- #1

- 462

- 1

- Thread starter radagast
- Start date

- #1

- 462

- 1

- #2

Hurkyl

Staff Emeritus

Science Advisor

Gold Member

- 14,916

- 19

Gaussian elimination isn't that bad, is it?

Important concerns include the size of your matrix, any qualitative properties it may have (such as sparse, symmetric, or banded), and what you want to do with the inverse.

For example, http://www.icos.ethz.ch/teaching/ml02/painless-conjugate-gradient.pdf [Broken] is well suited for solving the equation Ax=b when A is large sparse matrix, but it won't explicitly compute A inverse.

Important concerns include the size of your matrix, any qualitative properties it may have (such as sparse, symmetric, or banded), and what you want to do with the inverse.

For example, http://www.icos.ethz.ch/teaching/ml02/painless-conjugate-gradient.pdf [Broken] is well suited for solving the equation Ax=b when A is large sparse matrix, but it won't explicitly compute A inverse.

Last edited by a moderator:

- #3

chroot

Staff Emeritus

Science Advisor

Gold Member

- 10,226

- 34

- Warren

- #4

Integral

Staff Emeritus

Science Advisor

Gold Member

- 7,201

- 56

Rather then do a web search, I searched my book shelf. These routines are in Fortran (thats how old my book shelf is!) Not sure if this is the best format, but you will have to translate to a modern language anyway. I believe that the syntax is simple enough that you should (with fundamental programing knowledge) be able to figure out what is happening. I think I have edited out all of the OCR glitches, my software wanted to read an = as a z, so if there is a wierd z that makes no sense floating around replace it with an =.

These were lifted from__Elementary Numerical Analysis__ by Conte and de Boor.

I will do my best to answer any questions.

BTW

In the text preceding the inversion routine the authors say

INTEGER IPIVOT(N) , I,IP,J

REAL B(N) ,W(N,N) ,X(N), SUM

C****** I N PUT ******

C W, IPIVOT, N ARE AS ON OUTPUT FROM F ACT 0 R , APPLIED TO THE

C MATRIX A OF ORDER N.

C B IS AN N-VECTOR, GIVING THE RIGHT SIDE OF THE SYSTEM TO BE SOLVED.

C****** 0 U T PUT ******

C X IS THE N-VECTOR SATISFYING A*X . B.

C****** MET HOD ******

C ALGORITHM 4.4 IS USED, I.E., THE FACTORIZATION OF A CONTAINED IN

C W AND IPIVOT (AS GENERATED IN FACTOR) IS USED TO SOLVE A*X = B

C FOR X BY SOLVING TWO TRIANGULAR SYSTEMS.

C

IF (N .LE. 1) THEN

X(1) = B(1)/W(1,1)

RETURN

END IF

IP = IPIVOT(1)

X(1) = B(IP)

DO 15 I=2,N

SUM = 0.

DO 14 J=I,I-l

14 SUM = W(I,J)*X(J) + SUM

IP z IPIVOT(I)

15 X(I) = B(IP) - SUM

C

X(N) = X(N)/W(N,N)

DO 20 I=N-l,I,-1

SUM = 9.

DO 19 J=I+l,N

19 SUM = W(I,J)*X(J) + SUM

20 X(I) = (X(I) - SUM)/W(I,I)

RETURN

END

SUBROUTINE FACTOR ( W, N, D, IPIVOT, IFLAG )

INTEGER IFLAG,IPIVOT(N) , I,ISTAR,J,K

REAL D(N) ,W(N,N), AWIKOD,COLMAX,RATIO,ROWMAX,TEMP

**** INPUT ******

C W ARRAY OF SIZE (N,N) CONTAINING THE MATRIX A OF ORDER N TO BE

C FACTORED.

C N THE ORDER OF THE MATRIX

C***** WORK AREA ******

C D A REAL VECTOR OF LENGTH N, TO HOLD ROW SIZES

C***** 0UTPUT ******

C W ARRAY OF SIZE (N,N) CONTAINING THE LU FACTORIZATION OF P*A FOR

C SOME PERMUTATION MATRIX P SPECIFIED BY IPIVOT.

C IPIVOT INTEGER VECTOR OF LENGTH N INDICATING THAT ROW IPIVOT(K)

C WAS USED TO ELIMINATE X(K) , K-l,...,N .

C IFLAG AN INTEGER,

C = 1, IF AN EVEN NUMBER OF INTERCHANGES WAS CARRIED OUT,

C = -1,IF AN ODD NUMBER OF INTERCHANGES WAS CARRIED OUT,

C = 0, IF THE UPPER TRIANGULAR FACTOR HAS ONE OR MORE ZERO DIA-

C GONAL ENTRIES.

C THUS, DETERMINANT (A) - IFLAG*W(l,l)*...*W(N,N) .

C IF IFLAG .NE. 0, THEN THE LINEAR SYSTEM A*X 8 B CAN BE SOLVED FOR

C X BY A

C CALL SUBST (W, IPIVOT, B, N, X )

C**** METHOD ******

C THE PROGRAM FOLLOWS ALGORITHM 4.2, USING SCALED PARTIAL PIVOTING.

C

IFLAG = 1

C INITIALIZE IPIVOT, D

DO 9 I-l,N

IPIVOT(I) = I

ROWMAX = 0.

DO 5 J=l,N

5 ROWMAX = AMAXl(ROWMAX,ABS(W(I,J)))

IF (ROWMAX .Eq. 0.) THEN

IFLAG = 0

ROWMAX = 1.

END IF

9 D(I) = ROWMAX

IF (N .LE. 1) RETURN

C FACTORIZATION

DO 20 K=l,N-l

C DETERMINE PIVOT ROW, THE ROW ISTAR.

COLMAX = ABS(W(K,K)/D(K)

ISTAR = K

DO 13 I=K+l,N

AWIKOD = ABS(W(I,K))/D(I)

IF (AWIKOD .GT. COLMAX) THEN

COLMAX = AWIKOD

ISTAR = I

END IF

13 CONTINUE

IF (COLMAX .EO. 0.) THEN

IFLAG = 0

ELSE

IF (ISTAR .GT. K) THEN

C MAKE K THE PIVOT ROW BY INTERCHANGING IT WITH

C THE CHOSEN ROW ISTAR.

IFLAG = -IFLAG

I = IPIVOT(ISTAR)

IPIVOT(ISTAR) = IPIVOT(K)

IPIVOT(K) = I

TEMP = D(ISTAR)

D(ISTAR) = D(K)

D(K) = TEMP

DO 15 J=l,N

TEMP = W(ISTAR,J)

W(ISTAR,J) = W(K,J)

15 W(K,J) = TEMP

END IF

C ELIMINATE X(K) FROM ROWS K+l,...,N.

l6 DO 19 I=K+l,N

W(I,K) = W(I,K)/W(K,K)

RATIO = W(I,K)

DO 19 J=K=1,N

W(I,J) = W(I,J) - RATIO*W(K,J)

19 CONTINUE

END IF

20 CONTINUE

IF (W(N,N) .EQ. 0.) IFLAG -0

RETURN

END

C PROGRAM FOR CALCULATING THE INVERSE OF A GIVEN MATRIX

C CALLS FACT0R, SUBST.

PARAMETER NMAX=30,NMAXSQ=NMAX*NMAX

INTEGER I,IBEG,IFLAG,IPIVOT(NMAX) ,J,N,NSQ

REAL A(NMAXSQ) ,AINV(NMAXSQ) ,B(NMAX)

1 READ 501, N

501 FORMAT(I2)

IF (N .LT. 1 .OR. N .GT. NMAX) STOP

C READ IN MATRIX ROW BY ROW

NSQ = N*N

DO 10 I=I,N

18 READ 510, (A(J) ,J=I,NSQ,N)

510 FORMAT(5EI5.7)

C

CALL FACTOR ( A, N, B, IPIVOT, IFLAG )

IF (IFLAG .EQ. 0) THEN

PRINT 611

611 FORMAT('IMATRIX IS SINGULAR')

GO TO 1

END IF

DO 21 I=l,N

21 B(I) = 0.

IBEG = 1

DO 30 J=I,N

B(J) = 1.

CALL SUBST ( A, IPIVOT, B, N, AINV(IBEG) )

B(J) = 0.

38 IBEG = IBEG + N

PRINT 630

630 FORMAT('ITHE COMPUTED INVERSE IS '//)

DO 31 I=l,N

31 PRINT 631, I, (AINV(J) ,J-I,NSQ,N)

631 FORMAT('0ROW ',I2,8EI5.7/(7X,8EI5.7))

GO TO 1

END

These were lifted from

I will do my best to answer any questions.

BTW

In the text preceding the inversion routine the authors say

SUBROUTINE SUBST ( W, IPIVOT, B, N, X )..., we hasten to point out that there is usually no good reason for ever calculating the inverse.

INTEGER IPIVOT(N) , I,IP,J

REAL B(N) ,W(N,N) ,X(N), SUM

C****** I N PUT ******

C W, IPIVOT, N ARE AS ON OUTPUT FROM F ACT 0 R , APPLIED TO THE

C MATRIX A OF ORDER N.

C B IS AN N-VECTOR, GIVING THE RIGHT SIDE OF THE SYSTEM TO BE SOLVED.

C****** 0 U T PUT ******

C X IS THE N-VECTOR SATISFYING A*X . B.

C****** MET HOD ******

C ALGORITHM 4.4 IS USED, I.E., THE FACTORIZATION OF A CONTAINED IN

C W AND IPIVOT (AS GENERATED IN FACTOR) IS USED TO SOLVE A*X = B

C FOR X BY SOLVING TWO TRIANGULAR SYSTEMS.

C

IF (N .LE. 1) THEN

X(1) = B(1)/W(1,1)

RETURN

END IF

IP = IPIVOT(1)

X(1) = B(IP)

DO 15 I=2,N

SUM = 0.

DO 14 J=I,I-l

14 SUM = W(I,J)*X(J) + SUM

IP z IPIVOT(I)

15 X(I) = B(IP) - SUM

C

X(N) = X(N)/W(N,N)

DO 20 I=N-l,I,-1

SUM = 9.

DO 19 J=I+l,N

19 SUM = W(I,J)*X(J) + SUM

20 X(I) = (X(I) - SUM)/W(I,I)

RETURN

END

SUBROUTINE FACTOR ( W, N, D, IPIVOT, IFLAG )

INTEGER IFLAG,IPIVOT(N) , I,ISTAR,J,K

REAL D(N) ,W(N,N), AWIKOD,COLMAX,RATIO,ROWMAX,TEMP

**** INPUT ******

C W ARRAY OF SIZE (N,N) CONTAINING THE MATRIX A OF ORDER N TO BE

C FACTORED.

C N THE ORDER OF THE MATRIX

C***** WORK AREA ******

C D A REAL VECTOR OF LENGTH N, TO HOLD ROW SIZES

C***** 0UTPUT ******

C W ARRAY OF SIZE (N,N) CONTAINING THE LU FACTORIZATION OF P*A FOR

C SOME PERMUTATION MATRIX P SPECIFIED BY IPIVOT.

C IPIVOT INTEGER VECTOR OF LENGTH N INDICATING THAT ROW IPIVOT(K)

C WAS USED TO ELIMINATE X(K) , K-l,...,N .

C IFLAG AN INTEGER,

C = 1, IF AN EVEN NUMBER OF INTERCHANGES WAS CARRIED OUT,

C = -1,IF AN ODD NUMBER OF INTERCHANGES WAS CARRIED OUT,

C = 0, IF THE UPPER TRIANGULAR FACTOR HAS ONE OR MORE ZERO DIA-

C GONAL ENTRIES.

C THUS, DETERMINANT (A) - IFLAG*W(l,l)*...*W(N,N) .

C IF IFLAG .NE. 0, THEN THE LINEAR SYSTEM A*X 8 B CAN BE SOLVED FOR

C X BY A

C CALL SUBST (W, IPIVOT, B, N, X )

C**** METHOD ******

C THE PROGRAM FOLLOWS ALGORITHM 4.2, USING SCALED PARTIAL PIVOTING.

C

IFLAG = 1

C INITIALIZE IPIVOT, D

DO 9 I-l,N

IPIVOT(I) = I

ROWMAX = 0.

DO 5 J=l,N

5 ROWMAX = AMAXl(ROWMAX,ABS(W(I,J)))

IF (ROWMAX .Eq. 0.) THEN

IFLAG = 0

ROWMAX = 1.

END IF

9 D(I) = ROWMAX

IF (N .LE. 1) RETURN

C FACTORIZATION

DO 20 K=l,N-l

C DETERMINE PIVOT ROW, THE ROW ISTAR.

COLMAX = ABS(W(K,K)/D(K)

ISTAR = K

DO 13 I=K+l,N

AWIKOD = ABS(W(I,K))/D(I)

IF (AWIKOD .GT. COLMAX) THEN

COLMAX = AWIKOD

ISTAR = I

END IF

13 CONTINUE

IF (COLMAX .EO. 0.) THEN

IFLAG = 0

ELSE

IF (ISTAR .GT. K) THEN

C MAKE K THE PIVOT ROW BY INTERCHANGING IT WITH

C THE CHOSEN ROW ISTAR.

IFLAG = -IFLAG

I = IPIVOT(ISTAR)

IPIVOT(ISTAR) = IPIVOT(K)

IPIVOT(K) = I

TEMP = D(ISTAR)

D(ISTAR) = D(K)

D(K) = TEMP

DO 15 J=l,N

TEMP = W(ISTAR,J)

W(ISTAR,J) = W(K,J)

15 W(K,J) = TEMP

END IF

C ELIMINATE X(K) FROM ROWS K+l,...,N.

l6 DO 19 I=K+l,N

W(I,K) = W(I,K)/W(K,K)

RATIO = W(I,K)

DO 19 J=K=1,N

W(I,J) = W(I,J) - RATIO*W(K,J)

19 CONTINUE

END IF

20 CONTINUE

IF (W(N,N) .EQ. 0.) IFLAG -0

RETURN

END

C PROGRAM FOR CALCULATING THE INVERSE OF A GIVEN MATRIX

C CALLS FACT0R, SUBST.

PARAMETER NMAX=30,NMAXSQ=NMAX*NMAX

INTEGER I,IBEG,IFLAG,IPIVOT(NMAX) ,J,N,NSQ

REAL A(NMAXSQ) ,AINV(NMAXSQ) ,B(NMAX)

1 READ 501, N

501 FORMAT(I2)

IF (N .LT. 1 .OR. N .GT. NMAX) STOP

C READ IN MATRIX ROW BY ROW

NSQ = N*N

DO 10 I=I,N

18 READ 510, (A(J) ,J=I,NSQ,N)

510 FORMAT(5EI5.7)

C

CALL FACTOR ( A, N, B, IPIVOT, IFLAG )

IF (IFLAG .EQ. 0) THEN

PRINT 611

611 FORMAT('IMATRIX IS SINGULAR')

GO TO 1

END IF

DO 21 I=l,N

21 B(I) = 0.

IBEG = 1

DO 30 J=I,N

B(J) = 1.

CALL SUBST ( A, IPIVOT, B, N, AINV(IBEG) )

B(J) = 0.

38 IBEG = IBEG + N

PRINT 630

630 FORMAT('ITHE COMPUTED INVERSE IS '//)

DO 31 I=l,N

31 PRINT 631, I, (AINV(J) ,J-I,NSQ,N)

631 FORMAT('0ROW ',I2,8EI5.7/(7X,8EI5.7))

GO TO 1

END

Last edited:

- #5

chroot

Staff Emeritus

Science Advisor

Gold Member

- 10,226

- 34

Well, Fortran is pretty lame overall.

Numerical Recipes also has a lot of good material on the topic, in both fortran and C++. And, as quite a nice gift to the scientific community, the books are available in their entirety online:

http://ww.nr.com [Broken]

- Warren

Numerical Recipes also has a lot of good material on the topic, in both fortran and C++. And, as quite a nice gift to the scientific community, the books are available in their entirety online:

http://ww.nr.com [Broken]

- Warren

Last edited by a moderator:

- #6

Integral

Staff Emeritus

Science Advisor

Gold Member

- 7,201

- 56

The routines deserve a look inspite of the language, the above text generally presented very nice efficient routines, even if the language is archaic. :)

BTW: you may want to edit your link www usually works better then ww!

- #7

- 462

- 1

I was after the algorithm, I can always code it up into Java.

I throw and fire pottery for a hobby, and I've always had it in the back of my mind to write a glaze manipulation program. All the important properties of the ingredients are linearly translated into the glaze, so by using AX = B where A is a vector of ingredient percentages, X is the component/properties matrix, and B is the glaze component/properties vector, I can make changes in glaze properties (B), derive X

A mundane application, but it's a real pain when the thermal coefficient of expansion for the glaze and clay body don't match.

- #8

- 515

- 1

Originally posted by chroot

Numerical Recipes also has a lot of good material on the topic, in both fortran and C++. And, as quite a nice gift to the scientific community, the books are available in their entirety online:

http://www.nr.com

great site....

I'we wanted a good book on this ever since I studied this course in college (the profeesor was such a looser though I couldn't understand much from him).

My eternal gratitude....

PS: edited the above link....

- #9

mathwonk

Science Advisor

Homework Helper

- 11,011

- 1,208

aha! i got the center of the earth, or the beginning of time. the goal of all paleontologists!

- Last Post

- Replies
- 7

- Views
- 3K

- Last Post

- Replies
- 1

- Views
- 3K

- Last Post

- Replies
- 2

- Views
- 2K

- Last Post

- Replies
- 2

- Views
- 2K

- Last Post

- Replies
- 2

- Views
- 2K

- Last Post

- Replies
- 4

- Views
- 2K

- Last Post

- Replies
- 3

- Views
- 2K

- Last Post

- Replies
- 2

- Views
- 2K

- Last Post

- Replies
- 2

- Views
- 2K

- Replies
- 1

- Views
- 805