1 SUBROUTINE DTRCO(T,LDT,N,RCOND,Z,JOB) 2C***BEGIN PROLOGUE DTRCO 3C***DATE WRITTEN 780814 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000330 Modified array declarations. (JEC) 7C***CATEGORY NO. D2A3 8C***KEYWORDS CONDITION,DOUBLE PRECISION,FACTOR,LINEAR ALGEBRA,LINPACK, 9C MATRIX,TRIANGULAR 10C***AUTHOR MOLER, C. B., (U. OF NEW MEXICO) 11C***PURPOSE Estimates the condition of a double precision TRIANGULAR 12C matrix. 13C***DESCRIPTION 14C 15C DTRCO estimates the condition of a double precision triangular 16C matrix. 17C 18C On Entry 19C 20C T DOUBLE PRECISION(LDT,N) 21C T contains the triangular matrix. The zero 22C elements of the matrix are not referenced, and 23C the corresponding elements of the array can be 24C used to store other information. 25C 26C LDT INTEGER 27C LDT is the leading dimension of the array T. 28C 29C N INTEGER 30C N is the order of the system. 31C 32C JOB INTEGER 33C = 0 T is lower triangular. 34C = nonzero T is upper triangular. 35C 36C On Return 37C 38C RCOND DOUBLE PRECISION 39C an estimate of the reciprocal condition of T . 40C For the system T*X = B , relative perturbations 41C in T and B of size EPSILON may cause 42C relative perturbations in X of size EPSILON/RCOND . 43C If RCOND is so small that the logical expression 44C 1.0 + RCOND .EQ. 1.0 45C is true, then T may be singular to working 46C precision. In particular, RCOND is zero if 47C exact singularity is detected or the estimate 48C underflows. 49C 50C Z DOUBLE PRECISION(N) 51C a work vector whose contents are usually unimportant. 52C If T is close to a singular matrix, then Z is 53C an approximate null vector in the sense that 54C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . 55C 56C LINPACK. This version dated 08/14/78 . 57C Cleve Moler, University of New Mexico, Argonne National Lab. 58C 59C Subroutines and Functions 60C 61C BLAS DAXPY,DSCAL,DASUM 62C Fortran DABS,DMAX1,DSIGN 63C***REFERENCES DONGARRA J.J., BUNCH J.R., MOLER C.B., STEWART G.W., 64C *LINPACK USERS GUIDE*, SIAM, 1979. 65C***ROUTINES CALLED DASUM,DAXPY,DSCAL 66C***END PROLOGUE DTRCO 67 INTEGER LDT,N,JOB 68 DOUBLE PRECISION T(LDT,*),Z(*) 69 DOUBLE PRECISION RCOND 70C 71 DOUBLE PRECISION W,WK,WKM,EK 72 DOUBLE PRECISION TNORM,YNORM,S,SM,DASUM 73 INTEGER I1,J,J1,J2,K,KK,L 74 LOGICAL LOWER 75C***FIRST EXECUTABLE STATEMENT DTRCO 76 LOWER = JOB .EQ. 0 77C 78C COMPUTE 1-NORM OF T 79C 80 TNORM = 0.0D0 81 DO 10 J = 1, N 82 L = J 83 IF (LOWER) L = N + 1 - J 84 I1 = 1 85 IF (LOWER) I1 = J 86 TNORM = DMAX1(TNORM,DASUM(L,T(I1,J),1)) 87 10 CONTINUE 88C 89C RCOND = 1/(NORM(T)*(ESTIMATE OF NORM(INVERSE(T)))) . 90C ESTIMATE = NORM(Z)/NORM(Y) WHERE T*Z = Y AND TRANS(T)*Y = E . 91C TRANS(T) IS THE TRANSPOSE OF T . 92C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL 93C GROWTH IN THE ELEMENTS OF Y . 94C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. 95C 96C SOLVE TRANS(T)*Y = E 97C 98 EK = 1.0D0 99 DO 20 J = 1, N 100 Z(J) = 0.0D0 101 20 CONTINUE 102 DO 100 KK = 1, N 103 K = KK 104 IF (LOWER) K = N + 1 - KK 105 IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) 106 IF (DABS(EK-Z(K)) .LE. DABS(T(K,K))) GO TO 30 107 S = DABS(T(K,K))/DABS(EK-Z(K)) 108 CALL DSCAL(N,S,Z,1) 109 EK = S*EK 110 30 CONTINUE 111 WK = EK - Z(K) 112 WKM = -EK - Z(K) 113 S = DABS(WK) 114 SM = DABS(WKM) 115 IF (T(K,K) .EQ. 0.0D0) GO TO 40 116 WK = WK/T(K,K) 117 WKM = WKM/T(K,K) 118 GO TO 50 119 40 CONTINUE 120 WK = 1.0D0 121 WKM = 1.0D0 122 50 CONTINUE 123 IF (KK .EQ. N) GO TO 90 124 J1 = K + 1 125 IF (LOWER) J1 = 1 126 J2 = N 127 IF (LOWER) J2 = K - 1 128 DO 60 J = J1, J2 129 SM = SM + DABS(Z(J)+WKM*T(K,J)) 130 Z(J) = Z(J) + WK*T(K,J) 131 S = S + DABS(Z(J)) 132 60 CONTINUE 133 IF (S .GE. SM) GO TO 80 134 W = WKM - WK 135 WK = WKM 136 DO 70 J = J1, J2 137 Z(J) = Z(J) + W*T(K,J) 138 70 CONTINUE 139 80 CONTINUE 140 90 CONTINUE 141 Z(K) = WK 142 100 CONTINUE 143 S = 1.0D0/DASUM(N,Z,1) 144 CALL DSCAL(N,S,Z,1) 145C 146 YNORM = 1.0D0 147C 148C SOLVE T*Z = Y 149C 150 DO 130 KK = 1, N 151 K = N + 1 - KK 152 IF (LOWER) K = KK 153 IF (DABS(Z(K)) .LE. DABS(T(K,K))) GO TO 110 154 S = DABS(T(K,K))/DABS(Z(K)) 155 CALL DSCAL(N,S,Z,1) 156 YNORM = S*YNORM 157 110 CONTINUE 158 IF (T(K,K) .NE. 0.0D0) Z(K) = Z(K)/T(K,K) 159 IF (T(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 160 I1 = 1 161 IF (LOWER) I1 = K + 1 162 IF (KK .GE. N) GO TO 120 163 W = -Z(K) 164 CALL DAXPY(N-KK,W,T(I1,K),1,Z(I1),1) 165 120 CONTINUE 166 130 CONTINUE 167C MAKE ZNORM = 1.0 168 S = 1.0D0/DASUM(N,Z,1) 169 CALL DSCAL(N,S,Z,1) 170 YNORM = S*YNORM 171C 172 IF (TNORM .NE. 0.0D0) RCOND = YNORM/TNORM 173 IF (TNORM .EQ. 0.0D0) RCOND = 0.0D0 174 RETURN 175 END 176