home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / qzit.f < prev    next >
Text File  |  1996-09-28  |  12KB  |  355 lines
  1.       subroutine qzit(nm,n,a,b,eps1,matz,z,ierr)
  2. c
  3.       integer i,j,k,l,n,en,k1,k2,ld,ll,l1,na,nm,ish,itn,its,km1,lm1,
  4.      x        enm2,ierr,lor1,enorn
  5.       double precision a(nm,n),b(nm,n),z(nm,n)
  6.       double precision r,s,t,a1,a2,a3,ep,sh,u1,u2,u3,v1,v2,v3,ani,a11,
  7.      x       a12,a21,a22,a33,a34,a43,a44,bni,b11,b12,b22,b33,b34,
  8.      x       b44,epsa,epsb,eps1,anorm,bnorm,epslon
  9.       logical matz,notlas
  10. c
  11. c     this subroutine is the second step of the qz algorithm
  12. c     for solving generalized matrix eigenvalue problems,
  13. c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart,
  14. c     as modified in technical note nasa tn d-7305(1973) by ward.
  15. c
  16. c     this subroutine accepts a pair of real matrices, one of them
  17. c     in upper hessenberg form and the other in upper triangular form.
  18. c     it reduces the hessenberg matrix to quasi-triangular form using
  19. c     orthogonal transformations while maintaining the triangular form
  20. c     of the other matrix.  it is usually preceded by  qzhes  and
  21. c     followed by  qzval  and, possibly,  qzvec.
  22. c
  23. c     on input
  24. c
  25. c        nm must be set to the row dimension of two-dimensional
  26. c          array parameters as declared in the calling program
  27. c          dimension statement.
  28. c
  29. c        n is the order of the matrices.
  30. c
  31. c        a contains a real upper hessenberg matrix.
  32. c
  33. c        b contains a real upper triangular matrix.
  34. c
  35. c        eps1 is a tolerance used to determine negligible elements.
  36. c          eps1 = 0.0 (or negative) may be input, in which case an
  37. c          element will be neglected only if it is less than roundoff
  38. c          error times the norm of its matrix.  if the input eps1 is
  39. c          positive, then an element will be considered negligible
  40. c          if it is less than eps1 times the norm of its matrix.  a
  41. c          positive value of eps1 may result in faster execution,
  42. c          but less accurate results.
  43. c
  44. c        matz should be set to .true. if the right hand transformations
  45. c          are to be accumulated for later use in computing
  46. c          eigenvectors, and to .false. otherwise.
  47. c
  48. c        z contains, if matz has been set to .true., the
  49. c          transformation matrix produced in the reduction
  50. c          by  qzhes, if performed, or else the identity matrix.
  51. c          if matz has been set to .false., z is not referenced.
  52. c
  53. c     on output
  54. c
  55. c        a has been reduced to quasi-triangular form.  the elements
  56. c          below the first subdiagonal are still zero and no two
  57. c          consecutive subdiagonal elements are nonzero.
  58. c
  59. c        b is still in upper triangular form, although its elements
  60. c          have been altered.  the location b(n,1) is used to store
  61. c          eps1 times the norm of b for later use by  qzval  and  qzvec.
  62. c
  63. c        z contains the product of the right hand transformations
  64. c          (for both steps) if matz has been set to .true..
  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     questions and comments should be directed to burton s. garbow,
  72. c     mathematics and computer science div, argonne national laboratory
  73. c
  74. c     this version dated august 1983.
  75. c
  76. c     ------------------------------------------------------------------
  77. c
  78.       ierr = 0
  79. c     .......... compute epsa,epsb ..........
  80.       anorm = 0.0d0
  81.       bnorm = 0.0d0
  82. c
  83.       do 30 i = 1, n
  84.          ani = 0.0d0
  85.          if (i .ne. 1) ani = dabs(a(i,i-1))
  86.          bni = 0.0d0
  87. c
  88.          do 20 j = i, n
  89.             ani = ani + dabs(a(i,j))
  90.             bni = bni + dabs(b(i,j))
  91.    20    continue
  92. c
  93.          if (ani .gt. anorm) anorm = ani
  94.          if (bni .gt. bnorm) bnorm = bni
  95.    30 continue
  96. c
  97.       if (anorm .eq. 0.0d0) anorm = 1.0d0
  98.       if (bnorm .eq. 0.0d0) bnorm = 1.0d0
  99.       ep = eps1
  100.       if (ep .gt. 0.0d0) go to 50
  101. c     .......... use roundoff level if eps1 is zero ..........
  102.       ep = epslon(1.0d0)
  103.    50 epsa = ep * anorm
  104.       epsb = ep * bnorm
  105. c     .......... reduce a to quasi-triangular form, while
  106. c                keeping b triangular ..........
  107.       lor1 = 1
  108.       enorn = n
  109.       en = n
  110.       itn = 30*n
  111. c     .......... begin qz step ..........
  112.    60 if (en .le. 2) go to 1001
  113.       if (.not. matz) enorn = en
  114.       its = 0
  115.       na = en - 1
  116.       enm2 = na - 1
  117.    70 ish = 2
  118. c     .......... check for convergence or reducibility.
  119. c                for l=en step -1 until 1 do -- ..........
  120.       do 80 ll = 1, en
  121.          lm1 = en - ll
  122.          l = lm1 + 1
  123.          if (l .eq. 1) go to 95
  124.          if (dabs(a(l,lm1)) .le. epsa) go to 90
  125.    80 continue
  126. c
  127.    90 a(l,lm1) = 0.0d0
  128.       if (l .lt. na) go to 95
  129. c     .......... 1-by-1 or 2-by-2 block isolated ..........
  130.       en = lm1
  131.       go to 60
  132. c     .......... check for small top of b ..........
  133.    95 ld = l
  134.   100 l1 = l + 1
  135.       b11 = b(l,l)
  136.       if (dabs(b11) .gt. epsb) go to 120
  137.       b(l,l) = 0.0d0
  138.       s = dabs(a(l,l)) + dabs(a(l1,l))
  139.       u1 = a(l,l) / s
  140.       u2 = a(l1,l) / s
  141.       r = dsign(dsqrt(u1*u1+u2*u2),u1)
  142.       v1 = -(u1 + r) / r
  143.       v2 = -u2 / r
  144.       u2 = v2 / v1
  145. c
  146.       do 110 j = l, enorn
  147.          t = a(l,j) + u2 * a(l1,j)
  148.          a(l,j) = a(l,j) + t * v1
  149.          a(l1,j) = a(l1,j) + t * v2
  150.          t = b(l,j) + u2 * b(l1,j)
  151.          b(l,j) = b(l,j) + t * v1
  152.          b(l1,j) = b(l1,j) + t * v2
  153.   110 continue
  154. c
  155.       if (l .ne. 1) a(l,lm1) = -a(l,lm1)
  156.       lm1 = l
  157.       l = l1
  158.       go to 90
  159.   120 a11 = a(l,l) / b11
  160.       a21 = a(l1,l) / b11
  161.       if (ish .eq. 1) go to 140
  162. c     .......... iteration strategy ..........
  163.       if (itn .eq. 0) go to 1000
  164.       if (its .eq. 10) go to 155
  165. c     .......... determine type of shift ..........
  166.       b22 = b(l1,l1)
  167.       if (dabs(b22) .lt. epsb) b22 = epsb
  168.       b33 = b(na,na)
  169.       if (dabs(b33) .lt. epsb) b33 = epsb
  170.       b44 = b(en,en)
  171.       if (dabs(b44) .lt. epsb) b44 = epsb
  172.       a33 = a(na,na) / b33
  173.       a34 = a(na,en) / b44
  174.       a43 = a(en,na) / b33
  175.       a44 = a(en,en) / b44
  176.       b34 = b(na,en) / b44
  177.       t = 0.5d0 * (a43 * b34 - a33 - a44)
  178.       r = t * t + a34 * a43 - a33 * a44
  179.       if (r .lt. 0.0d0) go to 150
  180. c     .......... determine single shift zeroth column of a ..........
  181.       ish = 1
  182.       r = dsqrt(r)
  183.       sh = -t + r
  184.       s = -t - r
  185.       if (dabs(s-a44) .lt. dabs(sh-a44)) sh = s
  186. c     .......... look for two consecutive small
  187. c                sub-diagonal elements of a.
  188. c                for l=en-2 step -1 until ld do -- ..........
  189.       do 130 ll = ld, enm2
  190.          l = enm2 + ld - ll
  191.          if (l .eq. ld) go to 140
  192.          lm1 = l - 1
  193.          l1 = l + 1
  194.          t = a(l,l)
  195.          if (dabs(b(l,l)) .gt. epsb) t = t - sh * b(l,l)
  196.          if (dabs(a(l,lm1)) .le. dabs(t/a(l1,l)) * epsa) go to 100
  197.   130 continue
  198. c
  199.   140 a1 = a11 - sh
  200.       a2 = a21
  201.       if (l .ne. ld) a(l,lm1) = -a(l,lm1)
  202.       go to 160
  203. c     .......... determine double shift zeroth column of a ..........
  204.   150 a12 = a(l,l1) / b22
  205.       a22 = a(l1,l1) / b22
  206.       b12 = b(l,l1) / b22
  207.       a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11)
  208.      x     / a21 + a12 - a11 * b12
  209.       a2 = (a22 - a11) - a21 * b12 - (a33 - a11) - (a44 - a11)
  210.      x     + a43 * b34
  211.       a3 = a(l1+1,l1) / b22
  212.       go to 160
  213. c     .......... ad hoc shift ..........
  214.   155 a1 = 0.0d0
  215.       a2 = 1.0d0
  216.       a3 = 1.1605d0
  217.   160 its = its + 1
  218.       itn = itn - 1
  219.       if (.not. matz) lor1 = ld
  220. c     .......... main loop ..........
  221.       do 260 k = l, na
  222.          notlas = k .ne. na .and. ish .eq. 2
  223.          k1 = k + 1
  224.          k2 = k + 2
  225.          km1 = max0(k-1,l)
  226.          ll = min0(en,k1+ish)
  227.          if (notlas) go to 190
  228. c     .......... zero a(k+1,k-1) ..........
  229.          if (k .eq. l) go to 170
  230.          a1 = a(k,km1)
  231.          a2 = a(k1,km1)
  232.   170    s = dabs(a1) + dabs(a2)
  233.          if (s .eq. 0.0d0) go to 70
  234.          u1 = a1 / s
  235.          u2 = a2 / s
  236.          r = dsign(dsqrt(u1*u1+u2*u2),u1)
  237.          v1 = -(u1 + r) / r
  238.          v2 = -u2 / r
  239.          u2 = v2 / v1
  240. c
  241.          do 180 j = km1, enorn
  242.             t = a(k,j) + u2 * a(k1,j)
  243.             a(k,j) = a(k,j) + t * v1
  244.             a(k1,j) = a(k1,j) + t * v2
  245.             t = b(k,j) + u2 * b(k1,j)
  246.             b(k,j) = b(k,j) + t * v1
  247.             b(k1,j) = b(k1,j) + t * v2
  248.   180    continue
  249. c
  250.          if (k .ne. l) a(k1,km1) = 0.0d0
  251.          go to 240
  252. c     .......... zero a(k+1,k-1) and a(k+2,k-1) ..........
  253.   190    if (k .eq. l) go to 200
  254.          a1 = a(k,km1)
  255.          a2 = a(k1,km1)
  256.          a3 = a(k2,km1)
  257.   200    s = dabs(a1) + dabs(a2) + dabs(a3)
  258.          if (s .eq. 0.0d0) go to 260
  259.          u1 = a1 / s
  260.          u2 = a2 / s
  261.          u3 = a3 / s
  262.          r = dsign(dsqrt(u1*u1+u2*u2+u3*u3),u1)
  263.          v1 = -(u1 + r) / r
  264.          v2 = -u2 / r
  265.          v3 = -u3 / r
  266.          u2 = v2 / v1
  267.          u3 = v3 / v1
  268. c
  269.          do 210 j = km1, enorn
  270.             t = a(k,j) + u2 * a(k1,j) + u3 * a(k2,j)
  271.             a(k,j) = a(k,j) + t * v1
  272.             a(k1,j) = a(k1,j) + t * v2
  273.             a(k2,j) = a(k2,j) + t * v3
  274.             t = b(k,j) + u2 * b(k1,j) + u3 * b(k2,j)
  275.             b(k,j) = b(k,j) + t * v1
  276.             b(k1,j) = b(k1,j) + t * v2
  277.             b(k2,j) = b(k2,j) + t * v3
  278.   210    continue
  279. c
  280.          if (k .eq. l) go to 220
  281.          a(k1,km1) = 0.0d0
  282.          a(k2,km1) = 0.0d0
  283. c     .......... zero b(k+2,k+1) and b(k+2,k) ..........
  284.   220    s = dabs(b(k2,k2)) + dabs(b(k2,k1)) + dabs(b(k2,k))
  285.          if (s .eq. 0.0d0) go to 240
  286.          u1 = b(k2,k2) / s
  287.          u2 = b(k2,k1) / s
  288.          u3 = b(k2,k) / s
  289.          r = dsign(dsqrt(u1*u1+u2*u2+u3*u3),u1)
  290.          v1 = -(u1 + r) / r
  291.          v2 = -u2 / r
  292.          v3 = -u3 / r
  293.          u2 = v2 / v1
  294.          u3 = v3 / v1
  295. c
  296.          do 230 i = lor1, ll
  297.             t = a(i,k2) + u2 * a(i,k1) + u3 * a(i,k)
  298.             a(i,k2) = a(i,k2) + t * v1
  299.             a(i,k1) = a(i,k1) + t * v2
  300.             a(i,k) = a(i,k) + t * v3
  301.             t = b(i,k2) + u2 * b(i,k1) + u3 * b(i,k)
  302.             b(i,k2) = b(i,k2) + t * v1
  303.             b(i,k1) = b(i,k1) + t * v2
  304.             b(i,k) = b(i,k) + t * v3
  305.   230    continue
  306. c
  307.          b(k2,k) = 0.0d0
  308.          b(k2,k1) = 0.0d0
  309.          if (.not. matz) go to 240
  310. c
  311.          do 235 i = 1, n
  312.             t = z(i,k2) + u2 * z(i,k1) + u3 * z(i,k)
  313.             z(i,k2) = z(i,k2) + t * v1
  314.             z(i,k1) = z(i,k1) + t * v2
  315.             z(i,k) = z(i,k) + t * v3
  316.   235    continue
  317. c     .......... zero b(k+1,k) ..........
  318.   240    s = dabs(b(k1,k1)) + dabs(b(k1,k))
  319.          if (s .eq. 0.0d0) go to 260
  320.          u1 = b(k1,k1) / s
  321.          u2 = b(k1,k) / s
  322.          r = dsign(dsqrt(u1*u1+u2*u2),u1)
  323.          v1 = -(u1 + r) / r
  324.          v2 = -u2 / r
  325.          u2 = v2 / v1
  326. c
  327.          do 250 i = lor1, ll
  328.             t = a(i,k1) + u2 * a(i,k)
  329.             a(i,k1) = a(i,k1) + t * v1
  330.             a(i,k) = a(i,k) + t * v2
  331.             t = b(i,k1) + u2 * b(i,k)
  332.             b(i,k1) = b(i,k1) + t * v1
  333.             b(i,k) = b(i,k) + t * v2
  334.   250    continue
  335. c
  336.          b(k1,k) = 0.0d0
  337.          if (.not. matz) go to 260
  338. c
  339.          do 255 i = 1, n
  340.             t = z(i,k1) + u2 * z(i,k)
  341.             z(i,k1) = z(i,k1) + t * v1
  342.             z(i,k) = z(i,k) + t * v2
  343.   255    continue
  344. c
  345.   260 continue
  346. c     .......... end qz step ..........
  347.       go to 70
  348. c     .......... set error -- all eigenvalues have not
  349. c                converged after 30*n iterations ..........
  350.  1000 ierr = en
  351. c     .......... save epsb for use by qzval and qzvec ..........
  352.  1001 if (n .gt. 1) b(n,1) = epsb
  353.       return
  354.       end
  355.