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