home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / comlr2.f < prev    next >
Text File  |  1996-09-28  |  12KB  |  363 lines
  1.       subroutine comlr2(nm,n,low,igh,int,hr,hi,wr,wi,zr,zi,ierr)
  2. C  MESHED overflow control WITH vectors of isolated roots (10/19/89 BSG)
  3. C  MESHED overflow control WITH triangular multiply (10/30/89 BSG)
  4. c
  5.       integer i,j,k,l,m,n,en,ii,jj,ll,mm,nm,nn,igh,im1,ip1,
  6.      x        itn,its,low,mp1,enm1,iend,ierr
  7.       double precision hr(nm,n),hi(nm,n),wr(n),wi(n),zr(nm,n),zi(nm,n)
  8.       double precision si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,norm,tst1,tst2
  9.       integer int(igh)
  10. c
  11. c     this subroutine is a translation of the algol procedure comlr2,
  12. c     num. math. 16, 181-204(1970) by peters and wilkinson.
  13. c     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
  14. c
  15. c     this subroutine finds the eigenvalues and eigenvectors
  16. c     of a complex upper hessenberg matrix by the modified lr
  17. c     method.  the eigenvectors of a complex general matrix
  18. c     can also be found if  comhes  has been used to reduce
  19. c     this general matrix to hessenberg form.
  20. c
  21. c     on input
  22. c
  23. c        nm must be set to the row dimension of two-dimensional
  24. c          array parameters as declared in the calling program
  25. c          dimension statement.
  26. c
  27. c        n is the order of the matrix.
  28. c
  29. c        low and igh are integers determined by the balancing
  30. c          subroutine  cbal.  if  cbal  has not been used,
  31. c          set low=1, igh=n.
  32. c
  33. c        int contains information on the rows and columns interchanged
  34. c          in the reduction by  comhes, if performed.  only elements
  35. c          low through igh are used.  if the eigenvectors of the hessen-
  36. c          berg matrix are desired, set int(j)=j for these elements.
  37. c
  38. c        hr and hi contain the real and imaginary parts,
  39. c          respectively, of the complex upper hessenberg matrix.
  40. c          their lower triangles below the subdiagonal contain the
  41. c          multipliers which were used in the reduction by  comhes,
  42. c          if performed.  if the eigenvectors of the hessenberg
  43. c          matrix are desired, these elements must be set to zero.
  44. c
  45. c     on output
  46. c
  47. c        the upper hessenberg portions of hr and hi have been
  48. c          destroyed, but the location hr(1,1) contains the norm
  49. c          of the triangularized matrix.
  50. c
  51. c        wr and wi contain the real and imaginary parts,
  52. c          respectively, of the eigenvalues.  if an error
  53. c          exit is made, the eigenvalues should be correct
  54. c          for indices ierr+1,...,n.
  55. c
  56. c        zr and zi contain the real and imaginary parts,
  57. c          respectively, of the eigenvectors.  the eigenvectors
  58. c          are unnormalized.  if an error exit is made, none of
  59. c          the eigenvectors has been found.
  60. c
  61. c        ierr is set to
  62. c          zero       for normal return,
  63. c          j          if the limit of 30*n iterations is exhausted
  64. c                     while the j-th eigenvalue is being sought.
  65. c
  66. c     
  67. c     calls cdiv for complex division.
  68. c     calls csroot for complex square root.
  69. c
  70. c     questions and comments should be directed to burton s. garbow,
  71. c     mathematics and computer science div, argonne national laboratory
  72. c
  73. c     this version dated october 1989.
  74. c
  75. c     ------------------------------------------------------------------
  76. c
  77.       ierr = 0
  78. c     .......... initialize eigenvector matrix ..........
  79.       do 100 i = 1, n
  80. c
  81.          do 100 j = 1, n
  82.             zr(i,j) = 0.0d0
  83.             zi(i,j) = 0.0d0
  84.             if (i .eq. j) zr(i,j) = 1.0d0
  85.   100 continue
  86. c     .......... form the matrix of accumulated transformations
  87. c                from the information left by comhes ..........
  88.       iend = igh - low - 1
  89.       if (iend .le. 0) go to 180
  90. c     .......... for i=igh-1 step -1 until low+1 do -- ..........
  91.       do 160 ii = 1, iend
  92.          i = igh - ii
  93.          ip1 = i + 1
  94. c
  95.          do 120 k = ip1, igh
  96.             zr(k,i) = hr(k,i-1)
  97.             zi(k,i) = hi(k,i-1)
  98.   120    continue
  99. c
  100.          j = int(i)
  101.          if (i .eq. j) go to 160
  102. c
  103.          do 140 k = i, igh
  104.             zr(i,k) = zr(j,k)
  105.             zi(i,k) = zi(j,k)
  106.             zr(j,k) = 0.0d0
  107.             zi(j,k) = 0.0d0
  108.   140    continue
  109. c
  110.          zr(j,i) = 1.0d0
  111.   160 continue
  112. c     .......... store roots isolated by cbal ..........
  113.   180 do 200 i = 1, n
  114.          if (i .ge. low .and. i .le. igh) go to 200
  115.          wr(i) = hr(i,i)
  116.          wi(i) = hi(i,i)
  117.   200 continue
  118. c
  119.       en = igh
  120.       tr = 0.0d0
  121.       ti = 0.0d0
  122.       itn = 30*n
  123. c     .......... search for next eigenvalue ..........
  124.   220 if (en .lt. low) go to 680
  125.       its = 0
  126.       enm1 = en - 1
  127. c     .......... look for single small sub-diagonal element
  128. c                for l=en step -1 until low do -- ..........
  129.   240 do 260 ll = low, en
  130.          l = en + low - ll
  131.          if (l .eq. low) go to 300
  132.          tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1))
  133.      x            + dabs(hr(l,l)) + dabs(hi(l,l))
  134.          tst2 = tst1 + dabs(hr(l,l-1)) + dabs(hi(l,l-1))
  135.          if (tst2 .eq. tst1) go to 300
  136.   260 continue
  137. c     .......... form shift ..........
  138.   300 if (l .eq. en) go to 660
  139.       if (itn .eq. 0) go to 1000
  140.       if (its .eq. 10 .or. its .eq. 20) go to 320
  141.       sr = hr(en,en)
  142.       si = hi(en,en)
  143.       xr = hr(enm1,en) * hr(en,enm1) - hi(enm1,en) * hi(en,enm1)
  144.       xi = hr(enm1,en) * hi(en,enm1) + hi(enm1,en) * hr(en,enm1)
  145.       if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340
  146.       yr = (hr(enm1,enm1) - sr) / 2.0d0
  147.       yi = (hi(enm1,enm1) - si) / 2.0d0
  148.       call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi)
  149.       if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310
  150.       zzr = -zzr
  151.       zzi = -zzi
  152.   310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi)
  153.       sr = sr - xr
  154.       si = si - xi
  155.       go to 340
  156. c     .......... form exceptional shift ..........
  157.   320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2))
  158.       si = dabs(hi(en,enm1)) + dabs(hi(enm1,en-2))
  159. c
  160.   340 do 360 i = low, en
  161.          hr(i,i) = hr(i,i) - sr
  162.          hi(i,i) = hi(i,i) - si
  163.   360 continue
  164. c
  165.       tr = tr + sr
  166.       ti = ti + si
  167.       its = its + 1
  168.       itn = itn - 1
  169. c     .......... look for two consecutive small
  170. c                sub-diagonal elements ..........
  171.       xr = dabs(hr(enm1,enm1)) + dabs(hi(enm1,enm1))
  172.       yr = dabs(hr(en,enm1)) + dabs(hi(en,enm1))
  173.       zzr = dabs(hr(en,en)) + dabs(hi(en,en))
  174. c     .......... for m=en-1 step -1 until l do -- ..........
  175.       do 380 mm = l, enm1
  176.          m = enm1 + l - mm
  177.          if (m .eq. l) go to 420
  178.          yi = yr
  179.          yr = dabs(hr(m,m-1)) + dabs(hi(m,m-1))
  180.          xi = zzr
  181.          zzr = xr
  182.          xr = dabs(hr(m-1,m-1)) + dabs(hi(m-1,m-1))
  183.          tst1 = zzr / yi * (zzr + xr + xi)
  184.          tst2 = tst1 + yr
  185.          if (tst2 .eq. tst1) go to 420
  186.   380 continue
  187. c     .......... triangular decomposition h=l*r ..........
  188.   420 mp1 = m + 1
  189. c
  190.       do 520 i = mp1, en
  191.          im1 = i - 1
  192.          xr = hr(im1,im1)
  193.          xi = hi(im1,im1)
  194.          yr = hr(i,im1)
  195.          yi = hi(i,im1)
  196.          if (dabs(xr) + dabs(xi) .ge. dabs(yr) + dabs(yi)) go to 460
  197. c     .......... interchange rows of hr and hi ..........
  198.          do 440 j = im1, n
  199.             zzr = hr(im1,j)
  200.             hr(im1,j) = hr(i,j)
  201.             hr(i,j) = zzr
  202.             zzi = hi(im1,j)
  203.             hi(im1,j) = hi(i,j)
  204.             hi(i,j) = zzi
  205.   440    continue
  206. c
  207.          call cdiv(xr,xi,yr,yi,zzr,zzi)
  208.          wr(i) = 1.0d0
  209.          go to 480
  210.   460    call cdiv(yr,yi,xr,xi,zzr,zzi)
  211.          wr(i) = -1.0d0
  212.   480    hr(i,im1) = zzr
  213.          hi(i,im1) = zzi
  214. c
  215.          do 500 j = i, n
  216.             hr(i,j) = hr(i,j) - zzr * hr(im1,j) + zzi * hi(im1,j)
  217.             hi(i,j) = hi(i,j) - zzr * hi(im1,j) - zzi * hr(im1,j)
  218.   500    continue
  219. c
  220.   520 continue
  221. c     .......... composition r*l=h ..........
  222.       do 640 j = mp1, en
  223.          xr = hr(j,j-1)
  224.          xi = hi(j,j-1)
  225.          hr(j,j-1) = 0.0d0
  226.          hi(j,j-1) = 0.0d0
  227. c     .......... interchange columns of hr, hi, zr, and zi,
  228. c                if necessary ..........
  229.          if (wr(j) .le. 0.0d0) go to 580
  230. c
  231.          do 540 i = 1, j
  232.             zzr = hr(i,j-1)
  233.             hr(i,j-1) = hr(i,j)
  234.             hr(i,j) = zzr
  235.             zzi = hi(i,j-1)
  236.             hi(i,j-1) = hi(i,j)
  237.             hi(i,j) = zzi
  238.   540    continue
  239. c
  240.          do 560 i = low, igh
  241.             zzr = zr(i,j-1)
  242.             zr(i,j-1) = zr(i,j)
  243.             zr(i,j) = zzr
  244.             zzi = zi(i,j-1)
  245.             zi(i,j-1) = zi(i,j)
  246.             zi(i,j) = zzi
  247.   560    continue
  248. c
  249.   580    do 600 i = 1, j
  250.             hr(i,j-1) = hr(i,j-1) + xr * hr(i,j) - xi * hi(i,j)
  251.             hi(i,j-1) = hi(i,j-1) + xr * hi(i,j) + xi * hr(i,j)
  252.   600    continue
  253. c     .......... accumulate transformations ..........
  254.          do 620 i = low, igh
  255.             zr(i,j-1) = zr(i,j-1) + xr * zr(i,j) - xi * zi(i,j)
  256.             zi(i,j-1) = zi(i,j-1) + xr * zi(i,j) + xi * zr(i,j)
  257.   620    continue
  258. c
  259.   640 continue
  260. c
  261.       go to 240
  262. c     .......... a root found ..........
  263.   660 hr(en,en) = hr(en,en) + tr
  264.       wr(en) = hr(en,en)
  265.       hi(en,en) = hi(en,en) + ti
  266.       wi(en) = hi(en,en)
  267.       en = enm1
  268.       go to 220
  269. c     .......... all roots found.  backsubstitute to find
  270. c                vectors of upper triangular form ..........
  271.   680 norm = 0.0d0
  272. c
  273.       do 720 i = 1, n
  274. c
  275.          do 720 j = i, n
  276.             tr = dabs(hr(i,j)) + dabs(hi(i,j))
  277.             if (tr .gt. norm) norm = tr
  278.   720 continue
  279. c
  280.       hr(1,1) = norm
  281.       if (n .eq. 1 .or. norm .eq. 0.0d0) go to 1001
  282. c     .......... for en=n step -1 until 2 do -- ..........
  283.       do 800 nn = 2, n
  284.          en = n + 2 - nn
  285.          xr = wr(en)
  286.          xi = wi(en)
  287.          hr(en,en) = 1.0d0
  288.          hi(en,en) = 0.0d0
  289.          enm1 = en - 1
  290. c     .......... for i=en-1 step -1 until 1 do -- ..........
  291.          do 780 ii = 1, enm1
  292.             i = en - ii
  293.             zzr = 0.0d0
  294.             zzi = 0.0d0
  295.             ip1 = i + 1
  296. c
  297.             do 740 j = ip1, en
  298.                zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en)
  299.                zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en)
  300.   740       continue
  301. c
  302.             yr = xr - wr(i)
  303.             yi = xi - wi(i)
  304.             if (yr .ne. 0.0d0 .or. yi .ne. 0.0d0) go to 765
  305.                tst1 = norm
  306.                yr = tst1
  307.   760          yr = 0.01d0 * yr
  308.                tst2 = norm + yr
  309.                if (tst2 .gt. tst1) go to 760
  310.   765       continue
  311.             call cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en))
  312. c     .......... overflow control ..........
  313.             tr = dabs(hr(i,en)) + dabs(hi(i,en))
  314.             if (tr .eq. 0.0d0) go to 780
  315.             tst1 = tr
  316.             tst2 = tst1 + 1.0d0/tst1
  317.             if (tst2 .gt. tst1) go to 780
  318.             do 770 j = i, en
  319.                hr(j,en) = hr(j,en)/tr
  320.                hi(j,en) = hi(j,en)/tr
  321.   770       continue
  322. c
  323.   780    continue
  324. c
  325.   800 continue
  326. c     .......... end backsubstitution ..........
  327. c     .......... vectors of isolated roots ..........
  328.       do  840 i = 1, N
  329.          if (i .ge. low .and. i .le. igh) go to 840
  330. c
  331.          do 820 j = I, n
  332.             zr(i,j) = hr(i,j)
  333.             zi(i,j) = hi(i,j)
  334.   820    continue
  335. c
  336.   840 continue
  337. c     .......... multiply by transformation matrix to give
  338. c                vectors of original full matrix.
  339. c                for j=n step -1 until low do -- ..........
  340.       do 880 jj = low, N
  341.          j = n + low - jj
  342.          m = min0(j,igh)
  343. c
  344.          do 880 i = low, igh
  345.             zzr = 0.0d0
  346.             zzi = 0.0d0
  347. c
  348.             do 860 k = low, m
  349.                zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j)
  350.                zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j)
  351.   860       continue
  352. c
  353.             zr(i,j) = zzr
  354.             zi(i,j) = zzi
  355.   880 continue
  356. c
  357.       go to 1001
  358. c     .......... set error -- all eigenvalues have not
  359. c                converged after 30*n iterations ..........
  360.  1000 ierr = en
  361.  1001 return
  362.       end
  363.