![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to matfn6.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 matfn6.f CVS log | Up to [Research Unix] / researchv10dc / cmd / matlab / src |
1.1 ! root 1: SUBROUTINE MATFN6 ! 2: C ! 3: C EVALUATE UTILITY FUNCTIONS ! 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: INTEGER SEMI,ID(4),UNIFOR(4),NORMAL(4),SEED(4) ! 13: DOUBLE PRECISION EPS0,EPS,S,SR,SI,T ! 14: DOUBLE PRECISION FLOP,URAND ! 15: LOGICAL EQID ! 16: DATA SEMI/39/ ! 17: DATA UNIFOR/30,23,18,15/,NORMAL/23,24,27,22/,SEED/28,14,14,13/ ! 18: C ! 19: IF (DDT .EQ. 1) WRITE(WTE,100) FIN ! 20: 100 FORMAT(0X,'MATFN6',I4) ! 21: C FUNCTIONS/FIN ! 22: C MAGI DIAG SUM PROD USER EYE RAND ONES CHOP SIZE KRON TRIL TRIU ! 23: C 1 2 3 4 5 6 7 8 9 10 11-13 14 15 ! 24: L = LSTK(TOP) ! 25: M = MSTK(TOP) ! 26: N = NSTK(TOP) ! 27: GO TO (75,80,65,67,70,90,90,90,60,77,50,50,50,80,80),FIN ! 28: C ! 29: C KRONECKER PRODUCT ! 30: 50 IF (RHS .NE. 2) CALL ERROR(39) ! 31: IF (ERR .GT. 0) RETURN ! 32: TOP = TOP - 1 ! 33: L = LSTK(TOP) ! 34: MA = MSTK(TOP) ! 35: NA = NSTK(TOP) ! 36: LA = L + MAX0(M*N*MA*NA,M*N+MA*NA) ! 37: LB = LA + MA*NA ! 38: ERR = LB + M*N - LSTK(BOT) ! 39: IF (ERR .GT. 0) CALL ERROR(17) ! 40: IF (ERR .GT. 0) RETURN ! 41: C MOVE A AND B ABOVE RESULT ! 42: CALL WCOPY(MA*NA+M*N,STKR(L),STKI(L),1,STKR(LA),STKI(LA),1) ! 43: DO 54 JA = 1, NA ! 44: DO 53 J = 1, N ! 45: LJ = LB + (J-1)*M ! 46: DO 52 IA = 1, MA ! 47: C GET J-TH COLUMN OF B ! 48: CALL WCOPY(M,STKR(LJ),STKI(LJ),1,STKR(L),STKI(L),1) ! 49: C ADDRESS OF A(IA,JA) ! 50: LS = LA + IA-1 + (JA-1)*MA ! 51: DO 51 I = 1, M ! 52: C A(IA,JA) OP B(I,J) ! 53: IF (FIN .EQ. 11) CALL WMUL(STKR(LS),STKI(LS), ! 54: $ STKR(L),STKI(L),STKR(L),STKI(L)) ! 55: IF (FIN .EQ. 12) CALL WDIV(STKR(LS),STKI(LS), ! 56: $ STKR(L),STKI(L),STKR(L),STKI(L)) ! 57: IF (FIN .EQ. 13) CALL WDIV(STKR(L),STKI(L), ! 58: $ STKR(LS),STKI(LS),STKR(L),STKI(L)) ! 59: IF (ERR .GT. 0) RETURN ! 60: L = L + 1 ! 61: 51 CONTINUE ! 62: 52 CONTINUE ! 63: 53 CONTINUE ! 64: 54 CONTINUE ! 65: MSTK(TOP) = M*MA ! 66: NSTK(TOP) = N*NA ! 67: GO TO 99 ! 68: C ! 69: C CHOP ! 70: 60 EPS0 = 1.0D0 ! 71: 61 EPS0 = EPS0/2.0D0 ! 72: T = FLOP(1.0D0 + EPS0) ! 73: IF (T .GT. 1.0D0) GO TO 61 ! 74: EPS0 = 2.0D0*EPS0 ! 75: FLP(2) = IDINT(STKR(L)) ! 76: IF (SYM .NE. SEMI) WRITE(WTE,62) FLP(2) ! 77: 62 FORMAT(/0X,'CHOP ',I2,' PLACES.') ! 78: EPS = 1.0D0 ! 79: 63 EPS = EPS/2.0D0 ! 80: T = FLOP(1.0D0 + EPS) ! 81: IF (T .GT. 1.0D0) GO TO 63 ! 82: EPS = 2.0D0*EPS ! 83: T = STKR(VSIZE-4) ! 84: IF (T.LT.EPS .OR. T.EQ.EPS0) STKR(VSIZE-4) = EPS ! 85: MSTK(TOP) = 0 ! 86: GO TO 99 ! 87: C ! 88: C SUM ! 89: 65 SR = 0.0D0 ! 90: SI = 0.0D0 ! 91: MN = M*N ! 92: DO 66 I = 1, MN ! 93: LS = L+I-1 ! 94: SR = FLOP(SR+STKR(LS)) ! 95: SI = FLOP(SI+STKI(LS)) ! 96: 66 CONTINUE ! 97: GO TO 69 ! 98: C ! 99: C PROD ! 100: 67 SR = 1.0D0 ! 101: SI = 0.0D0 ! 102: MN = M*N ! 103: DO 68 I = 1, MN ! 104: LS = L+I-1 ! 105: CALL WMUL(STKR(LS),STKI(LS),SR,SI,SR,SI) ! 106: 68 CONTINUE ! 107: 69 STKR(L) = SR ! 108: STKI(L) = SI ! 109: MSTK(TOP) = 1 ! 110: NSTK(TOP) = 1 ! 111: GO TO 99 ! 112: C ! 113: C USER ! 114: 70 S = 0.0D0 ! 115: T = 0.0D0 ! 116: IF (RHS .LT. 2) GO TO 72 ! 117: IF (RHS .LT. 3) GO TO 71 ! 118: T = STKR(L) ! 119: TOP = TOP-1 ! 120: L = LSTK(TOP) ! 121: M = MSTK(TOP) ! 122: N = NSTK(TOP) ! 123: 71 S = STKR(L) ! 124: TOP = TOP-1 ! 125: L = LSTK(TOP) ! 126: M = MSTK(TOP) ! 127: N = NSTK(TOP) ! 128: 72 CALL USER(STKR(L),M,N,S,T) ! 129: CALL RSET(M*N,0.0D0,STKI(L),1) ! 130: MSTK(TOP) = M ! 131: NSTK(TOP) = N ! 132: GO TO 99 ! 133: C ! 134: C MAGIC ! 135: 75 N = MAX0(IDINT(STKR(L)),0) ! 136: IF (N .EQ. 2) N = 0 ! 137: IF (N .GT. 0) CALL MAGIC(STKR(L),N,N) ! 138: CALL RSET(N*N,0.0D0,STKI(L),1) ! 139: MSTK(TOP) = N ! 140: NSTK(TOP) = N ! 141: GO TO 99 ! 142: C ! 143: C SIZE ! 144: 77 STKR(L) = M ! 145: STKR(L+1) = N ! 146: STKI(L) = 0.0D0 ! 147: STKI(L+1) = 0.0D0 ! 148: MSTK(TOP) = 1 ! 149: NSTK(TOP) = 2 ! 150: IF (LHS .EQ. 1) GO TO 99 ! 151: NSTK(TOP) = 1 ! 152: TOP = TOP + 1 ! 153: LSTK(TOP) = L+1 ! 154: MSTK(TOP) = 1 ! 155: NSTK(TOP) = 1 ! 156: GO TO 99 ! 157: C ! 158: C DIAG, TRIU, TRIL ! 159: 80 K = 0 ! 160: IF (RHS .NE. 2) GO TO 81 ! 161: K = IDINT(STKR(L)) ! 162: TOP = TOP-1 ! 163: L = LSTK(TOP) ! 164: M = MSTK(TOP) ! 165: N = NSTK(TOP) ! 166: 81 IF (FIN .GE. 14) GO TO 85 ! 167: IF (M .EQ. 1 .OR. N .EQ. 1) GO TO 83 ! 168: IF (K.GE.0) MN=MIN0(M,N-K) ! 169: IF (K.LT.0) MN=MIN0(M+K,N) ! 170: MSTK(TOP) = MAX0(MN,0) ! 171: NSTK(TOP) = 1 ! 172: IF (MN .LE. 0) GO TO 99 ! 173: DO 82 I = 1, MN ! 174: IF (K.GE.0) LS = L+(I-1)+(I+K-1)*M ! 175: IF (K.LT.0) LS = L+(I-K-1)+(I-1)*M ! 176: LL = L+I-1 ! 177: STKR(LL) = STKR(LS) ! 178: STKI(LL) = STKI(LS) ! 179: 82 CONTINUE ! 180: GO TO 99 ! 181: 83 N = MAX0(M,N)+IABS(K) ! 182: ERR = L+N*N - LSTK(BOT) ! 183: IF (ERR .GT. 0) CALL ERROR(17) ! 184: IF (ERR .GT. 0) RETURN ! 185: MSTK(TOP) = N ! 186: NSTK(TOP) = N ! 187: DO 84 JB = 1, N ! 188: DO 84 IB = 1, N ! 189: J = N+1-JB ! 190: I = N+1-IB ! 191: SR = 0.0D0 ! 192: SI = 0.0D0 ! 193: IF (K.GE.0) LS = L+I-1 ! 194: IF (K.LT.0) LS = L+J-1 ! 195: LL = L+I-1+(J-1)*N ! 196: IF (J-I .EQ. K) SR = STKR(LS) ! 197: IF (J-I .EQ. K) SI = STKI(LS) ! 198: STKR(LL) = SR ! 199: STKI(LL) = SI ! 200: 84 CONTINUE ! 201: GO TO 99 ! 202: C ! 203: C TRIL, TRIU ! 204: 85 DO 87 J = 1, N ! 205: LD = L + J - K - 1 + (J-1)*M ! 206: IF (FIN .EQ. 14) LL = J - K - 1 ! 207: IF (FIN .EQ. 14) LS = LD - LL ! 208: IF (FIN .EQ. 15) LL = M - J + K ! 209: IF (FIN .EQ. 15) LS = LD + 1 ! 210: IF (LL .GT. 0) CALL WSET(LL,0.0D0,0.0D0,STKR(LS),STKI(LS),1) ! 211: 87 CONTINUE ! 212: GO TO 99 ! 213: C ! 214: C EYE, RAND, ONES ! 215: 90 IF (M.GT.1 .OR. RHS.EQ.0) GO TO 94 ! 216: IF (RHS .NE. 2) GO TO 91 ! 217: NN = IDINT(STKR(L)) ! 218: TOP = TOP-1 ! 219: L = LSTK(TOP) ! 220: N = NSTK(TOP) ! 221: 91 IF (FIN.NE.7 .OR. N.LT.4) GO TO 93 ! 222: DO 92 I = 1, 4 ! 223: LS = L+I-1 ! 224: ID(I) = IDINT(STKR(LS)) ! 225: 92 CONTINUE ! 226: IF (EQID(ID,UNIFOR).OR.EQID(ID,NORMAL)) GO TO 97 ! 227: IF (EQID(ID,SEED)) GO TO 98 ! 228: 93 IF (N .GT. 1) GO TO 94 ! 229: M = MAX0(IDINT(STKR(L)),0) ! 230: IF (RHS .EQ. 2) N = MAX0(NN,0) ! 231: IF (RHS .NE. 2) N = M ! 232: ERR = L+M*N - LSTK(BOT) ! 233: IF (ERR .GT. 0) CALL ERROR(17) ! 234: IF (ERR .GT. 0) RETURN ! 235: MSTK(TOP) = M ! 236: NSTK(TOP) = N ! 237: IF (M*N .EQ. 0) GO TO 99 ! 238: 94 DO 96 J = 1, N ! 239: DO 96 I = 1, M ! 240: LL = L+I-1+(J-1)*M ! 241: STKR(LL) = 0.0D0 ! 242: STKI(LL) = 0.0D0 ! 243: IF (I.EQ.J .OR. FIN.EQ.8) STKR(LL) = 1.0D0 ! 244: IF (FIN.EQ.7 .AND. RAN(2).EQ.0) STKR(LL) = FLOP(URAND(RAN(1))) ! 245: IF (FIN.NE.7 .OR. RAN(2).EQ.0) GO TO 96 ! 246: 95 SR = 2.0D0*URAND(RAN(1))-1.0D0 ! 247: SI = 2.0D0*URAND(RAN(1))-1.0D0 ! 248: T = SR*SR + SI*SI ! 249: IF (T .GT. 1.0D0) GO TO 95 ! 250: STKR(LL) = FLOP(SR*DSQRT(-2.0D0*DLOG(T)/T)) ! 251: 96 CONTINUE ! 252: GO TO 99 ! 253: C ! 254: C SWITCH UNIFORM AND NORMAL ! 255: 97 RAN(2) = ID(1) - UNIFOR(1) ! 256: MSTK(TOP) = 0 ! 257: GO TO 99 ! 258: C ! 259: C SEED ! 260: 98 IF (RHS .EQ. 2) RAN(1) = NN ! 261: STKR(L) = RAN(1) ! 262: MSTK(TOP) = 1 ! 263: IF (RHS .EQ. 2) MSTK(TOP) = 0 ! 264: NSTK(TOP) = 1 ! 265: GO TO 99 ! 266: C ! 267: 99 RETURN ! 268: END
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.