|
2329
|
1 SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) |
|
|
2 INTEGER LDA,N,IPVT(1) |
|
|
3 DOUBLE PRECISION A(LDA,1),Z(1) |
|
|
4 DOUBLE PRECISION RCOND |
|
|
5 C |
|
|
6 C DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION |
|
|
7 C AND ESTIMATES THE CONDITION OF THE MATRIX. |
|
|
8 C |
|
|
9 C IF RCOND IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER. |
|
|
10 C TO SOLVE A*X = B , FOLLOW DGECO BY DGESL. |
|
|
11 C TO COMPUTE INVERSE(A)*C , FOLLOW DGECO BY DGESL. |
|
|
12 C TO COMPUTE DETERMINANT(A) , FOLLOW DGECO BY DGEDI. |
|
|
13 C TO COMPUTE INVERSE(A) , FOLLOW DGECO BY DGEDI. |
|
|
14 C |
|
|
15 C ON ENTRY |
|
|
16 C |
|
|
17 C A DOUBLE PRECISION(LDA, N) |
|
|
18 C THE MATRIX TO BE FACTORED. |
|
|
19 C |
|
|
20 C LDA INTEGER |
|
|
21 C THE LEADING DIMENSION OF THE ARRAY A . |
|
|
22 C |
|
|
23 C N INTEGER |
|
|
24 C THE ORDER OF THE MATRIX A . |
|
|
25 C |
|
|
26 C ON RETURN |
|
|
27 C |
|
|
28 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS |
|
|
29 C WHICH WERE USED TO OBTAIN IT. |
|
|
30 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE |
|
|
31 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER |
|
|
32 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. |
|
|
33 C |
|
|
34 C IPVT INTEGER(N) |
|
|
35 C AN INTEGER VECTOR OF PIVOT INDICES. |
|
|
36 C |
|
|
37 C RCOND DOUBLE PRECISION |
|
|
38 C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . |
|
|
39 C FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS |
|
|
40 C IN A AND B OF SIZE EPSILON MAY CAUSE |
|
|
41 C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . |
|
|
42 C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION |
|
|
43 C 1.0 + RCOND .EQ. 1.0 |
|
|
44 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING |
|
|
45 C PRECISION. IN PARTICULAR, RCOND IS ZERO IF |
|
|
46 C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE |
|
|
47 C UNDERFLOWS. |
|
|
48 C |
|
|
49 C Z DOUBLE PRECISION(N) |
|
|
50 C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. |
|
|
51 C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS |
|
|
52 C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT |
|
|
53 C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . |
|
|
54 C |
|
|
55 C LINPACK. THIS VERSION DATED 08/14/78 . |
|
|
56 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. |
|
|
57 C |
|
|
58 C SUBROUTINES AND FUNCTIONS |
|
|
59 C |
|
|
60 C LINPACK DGEFA |
|
|
61 C BLAS DAXPY,DDOT,DSCAL,DASUM |
|
|
62 C FORTRAN DABS,DMAX1,DSIGN |
|
|
63 C |
|
|
64 C INTERNAL VARIABLES |
|
|
65 C |
|
|
66 DOUBLE PRECISION DDOT,EK,T,WK,WKM |
|
|
67 DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM |
|
|
68 INTEGER INFO,J,K,KB,KP1,L |
|
|
69 C |
|
|
70 C |
|
|
71 C COMPUTE 1-NORM OF A |
|
|
72 C |
|
|
73 ANORM = 0.0D0 |
|
|
74 DO 10 J = 1, N |
|
|
75 ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1)) |
|
|
76 10 CONTINUE |
|
|
77 C |
|
|
78 C FACTOR |
|
|
79 C |
|
|
80 CALL DGEFA(A,LDA,N,IPVT,INFO) |
|
|
81 C |
|
|
82 C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . |
|
|
83 C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . |
|
|
84 C TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE |
|
|
85 C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE |
|
|
86 C TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID |
|
|
87 C OVERFLOW. |
|
|
88 C |
|
|
89 C SOLVE TRANS(U)*W = E |
|
|
90 C |
|
|
91 EK = 1.0D0 |
|
|
92 DO 20 J = 1, N |
|
|
93 Z(J) = 0.0D0 |
|
|
94 20 CONTINUE |
|
|
95 DO 100 K = 1, N |
|
|
96 IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) |
|
|
97 IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 |
|
|
98 S = DABS(A(K,K))/DABS(EK-Z(K)) |
|
|
99 CALL DSCAL(N,S,Z,1) |
|
|
100 EK = S*EK |
|
|
101 30 CONTINUE |
|
|
102 WK = EK - Z(K) |
|
|
103 WKM = -EK - Z(K) |
|
|
104 S = DABS(WK) |
|
|
105 SM = DABS(WKM) |
|
|
106 IF (A(K,K) .EQ. 0.0D0) GO TO 40 |
|
|
107 WK = WK/A(K,K) |
|
|
108 WKM = WKM/A(K,K) |
|
|
109 GO TO 50 |
|
|
110 40 CONTINUE |
|
|
111 WK = 1.0D0 |
|
|
112 WKM = 1.0D0 |
|
|
113 50 CONTINUE |
|
|
114 KP1 = K + 1 |
|
|
115 IF (KP1 .GT. N) GO TO 90 |
|
|
116 DO 60 J = KP1, N |
|
|
117 SM = SM + DABS(Z(J)+WKM*A(K,J)) |
|
|
118 Z(J) = Z(J) + WK*A(K,J) |
|
|
119 S = S + DABS(Z(J)) |
|
|
120 60 CONTINUE |
|
|
121 IF (S .GE. SM) GO TO 80 |
|
|
122 T = WKM - WK |
|
|
123 WK = WKM |
|
|
124 DO 70 J = KP1, N |
|
|
125 Z(J) = Z(J) + T*A(K,J) |
|
|
126 70 CONTINUE |
|
|
127 80 CONTINUE |
|
|
128 90 CONTINUE |
|
|
129 Z(K) = WK |
|
|
130 100 CONTINUE |
|
|
131 S = 1.0D0/DASUM(N,Z,1) |
|
|
132 CALL DSCAL(N,S,Z,1) |
|
|
133 C |
|
|
134 C SOLVE TRANS(L)*Y = W |
|
|
135 C |
|
|
136 DO 120 KB = 1, N |
|
|
137 K = N + 1 - KB |
|
|
138 IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) |
|
|
139 IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 |
|
|
140 S = 1.0D0/DABS(Z(K)) |
|
|
141 CALL DSCAL(N,S,Z,1) |
|
|
142 110 CONTINUE |
|
|
143 L = IPVT(K) |
|
|
144 T = Z(L) |
|
|
145 Z(L) = Z(K) |
|
|
146 Z(K) = T |
|
|
147 120 CONTINUE |
|
|
148 S = 1.0D0/DASUM(N,Z,1) |
|
|
149 CALL DSCAL(N,S,Z,1) |
|
|
150 C |
|
|
151 YNORM = 1.0D0 |
|
|
152 C |
|
|
153 C SOLVE L*V = Y |
|
|
154 C |
|
|
155 DO 140 K = 1, N |
|
|
156 L = IPVT(K) |
|
|
157 T = Z(L) |
|
|
158 Z(L) = Z(K) |
|
|
159 Z(K) = T |
|
|
160 IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) |
|
|
161 IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 |
|
|
162 S = 1.0D0/DABS(Z(K)) |
|
|
163 CALL DSCAL(N,S,Z,1) |
|
|
164 YNORM = S*YNORM |
|
|
165 130 CONTINUE |
|
|
166 140 CONTINUE |
|
|
167 S = 1.0D0/DASUM(N,Z,1) |
|
|
168 CALL DSCAL(N,S,Z,1) |
|
|
169 YNORM = S*YNORM |
|
|
170 C |
|
|
171 C SOLVE U*Z = V |
|
|
172 C |
|
|
173 DO 160 KB = 1, N |
|
|
174 K = N + 1 - KB |
|
|
175 IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 |
|
|
176 S = DABS(A(K,K))/DABS(Z(K)) |
|
|
177 CALL DSCAL(N,S,Z,1) |
|
|
178 YNORM = S*YNORM |
|
|
179 150 CONTINUE |
|
|
180 IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) |
|
|
181 IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 |
|
|
182 T = -Z(K) |
|
|
183 CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) |
|
|
184 160 CONTINUE |
|
|
185 C MAKE ZNORM = 1.0 |
|
|
186 S = 1.0D0/DASUM(N,Z,1) |
|
|
187 CALL DSCAL(N,S,Z,1) |
|
|
188 YNORM = S*YNORM |
|
|
189 C |
|
|
190 IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM |
|
|
191 IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 |
|
|
192 RETURN |
|
|
193 END |
