0001:       SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
0002:      $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
0003: *
0004: *  -- LAPACK routine (version 3.2) --
0005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
0006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
0007: *     November 2006
0008: *
0009: *     .. Scalar Arguments ..
0010:       CHARACTER          HOWMNY, SIDE
0011:       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
0012: *     ..
0013: *     .. Array Arguments ..
0014:       LOGICAL            SELECT( * )
0015:       DOUBLE PRECISION   P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
0016:      $                   VR( LDVR, * ), WORK( * )
0017: *     ..
0018: *
0019: *
0020: *  Purpose
0021: *  =======
0022: *
0023: *  DTGEVC computes some or all of the right and/or left eigenvectors of
0024: *  a pair of real matrices (S,P), where S is a quasi-triangular matrix
0025: *  and P is upper triangular.  Matrix pairs of this type are produced by
0026: *  the generalized Schur factorization of a matrix pair (A,B):
0027: *
0028: *     A = Q*S*Z**T,  B = Q*P*Z**T
0029: *
0030: *  as computed by DGGHRD + DHGEQZ.
0031: *
0032: *  The right eigenvector x and the left eigenvector y of (S,P)
0033: *  corresponding to an eigenvalue w are defined by:
0034: *  
0035: *     S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
0036: *  
0037: *  where y**H denotes the conjugate tranpose of y.
0038: *  The eigenvalues are not input to this routine, but are computed
0039: *  directly from the diagonal blocks of S and P.
0040: *  
0041: *  This routine returns the matrices X and/or Y of right and left
0042: *  eigenvectors of (S,P), or the products Z*X and/or Q*Y,
0043: *  where Z and Q are input matrices.
0044: *  If Q and Z are the orthogonal factors from the generalized Schur
0045: *  factorization of a matrix pair (A,B), then Z*X and Q*Y
0046: *  are the matrices of right and left eigenvectors of (A,B).
0047: * 
0048: *  Arguments
0049: *  =========
0050: *
0051: *  SIDE    (input) CHARACTER*1
0052: *          = 'R': compute right eigenvectors only;
0053: *          = 'L': compute left eigenvectors only;
0054: *          = 'B': compute both right and left eigenvectors.
0055: *
0056: *  HOWMNY  (input) CHARACTER*1
0057: *          = 'A': compute all right and/or left eigenvectors;
0058: *          = 'B': compute all right and/or left eigenvectors,
0059: *                 backtransformed by the matrices in VR and/or VL;
0060: *          = 'S': compute selected right and/or left eigenvectors,
0061: *                 specified by the logical array SELECT.
0062: *
0063: *  SELECT  (input) LOGICAL array, dimension (N)
0064: *          If HOWMNY='S', SELECT specifies the eigenvectors to be
0065: *          computed.  If w(j) is a real eigenvalue, the corresponding
0066: *          real eigenvector is computed if SELECT(j) is .TRUE..
0067: *          If w(j) and w(j+1) are the real and imaginary parts of a
0068: *          complex eigenvalue, the corresponding complex eigenvector
0069: *          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
0070: *          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
0071: *          set to .FALSE..
0072: *          Not referenced if HOWMNY = 'A' or 'B'.
0073: *
0074: *  N       (input) INTEGER
0075: *          The order of the matrices S and P.  N >= 0.
0076: *
0077: *  S       (input) DOUBLE PRECISION array, dimension (LDS,N)
0078: *          The upper quasi-triangular matrix S from a generalized Schur
0079: *          factorization, as computed by DHGEQZ.
0080: *
0081: *  LDS     (input) INTEGER
0082: *          The leading dimension of array S.  LDS >= max(1,N).
0083: *
0084: *  P       (input) DOUBLE PRECISION array, dimension (LDP,N)
0085: *          The upper triangular matrix P from a generalized Schur
0086: *          factorization, as computed by DHGEQZ.
0087: *          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
0088: *          of S must be in positive diagonal form.
0089: *
0090: *  LDP     (input) INTEGER
0091: *          The leading dimension of array P.  LDP >= max(1,N).
0092: *
0093: *  VL      (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
0094: *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
0095: *          contain an N-by-N matrix Q (usually the orthogonal matrix Q
0096: *          of left Schur vectors returned by DHGEQZ).
0097: *          On exit, if SIDE = 'L' or 'B', VL contains:
0098: *          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
0099: *          if HOWMNY = 'B', the matrix Q*Y;
0100: *          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
0101: *                      SELECT, stored consecutively in the columns of
0102: *                      VL, in the same order as their eigenvalues.
0103: *
0104: *          A complex eigenvector corresponding to a complex eigenvalue
0105: *          is stored in two consecutive columns, the first holding the
0106: *          real part, and the second the imaginary part.
0107: *
0108: *          Not referenced if SIDE = 'R'.
0109: *
0110: *  LDVL    (input) INTEGER
0111: *          The leading dimension of array VL.  LDVL >= 1, and if
0112: *          SIDE = 'L' or 'B', LDVL >= N.
0113: *
0114: *  VR      (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
0115: *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
0116: *          contain an N-by-N matrix Z (usually the orthogonal matrix Z
0117: *          of right Schur vectors returned by DHGEQZ).
0118: *
0119: *          On exit, if SIDE = 'R' or 'B', VR contains:
0120: *          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
0121: *          if HOWMNY = 'B' or 'b', the matrix Z*X;
0122: *          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
0123: *                      specified by SELECT, stored consecutively in the
0124: *                      columns of VR, in the same order as their
0125: *                      eigenvalues.
0126: *
0127: *          A complex eigenvector corresponding to a complex eigenvalue
0128: *          is stored in two consecutive columns, the first holding the
0129: *          real part and the second the imaginary part.
0130: *          
0131: *          Not referenced if SIDE = 'L'.
0132: *
0133: *  LDVR    (input) INTEGER
0134: *          The leading dimension of the array VR.  LDVR >= 1, and if
0135: *          SIDE = 'R' or 'B', LDVR >= N.
0136: *
0137: *  MM      (input) INTEGER
0138: *          The number of columns in the arrays VL and/or VR. MM >= M.
0139: *
0140: *  M       (output) INTEGER
0141: *          The number of columns in the arrays VL and/or VR actually
0142: *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
0143: *          is set to N.  Each selected real eigenvector occupies one
0144: *          column and each selected complex eigenvector occupies two
0145: *          columns.
0146: *
0147: *  WORK    (workspace) DOUBLE PRECISION array, dimension (6*N)
0148: *
0149: *  INFO    (output) INTEGER
0150: *          = 0:  successful exit.
0151: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
0152: *          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
0153: *                eigenvalue.
0154: *
0155: *  Further Details
0156: *  ===============
0157: *
0158: *  Allocation of workspace:
0159: *  ---------- -- ---------
0160: *
0161: *     WORK( j ) = 1-norm of j-th column of A, above the diagonal
0162: *     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
0163: *     WORK( 2*N+1:3*N ) = real part of eigenvector
0164: *     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
0165: *     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
0166: *     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
0167: *
0168: *  Rowwise vs. columnwise solution methods:
0169: *  ------- --  ---------- -------- -------
0170: *
0171: *  Finding a generalized eigenvector consists basically of solving the
0172: *  singular triangular system
0173: *
0174: *   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)
0175: *
0176: *  Consider finding the i-th right eigenvector (assume all eigenvalues
0177: *  are real). The equation to be solved is:
0178: *       n                   i
0179: *  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
0180: *      k=j                 k=j
0181: *
0182: *  where  C = (A - w B)  (The components v(i+1:n) are 0.)
0183: *
0184: *  The "rowwise" method is:
0185: *
0186: *  (1)  v(i) := 1
0187: *  for j = i-1,. . .,1:
0188: *                          i
0189: *      (2) compute  s = - sum C(j,k) v(k)   and
0190: *                        k=j+1
0191: *
0192: *      (3) v(j) := s / C(j,j)
0193: *
0194: *  Step 2 is sometimes called the "dot product" step, since it is an
0195: *  inner product between the j-th row and the portion of the eigenvector
0196: *  that has been computed so far.
0197: *
0198: *  The "columnwise" method consists basically in doing the sums
0199: *  for all the rows in parallel.  As each v(j) is computed, the
0200: *  contribution of v(j) times the j-th column of C is added to the
0201: *  partial sums.  Since FORTRAN arrays are stored columnwise, this has
0202: *  the advantage that at each step, the elements of C that are accessed
0203: *  are adjacent to one another, whereas with the rowwise method, the
0204: *  elements accessed at a step are spaced LDS (and LDP) words apart.
0205: *
0206: *  When finding left eigenvectors, the matrix in question is the
0207: *  transpose of the one in storage, so the rowwise method then
0208: *  actually accesses columns of A and B at each step, and so is the
0209: *  preferred method.
0210: *
0211: *  =====================================================================
0212: *
0213: *     .. Parameters ..
0214:       DOUBLE PRECISION   ZERO, ONE, SAFETY
0215:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
0216:      $                   SAFETY = 1.0D+2 )
0217: *     ..
0218: *     .. Local Scalars ..
0219:       LOGICAL            COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
0220:      $                   ILBBAD, ILCOMP, ILCPLX, LSA, LSB
0221:       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
0222:      $                   J, JA, JC, JE, JR, JW, NA, NW
0223:       DOUBLE PRECISION   ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
0224:      $                   BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
0225:      $                   CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
0226:      $                   CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
0227:      $                   SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
0228:      $                   XSCALE
0229: *     ..
0230: *     .. Local Arrays ..
0231:       DOUBLE PRECISION   BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
0232:      $                   SUMP( 2, 2 )
0233: *     ..
0234: *     .. External Functions ..
0235:       LOGICAL            LSAME
0236:       DOUBLE PRECISION   DLAMCH
0237:       EXTERNAL           LSAME, DLAMCH
0238: *     ..
0239: *     .. External Subroutines ..
0240:       EXTERNAL           DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA
0241: *     ..
0242: *     .. Intrinsic Functions ..
0243:       INTRINSIC          ABS, MAX, MIN
0244: *     ..
0245: *     .. Executable Statements ..
0246: *
0247: *     Decode and Test the input parameters
0248: *
0249:       IF( LSAME( HOWMNY, 'A' ) ) THEN
0250:          IHWMNY = 1
0251:          ILALL = .TRUE.
0252:          ILBACK = .FALSE.
0253:       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
0254:          IHWMNY = 2
0255:          ILALL = .FALSE.
0256:          ILBACK = .FALSE.
0257:       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
0258:          IHWMNY = 3
0259:          ILALL = .TRUE.
0260:          ILBACK = .TRUE.
0261:       ELSE
0262:          IHWMNY = -1
0263:          ILALL = .TRUE.
0264:       END IF
0265: *
0266:       IF( LSAME( SIDE, 'R' ) ) THEN
0267:          ISIDE = 1
0268:          COMPL = .FALSE.
0269:          COMPR = .TRUE.
0270:       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
0271:          ISIDE = 2
0272:          COMPL = .TRUE.
0273:          COMPR = .FALSE.
0274:       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
0275:          ISIDE = 3
0276:          COMPL = .TRUE.
0277:          COMPR = .TRUE.
0278:       ELSE
0279:          ISIDE = -1
0280:       END IF
0281: *
0282:       INFO = 0
0283:       IF( ISIDE.LT.0 ) THEN
0284:          INFO = -1
0285:       ELSE IF( IHWMNY.LT.0 ) THEN
0286:          INFO = -2
0287:       ELSE IF( N.LT.0 ) THEN
0288:          INFO = -4
0289:       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
0290:          INFO = -6
0291:       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
0292:          INFO = -8
0293:       END IF
0294:       IF( INFO.NE.0 ) THEN
0295:          CALL XERBLA( 'DTGEVC', -INFO )
0296:          RETURN
0297:       END IF
0298: *
0299: *     Count the number of eigenvectors to be computed
0300: *
0301:       IF( .NOT.ILALL ) THEN
0302:          IM = 0
0303:          ILCPLX = .FALSE.
0304:          DO 10 J = 1, N
0305:             IF( ILCPLX ) THEN
0306:                ILCPLX = .FALSE.
0307:                GO TO 10
0308:             END IF
0309:             IF( J.LT.N ) THEN
0310:                IF( S( J+1, J ).NE.ZERO )
0311:      $            ILCPLX = .TRUE.
0312:             END IF
0313:             IF( ILCPLX ) THEN
0314:                IF( SELECT( J ) .OR. SELECT( J+1 ) )
0315:      $            IM = IM + 2
0316:             ELSE
0317:                IF( SELECT( J ) )
0318:      $            IM = IM + 1
0319:             END IF
0320:    10    CONTINUE
0321:       ELSE
0322:          IM = N
0323:       END IF
0324: *
0325: *     Check 2-by-2 diagonal blocks of A, B
0326: *
0327:       ILABAD = .FALSE.
0328:       ILBBAD = .FALSE.
0329:       DO 20 J = 1, N - 1
0330:          IF( S( J+1, J ).NE.ZERO ) THEN
0331:             IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
0332:      $          P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
0333:             IF( J.LT.N-1 ) THEN
0334:                IF( S( J+2, J+1 ).NE.ZERO )
0335:      $            ILABAD = .TRUE.
0336:             END IF
0337:          END IF
0338:    20 CONTINUE
0339: *
0340:       IF( ILABAD ) THEN
0341:          INFO = -5
0342:       ELSE IF( ILBBAD ) THEN
0343:          INFO = -7
0344:       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
0345:          INFO = -10
0346:       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
0347:          INFO = -12
0348:       ELSE IF( MM.LT.IM ) THEN
0349:          INFO = -13
0350:       END IF
0351:       IF( INFO.NE.0 ) THEN
0352:          CALL XERBLA( 'DTGEVC', -INFO )
0353:          RETURN
0354:       END IF
0355: *
0356: *     Quick return if possible
0357: *
0358:       M = IM
0359:       IF( N.EQ.0 )
0360:      $   RETURN
0361: *
0362: *     Machine Constants
0363: *
0364:       SAFMIN = DLAMCH( 'Safe minimum' )
0365:       BIG = ONE / SAFMIN
0366:       CALL DLABAD( SAFMIN, BIG )
0367:       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
0368:       SMALL = SAFMIN*N / ULP
0369:       BIG = ONE / SMALL
0370:       BIGNUM = ONE / ( SAFMIN*N )
0371: *
0372: *     Compute the 1-norm of each column of the strictly upper triangular
0373: *     part (i.e., excluding all elements belonging to the diagonal
0374: *     blocks) of A and B to check for possible overflow in the
0375: *     triangular solver.
0376: *
0377:       ANORM = ABS( S( 1, 1 ) )
0378:       IF( N.GT.1 )
0379:      $   ANORM = ANORM + ABS( S( 2, 1 ) )
0380:       BNORM = ABS( P( 1, 1 ) )
0381:       WORK( 1 ) = ZERO
0382:       WORK( N+1 ) = ZERO
0383: *
0384:       DO 50 J = 2, N
0385:          TEMP = ZERO
0386:          TEMP2 = ZERO
0387:          IF( S( J, J-1 ).EQ.ZERO ) THEN
0388:             IEND = J - 1
0389:          ELSE
0390:             IEND = J - 2
0391:          END IF
0392:          DO 30 I = 1, IEND
0393:             TEMP = TEMP + ABS( S( I, J ) )
0394:             TEMP2 = TEMP2 + ABS( P( I, J ) )
0395:    30    CONTINUE
0396:          WORK( J ) = TEMP
0397:          WORK( N+J ) = TEMP2
0398:          DO 40 I = IEND + 1, MIN( J+1, N )
0399:             TEMP = TEMP + ABS( S( I, J ) )
0400:             TEMP2 = TEMP2 + ABS( P( I, J ) )
0401:    40    CONTINUE
0402:          ANORM = MAX( ANORM, TEMP )
0403:          BNORM = MAX( BNORM, TEMP2 )
0404:    50 CONTINUE
0405: *
0406:       ASCALE = ONE / MAX( ANORM, SAFMIN )
0407:       BSCALE = ONE / MAX( BNORM, SAFMIN )
0408: *
0409: *     Left eigenvectors
0410: *
0411:       IF( COMPL ) THEN
0412:          IEIG = 0
0413: *
0414: *        Main loop over eigenvalues
0415: *
0416:          ILCPLX = .FALSE.
0417:          DO 220 JE = 1, N
0418: *
0419: *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
0420: *           (b) this would be the second of a complex pair.
0421: *           Check for complex eigenvalue, so as to be sure of which
0422: *           entry(-ies) of SELECT to look at.
0423: *
0424:             IF( ILCPLX ) THEN
0425:                ILCPLX = .FALSE.
0426:                GO TO 220
0427:             END IF
0428:             NW = 1
0429:             IF( JE.LT.N ) THEN
0430:                IF( S( JE+1, JE ).NE.ZERO ) THEN
0431:                   ILCPLX = .TRUE.
0432:                   NW = 2
0433:                END IF
0434:             END IF
0435:             IF( ILALL ) THEN
0436:                ILCOMP = .TRUE.
0437:             ELSE IF( ILCPLX ) THEN
0438:                ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 )
0439:             ELSE
0440:                ILCOMP = SELECT( JE )
0441:             END IF
0442:             IF( .NOT.ILCOMP )
0443:      $         GO TO 220
0444: *
0445: *           Decide if (a) singular pencil, (b) real eigenvalue, or
0446: *           (c) complex eigenvalue.
0447: *
0448:             IF( .NOT.ILCPLX ) THEN
0449:                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
0450:      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
0451: *
0452: *                 Singular matrix pencil -- return unit eigenvector
0453: *
0454:                   IEIG = IEIG + 1
0455:                   DO 60 JR = 1, N
0456:                      VL( JR, IEIG ) = ZERO
0457:    60             CONTINUE
0458:                   VL( IEIG, IEIG ) = ONE
0459:                   GO TO 220
0460:                END IF
0461:             END IF
0462: *
0463: *           Clear vector
0464: *
0465:             DO 70 JR = 1, NW*N
0466:                WORK( 2*N+JR ) = ZERO
0467:    70       CONTINUE
0468: *                                                 T
0469: *           Compute coefficients in  ( a A - b B )  y = 0
0470: *              a  is  ACOEF
0471: *              b  is  BCOEFR + i*BCOEFI
0472: *
0473:             IF( .NOT.ILCPLX ) THEN
0474: *
0475: *              Real eigenvalue
0476: *
0477:                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
0478:      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
0479:                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
0480:                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
0481:                ACOEF = SBETA*ASCALE
0482:                BCOEFR = SALFAR*BSCALE
0483:                BCOEFI = ZERO
0484: *
0485: *              Scale to avoid underflow
0486: *
0487:                SCALE = ONE
0488:                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
0489:                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
0490:      $               SMALL
0491:                IF( LSA )
0492:      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
0493:                IF( LSB )
0494:      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
0495:      $                    MIN( BNORM, BIG ) )
0496:                IF( LSA .OR. LSB ) THEN
0497:                   SCALE = MIN( SCALE, ONE /
0498:      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
0499:      $                    ABS( BCOEFR ) ) ) )
0500:                   IF( LSA ) THEN
0501:                      ACOEF = ASCALE*( SCALE*SBETA )
0502:                   ELSE
0503:                      ACOEF = SCALE*ACOEF
0504:                   END IF
0505:                   IF( LSB ) THEN
0506:                      BCOEFR = BSCALE*( SCALE*SALFAR )
0507:                   ELSE
0508:                      BCOEFR = SCALE*BCOEFR
0509:                   END IF
0510:                END IF
0511:                ACOEFA = ABS( ACOEF )
0512:                BCOEFA = ABS( BCOEFR )
0513: *
0514: *              First component is 1
0515: *
0516:                WORK( 2*N+JE ) = ONE
0517:                XMAX = ONE
0518:             ELSE
0519: *
0520: *              Complex eigenvalue
0521: *
0522:                CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
0523:      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
0524:      $                     BCOEFI )
0525:                BCOEFI = -BCOEFI
0526:                IF( BCOEFI.EQ.ZERO ) THEN
0527:                   INFO = JE
0528:                   RETURN
0529:                END IF
0530: *
0531: *              Scale to avoid over/underflow
0532: *
0533:                ACOEFA = ABS( ACOEF )
0534:                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
0535:                SCALE = ONE
0536:                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
0537:      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
0538:                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
0539:      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
0540:                IF( SAFMIN*ACOEFA.GT.ASCALE )
0541:      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
0542:                IF( SAFMIN*BCOEFA.GT.BSCALE )
0543:      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
0544:                IF( SCALE.NE.ONE ) THEN
0545:                   ACOEF = SCALE*ACOEF
0546:                   ACOEFA = ABS( ACOEF )
0547:                   BCOEFR = SCALE*BCOEFR
0548:                   BCOEFI = SCALE*BCOEFI
0549:                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
0550:                END IF
0551: *
0552: *              Compute first two components of eigenvector
0553: *
0554:                TEMP = ACOEF*S( JE+1, JE )
0555:                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
0556:                TEMP2I = -BCOEFI*P( JE, JE )
0557:                IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
0558:                   WORK( 2*N+JE ) = ONE
0559:                   WORK( 3*N+JE ) = ZERO
0560:                   WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
0561:                   WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
0562:                ELSE
0563:                   WORK( 2*N+JE+1 ) = ONE
0564:                   WORK( 3*N+JE+1 ) = ZERO
0565:                   TEMP = ACOEF*S( JE, JE+1 )
0566:                   WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
0567:      $                             S( JE+1, JE+1 ) ) / TEMP
0568:                   WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
0569:                END IF
0570:                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
0571:      $                ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
0572:             END IF
0573: *
0574:             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
0575: *
0576: *                                           T
0577: *           Triangular solve of  (a A - b B)  y = 0
0578: *
0579: *                                   T
0580: *           (rowwise in  (a A - b B) , or columnwise in (a A - b B) )
0581: *
0582:             IL2BY2 = .FALSE.
0583: *
0584:             DO 160 J = JE + NW, N
0585:                IF( IL2BY2 ) THEN
0586:                   IL2BY2 = .FALSE.
0587:                   GO TO 160
0588:                END IF
0589: *
0590:                NA = 1
0591:                BDIAG( 1 ) = P( J, J )
0592:                IF( J.LT.N ) THEN
0593:                   IF( S( J+1, J ).NE.ZERO ) THEN
0594:                      IL2BY2 = .TRUE.
0595:                      BDIAG( 2 ) = P( J+1, J+1 )
0596:                      NA = 2
0597:                   END IF
0598:                END IF
0599: *
0600: *              Check whether scaling is necessary for dot products
0601: *
0602:                XSCALE = ONE / MAX( ONE, XMAX )
0603:                TEMP = MAX( WORK( J ), WORK( N+J ),
0604:      $                ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) )
0605:                IF( IL2BY2 )
0606:      $            TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ),
0607:      $                   ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
0608:                IF( TEMP.GT.BIGNUM*XSCALE ) THEN
0609:                   DO 90 JW = 0, NW - 1
0610:                      DO 80 JR = JE, J - 1
0611:                         WORK( ( JW+2 )*N+JR ) = XSCALE*
0612:      $                     WORK( ( JW+2 )*N+JR )
0613:    80                CONTINUE
0614:    90             CONTINUE
0615:                   XMAX = XMAX*XSCALE
0616:                END IF
0617: *
0618: *              Compute dot products
0619: *
0620: *                    j-1
0621: *              SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
0622: *                    k=je
0623: *
0624: *              To reduce the op count, this is done as
0625: *
0626: *              _        j-1                  _        j-1
0627: *              a*conjg( sum  S(k,j)*x(k) ) - b*conjg( sum  P(k,j)*x(k) )
0628: *                       k=je                          k=je
0629: *
0630: *              which may cause underflow problems if A or B are close
0631: *              to underflow.  (E.g., less than SMALL.)
0632: *
0633: *
0634: *              A series of compiler directives to defeat vectorization
0635: *              for the next loop
0636: *
0637: *$PL$ CMCHAR=' '
0638: CDIR$          NEXTSCALAR
0639: C$DIR          SCALAR
0640: CDIR$          NEXT SCALAR
0641: CVD$L          NOVECTOR
0642: CDEC$          NOVECTOR
0643: CVD$           NOVECTOR
0644: *VDIR          NOVECTOR
0645: *VOCL          LOOP,SCALAR
0646: CIBM           PREFER SCALAR
0647: *$PL$ CMCHAR='*'
0648: *
0649:                DO 120 JW = 1, NW
0650: *
0651: *$PL$ CMCHAR=' '
0652: CDIR$             NEXTSCALAR
0653: C$DIR             SCALAR
0654: CDIR$             NEXT SCALAR
0655: CVD$L             NOVECTOR
0656: CDEC$             NOVECTOR
0657: CVD$              NOVECTOR
0658: *VDIR             NOVECTOR
0659: *VOCL             LOOP,SCALAR
0660: CIBM              PREFER SCALAR
0661: *$PL$ CMCHAR='*'
0662: *
0663:                   DO 110 JA = 1, NA
0664:                      SUMS( JA, JW ) = ZERO
0665:                      SUMP( JA, JW ) = ZERO
0666: *
0667:                      DO 100 JR = JE, J - 1
0668:                         SUMS( JA, JW ) = SUMS( JA, JW ) +
0669:      $                                   S( JR, J+JA-1 )*
0670:      $                                   WORK( ( JW+1 )*N+JR )
0671:                         SUMP( JA, JW ) = SUMP( JA, JW ) +
0672:      $                                   P( JR, J+JA-1 )*
0673:      $                                   WORK( ( JW+1 )*N+JR )
0674:   100                CONTINUE
0675:   110             CONTINUE
0676:   120          CONTINUE
0677: *
0678: *$PL$ CMCHAR=' '
0679: CDIR$          NEXTSCALAR
0680: C$DIR          SCALAR
0681: CDIR$          NEXT SCALAR
0682: CVD$L          NOVECTOR
0683: CDEC$          NOVECTOR
0684: CVD$           NOVECTOR
0685: *VDIR          NOVECTOR
0686: *VOCL          LOOP,SCALAR
0687: CIBM           PREFER SCALAR
0688: *$PL$ CMCHAR='*'
0689: *
0690:                DO 130 JA = 1, NA
0691:                   IF( ILCPLX ) THEN
0692:                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
0693:      $                              BCOEFR*SUMP( JA, 1 ) -
0694:      $                              BCOEFI*SUMP( JA, 2 )
0695:                      SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
0696:      $                              BCOEFR*SUMP( JA, 2 ) +
0697:      $                              BCOEFI*SUMP( JA, 1 )
0698:                   ELSE
0699:                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
0700:      $                              BCOEFR*SUMP( JA, 1 )
0701:                   END IF
0702:   130          CONTINUE
0703: *
0704: *                                  T
0705: *              Solve  ( a A - b B )  y = SUM(,)
0706: *              with scaling and perturbation of the denominator
0707: *
0708:                CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
0709:      $                      BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
0710:      $                      BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
0711:      $                      IINFO )
0712:                IF( SCALE.LT.ONE ) THEN
0713:                   DO 150 JW = 0, NW - 1
0714:                      DO 140 JR = JE, J - 1
0715:                         WORK( ( JW+2 )*N+JR ) = SCALE*
0716:      $                     WORK( ( JW+2 )*N+JR )
0717:   140                CONTINUE
0718:   150             CONTINUE
0719:                   XMAX = SCALE*XMAX
0720:                END IF
0721:                XMAX = MAX( XMAX, TEMP )
0722:   160       CONTINUE
0723: *
0724: *           Copy eigenvector to VL, back transforming if
0725: *           HOWMNY='B'.
0726: *
0727:             IEIG = IEIG + 1
0728:             IF( ILBACK ) THEN
0729:                DO 170 JW = 0, NW - 1
0730:                   CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
0731:      $                        WORK( ( JW+2 )*N+JE ), 1, ZERO,
0732:      $                        WORK( ( JW+4 )*N+1 ), 1 )
0733:   170          CONTINUE
0734:                CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
0735:      $                      LDVL )
0736:                IBEG = 1
0737:             ELSE
0738:                CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
0739:      $                      LDVL )
0740:                IBEG = JE
0741:             END IF
0742: *
0743: *           Scale eigenvector
0744: *
0745:             XMAX = ZERO
0746:             IF( ILCPLX ) THEN
0747:                DO 180 J = IBEG, N
0748:                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
0749:      $                   ABS( VL( J, IEIG+1 ) ) )
0750:   180          CONTINUE
0751:             ELSE
0752:                DO 190 J = IBEG, N
0753:                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
0754:   190          CONTINUE
0755:             END IF
0756: *
0757:             IF( XMAX.GT.SAFMIN ) THEN
0758:                XSCALE = ONE / XMAX
0759: *
0760:                DO 210 JW = 0, NW - 1
0761:                   DO 200 JR = IBEG, N
0762:                      VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
0763:   200             CONTINUE
0764:   210          CONTINUE
0765:             END IF
0766:             IEIG = IEIG + NW - 1
0767: *
0768:   220    CONTINUE
0769:       END IF
0770: *
0771: *     Right eigenvectors
0772: *
0773:       IF( COMPR ) THEN
0774:          IEIG = IM + 1
0775: *
0776: *        Main loop over eigenvalues
0777: *
0778:          ILCPLX = .FALSE.
0779:          DO 500 JE = N, 1, -1
0780: *
0781: *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
0782: *           (b) this would be the second of a complex pair.
0783: *           Check for complex eigenvalue, so as to be sure of which
0784: *           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
0785: *           or SELECT(JE-1).
0786: *           If this is a complex pair, the 2-by-2 diagonal block
0787: *           corresponding to the eigenvalue is in rows/columns JE-1:JE
0788: *
0789:             IF( ILCPLX ) THEN
0790:                ILCPLX = .FALSE.
0791:                GO TO 500
0792:             END IF
0793:             NW = 1
0794:             IF( JE.GT.1 ) THEN
0795:                IF( S( JE, JE-1 ).NE.ZERO ) THEN
0796:                   ILCPLX = .TRUE.
0797:                   NW = 2
0798:                END IF
0799:             END IF
0800:             IF( ILALL ) THEN
0801:                ILCOMP = .TRUE.
0802:             ELSE IF( ILCPLX ) THEN
0803:                ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
0804:             ELSE
0805:                ILCOMP = SELECT( JE )
0806:             END IF
0807:             IF( .NOT.ILCOMP )
0808:      $         GO TO 500
0809: *
0810: *           Decide if (a) singular pencil, (b) real eigenvalue, or
0811: *           (c) complex eigenvalue.
0812: *
0813:             IF( .NOT.ILCPLX ) THEN
0814:                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
0815:      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
0816: *
0817: *                 Singular matrix pencil -- unit eigenvector
0818: *
0819:                   IEIG = IEIG - 1
0820:                   DO 230 JR = 1, N
0821:                      VR( JR, IEIG ) = ZERO
0822:   230             CONTINUE
0823:                   VR( IEIG, IEIG ) = ONE
0824:                   GO TO 500
0825:                END IF
0826:             END IF
0827: *
0828: *           Clear vector
0829: *
0830:             DO 250 JW = 0, NW - 1
0831:                DO 240 JR = 1, N
0832:                   WORK( ( JW+2 )*N+JR ) = ZERO
0833:   240          CONTINUE
0834:   250       CONTINUE
0835: *
0836: *           Compute coefficients in  ( a A - b B ) x = 0
0837: *              a  is  ACOEF
0838: *              b  is  BCOEFR + i*BCOEFI
0839: *
0840:             IF( .NOT.ILCPLX ) THEN
0841: *
0842: *              Real eigenvalue
0843: *
0844:                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
0845:      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
0846:                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
0847:                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
0848:                ACOEF = SBETA*ASCALE
0849:                BCOEFR = SALFAR*BSCALE
0850:                BCOEFI = ZERO
0851: *
0852: *              Scale to avoid underflow
0853: *
0854:                SCALE = ONE
0855:                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
0856:                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
0857:      $               SMALL
0858:                IF( LSA )
0859:      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
0860:                IF( LSB )
0861:      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
0862:      $                    MIN( BNORM, BIG ) )
0863:                IF( LSA .OR. LSB ) THEN
0864:                   SCALE = MIN( SCALE, ONE /
0865:      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
0866:      $                    ABS( BCOEFR ) ) ) )
0867:                   IF( LSA ) THEN
0868:                      ACOEF = ASCALE*( SCALE*SBETA )
0869:                   ELSE
0870:                      ACOEF = SCALE*ACOEF
0871:                   END IF
0872:                   IF( LSB ) THEN
0873:                      BCOEFR = BSCALE*( SCALE*SALFAR )
0874:                   ELSE
0875:                      BCOEFR = SCALE*BCOEFR
0876:                   END IF
0877:                END IF
0878:                ACOEFA = ABS( ACOEF )
0879:                BCOEFA = ABS( BCOEFR )
0880: *
0881: *              First component is 1
0882: *
0883:                WORK( 2*N+JE ) = ONE
0884:                XMAX = ONE
0885: *
0886: *              Compute contribution from column JE of A and B to sum
0887: *              (See "Further Details", above.)
0888: *
0889:                DO 260 JR = 1, JE - 1
0890:                   WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
0891:      $                             ACOEF*S( JR, JE )
0892:   260          CONTINUE
0893:             ELSE
0894: *
0895: *              Complex eigenvalue
0896: *
0897:                CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
0898:      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
0899:      $                     BCOEFI )
0900:                IF( BCOEFI.EQ.ZERO ) THEN
0901:                   INFO = JE - 1
0902:                   RETURN
0903:                END IF
0904: *
0905: *              Scale to avoid over/underflow
0906: *
0907:                ACOEFA = ABS( ACOEF )
0908:                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
0909:                SCALE = ONE
0910:                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
0911:      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
0912:                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
0913:      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
0914:                IF( SAFMIN*ACOEFA.GT.ASCALE )
0915:      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
0916:                IF( SAFMIN*BCOEFA.GT.BSCALE )
0917:      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
0918:                IF( SCALE.NE.ONE ) THEN
0919:                   ACOEF = SCALE*ACOEF
0920:                   ACOEFA = ABS( ACOEF )
0921:                   BCOEFR = SCALE*BCOEFR
0922:                   BCOEFI = SCALE*BCOEFI
0923:                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
0924:                END IF
0925: *
0926: *              Compute first two components of eigenvector
0927: *              and contribution to sums
0928: *
0929:                TEMP = ACOEF*S( JE, JE-1 )
0930:                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
0931:                TEMP2I = -BCOEFI*P( JE, JE )
0932:                IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
0933:                   WORK( 2*N+JE ) = ONE
0934:                   WORK( 3*N+JE ) = ZERO
0935:                   WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
0936:                   WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
0937:                ELSE
0938:                   WORK( 2*N+JE-1 ) = ONE
0939:                   WORK( 3*N+JE-1 ) = ZERO
0940:                   TEMP = ACOEF*S( JE-1, JE )
0941:                   WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
0942:      $                             S( JE-1, JE-1 ) ) / TEMP
0943:                   WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
0944:                END IF
0945: *
0946:                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
0947:      $                ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
0948: *
0949: *              Compute contribution from columns JE and JE-1
0950: *              of A and B to the sums.
0951: *
0952:                CREALA = ACOEF*WORK( 2*N+JE-1 )
0953:                CIMAGA = ACOEF*WORK( 3*N+JE-1 )
0954:                CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
0955:      $                  BCOEFI*WORK( 3*N+JE-1 )
0956:                CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
0957:      $                  BCOEFR*WORK( 3*N+JE-1 )
0958:                CRE2A = ACOEF*WORK( 2*N+JE )
0959:                CIM2A = ACOEF*WORK( 3*N+JE )
0960:                CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
0961:                CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
0962:                DO 270 JR = 1, JE - 2
0963:                   WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
0964:      $                             CREALB*P( JR, JE-1 ) -
0965:      $                             CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
0966:                   WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
0967:      $                             CIMAGB*P( JR, JE-1 ) -
0968:      $                             CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
0969:   270          CONTINUE
0970:             END IF
0971: *
0972:             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
0973: *
0974: *           Columnwise triangular solve of  (a A - b B)  x = 0
0975: *
0976:             IL2BY2 = .FALSE.
0977:             DO 370 J = JE - NW, 1, -1
0978: *
0979: *              If a 2-by-2 block, is in position j-1:j, wait until
0980: *              next iteration to process it (when it will be j:j+1)
0981: *
0982:                IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
0983:                   IF( S( J, J-1 ).NE.ZERO ) THEN
0984:                      IL2BY2 = .TRUE.
0985:                      GO TO 370
0986:                   END IF
0987:                END IF
0988:                BDIAG( 1 ) = P( J, J )
0989:                IF( IL2BY2 ) THEN
0990:                   NA = 2
0991:                   BDIAG( 2 ) = P( J+1, J+1 )
0992:                ELSE
0993:                   NA = 1
0994:                END IF
0995: *
0996: *              Compute x(j) (and x(j+1), if 2-by-2 block)
0997: *
0998:                CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
0999:      $                      LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
1000:      $                      N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
1001:      $                      IINFO )
1002:                IF( SCALE.LT.ONE ) THEN
1003: *
1004:                   DO 290 JW = 0, NW - 1
1005:                      DO 280 JR = 1, JE
1006:                         WORK( ( JW+2 )*N+JR ) = SCALE*
1007:      $                     WORK( ( JW+2 )*N+JR )
1008:   280                CONTINUE
1009:   290             CONTINUE
1010:                END IF
1011:                XMAX = MAX( SCALE*XMAX, TEMP )
1012: *
1013:                DO 310 JW = 1, NW
1014:                   DO 300 JA = 1, NA
1015:                      WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
1016:   300             CONTINUE
1017:   310          CONTINUE
1018: *
1019: *              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1020: *
1021:                IF( J.GT.1 ) THEN
1022: *
1023: *                 Check whether scaling is necessary for sum.
1024: *
1025:                   XSCALE = ONE / MAX( ONE, XMAX )
1026:                   TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
1027:                   IF( IL2BY2 )
1028:      $               TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
1029:      $                      WORK( N+J+1 ) )
1030:                   TEMP = MAX( TEMP, ACOEFA, BCOEFA )
1031:                   IF( TEMP.GT.BIGNUM*XSCALE ) THEN
1032: *
1033:                      DO 330 JW = 0, NW - 1
1034:                         DO 320 JR = 1, JE
1035:                            WORK( ( JW+2 )*N+JR ) = XSCALE*
1036:      $                        WORK( ( JW+2 )*N+JR )
1037:   320                   CONTINUE
1038:   330                CONTINUE
1039:                      XMAX = XMAX*XSCALE
1040:                   END IF
1041: *
1042: *                 Compute the contributions of the off-diagonals of
1043: *                 column j (and j+1, if 2-by-2 block) of A and B to the
1044: *                 sums.
1045: *
1046: *
1047:                   DO 360 JA = 1, NA
1048:                      IF( ILCPLX ) THEN
1049:                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1050:                         CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
1051:                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
1052:      $                           BCOEFI*WORK( 3*N+J+JA-1 )
1053:                         CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
1054:      $                           BCOEFR*WORK( 3*N+J+JA-1 )
1055:                         DO 340 JR = 1, J - 1
1056:                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1057:      $                                      CREALA*S( JR, J+JA-1 ) +
1058:      $                                      CREALB*P( JR, J+JA-1 )
1059:                            WORK( 3*N+JR ) = WORK( 3*N+JR ) -
1060:      $                                      CIMAGA*S( JR, J+JA-1 ) +
1061:      $                                      CIMAGB*P( JR, J+JA-1 )
1062:   340                   CONTINUE
1063:                      ELSE
1064:                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1065:                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
1066:                         DO 350 JR = 1, J - 1
1067:                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1068:      $                                      CREALA*S( JR, J+JA-1 ) +
1069:      $                                      CREALB*P( JR, J+JA-1 )
1070:   350                   CONTINUE
1071:                      END IF
1072:   360             CONTINUE
1073:                END IF
1074: *
1075:                IL2BY2 = .FALSE.
1076:   370       CONTINUE
1077: *
1078: *           Copy eigenvector to VR, back transforming if
1079: *           HOWMNY='B'.
1080: *
1081:             IEIG = IEIG - NW
1082:             IF( ILBACK ) THEN
1083: *
1084:                DO 410 JW = 0, NW - 1
1085:                   DO 380 JR = 1, N
1086:                      WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
1087:      $                                       VR( JR, 1 )
1088:   380             CONTINUE
1089: *
1090: *                 A series of compiler directives to defeat
1091: *                 vectorization for the next loop
1092: *
1093: *
1094:                   DO 400 JC = 2, JE
1095:                      DO 390 JR = 1, N
1096:                         WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
1097:      $                     WORK( ( JW+2 )*N+JC )*VR( JR, JC )
1098:   390                CONTINUE
1099:   400             CONTINUE
1100:   410          CONTINUE
1101: *
1102:                DO 430 JW = 0, NW - 1
1103:                   DO 420 JR = 1, N
1104:                      VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
1105:   420             CONTINUE
1106:   430          CONTINUE
1107: *
1108:                IEND = N
1109:             ELSE
1110:                DO 450 JW = 0, NW - 1
1111:                   DO 440 JR = 1, N
1112:                      VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
1113:   440             CONTINUE
1114:   450          CONTINUE
1115: *
1116:                IEND = JE
1117:             END IF
1118: *
1119: *           Scale eigenvector
1120: *
1121:             XMAX = ZERO
1122:             IF( ILCPLX ) THEN
1123:                DO 460 J = 1, IEND
1124:                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
1125:      $                   ABS( VR( J, IEIG+1 ) ) )
1126:   460          CONTINUE
1127:             ELSE
1128:                DO 470 J = 1, IEND
1129:                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
1130:   470          CONTINUE
1131:             END IF
1132: *
1133:             IF( XMAX.GT.SAFMIN ) THEN
1134:                XSCALE = ONE / XMAX
1135:                DO 490 JW = 0, NW - 1
1136:                   DO 480 JR = 1, IEND
1137:                      VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
1138:   480             CONTINUE
1139:   490          CONTINUE
1140:             END IF
1141:   500    CONTINUE
1142:       END IF
1143: *
1144:       RETURN
1145: *
1146: *     End of DTGEVC
1147: *
1148:       END
1149: