      subroutine qzvec(nm,n,a,b,alfr,alfi,beta,z)
c
      integer i,j,k,m,n,en,ii,jj,na,nm,nn,isw,enm2
      double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n)
      double precision d,q,r,s,t,w,x,y,di,dr,ra,rr,sa,ti,tr,t1,t2,w1,x1,
     x       zz,z1,alfm,almi,almr,betm,epsb
c
c     this subroutine is the optional fourth step of the qz algorithm
c     for solving generalized matrix eigenvalue problems,
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
c
c     this subroutine accepts a pair of real matrices, one of them in
c     quasi-triangular form (in which each 2-by-2 block corresponds to
c     a pair of complex eigenvalues) and the other in upper triangular
c     form.  it computes the eigenvectors of the triangular problem and
c     transforms the results back to the original coordinate system.
c     it is usually preceded by  qzhes,  qzit, and  qzval.
c
c     on input
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c
c        n is the order of the matrices.
c
c        a contains a real upper quasi-triangular matrix.
c
c        b contains a real upper triangular matrix.  in addition,
c          location b(n,1) contains the tolerance quantity (epsb)
c          computed and saved in  qzit.
c
c        alfr, alfi, and beta  are vectors with components whose
c          ratios ((alfr+i*alfi)/beta) are the generalized
c          eigenvalues.  they are usually obtained from  qzval.
c
c        z contains the transformation matrix produced in the
c          reductions by  qzhes,  qzit, and  qzval, if performed.
c          if the eigenvectors of the triangular problem are
c          desired, z must contain the identity matrix.
c
c     on output
c
c        a is unaltered.  its subdiagonal elements provide information
c           about the storage of the complex eigenvectors.
c
c        b has been destroyed.
c
c        alfr, alfi, and beta are unaltered.
c
c        z contains the real and imaginary parts of the eigenvectors.
c          if alfi(i) .eq. 0.0, the i-th eigenvalue is real and
c            the i-th column of z contains its eigenvector.
c          if alfi(i) .ne. 0.0, the i-th eigenvalue is complex.
c            if alfi(i) .gt. 0.0, the eigenvalue is the first of
c              a complex pair and the i-th and (i+1)-th columns
c              of z contain its eigenvector.
c            if alfi(i) .lt. 0.0, the eigenvalue is the second of
c              a complex pair and the (i-1)-th and i-th columns
c              of z contain the conjugate of its eigenvector.
c          each eigenvector is normalized so that the modulus
c          of its largest component is 1.0 .
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      epsb = b(n,1)
      isw = 1
c     .......... for en=n step -1 until 1 do -- ..........
      do 800 nn = 1, n
         en = n + 1 - nn
         na = en - 1
         if (isw .eq. 2) go to 795
         if (alfi(en) .ne. 0.0d0) go to 710
c     .......... real vector ..........
         m = en
         b(en,en) = 1.0d0
         if (na .eq. 0) go to 800
         alfm = alfr(m)
         betm = beta(m)
c     .......... for i=en-1 step -1 until 1 do -- ..........
         do 700 ii = 1, na
            i = en - ii
            w = betm * a(i,i) - alfm * b(i,i)
            r = 0.0d0
c
            do 610 j = m, en
  610       r = r + (betm * a(i,j) - alfm * b(i,j)) * b(j,en)
c
            if (i .eq. 1 .or. isw .eq. 2) go to 630
            if (betm * a(i,i-1) .eq. 0.0d0) go to 630
            zz = w
            s = r
            go to 690
  630       m = i
            if (isw .eq. 2) go to 640
c     .......... real 1-by-1 block ..........
            t = w
            if (w .eq. 0.0d0) t = epsb
            b(i,en) = -r / t
            go to 700
c     .......... real 2-by-2 block ..........
  640       x = betm * a(i,i+1) - alfm * b(i,i+1)
            y = betm * a(i+1,i)
            q = w * zz - x * y
            t = (x * s - zz * r) / q
            b(i,en) = t
            if (dabs(x) .le. dabs(zz)) go to 650
            b(i+1,en) = (-r - w * t) / x
            go to 690
  650       b(i+1,en) = (-s - y * t) / zz
  690       isw = 3 - isw
  700    continue
c     .......... end real vector ..........
         go to 800
c     .......... complex vector ..........
  710    m = na
         almr = alfr(m)
         almi = alfi(m)
         betm = beta(m)
c     .......... last vector component chosen imaginary so that
c                eigenvector matrix is triangular ..........
         y = betm * a(en,na)
         b(na,na) = -almi * b(en,en) / y
         b(na,en) = (almr * b(en,en) - betm * a(en,en)) / y
         b(en,na) = 0.0d0
         b(en,en) = 1.0d0
         enm2 = na - 1
         if (enm2 .eq. 0) go to 795
c     .......... for i=en-2 step -1 until 1 do -- ..........
         do 790 ii = 1, enm2
            i = na - ii
            w = betm * a(i,i) - almr * b(i,i)
            w1 = -almi * b(i,i)
            ra = 0.0d0
            sa = 0.0d0
c
            do 760 j = m, en
               x = betm * a(i,j) - almr * b(i,j)
               x1 = -almi * b(i,j)
               ra = ra + x * b(j,na) - x1 * b(j,en)
               sa = sa + x * b(j,en) + x1 * b(j,na)
  760       continue
c
            if (i .eq. 1 .or. isw .eq. 2) go to 770
            if (betm * a(i,i-1) .eq. 0.0d0) go to 770
            zz = w
            z1 = w1
            r = ra
            s = sa
            isw = 2
            go to 790
  770       m = i
            if (isw .eq. 2) go to 780
c     .......... complex 1-by-1 block ..........
            tr = -ra
            ti = -sa
  773       dr = w
            di = w1
c     .......... complex divide (t1,t2) = (tr,ti) / (dr,di) ..........
  775       if (dabs(di) .gt. dabs(dr)) go to 777
            rr = di / dr
            d = dr + di * rr
            t1 = (tr + ti * rr) / d
            t2 = (ti - tr * rr) / d
            go to (787,782), isw
  777       rr = dr / di
            d = dr * rr + di
            t1 = (tr * rr + ti) / d
            t2 = (ti * rr - tr) / d
            go to (787,782), isw
c     .......... complex 2-by-2 block ..........
  780       x = betm * a(i,i+1) - almr * b(i,i+1)
            x1 = -almi * b(i,i+1)
            y = betm * a(i+1,i)
            tr = y * ra - w * r + w1 * s
            ti = y * sa - w * s - w1 * r
            dr = w * zz - w1 * z1 - x * y
            di = w * z1 + w1 * zz - x1 * y
            if (dr .eq. 0.0d0 .and. di .eq. 0.0d0) dr = epsb
            go to 775
  782       b(i+1,na) = t1
            b(i+1,en) = t2
            isw = 1
            if (dabs(y) .gt. dabs(w) + dabs(w1)) go to 785
            tr = -ra - x * b(i+1,na) + x1 * b(i+1,en)
            ti = -sa - x * b(i+1,en) - x1 * b(i+1,na)
            go to 773
  785       t1 = (-r - zz * b(i+1,na) + z1 * b(i+1,en)) / y
            t2 = (-s - zz * b(i+1,en) - z1 * b(i+1,na)) / y
  787       b(i,na) = t1
            b(i,en) = t2
  790    continue
c     .......... end complex vector ..........
  795    isw = 3 - isw
  800 continue
c     .......... end back substitution.
c                transform to original coordinate system.
c                for j=n step -1 until 1 do -- ..........
      do 880 jj = 1, n
         j = n + 1 - jj
c
         do 880 i = 1, n
            zz = 0.0d0
c
            do 860 k = 1, j
  860       zz = zz + z(i,k) * b(k,j)
c
            z(i,j) = zz
  880 continue
c     .......... normalize so that modulus of largest
c                component of each vector is 1.
c                (isw is 1 initially from before) ..........
      do 950 j = 1, n
         d = 0.0d0
         if (isw .eq. 2) go to 920
         if (alfi(j) .ne. 0.0d0) go to 945
c
         do 890 i = 1, n
            if (dabs(z(i,j)) .gt. d) d = dabs(z(i,j))
  890    continue
c
         do 900 i = 1, n
  900    z(i,j) = z(i,j) / d
c
         go to 950
c
  920    do 930 i = 1, n
            r = dabs(z(i,j-1)) + dabs(z(i,j))
            if (r .ne. 0.0d0) r = r * dsqrt((z(i,j-1)/r)**2
     x                                     +(z(i,j)/r)**2)
            if (r .gt. d) d = r
  930    continue
c
         do 940 i = 1, n
            z(i,j-1) = z(i,j-1) / d
            z(i,j) = z(i,j) / d
  940    continue
c
  945    isw = 3 - isw
  950 continue
c
      return
      end
