home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tql1.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  136 lines
  1.       subroutine tql1(n,d,e,ierr)
  2. c
  3.       integer i,j,l,m,n,ii,l1,l2,mml,ierr
  4.       double precision d(n),e(n)
  5.       double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag
  6. c
  7. c     this subroutine is a translation of the algol procedure tql1,
  8. c     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
  9. c     wilkinson.
  10. c     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
  11. c
  12. c     this subroutine finds the eigenvalues of a symmetric
  13. c     tridiagonal matrix by the ql method.
  14. c
  15. c     on input
  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
  22. c          in its last n-1 positions.  e(1) is arbitrary.
  23. c
  24. c      on output
  25. c
  26. c        d contains the eigenvalues in ascending order.  if an
  27. c          error exit is made, the eigenvalues are correct and
  28. c          ordered for indices 1,2,...ierr-1, but may not be
  29. c          the smallest eigenvalues.
  30. c
  31. c        e has been destroyed.
  32. c
  33. c        ierr is set to
  34. c          zero       for normal return,
  35. c          j          if the j-th eigenvalue has not been
  36. c                     determined after 30 iterations.
  37. c
  38. c     calls pythag for  dsqrt(a*a + b*b) .
  39. c
  40. c     questions and comments should be directed to burton s. garbow,
  41. c     mathematics and computer science div, argonne national laboratory
  42. c
  43. c     this version dated august 1983.
  44. c
  45. c     ------------------------------------------------------------------
  46. c
  47.       ierr = 0
  48.       if (n .eq. 1) go to 1001
  49. c
  50.       do 100 i = 2, n
  51.   100 e(i-1) = e(i)
  52. c
  53.       f = 0.0d0
  54.       tst1 = 0.0d0
  55.       e(n) = 0.0d0
  56. c
  57.       do 290 l = 1, n
  58.          j = 0
  59.          h = dabs(d(l)) + dabs(e(l))
  60.          if (tst1 .lt. h) tst1 = h
  61. c     .......... look for small sub-diagonal element ..........
  62.          do 110 m = l, n
  63.             tst2 = tst1 + dabs(e(m))
  64.             if (tst2 .eq. tst1) go to 120
  65. c     .......... e(n) is always zero, so there is no exit
  66. c                through the bottom of the loop ..........
  67.   110    continue
  68. c
  69.   120    if (m .eq. l) go to 210
  70.   130    if (j .eq. 30) go to 1000
  71.          j = j + 1
  72. c     .......... form shift ..........
  73.          l1 = l + 1
  74.          l2 = l1 + 1
  75.          g = d(l)
  76.          p = (d(l1) - g) / (2.0d0 * e(l))
  77.          r = pythag(p,1.0d0)
  78.          d(l) = e(l) / (p + dsign(r,p))
  79.          d(l1) = e(l) * (p + dsign(r,p))
  80.          dl1 = d(l1)
  81.          h = g - d(l)
  82.          if (l2 .gt. n) go to 145
  83. c
  84.          do 140 i = l2, n
  85.   140    d(i) = d(i) - h
  86. c
  87.   145    f = f + h
  88. c     .......... ql transformation ..........
  89.          p = d(m)
  90.          c = 1.0d0
  91.          c2 = c
  92.          el1 = e(l1)
  93.          s = 0.0d0
  94.          mml = m - l
  95. c     .......... for i=m-1 step -1 until l do -- ..........
  96.          do 200 ii = 1, mml
  97.             c3 = c2
  98.             c2 = c
  99.             s2 = s
  100.             i = m - ii
  101.             g = c * e(i)
  102.             h = c * p
  103.             r = pythag(p,e(i))
  104.             e(i+1) = s * r
  105.             s = e(i) / r
  106.             c = p / r
  107.             p = c * d(i) - s * g
  108.             d(i+1) = h + s * (c * g + s * d(i))
  109.   200    continue
  110. c
  111.          p = -s * s2 * c3 * el1 * e(l) / dl1
  112.          e(l) = s * p
  113.          d(l) = c * p
  114.          tst2 = tst1 + dabs(e(l))
  115.          if (tst2 .gt. tst1) go to 130
  116.   210    p = d(l) + f
  117. c     .......... order eigenvalues ..........
  118.          if (l .eq. 1) go to 250
  119. c     .......... for i=l step -1 until 2 do -- ..........
  120.          do 230 ii = 2, l
  121.             i = l + 2 - ii
  122.             if (p .ge. d(i-1)) go to 270
  123.             d(i) = d(i-1)
  124.   230    continue
  125. c
  126.   250    i = 1
  127.   270    d(i) = p
  128.   290 continue
  129. c
  130.       go to 1001
  131. c     .......... set error -- no convergence to an
  132. c                eigenvalue after 30 iterations ..........
  133.  1000 ierr = l
  134.  1001 return
  135.       end
  136.