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 / matfn1.for < prev    next >
Encoding:
Text File  |  1991-04-13  |  8.2 KB  |  298 lines
  1.       SUBROUTINE MATFN1
  2.       IMPLICIT NONE
  3. C
  4. C EVALUATE FUNCTIONS INVOLVING GAUSSIAN ELIMINATION
  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, J, K, KA, KB, L, L2, L3, LS, LL, LJ, LU, LI, LK
  12.       INTEGER M, M2, N, N2, NN, INFO
  13.       DOUBLE PRECISION DTR(2), DTI(2), SR, SI, RCOND, T, T0, T1, EPS
  14. C
  15.       DOUBLE PRECISION FLOP, WASUM, DFLOAT
  16. C
  17. C
  18.       IF (DDT.EQ.1) WRITE (WTE, 100) FIN
  19. 100   FORMAT (' MATFN1', I4)
  20. C
  21.       L = LSTK(TOP)
  22.       M = MSTK(TOP)
  23.       N = NSTK(TOP)
  24.       IF (FIN.EQ.-1) GO TO 10
  25.       IF (FIN.EQ.-2) GO TO 20
  26.       GO TO (30, 40, 50, 60, 70, 80, 85), FIN
  27. C
  28. C ***      MATRIX RIGHT DIVISION, A/A2
  29. 10    CONTINUE
  30.       L2 = LSTK(TOP+1)
  31.       M2 = MSTK(TOP+1)
  32.       N2 = NSTK(TOP+1)
  33.       IF (M2.NE.N2) CALL ERROR (20)
  34.       IF (ERR.GT.0) RETURN
  35.       IF (M*N.EQ.1) GO TO 16
  36.       IF (N.NE.N2) CALL ERROR (11)
  37.       IF (ERR.GT.0) RETURN
  38.       L3 = L2+M2*N2
  39.       ERR = L3+N2-LSTK(BOT)
  40.       IF (ERR.GT.0) CALL ERROR (17)
  41.       IF (ERR.GT.0) RETURN
  42.       CALL WGECO (STKR(L2), STKI(L2), M2, N2, BUF, RCOND,
  43.      .            STKR(L3), STKI(L3))
  44.       IF (RCOND.EQ.0.0D0) CALL ERROR (19)
  45.       IF (ERR.GT.0) RETURN
  46.       T = FLOP (1.0D0+RCOND)
  47.       IF (T.EQ.1.0D0 .AND. FUN.NE.21) WRITE (WTE, 11) RCOND
  48.       IF (T.EQ.1.0D0 .AND. FUN.NE.21 .AND. WIO.NE.0)
  49.      .       WRITE (WIO, 11) RCOND
  50. 11    FORMAT (' WARNING.', /,
  51.      .        ' MATRIX IS CLOSE TO SINGULAR OR BADLY SCALED.', /,
  52.      .        ' RESULTS MAY BE INACCURATE.  RCOND = ', 1PD13.4, /)
  53.       IF (T.EQ.1.0D0 .AND. FUN.EQ.21) WRITE (WTE, 12) RCOND
  54.       IF (T.EQ.1.0D0 .AND. FUN.EQ.21 .AND. WIO.NE.0)
  55.      .       WRITE (WIO, 12) RCOND
  56. 12    FORMAT (' WARNING.', /,
  57.      .        ' EIGENVECTORS ARE BADLY CONDITIONED.', /,
  58.      .        ' RESULTS MAY BE INACCURATE.  RCOND = ', 1PD13.4, /)
  59.       DO 15 I = 1, M
  60.         DO 13 J = 1, N
  61.           LS = L+I-1+(J-1)*M
  62.           LL = L3+J-1
  63.           STKR(LL) = STKR(LS)
  64.           STKI(LL) = -STKI(LS)
  65. 13      CONTINUE
  66.         CALL WGESL (STKR(L2), STKI(L2), M2, N2, BUF,
  67.      .              STKR(L3), STKI(L3), 1)
  68.         DO 14 J = 1, N
  69.           LL = L+I-1+(J-1)*M
  70.           LS = L3+J-1
  71.           STKR(LL) = STKR(LS)
  72.           STKI(LL) = -STKI(LS)
  73. 14      CONTINUE
  74. 15    CONTINUE
  75.       IF (FUN.NE.21) GO TO 99
  76. C
  77. C ***      CHECK FOR IMAGINARY ROUNDOFF IN MATRIX FUNCTIONS
  78.       SR = WASUM (N*N, STKR(L), STKR(L), 1)
  79.       SI = WASUM (N*N, STKI(L), STKI(L), 1)
  80.       EPS = STKR(VSIZE-4)
  81.       T = EPS*SR
  82.       IF (DDT.EQ.18) WRITE (WTE, 115) SR, SI, EPS, T
  83. 115   FORMAT (' SR, SI, EPS, T', 1P4D13.4)
  84.       IF (SI.LE.EPS*SR) CALL RSET (N*N, 0.0D0, STKI(L), 1)
  85.       GO TO 99
  86. C
  87. 16    CONTINUE
  88.       SR = STKR(L)
  89.       SI = STKI(L)
  90.       N = N2
  91.       M = N
  92.       MSTK(TOP) = N
  93.       NSTK(TOP) = N
  94.       CALL WCOPY (N*N, STKR(L2), STKI(L2), 1, STKR(L), STKI(L), 1)
  95.       GO TO 30
  96. C
  97. C ***      MATRIX LEFT DIVISION A BACKSLASH A2
  98. 20    CONTINUE
  99.       L2 = LSTK(TOP+1)
  100.       M2 = MSTK(TOP+1)
  101.       N2 = NSTK(TOP+1)
  102.       IF (M.NE.N) CALL ERROR (20)
  103.       IF (ERR.GT.0) RETURN
  104.       IF (M2*N2.EQ.1) GO TO 26
  105.       L3 = L2+M2*N2
  106.       ERR = L3+N-LSTK(BOT)
  107.       IF (ERR.GT.0) CALL ERROR (17)
  108.       IF (ERR.GT.0) RETURN
  109.       CALL WGECO (STKR(L), STKI(L), M, N, BUF, RCOND,
  110.      .            STKR(L3), STKI(L3))
  111.       IF (RCOND.EQ.0.0D0) CALL ERROR (19)
  112.       IF (ERR.GT.0) RETURN
  113.       T = FLOP (1.0D0+RCOND)
  114.       IF (T.EQ.1.0D0) WRITE (WTE, 11) RCOND
  115.       IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE (WIO, 11) RCOND
  116.       IF (M2.NE.N) CALL ERROR (12)
  117.       IF (ERR.GT.0) RETURN
  118.       DO 23 J = 1, N2
  119.         LJ = L2+(J-1)*M2
  120.         CALL WGESL (STKR(L), STKI(L), M, N, BUF, STKR(LJ), STKI(LJ), 0)
  121. 23    CONTINUE
  122.       NSTK(TOP) = N2
  123.       CALL WCOPY (M2*N2, STKR(L2), STKI(L2), 1, STKR(L), STKI(L), 1)
  124.       GO TO 99
  125. C
  126. 26    CONTINUE
  127.       SR = STKR(L2)
  128.       SI = STKI(L2)
  129.       GO TO 30
  130. C
  131. C ***      INV
  132. 30    CONTINUE
  133.       IF (M.NE.N) CALL ERROR (20)
  134.       IF (ERR.GT.0) RETURN
  135.       IF (DDT.EQ.17) GO TO 32
  136.       DO 31 J = 1, N
  137.       DO 31 I = 1, N
  138.         LS = L+I-1+(J-1)*N
  139.         T0 = STKR(LS)
  140.         T1 = FLOP (1.0D0/(DFLOAT (I+J-1)))
  141.         IF (T0.NE.T1) GO TO 32
  142. 31    CONTINUE
  143.       GO TO 72
  144. C
  145. 32    CONTINUE
  146.       L3 = L+N*N
  147.       ERR = L3+N-LSTK(BOT)
  148.       IF (ERR.GT.0) CALL ERROR (17)
  149.       IF (ERR.GT.0) RETURN
  150.       CALL WGECO (STKR(L), STKI(L), M, N, BUF, RCOND,
  151.      .            STKR(L3), STKI(L3))
  152.       IF (RCOND.EQ.0.0D0) CALL ERROR (19)
  153.       IF (ERR.GT.0) RETURN
  154.       T = FLOP (1.0D0+RCOND)
  155.       IF (T.EQ.1.0D0) WRITE (WTE, 11) RCOND
  156.       IF (T.EQ.1.0D0 .AND. WIO.NE.0) WRITE (WIO, 11) RCOND
  157.       CALL WGEDI (STKR(L), STKI(L), M, N, BUF, DTR, DTI,
  158.      .            STKR(L3), STKI(L3), 1)
  159.       IF (FIN.LT.0) CALL WSCAL (N*N, SR, SI, STKR(L), STKI(L), 1)
  160.       GO TO 99
  161. C
  162. C ***      DET
  163. 40    CONTINUE
  164.       IF (M.NE.N) CALL ERROR (20)
  165.       IF (ERR.GT.0) RETURN
  166.       CALL WGEFA (STKR(L), STKI(L), M, N, BUF, INFO)
  167.       CALL WGEDI (STKR(L), STKI(L), M, N, BUF, DTR, DTI, SR, SI, 10)
  168.       K = IDINT (DTR(2))
  169.       KA = IABS (K)+2
  170.       T = 1.0D0
  171.       DO 41 I = 1, KA
  172.         T = T/10.0D0
  173.         IF (T.EQ.0.0D0) GO TO 42
  174. 41    CONTINUE
  175.       STKR(L) = DTR(1)*10.D0**K
  176.       STKI(L) = DTI(1)*10.D0**K
  177.       MSTK(TOP) = 1
  178.       NSTK(TOP) = 1
  179.       GO TO 99
  180. C
  181. 42    CONTINUE
  182.       IF (DTI(1).EQ.0.0D0) WRITE (WTE, 43) DTR(1), K
  183.       IF (DTI(1).NE.0.0D0) WRITE (WTE, 44) DTR(1), DTI(1), K
  184. 43    FORMAT (' DET = ', F7.4, ' * 10**', I4)
  185. 44    FORMAT (' DET = ', F7.4, ' + ', F7.4, ' i  * 10**', I4)
  186.       STKR(L) = DTR(1)
  187.       STKI(L) = DTI(1)
  188.       STKR(L+1) = DTR(2)
  189.       STKI(L+1) = 0.0D0
  190.       MSTK(TOP) = 1
  191.       NSTK(TOP) = 2
  192.       GO TO 99
  193. C
  194. C ***      RCOND
  195. 50    CONTINUE
  196.       IF (M.NE.N) CALL ERROR (20)
  197.       IF (ERR.GT.0) RETURN
  198.       L3 = L+N*N
  199.       ERR = L3+N-LSTK(BOT)
  200.       IF (ERR.GT.0) CALL ERROR (17)
  201.       IF (ERR.GT.0) RETURN
  202.       CALL WGECO (STKR(L), STKI(L), M, N, BUF, RCOND,
  203.      .            STKR(L3), STKI(L3))
  204.       STKR(L) = RCOND
  205.       STKI(L) = 0.0D0
  206.       MSTK(TOP) = 1
  207.       NSTK(TOP) = 1
  208.       IF (LHS.EQ.1) GO TO 99
  209.       L = L+1
  210.       CALL WCOPY (N, STKR(L3), STKI(L3), 1, STKR(L), STKI(L), 1)
  211.       TOP = TOP+1
  212.       LSTK(TOP) = L
  213.       MSTK(TOP) = N
  214.       NSTK(TOP) = 1
  215.       GO TO 99
  216. C
  217. C ***      LU
  218. 60    CONTINUE
  219.       IF (M.NE.N) CALL ERROR (20)
  220.       IF (ERR.GT.0) RETURN
  221.       CALL WGEFA (STKR(L), STKI(L), M, N, BUF, INFO)
  222.       IF (LHS.NE.2) GO TO 99
  223.       NN = N*N
  224.       IF (TOP+1.GE.BOT) CALL ERROR (18)
  225.       IF (ERR.GT.0) RETURN
  226.       TOP = TOP+1
  227.       LSTK(TOP) = L+NN
  228.       MSTK(TOP) = N
  229.       NSTK(TOP) = N
  230.       ERR = L+NN+NN-LSTK(BOT)
  231.       IF (ERR.GT.0) CALL ERROR (17)
  232.       IF (ERR.GT.0) RETURN
  233.       DO 64 KB = 1, N
  234.         K = N+1-KB
  235.         DO 61 I = 1, N
  236.           LL = L+I-1+(K-1)*N
  237.           LU = LL+NN
  238.           IF (I.LE.K) STKR(LU) = STKR(LL)
  239.           IF (I.LE.K) STKI(LU) = STKI(LL)
  240.           IF (I.GT.K) STKR(LU) = 0.0D0
  241.           IF (I.GT.K) STKI(LU) = 0.0D0
  242.           IF (I.LT.K) STKR(LL) = 0.0D0
  243.           IF (I.LT.K) STKI(LL) = 0.0D0
  244.           IF (I.EQ.K) STKR(LL) = 1.0D0
  245.           IF (I.EQ.K) STKI(LL) = 0.0D0
  246.           IF (I.GT.K) STKR(LL) = -STKR(LL)
  247.           IF (I.GT.K) STKI(LL) = -STKI(LL)
  248. 61      CONTINUE
  249.         I = BUF(K)
  250.         IF (I.EQ.K) GO TO 64
  251.         LI = L+I-1+(K-1)*N
  252.         LK = L+K-1+(K-1)*N
  253.         CALL WSWAP (N-K+1, STKR(LI), STKI(LI), N,
  254.      .              STKR(LK), STKI(LK), N)
  255. 64    CONTINUE
  256.       GO TO 99
  257. C
  258. C ***      HILBERT
  259. 70    CONTINUE
  260.       N = IDINT (STKR(L))
  261.       MSTK(TOP) = N
  262.       NSTK(TOP) = N
  263. 72    CONTINUE
  264.       CALL HILBER (STKR(L), N, N)
  265.       CALL RSET (N*N, 0.0D0, STKI(L), 1)
  266.       IF (FIN.LT.0) CALL WSCAL (N*N, SR, SI, STKR(L), STKI(L), 1)
  267.       GO TO 99
  268. C
  269. C ***      CHOLESKY
  270. 80    CONTINUE
  271.       IF (M.NE.N) CALL ERROR (20)
  272.       IF (ERR.GT.0) RETURN
  273.       CALL WPOFA (STKR(L), STKI(L), M, N, ERR)
  274.       IF (ERR.NE.0) CALL ERROR (29)
  275.       IF (ERR.GT.0) RETURN
  276.       DO 81 J = 1, N
  277.         LL = L+J+(J-1)*M
  278.         CALL WSET (M-J, 0.0D0, 0.0D0, STKR(LL), STKI(LL), 1)
  279. 81    CONTINUE
  280.       GO TO 99
  281. C
  282. C ***      RREF
  283. 85    CONTINUE
  284.       IF (RHS.LT.2) GO TO 86
  285.       TOP = TOP-1
  286.       L = LSTK(TOP)
  287.       IF (MSTK(TOP).NE.M) CALL ERROR (5)
  288.       IF (ERR.GT.0) RETURN
  289.       N = N+NSTK(TOP)
  290. 86    CONTINUE
  291.       CALL RREF (STKR(L), STKI(L), M, M, N, STKR(VSIZE-4))
  292.       NSTK(TOP) = N
  293.       GO TO 99
  294. C
  295. 99    CONTINUE
  296.       RETURN
  297.       END
  298.