home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / eispack / qzhes.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  196 lines
  1.       subroutine qzhes(nm,n,a,b,matz,z)
  2. c
  3.       integer i,j,k,l,n,lb,l1,nm,nk1,nm1,nm2
  4.       double precision a(nm,n),b(nm,n),z(nm,n)
  5.       double precision r,s,t,u1,u2,v1,v2,rho
  6.       logical matz
  7. c
  8. c     this subroutine is the first step of the qz algorithm
  9. c     for solving generalized matrix eigenvalue problems,
  10. c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
  11. c
  12. c     this subroutine accepts a pair of real general matrices and
  13. c     reduces one of them to upper hessenberg form and the other
  14. c     to upper triangular form using orthogonal transformations.
  15. c     it is usually followed by  qzit,  qzval  and, possibly,  qzvec.
  16. c
  17. c     on input
  18. c
  19. c        nm must be set to the row dimension of two-dimensional
  20. c          array parameters as declared in the calling program
  21. c          dimension statement.
  22. c
  23. c        n is the order of the matrices.
  24. c
  25. c        a contains a real general matrix.
  26. c
  27. c        b contains a real general matrix.
  28. c
  29. c        matz should be set to .true. if the right hand transformations
  30. c          are to be accumulated for later use in computing
  31. c          eigenvectors, and to .false. otherwise.
  32. c
  33. c     on output
  34. c
  35. c        a has been reduced to upper hessenberg form.  the elements
  36. c          below the first subdiagonal have been set to zero.
  37. c
  38. c        b has been reduced to upper triangular form.  the elements
  39. c          below the main diagonal have been set to zero.
  40. c
  41. c        z contains the product of the right hand transformations if
  42. c          matz has been set to .true.  otherwise, z is not referenced.
  43. c
  44. c     questions and comments should be directed to burton s. garbow,
  45. c     mathematics and computer science div, argonne national laboratory
  46. c
  47. c     this version dated august 1983.
  48. c
  49. c     ------------------------------------------------------------------
  50. c
  51. c     .......... initialize z ..........
  52.       if (.not. matz) go to 10
  53. c
  54.       do 3 j = 1, n
  55. c
  56.          do 2 i = 1, n
  57.             z(i,j) = 0.0d0
  58.     2    continue
  59. c
  60.          z(j,j) = 1.0d0
  61.     3 continue
  62. c     .......... reduce b to upper triangular form ..........
  63.    10 if (n .le. 1) go to 170
  64.       nm1 = n - 1
  65. c
  66.       do 100 l = 1, nm1
  67.          l1 = l + 1
  68.          s = 0.0d0
  69. c
  70.          do 20 i = l1, n
  71.             s = s + dabs(b(i,l))
  72.    20    continue
  73. c
  74.          if (s .eq. 0.0d0) go to 100
  75.          s = s + dabs(b(l,l))
  76.          r = 0.0d0
  77. c
  78.          do 25 i = l, n
  79.             b(i,l) = b(i,l) / s
  80.             r = r + b(i,l)**2
  81.    25    continue
  82. c
  83.          r = dsign(dsqrt(r),b(l,l))
  84.          b(l,l) = b(l,l) + r
  85.          rho = r * b(l,l)
  86. c
  87.          do 50 j = l1, n
  88.             t = 0.0d0
  89. c
  90.             do 30 i = l, n
  91.                t = t + b(i,l) * b(i,j)
  92.    30       continue
  93. c
  94.             t = -t / rho
  95. c
  96.             do 40 i = l, n
  97.                b(i,j) = b(i,j) + t * b(i,l)
  98.    40       continue
  99. c
  100.    50    continue
  101. c
  102.          do 80 j = 1, n
  103.             t = 0.0d0
  104. c
  105.             do 60 i = l, n
  106.                t = t + b(i,l) * a(i,j)
  107.    60       continue
  108. c
  109.             t = -t / rho
  110. c
  111.             do 70 i = l, n
  112.                a(i,j) = a(i,j) + t * b(i,l)
  113.    70       continue
  114. c
  115.    80    continue
  116. c
  117.          b(l,l) = -s * r
  118. c
  119.          do 90 i = l1, n
  120.             b(i,l) = 0.0d0
  121.    90    continue
  122. c
  123.   100 continue
  124. c     .......... reduce a to upper hessenberg form, while
  125. c                keeping b triangular ..........
  126.       if (n .eq. 2) go to 170
  127.       nm2 = n - 2
  128. c
  129.       do 160 k = 1, nm2
  130.          nk1 = nm1 - k
  131. c     .......... for l=n-1 step -1 until k+1 do -- ..........
  132.          do 150 lb = 1, nk1
  133.             l = n - lb
  134.             l1 = l + 1
  135. c     .......... zero a(l+1,k) ..........
  136.             s = dabs(a(l,k)) + dabs(a(l1,k))
  137.             if (s .eq. 0.0d0) go to 150
  138.             u1 = a(l,k) / s
  139.             u2 = a(l1,k) / s
  140.             r = dsign(dsqrt(u1*u1+u2*u2),u1)
  141.             v1 =  -(u1 + r) / r
  142.             v2 = -u2 / r
  143.             u2 = v2 / v1
  144. c
  145.             do 110 j = k, n
  146.                t = a(l,j) + u2 * a(l1,j)
  147.                a(l,j) = a(l,j) + t * v1
  148.                a(l1,j) = a(l1,j) + t * v2
  149.   110       continue
  150. c
  151.             a(l1,k) = 0.0d0
  152. c
  153.             do 120 j = l, n
  154.                t = b(l,j) + u2 * b(l1,j)
  155.                b(l,j) = b(l,j) + t * v1
  156.                b(l1,j) = b(l1,j) + t * v2
  157.   120       continue
  158. c     .......... zero b(l+1,l) ..........
  159.             s = dabs(b(l1,l1)) + dabs(b(l1,l))
  160.             if (s .eq. 0.0d0) go to 150
  161.             u1 = b(l1,l1) / s
  162.             u2 = b(l1,l) / s
  163.             r = dsign(dsqrt(u1*u1+u2*u2),u1)
  164.             v1 =  -(u1 + r) / r
  165.             v2 = -u2 / r
  166.             u2 = v2 / v1
  167. c
  168.             do 130 i = 1, l1
  169.                t = b(i,l1) + u2 * b(i,l)
  170.                b(i,l1) = b(i,l1) + t * v1
  171.                b(i,l) = b(i,l) + t * v2
  172.   130       continue
  173. c
  174.             b(l1,l) = 0.0d0
  175. c
  176.             do 140 i = 1, n
  177.                t = a(i,l1) + u2 * a(i,l)
  178.                a(i,l1) = a(i,l1) + t * v1
  179.                a(i,l) = a(i,l) + t * v2
  180.   140       continue
  181. c
  182.             if (.not. matz) go to 150
  183. c
  184.             do 145 i = 1, n
  185.                t = z(i,l1) + u2 * z(i,l)
  186.                z(i,l1) = z(i,l1) + t * v1
  187.                z(i,l) = z(i,l) + t * v2
  188.   145       continue
  189. c
  190.   150    continue
  191. c
  192.   160 continue
  193. c
  194.   170 return
  195.       end
  196.