LAPACK: SRC/dtgevc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00002      $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          HOWMNY, SIDE
00011       INTEGER            INFO, LDP, LDS, LDVL, LDVR, M, MM, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            SELECT( * )
00015       DOUBLE PRECISION   P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00016      $                   VR( LDVR, * ), WORK( * )
00017 *     ..
00018 *
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DTGEVC computes some or all of the right and/or left eigenvectors of
00024 *  a pair of real matrices (S,P), where S is a quasi-triangular matrix
00025 *  and P is upper triangular.  Matrix pairs of this type are produced by
00026 *  the generalized Schur factorization of a matrix pair (A,B):
00027 *
00028 *     A = Q*S*Z**T,  B = Q*P*Z**T
00029 *
00030 *  as computed by DGGHRD + DHGEQZ.
00031 *
00032 *  The right eigenvector x and the left eigenvector y of (S,P)
00033 *  corresponding to an eigenvalue w are defined by:
00034 *  
00035 *     S*x = w*P*x,  (y**H)*S = w*(y**H)*P,
00036 *  
00037 *  where y**H denotes the conjugate tranpose of y.
00038 *  The eigenvalues are not input to this routine, but are computed
00039 *  directly from the diagonal blocks of S and P.
00040 *  
00041 *  This routine returns the matrices X and/or Y of right and left
00042 *  eigenvectors of (S,P), or the products Z*X and/or Q*Y,
00043 *  where Z and Q are input matrices.
00044 *  If Q and Z are the orthogonal factors from the generalized Schur
00045 *  factorization of a matrix pair (A,B), then Z*X and Q*Y
00046 *  are the matrices of right and left eigenvectors of (A,B).
00047 * 
00048 *  Arguments
00049 *  =========
00050 *
00051 *  SIDE    (input) CHARACTER*1
00052 *          = 'R': compute right eigenvectors only;
00053 *          = 'L': compute left eigenvectors only;
00054 *          = 'B': compute both right and left eigenvectors.
00055 *
00056 *  HOWMNY  (input) CHARACTER*1
00057 *          = 'A': compute all right and/or left eigenvectors;
00058 *          = 'B': compute all right and/or left eigenvectors,
00059 *                 backtransformed by the matrices in VR and/or VL;
00060 *          = 'S': compute selected right and/or left eigenvectors,
00061 *                 specified by the logical array SELECT.
00062 *
00063 *  SELECT  (input) LOGICAL array, dimension (N)
00064 *          If HOWMNY='S', SELECT specifies the eigenvectors to be
00065 *          computed.  If w(j) is a real eigenvalue, the corresponding
00066 *          real eigenvector is computed if SELECT(j) is .TRUE..
00067 *          If w(j) and w(j+1) are the real and imaginary parts of a
00068 *          complex eigenvalue, the corresponding complex eigenvector
00069 *          is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
00070 *          and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
00071 *          set to .FALSE..
00072 *          Not referenced if HOWMNY = 'A' or 'B'.
00073 *
00074 *  N       (input) INTEGER
00075 *          The order of the matrices S and P.  N >= 0.
00076 *
00077 *  S       (input) DOUBLE PRECISION array, dimension (LDS,N)
00078 *          The upper quasi-triangular matrix S from a generalized Schur
00079 *          factorization, as computed by DHGEQZ.
00080 *
00081 *  LDS     (input) INTEGER
00082 *          The leading dimension of array S.  LDS >= max(1,N).
00083 *
00084 *  P       (input) DOUBLE PRECISION array, dimension (LDP,N)
00085 *          The upper triangular matrix P from a generalized Schur
00086 *          factorization, as computed by DHGEQZ.
00087 *          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
00088 *          of S must be in positive diagonal form.
00089 *
00090 *  LDP     (input) INTEGER
00091 *          The leading dimension of array P.  LDP >= max(1,N).
00092 *
00093 *  VL      (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
00094 *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00095 *          contain an N-by-N matrix Q (usually the orthogonal matrix Q
00096 *          of left Schur vectors returned by DHGEQZ).
00097 *          On exit, if SIDE = 'L' or 'B', VL contains:
00098 *          if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
00099 *          if HOWMNY = 'B', the matrix Q*Y;
00100 *          if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
00101 *                      SELECT, stored consecutively in the columns of
00102 *                      VL, in the same order as their eigenvalues.
00103 *
00104 *          A complex eigenvector corresponding to a complex eigenvalue
00105 *          is stored in two consecutive columns, the first holding the
00106 *          real part, and the second the imaginary part.
00107 *
00108 *          Not referenced if SIDE = 'R'.
00109 *
00110 *  LDVL    (input) INTEGER
00111 *          The leading dimension of array VL.  LDVL >= 1, and if
00112 *          SIDE = 'L' or 'B', LDVL >= N.
00113 *
00114 *  VR      (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
00115 *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00116 *          contain an N-by-N matrix Z (usually the orthogonal matrix Z
00117 *          of right Schur vectors returned by DHGEQZ).
00118 *
00119 *          On exit, if SIDE = 'R' or 'B', VR contains:
00120 *          if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
00121 *          if HOWMNY = 'B' or 'b', the matrix Z*X;
00122 *          if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
00123 *                      specified by SELECT, stored consecutively in the
00124 *                      columns of VR, in the same order as their
00125 *                      eigenvalues.
00126 *
00127 *          A complex eigenvector corresponding to a complex eigenvalue
00128 *          is stored in two consecutive columns, the first holding the
00129 *          real part and the second the imaginary part.
00130 *          
00131 *          Not referenced if SIDE = 'L'.
00132 *
00133 *  LDVR    (input) INTEGER
00134 *          The leading dimension of the array VR.  LDVR >= 1, and if
00135 *          SIDE = 'R' or 'B', LDVR >= N.
00136 *
00137 *  MM      (input) INTEGER
00138 *          The number of columns in the arrays VL and/or VR. MM >= M.
00139 *
00140 *  M       (output) INTEGER
00141 *          The number of columns in the arrays VL and/or VR actually
00142 *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
00143 *          is set to N.  Each selected real eigenvector occupies one
00144 *          column and each selected complex eigenvector occupies two
00145 *          columns.
00146 *
00147 *  WORK    (workspace) DOUBLE PRECISION array, dimension (6*N)
00148 *
00149 *  INFO    (output) INTEGER
00150 *          = 0:  successful exit.
00151 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00152 *          > 0:  the 2-by-2 block (INFO:INFO+1) does not have a complex
00153 *                eigenvalue.
00154 *
00155 *  Further Details
00156 *  ===============
00157 *
00158 *  Allocation of workspace:
00159 *  ---------- -- ---------
00160 *
00161 *     WORK( j ) = 1-norm of j-th column of A, above the diagonal
00162 *     WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
00163 *     WORK( 2*N+1:3*N ) = real part of eigenvector
00164 *     WORK( 3*N+1:4*N ) = imaginary part of eigenvector
00165 *     WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
00166 *     WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
00167 *
00168 *  Rowwise vs. columnwise solution methods:
00169 *  ------- --  ---------- -------- -------
00170 *
00171 *  Finding a generalized eigenvector consists basically of solving the
00172 *  singular triangular system
00173 *
00174 *   (A - w B) x = 0     (for right) or:   (A - w B)**H y = 0  (for left)
00175 *
00176 *  Consider finding the i-th right eigenvector (assume all eigenvalues
00177 *  are real). The equation to be solved is:
00178 *       n                   i
00179 *  0 = sum  C(j,k) v(k)  = sum  C(j,k) v(k)     for j = i,. . .,1
00180 *      k=j                 k=j
00181 *
00182 *  where  C = (A - w B)  (The components v(i+1:n) are 0.)
00183 *
00184 *  The "rowwise" method is:
00185 *
00186 *  (1)  v(i) := 1
00187 *  for j = i-1,. . .,1:
00188 *                          i
00189 *      (2) compute  s = - sum C(j,k) v(k)   and
00190 *                        k=j+1
00191 *
00192 *      (3) v(j) := s / C(j,j)
00193 *
00194 *  Step 2 is sometimes called the "dot product" step, since it is an
00195 *  inner product between the j-th row and the portion of the eigenvector
00196 *  that has been computed so far.
00197 *
00198 *  The "columnwise" method consists basically in doing the sums
00199 *  for all the rows in parallel.  As each v(j) is computed, the
00200 *  contribution of v(j) times the j-th column of C is added to the
00201 *  partial sums.  Since FORTRAN arrays are stored columnwise, this has
00202 *  the advantage that at each step, the elements of C that are accessed
00203 *  are adjacent to one another, whereas with the rowwise method, the
00204 *  elements accessed at a step are spaced LDS (and LDP) words apart.
00205 *
00206 *  When finding left eigenvectors, the matrix in question is the
00207 *  transpose of the one in storage, so the rowwise method then
00208 *  actually accesses columns of A and B at each step, and so is the
00209 *  preferred method.
00210 *
00211 *  =====================================================================
00212 *
00213 *     .. Parameters ..
00214       DOUBLE PRECISION   ZERO, ONE, SAFETY
00215       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
00216      $                   SAFETY = 1.0D+2 )
00217 *     ..
00218 *     .. Local Scalars ..
00219       LOGICAL            COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
00220      $                   ILBBAD, ILCOMP, ILCPLX, LSA, LSB
00221       INTEGER            I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
00222      $                   J, JA, JC, JE, JR, JW, NA, NW
00223       DOUBLE PRECISION   ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
00224      $                   BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A,
00225      $                   CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA,
00226      $                   CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
00227      $                   SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
00228      $                   XSCALE
00229 *     ..
00230 *     .. Local Arrays ..
00231       DOUBLE PRECISION   BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
00232      $                   SUMP( 2, 2 )
00233 *     ..
00234 *     .. External Functions ..
00235       LOGICAL            LSAME
00236       DOUBLE PRECISION   DLAMCH
00237       EXTERNAL           LSAME, DLAMCH
00238 *     ..
00239 *     .. External Subroutines ..
00240       EXTERNAL           DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA
00241 *     ..
00242 *     .. Intrinsic Functions ..
00243       INTRINSIC          ABS, MAX, MIN
00244 *     ..
00245 *     .. Executable Statements ..
00246 *
00247 *     Decode and Test the input parameters
00248 *
00249       IF( LSAME( HOWMNY, 'A' ) ) THEN
00250          IHWMNY = 1
00251          ILALL = .TRUE.
00252          ILBACK = .FALSE.
00253       ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
00254          IHWMNY = 2
00255          ILALL = .FALSE.
00256          ILBACK = .FALSE.
00257       ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
00258          IHWMNY = 3
00259          ILALL = .TRUE.
00260          ILBACK = .TRUE.
00261       ELSE
00262          IHWMNY = -1
00263          ILALL = .TRUE.
00264       END IF
00265 *
00266       IF( LSAME( SIDE, 'R' ) ) THEN
00267          ISIDE = 1
00268          COMPL = .FALSE.
00269          COMPR = .TRUE.
00270       ELSE IF( LSAME( SIDE, 'L' ) ) THEN
00271          ISIDE = 2
00272          COMPL = .TRUE.
00273          COMPR = .FALSE.
00274       ELSE IF( LSAME( SIDE, 'B' ) ) THEN
00275          ISIDE = 3
00276          COMPL = .TRUE.
00277          COMPR = .TRUE.
00278       ELSE
00279          ISIDE = -1
00280       END IF
00281 *
00282       INFO = 0
00283       IF( ISIDE.LT.0 ) THEN
00284          INFO = -1
00285       ELSE IF( IHWMNY.LT.0 ) THEN
00286          INFO = -2
00287       ELSE IF( N.LT.0 ) THEN
00288          INFO = -4
00289       ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
00290          INFO = -6
00291       ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
00292          INFO = -8
00293       END IF
00294       IF( INFO.NE.0 ) THEN
00295          CALL XERBLA( 'DTGEVC', -INFO )
00296          RETURN
00297       END IF
00298 *
00299 *     Count the number of eigenvectors to be computed
00300 *
00301       IF( .NOT.ILALL ) THEN
00302          IM = 0
00303          ILCPLX = .FALSE.
00304          DO 10 J = 1, N
00305             IF( ILCPLX ) THEN
00306                ILCPLX = .FALSE.
00307                GO TO 10
00308             END IF
00309             IF( J.LT.N ) THEN
00310                IF( S( J+1, J ).NE.ZERO )
00311      $            ILCPLX = .TRUE.
00312             END IF
00313             IF( ILCPLX ) THEN
00314                IF( SELECT( J ) .OR. SELECT( J+1 ) )
00315      $            IM = IM + 2
00316             ELSE
00317                IF( SELECT( J ) )
00318      $            IM = IM + 1
00319             END IF
00320    10    CONTINUE
00321       ELSE
00322          IM = N
00323       END IF
00324 *
00325 *     Check 2-by-2 diagonal blocks of A, B
00326 *
00327       ILABAD = .FALSE.
00328       ILBBAD = .FALSE.
00329       DO 20 J = 1, N - 1
00330          IF( S( J+1, J ).NE.ZERO ) THEN
00331             IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
00332      $          P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
00333             IF( J.LT.N-1 ) THEN
00334                IF( S( J+2, J+1 ).NE.ZERO )
00335      $            ILABAD = .TRUE.
00336             END IF
00337          END IF
00338    20 CONTINUE
00339 *
00340       IF( ILABAD ) THEN
00341          INFO = -5
00342       ELSE IF( ILBBAD ) THEN
00343          INFO = -7
00344       ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
00345          INFO = -10
00346       ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
00347          INFO = -12
00348       ELSE IF( MM.LT.IM ) THEN
00349          INFO = -13
00350       END IF
00351       IF( INFO.NE.0 ) THEN
00352          CALL XERBLA( 'DTGEVC', -INFO )
00353          RETURN
00354       END IF
00355 *
00356 *     Quick return if possible
00357 *
00358       M = IM
00359       IF( N.EQ.0 )
00360      $   RETURN
00361 *
00362 *     Machine Constants
00363 *
00364       SAFMIN = DLAMCH( 'Safe minimum' )
00365       BIG = ONE / SAFMIN
00366       CALL DLABAD( SAFMIN, BIG )
00367       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00368       SMALL = SAFMIN*N / ULP
00369       BIG = ONE / SMALL
00370       BIGNUM = ONE / ( SAFMIN*N )
00371 *
00372 *     Compute the 1-norm of each column of the strictly upper triangular
00373 *     part (i.e., excluding all elements belonging to the diagonal
00374 *     blocks) of A and B to check for possible overflow in the
00375 *     triangular solver.
00376 *
00377       ANORM = ABS( S( 1, 1 ) )
00378       IF( N.GT.1 )
00379      $   ANORM = ANORM + ABS( S( 2, 1 ) )
00380       BNORM = ABS( P( 1, 1 ) )
00381       WORK( 1 ) = ZERO
00382       WORK( N+1 ) = ZERO
00383 *
00384       DO 50 J = 2, N
00385          TEMP = ZERO
00386          TEMP2 = ZERO
00387          IF( S( J, J-1 ).EQ.ZERO ) THEN
00388             IEND = J - 1
00389          ELSE
00390             IEND = J - 2
00391          END IF
00392          DO 30 I = 1, IEND
00393             TEMP = TEMP + ABS( S( I, J ) )
00394             TEMP2 = TEMP2 + ABS( P( I, J ) )
00395    30    CONTINUE
00396          WORK( J ) = TEMP
00397          WORK( N+J ) = TEMP2
00398          DO 40 I = IEND + 1, MIN( J+1, N )
00399             TEMP = TEMP + ABS( S( I, J ) )
00400             TEMP2 = TEMP2 + ABS( P( I, J ) )
00401    40    CONTINUE
00402          ANORM = MAX( ANORM, TEMP )
00403          BNORM = MAX( BNORM, TEMP2 )
00404    50 CONTINUE
00405 *
00406       ASCALE = ONE / MAX( ANORM, SAFMIN )
00407       BSCALE = ONE / MAX( BNORM, SAFMIN )
00408 *
00409 *     Left eigenvectors
00410 *
00411       IF( COMPL ) THEN
00412          IEIG = 0
00413 *
00414 *        Main loop over eigenvalues
00415 *
00416          ILCPLX = .FALSE.
00417          DO 220 JE = 1, N
00418 *
00419 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
00420 *           (b) this would be the second of a complex pair.
00421 *           Check for complex eigenvalue, so as to be sure of which
00422 *           entry(-ies) of SELECT to look at.
00423 *
00424             IF( ILCPLX ) THEN
00425                ILCPLX = .FALSE.
00426                GO TO 220
00427             END IF
00428             NW = 1
00429             IF( JE.LT.N ) THEN
00430                IF( S( JE+1, JE ).NE.ZERO ) THEN
00431                   ILCPLX = .TRUE.
00432                   NW = 2
00433                END IF
00434             END IF
00435             IF( ILALL ) THEN
00436                ILCOMP = .TRUE.
00437             ELSE IF( ILCPLX ) THEN
00438                ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 )
00439             ELSE
00440                ILCOMP = SELECT( JE )
00441             END IF
00442             IF( .NOT.ILCOMP )
00443      $         GO TO 220
00444 *
00445 *           Decide if (a) singular pencil, (b) real eigenvalue, or
00446 *           (c) complex eigenvalue.
00447 *
00448             IF( .NOT.ILCPLX ) THEN
00449                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
00450      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
00451 *
00452 *                 Singular matrix pencil -- return unit eigenvector
00453 *
00454                   IEIG = IEIG + 1
00455                   DO 60 JR = 1, N
00456                      VL( JR, IEIG ) = ZERO
00457    60             CONTINUE
00458                   VL( IEIG, IEIG ) = ONE
00459                   GO TO 220
00460                END IF
00461             END IF
00462 *
00463 *           Clear vector
00464 *
00465             DO 70 JR = 1, NW*N
00466                WORK( 2*N+JR ) = ZERO
00467    70       CONTINUE
00468 *                                                 T
00469 *           Compute coefficients in  ( a A - b B )  y = 0
00470 *              a  is  ACOEF
00471 *              b  is  BCOEFR + i*BCOEFI
00472 *
00473             IF( .NOT.ILCPLX ) THEN
00474 *
00475 *              Real eigenvalue
00476 *
00477                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
00478      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
00479                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
00480                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
00481                ACOEF = SBETA*ASCALE
00482                BCOEFR = SALFAR*BSCALE
00483                BCOEFI = ZERO
00484 *
00485 *              Scale to avoid underflow
00486 *
00487                SCALE = ONE
00488                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
00489                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
00490      $               SMALL
00491                IF( LSA )
00492      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00493                IF( LSB )
00494      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
00495      $                    MIN( BNORM, BIG ) )
00496                IF( LSA .OR. LSB ) THEN
00497                   SCALE = MIN( SCALE, ONE /
00498      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
00499      $                    ABS( BCOEFR ) ) ) )
00500                   IF( LSA ) THEN
00501                      ACOEF = ASCALE*( SCALE*SBETA )
00502                   ELSE
00503                      ACOEF = SCALE*ACOEF
00504                   END IF
00505                   IF( LSB ) THEN
00506                      BCOEFR = BSCALE*( SCALE*SALFAR )
00507                   ELSE
00508                      BCOEFR = SCALE*BCOEFR
00509                   END IF
00510                END IF
00511                ACOEFA = ABS( ACOEF )
00512                BCOEFA = ABS( BCOEFR )
00513 *
00514 *              First component is 1
00515 *
00516                WORK( 2*N+JE ) = ONE
00517                XMAX = ONE
00518             ELSE
00519 *
00520 *              Complex eigenvalue
00521 *
00522                CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
00523      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
00524      $                     BCOEFI )
00525                BCOEFI = -BCOEFI
00526                IF( BCOEFI.EQ.ZERO ) THEN
00527                   INFO = JE
00528                   RETURN
00529                END IF
00530 *
00531 *              Scale to avoid over/underflow
00532 *
00533                ACOEFA = ABS( ACOEF )
00534                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00535                SCALE = ONE
00536                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
00537      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
00538                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
00539      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
00540                IF( SAFMIN*ACOEFA.GT.ASCALE )
00541      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
00542                IF( SAFMIN*BCOEFA.GT.BSCALE )
00543      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
00544                IF( SCALE.NE.ONE ) THEN
00545                   ACOEF = SCALE*ACOEF
00546                   ACOEFA = ABS( ACOEF )
00547                   BCOEFR = SCALE*BCOEFR
00548                   BCOEFI = SCALE*BCOEFI
00549                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00550                END IF
00551 *
00552 *              Compute first two components of eigenvector
00553 *
00554                TEMP = ACOEF*S( JE+1, JE )
00555                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
00556                TEMP2I = -BCOEFI*P( JE, JE )
00557                IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
00558                   WORK( 2*N+JE ) = ONE
00559                   WORK( 3*N+JE ) = ZERO
00560                   WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
00561                   WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
00562                ELSE
00563                   WORK( 2*N+JE+1 ) = ONE
00564                   WORK( 3*N+JE+1 ) = ZERO
00565                   TEMP = ACOEF*S( JE, JE+1 )
00566                   WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
00567      $                             S( JE+1, JE+1 ) ) / TEMP
00568                   WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
00569                END IF
00570                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
00571      $                ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
00572             END IF
00573 *
00574             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00575 *
00576 *                                           T
00577 *           Triangular solve of  (a A - b B)  y = 0
00578 *
00579 *                                   T
00580 *           (rowwise in  (a A - b B) , or columnwise in (a A - b B) )
00581 *
00582             IL2BY2 = .FALSE.
00583 *
00584             DO 160 J = JE + NW, N
00585                IF( IL2BY2 ) THEN
00586                   IL2BY2 = .FALSE.
00587                   GO TO 160
00588                END IF
00589 *
00590                NA = 1
00591                BDIAG( 1 ) = P( J, J )
00592                IF( J.LT.N ) THEN
00593                   IF( S( J+1, J ).NE.ZERO ) THEN
00594                      IL2BY2 = .TRUE.
00595                      BDIAG( 2 ) = P( J+1, J+1 )
00596                      NA = 2
00597                   END IF
00598                END IF
00599 *
00600 *              Check whether scaling is necessary for dot products
00601 *
00602                XSCALE = ONE / MAX( ONE, XMAX )
00603                TEMP = MAX( WORK( J ), WORK( N+J ),
00604      $                ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) )
00605                IF( IL2BY2 )
00606      $            TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ),
00607      $                   ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) )
00608                IF( TEMP.GT.BIGNUM*XSCALE ) THEN
00609                   DO 90 JW = 0, NW - 1
00610                      DO 80 JR = JE, J - 1
00611                         WORK( ( JW+2 )*N+JR ) = XSCALE*
00612      $                     WORK( ( JW+2 )*N+JR )
00613    80                CONTINUE
00614    90             CONTINUE
00615                   XMAX = XMAX*XSCALE
00616                END IF
00617 *
00618 *              Compute dot products
00619 *
00620 *                    j-1
00621 *              SUM = sum  conjg( a*S(k,j) - b*P(k,j) )*x(k)
00622 *                    k=je
00623 *
00624 *              To reduce the op count, this is done as
00625 *
00626 *              _        j-1                  _        j-1
00627 *              a*conjg( sum  S(k,j)*x(k) ) - b*conjg( sum  P(k,j)*x(k) )
00628 *                       k=je                          k=je
00629 *
00630 *              which may cause underflow problems if A or B are close
00631 *              to underflow.  (E.g., less than SMALL.)
00632 *
00633 *
00634 *              A series of compiler directives to defeat vectorization
00635 *              for the next loop
00636 *
00637 *$PL$ CMCHAR=' '
00638 CDIR$          NEXTSCALAR
00639 C$DIR          SCALAR
00640 CDIR$          NEXT SCALAR
00641 CVD$L          NOVECTOR
00642 CDEC$          NOVECTOR
00643 CVD$           NOVECTOR
00644 *VDIR          NOVECTOR
00645 *VOCL          LOOP,SCALAR
00646 CIBM           PREFER SCALAR
00647 *$PL$ CMCHAR='*'
00648 *
00649                DO 120 JW = 1, NW
00650 *
00651 *$PL$ CMCHAR=' '
00652 CDIR$             NEXTSCALAR
00653 C$DIR             SCALAR
00654 CDIR$             NEXT SCALAR
00655 CVD$L             NOVECTOR
00656 CDEC$             NOVECTOR
00657 CVD$              NOVECTOR
00658 *VDIR             NOVECTOR
00659 *VOCL             LOOP,SCALAR
00660 CIBM              PREFER SCALAR
00661 *$PL$ CMCHAR='*'
00662 *
00663                   DO 110 JA = 1, NA
00664                      SUMS( JA, JW ) = ZERO
00665                      SUMP( JA, JW ) = ZERO
00666 *
00667                      DO 100 JR = JE, J - 1
00668                         SUMS( JA, JW ) = SUMS( JA, JW ) +
00669      $                                   S( JR, J+JA-1 )*
00670      $                                   WORK( ( JW+1 )*N+JR )
00671                         SUMP( JA, JW ) = SUMP( JA, JW ) +
00672      $                                   P( JR, J+JA-1 )*
00673      $                                   WORK( ( JW+1 )*N+JR )
00674   100                CONTINUE
00675   110             CONTINUE
00676   120          CONTINUE
00677 *
00678 *$PL$ CMCHAR=' '
00679 CDIR$          NEXTSCALAR
00680 C$DIR          SCALAR
00681 CDIR$          NEXT SCALAR
00682 CVD$L          NOVECTOR
00683 CDEC$          NOVECTOR
00684 CVD$           NOVECTOR
00685 *VDIR          NOVECTOR
00686 *VOCL          LOOP,SCALAR
00687 CIBM           PREFER SCALAR
00688 *$PL$ CMCHAR='*'
00689 *
00690                DO 130 JA = 1, NA
00691                   IF( ILCPLX ) THEN
00692                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
00693      $                              BCOEFR*SUMP( JA, 1 ) -
00694      $                              BCOEFI*SUMP( JA, 2 )
00695                      SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
00696      $                              BCOEFR*SUMP( JA, 2 ) +
00697      $                              BCOEFI*SUMP( JA, 1 )
00698                   ELSE
00699                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
00700      $                              BCOEFR*SUMP( JA, 1 )
00701                   END IF
00702   130          CONTINUE
00703 *
00704 *                                  T
00705 *              Solve  ( a A - b B )  y = SUM(,)
00706 *              with scaling and perturbation of the denominator
00707 *
00708                CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
00709      $                      BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
00710      $                      BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
00711      $                      IINFO )
00712                IF( SCALE.LT.ONE ) THEN
00713                   DO 150 JW = 0, NW - 1
00714                      DO 140 JR = JE, J - 1
00715                         WORK( ( JW+2 )*N+JR ) = SCALE*
00716      $                     WORK( ( JW+2 )*N+JR )
00717   140                CONTINUE
00718   150             CONTINUE
00719                   XMAX = SCALE*XMAX
00720                END IF
00721                XMAX = MAX( XMAX, TEMP )
00722   160       CONTINUE
00723 *
00724 *           Copy eigenvector to VL, back transforming if
00725 *           HOWMNY='B'.
00726 *
00727             IEIG = IEIG + 1
00728             IF( ILBACK ) THEN
00729                DO 170 JW = 0, NW - 1
00730                   CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
00731      $                        WORK( ( JW+2 )*N+JE ), 1, ZERO,
00732      $                        WORK( ( JW+4 )*N+1 ), 1 )
00733   170          CONTINUE
00734                CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
00735      $                      LDVL )
00736                IBEG = 1
00737             ELSE
00738                CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
00739      $                      LDVL )
00740                IBEG = JE
00741             END IF
00742 *
00743 *           Scale eigenvector
00744 *
00745             XMAX = ZERO
00746             IF( ILCPLX ) THEN
00747                DO 180 J = IBEG, N
00748                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
00749      $                   ABS( VL( J, IEIG+1 ) ) )
00750   180          CONTINUE
00751             ELSE
00752                DO 190 J = IBEG, N
00753                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
00754   190          CONTINUE
00755             END IF
00756 *
00757             IF( XMAX.GT.SAFMIN ) THEN
00758                XSCALE = ONE / XMAX
00759 *
00760                DO 210 JW = 0, NW - 1
00761                   DO 200 JR = IBEG, N
00762                      VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
00763   200             CONTINUE
00764   210          CONTINUE
00765             END IF
00766             IEIG = IEIG + NW - 1
00767 *
00768   220    CONTINUE
00769       END IF
00770 *
00771 *     Right eigenvectors
00772 *
00773       IF( COMPR ) THEN
00774          IEIG = IM + 1
00775 *
00776 *        Main loop over eigenvalues
00777 *
00778          ILCPLX = .FALSE.
00779          DO 500 JE = N, 1, -1
00780 *
00781 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
00782 *           (b) this would be the second of a complex pair.
00783 *           Check for complex eigenvalue, so as to be sure of which
00784 *           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
00785 *           or SELECT(JE-1).
00786 *           If this is a complex pair, the 2-by-2 diagonal block
00787 *           corresponding to the eigenvalue is in rows/columns JE-1:JE
00788 *
00789             IF( ILCPLX ) THEN
00790                ILCPLX = .FALSE.
00791                GO TO 500
00792             END IF
00793             NW = 1
00794             IF( JE.GT.1 ) THEN
00795                IF( S( JE, JE-1 ).NE.ZERO ) THEN
00796                   ILCPLX = .TRUE.
00797                   NW = 2
00798                END IF
00799             END IF
00800             IF( ILALL ) THEN
00801                ILCOMP = .TRUE.
00802             ELSE IF( ILCPLX ) THEN
00803                ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
00804             ELSE
00805                ILCOMP = SELECT( JE )
00806             END IF
00807             IF( .NOT.ILCOMP )
00808      $         GO TO 500
00809 *
00810 *           Decide if (a) singular pencil, (b) real eigenvalue, or
00811 *           (c) complex eigenvalue.
00812 *
00813             IF( .NOT.ILCPLX ) THEN
00814                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
00815      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
00816 *
00817 *                 Singular matrix pencil -- unit eigenvector
00818 *
00819                   IEIG = IEIG - 1
00820                   DO 230 JR = 1, N
00821                      VR( JR, IEIG ) = ZERO
00822   230             CONTINUE
00823                   VR( IEIG, IEIG ) = ONE
00824                   GO TO 500
00825                END IF
00826             END IF
00827 *
00828 *           Clear vector
00829 *
00830             DO 250 JW = 0, NW - 1
00831                DO 240 JR = 1, N
00832                   WORK( ( JW+2 )*N+JR ) = ZERO
00833   240          CONTINUE
00834   250       CONTINUE
00835 *
00836 *           Compute coefficients in  ( a A - b B ) x = 0
00837 *              a  is  ACOEF
00838 *              b  is  BCOEFR + i*BCOEFI
00839 *
00840             IF( .NOT.ILCPLX ) THEN
00841 *
00842 *              Real eigenvalue
00843 *
00844                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
00845      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
00846                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
00847                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
00848                ACOEF = SBETA*ASCALE
00849                BCOEFR = SALFAR*BSCALE
00850                BCOEFI = ZERO
00851 *
00852 *              Scale to avoid underflow
00853 *
00854                SCALE = ONE
00855                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
00856                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
00857      $               SMALL
00858                IF( LSA )
00859      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00860                IF( LSB )
00861      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
00862      $                    MIN( BNORM, BIG ) )
00863                IF( LSA .OR. LSB ) THEN
00864                   SCALE = MIN( SCALE, ONE /
00865      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
00866      $                    ABS( BCOEFR ) ) ) )
00867                   IF( LSA ) THEN
00868                      ACOEF = ASCALE*( SCALE*SBETA )
00869                   ELSE
00870                      ACOEF = SCALE*ACOEF
00871                   END IF
00872                   IF( LSB ) THEN
00873                      BCOEFR = BSCALE*( SCALE*SALFAR )
00874                   ELSE
00875                      BCOEFR = SCALE*BCOEFR
00876                   END IF
00877                END IF
00878                ACOEFA = ABS( ACOEF )
00879                BCOEFA = ABS( BCOEFR )
00880 *
00881 *              First component is 1
00882 *
00883                WORK( 2*N+JE ) = ONE
00884                XMAX = ONE
00885 *
00886 *              Compute contribution from column JE of A and B to sum
00887 *              (See "Further Details", above.)
00888 *
00889                DO 260 JR = 1, JE - 1
00890                   WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
00891      $                             ACOEF*S( JR, JE )
00892   260          CONTINUE
00893             ELSE
00894 *
00895 *              Complex eigenvalue
00896 *
00897                CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
00898      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
00899      $                     BCOEFI )
00900                IF( BCOEFI.EQ.ZERO ) THEN
00901                   INFO = JE - 1
00902                   RETURN
00903                END IF
00904 *
00905 *              Scale to avoid over/underflow
00906 *
00907                ACOEFA = ABS( ACOEF )
00908                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00909                SCALE = ONE
00910                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
00911      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
00912                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
00913      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
00914                IF( SAFMIN*ACOEFA.GT.ASCALE )
00915      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
00916                IF( SAFMIN*BCOEFA.GT.BSCALE )
00917      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
00918                IF( SCALE.NE.ONE ) THEN
00919                   ACOEF = SCALE*ACOEF
00920                   ACOEFA = ABS( ACOEF )
00921                   BCOEFR = SCALE*BCOEFR
00922                   BCOEFI = SCALE*BCOEFI
00923                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00924                END IF
00925 *
00926 *              Compute first two components of eigenvector
00927 *              and contribution to sums
00928 *
00929                TEMP = ACOEF*S( JE, JE-1 )
00930                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
00931                TEMP2I = -BCOEFI*P( JE, JE )
00932                IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
00933                   WORK( 2*N+JE ) = ONE
00934                   WORK( 3*N+JE ) = ZERO
00935                   WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
00936                   WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
00937                ELSE
00938                   WORK( 2*N+JE-1 ) = ONE
00939                   WORK( 3*N+JE-1 ) = ZERO
00940                   TEMP = ACOEF*S( JE-1, JE )
00941                   WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
00942      $                             S( JE-1, JE-1 ) ) / TEMP
00943                   WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
00944                END IF
00945 *
00946                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
00947      $                ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
00948 *
00949 *              Compute contribution from columns JE and JE-1
00950 *              of A and B to the sums.
00951 *
00952                CREALA = ACOEF*WORK( 2*N+JE-1 )
00953                CIMAGA = ACOEF*WORK( 3*N+JE-1 )
00954                CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
00955      $                  BCOEFI*WORK( 3*N+JE-1 )
00956                CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
00957      $                  BCOEFR*WORK( 3*N+JE-1 )
00958                CRE2A = ACOEF*WORK( 2*N+JE )
00959                CIM2A = ACOEF*WORK( 3*N+JE )
00960                CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
00961                CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
00962                DO 270 JR = 1, JE - 2
00963                   WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
00964      $                             CREALB*P( JR, JE-1 ) -
00965      $                             CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
00966                   WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
00967      $                             CIMAGB*P( JR, JE-1 ) -
00968      $                             CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
00969   270          CONTINUE
00970             END IF
00971 *
00972             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00973 *
00974 *           Columnwise triangular solve of  (a A - b B)  x = 0
00975 *
00976             IL2BY2 = .FALSE.
00977             DO 370 J = JE - NW, 1, -1
00978 *
00979 *              If a 2-by-2 block, is in position j-1:j, wait until
00980 *              next iteration to process it (when it will be j:j+1)
00981 *
00982                IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
00983                   IF( S( J, J-1 ).NE.ZERO ) THEN
00984                      IL2BY2 = .TRUE.
00985                      GO TO 370
00986                   END IF
00987                END IF
00988                BDIAG( 1 ) = P( J, J )
00989                IF( IL2BY2 ) THEN
00990                   NA = 2
00991                   BDIAG( 2 ) = P( J+1, J+1 )
00992                ELSE
00993                   NA = 1
00994                END IF
00995 *
00996 *              Compute x(j) (and x(j+1), if 2-by-2 block)
00997 *
00998                CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
00999      $                      LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
01000      $                      N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
01001      $                      IINFO )
01002                IF( SCALE.LT.ONE ) THEN
01003 *
01004                   DO 290 JW = 0, NW - 1
01005                      DO 280 JR = 1, JE
01006                         WORK( ( JW+2 )*N+JR ) = SCALE*
01007      $                     WORK( ( JW+2 )*N+JR )
01008   280                CONTINUE
01009   290             CONTINUE
01010                END IF
01011                XMAX = MAX( SCALE*XMAX, TEMP )
01012 *
01013                DO 310 JW = 1, NW
01014                   DO 300 JA = 1, NA
01015                      WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
01016   300             CONTINUE
01017   310          CONTINUE
01018 *
01019 *              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
01020 *
01021                IF( J.GT.1 ) THEN
01022 *
01023 *                 Check whether scaling is necessary for sum.
01024 *
01025                   XSCALE = ONE / MAX( ONE, XMAX )
01026                   TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
01027                   IF( IL2BY2 )
01028      $               TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
01029      $                      WORK( N+J+1 ) )
01030                   TEMP = MAX( TEMP, ACOEFA, BCOEFA )
01031                   IF( TEMP.GT.BIGNUM*XSCALE ) THEN
01032 *
01033                      DO 330 JW = 0, NW - 1
01034                         DO 320 JR = 1, JE
01035                            WORK( ( JW+2 )*N+JR ) = XSCALE*
01036      $                        WORK( ( JW+2 )*N+JR )
01037   320                   CONTINUE
01038   330                CONTINUE
01039                      XMAX = XMAX*XSCALE
01040                   END IF
01041 *
01042 *                 Compute the contributions of the off-diagonals of
01043 *                 column j (and j+1, if 2-by-2 block) of A and B to the
01044 *                 sums.
01045 *
01046 *
01047                   DO 360 JA = 1, NA
01048                      IF( ILCPLX ) THEN
01049                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
01050                         CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
01051                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
01052      $                           BCOEFI*WORK( 3*N+J+JA-1 )
01053                         CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
01054      $                           BCOEFR*WORK( 3*N+J+JA-1 )
01055                         DO 340 JR = 1, J - 1
01056                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
01057      $                                      CREALA*S( JR, J+JA-1 ) +
01058      $                                      CREALB*P( JR, J+JA-1 )
01059                            WORK( 3*N+JR ) = WORK( 3*N+JR ) -
01060      $                                      CIMAGA*S( JR, J+JA-1 ) +
01061      $                                      CIMAGB*P( JR, J+JA-1 )
01062   340                   CONTINUE
01063                      ELSE
01064                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
01065                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
01066                         DO 350 JR = 1, J - 1
01067                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
01068      $                                      CREALA*S( JR, J+JA-1 ) +
01069      $                                      CREALB*P( JR, J+JA-1 )
01070   350                   CONTINUE
01071                      END IF
01072   360             CONTINUE
01073                END IF
01074 *
01075                IL2BY2 = .FALSE.
01076   370       CONTINUE
01077 *
01078 *           Copy eigenvector to VR, back transforming if
01079 *           HOWMNY='B'.
01080 *
01081             IEIG = IEIG - NW
01082             IF( ILBACK ) THEN
01083 *
01084                DO 410 JW = 0, NW - 1
01085                   DO 380 JR = 1, N
01086                      WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
01087      $                                       VR( JR, 1 )
01088   380             CONTINUE
01089 *
01090 *                 A series of compiler directives to defeat
01091 *                 vectorization for the next loop
01092 *
01093 *
01094                   DO 400 JC = 2, JE
01095                      DO 390 JR = 1, N
01096                         WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
01097      $                     WORK( ( JW+2 )*N+JC )*VR( JR, JC )
01098   390                CONTINUE
01099   400             CONTINUE
01100   410          CONTINUE
01101 *
01102                DO 430 JW = 0, NW - 1
01103                   DO 420 JR = 1, N
01104                      VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
01105   420             CONTINUE
01106   430          CONTINUE
01107 *
01108                IEND = N
01109             ELSE
01110                DO 450 JW = 0, NW - 1
01111                   DO 440 JR = 1, N
01112                      VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
01113   440             CONTINUE
01114   450          CONTINUE
01115 *
01116                IEND = JE
01117             END IF
01118 *
01119 *           Scale eigenvector
01120 *
01121             XMAX = ZERO
01122             IF( ILCPLX ) THEN
01123                DO 460 J = 1, IEND
01124                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
01125      $                   ABS( VR( J, IEIG+1 ) ) )
01126   460          CONTINUE
01127             ELSE
01128                DO 470 J = 1, IEND
01129                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
01130   470          CONTINUE
01131             END IF
01132 *
01133             IF( XMAX.GT.SAFMIN ) THEN
01134                XSCALE = ONE / XMAX
01135                DO 490 JW = 0, NW - 1
01136                   DO 480 JR = 1, IEND
01137                      VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
01138   480             CONTINUE
01139   490          CONTINUE
01140             END IF
01141   500    CONTINUE
01142       END IF
01143 *
01144       RETURN
01145 *
01146 *     End of DTGEVC
01147 *
01148       END