home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / hqr.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  236 lines
  1.       subroutine hqr(nm,n,low,igh,h,wr,wi,ierr)
  2. C  RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG)
  3. c
  4.       integer i,j,k,l,m,n,en,ll,mm,na,nm,igh,itn,its,low,mp2,enm2,ierr
  5.       double precision h(nm,n),wr(n),wi(n)
  6.       double precision p,q,r,s,t,w,x,y,zz,norm,tst1,tst2
  7.       logical notlas
  8. c
  9. c     this subroutine is a translation of the algol procedure hqr,
  10. c     num. math. 14, 219-231(1970) by martin, peters, and wilkinson.
  11. c     handbook for auto. comp., vol.ii-linear algebra, 359-371(1971).
  12. c
  13. c     this subroutine finds the eigenvalues of a real
  14. c     upper hessenberg matrix by the qr method.
  15. c
  16. c     on input
  17. c
  18. c        nm must be set to the row dimension of two-dimensional
  19. c          array parameters as declared in the calling program
  20. c          dimension statement.
  21. c
  22. c        n is the order of the matrix.
  23. c
  24. c        low and igh are integers determined by the balancing
  25. c          subroutine  balanc.  if  balanc  has not been used,
  26. c          set low=1, igh=n.
  27. c
  28. c        h contains the upper hessenberg matrix.  information about
  29. c          the transformations used in the reduction to hessenberg
  30. c          form by  elmhes  or  orthes, if performed, is stored
  31. c          in the remaining triangle under the hessenberg matrix.
  32. c
  33. c     on output
  34. c
  35. c        h has been destroyed.  therefore, it must be saved
  36. c          before calling  hqr  if subsequent calculation and
  37. c          back transformation of eigenvectors is to be performed.
  38. c
  39. c        wr and wi contain the real and imaginary parts,
  40. c          respectively, of the eigenvalues.  the eigenvalues
  41. c          are unordered except that complex conjugate pairs
  42. c          of values appear consecutively with the eigenvalue
  43. c          having the positive imaginary part first.  if an
  44. c          error exit is made, the eigenvalues should be correct
  45. c          for indices ierr+1,...,n.
  46. c
  47. c        ierr is set to
  48. c          zero       for normal return,
  49. c          j          if the limit of 30*n iterations is exhausted
  50. c                     while the j-th eigenvalue is being sought.
  51. c
  52. c     questions and comments should be directed to burton s. garbow,
  53. c     mathematics and computer science div, argonne national laboratory
  54. c
  55. c     this version dated september 1989.
  56. c
  57. c     ------------------------------------------------------------------
  58. c
  59.       ierr = 0
  60.       norm = 0.0d0
  61.       k = 1
  62. c     .......... store roots isolated by balanc
  63. c                and compute matrix norm ..........
  64.       do 50 i = 1, n
  65. c
  66.          do 40 j = k, n
  67.    40    norm = norm + dabs(h(i,j))
  68. c
  69.          k = i
  70.          if (i .ge. low .and. i .le. igh) go to 50
  71.          wr(i) = h(i,i)
  72.          wi(i) = 0.0d0
  73.    50 continue
  74. c
  75.       en = igh
  76.       t = 0.0d0
  77.       itn = 30*n
  78. c     .......... search for next eigenvalues ..........
  79.    60 if (en .lt. low) go to 1001
  80.       its = 0
  81.       na = en - 1
  82.       enm2 = na - 1
  83. c     .......... look for single small sub-diagonal element
  84. c                for l=en step -1 until low do -- ..........
  85.    70 do 80 ll = low, en
  86.          l = en + low - ll
  87.          if (l .eq. low) go to 100
  88.          s = dabs(h(l-1,l-1)) + dabs(h(l,l))
  89.          if (s .eq. 0.0d0) s = norm
  90.          tst1 = s
  91.          tst2 = tst1 + dabs(h(l,l-1))
  92.          if (tst2 .eq. tst1) go to 100
  93.    80 continue
  94. c     .......... form shift ..........
  95.   100 x = h(en,en)
  96.       if (l .eq. en) go to 270
  97.       y = h(na,na)
  98.       w = h(en,na) * h(na,en)
  99.       if (l .eq. na) go to 280
  100.       if (itn .eq. 0) go to 1000
  101.       if (its .ne. 10 .and. its .ne. 20) go to 130
  102. c     .......... form exceptional shift ..........
  103.       t = t + x
  104. c
  105.       do 120 i = low, en
  106.   120 h(i,i) = h(i,i) - x
  107. c
  108.       s = dabs(h(en,na)) + dabs(h(na,enm2))
  109.       x = 0.75d0 * s
  110.       y = x
  111.       w = -0.4375d0 * s * s
  112.   130 its = its + 1
  113.       itn = itn - 1
  114. c     .......... look for two consecutive small
  115. c                sub-diagonal elements.
  116. c                for m=en-2 step -1 until l do -- ..........
  117.       do 140 mm = l, enm2
  118.          m = enm2 + l - mm
  119.          zz = h(m,m)
  120.          r = x - zz
  121.          s = y - zz
  122.          p = (r * s - w) / h(m+1,m) + h(m,m+1)
  123.          q = h(m+1,m+1) - zz - r - s
  124.          r = h(m+2,m+1)
  125.          s = dabs(p) + dabs(q) + dabs(r)
  126.          p = p / s
  127.          q = q / s
  128.          r = r / s
  129.          if (m .eq. l) go to 150
  130.          tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1)))
  131.          tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r))
  132.          if (tst2 .eq. tst1) go to 150
  133.   140 continue
  134. c
  135.   150 mp2 = m + 2
  136. c
  137.       do 160 i = mp2, en
  138.          h(i,i-2) = 0.0d0
  139.          if (i .eq. mp2) go to 160
  140.          h(i,i-3) = 0.0d0
  141.   160 continue
  142. c     .......... double qr step involving rows l to en and
  143. c                columns m to en ..........
  144.       do 260 k = m, na
  145.          notlas = k .ne. na
  146.          if (k .eq. m) go to 170
  147.          p = h(k,k-1)
  148.          q = h(k+1,k-1)
  149.          r = 0.0d0
  150.          if (notlas) r = h(k+2,k-1)
  151.          x = dabs(p) + dabs(q) + dabs(r)
  152.          if (x .eq. 0.0d0) go to 260
  153.          p = p / x
  154.          q = q / x
  155.          r = r / x
  156.   170    s = dsign(dsqrt(p*p+q*q+r*r),p)
  157.          if (k .eq. m) go to 180
  158.          h(k,k-1) = -s * x
  159.          go to 190
  160.   180    if (l .ne. m) h(k,k-1) = -h(k,k-1)
  161.   190    p = p + s
  162.          x = p / s
  163.          y = q / s
  164.          zz = r / s
  165.          q = q / p
  166.          r = r / p
  167.          if (notlas) go to 225
  168. c     .......... row modification ..........
  169.          do 200 j = k, EN
  170.             p = h(k,j) + q * h(k+1,j)
  171.             h(k,j) = h(k,j) - p * x
  172.             h(k+1,j) = h(k+1,j) - p * y
  173.   200    continue
  174. c
  175.          j = min0(en,k+3)
  176. c     .......... column modification ..........
  177.          do 210 i = L, j
  178.             p = x * h(i,k) + y * h(i,k+1)
  179.             h(i,k) = h(i,k) - p
  180.             h(i,k+1) = h(i,k+1) - p * q
  181.   210    continue
  182.          go to 255
  183.   225    continue
  184. c     .......... row modification ..........
  185.          do 230 j = k, EN
  186.             p = h(k,j) + q * h(k+1,j) + r * h(k+2,j)
  187.             h(k,j) = h(k,j) - p * x
  188.             h(k+1,j) = h(k+1,j) - p * y
  189.             h(k+2,j) = h(k+2,j) - p * zz
  190.   230    continue
  191. c
  192.          j = min0(en,k+3)
  193. c     .......... column modification ..........
  194.          do 240 i = L, j
  195.             p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2)
  196.             h(i,k) = h(i,k) - p
  197.             h(i,k+1) = h(i,k+1) - p * q
  198.             h(i,k+2) = h(i,k+2) - p * r
  199.   240    continue
  200.   255    continue
  201. c
  202.   260 continue
  203. c
  204.       go to 70
  205. c     .......... one root found ..........
  206.   270 wr(en) = x + t
  207.       wi(en) = 0.0d0
  208.       en = na
  209.       go to 60
  210. c     .......... two roots found ..........
  211.   280 p = (y - x) / 2.0d0
  212.       q = p * p + w
  213.       zz = dsqrt(dabs(q))
  214.       x = x + t
  215.       if (q .lt. 0.0d0) go to 320
  216. c     .......... real pair ..........
  217.       zz = p + dsign(zz,p)
  218.       wr(na) = x + zz
  219.       wr(en) = wr(na)
  220.       if (zz .ne. 0.0d0) wr(en) = x - w / zz
  221.       wi(na) = 0.0d0
  222.       wi(en) = 0.0d0
  223.       go to 330
  224. c     .......... complex pair ..........
  225.   320 wr(na) = x + p
  226.       wr(en) = x + p
  227.       wi(na) = zz
  228.       wi(en) = -zz
  229.   330 en = enm2
  230.       go to 60
  231. c     .......... set error -- all eigenvalues have not
  232. c                converged after 30*n iterations ..........
  233.  1000 ierr = en
  234.  1001 return
  235.       end
  236.