home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / qzvec.f < prev    next >
Text File  |  1996-09-28  |  9KB  |  249 lines
  1.       subroutine qzvec(nm,n,a,b,alfr,alfi,beta,z)
  2. c
  3.       integer i,j,k,m,n,en,ii,jj,na,nm,nn,isw,enm2
  4.       double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n)
  5.       double precision d,q,r,s,t,w,x,y,di,dr,ra,rr,sa,ti,tr,t1,t2,w1,x1,
  6.      x       zz,z1,alfm,almi,almr,betm,epsb
  7. c
  8. c     this subroutine is the optional fourth 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 matrices, one of them in
  13. c     quasi-triangular form (in which each 2-by-2 block corresponds to
  14. c     a pair of complex eigenvalues) and the other in upper triangular
  15. c     form.  it computes the eigenvectors of the triangular problem and
  16. c     transforms the results back to the original coordinate system.
  17. c     it is usually preceded by  qzhes,  qzit, and  qzval.
  18. c
  19. c     on input
  20. c
  21. c        nm must be set to the row dimension of two-dimensional
  22. c          array parameters as declared in the calling program
  23. c          dimension statement.
  24. c
  25. c        n is the order of the matrices.
  26. c
  27. c        a contains a real upper quasi-triangular matrix.
  28. c
  29. c        b contains a real upper triangular matrix.  in addition,
  30. c          location b(n,1) contains the tolerance quantity (epsb)
  31. c          computed and saved in  qzit.
  32. c
  33. c        alfr, alfi, and beta  are vectors with components whose
  34. c          ratios ((alfr+i*alfi)/beta) are the generalized
  35. c          eigenvalues.  they are usually obtained from  qzval.
  36. c
  37. c        z contains the transformation matrix produced in the
  38. c          reductions by  qzhes,  qzit, and  qzval, if performed.
  39. c          if the eigenvectors of the triangular problem are
  40. c          desired, z must contain the identity matrix.
  41. c
  42. c     on output
  43. c
  44. c        a is unaltered.  its subdiagonal elements provide information
  45. c           about the storage of the complex eigenvectors.
  46. c
  47. c        b has been destroyed.
  48. c
  49. c        alfr, alfi, and beta are unaltered.
  50. c
  51. c        z contains the real and imaginary parts of the eigenvectors.
  52. c          if alfi(i) .eq. 0.0, the i-th eigenvalue is real and
  53. c            the i-th column of z contains its eigenvector.
  54. c          if alfi(i) .ne. 0.0, the i-th eigenvalue is complex.
  55. c            if alfi(i) .gt. 0.0, the eigenvalue is the first of
  56. c              a complex pair and the i-th and (i+1)-th columns
  57. c              of z contain its eigenvector.
  58. c            if alfi(i) .lt. 0.0, the eigenvalue is the second of
  59. c              a complex pair and the (i-1)-th and i-th columns
  60. c              of z contain the conjugate of its eigenvector.
  61. c          each eigenvector is normalized so that the modulus
  62. c          of its largest component is 1.0 .
  63. c
  64. c     questions and comments should be directed to burton s. garbow,
  65. c     mathematics and computer science div, argonne national laboratory
  66. c
  67. c     this version dated august 1983.
  68. c
  69. c     ------------------------------------------------------------------
  70. c
  71.       epsb = b(n,1)
  72.       isw = 1
  73. c     .......... for en=n step -1 until 1 do -- ..........
  74.       do 800 nn = 1, n
  75.          en = n + 1 - nn
  76.          na = en - 1
  77.          if (isw .eq. 2) go to 795
  78.          if (alfi(en) .ne. 0.0d0) go to 710
  79. c     .......... real vector ..........
  80.          m = en
  81.          b(en,en) = 1.0d0
  82.          if (na .eq. 0) go to 800
  83.          alfm = alfr(m)
  84.          betm = beta(m)
  85. c     .......... for i=en-1 step -1 until 1 do -- ..........
  86.          do 700 ii = 1, na
  87.             i = en - ii
  88.             w = betm * a(i,i) - alfm * b(i,i)
  89.             r = 0.0d0
  90. c
  91.             do 610 j = m, en
  92.   610       r = r + (betm * a(i,j) - alfm * b(i,j)) * b(j,en)
  93. c
  94.             if (i .eq. 1 .or. isw .eq. 2) go to 630
  95.             if (betm * a(i,i-1) .eq. 0.0d0) go to 630
  96.             zz = w
  97.             s = r
  98.             go to 690
  99.   630       m = i
  100.             if (isw .eq. 2) go to 640
  101. c     .......... real 1-by-1 block ..........
  102.             t = w
  103.             if (w .eq. 0.0d0) t = epsb
  104.             b(i,en) = -r / t
  105.             go to 700
  106. c     .......... real 2-by-2 block ..........
  107.   640       x = betm * a(i,i+1) - alfm * b(i,i+1)
  108.             y = betm * a(i+1,i)
  109.             q = w * zz - x * y
  110.             t = (x * s - zz * r) / q
  111.             b(i,en) = t
  112.             if (dabs(x) .le. dabs(zz)) go to 650
  113.             b(i+1,en) = (-r - w * t) / x
  114.             go to 690
  115.   650       b(i+1,en) = (-s - y * t) / zz
  116.   690       isw = 3 - isw
  117.   700    continue
  118. c     .......... end real vector ..........
  119.          go to 800
  120. c     .......... complex vector ..........
  121.   710    m = na
  122.          almr = alfr(m)
  123.          almi = alfi(m)
  124.          betm = beta(m)
  125. c     .......... last vector component chosen imaginary so that
  126. c                eigenvector matrix is triangular ..........
  127.          y = betm * a(en,na)
  128.          b(na,na) = -almi * b(en,en) / y
  129.          b(na,en) = (almr * b(en,en) - betm * a(en,en)) / y
  130.          b(en,na) = 0.0d0
  131.          b(en,en) = 1.0d0
  132.          enm2 = na - 1
  133.          if (enm2 .eq. 0) go to 795
  134. c     .......... for i=en-2 step -1 until 1 do -- ..........
  135.          do 790 ii = 1, enm2
  136.             i = na - ii
  137.             w = betm * a(i,i) - almr * b(i,i)
  138.             w1 = -almi * b(i,i)
  139.             ra = 0.0d0
  140.             sa = 0.0d0
  141. c
  142.             do 760 j = m, en
  143.                x = betm * a(i,j) - almr * b(i,j)
  144.                x1 = -almi * b(i,j)
  145.                ra = ra + x * b(j,na) - x1 * b(j,en)
  146.                sa = sa + x * b(j,en) + x1 * b(j,na)
  147.   760       continue
  148. c
  149.             if (i .eq. 1 .or. isw .eq. 2) go to 770
  150.             if (betm * a(i,i-1) .eq. 0.0d0) go to 770
  151.             zz = w
  152.             z1 = w1
  153.             r = ra
  154.             s = sa
  155.             isw = 2
  156.             go to 790
  157.   770       m = i
  158.             if (isw .eq. 2) go to 780
  159. c     .......... complex 1-by-1 block ..........
  160.             tr = -ra
  161.             ti = -sa
  162.   773       dr = w
  163.             di = w1
  164. c     .......... complex divide (t1,t2) = (tr,ti) / (dr,di) ..........
  165.   775       if (dabs(di) .gt. dabs(dr)) go to 777
  166.             rr = di / dr
  167.             d = dr + di * rr
  168.             t1 = (tr + ti * rr) / d
  169.             t2 = (ti - tr * rr) / d
  170.             go to (787,782), isw
  171.   777       rr = dr / di
  172.             d = dr * rr + di
  173.             t1 = (tr * rr + ti) / d
  174.             t2 = (ti * rr - tr) / d
  175.             go to (787,782), isw
  176. c     .......... complex 2-by-2 block ..........
  177.   780       x = betm * a(i,i+1) - almr * b(i,i+1)
  178.             x1 = -almi * b(i,i+1)
  179.             y = betm * a(i+1,i)
  180.             tr = y * ra - w * r + w1 * s
  181.             ti = y * sa - w * s - w1 * r
  182.             dr = w * zz - w1 * z1 - x * y
  183.             di = w * z1 + w1 * zz - x1 * y
  184.             if (dr .eq. 0.0d0 .and. di .eq. 0.0d0) dr = epsb
  185.             go to 775
  186.   782       b(i+1,na) = t1
  187.             b(i+1,en) = t2
  188.             isw = 1
  189.             if (dabs(y) .gt. dabs(w) + dabs(w1)) go to 785
  190.             tr = -ra - x * b(i+1,na) + x1 * b(i+1,en)
  191.             ti = -sa - x * b(i+1,en) - x1 * b(i+1,na)
  192.             go to 773
  193.   785       t1 = (-r - zz * b(i+1,na) + z1 * b(i+1,en)) / y
  194.             t2 = (-s - zz * b(i+1,en) - z1 * b(i+1,na)) / y
  195.   787       b(i,na) = t1
  196.             b(i,en) = t2
  197.   790    continue
  198. c     .......... end complex vector ..........
  199.   795    isw = 3 - isw
  200.   800 continue
  201. c     .......... end back substitution.
  202. c                transform to original coordinate system.
  203. c                for j=n step -1 until 1 do -- ..........
  204.       do 880 jj = 1, n
  205.          j = n + 1 - jj
  206. c
  207.          do 880 i = 1, n
  208.             zz = 0.0d0
  209. c
  210.             do 860 k = 1, j
  211.   860       zz = zz + z(i,k) * b(k,j)
  212. c
  213.             z(i,j) = zz
  214.   880 continue
  215. c     .......... normalize so that modulus of largest
  216. c                component of each vector is 1.
  217. c                (isw is 1 initially from before) ..........
  218.       do 950 j = 1, n
  219.          d = 0.0d0
  220.          if (isw .eq. 2) go to 920
  221.          if (alfi(j) .ne. 0.0d0) go to 945
  222. c
  223.          do 890 i = 1, n
  224.             if (dabs(z(i,j)) .gt. d) d = dabs(z(i,j))
  225.   890    continue
  226. c
  227.          do 900 i = 1, n
  228.   900    z(i,j) = z(i,j) / d
  229. c
  230.          go to 950
  231. c
  232.   920    do 930 i = 1, n
  233.             r = dabs(z(i,j-1)) + dabs(z(i,j))
  234.             if (r .ne. 0.0d0) r = r * dsqrt((z(i,j-1)/r)**2
  235.      x                                     +(z(i,j)/r)**2)
  236.             if (r .gt. d) d = r
  237.   930    continue
  238. c
  239.          do 940 i = 1, n
  240.             z(i,j-1) = z(i,j-1) / d
  241.             z(i,j) = z(i,j) / d
  242.   940    continue
  243. c
  244.   945    isw = 3 - isw
  245.   950 continue
  246. c
  247.       return
  248.       end
  249.