home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / comlr.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  206 lines
  1.       subroutine comlr(nm,n,low,igh,hr,hi,wr,wi,ierr)
  2. c
  3.       integer i,j,l,m,n,en,ll,mm,nm,igh,im1,itn,its,low,mp1,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,tst1,tst2
  6. c
  7. c     this subroutine is a translation of the algol procedure comlr,
  8. c     num. math. 12, 369-376(1968) by martin and wilkinson.
  9. c     handbook for auto. comp., vol.ii-linear algebra, 396-403(1971).
  10. c
  11. c     this subroutine finds the eigenvalues of a complex
  12. c     upper hessenberg matrix by the modified lr method.
  13. c
  14. c     on input
  15. c
  16. c        nm must be set to the row dimension of two-dimensional
  17. c          array parameters as declared in the calling program
  18. c          dimension statement.
  19. c
  20. c        n is the order of the matrix.
  21. c
  22. c        low and igh are integers determined by the balancing
  23. c          subroutine  cbal.  if  cbal  has not been used,
  24. c          set low=1, igh=n.
  25. c
  26. c        hr and hi contain the real and imaginary parts,
  27. c          respectively, of the complex upper hessenberg matrix.
  28. c          their lower triangles below the subdiagonal contain the
  29. c          multipliers which were used in the reduction by  comhes,
  30. c          if performed.
  31. c
  32. c     on output
  33. c
  34. c        the upper hessenberg portions of hr and hi have been
  35. c          destroyed.  therefore, they must be saved before
  36. c          calling  comlr  if subsequent calculation of
  37. c          eigenvectors is to be performed.
  38. c
  39. c        wr and wi contain the real and imaginary parts,
  40. c          respectively, of the eigenvalues.  if an error
  41. c          exit is made, the eigenvalues should be correct
  42. c          for indices ierr+1,...,n.
  43. c
  44. c        ierr is set to
  45. c          zero       for normal return,
  46. c          j          if the limit of 30*n iterations is exhausted
  47. c                     while the j-th eigenvalue is being sought.
  48. c
  49. c     calls cdiv for complex division.
  50. c     calls csroot for complex square root.
  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 august 1983.
  56. c
  57. c     ------------------------------------------------------------------
  58. c
  59.       ierr = 0
  60. c     .......... store roots isolated by cbal ..........
  61.       do 200 i = 1, n
  62.          if (i .ge. low .and. i .le. igh) go to 200
  63.          wr(i) = hr(i,i)
  64.          wi(i) = hi(i,i)
  65.   200 continue
  66. c
  67.       en = igh
  68.       tr = 0.0d0
  69.       ti = 0.0d0
  70.       itn = 30*n
  71. c     .......... search for next eigenvalue ..........
  72.   220 if (en .lt. low) go to 1001
  73.       its = 0
  74.       enm1 = en - 1
  75. c     .......... look for single small sub-diagonal element
  76. c                for l=en step -1 until low d0 -- ..........
  77.   240 do 260 ll = low, en
  78.          l = en + low - ll
  79.          if (l .eq. low) go to 300
  80.          tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1))
  81.      x            + dabs(hr(l,l)) + dabs(hi(l,l))
  82.          tst2 = tst1 + dabs(hr(l,l-1)) + dabs(hi(l,l-1))
  83.          if (tst2 .eq. tst1) go to 300
  84.   260 continue
  85. c     .......... form shift ..........
  86.   300 if (l .eq. en) go to 660
  87.       if (itn .eq. 0) go to 1000
  88.       if (its .eq. 10 .or. its .eq. 20) go to 320
  89.       sr = hr(en,en)
  90.       si = hi(en,en)
  91.       xr = hr(enm1,en) * hr(en,enm1) - hi(enm1,en) * hi(en,enm1)
  92.       xi = hr(enm1,en) * hi(en,enm1) + hi(enm1,en) * hr(en,enm1)
  93.       if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340
  94.       yr = (hr(enm1,enm1) - sr) / 2.0d0
  95.       yi = (hi(enm1,enm1) - si) / 2.0d0
  96.       call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi)
  97.       if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310
  98.       zzr = -zzr
  99.       zzi = -zzi
  100.   310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi)
  101.       sr = sr - xr
  102.       si = si - xi
  103.       go to 340
  104. c     .......... form exceptional shift ..........
  105.   320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2))
  106.       si = dabs(hi(en,enm1)) + dabs(hi(enm1,en-2))
  107. c
  108.   340 do 360 i = low, en
  109.          hr(i,i) = hr(i,i) - sr
  110.          hi(i,i) = hi(i,i) - si
  111.   360 continue
  112. c
  113.       tr = tr + sr
  114.       ti = ti + si
  115.       its = its + 1
  116.       itn = itn - 1
  117. c     .......... look for two consecutive small
  118. c                sub-diagonal elements ..........
  119.       xr = dabs(hr(enm1,enm1)) + dabs(hi(enm1,enm1))
  120.       yr = dabs(hr(en,enm1)) + dabs(hi(en,enm1))
  121.       zzr = dabs(hr(en,en)) + dabs(hi(en,en))
  122. c     .......... for m=en-1 step -1 until l do -- ..........
  123.       do 380 mm = l, enm1
  124.          m = enm1 + l - mm
  125.          if (m .eq. l) go to 420
  126.          yi = yr
  127.          yr = dabs(hr(m,m-1)) + dabs(hi(m,m-1))
  128.          xi = zzr
  129.          zzr = xr
  130.          xr = dabs(hr(m-1,m-1)) + dabs(hi(m-1,m-1))
  131.          tst1 = zzr / yi * (zzr + xr + xi)
  132.          tst2 = tst1 + yr
  133.          if (tst2 .eq. tst1) go to 420
  134.   380 continue
  135. c     .......... triangular decomposition h=l*r ..........
  136.   420 mp1 = m + 1
  137. c
  138.       do 520 i = mp1, en
  139.          im1 = i - 1
  140.          xr = hr(im1,im1)
  141.          xi = hi(im1,im1)
  142.          yr = hr(i,im1)
  143.          yi = hi(i,im1)
  144.          if (dabs(xr) + dabs(xi) .ge. dabs(yr) + dabs(yi)) go to 460
  145. c     .......... interchange rows of hr and hi ..........
  146.          do 440 j = im1, en
  147.             zzr = hr(im1,j)
  148.             hr(im1,j) = hr(i,j)
  149.             hr(i,j) = zzr
  150.             zzi = hi(im1,j)
  151.             hi(im1,j) = hi(i,j)
  152.             hi(i,j) = zzi
  153.   440    continue
  154. c
  155.          call cdiv(xr,xi,yr,yi,zzr,zzi)
  156.          wr(i) = 1.0d0
  157.          go to 480
  158.   460    call cdiv(yr,yi,xr,xi,zzr,zzi)
  159.          wr(i) = -1.0d0
  160.   480    hr(i,im1) = zzr
  161.          hi(i,im1) = zzi
  162. c
  163.          do 500 j = i, en
  164.             hr(i,j) = hr(i,j) - zzr * hr(im1,j) + zzi * hi(im1,j)
  165.             hi(i,j) = hi(i,j) - zzr * hi(im1,j) - zzi * hr(im1,j)
  166.   500    continue
  167. c
  168.   520 continue
  169. c     .......... composition r*l=h ..........
  170.       do 640 j = mp1, en
  171.          xr = hr(j,j-1)
  172.          xi = hi(j,j-1)
  173.          hr(j,j-1) = 0.0d0
  174.          hi(j,j-1) = 0.0d0
  175. c     .......... interchange columns of hr and hi,
  176. c                if necessary ..........
  177.          if (wr(j) .le. 0.0d0) go to 580
  178. c
  179.          do 540 i = l, j
  180.             zzr = hr(i,j-1)
  181.             hr(i,j-1) = hr(i,j)
  182.             hr(i,j) = zzr
  183.             zzi = hi(i,j-1)
  184.             hi(i,j-1) = hi(i,j)
  185.             hi(i,j) = zzi
  186.   540    continue
  187. c
  188.   580    do 600 i = l, j
  189.             hr(i,j-1) = hr(i,j-1) + xr * hr(i,j) - xi * hi(i,j)
  190.             hi(i,j-1) = hi(i,j-1) + xr * hi(i,j) + xi * hr(i,j)
  191.   600    continue
  192. c
  193.   640 continue
  194. c
  195.       go to 240
  196. c     .......... a root found ..........
  197.   660 wr(en) = hr(en,en) + tr
  198.       wi(en) = hi(en,en) + ti
  199.       en = enm1
  200.       go to 220
  201. c     .......... set error -- all eigenvalues have not
  202. c                converged after 30*n iterations ..........
  203.  1000 ierr = en
  204.  1001 return
  205.       end
  206.