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