![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to wgeco.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10dc / cmd / matlab / src|
Return to wgeco.f CVS log | Up to [Research Unix] / researchv10dc / cmd / matlab / src |
1.1 ! root 1: SUBROUTINE WGECO(AR,AI,LDA,N,IPVT,RCOND,ZR,ZI) ! 2: INTEGER LDA,N,IPVT(1) ! 3: DOUBLE PRECISION AR(LDA,1),AI(LDA,1),ZR(1),ZI(1) ! 4: DOUBLE PRECISION RCOND ! 5: C ! 6: C WGECO FACTORS A DOUBLE-COMPLEX MATRIX BY GAUSSIAN ELIMINATION ! 7: C AND ESTIMATES THE CONDITION OF THE MATRIX. ! 8: C ! 9: C IF RCOND IS NOT NEEDED, WGEFA IS SLIGHTLY FASTER. ! 10: C TO SOLVE A*X = B , FOLLOW WGECO BY WGESL. ! 11: C TO COMPUTE INVERSE(A)*C , FOLLOW WGECO BY WGESL. ! 12: C TO COMPUTE DETERMINANT(A) , FOLLOW WGECO BY WGEDI. ! 13: C TO COMPUTE INVERSE(A) , FOLLOW WGECO BY WGEDI. ! 14: C ! 15: C ON ENTRY ! 16: C ! 17: C A DOUBLE-COMPLEX(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-COMPLEX(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 07/01/79 . ! 56: C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. ! 57: C ! 58: C SUBROUTINES AND FUNCTIONS ! 59: C ! 60: C LINPACK WGEFA ! 61: C BLAS WAXPY,WDOTC,WASUM ! 62: C FORTRAN DABS,DMAX1 ! 63: C ! 64: C INTERNAL VARIABLES ! 65: C ! 66: DOUBLE PRECISION WDOTCR,WDOTCI,EKR,EKI,TR,TI,WKR,WKI,WKMR,WKMI ! 67: DOUBLE PRECISION ANORM,S,WASUM,SM,YNORM,FLOP ! 68: INTEGER INFO,J,K,KB,KP1,L ! 69: C ! 70: DOUBLE PRECISION ZDUMR,ZDUMI ! 71: DOUBLE PRECISION CABS1 ! 72: CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI) ! 73: C ! 74: C COMPUTE 1-NORM OF A ! 75: C ! 76: ANORM = 0.0D0 ! 77: DO 10 J = 1, N ! 78: ANORM = DMAX1(ANORM,WASUM(N,AR(1,J),AI(1,J),1)) ! 79: 10 CONTINUE ! 80: C ! 81: C FACTOR ! 82: C ! 83: CALL WGEFA(AR,AI,LDA,N,IPVT,INFO) ! 84: C ! 85: C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . ! 86: C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND CTRANS(A)*Y = E . ! 87: C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A . ! 88: C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL ! 89: C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W = E . ! 90: C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. ! 91: C ! 92: C SOLVE CTRANS(U)*W = E ! 93: C ! 94: EKR = 1.0D0 ! 95: EKI = 0.0D0 ! 96: DO 20 J = 1, N ! 97: ZR(J) = 0.0D0 ! 98: ZI(J) = 0.0D0 ! 99: 20 CONTINUE ! 100: DO 110 K = 1, N ! 101: CALL WSIGN(EKR,EKI,-ZR(K),-ZI(K),EKR,EKI) ! 102: IF (CABS1(EKR-ZR(K),EKI-ZI(K)) ! 103: * .LE. CABS1(AR(K,K),AI(K,K))) GO TO 40 ! 104: S = CABS1(AR(K,K),AI(K,K)) ! 105: * /CABS1(EKR-ZR(K),EKI-ZI(K)) ! 106: CALL WRSCAL(N,S,ZR,ZI,1) ! 107: EKR = S*EKR ! 108: EKI = S*EKI ! 109: 40 CONTINUE ! 110: WKR = EKR - ZR(K) ! 111: WKI = EKI - ZI(K) ! 112: WKMR = -EKR - ZR(K) ! 113: WKMI = -EKI - ZI(K) ! 114: S = CABS1(WKR,WKI) ! 115: SM = CABS1(WKMR,WKMI) ! 116: IF (CABS1(AR(K,K),AI(K,K)) .EQ. 0.0D0) GO TO 50 ! 117: CALL WDIV(WKR,WKI,AR(K,K),-AI(K,K),WKR,WKI) ! 118: CALL WDIV(WKMR,WKMI,AR(K,K),-AI(K,K),WKMR,WKMI) ! 119: GO TO 60 ! 120: 50 CONTINUE ! 121: WKR = 1.0D0 ! 122: WKI = 0.0D0 ! 123: WKMR = 1.0D0 ! 124: WKMI = 0.0D0 ! 125: 60 CONTINUE ! 126: KP1 = K + 1 ! 127: IF (KP1 .GT. N) GO TO 100 ! 128: DO 70 J = KP1, N ! 129: CALL WMUL(WKMR,WKMI,AR(K,J),-AI(K,J),TR,TI) ! 130: SM = FLOP(SM + CABS1(ZR(J)+TR,ZI(J)+TI)) ! 131: CALL WAXPY(1,WKR,WKI,AR(K,J),-AI(K,J),1, ! 132: $ ZR(J),ZI(J),1) ! 133: S = FLOP(S + CABS1(ZR(J),ZI(J))) ! 134: 70 CONTINUE ! 135: IF (S .GE. SM) GO TO 90 ! 136: TR = WKMR - WKR ! 137: TI = WKMI - WKI ! 138: WKR = WKMR ! 139: WKI = WKMI ! 140: DO 80 J = KP1, N ! 141: CALL WAXPY(1,TR,TI,AR(K,J),-AI(K,J),1, ! 142: $ ZR(J),ZI(J),1) ! 143: 80 CONTINUE ! 144: 90 CONTINUE ! 145: 100 CONTINUE ! 146: ZR(K) = WKR ! 147: ZI(K) = WKI ! 148: 110 CONTINUE ! 149: S = 1.0D0/WASUM(N,ZR,ZI,1) ! 150: CALL WRSCAL(N,S,ZR,ZI,1) ! 151: C ! 152: C SOLVE CTRANS(L)*Y = W ! 153: C ! 154: DO 140 KB = 1, N ! 155: K = N + 1 - KB ! 156: IF (K .GE. N) GO TO 120 ! 157: ZR(K) = ZR(K) ! 158: * + WDOTCR(N-K,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),1) ! 159: ZI(K) = ZI(K) ! 160: * + WDOTCI(N-K,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1),1) ! 161: 120 CONTINUE ! 162: IF (CABS1(ZR(K),ZI(K)) .LE. 1.0D0) GO TO 130 ! 163: S = 1.0D0/CABS1(ZR(K),ZI(K)) ! 164: CALL WRSCAL(N,S,ZR,ZI,1) ! 165: 130 CONTINUE ! 166: L = IPVT(K) ! 167: TR = ZR(L) ! 168: TI = ZI(L) ! 169: ZR(L) = ZR(K) ! 170: ZI(L) = ZI(K) ! 171: ZR(K) = TR ! 172: ZI(K) = TI ! 173: 140 CONTINUE ! 174: S = 1.0D0/WASUM(N,ZR,ZI,1) ! 175: CALL WRSCAL(N,S,ZR,ZI,1) ! 176: C ! 177: YNORM = 1.0D0 ! 178: C ! 179: C SOLVE L*V = Y ! 180: C ! 181: DO 160 K = 1, N ! 182: L = IPVT(K) ! 183: TR = ZR(L) ! 184: TI = ZI(L) ! 185: ZR(L) = ZR(K) ! 186: ZI(L) = ZI(K) ! 187: ZR(K) = TR ! 188: ZI(K) = TI ! 189: IF (K .LT. N) ! 190: * CALL WAXPY(N-K,TR,TI,AR(K+1,K),AI(K+1,K),1,ZR(K+1),ZI(K+1), ! 191: * 1) ! 192: IF (CABS1(ZR(K),ZI(K)) .LE. 1.0D0) GO TO 150 ! 193: S = 1.0D0/CABS1(ZR(K),ZI(K)) ! 194: CALL WRSCAL(N,S,ZR,ZI,1) ! 195: YNORM = S*YNORM ! 196: 150 CONTINUE ! 197: 160 CONTINUE ! 198: S = 1.0D0/WASUM(N,ZR,ZI,1) ! 199: CALL WRSCAL(N,S,ZR,ZI,1) ! 200: YNORM = S*YNORM ! 201: C ! 202: C SOLVE U*Z = V ! 203: C ! 204: DO 200 KB = 1, N ! 205: K = N + 1 - KB ! 206: IF (CABS1(ZR(K),ZI(K)) ! 207: * .LE. CABS1(AR(K,K),AI(K,K))) GO TO 170 ! 208: S = CABS1(AR(K,K),AI(K,K)) ! 209: * /CABS1(ZR(K),ZI(K)) ! 210: CALL WRSCAL(N,S,ZR,ZI,1) ! 211: YNORM = S*YNORM ! 212: 170 CONTINUE ! 213: IF (CABS1(AR(K,K),AI(K,K)) .EQ. 0.0D0) GO TO 180 ! 214: CALL WDIV(ZR(K),ZI(K),AR(K,K),AI(K,K),ZR(K),ZI(K)) ! 215: 180 CONTINUE ! 216: IF (CABS1(AR(K,K),AI(K,K)) .NE. 0.0D0) GO TO 190 ! 217: ZR(K) = 1.0D0 ! 218: ZI(K) = 0.0D0 ! 219: 190 CONTINUE ! 220: TR = -ZR(K) ! 221: TI = -ZI(K) ! 222: CALL WAXPY(K-1,TR,TI,AR(1,K),AI(1,K),1,ZR(1),ZI(1),1) ! 223: 200 CONTINUE ! 224: C MAKE ZNORM = 1.0 ! 225: S = 1.0D0/WASUM(N,ZR,ZI,1) ! 226: CALL WRSCAL(N,S,ZR,ZI,1) ! 227: YNORM = S*YNORM ! 228: C ! 229: IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM ! 230: IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 ! 231: RETURN ! 232: END
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.