home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / imtql2.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  156 lines
  1.       subroutine imtql2(nm,n,d,e,z,ierr)
  2. c
  3.       integer i,j,k,l,m,n,ii,nm,mml,ierr
  4.       double precision d(n),e(n),z(nm,n)
  5.       double precision b,c,f,g,p,r,s,tst1,tst2,pythag
  6. c
  7. c     this subroutine is a translation of the algol procedure imtql2,
  8. c     num. math. 12, 377-383(1968) by martin and wilkinson,
  9. c     as modified in num. math. 15, 450(1970) by dubrulle.
  10. c     handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
  11. c
  12. c     this subroutine finds the eigenvalues and eigenvectors
  13. c     of a symmetric tridiagonal matrix by the implicit ql method.
  14. c     the eigenvectors of a full symmetric matrix can also
  15. c     be found if  tred2  has been used to reduce this
  16. c     full matrix to tridiagonal form.
  17. c
  18. c     on input
  19. c
  20. c        nm must be set to the row dimension of two-dimensional
  21. c          array parameters as declared in the calling program
  22. c          dimension statement.
  23. c
  24. c        n is the order of the matrix.
  25. c
  26. c        d contains the diagonal elements of the input matrix.
  27. c
  28. c        e contains the subdiagonal elements of the input matrix
  29. c          in its last n-1 positions.  e(1) is arbitrary.
  30. c
  31. c        z contains the transformation matrix produced in the
  32. c          reduction by  tred2, if performed.  if the eigenvectors
  33. c          of the tridiagonal matrix are desired, z must contain
  34. c          the identity matrix.
  35. c
  36. c      on output
  37. c
  38. c        d contains the eigenvalues in ascending order.  if an
  39. c          error exit is made, the eigenvalues are correct but
  40. c          unordered for indices 1,2,...,ierr-1.
  41. c
  42. c        e has been destroyed.
  43. c
  44. c        z contains orthonormal eigenvectors of the symmetric
  45. c          tridiagonal (or full) matrix.  if an error exit is made,
  46. c          z contains the eigenvectors associated with the stored
  47. c          eigenvalues.
  48. c
  49. c        ierr is set to
  50. c          zero       for normal return,
  51. c          j          if the j-th eigenvalue has not been
  52. c                     determined after 30 iterations.
  53. c
  54. c     calls pythag for  dsqrt(a*a + b*b) .
  55. c
  56. c     questions and comments should be directed to burton s. garbow,
  57. c     mathematics and computer science div, argonne national laboratory
  58. c
  59. c     this version dated august 1983.
  60. c
  61. c     ------------------------------------------------------------------
  62. c
  63.       ierr = 0
  64.       if (n .eq. 1) go to 1001
  65. c
  66.       do 100 i = 2, n
  67.   100 e(i-1) = e(i)
  68. c
  69.       e(n) = 0.0d0
  70. c
  71.       do 240 l = 1, n
  72.          j = 0
  73. c     .......... look for small sub-diagonal element ..........
  74.   105    do 110 m = l, n
  75.             if (m .eq. n) go to 120
  76.             tst1 = dabs(d(m)) + dabs(d(m+1))
  77.             tst2 = tst1 + dabs(e(m))
  78.             if (tst2 .eq. tst1) go to 120
  79.   110    continue
  80. c
  81.   120    p = d(l)
  82.          if (m .eq. l) go to 240
  83.          if (j .eq. 30) go to 1000
  84.          j = j + 1
  85. c     .......... form shift ..........
  86.          g = (d(l+1) - p) / (2.0d0 * e(l))
  87.          r = pythag(g,1.0d0)
  88.          g = d(m) - p + e(l) / (g + dsign(r,g))
  89.          s = 1.0d0
  90.          c = 1.0d0
  91.          p = 0.0d0
  92.          mml = m - l
  93. c     .......... for i=m-1 step -1 until l do -- ..........
  94.          do 200 ii = 1, mml
  95.             i = m - ii
  96.             f = s * e(i)
  97.             b = c * e(i)
  98.             r = pythag(f,g)
  99.             e(i+1) = r
  100.             if (r .eq. 0.0d0) go to 210
  101.             s = f / r
  102.             c = g / r
  103.             g = d(i+1) - p
  104.             r = (d(i) - g) * s + 2.0d0 * c * b
  105.             p = s * r
  106.             d(i+1) = g + p
  107.             g = c * r - b
  108. c     .......... form vector ..........
  109.             do 180 k = 1, n
  110.                f = z(k,i+1)
  111.                z(k,i+1) = s * z(k,i) + c * f
  112.                z(k,i) = c * z(k,i) - s * f
  113.   180       continue
  114. c
  115.   200    continue
  116. c
  117.          d(l) = d(l) - p
  118.          e(l) = g
  119.          e(m) = 0.0d0
  120.          go to 105
  121. c     .......... recover from underflow ..........
  122.   210    d(i+1) = d(i+1) - p
  123.          e(m) = 0.0d0
  124.          go to 105
  125.   240 continue
  126. c     .......... order eigenvalues and eigenvectors ..........
  127.       do 300 ii = 2, n
  128.          i = ii - 1
  129.          k = i
  130.          p = d(i)
  131. c
  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. c
  138.          if (k .eq. i) go to 300
  139.          d(k) = d(i)
  140.          d(i) = p
  141. c
  142.          do 280 j = 1, n
  143.             p = z(j,i)
  144.             z(j,i) = z(j,k)
  145.             z(j,k) = p
  146.   280    continue
  147. c
  148.   300 continue
  149. c
  150.       go to 1001
  151. c     .......... set error -- no convergence to an
  152. c                eigenvalue after 30 iterations ..........
  153.  1000 ierr = l
  154.  1001 return
  155.       end
  156.