SRC/dtgevc.f

   1 *> \brief \b DTGEVC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DTGEVC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgevc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgevc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgevc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
  22 *                          LDVL, VR, LDVR, MM, M, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          HOWMNY, SIDE
  26 *       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       LOGICAL            SELECT( * )
  30 *       DOUBLE PRECISION   P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
  31 *      $                   VR( LDVR, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> DTGEVC computes some or all of the right and/or left eigenvectors of
  42 *> a pair of real matrices (S,P), where S is a quasi-triangular matrix
  43 *> and P is upper triangular.  Matrix pairs of this type are produced by
  44 *> the generalized Schur factorization of a matrix pair (A,B):
  45 *>
  46 *>    A = Q*S*Z**T,  B = Q*P*Z**T
  47 *>
  48 *> as computed by DGGHRD + DHGEQZ.
  49 *>
  50 *> The right eigenvector x and the left eigenvector y of (S,P)
  51 *> corresponding to an eigenvalue w are defined by:
  52 *>
  53 *>    S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
  54 *>
  55 *> where y**H denotes the conjugate tranpose of y.
  56 *> The eigenvalues are not input to this routine, but are computed
  57 *> directly from the diagonal blocks of S and P.
  58 *>
  59 *> This routine returns the matrices X and/or Y of right and left
  60 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
  61 *> where Z and Q are input matrices.
  62 *> If Q and Z are the orthogonal factors from the generalized Schur
  63 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
  64 *> are the matrices of right and left eigenvectors of (A,B).
  65 *>
  66 *> \endverbatim
  67 *
  68 *  Arguments:
  69 *  ==========
  70 *
  71 *> \param[in] SIDE
  72 *> \verbatim
  73 *>          SIDE is CHARACTER*1
  74 *>          = 'R': compute right eigenvectors only;
  75 *>          = 'L': compute left eigenvectors only;
  76 *>          = 'B': compute both right and left eigenvectors.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] HOWMNY
  80 *> \verbatim
  81 *>          HOWMNY is CHARACTER*1
  82 *>          = 'A': compute all right and/or left eigenvectors;
  83 *>          = 'B': compute all right and/or left eigenvectors,
  84 *>                 backtransformed by the matrices in VR and/or VL;
  85 *>          = 'S': compute selected right and/or left eigenvectors,
  86 *>                 specified by the logical array SELECT.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] SELECT
  90 *> \verbatim
  91 *>          SELECT is LOGICAL array, dimension (N)
  92 *>          If HOWMNY='S', SELECT specifies the eigenvectors to be
  93 *>          computed.  If w(j) is a real eigenvalue, the corresponding
  94 *>          real eigenvector is computed if SELECT(j) is .TRUE..
  95 *>          If w(j) and w(j+1) are the real and imaginary parts of a
  96 *>          complex eigenvalue, the corresponding complex eigenvector
  97 *>          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
  98 *>          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
  99 *>          set to .FALSE..
 100 *>          Not referenced if HOWMNY = 'A' or 'B'.
 101 *> \endverbatim
 102 *>
 103 *> \param[in] N
 104 *> \verbatim
 105 *>          N is INTEGER
 106 *>          The order of the matrices S and P.  N >= 0.
 107 *> \endverbatim
 108 *>
 109 *> \param[in] S
 110 *> \verbatim
 111 *>          S is DOUBLE PRECISION array, dimension (LDS,N)
 112 *>          The upper quasi-triangular matrix S from a generalized Schur
 113 *>          factorization, as computed by DHGEQZ.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] LDS
 117 *> \verbatim
 118 *>          LDS is INTEGER
 119 *>          The leading dimension of array S.  LDS >= max(1,N).
 120 *> \endverbatim
 121 *>
 122 *> \param[in] P
 123 *> \verbatim
 124 *>          P is DOUBLE PRECISION array, dimension (LDP,N)
 125 *>          The upper triangular matrix P from a generalized Schur
 126 *>          factorization, as computed by DHGEQZ.
 127 *>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
 128 *>          of S must be in positive diagonal form.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LDP
 132 *> \verbatim
 133 *>          LDP is INTEGER
 134 *>          The leading dimension of array P.  LDP >= max(1,N).
 135 *> \endverbatim
 136 *>
 137 *> \param[in,out] VL
 138 *> \verbatim
 139 *>          VL is DOUBLE PRECISION array, dimension (LDVL,MM)
 140 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
 141 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
 142 *>          of left Schur vectors returned by DHGEQZ).
 143 *>          On exit, if SIDE = 'L' or 'B', VL contains:
 144 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
 145 *>          if HOWMNY = 'B', the matrix Q*Y;
 146 *>          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
 147 *>                      SELECT, stored consecutively in the columns of
 148 *>                      VL, in the same order as their eigenvalues.
 149 *>
 150 *>          A complex eigenvector corresponding to a complex eigenvalue
 151 *>          is stored in two consecutive columns, the first holding the
 152 *>          real part, and the second the imaginary part.
 153 *>
 154 *>          Not referenced if SIDE = 'R'.
 155 *> \endverbatim
 156 *>
 157 *> \param[in] LDVL
 158 *> \verbatim
 159 *>          LDVL is INTEGER
 160 *>          The leading dimension of array VL.  LDVL >= 1, and if
 161 *>          SIDE = 'L' or 'B', LDVL >= N.
 162 *> \endverbatim
 163 *>
 164 *> \param[in,out] VR
 165 *> \verbatim
 166 *>          VR is DOUBLE PRECISION array, dimension (LDVR,MM)
 167 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
 168 *>          contain an N-by-N matrix Z (usually the orthogonal matrix Z
 169 *>          of right Schur vectors returned by DHGEQZ).
 170 *>
 171 *>          On exit, if SIDE = 'R' or 'B', VR contains:
 172 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
 173 *>          if HOWMNY = 'B' or 'b', the matrix Z*X;
 174 *>          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
 175 *>                      specified by SELECT, stored consecutively in the
 176 *>                      columns of VR, in the same order as their
 177 *>                      eigenvalues.
 178 *>
 179 *>          A complex eigenvector corresponding to a complex eigenvalue
 180 *>          is stored in two consecutive columns, the first holding the
 181 *>          real part and the second the imaginary part.
 182 *>
 183 *>          Not referenced if SIDE = 'L'.
 184 *> \endverbatim
 185 *>
 186 *> \param[in] LDVR
 187 *> \verbatim
 188 *>          LDVR is INTEGER
 189 *>          The leading dimension of the array VR.  LDVR >= 1, and if
 190 *>          SIDE = 'R' or 'B', LDVR >= N.
 191 *> \endverbatim
 192 *>
 193 *> \param[in] MM
 194 *> \verbatim
 195 *>          MM is INTEGER
 196 *>          The number of columns in the arrays VL and/or VR. MM >= M.
 197 *> \endverbatim
 198 *>
 199 *> \param[out] M
 200 *> \verbatim
 201 *>          M is INTEGER
 202 *>          The number of columns in the arrays VL and/or VR actually
 203 *>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
 204 *>          is set to N.  Each selected real eigenvector occupies one
 205 *>          column and each selected complex eigenvector occupies two
 206 *>          columns.
 207 *> \endverbatim
 208 *>
 209 *> \param[out] WORK
 210 *> \verbatim
 211 *>          WORK is DOUBLE PRECISION array, dimension (6*N)
 212 *> \endverbatim
 213 *>
 214 *> \param[out] INFO
 215 *> \verbatim
 216 *>          INFO is INTEGER
 217 *>          = 0:  successful exit.
 218 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 219 *>          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
 220 *>                eigenvalue.
 221 *> \endverbatim
 222 *
 223 *  Authors:
 224 *  ========
 225 *
 226 *> \author Univ. of Tennessee
 227 *> \author Univ. of California Berkeley
 228 *> \author Univ. of Colorado Denver
 229 *> \author NAG Ltd.
 230 *
 231 *> \date November 2011
 232 *
 233 *> \ingroup doubleGEcomputational
 234 *
 235 *> \par Further Details:
 236 *  =====================
 237 *>
 238 *> \verbatim
 239 *>
 240 *>  Allocation of workspace:
 241 *>  ---------- -- ---------
 242 *>
 243 *>     WORK( j ) = 1-norm of j-th column of A, above the diagonal
 244 *>     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
 245 *>     WORK( 2*N+1:3*N ) = real part of eigenvector
 246 *>     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
 247 *>     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
 248 *>     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
 249 *>
 250 *>  Rowwise vs. columnwise solution methods:
 251 *>  ------- --  ---------- -------- -------
 252 *>
 253 *>  Finding a generalized eigenvector consists basically of solving the
 254 *>  singular triangular system
 255 *>
 256 *>   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)
 257 *>
 258 *>  Consider finding the i-th right eigenvector (assume all eigenvalues
 259 *>  are real). The equation to be solved is:
 260 *>       n                   i
 261 *>  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
 262 *>      k=j                 k=j
 263 *>
 264 *>  where  C = (A - w B)  (The components v(i+1:n) are 0.)
 265 *>
 266 *>  The "rowwise" method is:
 267 *>
 268 *>  (1)  v(i) := 1
 269 *>  for j = i-1,. . .,1:
 270 *>                          i
 271 *>      (2) compute  s = - sum C(j,k) v(k)   and
 272 *>                        k=j+1
 273 *>
 274 *>      (3) v(j) := s / C(j,j)
 275 *>
 276 *>  Step 2 is sometimes called the "dot product" step, since it is an
 277 *>  inner product between the j-th row and the portion of the eigenvector
 278 *>  that has been computed so far.
 279 *>
 280 *>  The "columnwise" method consists basically in doing the sums
 281 *>  for all the rows in parallel.  As each v(j) is computed, the
 282 *>  contribution of v(j) times the j-th column of C is added to the
 283 *>  partial sums.  Since FORTRAN arrays are stored columnwise, this has
 284 *>  the advantage that at each step, the elements of C that are accessed
 285 *>  are adjacent to one another, whereas with the rowwise method, the
 286 *>  elements accessed at a step are spaced LDS (and LDP) words apart.
 287 *>
 288 *>  When finding left eigenvectors, the matrix in question is the
 289 *>  transpose of the one in storage, so the rowwise method then
 290 *>  actually accesses columns of A and B at each step, and so is the
 291 *>  preferred method.
 292 *> \endverbatim
 293 *>
 294 *  =====================================================================
 295       SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
 296      $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
 297 *
 298 *  -- LAPACK computational routine (version 3.4.0) --
 299 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 300 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 301 *     November 2011
 302 *
 303 *     .. Scalar Arguments ..
 304       CHARACTER          HOWMNY, SIDE
 305       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
 306 *     ..
 307 *     .. Array Arguments ..
 308       LOGICAL            SELECT( * )
 309       DOUBLE PRECISION   P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
 310      $                   VR( LDVR, * ), WORK( * )
 311 *     ..
 312 *
 313 *
 314 *  =====================================================================
 315 *
 316 *     .. Parameters ..
 317       DOUBLE PRECISION   ZERO, ONE, SAFETY
 318       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
 319      $                   SAFETY = 1.0D+2 )
 320 *     ..
 321 *     .. Local Scalars ..
 322       LOGICAL            COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
 323      $                   ILBBAD, ILCOMP, ILCPLX, LSA, LSB
 324       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
 325      $                   J, JA, JC, JE, JR, JW, NA, NW
 326       DOUBLE PRECISION   ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
 327      $                   BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
 328      $                   CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
 329      $                   CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
 330      $                   SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
 331      $                   XSCALE
 332 *     ..
 333 *     .. Local Arrays ..
 334       DOUBLE PRECISION   BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
 335      $                   SUMP( 2, 2 )
 336 *     ..
 337 *     .. External Functions ..
 338       LOGICAL            LSAME
 339       DOUBLE PRECISION   DLAMCH
 340       EXTERNAL           LSAME, DLAMCH
 341 *     ..
 342 *     .. External Subroutines ..
 343       EXTERNAL           DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA
 344 *     ..
 345 *     .. Intrinsic Functions ..
 346       INTRINSIC          ABS, MAX, MIN
 347 *     ..
 348 *     .. Executable Statements ..
 349 *
 350 *     Decode and Test the input parameters
 351 *
 352       IF( LSAME( HOWMNY, 'A' ) ) THEN
 353          IHWMNY = 1
 354          ILALL = .TRUE.
 355          ILBACK = .FALSE.
 356       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
 357          IHWMNY = 2
 358          ILALL = .FALSE.
 359          ILBACK = .FALSE.
 360       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
 361          IHWMNY = 3
 362          ILALL = .TRUE.
 363          ILBACK = .TRUE.
 364       ELSE
 365          IHWMNY = -1
 366          ILALL = .TRUE.
 367       END IF
 368 *
 369       IF( LSAME( SIDE, 'R' ) ) THEN
 370          ISIDE = 1
 371          COMPL = .FALSE.
 372          COMPR = .TRUE.
 373       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
 374          ISIDE = 2
 375          COMPL = .TRUE.
 376          COMPR = .FALSE.
 377       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
 378          ISIDE = 3
 379          COMPL = .TRUE.
 380          COMPR = .TRUE.
 381       ELSE
 382          ISIDE = -1
 383       END IF
 384 *
 385       INFO = 0
 386       IF( ISIDE.LT.0 ) THEN
 387          INFO = -1
 388       ELSE IF( IHWMNY.LT.0 ) THEN
 389          INFO = -2
 390       ELSE IF( N.LT.0 ) THEN
 391          INFO = -4
 392       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
 393          INFO = -6
 394       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
 395          INFO = -8
 396       END IF
 397       IF( INFO.NE.0 ) THEN
 398          CALL XERBLA( 'DTGEVC', -INFO )
 399          RETURN
 400       END IF
 401 *
 402 *     Count the number of eigenvectors to be computed
 403 *
 404       IF( .NOT.ILALL ) THEN
 405          IM = 0
 406          ILCPLX = .FALSE.
 407          DO 10 J = 1, N
 408             IF( ILCPLX ) THEN
 409                ILCPLX = .FALSE.
 410                GO TO 10
 411             END IF
 412             IF( J.LT.N ) THEN
 413                IF( S( J+1, J ).NE.ZERO )
 414      $            ILCPLX = .TRUE.
 415             END IF
 416             IF( ILCPLX ) THEN
 417                IF( SELECT( J ) .OR. SELECT( J+1 ) )
 418      $            IM = IM + 2
 419             ELSE
 420                IF( SELECT( J ) )
 421      $            IM = IM + 1
 422             END IF
 423    10    CONTINUE
 424       ELSE
 425          IM = N
 426       END IF
 427 *
 428 *     Check 2-by-2 diagonal blocks of A, B
 429 *
 430       ILABAD = .FALSE.
 431       ILBBAD = .FALSE.
 432       DO 20 J = 1, N - 1
 433          IF( S( J+1, J ).NE.ZERO ) THEN
 434             IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
 435      $          P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
 436             IF( J.LT.N-1 ) THEN
 437                IF( S( J+2, J+1 ).NE.ZERO )
 438      $            ILABAD = .TRUE.
 439             END IF
 440          END IF
 441    20 CONTINUE
 442 *
 443       IF( ILABAD ) THEN
 444          INFO = -5
 445       ELSE IF( ILBBAD ) THEN
 446          INFO = -7
 447       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
 448          INFO = -10
 449       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
 450          INFO = -12
 451       ELSE IF( MM.LT.IM ) THEN
 452          INFO = -13
 453       END IF
 454       IF( INFO.NE.0 ) THEN
 455          CALL XERBLA( 'DTGEVC', -INFO )
 456          RETURN
 457       END IF
 458 *
 459 *     Quick return if possible
 460 *
 461       M = IM
 462       IF( N.EQ.0 )
 463      $   RETURN
 464 *
 465 *     Machine Constants
 466 *
 467       SAFMIN = DLAMCH( 'Safe minimum' )
 468       BIG = ONE / SAFMIN
 469       CALL DLABAD( SAFMIN, BIG )
 470       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
 471       SMALL = SAFMIN*N / ULP
 472       BIG = ONE / SMALL
 473       BIGNUM = ONE / ( SAFMIN*N )
 474 *
 475 *     Compute the 1-norm of each column of the strictly upper triangular
 476 *     part (i.e., excluding all elements belonging to the diagonal
 477 *     blocks) of A and B to check for possible overflow in the
 478 *     triangular solver.
 479 *
 480       ANORM = ABS( S( 1, 1 ) )
 481       IF( N.GT.1 )
 482      $   ANORM = ANORM + ABS( S( 2, 1 ) )
 483       BNORM = ABS( P( 1, 1 ) )
 484       WORK( 1 ) = ZERO
 485       WORK( N+1 ) = ZERO
 486 *
 487       DO 50 J = 2, N
 488          TEMP = ZERO
 489          TEMP2 = ZERO
 490          IF( S( J, J-1 ).EQ.ZERO ) THEN
 491             IEND = J - 1
 492          ELSE
 493             IEND = J - 2
 494          END IF
 495          DO 30 I = 1, IEND
 496             TEMP = TEMP + ABS( S( I, J ) )
 497             TEMP2 = TEMP2 + ABS( P( I, J ) )
 498    30    CONTINUE
 499          WORK( J ) = TEMP
 500          WORK( N+J ) = TEMP2
 501          DO 40 I = IEND + 1, MIN( J+1, N )
 502             TEMP = TEMP + ABS( S( I, J ) )
 503             TEMP2 = TEMP2 + ABS( P( I, J ) )
 504    40    CONTINUE
 505          ANORM = MAX( ANORM, TEMP )
 506          BNORM = MAX( BNORM, TEMP2 )
 507    50 CONTINUE
 508 *
 509       ASCALE = ONE / MAX( ANORM, SAFMIN )
 510       BSCALE = ONE / MAX( BNORM, SAFMIN )
 511 *
 512 *     Left eigenvectors
 513 *
 514       IF( COMPL ) THEN
 515          IEIG = 0
 516 *
 517 *        Main loop over eigenvalues
 518 *
 519          ILCPLX = .FALSE.
 520          DO 220 JE = 1, N
 521 *
 522 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
 523 *           (b) this would be the second of a complex pair.
 524 *           Check for complex eigenvalue, so as to be sure of which
 525 *           entry(-ies) of SELECT to look at.
 526 *
 527             IF( ILCPLX ) THEN
 528                ILCPLX = .FALSE.
 529                GO TO 220
 530             END IF
 531             NW = 1
 532             IF( JE.LT.N ) THEN
 533                IF( S( JE+1, JE ).NE.ZERO ) THEN
 534                   ILCPLX = .TRUE.
 535                   NW = 2
 536                END IF
 537             END IF
 538             IF( ILALL ) THEN
 539                ILCOMP = .TRUE.
 540             ELSE IF( ILCPLX ) THEN
 541                ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 )
 542             ELSE
 543                ILCOMP = SELECT( JE )
 544             END IF
 545             IF( .NOT.ILCOMP )
 546      $         GO TO 220
 547 *
 548 *           Decide if (a) singular pencil, (b) real eigenvalue, or
 549 *           (c) complex eigenvalue.
 550 *
 551             IF( .NOT.ILCPLX ) THEN
 552                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
 553      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
 554 *
 555 *                 Singular matrix pencil -- return unit eigenvector
 556 *
 557                   IEIG = IEIG + 1
 558                   DO 60 JR = 1, N
 559                      VL( JR, IEIG ) = ZERO
 560    60             CONTINUE
 561                   VL( IEIG, IEIG ) = ONE
 562                   GO TO 220
 563                END IF
 564             END IF
 565 *
 566 *           Clear vector
 567 *
 568             DO 70 JR = 1, NW*N
 569                WORK( 2*N+JR ) = ZERO
 570    70       CONTINUE
 571 *                                                 T
 572 *           Compute coefficients in  ( a A - b B )  y = 0
 573 *              a  is  ACOEF
 574 *              b  is  BCOEFR + i*BCOEFI
 575 *
 576             IF( .NOT.ILCPLX ) THEN
 577 *
 578 *              Real eigenvalue
 579 *
 580                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
 581      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
 582                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
 583                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
 584                ACOEF = SBETA*ASCALE
 585                BCOEFR = SALFAR*BSCALE
 586                BCOEFI = ZERO
 587 *
 588 *              Scale to avoid underflow
 589 *
 590                SCALE = ONE
 591                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
 592                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
 593      $               SMALL
 594                IF( LSA )
 595      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
 596                IF( LSB )
 597      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
 598      $                    MIN( BNORM, BIG ) )
 599                IF( LSA .OR. LSB ) THEN
 600                   SCALE = MIN( SCALE, ONE /
 601      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
 602      $                    ABS( BCOEFR ) ) ) )
 603                   IF( LSA ) THEN
 604                      ACOEF = ASCALE*( SCALE*SBETA )
 605                   ELSE
 606                      ACOEF = SCALE*ACOEF
 607                   END IF
 608                   IF( LSB ) THEN
 609                      BCOEFR = BSCALE*( SCALE*SALFAR )
 610                   ELSE
 611                      BCOEFR = SCALE*BCOEFR
 612                   END IF
 613                END IF
 614                ACOEFA = ABS( ACOEF )
 615                BCOEFA = ABS( BCOEFR )
 616 *
 617 *              First component is 1
 618 *
 619                WORK( 2*N+JE ) = ONE
 620                XMAX = ONE
 621             ELSE
 622 *
 623 *              Complex eigenvalue
 624 *
 625                CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
 626      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
 627      $                     BCOEFI )
 628                BCOEFI = -BCOEFI
 629                IF( BCOEFI.EQ.ZERO ) THEN
 630                   INFO = JE
 631                   RETURN
 632                END IF
 633 *
 634 *              Scale to avoid over/underflow
 635 *
 636                ACOEFA = ABS( ACOEF )
 637                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
 638                SCALE = ONE
 639                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
 640      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
 641                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
 642      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
 643                IF( SAFMIN*ACOEFA.GT.ASCALE )
 644      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
 645                IF( SAFMIN*BCOEFA.GT.BSCALE )
 646      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
 647                IF( SCALE.NE.ONE ) THEN
 648                   ACOEF = SCALE*ACOEF
 649                   ACOEFA = ABS( ACOEF )
 650                   BCOEFR = SCALE*BCOEFR
 651                   BCOEFI = SCALE*BCOEFI
 652                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
 653                END IF
 654 *
 655 *              Compute first two components of eigenvector
 656 *
 657                TEMP = ACOEF*S( JE+1, JE )
 658                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
 659                TEMP2I = -BCOEFI*P( JE, JE )
 660                IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
 661                   WORK( 2*N+JE ) = ONE
 662                   WORK( 3*N+JE ) = ZERO
 663                   WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
 664                   WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
 665                ELSE
 666                   WORK( 2*N+JE+1 ) = ONE
 667                   WORK( 3*N+JE+1 ) = ZERO
 668                   TEMP = ACOEF*S( JE, JE+1 )
 669                   WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
 670      $                             S( JE+1, JE+1 ) ) / TEMP
 671                   WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
 672                END IF
 673                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
 674      $                ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
 675             END IF
 676 *
 677             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
 678 *
 679 *                                           T
 680 *           Triangular solve of  (a A - b B)  y = 0
 681 *
 682 *                                   T
 683 *           (rowwise in  (a A - b B) , or columnwise in (a A - b B) )
 684 *
 685             IL2BY2 = .FALSE.
 686 *
 687             DO 160 J = JE + NW, N
 688                IF( IL2BY2 ) THEN
 689                   IL2BY2 = .FALSE.
 690                   GO TO 160
 691                END IF
 692 *
 693                NA = 1
 694                BDIAG( 1 ) = P( J, J )
 695                IF( J.LT.N ) THEN
 696                   IF( S( J+1, J ).NE.ZERO ) THEN
 697                      IL2BY2 = .TRUE.
 698                      BDIAG( 2 ) = P( J+1, J+1 )
 699                      NA = 2
 700                   END IF
 701                END IF
 702 *
 703 *              Check whether scaling is necessary for dot products
 704 *
 705                XSCALE = ONE / MAX( ONE, XMAX )
 706                TEMP = MAX( WORK( J ), WORK( N+J ),
 707      $                ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) )
 708                IF( IL2BY2 )
 709      $            TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ),
 710      $                   ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
 711                IF( TEMP.GT.BIGNUM*XSCALE ) THEN
 712                   DO 90 JW = 0, NW - 1
 713                      DO 80 JR = JE, J - 1
 714                         WORK( ( JW+2 )*N+JR ) = XSCALE*
 715      $                     WORK( ( JW+2 )*N+JR )
 716    80                CONTINUE
 717    90             CONTINUE
 718                   XMAX = XMAX*XSCALE
 719                END IF
 720 *
 721 *              Compute dot products
 722 *
 723 *                    j-1
 724 *              SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
 725 *                    k=je
 726 *
 727 *              To reduce the op count, this is done as
 728 *
 729 *              _        j-1                  _        j-1
 730 *              a*conjg( sum  S(k,j)*x(k) ) - b*conjg( sum  P(k,j)*x(k) )
 731 *                       k=je                          k=je
 732 *
 733 *              which may cause underflow problems if A or B are close
 734 *              to underflow.  (E.g., less than SMALL.)
 735 *
 736 *
 737                DO 120 JW = 1, NW
 738                   DO 110 JA = 1, NA
 739                      SUMS( JA, JW ) = ZERO
 740                      SUMP( JA, JW ) = ZERO
 741 *
 742                      DO 100 JR = JE, J - 1
 743                         SUMS( JA, JW ) = SUMS( JA, JW ) +
 744      $                                   S( JR, J+JA-1 )*
 745      $                                   WORK( ( JW+1 )*N+JR )
 746                         SUMP( JA, JW ) = SUMP( JA, JW ) +
 747      $                                   P( JR, J+JA-1 )*
 748      $                                   WORK( ( JW+1 )*N+JR )
 749   100                CONTINUE
 750   110             CONTINUE
 751   120          CONTINUE
 752 *
 753                DO 130 JA = 1, NA
 754                   IF( ILCPLX ) THEN
 755                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
 756      $                              BCOEFR*SUMP( JA, 1 ) -
 757      $                              BCOEFI*SUMP( JA, 2 )
 758                      SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
 759      $                              BCOEFR*SUMP( JA, 2 ) +
 760      $                              BCOEFI*SUMP( JA, 1 )
 761                   ELSE
 762                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
 763      $                              BCOEFR*SUMP( JA, 1 )
 764                   END IF
 765   130          CONTINUE
 766 *
 767 *                                  T
 768 *              Solve  ( a A - b B )  y = SUM(,)
 769 *              with scaling and perturbation of the denominator
 770 *
 771                CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
 772      $                      BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
 773      $                      BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
 774      $                      IINFO )
 775                IF( SCALE.LT.ONE ) THEN
 776                   DO 150 JW = 0, NW - 1
 777                      DO 140 JR = JE, J - 1
 778                         WORK( ( JW+2 )*N+JR ) = SCALE*
 779      $                     WORK( ( JW+2 )*N+JR )
 780   140                CONTINUE
 781   150             CONTINUE
 782                   XMAX = SCALE*XMAX
 783                END IF
 784                XMAX = MAX( XMAX, TEMP )
 785   160       CONTINUE
 786 *
 787 *           Copy eigenvector to VL, back transforming if
 788 *           HOWMNY='B'.
 789 *
 790             IEIG = IEIG + 1
 791             IF( ILBACK ) THEN
 792                DO 170 JW = 0, NW - 1
 793                   CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
 794      $                        WORK( ( JW+2 )*N+JE ), 1, ZERO,
 795      $                        WORK( ( JW+4 )*N+1 ), 1 )
 796   170          CONTINUE
 797                CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
 798      $                      LDVL )
 799                IBEG = 1
 800             ELSE
 801                CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
 802      $                      LDVL )
 803                IBEG = JE
 804             END IF
 805 *
 806 *           Scale eigenvector
 807 *
 808             XMAX = ZERO
 809             IF( ILCPLX ) THEN
 810                DO 180 J = IBEG, N
 811                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
 812      $                   ABS( VL( J, IEIG+1 ) ) )
 813   180          CONTINUE
 814             ELSE
 815                DO 190 J = IBEG, N
 816                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
 817   190          CONTINUE
 818             END IF
 819 *
 820             IF( XMAX.GT.SAFMIN ) THEN
 821                XSCALE = ONE / XMAX
 822 *
 823                DO 210 JW = 0, NW - 1
 824                   DO 200 JR = IBEG, N
 825                      VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
 826   200             CONTINUE
 827   210          CONTINUE
 828             END IF
 829             IEIG = IEIG + NW - 1
 830 *
 831   220    CONTINUE
 832       END IF
 833 *
 834 *     Right eigenvectors
 835 *
 836       IF( COMPR ) THEN
 837          IEIG = IM + 1
 838 *
 839 *        Main loop over eigenvalues
 840 *
 841          ILCPLX = .FALSE.
 842          DO 500 JE = N, 1, -1
 843 *
 844 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
 845 *           (b) this would be the second of a complex pair.
 846 *           Check for complex eigenvalue, so as to be sure of which
 847 *           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
 848 *           or SELECT(JE-1).
 849 *           If this is a complex pair, the 2-by-2 diagonal block
 850 *           corresponding to the eigenvalue is in rows/columns JE-1:JE
 851 *
 852             IF( ILCPLX ) THEN
 853                ILCPLX = .FALSE.
 854                GO TO 500
 855             END IF
 856             NW = 1
 857             IF( JE.GT.1 ) THEN
 858                IF( S( JE, JE-1 ).NE.ZERO ) THEN
 859                   ILCPLX = .TRUE.
 860                   NW = 2
 861                END IF
 862             END IF
 863             IF( ILALL ) THEN
 864                ILCOMP = .TRUE.
 865             ELSE IF( ILCPLX ) THEN
 866                ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
 867             ELSE
 868                ILCOMP = SELECT( JE )
 869             END IF
 870             IF( .NOT.ILCOMP )
 871      $         GO TO 500
 872 *
 873 *           Decide if (a) singular pencil, (b) real eigenvalue, or
 874 *           (c) complex eigenvalue.
 875 *
 876             IF( .NOT.ILCPLX ) THEN
 877                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
 878      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
 879 *
 880 *                 Singular matrix pencil -- unit eigenvector
 881 *
 882                   IEIG = IEIG - 1
 883                   DO 230 JR = 1, N
 884                      VR( JR, IEIG ) = ZERO
 885   230             CONTINUE
 886                   VR( IEIG, IEIG ) = ONE
 887                   GO TO 500
 888                END IF
 889             END IF
 890 *
 891 *           Clear vector
 892 *
 893             DO 250 JW = 0, NW - 1
 894                DO 240 JR = 1, N
 895                   WORK( ( JW+2 )*N+JR ) = ZERO
 896   240          CONTINUE
 897   250       CONTINUE
 898 *
 899 *           Compute coefficients in  ( a A - b B ) x = 0
 900 *              a  is  ACOEF
 901 *              b  is  BCOEFR + i*BCOEFI
 902 *
 903             IF( .NOT.ILCPLX ) THEN
 904 *
 905 *              Real eigenvalue
 906 *
 907                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
 908      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
 909                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
 910                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
 911                ACOEF = SBETA*ASCALE
 912                BCOEFR = SALFAR*BSCALE
 913                BCOEFI = ZERO
 914 *
 915 *              Scale to avoid underflow
 916 *
 917                SCALE = ONE
 918                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
 919                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
 920      $               SMALL
 921                IF( LSA )
 922      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
 923                IF( LSB )
 924      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
 925      $                    MIN( BNORM, BIG ) )
 926                IF( LSA .OR. LSB ) THEN
 927                   SCALE = MIN( SCALE, ONE /
 928      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
 929      $                    ABS( BCOEFR ) ) ) )
 930                   IF( LSA ) THEN
 931                      ACOEF = ASCALE*( SCALE*SBETA )
 932                   ELSE
 933                      ACOEF = SCALE*ACOEF
 934                   END IF
 935                   IF( LSB ) THEN
 936                      BCOEFR = BSCALE*( SCALE*SALFAR )
 937                   ELSE
 938                      BCOEFR = SCALE*BCOEFR
 939                   END IF
 940                END IF
 941                ACOEFA = ABS( ACOEF )
 942                BCOEFA = ABS( BCOEFR )
 943 *
 944 *              First component is 1
 945 *
 946                WORK( 2*N+JE ) = ONE
 947                XMAX = ONE
 948 *
 949 *              Compute contribution from column JE of A and B to sum
 950 *              (See "Further Details", above.)
 951 *
 952                DO 260 JR = 1, JE - 1
 953                   WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
 954      $                             ACOEF*S( JR, JE )
 955   260          CONTINUE
 956             ELSE
 957 *
 958 *              Complex eigenvalue
 959 *
 960                CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
 961      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
 962      $                     BCOEFI )
 963                IF( BCOEFI.EQ.ZERO ) THEN
 964                   INFO = JE - 1
 965                   RETURN
 966                END IF
 967 *
 968 *              Scale to avoid over/underflow
 969 *
 970                ACOEFA = ABS( ACOEF )
 971                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
 972                SCALE = ONE
 973                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
 974      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
 975                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
 976      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
 977                IF( SAFMIN*ACOEFA.GT.ASCALE )
 978      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
 979                IF( SAFMIN*BCOEFA.GT.BSCALE )
 980      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
 981                IF( SCALE.NE.ONE ) THEN
 982                   ACOEF = SCALE*ACOEF
 983                   ACOEFA = ABS( ACOEF )
 984                   BCOEFR = SCALE*BCOEFR
 985                   BCOEFI = SCALE*BCOEFI
 986                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
 987                END IF
 988 *
 989 *              Compute first two components of eigenvector
 990 *              and contribution to sums
 991 *
 992                TEMP = ACOEF*S( JE, JE-1 )
 993                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
 994                TEMP2I = -BCOEFI*P( JE, JE )
 995                IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
 996                   WORK( 2*N+JE ) = ONE
 997                   WORK( 3*N+JE ) = ZERO
 998                   WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
 999                   WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
1000                ELSE
1001                   WORK( 2*N+JE-1 ) = ONE
1002                   WORK( 3*N+JE-1 ) = ZERO
1003                   TEMP = ACOEF*S( JE-1, JE )
1004                   WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
1005      $                             S( JE-1, JE-1 ) ) / TEMP
1006                   WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
1007                END IF
1008 *
1009                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
1010      $                ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
1011 *
1012 *              Compute contribution from columns JE and JE-1
1013 *              of A and B to the sums.
1014 *
1015                CREALA = ACOEF*WORK( 2*N+JE-1 )
1016                CIMAGA = ACOEF*WORK( 3*N+JE-1 )
1017                CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
1018      $                  BCOEFI*WORK( 3*N+JE-1 )
1019                CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
1020      $                  BCOEFR*WORK( 3*N+JE-1 )
1021                CRE2A = ACOEF*WORK( 2*N+JE )
1022                CIM2A = ACOEF*WORK( 3*N+JE )
1023                CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
1024                CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
1025                DO 270 JR = 1, JE - 2
1026                   WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
1027      $                             CREALB*P( JR, JE-1 ) -
1028      $                             CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
1029                   WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
1030      $                             CIMAGB*P( JR, JE-1 ) -
1031      $                             CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
1032   270          CONTINUE
1033             END IF
1034 *
1035             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
1036 *
1037 *           Columnwise triangular solve of  (a A - b B)  x = 0
1038 *
1039             IL2BY2 = .FALSE.
1040             DO 370 J = JE - NW, 1, -1
1041 *
1042 *              If a 2-by-2 block, is in position j-1:j, wait until
1043 *              next iteration to process it (when it will be j:j+1)
1044 *
1045                IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
1046                   IF( S( J, J-1 ).NE.ZERO ) THEN
1047                      IL2BY2 = .TRUE.
1048                      GO TO 370
1049                   END IF
1050                END IF
1051                BDIAG( 1 ) = P( J, J )
1052                IF( IL2BY2 ) THEN
1053                   NA = 2
1054                   BDIAG( 2 ) = P( J+1, J+1 )
1055                ELSE
1056                   NA = 1
1057                END IF
1058 *
1059 *              Compute x(j) (and x(j+1), if 2-by-2 block)
1060 *
1061                CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
1062      $                      LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
1063      $                      N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
1064      $                      IINFO )
1065                IF( SCALE.LT.ONE ) THEN
1066 *
1067                   DO 290 JW = 0, NW - 1
1068                      DO 280 JR = 1, JE
1069                         WORK( ( JW+2 )*N+JR ) = SCALE*
1070      $                     WORK( ( JW+2 )*N+JR )
1071   280                CONTINUE
1072   290             CONTINUE
1073                END IF
1074                XMAX = MAX( SCALE*XMAX, TEMP )
1075 *
1076                DO 310 JW = 1, NW
1077                   DO 300 JA = 1, NA
1078                      WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
1079   300             CONTINUE
1080   310          CONTINUE
1081 *
1082 *              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1083 *
1084                IF( J.GT.1 ) THEN
1085 *
1086 *                 Check whether scaling is necessary for sum.
1087 *
1088                   XSCALE = ONE / MAX( ONE, XMAX )
1089                   TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
1090                   IF( IL2BY2 )
1091      $               TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
1092      $                      WORK( N+J+1 ) )
1093                   TEMP = MAX( TEMP, ACOEFA, BCOEFA )
1094                   IF( TEMP.GT.BIGNUM*XSCALE ) THEN
1095 *
1096                      DO 330 JW = 0, NW - 1
1097                         DO 320 JR = 1, JE
1098                            WORK( ( JW+2 )*N+JR ) = XSCALE*
1099      $                        WORK( ( JW+2 )*N+JR )
1100   320                   CONTINUE
1101   330                CONTINUE
1102                      XMAX = XMAX*XSCALE
1103                   END IF
1104 *
1105 *                 Compute the contributions of the off-diagonals of
1106 *                 column j (and j+1, if 2-by-2 block) of A and B to the
1107 *                 sums.
1108 *
1109 *
1110                   DO 360 JA = 1, NA
1111                      IF( ILCPLX ) THEN
1112                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1113                         CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
1114                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
1115      $                           BCOEFI*WORK( 3*N+J+JA-1 )
1116                         CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
1117      $                           BCOEFR*WORK( 3*N+J+JA-1 )
1118                         DO 340 JR = 1, J - 1
1119                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1120      $                                      CREALA*S( JR, J+JA-1 ) +
1121      $                                      CREALB*P( JR, J+JA-1 )
1122                            WORK( 3*N+JR ) = WORK( 3*N+JR ) -
1123      $                                      CIMAGA*S( JR, J+JA-1 ) +
1124      $                                      CIMAGB*P( JR, J+JA-1 )
1125   340                   CONTINUE
1126                      ELSE
1127                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1128                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
1129                         DO 350 JR = 1, J - 1
1130                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1131      $                                      CREALA*S( JR, J+JA-1 ) +
1132      $                                      CREALB*P( JR, J+JA-1 )
1133   350                   CONTINUE
1134                      END IF
1135   360             CONTINUE
1136                END IF
1137 *
1138                IL2BY2 = .FALSE.
1139   370       CONTINUE
1140 *
1141 *           Copy eigenvector to VR, back transforming if
1142 *           HOWMNY='B'.
1143 *
1144             IEIG = IEIG - NW
1145             IF( ILBACK ) THEN
1146 *
1147                DO 410 JW = 0, NW - 1
1148                   DO 380 JR = 1, N
1149                      WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
1150      $                                       VR( JR, 1 )
1151   380             CONTINUE
1152 *
1153 *                 A series of compiler directives to defeat
1154 *                 vectorization for the next loop
1155 *
1156 *
1157                   DO 400 JC = 2, JE
1158                      DO 390 JR = 1, N
1159                         WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
1160      $                     WORK( ( JW+2 )*N+JC )*VR( JR, JC )
1161   390                CONTINUE
1162   400             CONTINUE
1163   410          CONTINUE
1164 *
1165                DO 430 JW = 0, NW - 1
1166                   DO 420 JR = 1, N
1167                      VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
1168   420             CONTINUE
1169   430          CONTINUE
1170 *
1171                IEND = N
1172             ELSE
1173                DO 450 JW = 0, NW - 1
1174                   DO 440 JR = 1, N
1175                      VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
1176   440             CONTINUE
1177   450          CONTINUE
1178 *
1179                IEND = JE
1180             END IF
1181 *
1182 *           Scale eigenvector
1183 *
1184             XMAX = ZERO
1185             IF( ILCPLX ) THEN
1186                DO 460 J = 1, IEND
1187                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
1188      $                   ABS( VR( J, IEIG+1 ) ) )
1189   460          CONTINUE
1190             ELSE
1191                DO 470 J = 1, IEND
1192                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
1193   470          CONTINUE
1194             END IF
1195 *
1196             IF( XMAX.GT.SAFMIN ) THEN
1197                XSCALE = ONE / XMAX
1198                DO 490 JW = 0, NW - 1
1199                   DO 480 JR = 1, IEND
1200                      VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
1201   480             CONTINUE
1202   490          CONTINUE
1203             END IF
1204   500    CONTINUE
1205       END IF
1206 *
1207       RETURN
1208 *
1209 *     End of DTGEVC
1210 *
1211       END