home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / qzval.f < prev    next >
Text File  |  1996-09-28  |  10KB  |  278 lines
  1.       subroutine qzval(nm,n,a,b,alfr,alfi,beta,matz,z)
  2. c
  3.       integer i,j,n,en,na,nm,nn,isw
  4.       double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n)
  5.       double precision c,d,e,r,s,t,an,a1,a2,bn,cq,cz,di,dr,ei,ti,tr,u1,
  6.      x       u2,v1,v2,a1i,a11,a12,a2i,a21,a22,b11,b12,b22,sqi,sqr,
  7.      x       ssi,ssr,szi,szr,a11i,a11r,a12i,a12r,a22i,a22r,epsb
  8.       logical matz
  9. c
  10. c     this subroutine is the third step of the qz algorithm
  11. c     for solving generalized matrix eigenvalue problems,
  12. c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
  13. c
  14. c     this subroutine accepts a pair of real matrices, one of them
  15. c     in quasi-triangular form and the other in upper triangular form.
  16. c     it reduces the quasi-triangular matrix further, so that any
  17. c     remaining 2-by-2 blocks correspond to pairs of complex
  18. c     eigenvalues, and returns quantities whose ratios give the
  19. c     generalized eigenvalues.  it is usually preceded by  qzhes
  20. c     and  qzit  and may be followed by  qzvec.
  21. c
  22. c     on input
  23. c
  24. c        nm must be set to the row dimension of two-dimensional
  25. c          array parameters as declared in the calling program
  26. c          dimension statement.
  27. c
  28. c        n is the order of the matrices.
  29. c
  30. c        a contains a real upper quasi-triangular matrix.
  31. c
  32. c        b contains a real upper triangular matrix.  in addition,
  33. c          location b(n,1) contains the tolerance quantity (epsb)
  34. c          computed and saved in  qzit.
  35. c
  36. c        matz should be set to .true. if the right hand transformations
  37. c          are to be accumulated for later use in computing
  38. c          eigenvectors, and to .false. otherwise.
  39. c
  40. c        z contains, if matz has been set to .true., the
  41. c          transformation matrix produced in the reductions by qzhes
  42. c          and qzit, if performed, or else the identity matrix.
  43. c          if matz has been set to .false., z is not referenced.
  44. c
  45. c     on output
  46. c
  47. c        a has been reduced further to a quasi-triangular matrix
  48. c          in which all nonzero subdiagonal elements correspond to
  49. c          pairs of complex eigenvalues.
  50. c
  51. c        b is still in upper triangular form, although its elements
  52. c          have been altered.  b(n,1) is unaltered.
  53. c
  54. c        alfr and alfi contain the real and imaginary parts of the
  55. c          diagonal elements of the triangular matrix that would be
  56. c          obtained if a were reduced completely to triangular form
  57. c          by unitary transformations.  non-zero values of alfi occur
  58. c          in pairs, the first member positive and the second negative.
  59. c
  60. c        beta contains the diagonal elements of the corresponding b,
  61. c          normalized to be real and non-negative.  the generalized
  62. c          eigenvalues are then the ratios ((alfr+i*alfi)/beta).
  63. c
  64. c        z contains the product of the right hand transformations
  65. c          (for all three steps) if matz has been set to .true.
  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.       epsb = b(n,1)
  75.       isw = 1
  76. c     .......... find eigenvalues of quasi-triangular matrices.
  77. c                for en=n step -1 until 1 do -- ..........
  78.       do 510 nn = 1, n
  79.          en = n + 1 - nn
  80.          na = en - 1
  81.          if (isw .eq. 2) go to 505
  82.          if (en .eq. 1) go to 410
  83.          if (a(en,na) .ne. 0.0d0) go to 420
  84. c     .......... 1-by-1 block, one real root ..........
  85.   410    alfr(en) = a(en,en)
  86.          if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en)
  87.          beta(en) = dabs(b(en,en))
  88.          alfi(en) = 0.0d0
  89.          go to 510
  90. c     .......... 2-by-2 block ..........
  91.   420    if (dabs(b(na,na)) .le. epsb) go to 455
  92.          if (dabs(b(en,en)) .gt. epsb) go to 430
  93.          a1 = a(en,en)
  94.          a2 = a(en,na)
  95.          bn = 0.0d0
  96.          go to 435
  97.   430    an = dabs(a(na,na)) + dabs(a(na,en)) + dabs(a(en,na))
  98.      x      + dabs(a(en,en))
  99.          bn = dabs(b(na,na)) + dabs(b(na,en)) + dabs(b(en,en))
  100.          a11 = a(na,na) / an
  101.          a12 = a(na,en) / an
  102.          a21 = a(en,na) / an
  103.          a22 = a(en,en) / an
  104.          b11 = b(na,na) / bn
  105.          b12 = b(na,en) / bn
  106.          b22 = b(en,en) / bn
  107.          e = a11 / b11
  108.          ei = a22 / b22
  109.          s = a21 / (b11 * b22)
  110.          t = (a22 - e * b22) / b22
  111.          if (dabs(e) .le. dabs(ei)) go to 431
  112.          e = ei
  113.          t = (a11 - e * b11) / b11
  114.   431    c = 0.5d0 * (t - s * b12)
  115.          d = c * c + s * (a12 - e * b12)
  116.          if (d .lt. 0.0d0) go to 480
  117. c     .......... two real roots.
  118. c                zero both a(en,na) and b(en,na) ..........
  119.          e = e + (c + dsign(dsqrt(d),c))
  120.          a11 = a11 - e * b11
  121.          a12 = a12 - e * b12
  122.          a22 = a22 - e * b22
  123.          if (dabs(a11) + dabs(a12) .lt.
  124.      x       dabs(a21) + dabs(a22)) go to 432
  125.          a1 = a12
  126.          a2 = a11
  127.          go to 435
  128.   432    a1 = a22
  129.          a2 = a21
  130. c     .......... choose and apply real z ..........
  131.   435    s = dabs(a1) + dabs(a2)
  132.          u1 = a1 / s
  133.          u2 = a2 / s
  134.          r = dsign(dsqrt(u1*u1+u2*u2),u1)
  135.          v1 = -(u1 + r) / r
  136.          v2 = -u2 / r
  137.          u2 = v2 / v1
  138. c
  139.          do 440 i = 1, en
  140.             t = a(i,en) + u2 * a(i,na)
  141.             a(i,en) = a(i,en) + t * v1
  142.             a(i,na) = a(i,na) + t * v2
  143.             t = b(i,en) + u2 * b(i,na)
  144.             b(i,en) = b(i,en) + t * v1
  145.             b(i,na) = b(i,na) + t * v2
  146.   440    continue
  147. c
  148.          if (.not. matz) go to 450
  149. c
  150.          do 445 i = 1, n
  151.             t = z(i,en) + u2 * z(i,na)
  152.             z(i,en) = z(i,en) + t * v1
  153.             z(i,na) = z(i,na) + t * v2
  154.   445    continue
  155. c
  156.   450    if (bn .eq. 0.0d0) go to 475
  157.          if (an .lt. dabs(e) * bn) go to 455
  158.          a1 = b(na,na)
  159.          a2 = b(en,na)
  160.          go to 460
  161.   455    a1 = a(na,na)
  162.          a2 = a(en,na)
  163. c     .......... choose and apply real q ..........
  164.   460    s = dabs(a1) + dabs(a2)
  165.          if (s .eq. 0.0d0) go to 475
  166.          u1 = a1 / s
  167.          u2 = a2 / s
  168.          r = dsign(dsqrt(u1*u1+u2*u2),u1)
  169.          v1 = -(u1 + r) / r
  170.          v2 = -u2 / r
  171.          u2 = v2 / v1
  172. c
  173.          do 470 j = na, n
  174.             t = a(na,j) + u2 * a(en,j)
  175.             a(na,j) = a(na,j) + t * v1
  176.             a(en,j) = a(en,j) + t * v2
  177.             t = b(na,j) + u2 * b(en,j)
  178.             b(na,j) = b(na,j) + t * v1
  179.             b(en,j) = b(en,j) + t * v2
  180.   470    continue
  181. c
  182.   475    a(en,na) = 0.0d0
  183.          b(en,na) = 0.0d0
  184.          alfr(na) = a(na,na)
  185.          alfr(en) = a(en,en)
  186.          if (b(na,na) .lt. 0.0d0) alfr(na) = -alfr(na)
  187.          if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en)
  188.          beta(na) = dabs(b(na,na))
  189.          beta(en) = dabs(b(en,en))
  190.          alfi(en) = 0.0d0
  191.          alfi(na) = 0.0d0
  192.          go to 505
  193. c     .......... two complex roots ..........
  194.   480    e = e + c
  195.          ei = dsqrt(-d)
  196.          a11r = a11 - e * b11
  197.          a11i = ei * b11
  198.          a12r = a12 - e * b12
  199.          a12i = ei * b12
  200.          a22r = a22 - e * b22
  201.          a22i = ei * b22
  202.          if (dabs(a11r) + dabs(a11i) + dabs(a12r) + dabs(a12i) .lt.
  203.      x       dabs(a21) + dabs(a22r) + dabs(a22i)) go to 482
  204.          a1 = a12r
  205.          a1i = a12i
  206.          a2 = -a11r
  207.          a2i = -a11i
  208.          go to 485
  209.   482    a1 = a22r
  210.          a1i = a22i
  211.          a2 = -a21
  212.          a2i = 0.0d0
  213. c     .......... choose complex z ..........
  214.   485    cz = dsqrt(a1*a1+a1i*a1i)
  215.          if (cz .eq. 0.0d0) go to 487
  216.          szr = (a1 * a2 + a1i * a2i) / cz
  217.          szi = (a1 * a2i - a1i * a2) / cz
  218.          r = dsqrt(cz*cz+szr*szr+szi*szi)
  219.          cz = cz / r
  220.          szr = szr / r
  221.          szi = szi / r
  222.          go to 490
  223.   487    szr = 1.0d0
  224.          szi = 0.0d0
  225.   490    if (an .lt. (dabs(e) + ei) * bn) go to 492
  226.          a1 = cz * b11 + szr * b12
  227.          a1i = szi * b12
  228.          a2 = szr * b22
  229.          a2i = szi * b22
  230.          go to 495
  231.   492    a1 = cz * a11 + szr * a12
  232.          a1i = szi * a12
  233.          a2 = cz * a21 + szr * a22
  234.          a2i = szi * a22
  235. c     .......... choose complex q ..........
  236.   495    cq = dsqrt(a1*a1+a1i*a1i)
  237.          if (cq .eq. 0.0d0) go to 497
  238.          sqr = (a1 * a2 + a1i * a2i) / cq
  239.          sqi = (a1 * a2i - a1i * a2) / cq
  240.          r = dsqrt(cq*cq+sqr*sqr+sqi*sqi)
  241.          cq = cq / r
  242.          sqr = sqr / r
  243.          sqi = sqi / r
  244.          go to 500
  245.   497    sqr = 1.0d0
  246.          sqi = 0.0d0
  247. c     .......... compute diagonal elements that would result
  248. c                if transformations were applied ..........
  249.   500    ssr = sqr * szr + sqi * szi
  250.          ssi = sqr * szi - sqi * szr
  251.          i = 1
  252.          tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21
  253.      x      + ssr * a22
  254.          ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22
  255.          dr = cq * cz * b11 + cq * szr * b12 + ssr * b22
  256.          di = cq * szi * b12 + ssi * b22
  257.          go to 503
  258.   502    i = 2
  259.          tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21
  260.      x      + cq * cz * a22
  261.          ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21
  262.          dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22
  263.          di = -ssi * b11 - sqi * cz * b12
  264.   503    t = ti * dr - tr * di
  265.          j = na
  266.          if (t .lt. 0.0d0) j = en
  267.          r = dsqrt(dr*dr+di*di)
  268.          beta(j) = bn * r
  269.          alfr(j) = an * (tr * dr + ti * di) / r
  270.          alfi(j) = an * t / r
  271.          if (i .eq. 1) go to 502
  272.   505    isw = 3 - isw
  273.   510 continue
  274.       b(n,1) = epsb
  275. c
  276.       return
  277.       end
  278.