home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / otqlrat.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  131 lines
  1.       subroutine tqlrat(n,d,e2,ierr)
  2. c
  3.       integer i,j,l,m,n,ii,l1,mml,ierr
  4.       double precision d(n),e2(n)
  5.       double precision b,c,f,g,h,p,r,s,t,epslon,pythag
  6. c
  7. c     this subroutine is a translation of the algol procedure tqlrat,
  8. c     algorithm 464, comm. acm 16, 689(1973) by reinsch.
  9. c
  10. c     this subroutine finds the eigenvalues of a symmetric
  11. c     tridiagonal matrix by the rational ql method.
  12. c
  13. c     on input
  14. c
  15. c        n is the order of the matrix.
  16. c
  17. c        d contains the diagonal elements of the input matrix.
  18. c
  19. c        e2 contains the squares of the subdiagonal elements of the
  20. c          input matrix in its last n-1 positions.  e2(1) is arbitrary.
  21. c
  22. c      on output
  23. c
  24. c        d contains the eigenvalues in ascending order.  if an
  25. c          error exit is made, the eigenvalues are correct and
  26. c          ordered for indices 1,2,...ierr-1, but may not be
  27. c          the smallest eigenvalues.
  28. c
  29. c        e2 has been destroyed.
  30. c
  31. c        ierr is set to
  32. c          zero       for normal return,
  33. c          j          if the j-th eigenvalue has not been
  34. c                     determined after 30 iterations.
  35. c
  36. c     calls pythag for  dsqrt(a*a + b*b) .
  37. c
  38. c     questions and comments should be directed to burton s. garbow,
  39. c     mathematics and computer science div, argonne national laboratory
  40. c
  41. c     this version dated august 1983.
  42. c
  43. c     ------------------------------------------------------------------
  44. c
  45.       ierr = 0
  46.       if (n .eq. 1) go to 1001
  47. c
  48.       do 100 i = 2, n
  49.   100 e2(i-1) = e2(i)
  50. c
  51.       f = 0.0d0
  52.       t = 0.0d0
  53.       e2(n) = 0.0d0
  54. c
  55.       do 290 l = 1, n
  56.          j = 0
  57.          h = dabs(d(l)) + dsqrt(e2(l))
  58.          if (t .gt. h) go to 105
  59.          t = h
  60.          b = epslon(t)
  61.          c = b * b
  62. c     .......... look for small squared sub-diagonal element ..........
  63.   105    do 110 m = l, n
  64.             if (e2(m) .le. c) go to 120
  65. c     .......... e2(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.          s = dsqrt(e2(l))
  75.          g = d(l)
  76.          p = (d(l1) - g) / (2.0d0 * s)
  77.          r = pythag(p,1.0d0)
  78.          d(l) = s / (p + dsign(r,p))
  79.          h = g - d(l)
  80. c
  81.          do 140 i = l1, n
  82.   140    d(i) = d(i) - h
  83. c
  84.          f = f + h
  85. c     .......... rational ql transformation ..........
  86.          g = d(m)
  87.          if (g .eq. 0.0d0) g = b
  88.          h = g
  89.          s = 0.0d0
  90.          mml = m - l
  91. c     .......... for i=m-1 step -1 until l do -- ..........
  92.          do 200 ii = 1, mml
  93.             i = m - ii
  94.             p = g * h
  95.             r = p + e2(i)
  96.             e2(i+1) = s * r
  97.             s = e2(i) / r
  98.             d(i+1) = h + s * (h + d(i))
  99.             g = d(i) - e2(i) / g
  100.             if (g .eq. 0.0d0) g = b
  101.             h = g * p / r
  102.   200    continue
  103. c
  104.          e2(l) = s * g
  105.          d(l) = h
  106. c     .......... guard against underflow in convergence test ..........
  107.          if (h .eq. 0.0d0) go to 210
  108.          if (dabs(e2(l)) .le. dabs(c/h)) go to 210
  109.          e2(l) = h * e2(l)
  110.          if (e2(l) .ne. 0.0d0) go to 130
  111.   210    p = d(l) + f
  112. c     .......... order eigenvalues ..........
  113.          if (l .eq. 1) go to 250
  114. c     .......... for i=l step -1 until 2 do -- ..........
  115.          do 230 ii = 2, l
  116.             i = l + 2 - ii
  117.             if (p .ge. d(i-1)) go to 270
  118.             d(i) = d(i-1)
  119.   230    continue
  120. c
  121.   250    i = 1
  122.   270    d(i) = p
  123.   290 continue
  124. c
  125.       go to 1001
  126. c     .......... set error -- no convergence to an
  127. c                eigenvalue after 30 iterations ..........
  128.  1000 ierr = l
  129.  1001 return
  130.       end
  131.