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 / matfn2.for < prev    next >
Encoding:
Text File  |  1991-04-13  |  7.5 KB  |  258 lines
  1.       SUBROUTINE MATFN2
  2.       IMPLICIT NONE
  3. C
  4. C EVALUATE ELEMENTARY FUNCTIONS & FUNCTIONS INVOLVING EIGENVALUES/EIGENVECTORS
  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, L, L1, L2, LD, LE, LW, LS, LL, LJ
  12.       INTEGER M, N, NN, INC, JOB
  13.       DOUBLE PRECISION TR, TI, SR, SI, POWR, POWI
  14.       LOGICAL HERM, SCHUR, VECT, HESS
  15. C
  16.       DOUBLE PRECISION PYTHAG, FLOP, ROUND
  17. C
  18. C
  19.       IF (DDT.EQ.1) WRITE (WTE, 100) FIN
  20. 100   FORMAT (' MATFN2', I4)
  21. C
  22. C FUNCTIONS/FIN
  23. C     **   SIN  COS ATAN  EXP  SQRT LOG
  24. C      0    1    2    3    4    5    6
  25. C    EIG  SCHU HESS POLY ROOT
  26. C     11   12   13   14   15
  27. C    ABS  ROUN REAL IMAG CONJ
  28. C     21   22   23   24   25
  29. C
  30.       IF (FIN.NE.0) GO TO 05
  31.       L = LSTK(TOP+1)
  32.       POWR = STKR(L)
  33.       POWI = STKI(L)
  34. 05    CONTINUE
  35.       L = LSTK(TOP)
  36.       M = MSTK(TOP)
  37.       N = NSTK(TOP)
  38.       IF (FIN.GE.11 .AND. FIN.LE.13) GO TO 10
  39.       IF (FIN.EQ.14 .AND. (M.EQ.1 .OR. N.EQ.1)) GO TO 50
  40.       IF (FIN.EQ.14) GO TO 10
  41.       IF (FIN.EQ.15) GO TO 60
  42.       IF (FIN.GT.20) GO TO 40
  43.       IF (M.EQ.1 .OR. N.EQ.1) GO TO 40
  44. C
  45. C ***      EIGENVALUES AND VECTORS
  46. 10    CONTINUE
  47.       IF (M.NE.N) CALL ERROR (20)
  48.       IF (ERR.GT.0) RETURN
  49.       SCHUR = FIN.EQ.12
  50.       HESS = FIN.EQ.13
  51.       VECT = LHS.EQ.2 .OR. FIN.LT.10
  52.       NN = N*N
  53.       L2 = L+NN
  54.       LD = L2+NN
  55.       LE = LD+N
  56.       LW = LE+N
  57.       ERR = LW+N-LSTK(BOT)
  58.       IF (ERR.GT.0) CALL ERROR (17)
  59.       IF (ERR.GT.0) RETURN
  60.       CALL WCOPY (NN, STKR(L), STKI(L), 1, STKR(L2), STKI(L2), 1)
  61. C
  62. C ***      CHECK IF HERMITIAN
  63.       DO 16 J = 1, N
  64.         DO 15 I = 1, J
  65.           LS = L+I-1+(J-1)*N
  66.           LL = L+(I-1)*N+J-1
  67.           HERM = STKR(LL).EQ.STKR(LS) .AND. STKI(LL).EQ.-STKI(LS)
  68.           IF (.NOT.HERM) GO TO 30
  69. 15      CONTINUE
  70. 16    CONTINUE
  71. C
  72. C ***      HERMITIAN EIGENVALUE PROBLEM
  73.       CALL WSET (NN, 0.0D0, 0.0D0, STKR(L), STKI(L), 1)
  74.       CALL WSET (N, 1.0D0, 0.0D0, STKR(L), STKI(L), N+1)
  75.       CALL WSET (N, 0.0D0, 0.0D0, STKI(LD), STKI(LE), 1)
  76.       JOB = 0
  77.       IF (VECT) JOB = 1
  78.       CALL HTRIDI (N, N, STKR(L2), STKI(L2), STKR(LD), STKR(LE),
  79.      .             STKR(LE), STKR(LW))
  80.       IF (.NOT.HESS) CALL IMTQL2 (N, N, STKR(LD), STKR(LE), STKR(L),
  81.      .                            ERR, JOB)
  82.       IF (ERR.GT.0) CALL ERROR (24)
  83.       IF (ERR.GT.0) RETURN
  84.       IF (JOB.NE.0) CALL HTRIBK (N, N, STKR(L2), STKI(L2), STKR(LW),
  85.      .                           N, STKR(L), STKI(L))
  86.       GO TO 31
  87. C
  88. C ***      NON-HERMITIAN EIGENVALUE PROBLEM
  89. 30    CONTINUE
  90.       CALL CORTH (N, N, 1, N, STKR(L2), STKI(L2), STKR(LW), STKI(LW))
  91.       IF (.NOT.VECT .AND. HESS) GO TO 31
  92.       JOB = 0
  93.       IF (VECT) JOB = 2
  94.       IF (VECT.AND.SCHUR) JOB = 1
  95.       IF (HESS) JOB = 3
  96.       CALL COMQR3 (N, N, 1, N, STKR(LW), STKI(LW), STKR(L2), STKI(L2),
  97.      .             STKR(LD), STKI(LD), STKR(L), STKI(L), ERR, JOB)
  98.       IF (ERR.GT.0) CALL ERROR (24)
  99.       IF (ERR.GT.0) RETURN
  100. C
  101. C ***      VECTORS
  102. 31    CONTINUE
  103.       IF (.NOT.VECT) GO TO 34
  104.       IF (TOP+1.GE.BOT) CALL ERROR (18)
  105.       IF (ERR.GT.0) RETURN
  106.       TOP = TOP+1
  107.       LSTK(TOP) = L2
  108.       MSTK(TOP) = N
  109.       NSTK(TOP) = N
  110. C
  111. C ***      DIAGONAL OF VALUES OR CANONICAL FORMS
  112. 34    CONTINUE
  113.       IF (.NOT.VECT .AND. .NOT.SCHUR .AND. .NOT.HESS) GO TO 37
  114.       DO 36 J = 1, N
  115.         LJ = L2+(J-1)*N
  116.         IF (SCHUR .AND. (.NOT.HERM)) LJ = LJ+J
  117.         IF (HESS .AND. (.NOT.HERM)) LJ = LJ+J+1
  118.         LL = L2+J*N-LJ
  119.         CALL WSET (LL, 0.0D0, 0.0D0, STKR(LJ), STKI(LJ), 1)
  120. 36    CONTINUE
  121.       IF (.NOT.HESS .OR. HERM) CALL WCOPY (N, STKR(LD), STKI(LD),
  122.      .                                     1, STKR(L2), STKI(L2), N+1)
  123.       LL = L2+1
  124.       IF (HESS.AND.HERM) CALL WCOPY (N-1, STKR(LE+1), STKI(LE+1),
  125.      .                               1, STKR(LL), STKI(LL), N+1)
  126.       LL = L2+N
  127.       IF (HESS.AND.HERM) CALL WCOPY (N-1, STKR(LE+1), STKI(LE+1),
  128.      .                               1, STKR(LL), STKI(LL), N+1)
  129.       IF (FIN.LT.10) GO TO 42
  130.       IF (VECT .OR. .NOT.(SCHUR.OR.HESS)) GO TO 99
  131.       CALL WCOPY (NN, STKR(L2), STKI(L2), 1, STKR(L), STKI(L), 1)
  132.       GO TO 99
  133. C
  134. C ***      VECTOR OF EIGENVALUES
  135. 37    CONTINUE
  136.       IF (FIN.EQ.14) GO TO 52
  137.       CALL WCOPY (N, STKR(LD), STKI(LD), 1, STKR(L), STKI(L), 1)
  138.       NSTK(TOP) = 1
  139.       GO TO 99
  140. C
  141. C ***      ELEMENTARY FUNCTIONS
  142. C ***        FOR MATRICES.. X,D = EIG(A), FUN(A) = X*FUN(D)/X
  143. 40    CONTINUE
  144.       INC = 1
  145.       N = M*N
  146.       L2 = L
  147.       GO TO 44
  148. C
  149. 42    CONTINUE
  150.       INC = N+1
  151. 44    CONTINUE
  152.       DO 46 J = 1, N
  153.         LS = L2+(J-1)*INC
  154.         SR = STKR(LS)
  155.         SI = STKI(LS)
  156.         TI = 0.0D0
  157.         IF (FIN.NE.0) GO TO 45
  158.         CALL WLOG (SR, SI, SR, SI)
  159.         CALL WMUL (SR, SI, POWR, POWI, SR, SI)
  160.         TR = DEXP (SR)*DCOS (SI)
  161.         TI = DEXP (SR)*DSIN (SI)
  162. 45      CONTINUE
  163.         IF (FIN.EQ.1) TR = DSIN (SR)*DCOSH (SI)
  164.         IF (FIN.EQ.1) TI = DCOS (SR)*DSINH (SI)
  165.         IF (FIN.EQ.2) TR = DCOS (SR)*DCOSH (SI)
  166.         IF (FIN.EQ.2) TI = -DSIN (SR)*DSINH (SI)
  167.         IF (FIN.EQ.3) CALL WATAN (SR, SI, TR, TI)
  168.         IF (FIN.EQ.4) TR = DEXP (SR)*DCOS (SI)
  169.         IF (FIN.EQ.4) TI = DEXP (SR)*DSIN (SI)
  170.         IF (FIN.EQ.5) CALL WSQRT (SR, SI, TR, TI)
  171.         IF (FIN.EQ.6) CALL WLOG (SR, SI, TR, TI)
  172.         IF (FIN.EQ.21) TR = PYTHAG (SR, SI)
  173.         IF (FIN.EQ.22) TR = ROUND (SR)
  174.         IF (FIN.EQ.23) TR = SR
  175.         IF (FIN.EQ.24) TR = SI
  176.         IF (FIN.EQ.25) TR = SR
  177.         IF (FIN.EQ.25) TI = -SI
  178.         IF (ERR.GT.0) RETURN
  179.         STKR(LS) = FLOP (TR)
  180.         STKI(LS) = 0.0D0
  181.         IF (TI.NE.0.0D0) STKI(LS) = FLOP (TI)
  182. 46    CONTINUE
  183.       IF (INC.EQ.1) GO TO 99
  184.       DO 48 J = 1, N
  185.         LS = L2+(J-1)*INC
  186.         SR = STKR(LS)
  187.         SI = STKI(LS)
  188.         LS = L+(J-1)*N
  189.         LL = L2+(J-1)*N
  190.         CALL WCOPY (N, STKR(LS), STKI(LS), 1, STKR(LL), STKI(LL), 1)
  191.         CALL WSCAL (N, SR, SI, STKR(LS), STKI(LS), 1)
  192. 48    CONTINUE
  193. C
  194. C ***      SIGNAL MATFN1 TO DIVIDE BY EIGENVECTORS
  195.       FUN = 21
  196.       FIN = -1
  197.       TOP = TOP-1
  198.       GO TO 99
  199. C
  200. C ***      POLY:  FORM POLYNOMIAL WITH GIVEN VECTOR AS ROOTS
  201. 50    CONTINUE
  202.       N = MAX0 (M, N)
  203.       LD = L+N+1
  204.       CALL WCOPY (N, STKR(L), STKI(L), 1, STKR(LD), STKI(LD), 1)
  205. C
  206. C ***      FORM CHARACTERISTIC POLYNOMIAL
  207. 52    CONTINUE
  208.       CALL WSET (N+1, 0.0D0, 0.0D0, STKR(L), STKI(L), 1)
  209.       STKR(L) = 1.0D0
  210.       DO 56 J = 1, N
  211.         CALL WAXPY (J, -STKR(LD), -STKI(LD), STKR(L), STKI(L),
  212.      .              -1, STKR(L+1), STKI(L+1), -1)
  213.         LD = LD+1
  214. 56    CONTINUE
  215.       MSTK(TOP) = N+1
  216.       NSTK(TOP) = 1
  217.       GO TO 99
  218. C
  219. C ***      ROOTS
  220. 60    CONTINUE
  221.       LL = L+M*N
  222.       STKR(LL) = -1.0D0
  223.       STKI(LL) = 0.0D0
  224.       K = -1
  225. 61    CONTINUE
  226.       K = K+1
  227.       L1 = L+K
  228.       IF (DABS (STKR(L1))+DABS (STKI(L1)).EQ.0.0D0) GO TO 61
  229.       N = MAX0 (M*N-K-1,  0)
  230.       IF (N.LE.0) GO TO 65
  231.       L2 = L1+N+1
  232.       LW = L2+N*N
  233.       ERR = LW+N-LSTK(BOT)
  234.       IF (ERR.GT.0) CALL ERROR (17)
  235.       IF (ERR.GT.0) RETURN
  236.       CALL WSET (N*N+N, 0.0D0, 0.0D0, STKR(L2), STKI(L2), 1)
  237.       DO 64 J = 1, N
  238.         LL = L2+J+(J-1)*N
  239.         STKR(LL) = 1.0D0
  240.         LS = L1+J
  241.         LL = L2+(J-1)*N
  242.         CALL WDIV (-STKR(LS), -STKI(LS), STKR(L1), STKI(L1),
  243.      .             STKR(LL), STKI(LL))
  244.         IF (ERR.GT.0) RETURN
  245. 64    CONTINUE
  246.       CALL COMQR3 (N, N, 1, N, STKR(LW), STKI(LW), STKR(L2), STKI(L2),
  247.      .             STKR(L), STKI(L), TR, TI, ERR, 0)
  248.       IF (ERR.GT.0) CALL ERROR (24)
  249.       IF (ERR.GT.0) RETURN
  250. 65    CONTINUE
  251.       MSTK(TOP) = N
  252.       NSTK(TOP) = 1
  253.       GO TO 99
  254. C
  255. 99    CONTINUE
  256.       RETURN
  257.       END
  258.