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 / IMTQL2.FOR < prev    next >
Encoding:
Text File  |  1991-04-13  |  4.3 KB  |  157 lines
  1.       SUBROUTINE IMTQL2 (NM, N, D, E, Z, IERR, JOB)
  2.       IMPLICIT NONE
  3. C
  4. C A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, NUM. MATH. 12, 377-383 (1968)
  5. C   BY MARTIN AND WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450 (1970) BY
  6. C   DUBRULLE.  HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248 (1971).
  7. C
  8. C FINDS THE EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC TRIDIAGONAL MATRIX
  9. C   BY THE IMPLICIT QL METHOD.  THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX
  10. C   CAN ALSO BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS FULL MATRIX TO
  11. C   TRIDIAGONAL FORM.
  12. C
  13. C ON INPUT:
  14. C   NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  15. C      AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
  16. C
  17. C   N IS THE ORDER OF THE MATRIX.
  18. C
  19. C   D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  20. C
  21. C   E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX IN ITS LAST
  22. C      N-1 POSITIONS.  E(1) IS ARBITRARY.
  23. C
  24. C   Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION BY TRED2,
  25. C      IF PERFORMED.  IF THE EIGENVECTORS OF THE TRIDIAGONAL MATRIX ARE
  26. C      DESIRED, Z MUST CONTAIN THE IDENTITY MATRIX.
  27. C
  28. C ON OUTPUT:
  29. C   D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN ERROR EXIT IS MADE,
  30. C      THE EIGENVALUES ARE CORRECT BUT UNORDERED FOR INDICES 1,2,...,IERR-1.
  31. C
  32. C   E HAS BEEN DESTROYED.
  33. C
  34. C   Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC TRIDIAGONAL (OR
  35. C      FULL) MATRIX.  IF AN ERROR EXIT IS MADE, Z CONTAINS THE EIGENVECTORS
  36. C      ASSOCIATED WITH THE STORED EIGENVALUES.
  37. C
  38. C   IERR = ZERO FOR NORMAL RETURN,
  39. C          J    IF J-TH EIGENVALUE HAS NOT BEEN DETERMINED AFTER 30 ITERATIONS.
  40. C
  41. C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  42. C   APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  43. C
  44.       INTEGER NM, N, IERR, JOB
  45.       DOUBLE PRECISION D(N), E(N), Z(NM,N)
  46. C
  47.       INTEGER I, J, K, L, M, II, MML
  48.       DOUBLE PRECISION B, C, F, G, P, R, S
  49. C
  50.       DOUBLE PRECISION FLOP
  51. C
  52. C
  53.       IERR = 0
  54.       IF (N.EQ.1) GO TO 1001
  55.       DO 100 I = 2, N
  56.         E(I-1) = E(I)
  57. 100   CONTINUE
  58. C
  59.       E(N) = 0.0D0
  60.       DO 240 L = 1, N
  61.         J = 0
  62. C
  63. C ***      LOOK FOR SMALL SUB-DIAGONAL ELEMENT
  64. 105     CONTINUE
  65.         DO 110 M = L, N
  66.           IF (M.EQ.N) GO TO 120
  67.           P = FLOP (DABS (D(M))+DABS (D(M+1)))
  68.           S = FLOP (P+DABS (E(M)))
  69.           IF (P.EQ.S) GO TO 120
  70. 110     CONTINUE
  71. 120     CONTINUE
  72.         P = D(L)
  73.         IF (M.EQ.L) GO TO 240
  74.         IF (J.EQ.30) GO TO 1000
  75.         J = J+1
  76. C
  77. C ***      FORM SHIFT
  78.         G = FLOP ((D(L+1)-P)/(2.0D0*E(L)))
  79.         R = FLOP (DSQRT (G*G+1.0D0))
  80.         G = FLOP (D(M)-P+E(L)/(G+DSIGN (R, G)))
  81.         S = 1.0D0
  82.         C = 1.0D0
  83.         P = 0.0D0
  84.         MML = M-L
  85. C
  86. C ***      FOR I = M-1 STEP -1 UNTIL L DO
  87.         DO 200 II = 1, MML
  88.           I = M-II
  89.           F = FLOP (S*E(I))
  90.           B = FLOP (C*E(I))
  91.           IF (DABS (F).LT.DABS (G)) GO TO 150
  92.           C = FLOP (G/F)
  93.           R = FLOP (DSQRT (C*C+1.0D0))
  94.           E(I+1) = FLOP (F*R)
  95.           S = FLOP (1.0D0/R)
  96.           C = FLOP (C*S)
  97.           GO TO 160
  98. C
  99. 150       CONTINUE
  100.           S = FLOP (F/G)
  101.           R = FLOP (DSQRT (S*S+1.0D0))
  102.           E(I+1) = FLOP (G*R)
  103.           C = FLOP (1.0D0/R)
  104.           S = FLOP (S*C)
  105. 160       CONTINUE
  106.           G = FLOP (D(I+1)-P)
  107.           R = FLOP ((D(I)-G)*S+2.0D0*C*B)
  108.           P = FLOP (S*R)
  109.           D(I+1) = G+P
  110.           G = FLOP (C*R-B)
  111.           IF (JOB.EQ.0) GO TO 185
  112. C
  113. C ***      FORM VECTOR
  114.           DO 180 K = 1, N
  115.             F = Z(K,I+1)
  116.             Z(K,I+1) = FLOP (S*Z(K,I)+C*F)
  117.             Z(K,I) = FLOP (C*Z(K,I)-S*F)
  118. 180       CONTINUE
  119. 185       CONTINUE
  120. 200     CONTINUE
  121.         D(L) = FLOP (D(L)-P)
  122.         E(L) = G
  123.         E(M) = 0.0D0
  124.         GO TO 105
  125. 240   CONTINUE
  126. C
  127. C ***      ORDER EIGENVALUES AND EIGENVECTORS
  128.       DO 300 II = 2, N
  129.         I = II-1
  130.         K = I
  131.         P = D(I)
  132.         DO 260 J = II, N
  133.           IF (D(J).GE.P) GO TO 260
  134.           K = J
  135.           P = D(J)
  136. 260     CONTINUE
  137.         IF (K.EQ.I) GO TO 300
  138.         D(K) = D(I)
  139.         D(I) = P
  140.         IF (JOB.EQ.0) GO TO 285
  141.         DO 280 J = 1, N
  142.           P = Z(J,I)
  143.           Z(J,I) = Z(J,K)
  144.           Z(J,K) = P
  145. 280     CONTINUE
  146. 285     CONTINUE
  147. 300   CONTINUE
  148.       GO TO 1001
  149. C
  150. C ***      SET ERROR -- NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS
  151. 1000  CONTINUE
  152.       IERR = L
  153. C
  154. 1001  CONTINUE
  155.       RETURN
  156.       END
  157.