| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141 |
- *DECK DGEDI
- SUBROUTINE DGEDI (A, LDA, N, IPVT, DET, WORK, JOB)
- C***BEGIN PROLOGUE DGEDI
- C***PURPOSE Compute the determinant and inverse of a matrix using the
- C factors computed by DGECO or DGEFA.
- C***LIBRARY SLATEC (LINPACK)
- C***CATEGORY D3A1, D2A1
- C***TYPE DOUBLE PRECISION (SGEDI-S, DGEDI-D, CGEDI-C)
- C***KEYWORDS DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX
- C***AUTHOR Moler, C. B., (U. of New Mexico)
- C***DESCRIPTION
- C
- C DGEDI computes the determinant and inverse of a matrix
- C using the factors computed by DGECO or DGEFA.
- C
- C On Entry
- C
- C A DOUBLE PRECISION(LDA, N)
- C the output from DGECO or DGEFA.
- C
- C LDA INTEGER
- C the leading dimension of the array A .
- C
- C N INTEGER
- C the order of the matrix A .
- C
- C IPVT INTEGER(N)
- C the pivot vector from DGECO or DGEFA.
- C
- C WORK DOUBLE PRECISION(N)
- C work vector. Contents destroyed.
- C
- C JOB INTEGER
- C = 11 both determinant and inverse.
- C = 01 inverse only.
- C = 10 determinant only.
- C
- C On Return
- C
- C A inverse of original matrix if requested.
- C Otherwise unchanged.
- C
- C DET DOUBLE PRECISION(2)
- C determinant of original matrix if requested.
- C Otherwise not referenced.
- C Determinant = DET(1) * 10.0**DET(2)
- C with 1.0 .LE. ABS(DET(1)) .LT. 10.0
- C or DET(1) .EQ. 0.0 .
- C
- C Error Condition
- C
- C A division by zero will occur if the input factor contains
- C a zero on the diagonal and the inverse is requested.
- C It will not occur if the subroutines are called correctly
- C and if DGECO has set RCOND .GT. 0.0 or DGEFA has set
- C INFO .EQ. 0 .
- C
- C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
- C Stewart, LINPACK Users' Guide, SIAM, 1979.
- C***ROUTINES CALLED DAXPY, DSCAL, DSWAP
- C***REVISION HISTORY (YYMMDD)
- C 780814 DATE WRITTEN
- C 890531 Changed all specific intrinsics to generic. (WRB)
- C 890831 Modified array declarations. (WRB)
- C 890831 REVISION DATE from Version 3.2
- C 891214 Prologue converted to Version 4.0 format. (BAB)
- C 900326 Removed duplicate information from DESCRIPTION section.
- C (WRB)
- C 920501 Reformatted the REFERENCES section. (WRB)
- C***END PROLOGUE DGEDI
- INTEGER LDA,N,IPVT(*),JOB
- DOUBLE PRECISION A(LDA,*),DET(2),WORK(*)
- C
- DOUBLE PRECISION T
- DOUBLE PRECISION TEN
- INTEGER I,J,K,KB,KP1,L,NM1
- C***FIRST EXECUTABLE STATEMENT DGEDI
- C
- C COMPUTE DETERMINANT
- C
- IF (JOB/10 .EQ. 0) GO TO 70
- DET(1) = 1.0D0
- DET(2) = 0.0D0
- TEN = 10.0D0
- DO 50 I = 1, N
- IF (IPVT(I) .NE. I) DET(1) = -DET(1)
- DET(1) = A(I,I)*DET(1)
- IF (DET(1) .EQ. 0.0D0) GO TO 60
- 10 IF (ABS(DET(1)) .GE. 1.0D0) GO TO 20
- DET(1) = TEN*DET(1)
- DET(2) = DET(2) - 1.0D0
- GO TO 10
- 20 CONTINUE
- 30 IF (ABS(DET(1)) .LT. TEN) GO TO 40
- DET(1) = DET(1)/TEN
- DET(2) = DET(2) + 1.0D0
- GO TO 30
- 40 CONTINUE
- 50 CONTINUE
- 60 CONTINUE
- 70 CONTINUE
- C
- C COMPUTE INVERSE(U)
- C
- IF (MOD(JOB,10) .EQ. 0) GO TO 150
- DO 100 K = 1, N
- A(K,K) = 1.0D0/A(K,K)
- T = -A(K,K)
- CALL DSCAL(K-1,T,A(1,K),1)
- KP1 = K + 1
- IF (N .LT. KP1) GO TO 90
- DO 80 J = KP1, N
- T = A(K,J)
- A(K,J) = 0.0D0
- CALL DAXPY(K,T,A(1,K),1,A(1,J),1)
- 80 CONTINUE
- 90 CONTINUE
- 100 CONTINUE
- C
- C FORM INVERSE(U)*INVERSE(L)
- C
- NM1 = N - 1
- IF (NM1 .LT. 1) GO TO 140
- DO 130 KB = 1, NM1
- K = N - KB
- KP1 = K + 1
- DO 110 I = KP1, N
- WORK(I) = A(I,K)
- A(I,K) = 0.0D0
- 110 CONTINUE
- DO 120 J = KP1, N
- T = WORK(J)
- CALL DAXPY(N,T,A(1,J),1,A(1,K),1)
- 120 CONTINUE
- L = IPVT(K)
- IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1)
- 130 CONTINUE
- 140 CONTINUE
- 150 CONTINUE
- RETURN
- END
|
