![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to matfn1.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 matfn1.f CVS log | Up to [Research Unix] / researchv10dc / cmd / matlab / src |
1.1 ! root 1: SUBROUTINE MATFN1 ! 2: C ! 3: C EVALUATE FUNCTIONS INVOLVING GAUSSIAN ELIMINATION ! 4: C ! 5: DOUBLE PRECISION STKR(50005),STKI(50005) ! 6: INTEGER IDSTK(4,48),LSTK(48),MSTK(48),NSTK(48),VSIZE,LSIZE,BOT,TOP ! 7: INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),RIO,WIO,RTE,WTE,HIO ! 8: INTEGER SYM,SYN(4),BUF(256),CHAR,FLP(2),FIN,FUN,LHS,RHS,RAN(2) ! 9: COMMON /VSTK/ STKR,STKI,IDSTK,LSTK,MSTK,NSTK,VSIZE,LSIZE,BOT,TOP ! 10: COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,RIO,WIO,RTE,WTE,HIO ! 11: COMMON /COM/ SYM,SYN,BUF,CHAR,FLP,FIN,FUN,LHS,RHS,RAN ! 12: DOUBLE PRECISION DTR(2),DTI(2),SR,SI,RCOND,T,T0,T1,FLOP,EPS,WASUM ! 13: C ! 14: IF (DDT .EQ. 1) WRITE(WTE,100) FIN ! 15: 100 FORMAT(0X,'MATFN1',I4) ! 16: C ! 17: L = LSTK(TOP) ! 18: M = MSTK(TOP) ! 19: N = NSTK(TOP) ! 20: IF (FIN .EQ. -1) GO TO 10 ! 21: IF (FIN .EQ. -2) GO TO 20 ! 22: GO TO (30,40,50,60,70,80,85),FIN ! 23: C ! 24: C MATRIX RIGHT DIVISION, A/A2 ! 25: 10 L2 = LSTK(TOP+1) ! 26: M2 = MSTK(TOP+1) ! 27: N2 = NSTK(TOP+1) ! 28: IF (M2 .NE. N2) CALL ERROR(20) ! 29: IF (ERR .GT. 0) RETURN ! 30: IF (M*N .EQ. 1) GO TO 16 ! 31: IF (N .NE. N2) CALL ERROR(11) ! 32: IF (ERR .GT. 0) RETURN ! 33: L3 = L2 + M2*N2 ! 34: ERR = L3+N2 - LSTK(BOT) ! 35: IF (ERR .GT. 0) CALL ERROR(17) ! 36: IF (ERR .GT. 0) RETURN ! 37: CALL WGECO(STKR(L2),STKI(L2),M2,N2,BUF,RCOND,STKR(L3),STKI(L3)) ! 38: IF (RCOND .EQ. 0.0D0) CALL ERROR(19) ! 39: IF (ERR .GT. 0) RETURN ! 40: T = FLOP(1.0D0 + RCOND) ! 41: IF (T.EQ.1.0D0 .AND. FUN.NE.21) WRITE(WTE,11) RCOND ! 42: IF (T.EQ.1.0D0 .AND. FUN.NE.21 .AND. WIO.NE.0) WRITE(WIO,11) RCOND ! 43: 11 FORMAT(0X,'WARNING.' ! 44: $ /0X,'MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.' ! 45: $ /0X,'RESULTS MAY BE INACCURATE. RCOND =', 1PD13.4/) ! 46: IF (T.EQ.1.0D0 .AND. FUN.EQ.21) WRITE(WTE,12) RCOND ! 47: IF (T.EQ.1.0D0 .AND. FUN.EQ.21 .AND. WIO.NE.0) WRITE(WIO,12) RCOND ! 48: 12 FORMAT(0X,'WARNING.' ! 49: $ /0X,'EIGENVECTORS ARE BADLY CONDITIONED.' ! 50: $ /0X,'RESULTS MAY BE INACCURATE. RCOND =', 1PD13.4/) ! 51: DO 15 I = 1, M ! 52: DO 13 J = 1, N ! 53: LS = L+I-1+(J-1)*M ! 54: LL = L3+J-1 ! 55: STKR(LL) = STKR(LS) ! 56: STKI(LL) = -STKI(LS) ! 57: 13 CONTINUE ! 58: CALL WGESL(STKR(L2),STKI(L2),M2,N2,BUF,STKR(L3),STKI(L3),1) ! 59: DO 14 J = 1, N ! 60: LL = L+I-1+(J-1)*M ! 61: LS = L3+J-1 ! 62: STKR(LL) = STKR(LS) ! 63: STKI(LL) = -STKI(LS) ! 64: 14 CONTINUE ! 65: 15 CONTINUE ! 66: IF (FUN .NE. 21) GO TO 99 ! 67: C ! 68: C CHECK FOR IMAGINARY ROUNDOFF IN MATRIX FUNCTIONS ! 69: SR = WASUM(N*N,STKR(L),STKR(L),1) ! 70: SI = WASUM(N*N,STKI(L),STKI(L),1) ! 71: EPS = STKR(VSIZE-4) ! 72: T = EPS*SR ! 73: IF (DDT .EQ. 18) WRITE(WTE,115) SR,SI,EPS,T ! 74: 115 FORMAT(0X,'SR,SI,EPS,T',1P4D13.4) ! 75: IF (SI .LE. EPS*SR) CALL RSET(N*N,0.0D0,STKI(L),1) ! 76: GO TO 99 ! 77: C ! 78: 16 SR = STKR(L) ! 79: SI = STKI(L) ! 80: N = N2 ! 81: M = N ! 82: MSTK(TOP) = N ! 83: NSTK(TOP) = N ! 84: CALL WCOPY(N*N,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1) ! 85: GO TO 30 ! 86: C ! 87: C MATRIX LEFT DIVISION A BACKSLASH A2 ! 88: 20 L2 = LSTK(TOP+1) ! 89: M2 = MSTK(TOP+1) ! 90: N2 = NSTK(TOP+1) ! 91: IF (M .NE. N) CALL ERROR(20) ! 92: IF (ERR .GT. 0) RETURN ! 93: IF (M2*N2 .EQ. 1) GO TO 26 ! 94: L3 = L2 + M2*N2 ! 95: ERR = L3+N - LSTK(BOT) ! 96: IF (ERR .GT. 0) CALL ERROR(17) ! 97: IF (ERR .GT. 0) RETURN ! 98: CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3)) ! 99: IF (RCOND .EQ. 0.0D0) CALL ERROR(19) ! 100: IF (ERR .GT. 0) RETURN ! 101: T = FLOP(1.0D0 + RCOND) ! 102: IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND ! 103: IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND ! 104: IF (M2 .NE. N) CALL ERROR(12) ! 105: IF (ERR .GT. 0) RETURN ! 106: DO 23 J = 1, N2 ! 107: LJ = L2+(J-1)*M2 ! 108: CALL WGESL(STKR(L),STKI(L),M,N,BUF,STKR(LJ),STKI(LJ),0) ! 109: 23 CONTINUE ! 110: NSTK(TOP) = N2 ! 111: CALL WCOPY(M2*N2,STKR(L2),STKI(L2),1,STKR(L),STKI(L),1) ! 112: GO TO 99 ! 113: 26 SR = STKR(L2) ! 114: SI = STKI(L2) ! 115: GO TO 30 ! 116: C ! 117: C INV ! 118: C ! 119: 30 IF (M .NE. N) CALL ERROR(20) ! 120: IF (ERR .GT. 0) RETURN ! 121: IF (DDT .EQ. 17) GO TO 32 ! 122: DO 31 J = 1, N ! 123: DO 31 I = 1, N ! 124: LS = L+I-1+(J-1)*N ! 125: T0 = STKR(LS) ! 126: T1 = I+J-1 ! 127: T1 = FLOP(1.0D0/T1) ! 128: IF (T0 .NE. T1) GO TO 32 ! 129: 31 CONTINUE ! 130: GO TO 72 ! 131: 32 L3 = L + N*N ! 132: ERR = L3+N - LSTK(BOT) ! 133: IF (ERR .GT. 0) CALL ERROR(17) ! 134: IF (ERR .GT. 0) RETURN ! 135: CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3)) ! 136: IF (RCOND .EQ. 0.0D0) CALL ERROR(19) ! 137: IF (ERR .GT. 0) RETURN ! 138: T = FLOP(1.0D0 + RCOND) ! 139: IF (T .EQ. 1.0D0) WRITE(WTE,11) RCOND ! 140: IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE(WIO,11) RCOND ! 141: CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,STKR(L3),STKI(L3),1) ! 142: IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1) ! 143: GO TO 99 ! 144: C ! 145: C DET ! 146: C ! 147: 40 IF (M .NE. N) CALL ERROR(20) ! 148: IF (ERR .GT. 0) RETURN ! 149: CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO) ! 150: CALL WGEDI(STKR(L),STKI(L),M,N,BUF,DTR,DTI,SR,SI,10) ! 151: K = IDINT(DTR(2)) ! 152: KA = IABS(K)+2 ! 153: T = 1.0D0 ! 154: DO 41 I = 1, KA ! 155: T = T/10.0D0 ! 156: IF (T .EQ. 0.0D0) GO TO 42 ! 157: 41 CONTINUE ! 158: STKR(L) = DTR(1)*10.D0**K ! 159: STKI(L) = DTI(1)*10.D0**K ! 160: MSTK(TOP) = 1 ! 161: NSTK(TOP) = 1 ! 162: GO TO 99 ! 163: 42 IF (DTI(1) .EQ. 0.0D0) WRITE(WTE,43) DTR(1),K ! 164: IF (DTI(1) .NE. 0.0D0) WRITE(WTE,44) DTR(1),DTI(1),K ! 165: 43 FORMAT(0X,'DET = ',F7.4,7H * 10**,I4) ! 166: 44 FORMAT(0X,'DET = ',F7.4,' + ',F7.4,' i ',7H * 10**,I4) ! 167: STKR(L) = DTR(1) ! 168: STKI(L) = DTI(1) ! 169: STKR(L+1) = DTR(2) ! 170: STKI(L+1) = 0.0D0 ! 171: MSTK(TOP) = 1 ! 172: NSTK(TOP) = 2 ! 173: GO TO 99 ! 174: C ! 175: C RCOND ! 176: C ! 177: 50 IF (M .NE. N) CALL ERROR(20) ! 178: IF (ERR .GT. 0) RETURN ! 179: L3 = L + N*N ! 180: ERR = L3+N - LSTK(BOT) ! 181: IF (ERR .GT. 0) CALL ERROR(17) ! 182: IF (ERR .GT. 0) RETURN ! 183: CALL WGECO(STKR(L),STKI(L),M,N,BUF,RCOND,STKR(L3),STKI(L3)) ! 184: STKR(L) = RCOND ! 185: STKI(L) = 0.0D0 ! 186: MSTK(TOP) = 1 ! 187: NSTK(TOP) = 1 ! 188: IF (LHS .EQ. 1) GO TO 99 ! 189: L = L + 1 ! 190: CALL WCOPY(N,STKR(L3),STKI(L3),1,STKR(L),STKI(L),1) ! 191: TOP = TOP + 1 ! 192: LSTK(TOP) = L ! 193: MSTK(TOP) = N ! 194: NSTK(TOP) = 1 ! 195: GO TO 99 ! 196: C ! 197: C LU ! 198: C ! 199: 60 IF (M .NE. N) CALL ERROR(20) ! 200: IF (ERR .GT. 0) RETURN ! 201: CALL WGEFA(STKR(L),STKI(L),M,N,BUF,INFO) ! 202: IF (LHS .NE. 2) GO TO 99 ! 203: NN = N*N ! 204: IF (TOP+1 .GE. BOT) CALL ERROR(18) ! 205: IF (ERR .GT. 0) RETURN ! 206: TOP = TOP+1 ! 207: LSTK(TOP) = L + NN ! 208: MSTK(TOP) = N ! 209: NSTK(TOP) = N ! 210: ERR = L+NN+NN - LSTK(BOT) ! 211: IF (ERR .GT. 0) CALL ERROR(17) ! 212: IF (ERR .GT. 0) RETURN ! 213: DO 64 KB = 1, N ! 214: K = N+1-KB ! 215: DO 61 I = 1, N ! 216: LL = L+I-1+(K-1)*N ! 217: LU = LL + NN ! 218: IF (I .LE. K) STKR(LU) = STKR(LL) ! 219: IF (I .LE. K) STKI(LU) = STKI(LL) ! 220: IF (I .GT. K) STKR(LU) = 0.0D0 ! 221: IF (I .GT. K) STKI(LU) = 0.0D0 ! 222: IF (I .LT. K) STKR(LL) = 0.0D0 ! 223: IF (I .LT. K) STKI(LL) = 0.0D0 ! 224: IF (I .EQ. K) STKR(LL) = 1.0D0 ! 225: IF (I .EQ. K) STKI(LL) = 0.0D0 ! 226: IF (I .GT. K) STKR(LL) = -STKR(LL) ! 227: IF (I .GT. K) STKI(LL) = -STKI(LL) ! 228: 61 CONTINUE ! 229: I = BUF(K) ! 230: IF (I .EQ. K) GO TO 64 ! 231: LI = L+I-1+(K-1)*N ! 232: LK = L+K-1+(K-1)*N ! 233: CALL WSWAP(N-K+1,STKR(LI),STKI(LI),N,STKR(LK),STKI(LK),N) ! 234: 64 CONTINUE ! 235: GO TO 99 ! 236: C ! 237: C HILBERT ! 238: 70 N = IDINT(STKR(L)) ! 239: MSTK(TOP) = N ! 240: NSTK(TOP) = N ! 241: 72 CALL HILBER(STKR(L),N,N) ! 242: CALL RSET(N*N,0.0D0,STKI(L),1) ! 243: IF (FIN .LT. 0) CALL WSCAL(N*N,SR,SI,STKR(L),STKI(L),1) ! 244: GO TO 99 ! 245: C ! 246: C CHOLESKY ! 247: 80 IF (M .NE. N) CALL ERROR(20) ! 248: IF (ERR .GT. 0) RETURN ! 249: CALL WPOFA(STKR(L),STKI(L),M,N,ERR) ! 250: IF (ERR .NE. 0) CALL ERROR(29) ! 251: IF (ERR .GT. 0) RETURN ! 252: DO 81 J = 1, N ! 253: LL = L+J+(J-1)*M ! 254: CALL WSET(M-J,0.0D0,0.0D0,STKR(LL),STKI(LL),1) ! 255: 81 CONTINUE ! 256: GO TO 99 ! 257: C ! 258: C RREF ! 259: 85 IF (RHS .LT. 2) GO TO 86 ! 260: TOP = TOP-1 ! 261: L = LSTK(TOP) ! 262: IF (MSTK(TOP) .NE. M) CALL ERROR(5) ! 263: IF (ERR .GT. 0) RETURN ! 264: N = N + NSTK(TOP) ! 265: 86 CALL RREF(STKR(L),STKI(L),M,M,N,STKR(VSIZE-4)) ! 266: NSTK(TOP) = N ! 267: GO TO 99 ! 268: C ! 269: 99 RETURN ! 270: END
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.