home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / comqr.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  223 lines
  1.       subroutine comqr(nm,n,low,igh,hr,hi,wr,wi,ierr)
  2. c
  3.       integer i,j,l,n,en,ll,nm,igh,itn,its,low,lp1,enm1,ierr
  4.       double precision hr(nm,n),hi(nm,n),wr(n),wi(n)
  5.       double precision si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,norm,tst1,tst2,
  6.      x       pythag
  7. c
  8. c     this subroutine is a translation of a unitary analogue of the
  9. c     algol procedure  comlr, num. math. 12, 369-376(1968) by martin
  10. c     and wilkinson.
  11. c     handbook for auto. comp., vol.ii-linear algebra, 396-403(1971).
  12. c     the unitary analogue substitutes the qr algorithm of francis
  13. c     (comp. jour. 4, 332-345(1962)) for the lr algorithm.
  14. c
  15. c     this subroutine finds the eigenvalues of a complex
  16. c     upper hessenberg matrix by the qr method.
  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        low and igh are integers determined by the balancing
  27. c          subroutine  cbal.  if  cbal  has not been used,
  28. c          set low=1, igh=n.
  29. c
  30. c        hr and hi contain the real and imaginary parts,
  31. c          respectively, of the complex upper hessenberg matrix.
  32. c          their lower triangles below the subdiagonal contain
  33. c          information about the unitary transformations used in
  34. c          the reduction by  corth, if performed.
  35. c
  36. c     on output
  37. c
  38. c        the upper hessenberg portions of hr and hi have been
  39. c          destroyed.  therefore, they must be saved before
  40. c          calling  comqr  if subsequent calculation of
  41. c          eigenvectors is to be performed.
  42. c
  43. c        wr and wi contain the real and imaginary parts,
  44. c          respectively, of the eigenvalues.  if an error
  45. c          exit is made, the eigenvalues should be correct
  46. c          for indices ierr+1,...,n.
  47. c
  48. c        ierr is set to
  49. c          zero       for normal return,
  50. c          j          if the limit of 30*n iterations is exhausted
  51. c                     while the j-th eigenvalue is being sought.
  52. c
  53. c     calls cdiv for complex division.
  54. c     calls csroot for complex square root.
  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.       if (low .eq. igh) go to 180
  66. c     .......... create real subdiagonal elements ..........
  67.       l = low + 1
  68. c
  69.       do 170 i = l, igh
  70.          ll = min0(i+1,igh)
  71.          if (hi(i,i-1) .eq. 0.0d0) go to 170
  72.          norm = pythag(hr(i,i-1),hi(i,i-1))
  73.          yr = hr(i,i-1) / norm
  74.          yi = hi(i,i-1) / norm
  75.          hr(i,i-1) = norm
  76.          hi(i,i-1) = 0.0d0
  77. c
  78.          do 155 j = i, igh
  79.             si = yr * hi(i,j) - yi * hr(i,j)
  80.             hr(i,j) = yr * hr(i,j) + yi * hi(i,j)
  81.             hi(i,j) = si
  82.   155    continue
  83. c
  84.          do 160 j = low, ll
  85.             si = yr * hi(j,i) + yi * hr(j,i)
  86.             hr(j,i) = yr * hr(j,i) - yi * hi(j,i)
  87.             hi(j,i) = si
  88.   160    continue
  89. c
  90.   170 continue
  91. c     .......... store roots isolated by cbal ..........
  92.   180 do 200 i = 1, n
  93.          if (i .ge. low .and. i .le. igh) go to 200
  94.          wr(i) = hr(i,i)
  95.          wi(i) = hi(i,i)
  96.   200 continue
  97. c
  98.       en = igh
  99.       tr = 0.0d0
  100.       ti = 0.0d0
  101.       itn = 30*n
  102. c     .......... search for next eigenvalue ..........
  103.   220 if (en .lt. low) go to 1001
  104.       its = 0
  105.       enm1 = en - 1
  106. c     .......... look for single small sub-diagonal element
  107. c                for l=en step -1 until low d0 -- ..........
  108.   240 do 260 ll = low, en
  109.          l = en + low - ll
  110.          if (l .eq. low) go to 300
  111.          tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1))
  112.      x            + dabs(hr(l,l)) + dabs(hi(l,l))
  113.          tst2 = tst1 + dabs(hr(l,l-1))
  114.          if (tst2 .eq. tst1) go to 300
  115.   260 continue
  116. c     .......... form shift ..........
  117.   300 if (l .eq. en) go to 660
  118.       if (itn .eq. 0) go to 1000
  119.       if (its .eq. 10 .or. its .eq. 20) go to 320
  120.       sr = hr(en,en)
  121.       si = hi(en,en)
  122.       xr = hr(enm1,en) * hr(en,enm1)
  123.       xi = hi(enm1,en) * hr(en,enm1)
  124.       if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340
  125.       yr = (hr(enm1,enm1) - sr) / 2.0d0
  126.       yi = (hi(enm1,enm1) - si) / 2.0d0
  127.       call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi)
  128.       if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310
  129.       zzr = -zzr
  130.       zzi = -zzi
  131.   310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi)
  132.       sr = sr - xr
  133.       si = si - xi
  134.       go to 340
  135. c     .......... form exceptional shift ..........
  136.   320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2))
  137.       si = 0.0d0
  138. c
  139.   340 do 360 i = low, en
  140.          hr(i,i) = hr(i,i) - sr
  141.          hi(i,i) = hi(i,i) - si
  142.   360 continue
  143. c
  144.       tr = tr + sr
  145.       ti = ti + si
  146.       its = its + 1
  147.       itn = itn - 1
  148. c     .......... reduce to triangle (rows) ..........
  149.       lp1 = l + 1
  150. c
  151.       do 500 i = lp1, en
  152.          sr = hr(i,i-1)
  153.          hr(i,i-1) = 0.0d0
  154.          norm = pythag(pythag(hr(i-1,i-1),hi(i-1,i-1)),sr)
  155.          xr = hr(i-1,i-1) / norm
  156.          wr(i-1) = xr
  157.          xi = hi(i-1,i-1) / norm
  158.          wi(i-1) = xi
  159.          hr(i-1,i-1) = norm
  160.          hi(i-1,i-1) = 0.0d0
  161.          hi(i,i-1) = sr / norm
  162. c
  163.          do 490 j = i, en
  164.             yr = hr(i-1,j)
  165.             yi = hi(i-1,j)
  166.             zzr = hr(i,j)
  167.             zzi = hi(i,j)
  168.             hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr
  169.             hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi
  170.             hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr
  171.             hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi
  172.   490    continue
  173. c
  174.   500 continue
  175. c
  176.       si = hi(en,en)
  177.       if (si .eq. 0.0d0) go to 540
  178.       norm = pythag(hr(en,en),si)
  179.       sr = hr(en,en) / norm
  180.       si = si / norm
  181.       hr(en,en) = norm
  182.       hi(en,en) = 0.0d0
  183. c     .......... inverse operation (columns) ..........
  184.   540 do 600 j = lp1, en
  185.          xr = wr(j-1)
  186.          xi = wi(j-1)
  187. c
  188.          do 580 i = l, j
  189.             yr = hr(i,j-1)
  190.             yi = 0.0d0
  191.             zzr = hr(i,j)
  192.             zzi = hi(i,j)
  193.             if (i .eq. j) go to 560
  194.             yi = hi(i,j-1)
  195.             hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi
  196.   560       hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr
  197.             hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr
  198.             hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi
  199.   580    continue
  200. c
  201.   600 continue
  202. c
  203.       if (si .eq. 0.0d0) go to 240
  204. c
  205.       do 630 i = l, en
  206.          yr = hr(i,en)
  207.          yi = hi(i,en)
  208.          hr(i,en) = sr * yr - si * yi
  209.          hi(i,en) = sr * yi + si * yr
  210.   630 continue
  211. c
  212.       go to 240
  213. c     .......... a root found ..........
  214.   660 wr(en) = hr(en,en) + tr
  215.       wi(en) = hi(en,en) + ti
  216.       en = enm1
  217.       go to 220
  218. c     .......... set error -- all eigenvalues have not
  219. c                converged after 30*n iterations ..........
  220.  1000 ierr = en
  221.  1001 return
  222.       end
  223.