| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145 |
- *DECK STRDI
- SUBROUTINE STRDI (T, LDT, N, DET, JOB, INFO)
- C***BEGIN PROLOGUE STRDI
- C***PURPOSE Compute the determinant and inverse of a triangular matrix.
- C***LIBRARY SLATEC (LINPACK)
- C***CATEGORY D2A3, D3A3
- C***TYPE SINGLE PRECISION (STRDI-S, DTRDI-D, CTRDI-C)
- C***KEYWORDS DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX,
- C TRIANGULAR
- C***AUTHOR Moler, C. B., (U. of New Mexico)
- C***DESCRIPTION
- C
- C STRDI computes the determinant and inverse of a real
- C triangular matrix.
- C
- C On Entry
- C
- C T REAL(LDT,N)
- C T contains the triangular matrix. The zero
- C elements of the matrix are not referenced, and
- C the corresponding elements of the array can be
- C used to store other information.
- C
- C LDT INTEGER
- C LDT is the leading dimension of the array T.
- C
- C N INTEGER
- C N is the order of the system.
- C
- C JOB INTEGER
- C = 010 no det, inverse of lower triangular.
- C = 011 no det, inverse of upper triangular.
- C = 100 det, no inverse.
- C = 110 det, inverse of lower triangular.
- C = 111 det, inverse of upper triangular.
- C
- C On Return
- C
- C T inverse of original matrix if requested.
- C Otherwise unchanged.
- C
- C DET REAL(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 INFO INTEGER
- C INFO contains zero if the system is nonsingular
- C and the inverse is requested.
- C Otherwise INFO contains the index of
- C a zero diagonal element of T.
- 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 SAXPY, SSCAL
- C***REVISION HISTORY (YYMMDD)
- C 780814 DATE WRITTEN
- 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 STRDI
- INTEGER LDT,N,JOB,INFO
- REAL T(LDT,*),DET(2)
- C
- REAL TEMP
- REAL TEN
- INTEGER I,J,K,KB,KM1,KP1
- C***FIRST EXECUTABLE STATEMENT STRDI
- C
- C COMPUTE DETERMINANT
- C
- IF (JOB/100 .EQ. 0) GO TO 70
- DET(1) = 1.0E0
- DET(2) = 0.0E0
- TEN = 10.0E0
- DO 50 I = 1, N
- DET(1) = T(I,I)*DET(1)
- IF (DET(1) .EQ. 0.0E0) GO TO 60
- 10 IF (ABS(DET(1)) .GE. 1.0E0) GO TO 20
- DET(1) = TEN*DET(1)
- DET(2) = DET(2) - 1.0E0
- 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.0E0
- GO TO 30
- 40 CONTINUE
- 50 CONTINUE
- 60 CONTINUE
- 70 CONTINUE
- C
- C COMPUTE INVERSE OF UPPER TRIANGULAR
- C
- IF (MOD(JOB/10,10) .EQ. 0) GO TO 170
- IF (MOD(JOB,10) .EQ. 0) GO TO 120
- DO 100 K = 1, N
- INFO = K
- IF (T(K,K) .EQ. 0.0E0) GO TO 110
- T(K,K) = 1.0E0/T(K,K)
- TEMP = -T(K,K)
- CALL SSCAL(K-1,TEMP,T(1,K),1)
- KP1 = K + 1
- IF (N .LT. KP1) GO TO 90
- DO 80 J = KP1, N
- TEMP = T(K,J)
- T(K,J) = 0.0E0
- CALL SAXPY(K,TEMP,T(1,K),1,T(1,J),1)
- 80 CONTINUE
- 90 CONTINUE
- 100 CONTINUE
- INFO = 0
- 110 CONTINUE
- GO TO 160
- 120 CONTINUE
- C
- C COMPUTE INVERSE OF LOWER TRIANGULAR
- C
- DO 150 KB = 1, N
- K = N + 1 - KB
- INFO = K
- IF (T(K,K) .EQ. 0.0E0) GO TO 180
- T(K,K) = 1.0E0/T(K,K)
- TEMP = -T(K,K)
- IF (K .NE. N) CALL SSCAL(N-K,TEMP,T(K+1,K),1)
- KM1 = K - 1
- IF (KM1 .LT. 1) GO TO 140
- DO 130 J = 1, KM1
- TEMP = T(K,J)
- T(K,J) = 0.0E0
- CALL SAXPY(N-K+1,TEMP,T(K,K),1,T(K,J),1)
- 130 CONTINUE
- 140 CONTINUE
- 150 CONTINUE
- INFO = 0
- 160 CONTINUE
- 170 CONTINUE
- 180 CONTINUE
- RETURN
- END
|
