minv.f
164 lines
| 3.6 KiB
| text/x-fortran
|
FortranFixedLexer
| r0 | C $Id: minv.f 3304 2011-01-17 15:25:59Z brideout $ | |||
| C | ||||
| SUBROUTINE MINV(A,N,D,L,M) | ||||
| C | ||||
| C Inverts general matrix A (overwrites A) using standard | ||||
| C gauss-jordan method. The determinant is also calculated. a | ||||
| C determinant of zero indicates that the matrix is singular. | ||||
| C | ||||
| C Input: | ||||
| C n - order of matrix a | ||||
| C l - work vector of length n | ||||
| C m - work vector of length n | ||||
| C | ||||
| C Input, Output: | ||||
| C a - input matrix, destroyed in computation and replaced by | ||||
| C resultant inverse. | ||||
| C | ||||
| C Output: | ||||
| C d - resultant determinant | ||||
| C | ||||
| C .. Scalar Arguments .. | ||||
| DOUBLE PRECISION D | ||||
| INTEGER N | ||||
| C .. | ||||
| C .. Array Arguments .. | ||||
| DOUBLE PRECISION A(*) | ||||
| INTEGER L(*),M(*) | ||||
| C .. | ||||
| C .. Local Scalars .. | ||||
| DOUBLE PRECISION BIGA,HOLD | ||||
| INTEGER I,IJ,IK,IZ,J,JI,JK,JP,JQ,JR,K,KI,KJ,KK,NK | ||||
| C .. | ||||
| C .. Intrinsic Functions .. | ||||
| INTRINSIC ABS | ||||
| C .. | ||||
| C search for largest element | ||||
| D = 1.0D0 | ||||
| NK = -N | ||||
| DO 90 K = 1,N | ||||
| NK = NK + N | ||||
| L(K) = K | ||||
| M(K) = K | ||||
| KK = NK + K | ||||
| BIGA = A(KK) | ||||
| DO 20 J = K,N | ||||
| IZ = N*(J-1) | ||||
| DO 10 I = K,N | ||||
| IJ = IZ + I | ||||
| IF (ABS(BIGA).LT.ABS(A(IJ))) THEN | ||||
| BIGA = A(IJ) | ||||
| L(K) = I | ||||
| M(K) = J | ||||
| END IF | ||||
| 10 CONTINUE | ||||
| 20 CONTINUE | ||||
| C | ||||
| C interchange rows | ||||
| C | ||||
| J = L(K) | ||||
| IF (J.GT.K) THEN | ||||
| KI = K - N | ||||
| DO 30 I = 1,N | ||||
| KI = KI + N | ||||
| HOLD = -A(KI) | ||||
| JI = KI - K + J | ||||
| A(KI) = A(JI) | ||||
| A(JI) = HOLD | ||||
| 30 CONTINUE | ||||
| END IF | ||||
| C | ||||
| C interchange columns | ||||
| C | ||||
| I = M(K) | ||||
| IF (I.GT.K) THEN | ||||
| JP = N*(I-1) | ||||
| DO 40 J = 1,N | ||||
| JK = NK + J | ||||
| JI = JP + J | ||||
| HOLD = -A(JK) | ||||
| A(JK) = A(JI) | ||||
| A(JI) = HOLD | ||||
| 40 CONTINUE | ||||
| END IF | ||||
| C | ||||
| C divide column by minus pivot (value of pivot element is | ||||
| C contained in biga) | ||||
| C | ||||
| IF (BIGA.EQ.0) GO TO 130 | ||||
| DO 50 I = 1,N | ||||
| IF (I.NE.K) THEN | ||||
| IK = NK + I | ||||
| A(IK) = A(IK)/(-BIGA) | ||||
| END IF | ||||
| 50 CONTINUE | ||||
| C | ||||
| C reduce matrix | ||||
| C | ||||
| DO 70 I = 1,N | ||||
| IK = NK + I | ||||
| HOLD = A(IK) | ||||
| IJ = I - N | ||||
| DO 60 J = 1,N | ||||
| IJ = IJ + N | ||||
| IF (I.NE.K) THEN | ||||
| IF (J.NE.K) THEN | ||||
| KJ = IJ - I + K | ||||
| A(IJ) = HOLD*A(KJ) + A(IJ) | ||||
| END IF | ||||
| END IF | ||||
| 60 CONTINUE | ||||
| 70 CONTINUE | ||||
| C | ||||
| C divide row by pivot | ||||
| C | ||||
| KJ = K - N | ||||
| DO 80 J = 1,N | ||||
| KJ = KJ + N | ||||
| IF (J.NE.K) A(KJ) = A(KJ)/BIGA | ||||
| 80 CONTINUE | ||||
| C | ||||
| C product of pivots | ||||
| C | ||||
| D = D*BIGA | ||||
| C | ||||
| C replace pivot by reciprocal | ||||
| C | ||||
| A(KK) = 1.0D0/BIGA | ||||
| 90 CONTINUE | ||||
| C | ||||
| C final row and column interchange | ||||
| C | ||||
| K = N | ||||
| 100 CONTINUE | ||||
| K = (K-1) | ||||
| IF (K.LE.0) GO TO 140 | ||||
| I = L(K) | ||||
| IF (I.GT.K) THEN | ||||
| JQ = N*(K-1) | ||||
| JR = N*(I-1) | ||||
| DO 110 J = 1,N | ||||
| JK = JQ + J | ||||
| HOLD = A(JK) | ||||
| JI = JR + J | ||||
| A(JK) = -A(JI) | ||||
| A(JI) = HOLD | ||||
| 110 CONTINUE | ||||
| END IF | ||||
| J = M(K) | ||||
| IF (J.GT.K) THEN | ||||
| KI = K - N | ||||
| DO 120 I = 1,N | ||||
| KI = KI + N | ||||
| HOLD = A(KI) | ||||
| JI = KI - K + J | ||||
| A(KI) = -A(JI) | ||||
| A(JI) = HOLD | ||||
| 120 CONTINUE | ||||
| END IF | ||||
| GO TO 100 | ||||
| 130 CONTINUE | ||||
| D = 0.0D0 | ||||
| 140 CONTINUE | ||||
| RETURN | ||||
| END | ||||
