C/ FIGURE 3.6.1
      SUBROUTINE DEGENP(A,B,N,LAM0,POLY)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                              DECLARATIONS FOR ARGUMENTS
      DOUBLE PRECISION A(N,N),B(N,N),LAM0,POLY(N+1)
      INTEGER N
C                              DECLARATIONS FOR LOCAL VARIABLES
      DOUBLE PRECISION DET(N+1,N+1)
C
C  SUBROUTINE DEGENP CALLS HESSQZ TO REDUCE THE GENERALIZED EIGENVALUE 
C     PROBLEM
C
C                A*X = LAMBDA*B*X
C
C  TO A SIMILAR PROBLEM WITH THE SAME EIGENVALUES, WHERE A IS UPPER
C  HESSENBERG AND B IS UPPER TRIANGULAR, THEN CALCULATES THE COEFFICIENTS 
C  OF THE POLYNOMIAL DET(A-LAMBDA*B).  THE ROOTS OF THIS POLYNOMIAL WILL 
C  BE THE EIGENVALUES OF THE GENERALIZED PROBLEM. 
C
C  ARGUMENTS
C
C             ON INPUT                          ON OUTPUT
C             --------                          ---------
C
C    A      - THE N BY N A MATRIX.              A IS UPPER HESSENBERG. 
C
C    B      - THE N BY N B MATRIX.              B IS UPPER TRIANGULAR.  
C
C    N      - THE SIZE OF MATRICES A AND B.
C
C    LAM0   - A SCALAR, SEE POLY. (SET 
C             LAM0=0 TO GET THE USUAL
C             POLYNOMIAL COEFFICIENTS)
C
C    POLY   -                                   VECTOR OF LENGTH N+1,
C                                               CONTAINING POLYNOMIAL
C                                               COEFFICIENTS. 
C                                               DET(A-LAMBDA*B)=
C                                               SUM FROM I=0 TO N OF
C                                               POLY(I+1)*
C                                                   (LAMBDA-LAM0)**I
C
C-----------------------------------------------------------------------
C
C                        CALL HESSQZ TO REDUCE A TO UPPER HESSENBERG AND
C                        B TO UPPER TRIANGULAR FORM
      CALL HESSQZ(A,B,N)
C                        DEGENP USES A RECURRENCE RELATION TO CALCULATE
C                        DET(A-LAMBDA*B) AND ALL N DERIVATIVES AT 
C                        LAMBDA = LAM0, FROM WHICH THE POLYNOMIAL
C                        COEFFICIENTS CAN BE FOUND. 
      DET(1,N+1) = 1.0
      DO 5 I=2,N+1
         DET(I,N+1) = 0.0
    5 CONTINUE
      DET(1,N) = A(N,N)-LAM0*B(N,N)
      DET(2,N) = -B(N,N)
      DO 10 I=3,N+1
        DET(I,N) = 0.0
   10 CONTINUE
      DO 30 K=N-1,1,-1
         DET(1,K) = (A(K,K)-LAM0*B(K,K))*DET(1,K+1)
         DO 15 I=1,N
            DET(I+1,K) = (A(K,K)-LAM0*B(K,K))*DET(I+1,K+1)
     &                               - B(K,K)*DET(I,K+1)
   15    CONTINUE
         FACT = 1.0
         DO 25 J=K+1,N
            FACT = -FACT*A(J,J-1)
            IF (A(K,J).EQ.0.0 .AND. B(K,J).EQ.0.0) GO TO 25
            DET(1,K) = DET(1,K) + FACT*(A(K,J)-LAM0*B(K,J))*DET(1,J+1)
            DO 20 I=1,N
               DET(I+1,K) = DET(I+1,K) + FACT*
     &         ((A(K,J)-LAM0*B(K,J))*DET(I+1,J+1) - B(K,J)*DET(I,J+1))
   20       CONTINUE
   25    CONTINUE
   30 CONTINUE
      DO 35 I=1,N+1
         POLY(I) = DET(I,1)
   35 CONTINUE
      RETURN
      END


      SUBROUTINE HESSQZ(A,B,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                              DECLARATIONS FOR ARGUMENTS
      DOUBLE PRECISION A(N,N),B(N,N)
      INTEGER N
C
C  SUBROUTINE HESSQZ REDUCES THE GENERALIZED EIGENVALUE PROBLEM
C
C                A*X = LAMBDA*B*X
C
C  TO A SIMILAR PROBLEM WITH THE SAME EIGENVALUES, WHERE A IS 
C  UPPER HESSENBERG AND B IS UPPER TRIANGULAR. 
C
C  ARGUMENTS
C
C             ON INPUT                          ON OUTPUT
C             --------                          ---------
C
C    A      - THE N BY N A MATRIX.              A IS UPPER HESSENBERG. 
C
C    B      - THE N BY N B MATRIX.              B IS UPPER TRIANGULAR.  
C
C    N      - THE SIZE OF MATRICES A AND B.
C
C-----------------------------------------------------------------------
C                              PREMULTIPLY A AND B BY ORTHOGONAL MATRIX
C                              (PRODUCT OF GIVENS MATRICES) Q, SUCH   
C                              THAT QB IS UPPER TRIANGULAR.
      DO 20 I=1,N-1
         DO 15 J=I+1,N
            IF (B(J,I).EQ.0.0) GO TO 15
            DEN = SQRT(B(I,I)**2+B(J,I)**2)
            S = -B(J,I)/DEN
            C =  B(I,I)/DEN
            DO 5 K=I,N
               BIK = C*B(I,K)-S*B(J,K) 
               BJK = S*B(I,K)+C*B(J,K)
               B(I,K) = BIK
               B(J,K) = BJK
    5       CONTINUE
            DO 10 K=1,N
               AIK = C*A(I,K)-S*A(J,K) 
               AJK = S*A(I,K)+C*A(J,K)
               A(I,K) = AIK
               A(J,K) = AJK
   10       CONTINUE
   15    CONTINUE
   20 CONTINUE
C                              PREMULTIPLY A AND B BY ORTHOGONAL MATRIX
C                              Q, AND POSTMULTIPLY BY ORTHOGONAL MATRIX
C                              Z, SUCH THAT QAZ IS UPPER HESSENBERG AND  
C                              QBZ IS STILL UPPER TRIANGULAR      
      DO 50 I=1,N-2
         DO 45 J=N,I+2,-1
            IF (A(J,I).EQ.0.0) GO TO 45
C                              PREMULTIPLY A TO ZERO A(J,I)   
            DEN = SQRT(A(J-1,I)**2+A(J,I)**2)
            S = -A(J,I)/DEN 
            C =  A(J-1,I)/DEN
            DO 25 K=I,N
               A1K = C*A(J-1,K) - S*A(J,K)
               A2K = S*A(J-1,K) + C*A(J,K)
               A(J-1,K) = A1K
               A(J,K) = A2K
   25       CONTINUE
C                              PREMULTIPLY B BY SAME MATRIX, CREATING 
C                              NEW NONZERO B(J,J-1)
            DO 30 K=J-1,N
               B1K = C*B(J-1,K) - S*B(J,K)
               B2K = S*B(J-1,K) + C*B(J,K)
               B(J-1,K) = B1K
               B(J,K) = B2K
   30       CONTINUE
            IF (B(J,J-1).EQ.0.0) GO TO 45
C                              POSTMULTIPLY B TO ZERO B(J,J-1)     
            DEN = SQRT(B(J,J-1)**2+B(J,J)**2)
            S = -B(J,J-1)/DEN
            C =  B(J,J)/DEN
            DO 35 K=1,J
               BK1 = C*B(K,J-1) + S*B(K,J)
               BK2 = -S*B(K,J-1) + C*B(K,J)  
               B(K,J-1) = BK1
               B(K,J) = BK2
   35       CONTINUE
C                              POSTMULTIPLY A BY SAME MATRIX
            DO 40 K=1,N
               AK1 = C*A(K,J-1) + S*A(K,J)
               AK2 = -S*A(K,J-1) + C*A(K,J)  
               A(K,J-1) = AK1
               A(K,J) = AK2
   40       CONTINUE
   45    CONTINUE
   50 CONTINUE
      RETURN
      END


      SUBROUTINE DEGEN(A,B,N,EIG,NEIG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                              DECLARATIONS FOR ARGUMENTS
      COMPLEX*16 EIG(N)
      DOUBLE PRECISION A(N,N),B(N,N)
      INTEGER N,NEIG
C                              DECLARATIONS FOR LOCAL VARIABLES
      DOUBLE PRECISION ALFR(N),ALFI(N),BETA(N)
C
C  SUBROUTINE DEGEN FINDS ALL EIGENVALUES OF THE GENERALIZED 
C  EIGENVALUE PROBLEM
C
C                A*X = LAMBDA*B*X
C
C  ARGUMENTS
C
C             ON INPUT                          ON OUTPUT
C             --------                          ---------
C
C    A      - THE N BY N A MATRIX.              DESTROYED. 
C
C    B      - THE N BY N B MATRIX.              DESTROYED.  
C
C    N      - THE SIZE OF MATRICES A AND B.
C
C    EIG    -                                   A COMPLEX N-VECTOR
C                                               CONTAINING THE  
C                                               EIGENVALUES IN ITS
C                                               FIRST NEIG ELEMENTS
C
C    NEIG   -                                   NUMBER OF EIGENVALUES
C                                               (NEIG .LE. N)
C
C-----------------------------------------------------------------------
C                              EPS = MACHINE FLOATING POINT RELATIVE
C                                    PRECISION
C *****************************
      DATA EPS/2.D-16/
C *****************************
      EPS1 = 100*EPS
C                              REDUCE A TO UPPER HESSENBERG, B TO
C                              UPPER TRIANGULAR FORM
      CALL HESSQZ(A,B,N)
C                              REDUCE A TO QUASITRIANGULAR FORM
      CALL QZIT(N,A,B,EPS1,IERR)
      IF (IERR.NE.0) GO TO 20
C                              EXTRACT EIGENVALUES
      CALL QZVAL(N,A,B,ALFR,ALFI,BETA)
      NEIG = 0
      EPSB = B(N,1)
      DO 10 I=1,N
         IF (BETA(I).GT.EPSB) THEN
            NEIG = NEIG+1
            EIG(NEIG) = (ALFR(I) + ALFI(I)*CMPLX(0.0,1.0))/BETA(I)
         ENDIF 
   10 CONTINUE 
      RETURN
C                             QZ ITERATION DOES NOT CONVERGE
   20 PRINT 30
   30 FORMAT (' ***** QZ ITERATION DOES NOT CONVERGE *****')    
      RETURN
      END


      subroutine qzit(n,a,b,eps1,ierr)
c
      integer i,j,k,l,n,en,k1,k2,ld,ll,l1,na,ish,itn,its,km1,lm1,
     x        enm2,ierr,lor1,enorn
      double precision a(n,n),b(n,n)
      double precision r,s,t,a1,a2,a3,ep,sh,u1,u2,u3,v1,v2,v3,ani,a11,
     x       a12,a21,a22,a33,a34,a43,a44,bni,b11,b12,b22,b33,b34,
     x       b44,epsa,epsb,eps1,anorm,bnorm
      logical notlas
c
c     this subroutine is the second 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     as modified in technical note nasa tn d-7305(1973) by ward.
c
c     this subroutine accepts a pair of real matrices, one of them
c     in upper hessenberg form and the other in upper triangular form.
c     it reduces the hessenberg matrix to quasi-triangular form using
c     orthogonal transformations while maintaining the triangular form
c     of the other matrix. 
c
c     on input
c
c        n is the order of the matrices.
c
c        a contains a real upper hessenberg matrix.
c
c        b contains a real upper triangular matrix.
c
c        eps1 is a tolerance used to determine negligible elements.
c
c     on output
c
c        a has been reduced to quasi-triangular form.  the elements
c          below the first subdiagonal are still zero and no two
c          consecutive subdiagonal elements are nonzero.
c
c        b is still in upper triangular form, although its elements
c          have been altered.  the location b(n,1) is used to store
c          eps1 times the norm of b for later use by  qzval  and  qzvec.
c
c        ierr is set to
c          zero       for normal return,
c          j          if the limit of 30*n iterations is exhausted
c                     while the j-th eigenvalue is being sought.
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     ------------------------------------------------------------------
      ierr = 0
c     .......... compute epsa,epsb ..........
      anorm = 0.0d0
      bnorm = 0.0d0
c
      do 30 i = 1, n
         ani = 0.0d0
         if (i .ne. 1) ani = abs(a(i,i-1))
         bni = 0.0d0
c
         do 20 j = i, n
            ani = ani + abs(a(i,j))
            bni = bni + abs(b(i,j))
   20    continue
c
         if (ani .gt. anorm) anorm = ani
         if (bni .gt. bnorm) bnorm = bni
   30 continue
c
      if (anorm .eq. 0.0d0) anorm = 1.0d0
      if (bnorm .eq. 0.0d0) bnorm = 1.0d0
      ep = eps1
   50 epsa = ep * anorm
      epsb = ep * bnorm
c     .......... reduce a to quasi-triangular form, while
c                keeping b triangular ..........
      lor1 = 1
      enorn = n
      en = n
      itn = 30*n
c     .......... begin qz step ..........
   60 if (en .le. 2) go to 1001
      enorn = en
      its = 0
      na = en - 1
      enm2 = na - 1
   70 ish = 2
c     .......... check for convergence or reducibility.
c                for l=en step -1 until 1 do -- ..........
      do 80 ll = 1, en
         lm1 = en - ll
         l = lm1 + 1
         if (l .eq. 1) go to 95
         if (abs(a(l,lm1)) .le. epsa) go to 90
   80 continue
c
   90 a(l,lm1) = 0.0d0
      if (l .lt. na) go to 95
c     .......... 1-by-1 or 2-by-2 block isolated ..........
      en = lm1
      go to 60
c     .......... check for small top of b ..........
   95 ld = l
  100 l1 = l + 1
      b11 = b(l,l)
      if (abs(b11) .gt. epsb) go to 120
      b(l,l) = 0.0d0
      s = abs(a(l,l)) + abs(a(l1,l))
      u1 = a(l,l) / s
      u2 = a(l1,l) / s
      r = sign(sqrt(u1*u1+u2*u2),u1)
      v1 = -(u1 + r) / r
      v2 = -u2 / r
      u2 = v2 / v1
c
      do 110 j = l, enorn
         t = a(l,j) + u2 * a(l1,j)
         a(l,j) = a(l,j) + t * v1
         a(l1,j) = a(l1,j) + t * v2
         t = b(l,j) + u2 * b(l1,j)
         b(l,j) = b(l,j) + t * v1
         b(l1,j) = b(l1,j) + t * v2
  110 continue
c
      if (l .ne. 1) a(l,lm1) = -a(l,lm1)
      lm1 = l
      l = l1
      go to 90
  120 a11 = a(l,l) / b11
      a21 = a(l1,l) / b11
      if (ish .eq. 1) go to 140
c     .......... iteration strategy ..........
      if (itn .eq. 0) go to 1000
      if (its .eq. 10) go to 155
c     .......... determine type of shift ..........
      b22 = b(l1,l1)
      if (abs(b22) .lt. epsb) b22 = epsb
      b33 = b(na,na)
      if (abs(b33) .lt. epsb) b33 = epsb
      b44 = b(en,en)
      if (abs(b44) .lt. epsb) b44 = epsb
      a33 = a(na,na) / b33
      a34 = a(na,en) / b44
      a43 = a(en,na) / b33
      a44 = a(en,en) / b44
      b34 = b(na,en) / b44
      t = 0.5d0 * (a43 * b34 - a33 - a44)
      r = t * t + a34 * a43 - a33 * a44
      if (r .lt. 0.0d0) go to 150
c     .......... determine single shift zeroth column of a ..........
      ish = 1
      r = sqrt(r)
      sh = -t + r
      s = -t - r
      if (abs(s-a44) .lt. abs(sh-a44)) sh = s
c     .......... look for two consecutive small
c                sub-diagonal elements of a.
c                for l=en-2 step -1 until ld do -- ..........
      do 130 ll = ld, enm2
         l = enm2 + ld - ll
         if (l .eq. ld) go to 140
         lm1 = l - 1
         l1 = l + 1
         t = a(l,l)
         if (abs(b(l,l)) .gt. epsb) t = t - sh * b(l,l)
         if (abs(a(l,lm1)) .le. abs(t/a(l1,l)) * epsa) go to 100
  130 continue
c
  140 a1 = a11 - sh
      a2 = a21
      if (l .ne. ld) a(l,lm1) = -a(l,lm1)
      go to 160
c     .......... determine double shift zeroth column of a ..........
  150 a12 = a(l,l1) / b22
      a22 = a(l1,l1) / b22
      b12 = b(l,l1) / b22
      a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11)
     x     / a21 + a12 - a11 * b12
      a2 = (a22 - a11) - a21 * b12 - (a33 - a11) - (a44 - a11)
     x     + a43 * b34
      a3 = a(l1+1,l1) / b22
      go to 160
c     .......... ad hoc shift ..........
  155 a1 = 0.0d0
      a2 = 1.0d0
      a3 = 1.1605d0
  160 its = its + 1
      itn = itn - 1
      lor1 = ld
c     .......... main loop ..........
      do 260 k = l, na
         notlas = k .ne. na .and. ish .eq. 2
         k1 = k + 1
         k2 = k + 2
         km1 = max(k-1,l)
         ll = min(en,k1+ish)
         if (notlas) go to 190
c     .......... zero a(k+1,k-1) ..........
         if (k .eq. l) go to 170
         a1 = a(k,km1)
         a2 = a(k1,km1)
  170    s = abs(a1) + abs(a2)
         if (s .eq. 0.0d0) go to 70
         u1 = a1 / s
         u2 = a2 / s
         r = sign(sqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 180 j = km1, enorn
            t = a(k,j) + u2 * a(k1,j)
            a(k,j) = a(k,j) + t * v1
            a(k1,j) = a(k1,j) + t * v2
            t = b(k,j) + u2 * b(k1,j)
            b(k,j) = b(k,j) + t * v1
            b(k1,j) = b(k1,j) + t * v2
  180    continue
c
         if (k .ne. l) a(k1,km1) = 0.0d0
         go to 240
c     .......... zero a(k+1,k-1) and a(k+2,k-1) ..........
  190    if (k .eq. l) go to 200
         a1 = a(k,km1)
         a2 = a(k1,km1)
         a3 = a(k2,km1)
  200    s = abs(a1) + abs(a2) + abs(a3)
         if (s .eq. 0.0d0) go to 260
         u1 = a1 / s
         u2 = a2 / s
         u3 = a3 / s
         r = sign(sqrt(u1*u1+u2*u2+u3*u3),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         v3 = -u3 / r
         u2 = v2 / v1
         u3 = v3 / v1
c
         do 210 j = km1, enorn
            t = a(k,j) + u2 * a(k1,j) + u3 * a(k2,j)
            a(k,j) = a(k,j) + t * v1
            a(k1,j) = a(k1,j) + t * v2
            a(k2,j) = a(k2,j) + t * v3
            t = b(k,j) + u2 * b(k1,j) + u3 * b(k2,j)
            b(k,j) = b(k,j) + t * v1
            b(k1,j) = b(k1,j) + t * v2
            b(k2,j) = b(k2,j) + t * v3
  210    continue
c
         if (k .eq. l) go to 220
         a(k1,km1) = 0.0d0
         a(k2,km1) = 0.0d0
c     .......... zero b(k+2,k+1) and b(k+2,k) ..........
  220    s = abs(b(k2,k2)) + abs(b(k2,k1)) + abs(b(k2,k))
         if (s .eq. 0.0d0) go to 240
         u1 = b(k2,k2) / s
         u2 = b(k2,k1) / s
         u3 = b(k2,k) / s
         r = sign(sqrt(u1*u1+u2*u2+u3*u3),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         v3 = -u3 / r
         u2 = v2 / v1
         u3 = v3 / v1
c
         do 230 i = lor1, ll
            t = a(i,k2) + u2 * a(i,k1) + u3 * a(i,k)
            a(i,k2) = a(i,k2) + t * v1
            a(i,k1) = a(i,k1) + t * v2
            a(i,k) = a(i,k) + t * v3
            t = b(i,k2) + u2 * b(i,k1) + u3 * b(i,k)
            b(i,k2) = b(i,k2) + t * v1
            b(i,k1) = b(i,k1) + t * v2
            b(i,k) = b(i,k) + t * v3
  230    continue
c
         b(k2,k) = 0.0d0
         b(k2,k1) = 0.0d0
c     .......... zero b(k+1,k) ..........
  240    s = abs(b(k1,k1)) + abs(b(k1,k))
         if (s .eq. 0.0d0) go to 260
         u1 = b(k1,k1) / s
         u2 = b(k1,k) / s
         r = sign(sqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 250 i = lor1, ll
            t = a(i,k1) + u2 * a(i,k)
            a(i,k1) = a(i,k1) + t * v1
            a(i,k) = a(i,k) + t * v2
            t = b(i,k1) + u2 * b(i,k)
            b(i,k1) = b(i,k1) + t * v1
            b(i,k) = b(i,k) + t * v2
  250    continue
c
         b(k1,k) = 0.0d0
  260 continue
c     .......... end qz step ..........
      go to 70
c     .......... set error -- all eigenvalues have not
c                converged after 30*n iterations ..........
 1000 ierr = en
c     .......... save epsb for use by qzval and qzvec ..........
 1001 if (n .gt. 1) b(n,1) = epsb
      return
      end


      subroutine qzval(n,a,b,alfr,alfi,beta)
c
      integer i,j,n,en,na,nn,isw
      double precision a(n,n),b(n,n),alfr(n),alfi(n),beta(n)
      double precision c,d,e,r,s,t,an,a1,a2,bn,cq,cz,di,dr,ei,ti,tr,u1,
     x       u2,v1,v2,a1i,a11,a12,a2i,a21,a22,b11,b12,b22,sqi,sqr,
     x       ssi,ssr,szi,szr,a11i,a11r,a12i,a12r,a22i,a22r,epsb
c
c     this subroutine is the third 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
c     in quasi-triangular form and the other in upper triangular form.
c     it reduces the quasi-triangular matrix further, so that any
c     remaining 2-by-2 blocks correspond to pairs of complex
c     eigenvalues, and returns quantities whose ratios give the
c     generalized eigenvalues.  
c
c     on input
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     on output
c
c        a has been reduced further to a quasi-triangular matrix
c          in which all nonzero subdiagonal elements correspond to
c          pairs of complex eigenvalues.
c
c        b is still in upper triangular form, although its elements
c          have been altered.  b(n,1) is unaltered.
c
c        alfr and alfi contain the real and imaginary parts of the
c          diagonal elements of the triangular matrix that would be
c          obtained if a were reduced completely to triangular form
c          by unitary transformations.  non-zero values of alfi occur
c          in pairs, the first member positive and the second negative.
c
c        beta contains the diagonal elements of the corresponding b,
c          normalized to be real and non-negative.  the generalized
c          eigenvalues are then the ratios ((alfr+i*alfi)/beta).
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     .......... find eigenvalues of quasi-triangular matrices.
c                for en=n step -1 until 1 do -- ..........
      do 510 nn = 1, n
         en = n + 1 - nn
         na = en - 1
         if (isw .eq. 2) go to 505
         if (en .eq. 1) go to 410
         if (a(en,na) .ne. 0.0d0) go to 420
c     .......... 1-by-1 block, one real root ..........
  410    alfr(en) = a(en,en)
         if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en)
         beta(en) = abs(b(en,en))
         alfi(en) = 0.0d0
         go to 510
c     .......... 2-by-2 block ..........
  420    if (abs(b(na,na)) .le. epsb) go to 455
         if (abs(b(en,en)) .gt. epsb) go to 430
         a1 = a(en,en)
         a2 = a(en,na)
         bn = 0.0d0
         go to 435
  430    an = abs(a(na,na)) + abs(a(na,en)) + abs(a(en,na))
     x      + abs(a(en,en))
         bn = abs(b(na,na)) + abs(b(na,en)) + abs(b(en,en))
         a11 = a(na,na) / an
         a12 = a(na,en) / an
         a21 = a(en,na) / an
         a22 = a(en,en) / an
         b11 = b(na,na) / bn
         b12 = b(na,en) / bn
         b22 = b(en,en) / bn
         e = a11 / b11
         ei = a22 / b22
         s = a21 / (b11 * b22)
         t = (a22 - e * b22) / b22
         if (abs(e) .le. abs(ei)) go to 431
         e = ei
         t = (a11 - e * b11) / b11
  431    c = 0.5d0 * (t - s * b12)
         d = c * c + s * (a12 - e * b12)
         if (d .lt. 0.0d0) go to 480
c     .......... two real roots.
c                zero both a(en,na) and b(en,na) ..........
         e = e + (c + sign(sqrt(d),c))
         a11 = a11 - e * b11
         a12 = a12 - e * b12
         a22 = a22 - e * b22
         if (abs(a11) + abs(a12) .lt.
     x       abs(a21) + abs(a22)) go to 432
         a1 = a12
         a2 = a11
         go to 435
  432    a1 = a22
         a2 = a21
c     .......... choose and apply real z ..........
  435    s = abs(a1) + abs(a2)
         u1 = a1 / s
         u2 = a2 / s
         r = sign(sqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 440 i = 1, en
            t = a(i,en) + u2 * a(i,na)
            a(i,en) = a(i,en) + t * v1
            a(i,na) = a(i,na) + t * v2
            t = b(i,en) + u2 * b(i,na)
            b(i,en) = b(i,en) + t * v1
            b(i,na) = b(i,na) + t * v2
  440    continue
c
         if (bn .eq. 0.0d0) go to 475
         if (an .lt. abs(e) * bn) go to 455
         a1 = b(na,na)
         a2 = b(en,na)
         go to 460
  455    a1 = a(na,na)
         a2 = a(en,na)
c     .......... choose and apply real q ..........
  460    s = abs(a1) + abs(a2)
         if (s .eq. 0.0d0) go to 475
         u1 = a1 / s
         u2 = a2 / s
         r = sign(sqrt(u1*u1+u2*u2),u1)
         v1 = -(u1 + r) / r
         v2 = -u2 / r
         u2 = v2 / v1
c
         do 470 j = na, n
            t = a(na,j) + u2 * a(en,j)
            a(na,j) = a(na,j) + t * v1
            a(en,j) = a(en,j) + t * v2
            t = b(na,j) + u2 * b(en,j)
            b(na,j) = b(na,j) + t * v1
            b(en,j) = b(en,j) + t * v2
  470    continue
c
  475    a(en,na) = 0.0d0
         b(en,na) = 0.0d0
         alfr(na) = a(na,na)
         alfr(en) = a(en,en)
         if (b(na,na) .lt. 0.0d0) alfr(na) = -alfr(na)
         if (b(en,en) .lt. 0.0d0) alfr(en) = -alfr(en)
         beta(na) = abs(b(na,na))
         beta(en) = abs(b(en,en))
         alfi(en) = 0.0d0
         alfi(na) = 0.0d0
         go to 505
c     .......... two complex roots ..........
  480    e = e + c
         ei = sqrt(-d)
         a11r = a11 - e * b11
         a11i = ei * b11
         a12r = a12 - e * b12
         a12i = ei * b12
         a22r = a22 - e * b22
         a22i = ei * b22
         if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) .lt.
     x       abs(a21) + abs(a22r) + abs(a22i)) go to 482
         a1 = a12r
         a1i = a12i
         a2 = -a11r
         a2i = -a11i
         go to 485
  482    a1 = a22r
         a1i = a22i
         a2 = -a21
         a2i = 0.0d0
c     .......... choose complex z ..........
  485    cz = sqrt(a1*a1+a1i*a1i)
         if (cz .eq. 0.0d0) go to 487
         szr = (a1 * a2 + a1i * a2i) / cz
         szi = (a1 * a2i - a1i * a2) / cz
         r = sqrt(cz*cz+szr*szr+szi*szi)
         cz = cz / r
         szr = szr / r
         szi = szi / r
         go to 490
  487    szr = 1.0d0
         szi = 0.0d0
  490    if (an .lt. (abs(e) + ei) * bn) go to 492
         a1 = cz * b11 + szr * b12
         a1i = szi * b12
         a2 = szr * b22
         a2i = szi * b22
         go to 495
  492    a1 = cz * a11 + szr * a12
         a1i = szi * a12
         a2 = cz * a21 + szr * a22
         a2i = szi * a22
c     .......... choose complex q ..........
  495    cq = sqrt(a1*a1+a1i*a1i)
         if (cq .eq. 0.0d0) go to 497
         sqr = (a1 * a2 + a1i * a2i) / cq
         sqi = (a1 * a2i - a1i * a2) / cq
         r = sqrt(cq*cq+sqr*sqr+sqi*sqi)
         cq = cq / r
         sqr = sqr / r
         sqi = sqi / r
         go to 500
  497    sqr = 1.0d0
         sqi = 0.0d0
c     .......... compute diagonal elements that would result
c                if transformations were applied ..........
  500    ssr = sqr * szr + sqi * szi
         ssi = sqr * szi - sqi * szr
         i = 1
         tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21
     x      + ssr * a22
         ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22
         dr = cq * cz * b11 + cq * szr * b12 + ssr * b22
         di = cq * szi * b12 + ssi * b22
         go to 503
  502    i = 2
         tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21
     x      + cq * cz * a22
         ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21
         dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22
         di = -ssi * b11 - sqi * cz * b12
  503    t = ti * dr - tr * di
         j = na
         if (t .lt. 0.0d0) j = en
         r = sqrt(dr*dr+di*di)
         beta(j) = bn * r
         alfr(j) = an * (tr * dr + ti * di) / r
         alfi(j) = an * t / r
         if (i .eq. 1) go to 502
  505    isw = 3 - isw
  510 continue
      b(n,1) = epsb
c
      return
      end

