home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / imtqlv.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  150 lines
  1.       subroutine imtqlv(n,d,e,e2,w,ind,ierr,rv1)
  2. c
  3.       integer i,j,k,l,m,n,ii,mml,tag,ierr
  4.       double precision d(n),e(n),e2(n),w(n),rv1(n)
  5.       double precision b,c,f,g,p,r,s,tst1,tst2,pythag
  6.       integer ind(n)
  7. c
  8. c     this subroutine is a variant of  imtql1  which is a translation of
  9. c     algol procedure imtql1, num. math. 12, 377-383(1968) by martin and
  10. c     wilkinson, as modified in num. math. 15, 450(1970) by dubrulle.
  11. c     handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
  12. c
  13. c     this subroutine finds the eigenvalues of a symmetric tridiagonal
  14. c     matrix by the implicit ql method and associates with them
  15. c     their corresponding submatrix indices.
  16. c
  17. c     on input
  18. c
  19. c        n is the order of the matrix.
  20. c
  21. c        d contains the diagonal elements of the input matrix.
  22. c
  23. c        e contains the subdiagonal elements of the input matrix
  24. c          in its last n-1 positions.  e(1) is arbitrary.
  25. c
  26. c        e2 contains the squares of the corresponding elements of e.
  27. c          e2(1) is arbitrary.
  28. c
  29. c     on output
  30. c
  31. c        d and e are unaltered.
  32. c
  33. c        elements of e2, corresponding to elements of e regarded
  34. c          as negligible, have been replaced by zero causing the
  35. c          matrix to split into a direct sum of submatrices.
  36. c          e2(1) is also set to zero.
  37. c
  38. c        w contains the eigenvalues in ascending order.  if an
  39. c          error exit is made, the eigenvalues are correct and
  40. c          ordered for indices 1,2,...ierr-1, but may not be
  41. c          the smallest eigenvalues.
  42. c
  43. c        ind contains the submatrix indices associated with the
  44. c          corresponding eigenvalues in w -- 1 for eigenvalues
  45. c          belonging to the first submatrix from the top,
  46. c          2 for those belonging to the second submatrix, etc..
  47. c
  48. c        ierr is set to
  49. c          zero       for normal return,
  50. c          j          if the j-th eigenvalue has not been
  51. c                     determined after 30 iterations.
  52. c
  53. c        rv1 is a temporary storage array.
  54. c
  55. c     calls pythag for  dsqrt(a*a + b*b) .
  56. c
  57. c     questions and comments should be directed to burton s. garbow,
  58. c     mathematics and computer science div, argonne national laboratory
  59. c
  60. c     this version dated august 1983.
  61. c
  62. c     ------------------------------------------------------------------
  63. c
  64.       ierr = 0
  65.       k = 0
  66.       tag = 0
  67. c
  68.       do 100 i = 1, n
  69.          w(i) = d(i)
  70.          if (i .ne. 1) rv1(i-1) = e(i)
  71.   100 continue
  72. c
  73.       e2(1) = 0.0d0
  74.       rv1(n) = 0.0d0
  75. c
  76.       do 290 l = 1, n
  77.          j = 0
  78. c     .......... look for small sub-diagonal element ..........
  79.   105    do 110 m = l, n
  80.             if (m .eq. n) go to 120
  81.             tst1 = dabs(w(m)) + dabs(w(m+1))
  82.             tst2 = tst1 + dabs(rv1(m))
  83.             if (tst2 .eq. tst1) go to 120
  84. c     .......... guard against underflowed element of e2 ..........
  85.             if (e2(m+1) .eq. 0.0d0) go to 125
  86.   110    continue
  87. c
  88.   120    if (m .le. k) go to 130
  89.          if (m .ne. n) e2(m+1) = 0.0d0
  90.   125    k = m
  91.          tag = tag + 1
  92.   130    p = w(l)
  93.          if (m .eq. l) go to 215
  94.          if (j .eq. 30) go to 1000
  95.          j = j + 1
  96. c     .......... form shift ..........
  97.          g = (w(l+1) - p) / (2.0d0 * rv1(l))
  98.          r = pythag(g,1.0d0)
  99.          g = w(m) - p + rv1(l) / (g + dsign(r,g))
  100.          s = 1.0d0
  101.          c = 1.0d0
  102.          p = 0.0d0
  103.          mml = m - l
  104. c     .......... for i=m-1 step -1 until l do -- ..........
  105.          do 200 ii = 1, mml
  106.             i = m - ii
  107.             f = s * rv1(i)
  108.             b = c * rv1(i)
  109.             r = pythag(f,g)
  110.             rv1(i+1) = r
  111.             if (r .eq. 0.0d0) go to 210
  112.             s = f / r
  113.             c = g / r
  114.             g = w(i+1) - p
  115.             r = (w(i) - g) * s + 2.0d0 * c * b
  116.             p = s * r
  117.             w(i+1) = g + p
  118.             g = c * r - b
  119.   200    continue
  120. c
  121.          w(l) = w(l) - p
  122.          rv1(l) = g
  123.          rv1(m) = 0.0d0
  124.          go to 105
  125. c     .......... recover from underflow ..........
  126.   210    w(i+1) = w(i+1) - p
  127.          rv1(m) = 0.0d0
  128.          go to 105
  129. c     .......... order eigenvalues ..........
  130.   215    if (l .eq. 1) go to 250
  131. c     .......... for i=l step -1 until 2 do -- ..........
  132.          do 230 ii = 2, l
  133.             i = l + 2 - ii
  134.             if (p .ge. w(i-1)) go to 270
  135.             w(i) = w(i-1)
  136.             ind(i) = ind(i-1)
  137.   230    continue
  138. c
  139.   250    i = 1
  140.   270    w(i) = p
  141.          ind(i) = tag
  142.   290 continue
  143. c
  144.       go to 1001
  145. c     .......... set error -- no convergence to an
  146. c                eigenvalue after 30 iterations ..........
  147.  1000 ierr = l
  148.  1001 return
  149.       end
  150.