home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / hqr2.f < prev    next >
Text File  |  1996-09-28  |  14KB  |  450 lines
  1.       subroutine hqr2(nm,n,low,igh,h,wr,wi,z,ierr)
  2. c
  3.       integer i,j,k,l,m,n,en,ii,jj,ll,mm,na,nm,nn,
  4.      x        igh,itn,its,low,mp2,enm2,ierr
  5.       double precision h(nm,n),wr(n),wi(n),z(nm,n)
  6.       double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2
  7.       logical notlas
  8. c
  9. c     this subroutine is a translation of the algol procedure hqr2,
  10. c     num. math. 16, 181-204(1970) by peters and wilkinson.
  11. c     handbook for auto. comp., vol.ii-linear algebra, 372-395(1971).
  12. c
  13. c     this subroutine finds the eigenvalues and eigenvectors
  14. c     of a real upper hessenberg matrix by the qr method.  the
  15. c     eigenvectors of a real general matrix can also be found
  16. c     if  elmhes  and  eltran  or  orthes  and  ortran  have
  17. c     been used to reduce this general matrix to hessenberg form
  18. c     and to accumulate the similarity transformations.
  19. c
  20. c     on input
  21. c
  22. c        nm must be set to the row dimension of two-dimensional
  23. c          array parameters as declared in the calling program
  24. c          dimension statement.
  25. c
  26. c        n is the order of the matrix.
  27. c
  28. c        low and igh are integers determined by the balancing
  29. c          subroutine  balanc.  if  balanc  has not been used,
  30. c          set low=1, igh=n.
  31. c
  32. c        h contains the upper hessenberg matrix.
  33. c
  34. c        z contains the transformation matrix produced by  eltran
  35. c          after the reduction by  elmhes, or by  ortran  after the
  36. c          reduction by  orthes, if performed.  if the eigenvectors
  37. c          of the hessenberg matrix are desired, z must contain the
  38. c          identity matrix.
  39. c
  40. c     on output
  41. c
  42. c        h has been destroyed.
  43. c
  44. c        wr and wi contain the real and imaginary parts,
  45. c          respectively, of the eigenvalues.  the eigenvalues
  46. c          are unordered except that complex conjugate pairs
  47. c          of values appear consecutively with the eigenvalue
  48. c          having the positive imaginary part first.  if an
  49. c          error exit is made, the eigenvalues should be correct
  50. c          for indices ierr+1,...,n.
  51. c
  52. c        z contains the real and imaginary parts of the eigenvectors.
  53. c          if the i-th eigenvalue is real, the i-th column of z
  54. c          contains its eigenvector.  if the i-th eigenvalue is complex
  55. c          with positive imaginary part, the i-th and (i+1)-th
  56. c          columns of z contain the real and imaginary parts of its
  57. c          eigenvector.  the eigenvectors are unnormalized.  if an
  58. c          error exit is made, none of the eigenvectors has been found.
  59. c
  60. c        ierr is set to
  61. c          zero       for normal return,
  62. c          j          if the limit of 30*n iterations is exhausted
  63. c                     while the j-th eigenvalue is being sought.
  64. c
  65. c     calls cdiv for complex division.
  66. c
  67. c     questions and comments should be directed to burton s. garbow,
  68. c     mathematics and computer science div, argonne national laboratory
  69. c
  70. c     this version dated august 1983.
  71. c
  72. c     ------------------------------------------------------------------
  73. c
  74.       ierr = 0
  75.       norm = 0.0d0
  76.       k = 1
  77. c     .......... store roots isolated by balanc
  78. c                and compute matrix norm ..........
  79.       do 50 i = 1, n
  80. c
  81.          do 40 j = k, n
  82.    40    norm = norm + dabs(h(i,j))
  83. c
  84.          k = i
  85.          if (i .ge. low .and. i .le. igh) go to 50
  86.          wr(i) = h(i,i)
  87.          wi(i) = 0.0d0
  88.    50 continue
  89. c
  90.       en = igh
  91.       t = 0.0d0
  92.       itn = 30*n
  93. c     .......... search for next eigenvalues ..........
  94.    60 if (en .lt. low) go to 340
  95.       its = 0
  96.       na = en - 1
  97.       enm2 = na - 1
  98. c     .......... look for single small sub-diagonal element
  99. c                for l=en step -1 until low do -- ..........
  100.    70 do 80 ll = low, en
  101.          l = en + low - ll
  102.          if (l .eq. low) go to 100
  103.          s = dabs(h(l-1,l-1)) + dabs(h(l,l))
  104.          if (s .eq. 0.0d0) s = norm
  105.          tst1 = s
  106.          tst2 = tst1 + dabs(h(l,l-1))
  107.          if (tst2 .eq. tst1) go to 100
  108.    80 continue
  109. c     .......... form shift ..........
  110.   100 x = h(en,en)
  111.       if (l .eq. en) go to 270
  112.       y = h(na,na)
  113.       w = h(en,na) * h(na,en)
  114.       if (l .eq. na) go to 280
  115.       if (itn .eq. 0) go to 1000
  116.       if (its .ne. 10 .and. its .ne. 20) go to 130
  117. c     .......... form exceptional shift ..........
  118.       t = t + x
  119. c
  120.       do 120 i = low, en
  121.   120 h(i,i) = h(i,i) - x
  122. c
  123.       s = dabs(h(en,na)) + dabs(h(na,enm2))
  124.       x = 0.75d0 * s
  125.       y = x
  126.       w = -0.4375d0 * s * s
  127.   130 its = its + 1
  128.       itn = itn - 1
  129. c     .......... look for two consecutive small
  130. c                sub-diagonal elements.
  131. c                for m=en-2 step -1 until l do -- ..........
  132.       do 140 mm = l, enm2
  133.          m = enm2 + l - mm
  134.          zz = h(m,m)
  135.          r = x - zz
  136.          s = y - zz
  137.          p = (r * s - w) / h(m+1,m) + h(m,m+1)
  138.          q = h(m+1,m+1) - zz - r - s
  139.          r = h(m+2,m+1)
  140.          s = dabs(p) + dabs(q) + dabs(r)
  141.          p = p / s
  142.          q = q / s
  143.          r = r / s
  144.          if (m .eq. l) go to 150
  145.          tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1)))
  146.          tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r))
  147.          if (tst2 .eq. tst1) go to 150
  148.   140 continue
  149. c
  150.   150 mp2 = m + 2
  151. c
  152.       do 160 i = mp2, en
  153.          h(i,i-2) = 0.0d0
  154.          if (i .eq. mp2) go to 160
  155.          h(i,i-3) = 0.0d0
  156.   160 continue
  157. c     .......... double qr step involving rows l to en and
  158. c                columns m to en ..........
  159.       do 260 k = m, na
  160.          notlas = k .ne. na
  161.          if (k .eq. m) go to 170
  162.          p = h(k,k-1)
  163.          q = h(k+1,k-1)
  164.          r = 0.0d0
  165.          if (notlas) r = h(k+2,k-1)
  166.          x = dabs(p) + dabs(q) + dabs(r)
  167.          if (x .eq. 0.0d0) go to 260
  168.          p = p / x
  169.          q = q / x
  170.          r = r / x
  171.   170    s = dsign(dsqrt(p*p+q*q+r*r),p)
  172.          if (k .eq. m) go to 180
  173.          h(k,k-1) = -s * x
  174.          go to 190
  175.   180    if (l .ne. m) h(k,k-1) = -h(k,k-1)
  176.   190    p = p + s
  177.          x = p / s
  178.          y = q / s
  179.          zz = r / s
  180.          q = q / p
  181.          r = r / p
  182.          if (notlas) go to 225
  183. c     .......... row modification ..........
  184.          do 200 j = k, n
  185.             p = h(k,j) + q * h(k+1,j)
  186.             h(k,j) = h(k,j) - p * x
  187.             h(k+1,j) = h(k+1,j) - p * y
  188.   200    continue
  189. c
  190.          j = min0(en,k+3)
  191. c     .......... column modification ..........
  192.          do 210 i = 1, j
  193.             p = x * h(i,k) + y * h(i,k+1)
  194.             h(i,k) = h(i,k) - p
  195.             h(i,k+1) = h(i,k+1) - p * q
  196.   210    continue
  197. c     .......... accumulate transformations ..........
  198.          do 220 i = low, igh
  199.             p = x * z(i,k) + y * z(i,k+1)
  200.             z(i,k) = z(i,k) - p
  201.             z(i,k+1) = z(i,k+1) - p * q
  202.   220    continue
  203.          go to 255
  204.   225    continue
  205. c     .......... row modification ..........
  206.          do 230 j = k, n
  207.             p = h(k,j) + q * h(k+1,j) + r * h(k+2,j)
  208.             h(k,j) = h(k,j) - p * x
  209.             h(k+1,j) = h(k+1,j) - p * y
  210.             h(k+2,j) = h(k+2,j) - p * zz
  211.   230    continue
  212. c
  213.          j = min0(en,k+3)
  214. c     .......... column modification ..........
  215.          do 240 i = 1, j
  216.             p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2)
  217.             h(i,k) = h(i,k) - p
  218.             h(i,k+1) = h(i,k+1) - p * q
  219.             h(i,k+2) = h(i,k+2) - p * r
  220.   240    continue
  221. c     .......... accumulate transformations ..........
  222.          do 250 i = low, igh
  223.             p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2)
  224.             z(i,k) = z(i,k) - p
  225.             z(i,k+1) = z(i,k+1) - p * q
  226.             z(i,k+2) = z(i,k+2) - p * r
  227.   250    continue
  228.   255    continue
  229. c
  230.   260 continue
  231. c
  232.       go to 70
  233. c     .......... one root found ..........
  234.   270 h(en,en) = x + t
  235.       wr(en) = h(en,en)
  236.       wi(en) = 0.0d0
  237.       en = na
  238.       go to 60
  239. c     .......... two roots found ..........
  240.   280 p = (y - x) / 2.0d0
  241.       q = p * p + w
  242.       zz = dsqrt(dabs(q))
  243.       h(en,en) = x + t
  244.       x = h(en,en)
  245.       h(na,na) = y + t
  246.       if (q .lt. 0.0d0) go to 320
  247. c     .......... real pair ..........
  248.       zz = p + dsign(zz,p)
  249.       wr(na) = x + zz
  250.       wr(en) = wr(na)
  251.       if (zz .ne. 0.0d0) wr(en) = x - w / zz
  252.       wi(na) = 0.0d0
  253.       wi(en) = 0.0d0
  254.       x = h(en,na)
  255.       s = dabs(x) + dabs(zz)
  256.       p = x / s
  257.       q = zz / s
  258.       r = dsqrt(p*p+q*q)
  259.       p = p / r
  260.       q = q / r
  261. c     .......... row modification ..........
  262.       do 290 j = na, n
  263.          zz = h(na,j)
  264.          h(na,j) = q * zz + p * h(en,j)
  265.          h(en,j) = q * h(en,j) - p * zz
  266.   290 continue
  267. c     .......... column modification ..........
  268.       do 300 i = 1, en
  269.          zz = h(i,na)
  270.          h(i,na) = q * zz + p * h(i,en)
  271.          h(i,en) = q * h(i,en) - p * zz
  272.   300 continue
  273. c     .......... accumulate transformations ..........
  274.       do 310 i = low, igh
  275.          zz = z(i,na)
  276.          z(i,na) = q * zz + p * z(i,en)
  277.          z(i,en) = q * z(i,en) - p * zz
  278.   310 continue
  279. c
  280.       go to 330
  281. c     .......... complex pair ..........
  282.   320 wr(na) = x + p
  283.       wr(en) = x + p
  284.       wi(na) = zz
  285.       wi(en) = -zz
  286.   330 en = enm2
  287.       go to 60
  288. c     .......... all roots found.  backsubstitute to find
  289. c                vectors of upper triangular form ..........
  290.   340 if (norm .eq. 0.0d0) go to 1001
  291. c     .......... for en=n step -1 until 1 do -- ..........
  292.       do 800 nn = 1, n
  293.          en = n + 1 - nn
  294.          p = wr(en)
  295.          q = wi(en)
  296.          na = en - 1
  297.          if (q) 710, 600, 800
  298. c     .......... real vector ..........
  299.   600    m = en
  300.          h(en,en) = 1.0d0
  301.          if (na .eq. 0) go to 800
  302. c     .......... for i=en-1 step -1 until 1 do -- ..........
  303.          do 700 ii = 1, na
  304.             i = en - ii
  305.             w = h(i,i) - p
  306.             r = 0.0d0
  307. c
  308.             do 610 j = m, en
  309.   610       r = r + h(i,j) * h(j,en)
  310. c
  311.             if (wi(i) .ge. 0.0d0) go to 630
  312.             zz = w
  313.             s = r
  314.             go to 700
  315.   630       m = i
  316.             if (wi(i) .ne. 0.0d0) go to 640
  317.             t = w
  318.             if (t .ne. 0.0d0) go to 635
  319.                tst1 = norm
  320.                t = tst1
  321.   632          t = 0.01d0 * t
  322.                tst2 = norm + t
  323.                if (tst2 .gt. tst1) go to 632
  324.   635       h(i,en) = -r / t
  325.             go to 680
  326. c     .......... solve real equations ..........
  327.   640       x = h(i,i+1)
  328.             y = h(i+1,i)
  329.             q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i)
  330.             t = (x * s - zz * r) / q
  331.             h(i,en) = t
  332.             if (dabs(x) .le. dabs(zz)) go to 650
  333.             h(i+1,en) = (-r - w * t) / x
  334.             go to 680
  335.   650       h(i+1,en) = (-s - y * t) / zz
  336. c
  337. c     .......... overflow control ..........
  338.   680       t = dabs(h(i,en))
  339.             if (t .eq. 0.0d0) go to 700
  340.             tst1 = t
  341.             tst2 = tst1 + 1.0d0/tst1
  342.             if (tst2 .gt. tst1) go to 700
  343.             do 690 j = i, en
  344.                h(j,en) = h(j,en)/t
  345.   690       continue
  346. c
  347.   700    continue
  348. c     .......... end real vector ..........
  349.          go to 800
  350. c     .......... complex vector ..........
  351.   710    m = na
  352. c     .......... last vector component chosen imaginary so that
  353. c                eigenvector matrix is triangular ..........
  354.          if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720
  355.          h(na,na) = q / h(en,na)
  356.          h(na,en) = -(h(en,en) - p) / h(en,na)
  357.          go to 730
  358.   720    call cdiv(0.0d0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en))
  359.   730    h(en,na) = 0.0d0
  360.          h(en,en) = 1.0d0
  361.          enm2 = na - 1
  362.          if (enm2 .eq. 0) go to 800
  363. c     .......... for i=en-2 step -1 until 1 do -- ..........
  364.          do 795 ii = 1, enm2
  365.             i = na - ii
  366.             w = h(i,i) - p
  367.             ra = 0.0d0
  368.             sa = 0.0d0
  369. c
  370.             do 760 j = m, en
  371.                ra = ra + h(i,j) * h(j,na)
  372.                sa = sa + h(i,j) * h(j,en)
  373.   760       continue
  374. c
  375.             if (wi(i) .ge. 0.0d0) go to 770
  376.             zz = w
  377.             r = ra
  378.             s = sa
  379.             go to 795
  380.   770       m = i
  381.             if (wi(i) .ne. 0.0d0) go to 780
  382.             call cdiv(-ra,-sa,w,q,h(i,na),h(i,en))
  383.             go to 790
  384. c     .......... solve complex equations ..........
  385.   780       x = h(i,i+1)
  386.             y = h(i+1,i)
  387.             vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q
  388.             vi = (wr(i) - p) * 2.0d0 * q
  389.             if (vr .ne. 0.0d0 .or. vi .ne. 0.0d0) go to 784
  390.                tst1 = norm * (dabs(w) + dabs(q) + dabs(x)
  391.      x                      + dabs(y) + dabs(zz))
  392.                vr = tst1
  393.   783          vr = 0.01d0 * vr
  394.                tst2 = tst1 + vr
  395.                if (tst2 .gt. tst1) go to 783
  396.   784       call cdiv(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi,
  397.      x                h(i,na),h(i,en))
  398.             if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785
  399.             h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x
  400.             h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x
  401.             go to 790
  402.   785       call cdiv(-r-y*h(i,na),-s-y*h(i,en),zz,q,
  403.      x                h(i+1,na),h(i+1,en))
  404. c
  405. c     .......... overflow control ..........
  406.   790       t = dmax1(dabs(h(i,na)), dabs(h(i,en)))
  407.             if (t .eq. 0.0d0) go to 795
  408.             tst1 = t
  409.             tst2 = tst1 + 1.0d0/tst1
  410.             if (tst2 .gt. tst1) go to 795
  411.             do 792 j = i, en
  412.                h(j,na) = h(j,na)/t
  413.                h(j,en) = h(j,en)/t
  414.   792       continue
  415. c
  416.   795    continue
  417. c     .......... end complex vector ..........
  418.   800 continue
  419. c     .......... end back substitution.
  420. c                vectors of isolated roots ..........
  421.       do 840 i = 1, n
  422.          if (i .ge. low .and. i .le. igh) go to 840
  423. c
  424.          do 820 j = i, n
  425.   820    z(i,j) = h(i,j)
  426. c
  427.   840 continue
  428. c     .......... multiply by transformation matrix to give
  429. c                vectors of original full matrix.
  430. c                for j=n step -1 until low do -- ..........
  431.       do 880 jj = low, n
  432.          j = n + low - jj
  433.          m = min0(j,igh)
  434. c
  435.          do 880 i = low, igh
  436.             zz = 0.0d0
  437. c
  438.             do 860 k = low, m
  439.   860       zz = zz + z(i,k) * h(k,j)
  440. c
  441.             z(i,j) = zz
  442.   880 continue
  443. c
  444.       go to 1001
  445. c     .......... set error -- all eigenvalues have not
  446. c                converged after 30*n iterations ..........
  447.  1000 ierr = en
  448.  1001 return
  449.       end
  450.