LAPACK 3.3.1
Linear Algebra PACKage

stgevc.f

Go to the documentation of this file.
00001       SUBROUTINE STGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
00002      $                   LDVL, VR, LDVR, MM, M, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
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       REAL               P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
00016      $                   VR( LDVR, * ), WORK( * )
00017 *     ..
00018 *
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  STGEVC 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 SGGHRD + SHGEQZ.
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) REAL array, dimension (LDS,N)
00078 *          The upper quasi-triangular matrix S from a generalized Schur
00079 *          factorization, as computed by SHGEQZ.
00080 *
00081 *  LDS     (input) INTEGER
00082 *          The leading dimension of array S.  LDS >= max(1,N).
00083 *
00084 *  P       (input) REAL array, dimension (LDP,N)
00085 *          The upper triangular matrix P from a generalized Schur
00086 *          factorization, as computed by SHGEQZ.
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) REAL 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 SHGEQZ).
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) REAL 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 SHGEQZ).
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) REAL 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       REAL               ZERO, ONE, SAFETY
00215       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0,
00216      $                   SAFETY = 1.0E+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       REAL               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       REAL               BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
00232      $                   SUMP( 2, 2 )
00233 *     ..
00234 *     .. External Functions ..
00235       LOGICAL            LSAME
00236       REAL               SLAMCH
00237       EXTERNAL           LSAME, SLAMCH
00238 *     ..
00239 *     .. External Subroutines ..
00240       EXTERNAL           SGEMV, SLABAD, SLACPY, SLAG2, SLALN2, 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( 'STGEVC', -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( 'STGEVC', -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 = SLAMCH( 'Safe minimum' )
00365       BIG = ONE / SAFMIN
00366       CALL SLABAD( SAFMIN, BIG )
00367       ULP = SLAMCH( 'Epsilon' )*SLAMCH( '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 SLAG2( 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                DO 120 JW = 1, NW
00635                   DO 110 JA = 1, NA
00636                      SUMS( JA, JW ) = ZERO
00637                      SUMP( JA, JW ) = ZERO
00638 *
00639                      DO 100 JR = JE, J - 1
00640                         SUMS( JA, JW ) = SUMS( JA, JW ) +
00641      $                                   S( JR, J+JA-1 )*
00642      $                                   WORK( ( JW+1 )*N+JR )
00643                         SUMP( JA, JW ) = SUMP( JA, JW ) +
00644      $                                   P( JR, J+JA-1 )*
00645      $                                   WORK( ( JW+1 )*N+JR )
00646   100                CONTINUE
00647   110             CONTINUE
00648   120          CONTINUE
00649 *
00650                DO 130 JA = 1, NA
00651                   IF( ILCPLX ) THEN
00652                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
00653      $                              BCOEFR*SUMP( JA, 1 ) -
00654      $                              BCOEFI*SUMP( JA, 2 )
00655                      SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
00656      $                              BCOEFR*SUMP( JA, 2 ) +
00657      $                              BCOEFI*SUMP( JA, 1 )
00658                   ELSE
00659                      SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
00660      $                              BCOEFR*SUMP( JA, 1 )
00661                   END IF
00662   130          CONTINUE
00663 *
00664 *                                  T
00665 *              Solve  ( a A - b B )  y = SUM(,)
00666 *              with scaling and perturbation of the denominator
00667 *
00668                CALL SLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
00669      $                      BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
00670      $                      BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
00671      $                      IINFO )
00672                IF( SCALE.LT.ONE ) THEN
00673                   DO 150 JW = 0, NW - 1
00674                      DO 140 JR = JE, J - 1
00675                         WORK( ( JW+2 )*N+JR ) = SCALE*
00676      $                     WORK( ( JW+2 )*N+JR )
00677   140                CONTINUE
00678   150             CONTINUE
00679                   XMAX = SCALE*XMAX
00680                END IF
00681                XMAX = MAX( XMAX, TEMP )
00682   160       CONTINUE
00683 *
00684 *           Copy eigenvector to VL, back transforming if
00685 *           HOWMNY='B'.
00686 *
00687             IEIG = IEIG + 1
00688             IF( ILBACK ) THEN
00689                DO 170 JW = 0, NW - 1
00690                   CALL SGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
00691      $                        WORK( ( JW+2 )*N+JE ), 1, ZERO,
00692      $                        WORK( ( JW+4 )*N+1 ), 1 )
00693   170          CONTINUE
00694                CALL SLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
00695      $                      LDVL )
00696                IBEG = 1
00697             ELSE
00698                CALL SLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
00699      $                      LDVL )
00700                IBEG = JE
00701             END IF
00702 *
00703 *           Scale eigenvector
00704 *
00705             XMAX = ZERO
00706             IF( ILCPLX ) THEN
00707                DO 180 J = IBEG, N
00708                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
00709      $                   ABS( VL( J, IEIG+1 ) ) )
00710   180          CONTINUE
00711             ELSE
00712                DO 190 J = IBEG, N
00713                   XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
00714   190          CONTINUE
00715             END IF
00716 *
00717             IF( XMAX.GT.SAFMIN ) THEN
00718                XSCALE = ONE / XMAX
00719 *
00720                DO 210 JW = 0, NW - 1
00721                   DO 200 JR = IBEG, N
00722                      VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
00723   200             CONTINUE
00724   210          CONTINUE
00725             END IF
00726             IEIG = IEIG + NW - 1
00727 *
00728   220    CONTINUE
00729       END IF
00730 *
00731 *     Right eigenvectors
00732 *
00733       IF( COMPR ) THEN
00734          IEIG = IM + 1
00735 *
00736 *        Main loop over eigenvalues
00737 *
00738          ILCPLX = .FALSE.
00739          DO 500 JE = N, 1, -1
00740 *
00741 *           Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
00742 *           (b) this would be the second of a complex pair.
00743 *           Check for complex eigenvalue, so as to be sure of which
00744 *           entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
00745 *           or SELECT(JE-1).
00746 *           If this is a complex pair, the 2-by-2 diagonal block
00747 *           corresponding to the eigenvalue is in rows/columns JE-1:JE
00748 *
00749             IF( ILCPLX ) THEN
00750                ILCPLX = .FALSE.
00751                GO TO 500
00752             END IF
00753             NW = 1
00754             IF( JE.GT.1 ) THEN
00755                IF( S( JE, JE-1 ).NE.ZERO ) THEN
00756                   ILCPLX = .TRUE.
00757                   NW = 2
00758                END IF
00759             END IF
00760             IF( ILALL ) THEN
00761                ILCOMP = .TRUE.
00762             ELSE IF( ILCPLX ) THEN
00763                ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
00764             ELSE
00765                ILCOMP = SELECT( JE )
00766             END IF
00767             IF( .NOT.ILCOMP )
00768      $         GO TO 500
00769 *
00770 *           Decide if (a) singular pencil, (b) real eigenvalue, or
00771 *           (c) complex eigenvalue.
00772 *
00773             IF( .NOT.ILCPLX ) THEN
00774                IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
00775      $             ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
00776 *
00777 *                 Singular matrix pencil -- unit eigenvector
00778 *
00779                   IEIG = IEIG - 1
00780                   DO 230 JR = 1, N
00781                      VR( JR, IEIG ) = ZERO
00782   230             CONTINUE
00783                   VR( IEIG, IEIG ) = ONE
00784                   GO TO 500
00785                END IF
00786             END IF
00787 *
00788 *           Clear vector
00789 *
00790             DO 250 JW = 0, NW - 1
00791                DO 240 JR = 1, N
00792                   WORK( ( JW+2 )*N+JR ) = ZERO
00793   240          CONTINUE
00794   250       CONTINUE
00795 *
00796 *           Compute coefficients in  ( a A - b B ) x = 0
00797 *              a  is  ACOEF
00798 *              b  is  BCOEFR + i*BCOEFI
00799 *
00800             IF( .NOT.ILCPLX ) THEN
00801 *
00802 *              Real eigenvalue
00803 *
00804                TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
00805      $                ABS( P( JE, JE ) )*BSCALE, SAFMIN )
00806                SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
00807                SBETA = ( TEMP*P( JE, JE ) )*BSCALE
00808                ACOEF = SBETA*ASCALE
00809                BCOEFR = SALFAR*BSCALE
00810                BCOEFI = ZERO
00811 *
00812 *              Scale to avoid underflow
00813 *
00814                SCALE = ONE
00815                LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
00816                LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
00817      $               SMALL
00818                IF( LSA )
00819      $            SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
00820                IF( LSB )
00821      $            SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
00822      $                    MIN( BNORM, BIG ) )
00823                IF( LSA .OR. LSB ) THEN
00824                   SCALE = MIN( SCALE, ONE /
00825      $                    ( SAFMIN*MAX( ONE, ABS( ACOEF ),
00826      $                    ABS( BCOEFR ) ) ) )
00827                   IF( LSA ) THEN
00828                      ACOEF = ASCALE*( SCALE*SBETA )
00829                   ELSE
00830                      ACOEF = SCALE*ACOEF
00831                   END IF
00832                   IF( LSB ) THEN
00833                      BCOEFR = BSCALE*( SCALE*SALFAR )
00834                   ELSE
00835                      BCOEFR = SCALE*BCOEFR
00836                   END IF
00837                END IF
00838                ACOEFA = ABS( ACOEF )
00839                BCOEFA = ABS( BCOEFR )
00840 *
00841 *              First component is 1
00842 *
00843                WORK( 2*N+JE ) = ONE
00844                XMAX = ONE
00845 *
00846 *              Compute contribution from column JE of A and B to sum
00847 *              (See "Further Details", above.)
00848 *
00849                DO 260 JR = 1, JE - 1
00850                   WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
00851      $                             ACOEF*S( JR, JE )
00852   260          CONTINUE
00853             ELSE
00854 *
00855 *              Complex eigenvalue
00856 *
00857                CALL SLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
00858      $                     SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
00859      $                     BCOEFI )
00860                IF( BCOEFI.EQ.ZERO ) THEN
00861                   INFO = JE - 1
00862                   RETURN
00863                END IF
00864 *
00865 *              Scale to avoid over/underflow
00866 *
00867                ACOEFA = ABS( ACOEF )
00868                BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00869                SCALE = ONE
00870                IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
00871      $            SCALE = ( SAFMIN / ULP ) / ACOEFA
00872                IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
00873      $            SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
00874                IF( SAFMIN*ACOEFA.GT.ASCALE )
00875      $            SCALE = ASCALE / ( SAFMIN*ACOEFA )
00876                IF( SAFMIN*BCOEFA.GT.BSCALE )
00877      $            SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
00878                IF( SCALE.NE.ONE ) THEN
00879                   ACOEF = SCALE*ACOEF
00880                   ACOEFA = ABS( ACOEF )
00881                   BCOEFR = SCALE*BCOEFR
00882                   BCOEFI = SCALE*BCOEFI
00883                   BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
00884                END IF
00885 *
00886 *              Compute first two components of eigenvector
00887 *              and contribution to sums
00888 *
00889                TEMP = ACOEF*S( JE, JE-1 )
00890                TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
00891                TEMP2I = -BCOEFI*P( JE, JE )
00892                IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
00893                   WORK( 2*N+JE ) = ONE
00894                   WORK( 3*N+JE ) = ZERO
00895                   WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
00896                   WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
00897                ELSE
00898                   WORK( 2*N+JE-1 ) = ONE
00899                   WORK( 3*N+JE-1 ) = ZERO
00900                   TEMP = ACOEF*S( JE-1, JE )
00901                   WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
00902      $                             S( JE-1, JE-1 ) ) / TEMP
00903                   WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
00904                END IF
00905 *
00906                XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
00907      $                ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
00908 *
00909 *              Compute contribution from columns JE and JE-1
00910 *              of A and B to the sums.
00911 *
00912                CREALA = ACOEF*WORK( 2*N+JE-1 )
00913                CIMAGA = ACOEF*WORK( 3*N+JE-1 )
00914                CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
00915      $                  BCOEFI*WORK( 3*N+JE-1 )
00916                CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
00917      $                  BCOEFR*WORK( 3*N+JE-1 )
00918                CRE2A = ACOEF*WORK( 2*N+JE )
00919                CIM2A = ACOEF*WORK( 3*N+JE )
00920                CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
00921                CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
00922                DO 270 JR = 1, JE - 2
00923                   WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
00924      $                             CREALB*P( JR, JE-1 ) -
00925      $                             CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
00926                   WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
00927      $                             CIMAGB*P( JR, JE-1 ) -
00928      $                             CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
00929   270          CONTINUE
00930             END IF
00931 *
00932             DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
00933 *
00934 *           Columnwise triangular solve of  (a A - b B)  x = 0
00935 *
00936             IL2BY2 = .FALSE.
00937             DO 370 J = JE - NW, 1, -1
00938 *
00939 *              If a 2-by-2 block, is in position j-1:j, wait until
00940 *              next iteration to process it (when it will be j:j+1)
00941 *
00942                IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
00943                   IF( S( J, J-1 ).NE.ZERO ) THEN
00944                      IL2BY2 = .TRUE.
00945                      GO TO 370
00946                   END IF
00947                END IF
00948                BDIAG( 1 ) = P( J, J )
00949                IF( IL2BY2 ) THEN
00950                   NA = 2
00951                   BDIAG( 2 ) = P( J+1, J+1 )
00952                ELSE
00953                   NA = 1
00954                END IF
00955 *
00956 *              Compute x(j) (and x(j+1), if 2-by-2 block)
00957 *
00958                CALL SLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
00959      $                      LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
00960      $                      N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
00961      $                      IINFO )
00962                IF( SCALE.LT.ONE ) THEN
00963 *
00964                   DO 290 JW = 0, NW - 1
00965                      DO 280 JR = 1, JE
00966                         WORK( ( JW+2 )*N+JR ) = SCALE*
00967      $                     WORK( ( JW+2 )*N+JR )
00968   280                CONTINUE
00969   290             CONTINUE
00970                END IF
00971                XMAX = MAX( SCALE*XMAX, TEMP )
00972 *
00973                DO 310 JW = 1, NW
00974                   DO 300 JA = 1, NA
00975                      WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
00976   300             CONTINUE
00977   310          CONTINUE
00978 *
00979 *              w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
00980 *
00981                IF( J.GT.1 ) THEN
00982 *
00983 *                 Check whether scaling is necessary for sum.
00984 *
00985                   XSCALE = ONE / MAX( ONE, XMAX )
00986                   TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
00987                   IF( IL2BY2 )
00988      $               TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
00989      $                      WORK( N+J+1 ) )
00990                   TEMP = MAX( TEMP, ACOEFA, BCOEFA )
00991                   IF( TEMP.GT.BIGNUM*XSCALE ) THEN
00992 *
00993                      DO 330 JW = 0, NW - 1
00994                         DO 320 JR = 1, JE
00995                            WORK( ( JW+2 )*N+JR ) = XSCALE*
00996      $                        WORK( ( JW+2 )*N+JR )
00997   320                   CONTINUE
00998   330                CONTINUE
00999                      XMAX = XMAX*XSCALE
01000                   END IF
01001 *
01002 *                 Compute the contributions of the off-diagonals of
01003 *                 column j (and j+1, if 2-by-2 block) of A and B to the
01004 *                 sums.
01005 *
01006 *
01007                   DO 360 JA = 1, NA
01008                      IF( ILCPLX ) THEN
01009                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
01010                         CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
01011                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
01012      $                           BCOEFI*WORK( 3*N+J+JA-1 )
01013                         CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
01014      $                           BCOEFR*WORK( 3*N+J+JA-1 )
01015                         DO 340 JR = 1, J - 1
01016                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
01017      $                                      CREALA*S( JR, J+JA-1 ) +
01018      $                                      CREALB*P( JR, J+JA-1 )
01019                            WORK( 3*N+JR ) = WORK( 3*N+JR ) -
01020      $                                      CIMAGA*S( JR, J+JA-1 ) +
01021      $                                      CIMAGB*P( JR, J+JA-1 )
01022   340                   CONTINUE
01023                      ELSE
01024                         CREALA = ACOEF*WORK( 2*N+J+JA-1 )
01025                         CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
01026                         DO 350 JR = 1, J - 1
01027                            WORK( 2*N+JR ) = WORK( 2*N+JR ) -
01028      $                                      CREALA*S( JR, J+JA-1 ) +
01029      $                                      CREALB*P( JR, J+JA-1 )
01030   350                   CONTINUE
01031                      END IF
01032   360             CONTINUE
01033                END IF
01034 *
01035                IL2BY2 = .FALSE.
01036   370       CONTINUE
01037 *
01038 *           Copy eigenvector to VR, back transforming if
01039 *           HOWMNY='B'.
01040 *
01041             IEIG = IEIG - NW
01042             IF( ILBACK ) THEN
01043 *
01044                DO 410 JW = 0, NW - 1
01045                   DO 380 JR = 1, N
01046                      WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
01047      $                                       VR( JR, 1 )
01048   380             CONTINUE
01049 *
01050 *                 A series of compiler directives to defeat
01051 *                 vectorization for the next loop
01052 *
01053 *
01054                   DO 400 JC = 2, JE
01055                      DO 390 JR = 1, N
01056                         WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
01057      $                     WORK( ( JW+2 )*N+JC )*VR( JR, JC )
01058   390                CONTINUE
01059   400             CONTINUE
01060   410          CONTINUE
01061 *
01062                DO 430 JW = 0, NW - 1
01063                   DO 420 JR = 1, N
01064                      VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
01065   420             CONTINUE
01066   430          CONTINUE
01067 *
01068                IEND = N
01069             ELSE
01070                DO 450 JW = 0, NW - 1
01071                   DO 440 JR = 1, N
01072                      VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
01073   440             CONTINUE
01074   450          CONTINUE
01075 *
01076                IEND = JE
01077             END IF
01078 *
01079 *           Scale eigenvector
01080 *
01081             XMAX = ZERO
01082             IF( ILCPLX ) THEN
01083                DO 460 J = 1, IEND
01084                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
01085      $                   ABS( VR( J, IEIG+1 ) ) )
01086   460          CONTINUE
01087             ELSE
01088                DO 470 J = 1, IEND
01089                   XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
01090   470          CONTINUE
01091             END IF
01092 *
01093             IF( XMAX.GT.SAFMIN ) THEN
01094                XSCALE = ONE / XMAX
01095                DO 490 JW = 0, NW - 1
01096                   DO 480 JR = 1, IEND
01097                      VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
01098   480             CONTINUE
01099   490          CONTINUE
01100             END IF
01101   500    CONTINUE
01102       END IF
01103 *
01104       RETURN
01105 *
01106 *     End of STGEVC
01107 *
01108       END
 All Files Functions