dtgevc.f

Go to the documentation of this file.

*> \brief \b DTGEVC
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DTGEVC + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgevc.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgevc.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgevc.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
*                          LDVL, VR, LDVR, MM, M, WORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          HOWMNY, SIDE
*       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
*       ..
*       .. Array Arguments ..
*       LOGICAL            SELECT( * )
*       DOUBLE PRECISION   P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
*      $                   VR( LDVR, * ), WORK( * )
*       ..
*
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DTGEVC computes some or all of the right and/or left eigenvectors of
*> a pair of real matrices (S,P), where S is a quasi-triangular matrix
*> and P is upper triangular.  Matrix pairs of this type are produced by
*> the generalized Schur factorization of a matrix pair (A,B):
*>
*>    A = Q*S*Z**T,  B = Q*P*Z**T
*>
*> as computed by DGGHRD + DHGEQZ.
*>
*> The right eigenvector x and the left eigenvector y of (S,P)
*> corresponding to an eigenvalue w are defined by:
*>
*>    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
*>
*> where y**H denotes the conjugate transpose of y.
*> The eigenvalues are not input to this routine, but are computed
*> directly from the diagonal blocks of S and P.
*>
*> This routine returns the matrices X and/or Y of right and left
*> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
*> where Z and Q are input matrices.
*> If Q and Z are the orthogonal factors from the generalized Schur
*> factorization of a matrix pair (A,B), then Z*X and Q*Y
*> are the matrices of right and left eigenvectors of (A,B).
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] SIDE
*> \verbatim
*>          SIDE is CHARACTER*1
*>          = 'R': compute right eigenvectors only;
*>          = 'L': compute left eigenvectors only;
*>          = 'B': compute both right and left eigenvectors.
*> \endverbatim
*>
*> \param[in] HOWMNY
*> \verbatim
*>          HOWMNY is CHARACTER*1
*>          = 'A': compute all right and/or left eigenvectors;
*>          = 'B': compute all right and/or left eigenvectors,
*>                 backtransformed by the matrices in VR and/or VL;
*>          = 'S': compute selected right and/or left eigenvectors,
*>                 specified by the logical array SELECT.
*> \endverbatim
*>
*> \param[in] SELECT
*> \verbatim
*>          SELECT is LOGICAL array, dimension (N)
*>          If HOWMNY='S', SELECT specifies the eigenvectors to be
*>          computed.  If w(j) is a real eigenvalue, the corresponding
*>          real eigenvector is computed if SELECT(j) is .TRUE..
*>          If w(j) and w(j+1) are the real and imaginary parts of a
*>          complex eigenvalue, the corresponding complex eigenvector
*>          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
*>          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
*>          set to .FALSE..
*>          Not referenced if HOWMNY = 'A' or 'B'.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices S and P.  N >= 0.
*> \endverbatim
*>
*> \param[in] S
*> \verbatim
*>          S is DOUBLE PRECISION array, dimension (LDS,N)
*>          The upper quasi-triangular matrix S from a generalized Schur
*>          factorization, as computed by DHGEQZ.
*> \endverbatim
*>
*> \param[in] LDS
*> \verbatim
*>          LDS is INTEGER
*>          The leading dimension of array S.  LDS >= max(1,N).
*> \endverbatim
*>
*> \param[in] P
*> \verbatim
*>          P is DOUBLE PRECISION array, dimension (LDP,N)
*>          The upper triangular matrix P from a generalized Schur
*>          factorization, as computed by DHGEQZ.
*>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
*>          of S must be in positive diagonal form.
*> \endverbatim
*>
*> \param[in] LDP
*> \verbatim
*>          LDP is INTEGER
*>          The leading dimension of array P.  LDP >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] VL
*> \verbatim
*>          VL is DOUBLE PRECISION array, dimension (LDVL,MM)
*>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
*>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
*>          of left Schur vectors returned by DHGEQZ).
*>          On exit, if SIDE = 'L' or 'B', VL contains:
*>          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
*>          if HOWMNY = 'B', the matrix Q*Y;
*>          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
*>                      SELECT, stored consecutively in the columns of
*>                      VL, in the same order as their eigenvalues.
*>
*>          A complex eigenvector corresponding to a complex eigenvalue
*>          is stored in two consecutive columns, the first holding the
*>          real part, and the second the imaginary part.
*>
*>          Not referenced if SIDE = 'R'.
*> \endverbatim
*>
*> \param[in] LDVL
*> \verbatim
*>          LDVL is INTEGER
*>          The leading dimension of array VL.  LDVL >= 1, and if
*>          SIDE = 'L' or 'B', LDVL >= N.
*> \endverbatim
*>
*> \param[in,out] VR
*> \verbatim
*>          VR is DOUBLE PRECISION array, dimension (LDVR,MM)
*>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
*>          contain an N-by-N matrix Z (usually the orthogonal matrix Z
*>          of right Schur vectors returned by DHGEQZ).
*>
*>          On exit, if SIDE = 'R' or 'B', VR contains:
*>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
*>          if HOWMNY = 'B' or 'b', the matrix Z*X;
*>          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
*>                      specified by SELECT, stored consecutively in the
*>                      columns of VR, in the same order as their
*>                      eigenvalues.
*>
*>          A complex eigenvector corresponding to a complex eigenvalue
*>          is stored in two consecutive columns, the first holding the
*>          real part and the second the imaginary part.
*>
*>          Not referenced if SIDE = 'L'.
*> \endverbatim
*>
*> \param[in] LDVR
*> \verbatim
*>          LDVR is INTEGER
*>          The leading dimension of the array VR.  LDVR >= 1, and if
*>          SIDE = 'R' or 'B', LDVR >= N.
*> \endverbatim
*>
*> \param[in] MM
*> \verbatim
*>          MM is INTEGER
*>          The number of columns in the arrays VL and/or VR. MM >= M.
*> \endverbatim
*>
*> \param[out] M
*> \verbatim
*>          M is INTEGER
*>          The number of columns in the arrays VL and/or VR actually
*>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
*>          is set to N.  Each selected real eigenvector occupies one
*>          column and each selected complex eigenvector occupies two
*>          columns.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (6*N)
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit.
*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
*>          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
*>                eigenvalue.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tgevc
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  Allocation of workspace:
*>  ---------- -- ---------
*>
*>     WORK( j ) = 1-norm of j-th column of A, above the diagonal
*>     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
*>     WORK( 2*N+1:3*N ) = real part of eigenvector
*>     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
*>     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
*>     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
*>
*>  Rowwise vs. columnwise solution methods:
*>  ------- --  ---------- -------- -------
*>
*>  Finding a generalized eigenvector consists basically of solving the
*>  singular triangular system
*>
*>   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)
*>
*>  Consider finding the i-th right eigenvector (assume all eigenvalues
*>  are real). The equation to be solved is:
*>       n                   i
*>  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
*>      k=j                 k=j
*>
*>  where  C = (A - w B)  (The components v(i+1:n) are 0.)
*>
*>  The "rowwise" method is:
*>
*>  (1)  v(i) := 1
*>  for j = i-1,. . .,1:
*>                          i
*>      (2) compute  s = - sum C(j,k) v(k)   and
*>                        k=j+1
*>
*>      (3) v(j) := s / C(j,j)
*>
*>  Step 2 is sometimes called the "dot product" step, since it is an
*>  inner product between the j-th row and the portion of the eigenvector
*>  that has been computed so far.
*>
*>  The "columnwise" method consists basically in doing the sums
*>  for all the rows in parallel.  As each v(j) is computed, the
*>  contribution of v(j) times the j-th column of C is added to the
*>  partial sums.  Since FORTRAN arrays are stored columnwise, this has
*>  the advantage that at each step, the elements of C that are accessed
*>  are adjacent to one another, whereas with the rowwise method, the
*>  elements accessed at a step are spaced LDS (and LDP) words apart.
*>
*>  When finding left eigenvectors, the matrix in question is the
*>  transpose of the one in storage, so the rowwise method then
*>  actually accesses columns of A and B at each step, and so is the
*>  preferred method.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE dtgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
     $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          HOWMNY, SIDE
      INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
*     ..
*     .. Array Arguments ..
      LOGICAL            SELECT( * )
      DOUBLE PRECISION   P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
     $                   vr( ldvr, * ), work( * )
*     ..
*
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, SAFETY
      parameter( zero = 0.0d+0, one = 1.0d+0,
     $                   safety = 1.0d+2 )
*     ..
*     .. Local Scalars ..
      LOGICAL            COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
     $                   ilbbad, ilcomp, ilcplx, lsa, lsb
      INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
     $                   j, ja, jc, je, jr, jw, na, nw
      DOUBLE PRECISION   ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
     $                   bcoefr, big, bignum, bnorm, bscale, cim2a,
     $                   cim2b, cimaga, cimagb, cre2a, cre2b, creala,
     $                   crealb, dmin, safmin, salfar, sbeta, scale,
     $                   small, temp, temp2, temp2i, temp2r, ulp, xmax,
     $                   xscale
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
     $                   sump( 2, 2 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           lsame, dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dgemv, dlacpy, dlag2, dlaln2, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min
*     ..
*     .. Executable Statements ..
*
*     Decode and Test the input parameters
*
      IF( lsame( howmny, 'A' ) ) THEN
         ihwmny = 1
         ilall = .true.
         ilback = .false.
      ELSE IF( lsame( howmny, 'S' ) ) THEN
         ihwmny = 2
         ilall = .false.
         ilback = .false.
      ELSE IF( lsame( howmny, 'B' ) ) THEN
         ihwmny = 3
         ilall = .true.
         ilback = .true.
      ELSE
         ihwmny = -1
         ilall = .true.
      END IF
*
      IF( lsame( side, 'R' ) ) THEN
         iside = 1
         compl = .false.
         compr = .true.
      ELSE IF( lsame( side, 'L' ) ) THEN
         iside = 2
         compl = .true.
         compr = .false.
      ELSE IF( lsame( side, 'B' ) ) THEN
         iside = 3
         compl = .true.
         compr = .true.
      ELSE
         iside = -1
      END IF
*
      info = 0
      IF( iside.LT.0 ) THEN
         info = -1
      ELSE IF( ihwmny.LT.0 ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( lds.LT.max( 1, n ) ) THEN
         info = -6
      ELSE IF( ldp.LT.max( 1, n ) ) THEN
         info = -8
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DTGEVC', -info )
         RETURN
      END IF
*
*     Count the number of eigenvectors to be computed
*
      IF( .NOT.ilall ) THEN
         im = 0
         ilcplx = .false.
         DO 10 j = 1, n
            IF( ilcplx ) THEN
               ilcplx = .false.
               GO TO 10
            END IF
            IF( j.LT.n ) THEN
               IF( s( j+1, j ).NE.zero )
     $            ilcplx = .true.
            END IF
            IF( ilcplx ) THEN
               IF( SELECT( j ) .OR. SELECT( j+1 ) )
     $            im = im + 2
            ELSE
               IF( SELECT( j ) )
     $            im = im + 1
            END IF
   10    CONTINUE
      ELSE
         im = n
      END IF
*
*     Check 2-by-2 diagonal blocks of A, B
*
      ilabad = .false.
      ilbbad = .false.
      DO 20 j = 1, n - 1
         IF( s( j+1, j ).NE.zero ) THEN
            IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
     $          p( j, j+1 ).NE.zero )ilbbad = .true.
            IF( j.LT.n-1 ) THEN
               IF( s( j+2, j+1 ).NE.zero )
     $            ilabad = .true.
            END IF
         END IF
   20 CONTINUE
*
      IF( ilabad ) THEN
         info = -5
      ELSE IF( ilbbad ) THEN
         info = -7
      ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
         info = -10
      ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
         info = -12
      ELSE IF( mm.LT.im ) THEN
         info = -13
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DTGEVC', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      m = im
      IF( n.EQ.0 )
     $   RETURN
*
*     Machine Constants
*
      safmin = dlamch( 'Safe minimum' )
      big = one / safmin
      ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
      small = safmin*n / ulp
      big = one / small
      bignum = one / ( safmin*n )
*
*     Compute the 1-norm of each column of the strictly upper triangular
*     part (i.e., excluding all elements belonging to the diagonal
*     blocks) of A and B to check for possible overflow in the
*     triangular solver.
*
      anorm = abs( s( 1, 1 ) )
      IF( n.GT.1 )
     $   anorm = anorm + abs( s( 2, 1 ) )
      bnorm = abs( p( 1, 1 ) )
      work( 1 ) = zero
      work( n+1 ) = zero
*
      DO 50 j = 2, n
         temp = zero
         temp2 = zero
         IF( s( j, j-1 ).EQ.zero ) THEN
            iend = j - 1
         ELSE
            iend = j - 2
         END IF
         DO 30 i = 1, iend
            temp = temp + abs( s( i, j ) )
            temp2 = temp2 + abs( p( i, j ) )
   30    CONTINUE
         work( j ) = temp
         work( n+j ) = temp2
         DO 40 i = iend + 1, min( j+1, n )
            temp = temp + abs( s( i, j ) )
            temp2 = temp2 + abs( p( i, j ) )
   40    CONTINUE
         anorm = max( anorm, temp )
         bnorm = max( bnorm, temp2 )
   50 CONTINUE
*
      ascale = one / max( anorm, safmin )
      bscale = one / max( bnorm, safmin )
*
*     Left eigenvectors
*
      IF( compl ) THEN
         ieig = 0
*
*        Main loop over eigenvalues
*
         ilcplx = .false.
         DO 220 je = 1, n
*
*           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
*           (b) this would be the second of a complex pair.
*           Check for complex eigenvalue, so as to be sure of which
*           entry(-ies) of SELECT to look at.
*
            IF( ilcplx ) THEN
               ilcplx = .false.
               GO TO 220
            END IF
            nw = 1
            IF( je.LT.n ) THEN
               IF( s( je+1, je ).NE.zero ) THEN
                  ilcplx = .true.
                  nw = 2
               END IF
            END IF
            IF( ilall ) THEN
               ilcomp = .true.
            ELSE IF( ilcplx ) THEN
               ilcomp = SELECT( je ) .OR. SELECT( je+1 )
            ELSE
               ilcomp = SELECT( je )
            END IF
            IF( .NOT.ilcomp )
     $         GO TO 220
*
*           Decide if (a) singular pencil, (b) real eigenvalue, or
*           (c) complex eigenvalue.
*
            IF( .NOT.ilcplx ) THEN
               IF( abs( s( je, je ) ).LE.safmin .AND.
     $             abs( p( je, je ) ).LE.safmin ) THEN
*
*                 Singular matrix pencil -- return unit eigenvector
*
                  ieig = ieig + 1
                  DO 60 jr = 1, n
                     vl( jr, ieig ) = zero
   60             CONTINUE
                  vl( ieig, ieig ) = one
                  GO TO 220
               END IF
            END IF
*
*           Clear vector
*
            DO 70 jr = 1, nw*n
               work( 2*n+jr ) = zero
   70       CONTINUE
*                                                 T
*           Compute coefficients in  ( a A - b B )  y = 0
*              a  is  ACOEF
*              b  is  BCOEFR + i*BCOEFI
*
            IF( .NOT.ilcplx ) THEN
*
*              Real eigenvalue
*
               temp = one / max( abs( s( je, je ) )*ascale,
     $                abs( p( je, je ) )*bscale, safmin )
               salfar = ( temp*s( je, je ) )*ascale
               sbeta = ( temp*p( je, je ) )*bscale
               acoef = sbeta*ascale
               bcoefr = salfar*bscale
               bcoefi = zero
*
*              Scale to avoid underflow
*
               scale = one
               lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
               lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
     $               small
               IF( lsa )
     $            scale = ( small / abs( sbeta ) )*min( anorm, big )
               IF( lsb )
     $            scale = max( scale, ( small / abs( salfar ) )*
     $                    min( bnorm, big ) )
               IF( lsa .OR. lsb ) THEN
                  scale = min( scale, one /
     $                    ( safmin*max( one, abs( acoef ),
     $                    abs( bcoefr ) ) ) )
                  IF( lsa ) THEN
                     acoef = ascale*( scale*sbeta )
                  ELSE
                     acoef = scale*acoef
                  END IF
                  IF( lsb ) THEN
                     bcoefr = bscale*( scale*salfar )
                  ELSE
                     bcoefr = scale*bcoefr
                  END IF
               END IF
               acoefa = abs( acoef )
               bcoefa = abs( bcoefr )
*
*              First component is 1
*
               work( 2*n+je ) = one
               xmax = one
            ELSE
*
*              Complex eigenvalue
*
               CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
     $                     safmin*safety, acoef, temp, bcoefr, temp2,
     $                     bcoefi )
               bcoefi = -bcoefi
               IF( bcoefi.EQ.zero ) THEN
                  info = je
                  RETURN
               END IF
*
*              Scale to avoid over/underflow
*
               acoefa = abs( acoef )
               bcoefa = abs( bcoefr ) + abs( bcoefi )
               scale = one
               IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
     $            scale = ( safmin / ulp ) / acoefa
               IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
     $            scale = max( scale, ( safmin / ulp ) / bcoefa )
               IF( safmin*acoefa.GT.ascale )
     $            scale = ascale / ( safmin*acoefa )
               IF( safmin*bcoefa.GT.bscale )
     $            scale = min( scale, bscale / ( safmin*bcoefa ) )
               IF( scale.NE.one ) THEN
                  acoef = scale*acoef
                  acoefa = abs( acoef )
                  bcoefr = scale*bcoefr
                  bcoefi = scale*bcoefi
                  bcoefa = abs( bcoefr ) + abs( bcoefi )
               END IF
*
*              Compute first two components of eigenvector
*
               temp = acoef*s( je+1, je )
               temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
               temp2i = -bcoefi*p( je, je )
               IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
                  work( 2*n+je ) = one
                  work( 3*n+je ) = zero
                  work( 2*n+je+1 ) = -temp2r / temp
                  work( 3*n+je+1 ) = -temp2i / temp
               ELSE
                  work( 2*n+je+1 ) = one
                  work( 3*n+je+1 ) = zero
                  temp = acoef*s( je, je+1 )
                  work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
     $                             s( je+1, je+1 ) ) / temp
                  work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
               END IF
               xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
     $                abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
            END IF
*
            dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
*
*                                           T
*           Triangular solve of  (a A - b B)  y = 0
*
*                                   T
*           (rowwise in  (a A - b B) , or columnwise in (a A - b B) )
*
            il2by2 = .false.
*
            DO 160 j = je + nw, n
               IF( il2by2 ) THEN
                  il2by2 = .false.
                  GO TO 160
               END IF
*
               na = 1
               bdiag( 1 ) = p( j, j )
               IF( j.LT.n ) THEN
                  IF( s( j+1, j ).NE.zero ) THEN
                     il2by2 = .true.
                     bdiag( 2 ) = p( j+1, j+1 )
                     na = 2
                  END IF
               END IF
*
*              Check whether scaling is necessary for dot products
*
               xscale = one / max( one, xmax )
               temp = max( work( j ), work( n+j ),
     $                acoefa*work( j )+bcoefa*work( n+j ) )
               IF( il2by2 )
     $            temp = max( temp, work( j+1 ), work( n+j+1 ),
     $                   acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
               IF( temp.GT.bignum*xscale ) THEN
                  DO 90 jw = 0, nw - 1
                     DO 80 jr = je, j - 1
                        work( ( jw+2 )*n+jr ) = xscale*
     $                     work( ( jw+2 )*n+jr )
   80                CONTINUE
   90             CONTINUE
                  xmax = xmax*xscale
               END IF
*
*              Compute dot products
*
*                    j-1
*              SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
*                    k=je
*
*              To reduce the op count, this is done as
*
*              _        j-1                  _        j-1
*              a*conjg( sum  S(k,j)*x(k) ) - b*conjg( sum  P(k,j)*x(k) )
*                       k=je                          k=je
*
*              which may cause underflow problems if A or B are close
*              to underflow.  (E.g., less than SMALL.)
*
*
               DO 120 jw = 1, nw
                  DO 110 ja = 1, na
                     sums( ja, jw ) = zero
                     sump( ja, jw ) = zero
*
                     DO 100 jr = je, j - 1
                        sums( ja, jw ) = sums( ja, jw ) +
     $                                   s( jr, j+ja-1 )*
     $                                   work( ( jw+1 )*n+jr )
                        sump( ja, jw ) = sump( ja, jw ) +
     $                                   p( jr, j+ja-1 )*
     $                                   work( ( jw+1 )*n+jr )
  100                CONTINUE
  110             CONTINUE
  120          CONTINUE
*
               DO 130 ja = 1, na
                  IF( ilcplx ) THEN
                     sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
     $                              bcoefr*sump( ja, 1 ) -
     $                              bcoefi*sump( ja, 2 )
                     sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
     $                              bcoefr*sump( ja, 2 ) +
     $                              bcoefi*sump( ja, 1 )
                  ELSE
                     sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
     $                              bcoefr*sump( ja, 1 )
                  END IF
  130          CONTINUE
*
*                                  T
*              Solve  ( a A - b B )  y = SUM(,)
*              with scaling and perturbation of the denominator
*
               CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
     $                      bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
     $                      bcoefi, work( 2*n+j ), n, scale, temp,
     $                      iinfo )
               IF( scale.LT.one ) THEN
                  DO 150 jw = 0, nw - 1
                     DO 140 jr = je, j - 1
                        work( ( jw+2 )*n+jr ) = scale*
     $                     work( ( jw+2 )*n+jr )
  140                CONTINUE
  150             CONTINUE
                  xmax = scale*xmax
               END IF
               xmax = max( xmax, temp )
  160       CONTINUE
*
*           Copy eigenvector to VL, back transforming if
*           HOWMNY='B'.
*
            ieig = ieig + 1
            IF( ilback ) THEN
               DO 170 jw = 0, nw - 1
                  CALL dgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
     $                        work( ( jw+2 )*n+je ), 1, zero,
     $                        work( ( jw+4 )*n+1 ), 1 )
  170          CONTINUE
               CALL dlacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
     $                      ldvl )
               ibeg = 1
            ELSE
               CALL dlacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
     $                      ldvl )
               ibeg = je
            END IF
*
*           Scale eigenvector
*
            xmax = zero
            IF( ilcplx ) THEN
               DO 180 j = ibeg, n
                  xmax = max( xmax, abs( vl( j, ieig ) )+
     $                   abs( vl( j, ieig+1 ) ) )
  180          CONTINUE
            ELSE
               DO 190 j = ibeg, n
                  xmax = max( xmax, abs( vl( j, ieig ) ) )
  190          CONTINUE
            END IF
*
            IF( xmax.GT.safmin ) THEN
               xscale = one / xmax
*
               DO 210 jw = 0, nw - 1
                  DO 200 jr = ibeg, n
                     vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
  200             CONTINUE
  210          CONTINUE
            END IF
            ieig = ieig + nw - 1
*
  220    CONTINUE
      END IF
*
*     Right eigenvectors
*
      IF( compr ) THEN
         ieig = im + 1
*
*        Main loop over eigenvalues
*
         ilcplx = .false.
         DO 500 je = n, 1, -1
*
*           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
*           (b) this would be the second of a complex pair.
*           Check for complex eigenvalue, so as to be sure of which
*           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
*           or SELECT(JE-1).
*           If this is a complex pair, the 2-by-2 diagonal block
*           corresponding to the eigenvalue is in rows/columns JE-1:JE
*
            IF( ilcplx ) THEN
               ilcplx = .false.
               GO TO 500
            END IF
            nw = 1
            IF( je.GT.1 ) THEN
               IF( s( je, je-1 ).NE.zero ) THEN
                  ilcplx = .true.
                  nw = 2
               END IF
            END IF
            IF( ilall ) THEN
               ilcomp = .true.
            ELSE IF( ilcplx ) THEN
               ilcomp = SELECT( je ) .OR. SELECT( je-1 )
            ELSE
               ilcomp = SELECT( je )
            END IF
            IF( .NOT.ilcomp )
     $         GO TO 500
*
*           Decide if (a) singular pencil, (b) real eigenvalue, or
*           (c) complex eigenvalue.
*
            IF( .NOT.ilcplx ) THEN
               IF( abs( s( je, je ) ).LE.safmin .AND.
     $             abs( p( je, je ) ).LE.safmin ) THEN
*
*                 Singular matrix pencil -- unit eigenvector
*
                  ieig = ieig - 1
                  DO 230 jr = 1, n
                     vr( jr, ieig ) = zero
  230             CONTINUE
                  vr( ieig, ieig ) = one
                  GO TO 500
               END IF
            END IF
*
*           Clear vector
*
            DO 250 jw = 0, nw - 1
               DO 240 jr = 1, n
                  work( ( jw+2 )*n+jr ) = zero
  240          CONTINUE
  250       CONTINUE
*
*           Compute coefficients in  ( a A - b B ) x = 0
*              a  is  ACOEF
*              b  is  BCOEFR + i*BCOEFI
*
            IF( .NOT.ilcplx ) THEN
*
*              Real eigenvalue
*
               temp = one / max( abs( s( je, je ) )*ascale,
     $                abs( p( je, je ) )*bscale, safmin )
               salfar = ( temp*s( je, je ) )*ascale
               sbeta = ( temp*p( je, je ) )*bscale
               acoef = sbeta*ascale
               bcoefr = salfar*bscale
               bcoefi = zero
*
*              Scale to avoid underflow
*
               scale = one
               lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
               lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
     $               small
               IF( lsa )
     $            scale = ( small / abs( sbeta ) )*min( anorm, big )
               IF( lsb )
     $            scale = max( scale, ( small / abs( salfar ) )*
     $                    min( bnorm, big ) )
               IF( lsa .OR. lsb ) THEN
                  scale = min( scale, one /
     $                    ( safmin*max( one, abs( acoef ),
     $                    abs( bcoefr ) ) ) )
                  IF( lsa ) THEN
                     acoef = ascale*( scale*sbeta )
                  ELSE
                     acoef = scale*acoef
                  END IF
                  IF( lsb ) THEN
                     bcoefr = bscale*( scale*salfar )
                  ELSE
                     bcoefr = scale*bcoefr
                  END IF
               END IF
               acoefa = abs( acoef )
               bcoefa = abs( bcoefr )
*
*              First component is 1
*
               work( 2*n+je ) = one
               xmax = one
*
*              Compute contribution from column JE of A and B to sum
*              (See "Further Details", above.)
*
               DO 260 jr = 1, je - 1
                  work( 2*n+jr ) = bcoefr*p( jr, je ) -
     $                             acoef*s( jr, je )
  260          CONTINUE
            ELSE
*
*              Complex eigenvalue
*
               CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
     $                     safmin*safety, acoef, temp, bcoefr, temp2,
     $                     bcoefi )
               IF( bcoefi.EQ.zero ) THEN
                  info = je - 1
                  RETURN
               END IF
*
*              Scale to avoid over/underflow
*
               acoefa = abs( acoef )
               bcoefa = abs( bcoefr ) + abs( bcoefi )
               scale = one
               IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
     $            scale = ( safmin / ulp ) / acoefa
               IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
     $            scale = max( scale, ( safmin / ulp ) / bcoefa )
               IF( safmin*acoefa.GT.ascale )
     $            scale = ascale / ( safmin*acoefa )
               IF( safmin*bcoefa.GT.bscale )
     $            scale = min( scale, bscale / ( safmin*bcoefa ) )
               IF( scale.NE.one ) THEN
                  acoef = scale*acoef
                  acoefa = abs( acoef )
                  bcoefr = scale*bcoefr
                  bcoefi = scale*bcoefi
                  bcoefa = abs( bcoefr ) + abs( bcoefi )
               END IF
*
*              Compute first two components of eigenvector
*              and contribution to sums
*
               temp = acoef*s( je, je-1 )
               temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
               temp2i = -bcoefi*p( je, je )
               IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
                  work( 2*n+je ) = one
                  work( 3*n+je ) = zero
                  work( 2*n+je-1 ) = -temp2r / temp
                  work( 3*n+je-1 ) = -temp2i / temp
               ELSE
                  work( 2*n+je-1 ) = one
                  work( 3*n+je-1 ) = zero
                  temp = acoef*s( je-1, je )
                  work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
     $                             s( je-1, je-1 ) ) / temp
                  work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
               END IF
*
               xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
     $                abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
*
*              Compute contribution from columns JE and JE-1
*              of A and B to the sums.
*
               creala = acoef*work( 2*n+je-1 )
               cimaga = acoef*work( 3*n+je-1 )
               crealb = bcoefr*work( 2*n+je-1 ) -
     $                  bcoefi*work( 3*n+je-1 )
               cimagb = bcoefi*work( 2*n+je-1 ) +
     $                  bcoefr*work( 3*n+je-1 )
               cre2a = acoef*work( 2*n+je )
               cim2a = acoef*work( 3*n+je )
               cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
               cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
               DO 270 jr = 1, je - 2
                  work( 2*n+jr ) = -creala*s( jr, je-1 ) +
     $                             crealb*p( jr, je-1 ) -
     $                             cre2a*s( jr, je ) + cre2b*p( jr, je )
                  work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
     $                             cimagb*p( jr, je-1 ) -
     $                             cim2a*s( jr, je ) + cim2b*p( jr, je )
  270          CONTINUE
            END IF
*
            dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
*
*           Columnwise triangular solve of  (a A - b B)  x = 0
*
            il2by2 = .false.
            DO 370 j = je - nw, 1, -1
*
*              If a 2-by-2 block, is in position j-1:j, wait until
*              next iteration to process it (when it will be j:j+1)
*
               IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
                  IF( s( j, j-1 ).NE.zero ) THEN
                     il2by2 = .true.
                     GO TO 370
                  END IF
               END IF
               bdiag( 1 ) = p( j, j )
               IF( il2by2 ) THEN
                  na = 2
                  bdiag( 2 ) = p( j+1, j+1 )
               ELSE
                  na = 1
               END IF
*
*              Compute x(j) (and x(j+1), if 2-by-2 block)
*
               CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
     $                      lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
     $                      n, bcoefr, bcoefi, sum, 2, scale, temp,
     $                      iinfo )
               IF( scale.LT.one ) THEN
*
                  DO 290 jw = 0, nw - 1
                     DO 280 jr = 1, je
                        work( ( jw+2 )*n+jr ) = scale*
     $                     work( ( jw+2 )*n+jr )
  280                CONTINUE
  290             CONTINUE
               END IF
               xmax = max( scale*xmax, temp )
*
               DO 310 jw = 1, nw
                  DO 300 ja = 1, na
                     work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
  300             CONTINUE
  310          CONTINUE
*
*              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
*
               IF( j.GT.1 ) THEN
*
*                 Check whether scaling is necessary for sum.
*
                  xscale = one / max( one, xmax )
                  temp = acoefa*work( j ) + bcoefa*work( n+j )
                  IF( il2by2 )
     $               temp = max( temp, acoefa*work( j+1 )+bcoefa*
     $                      work( n+j+1 ) )
                  temp = max( temp, acoefa, bcoefa )
                  IF( temp.GT.bignum*xscale ) THEN
*
                     DO 330 jw = 0, nw - 1
                        DO 320 jr = 1, je
                           work( ( jw+2 )*n+jr ) = xscale*
     $                        work( ( jw+2 )*n+jr )
  320                   CONTINUE
  330                CONTINUE
                     xmax = xmax*xscale
                  END IF
*
*                 Compute the contributions of the off-diagonals of
*                 column j (and j+1, if 2-by-2 block) of A and B to the
*                 sums.
*
*
                  DO 360 ja = 1, na
                     IF( ilcplx ) THEN
                        creala = acoef*work( 2*n+j+ja-1 )
                        cimaga = acoef*work( 3*n+j+ja-1 )
                        crealb = bcoefr*work( 2*n+j+ja-1 ) -
     $                           bcoefi*work( 3*n+j+ja-1 )
                        cimagb = bcoefi*work( 2*n+j+ja-1 ) +
     $                           bcoefr*work( 3*n+j+ja-1 )
                        DO 340 jr = 1, j - 1
                           work( 2*n+jr ) = work( 2*n+jr ) -
     $                                      creala*s( jr, j+ja-1 ) +
     $                                      crealb*p( jr, j+ja-1 )
                           work( 3*n+jr ) = work( 3*n+jr ) -
     $                                      cimaga*s( jr, j+ja-1 ) +
     $                                      cimagb*p( jr, j+ja-1 )
  340                   CONTINUE
                     ELSE
                        creala = acoef*work( 2*n+j+ja-1 )
                        crealb = bcoefr*work( 2*n+j+ja-1 )
                        DO 350 jr = 1, j - 1
                           work( 2*n+jr ) = work( 2*n+jr ) -
     $                                      creala*s( jr, j+ja-1 ) +
     $                                      crealb*p( jr, j+ja-1 )
  350                   CONTINUE
                     END IF
  360             CONTINUE
               END IF
*
               il2by2 = .false.
  370       CONTINUE
*
*           Copy eigenvector to VR, back transforming if
*           HOWMNY='B'.
*
            ieig = ieig - nw
            IF( ilback ) THEN
*
               DO 410 jw = 0, nw - 1
                  DO 380 jr = 1, n
                     work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
     $                                       vr( jr, 1 )
  380             CONTINUE
*
*                 A series of compiler directives to defeat
*                 vectorization for the next loop
*
*
                  DO 400 jc = 2, je
                     DO 390 jr = 1, n
                        work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
     $                     work( ( jw+2 )*n+jc )*vr( jr, jc )
  390                CONTINUE
  400             CONTINUE
  410          CONTINUE
*
               DO 430 jw = 0, nw - 1
                  DO 420 jr = 1, n
                     vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
  420             CONTINUE
  430          CONTINUE
*
               iend = n
            ELSE
               DO 450 jw = 0, nw - 1
                  DO 440 jr = 1, n
                     vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
  440             CONTINUE
  450          CONTINUE
*
               iend = je
            END IF
*
*           Scale eigenvector
*
            xmax = zero
            IF( ilcplx ) THEN
               DO 460 j = 1, iend
                  xmax = max( xmax, abs( vr( j, ieig ) )+
     $                   abs( vr( j, ieig+1 ) ) )
  460          CONTINUE
            ELSE
               DO 470 j = 1, iend
                  xmax = max( xmax, abs( vr( j, ieig ) ) )
  470          CONTINUE
            END IF
*
            IF( xmax.GT.safmin ) THEN
               xscale = one / xmax
               DO 490 jw = 0, nw - 1
                  DO 480 jr = 1, iend
                     vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
  480             CONTINUE
  490          CONTINUE
            END IF
  500    CONTINUE
      END IF
*
      RETURN
*
*     End of DTGEVC
*
      END