home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / libs / matlab.lzh / MATLAB / MATLAB.LZH / Source / MatLab / matfn6.for < prev    next >
Encoding:
Text File  |  1991-04-13  |  7.5 KB  |  304 lines
  1.       SUBROUTINE MATFN6
  2.       IMPLICIT NONE
  3. C
  4. C EVALUATE UTILITY FUNCTIONS
  5. C
  6.       INCLUDE MATLAB$KOM:SIZEPARMS.INC
  7.       INCLUDE MATLAB$KOM:VSTK.KOM
  8.       INCLUDE MATLAB$KOM:IOP.KOM
  9.       INCLUDE MATLAB$KOM:COM.KOM
  10. C
  11.       INTEGER I, IA, IB, J, JA, JB, K, M, MA, MN, N, NA, NN
  12.       INTEGER L, LA, LB, LD, LJ, LL, LS
  13.       INTEGER SEMI, ID(4), UNIFOR(4), NORMAL(4), SEED(4)
  14.       DOUBLE PRECISION EPS0, EPS, S, SR, SI, T
  15. C
  16.       DOUBLE PRECISION FLOP, URAND
  17.       LOGICAL EQID
  18. C
  19.       DATA SEMI / 39 /, UNIFOR / 30, 23, 18, 15 /
  20.       DATA NORMAL / 23, 24, 27, 22 /, SEED / 28, 14, 14, 13 /
  21. C
  22. C
  23.       IF (DDT.EQ.1) WRITE (WTE, 100) FIN
  24. 100   FORMAT (' MATFN6', I4)
  25. C
  26. C FUNCTIONS/FIN
  27. C   MAGI DIAG SUM  PROD USER EYE  RAND ONES CHOP SIZE KRON  TRIL TRIU
  28. C     1    2    3    4    5    6    7    8    9   10  11-13  14   15
  29. C
  30.       L = LSTK(TOP)
  31.       M = MSTK(TOP)
  32.       N = NSTK(TOP)
  33.       GO TO (75, 80, 65, 67, 70, 90, 90, 90, 60, 77,
  34.      .       50, 50, 50, 80, 80), FIN
  35. C
  36. C ***      KRONECKER PRODUCT
  37. 50    CONTINUE
  38.       IF (RHS.NE.2) CALL ERROR (39)
  39.       IF (ERR.GT.0) RETURN
  40.       TOP = TOP-1
  41.       L = LSTK(TOP)
  42.       MA = MSTK(TOP)
  43.       NA = NSTK(TOP)
  44.       LA = L+MAX0 (M*N*MA*NA, M*N+MA*NA)
  45.       LB = LA+MA*NA
  46.       ERR = LB+M*N-LSTK(BOT)
  47.       IF (ERR.GT.0) CALL ERROR (17)
  48.       IF (ERR.GT.0) RETURN
  49. C ***      MOVE A AND B ABOVE RESULT
  50.       CALL WCOPY (MA*NA+M*N, STKR(L), STKI(L),
  51.      .            1, STKR(LA), STKI(LA), 1)
  52.       DO 54 JA = 1, NA
  53.         DO 53 J = 1, N
  54.           LJ = LB+(J-1)*M
  55.           DO 52 IA = 1, MA
  56. C ***            GET J-TH COLUMN OF B
  57.             CALL WCOPY (M, STKR(LJ), STKI(LJ), 1, STKR(L), STKI(L), 1)
  58. C ***            ADDRESS OF A(IA,JA)
  59.             LS = LA+IA-1+(JA-1)*MA
  60.             DO 51 I = 1, M
  61. C ***            A(IA,JA) OP B(I,J)
  62.               IF (FIN.EQ.11) CALL WMUL (STKR(LS), STKI(LS), STKR(L),
  63.      .                                  STKI(L), STKR(L), STKI(L))
  64.               IF (FIN.EQ.12) CALL WDIV (STKR(LS), STKI(LS), STKR(L),
  65.      .                                  STKI(L), STKR(L), STKI(L))
  66.               IF (FIN.EQ.13) CALL WDIV (STKR(L), STKI(L), STKR(LS),
  67.      .                                  STKI(LS), STKR(L), STKI(L))
  68.               IF (ERR.GT.0) RETURN
  69.               L = L+1
  70. 51          CONTINUE
  71. 52        CONTINUE
  72. 53      CONTINUE
  73. 54    CONTINUE
  74.       MSTK(TOP) = M*MA
  75.       NSTK(TOP) = N*NA
  76.       GO TO 99
  77. C
  78. C ***      CHOP
  79. 60    CONTINUE
  80.       EPS0 = 1.0D0
  81. 61    CONTINUE
  82.       EPS0 = EPS0/2.0D0
  83.       T = FLOP (1.0D0+EPS0)
  84.       IF (T.GT.1.0D0) GO TO 61
  85.       EPS0 = 2.0D0*EPS0
  86.       FLP(2) = IDINT (STKR(L))
  87.       IF (SYM.NE.SEMI) WRITE (WTE, 62) FLP(2)
  88. 62    FORMAT (/, ' CHOP ', I2, ' PLACES.')
  89.       EPS = 1.0D0
  90. 63    CONTINUE
  91.       EPS = EPS/2.0D0
  92.       T = FLOP (1.0D0+EPS)
  93.       IF (T.GT.1.0D0) GO TO 63
  94.       EPS = 2.0D0*EPS
  95.       T = STKR(VSIZE-4)
  96.       IF (T.LT.EPS .OR. T.EQ.EPS0) STKR(VSIZE-4) = EPS
  97.       MSTK(TOP) = 0
  98.       GO TO 99
  99. C
  100. C ***      SUM
  101. 65    CONTINUE
  102.       SR = 0.0D0
  103.       SI = 0.0D0
  104.       MN = M*N
  105.       DO 66 I = 1, MN
  106.         LS = L+I-1
  107.         SR = FLOP (SR+STKR(LS))
  108.         SI = FLOP (SI+STKI(LS))
  109. 66    CONTINUE
  110.       GO TO 69
  111. C
  112. C ***      PROD
  113. 67    CONTINUE
  114.       SR = 1.0D0
  115.       SI = 0.0D0
  116.       MN = M*N
  117.       DO 68 I = 1, MN
  118.         LS = L+I-1
  119.         CALL WMUL (STKR(LS), STKI(LS), SR, SI, SR, SI)
  120. 68    CONTINUE
  121. 69    CONTINUE
  122.       STKR(L) = SR
  123.       STKI(L) = SI
  124.       MSTK(TOP) = 1
  125.       NSTK(TOP) = 1
  126.       GO TO 99
  127. C
  128. C ***      USER
  129. 70    CONTINUE
  130.       S = 0.0D0
  131.       T = 0.0D0
  132.       IF (RHS.LT.2) GO TO 72
  133.       IF (RHS.LT.3) GO TO 71
  134.       T = STKR(L)
  135.       TOP = TOP-1
  136.       L = LSTK(TOP)
  137.       M = MSTK(TOP)
  138.       N = NSTK(TOP)
  139. 71    CONTINUE
  140.       S = STKR(L)
  141.       TOP = TOP-1
  142.       L = LSTK(TOP)
  143.       M = MSTK(TOP)
  144.       N = NSTK(TOP)
  145. 72    CONTINUE
  146.       CALL USER (STKR(L), M, N, S, T)
  147.       CALL RSET (M*N, 0.0D0, STKI(L), 1)
  148.       MSTK(TOP) = M
  149.       NSTK(TOP) = N
  150.       GO TO 99
  151. C
  152. C ***      MAGIC
  153. 75    CONTINUE
  154.       N = MAX0 (IDINT (STKR(L)), 0)
  155.       IF (N.EQ.2) N = 0
  156.       IF (N.GT.0) CALL MAGIC (STKR(L), N, N)
  157.       CALL RSET (N*N, 0.0D0, STKI(L), 1)
  158.       MSTK(TOP) = N
  159.       NSTK(TOP) = N
  160.       GO TO 99
  161. C
  162. C ***      SIZE
  163. 77    CONTINUE
  164.       STKR(L) = M
  165.       STKR(L+1) = N
  166.       STKI(L) = 0.0D0
  167.       STKI(L+1) = 0.0D0
  168.       MSTK(TOP) = 1
  169.       NSTK(TOP) = 2
  170.       IF (LHS.EQ.1) GO TO 99
  171.       NSTK(TOP) = 1
  172.       TOP = TOP+1
  173.       LSTK(TOP) = L+1
  174.       MSTK(TOP) = 1
  175.       NSTK(TOP) = 1
  176.       GO TO 99
  177. C
  178. C ***      DIAG, TRIU, TRIL
  179. 80    CONTINUE
  180.       K = 0
  181.       IF (RHS.NE.2) GO TO 81
  182.       K = IDINT (STKR(L))
  183.       TOP = TOP-1
  184.       L = LSTK(TOP)
  185.       M = MSTK(TOP)
  186.       N = NSTK(TOP)
  187. 81    CONTINUE
  188.       IF (FIN.GE.14) GO TO 85
  189.       IF (M.EQ.1 .OR. N.EQ.1) GO TO 83
  190.       IF (K.GE.0) MN = MIN0 (M, N-K)
  191.       IF (K.LT.0) MN = MIN0 (M+K, N)
  192.       MSTK(TOP) = MAX0 (MN, 0)
  193.       NSTK(TOP) = 1
  194.       IF (MN.LE.0) GO TO 99
  195.       DO 82 I = 1, MN
  196.         IF (K.GE.0) LS = L+(I-1)+(I+K-1)*M
  197.         IF (K.LT.0) LS = L+(I-K-1)+(I-1)*M
  198.         LL = L+I-1
  199.         STKR(LL) = STKR(LS)
  200.         STKI(LL) = STKI(LS)
  201. 82    CONTINUE
  202.       GO TO 99
  203. C
  204. 83    CONTINUE
  205.       N = MAX0 (M, N)+IABS (K)
  206.       ERR = L+N*N-LSTK(BOT)
  207.       IF (ERR.GT.0) CALL ERROR (17)
  208.       IF (ERR.GT.0) RETURN
  209.       MSTK(TOP) = N
  210.       NSTK(TOP) = N
  211.       DO 84 JB = 1, N
  212.       DO 84 IB = 1, N
  213.         J = N+1-JB
  214.         I = N+1-IB
  215.         SR = 0.0D0
  216.         SI = 0.0D0
  217.         IF (K.GE.0) LS = L+I-1
  218.         IF (K.LT.0) LS = L+J-1
  219.         LL = L+I-1+(J-1)*N
  220.         IF (J-I.EQ.K) SR = STKR(LS)
  221.         IF (J-I.EQ.K) SI = STKI(LS)
  222.         STKR(LL) = SR
  223.         STKI(LL) = SI
  224. 84    CONTINUE
  225.       GO TO 99
  226. C
  227. C ***      TRIL, TRIU
  228. 85    CONTINUE
  229.       DO 87 J = 1, N
  230.         LD = L+J-K-1+(J-1)*M
  231.         IF (FIN.EQ.14) LL = J-K-1
  232.         IF (FIN.EQ.14) LS = LD-LL
  233.         IF (FIN.EQ.15) LL = M-J+K
  234.         IF (FIN.EQ.15) LS = LD+1
  235.         IF (LL.GT.0) CALL WSET (LL, 0.0D0, 0.0D0, STKR(LS),
  236.      .                          STKI(LS), 1)
  237. 87    CONTINUE
  238.       GO TO 99
  239. C
  240. C ***      EYE, RAND, ONES
  241. 90    CONTINUE
  242.       IF (M.GT.1 .OR. RHS.EQ.0) GO TO 94
  243.       IF (RHS.NE.2) GO TO 91
  244.       NN = IDINT (STKR(L))
  245.       TOP = TOP-1
  246.       L = LSTK(TOP)
  247.       N = NSTK(TOP)
  248. 91    CONTINUE
  249.       IF (FIN.NE.7 .OR. N.LT.4) GO TO 93
  250.       DO 92 I = 1, 4
  251.         LS = L+I-1
  252.         ID(I) = IDINT (STKR(LS))
  253. 92    CONTINUE
  254.       IF (EQID (ID, UNIFOR) .OR. EQID (ID, NORMAL)) GO TO 97
  255.       IF (EQID (ID, SEED)) GO TO 98
  256. 93    CONTINUE
  257.       IF (N.GT.1) GO TO 94
  258.       M = MAX0 (IDINT (STKR(L)), 0)
  259.       IF (RHS.EQ.2) N = MAX0 (NN, 0)
  260.       IF (RHS.NE.2) N = M
  261.       ERR = L+M*N-LSTK(BOT)
  262.       IF (ERR.GT.0) CALL ERROR (17)
  263.       IF (ERR.GT.0) RETURN
  264.       MSTK(TOP) = M
  265.       NSTK(TOP) = N
  266.       IF (M*N.EQ.0) GO TO 99
  267. 94    CONTINUE
  268.       DO 96 J = 1, N
  269.       DO 96 I = 1, M
  270.         LL = L+I-1+(J-1)*M
  271.         STKR(LL) = 0.0D0
  272.         STKI(LL) = 0.0D0
  273.         IF (I.EQ.J .OR. FIN.EQ.8) STKR(LL) = 1.0D0
  274.         IF (FIN.EQ.7 .AND. RAN(2).EQ.0)
  275.      .                         STKR(LL) = FLOP (URAND (RAN(1)))
  276.         IF (FIN.NE.7 .OR. RAN(2).EQ.0) GO TO 96
  277. 95      CONTINUE
  278.         SR = 2.0D0*URAND (RAN(1))-1.0D0
  279.         SI = 2.0D0*URAND (RAN(1))-1.0D0
  280.         T = SR*SR+SI*SI
  281.         IF (T.GT.1.0D0) GO TO 95
  282.         STKR(LL) = FLOP (SR*DSQRT (-2.0D0*DLOG (T)/T))
  283. 96    CONTINUE
  284.       GO TO 99
  285. C
  286. C ***      SWITCH UNIFORM AND NORMAL
  287. 97    CONTINUE
  288.       RAN(2) = ID(1)-UNIFOR(1)
  289.       MSTK(TOP) = 0
  290.       GO TO 99
  291. C
  292. C ***      SEED
  293. 98    CONTINUE
  294.       IF (RHS.EQ.2) RAN(1) = NN
  295.       STKR(L) = RAN(1)
  296.       MSTK(TOP) = 1
  297.       IF (RHS.EQ.2) MSTK(TOP) = 0
  298.       NSTK(TOP) = 1
  299.       GO TO 99
  300. C
  301. 99    CONTINUE
  302.       RETURN
  303.       END
  304.