home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tqlrat.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  147 lines
  1. **** for old version, "send otqlrat from eispack"
  2. ** From dana!moler Tue, 1 Sep 87 10:15:40 PDT
  3. ** New TQLRAT
  4.       SUBROUTINE TQLRAT(N,D,E2,IERR)
  5. C
  6.       INTEGER I,J,L,M,N,II,L1,MML,IERR
  7.       DOUBLE PRECISION D(N),E2(N)
  8.       DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG
  9. C
  10. C     This subroutine is a translation of the Algol procedure tqlrat,
  11. C     Algorithm 464, Comm. ACM 16, 689(1973) by Reinsch.
  12. C
  13. C     This subroutine finds the eigenvalues of a symmetric
  14. C     tridiagonal matrix by the rational QL method.
  15. C
  16. C     On input
  17. C
  18. C        N is the order of the matrix.
  19. C
  20. C        D contains the diagonal elements of the input matrix.
  21. C
  22. C        E2 contains the squares of the subdiagonal elements of the
  23. C          input matrix in its last N-1 positions.  E2(1) is arbitrary.
  24. C
  25. C      On output
  26. C
  27. C        D contains the eigenvalues in ascending order.  If an
  28. C          error exit is made, the eigenvalues are correct and
  29. C          ordered for indices 1,2,...IERR-1, but may not be
  30. C          the smallest eigenvalues.
  31. C
  32. C        E2 has been destroyed.
  33. C
  34. C        IERR is set to
  35. C          zero       for normal return,
  36. C          J          if the J-th eigenvalue has not been
  37. C                     determined after 30 iterations.
  38. C
  39. C     Calls PYTHAG for  DSQRT(A*A + B*B) .
  40. C
  41. C     Questions and comments should be directed to Burton S. Garbow,
  42. C     Mathematics and Computer Science Div, Argonne National Laboratory
  43. C
  44. C     This version dated August 1987.
  45. C     Modified by C. Moler to fix underflow/overflow difficulties,
  46. C     especially on the VAX and other machines where epslon(1.0d0)**2
  47. C     nearly underflows.  See the loop involving statement 102 and
  48. C     the two statements just before statement 200.
  49. C
  50. C     ------------------------------------------------------------------
  51. C
  52.       IERR = 0
  53.       IF (N .EQ. 1) GO TO 1001
  54. C
  55.       DO 100 I = 2, N
  56.   100 E2(I-1) = E2(I)
  57. C
  58.       F = 0.0D0
  59.       T = 0.0D0
  60.       E2(N) = 0.0D0
  61. C
  62.       DO 290 L = 1, N
  63.          J = 0
  64.          H = DABS(D(L)) + DSQRT(E2(L))
  65.          IF (T .GT. H) GO TO 105
  66.          T = H
  67.          B = EPSLON(T)
  68.          C = B * B
  69.          if (c .ne. 0.0d0) go to 105
  70. C        Spliting tolerance underflowed.  Look for larger value.
  71.          do 102 i = l, n
  72.             h = dabs(d(i)) + dsqrt(e2(i))
  73.             if (h .gt. t) t = h
  74.   102    continue
  75.          b = epslon(t)
  76.          c = b * b
  77. C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
  78.   105    DO 110 M = L, N
  79.             IF (E2(M) .LE. C) GO TO 120
  80. C     .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
  81. C                THROUGH THE BOTTOM OF THE LOOP ..........
  82.   110    CONTINUE
  83. C
  84.   120    IF (M .EQ. L) GO TO 210
  85.   130    IF (J .EQ. 30) GO TO 1000
  86.          J = J + 1
  87. C     .......... FORM SHIFT ..........
  88.          L1 = L + 1
  89.          S = DSQRT(E2(L))
  90.          G = D(L)
  91.          P = (D(L1) - G) / (2.0D0 * S)
  92.          R = PYTHAG(P,1.0D0)
  93.          D(L) = S / (P + DSIGN(R,P))
  94.          H = G - D(L)
  95. C
  96.          DO 140 I = L1, N
  97.   140    D(I) = D(I) - H
  98. C
  99.          F = F + H
  100. C     .......... RATIONAL QL TRANSFORMATION ..........
  101.          G = D(M)
  102.          IF (G .EQ. 0.0D0) G = B
  103.          H = G
  104.          S = 0.0D0
  105.          MML = M - L
  106. C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
  107.          DO 200 II = 1, MML
  108.             I = M - II
  109.             P = G * H
  110.             R = P + E2(I)
  111.             E2(I+1) = S * R
  112.             S = E2(I) / R
  113.             D(I+1) = H + S * (H + D(I))
  114.             G = D(I) - E2(I) / G
  115. C           Avoid division by zero on next pass
  116.             if (g .eq. 0.0d0) g = epslon(d(i))
  117.             h = g * (p / r)
  118.   200    CONTINUE
  119. C
  120.          E2(L) = S * G
  121.          D(L) = H
  122. C     .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
  123.          IF (H .EQ. 0.0D0) GO TO 210
  124.          IF (DABS(E2(L)) .LE. DABS(C/H)) GO TO 210
  125.          E2(L) = H * E2(L)
  126.          IF (E2(L) .NE. 0.0D0) GO TO 130
  127.   210    P = D(L) + F
  128. C     .......... ORDER EIGENVALUES ..........
  129.          IF (L .EQ. 1) GO TO 250
  130. C     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
  131.          DO 230 II = 2, L
  132.             I = L + 2 - II
  133.             IF (P .GE. D(I-1)) GO TO 270
  134.             D(I) = D(I-1)
  135.   230    CONTINUE
  136. C
  137.   250    I = 1
  138.   270    D(I) = P
  139.   290 CONTINUE
  140. C
  141.       GO TO 1001
  142. C     .......... SET ERROR -- NO CONVERGENCE TO AN
  143. C                EIGENVALUE AFTER 30 ITERATIONS ..........
  144.  1000 IERR = L
  145.  1001 RETURN
  146.       END
  147.