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