LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
strevc.f
Go to the documentation of this file.
00001 *> \brief \b STREVC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STREVC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strevc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strevc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strevc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00022 *                          LDVR, MM, M, WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          HOWMNY, SIDE
00026 *       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       LOGICAL            SELECT( * )
00030 *       REAL               T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00031 *      $                   WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> STREVC computes some or all of the right and/or left eigenvectors of
00041 *> a real upper quasi-triangular matrix T.
00042 *> Matrices of this type are produced by the Schur factorization of
00043 *> a real general matrix:  A = Q*T*Q**T, as computed by SHSEQR.
00044 *> 
00045 *> The right eigenvector x and the left eigenvector y of T corresponding
00046 *> to an eigenvalue w are defined by:
00047 *> 
00048 *>    T*x = w*x,     (y**T)*T = w*(y**T)
00049 *> 
00050 *> where y**T denotes the transpose of y.
00051 *> The eigenvalues are not input to this routine, but are read directly
00052 *> from the diagonal blocks of T.
00053 *> 
00054 *> This routine returns the matrices X and/or Y of right and left
00055 *> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
00056 *> input matrix.  If Q is the orthogonal factor that reduces a matrix
00057 *> A to Schur form T, then Q*X and Q*Y are the matrices of right and
00058 *> left eigenvectors of A.
00059 *> \endverbatim
00060 *
00061 *  Arguments:
00062 *  ==========
00063 *
00064 *> \param[in] SIDE
00065 *> \verbatim
00066 *>          SIDE is CHARACTER*1
00067 *>          = 'R':  compute right eigenvectors only;
00068 *>          = 'L':  compute left eigenvectors only;
00069 *>          = 'B':  compute both right and left eigenvectors.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] HOWMNY
00073 *> \verbatim
00074 *>          HOWMNY is CHARACTER*1
00075 *>          = 'A':  compute all right and/or left eigenvectors;
00076 *>          = 'B':  compute all right and/or left eigenvectors,
00077 *>                  backtransformed by the matrices in VR and/or VL;
00078 *>          = 'S':  compute selected right and/or left eigenvectors,
00079 *>                  as indicated by the logical array SELECT.
00080 *> \endverbatim
00081 *>
00082 *> \param[in,out] SELECT
00083 *> \verbatim
00084 *>          SELECT is LOGICAL array, dimension (N)
00085 *>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
00086 *>          computed.
00087 *>          If w(j) is a real eigenvalue, the corresponding real
00088 *>          eigenvector is computed if SELECT(j) is .TRUE..
00089 *>          If w(j) and w(j+1) are the real and imaginary parts of a
00090 *>          complex eigenvalue, the corresponding complex eigenvector is
00091 *>          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
00092 *>          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
00093 *>          .FALSE..
00094 *>          Not referenced if HOWMNY = 'A' or 'B'.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] N
00098 *> \verbatim
00099 *>          N is INTEGER
00100 *>          The order of the matrix T. N >= 0.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] T
00104 *> \verbatim
00105 *>          T is REAL array, dimension (LDT,N)
00106 *>          The upper quasi-triangular matrix T in Schur canonical form.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LDT
00110 *> \verbatim
00111 *>          LDT is INTEGER
00112 *>          The leading dimension of the array T. LDT >= max(1,N).
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] VL
00116 *> \verbatim
00117 *>          VL is REAL array, dimension (LDVL,MM)
00118 *>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00119 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
00120 *>          of Schur vectors returned by SHSEQR).
00121 *>          On exit, if SIDE = 'L' or 'B', VL contains:
00122 *>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
00123 *>          if HOWMNY = 'B', the matrix Q*Y;
00124 *>          if HOWMNY = 'S', the left eigenvectors of T specified by
00125 *>                           SELECT, stored consecutively in the columns
00126 *>                           of VL, in the same order as their
00127 *>                           eigenvalues.
00128 *>          A complex eigenvector corresponding to a complex eigenvalue
00129 *>          is stored in two consecutive columns, the first holding the
00130 *>          real part, and the second the imaginary part.
00131 *>          Not referenced if SIDE = 'R'.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] LDVL
00135 *> \verbatim
00136 *>          LDVL is INTEGER
00137 *>          The leading dimension of the array VL.  LDVL >= 1, and if
00138 *>          SIDE = 'L' or 'B', LDVL >= N.
00139 *> \endverbatim
00140 *>
00141 *> \param[in,out] VR
00142 *> \verbatim
00143 *>          VR is REAL array, dimension (LDVR,MM)
00144 *>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00145 *>          contain an N-by-N matrix Q (usually the orthogonal matrix Q
00146 *>          of Schur vectors returned by SHSEQR).
00147 *>          On exit, if SIDE = 'R' or 'B', VR contains:
00148 *>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
00149 *>          if HOWMNY = 'B', the matrix Q*X;
00150 *>          if HOWMNY = 'S', the right eigenvectors of T specified by
00151 *>                           SELECT, stored consecutively in the columns
00152 *>                           of VR, in the same order as their
00153 *>                           eigenvalues.
00154 *>          A complex eigenvector corresponding to a complex eigenvalue
00155 *>          is stored in two consecutive columns, the first holding the
00156 *>          real part and the second the imaginary part.
00157 *>          Not referenced if SIDE = 'L'.
00158 *> \endverbatim
00159 *>
00160 *> \param[in] LDVR
00161 *> \verbatim
00162 *>          LDVR is INTEGER
00163 *>          The leading dimension of the array VR.  LDVR >= 1, and if
00164 *>          SIDE = 'R' or 'B', LDVR >= N.
00165 *> \endverbatim
00166 *>
00167 *> \param[in] MM
00168 *> \verbatim
00169 *>          MM is INTEGER
00170 *>          The number of columns in the arrays VL and/or VR. MM >= M.
00171 *> \endverbatim
00172 *>
00173 *> \param[out] M
00174 *> \verbatim
00175 *>          M is INTEGER
00176 *>          The number of columns in the arrays VL and/or VR actually
00177 *>          used to store the eigenvectors.
00178 *>          If HOWMNY = 'A' or 'B', M is set to N.
00179 *>          Each selected real eigenvector occupies one column and each
00180 *>          selected complex eigenvector occupies two columns.
00181 *> \endverbatim
00182 *>
00183 *> \param[out] WORK
00184 *> \verbatim
00185 *>          WORK is REAL array, dimension (3*N)
00186 *> \endverbatim
00187 *>
00188 *> \param[out] INFO
00189 *> \verbatim
00190 *>          INFO is INTEGER
00191 *>          = 0:  successful exit
00192 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00193 *> \endverbatim
00194 *
00195 *  Authors:
00196 *  ========
00197 *
00198 *> \author Univ. of Tennessee 
00199 *> \author Univ. of California Berkeley 
00200 *> \author Univ. of Colorado Denver 
00201 *> \author NAG Ltd. 
00202 *
00203 *> \date November 2011
00204 *
00205 *> \ingroup realOTHERcomputational
00206 *
00207 *> \par Further Details:
00208 *  =====================
00209 *>
00210 *> \verbatim
00211 *>
00212 *>  The algorithm used in this program is basically backward (forward)
00213 *>  substitution, with scaling to make the the code robust against
00214 *>  possible overflow.
00215 *>
00216 *>  Each eigenvector is normalized so that the element of largest
00217 *>  magnitude has magnitude 1; here the magnitude of a complex number
00218 *>  (x,y) is taken to be |x| + |y|.
00219 *> \endverbatim
00220 *>
00221 *  =====================================================================
00222       SUBROUTINE STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00223      $                   LDVR, MM, M, WORK, INFO )
00224 *
00225 *  -- LAPACK computational routine (version 3.4.0) --
00226 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00227 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00228 *     November 2011
00229 *
00230 *     .. Scalar Arguments ..
00231       CHARACTER          HOWMNY, SIDE
00232       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
00233 *     ..
00234 *     .. Array Arguments ..
00235       LOGICAL            SELECT( * )
00236       REAL               T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00237      $                   WORK( * )
00238 *     ..
00239 *
00240 *  =====================================================================
00241 *
00242 *     .. Parameters ..
00243       REAL               ZERO, ONE
00244       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00245 *     ..
00246 *     .. Local Scalars ..
00247       LOGICAL            ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
00248       INTEGER            I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
00249       REAL               BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
00250      $                   SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
00251      $                   XNORM
00252 *     ..
00253 *     .. External Functions ..
00254       LOGICAL            LSAME
00255       INTEGER            ISAMAX
00256       REAL               SDOT, SLAMCH
00257       EXTERNAL           LSAME, ISAMAX, SDOT, SLAMCH
00258 *     ..
00259 *     .. External Subroutines ..
00260       EXTERNAL           SAXPY, SCOPY, SGEMV, SLABAD, SLALN2, SSCAL,
00261      $                   XERBLA
00262 *     ..
00263 *     .. Intrinsic Functions ..
00264       INTRINSIC          ABS, MAX, SQRT
00265 *     ..
00266 *     .. Local Arrays ..
00267       REAL               X( 2, 2 )
00268 *     ..
00269 *     .. Executable Statements ..
00270 *
00271 *     Decode and test the input parameters
00272 *
00273       BOTHV = LSAME( SIDE, 'B' )
00274       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00275       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00276 *
00277       ALLV = LSAME( HOWMNY, 'A' )
00278       OVER = LSAME( HOWMNY, 'B' )
00279       SOMEV = LSAME( HOWMNY, 'S' )
00280 *
00281       INFO = 0
00282       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00283          INFO = -1
00284       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
00285          INFO = -2
00286       ELSE IF( N.LT.0 ) THEN
00287          INFO = -4
00288       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00289          INFO = -6
00290       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00291          INFO = -8
00292       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00293          INFO = -10
00294       ELSE
00295 *
00296 *        Set M to the number of columns required to store the selected
00297 *        eigenvectors, standardize the array SELECT if necessary, and
00298 *        test MM.
00299 *
00300          IF( SOMEV ) THEN
00301             M = 0
00302             PAIR = .FALSE.
00303             DO 10 J = 1, N
00304                IF( PAIR ) THEN
00305                   PAIR = .FALSE.
00306                   SELECT( J ) = .FALSE.
00307                ELSE
00308                   IF( J.LT.N ) THEN
00309                      IF( T( J+1, J ).EQ.ZERO ) THEN
00310                         IF( SELECT( J ) )
00311      $                     M = M + 1
00312                      ELSE
00313                         PAIR = .TRUE.
00314                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
00315                            SELECT( J ) = .TRUE.
00316                            M = M + 2
00317                         END IF
00318                      END IF
00319                   ELSE
00320                      IF( SELECT( N ) )
00321      $                  M = M + 1
00322                   END IF
00323                END IF
00324    10       CONTINUE
00325          ELSE
00326             M = N
00327          END IF
00328 *
00329          IF( MM.LT.M ) THEN
00330             INFO = -11
00331          END IF
00332       END IF
00333       IF( INFO.NE.0 ) THEN
00334          CALL XERBLA( 'STREVC', -INFO )
00335          RETURN
00336       END IF
00337 *
00338 *     Quick return if possible.
00339 *
00340       IF( N.EQ.0 )
00341      $   RETURN
00342 *
00343 *     Set the constants to control overflow.
00344 *
00345       UNFL = SLAMCH( 'Safe minimum' )
00346       OVFL = ONE / UNFL
00347       CALL SLABAD( UNFL, OVFL )
00348       ULP = SLAMCH( 'Precision' )
00349       SMLNUM = UNFL*( N / ULP )
00350       BIGNUM = ( ONE-ULP ) / SMLNUM
00351 *
00352 *     Compute 1-norm of each column of strictly upper triangular
00353 *     part of T to control overflow in triangular solver.
00354 *
00355       WORK( 1 ) = ZERO
00356       DO 30 J = 2, N
00357          WORK( J ) = ZERO
00358          DO 20 I = 1, J - 1
00359             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
00360    20    CONTINUE
00361    30 CONTINUE
00362 *
00363 *     Index IP is used to specify the real or complex eigenvalue:
00364 *       IP = 0, real eigenvalue,
00365 *            1, first of conjugate complex pair: (wr,wi)
00366 *           -1, second of conjugate complex pair: (wr,wi)
00367 *
00368       N2 = 2*N
00369 *
00370       IF( RIGHTV ) THEN
00371 *
00372 *        Compute right eigenvectors.
00373 *
00374          IP = 0
00375          IS = M
00376          DO 140 KI = N, 1, -1
00377 *
00378             IF( IP.EQ.1 )
00379      $         GO TO 130
00380             IF( KI.EQ.1 )
00381      $         GO TO 40
00382             IF( T( KI, KI-1 ).EQ.ZERO )
00383      $         GO TO 40
00384             IP = -1
00385 *
00386    40       CONTINUE
00387             IF( SOMEV ) THEN
00388                IF( IP.EQ.0 ) THEN
00389                   IF( .NOT.SELECT( KI ) )
00390      $               GO TO 130
00391                ELSE
00392                   IF( .NOT.SELECT( KI-1 ) )
00393      $               GO TO 130
00394                END IF
00395             END IF
00396 *
00397 *           Compute the KI-th eigenvalue (WR,WI).
00398 *
00399             WR = T( KI, KI )
00400             WI = ZERO
00401             IF( IP.NE.0 )
00402      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
00403      $              SQRT( ABS( T( KI-1, KI ) ) )
00404             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00405 *
00406             IF( IP.EQ.0 ) THEN
00407 *
00408 *              Real right eigenvector
00409 *
00410                WORK( KI+N ) = ONE
00411 *
00412 *              Form right-hand side
00413 *
00414                DO 50 K = 1, KI - 1
00415                   WORK( K+N ) = -T( K, KI )
00416    50          CONTINUE
00417 *
00418 *              Solve the upper quasi-triangular system:
00419 *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
00420 *
00421                JNXT = KI - 1
00422                DO 60 J = KI - 1, 1, -1
00423                   IF( J.GT.JNXT )
00424      $               GO TO 60
00425                   J1 = J
00426                   J2 = J
00427                   JNXT = J - 1
00428                   IF( J.GT.1 ) THEN
00429                      IF( T( J, J-1 ).NE.ZERO ) THEN
00430                         J1 = J - 1
00431                         JNXT = J - 2
00432                      END IF
00433                   END IF
00434 *
00435                   IF( J1.EQ.J2 ) THEN
00436 *
00437 *                    1-by-1 diagonal block
00438 *
00439                      CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00440      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00441      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00442 *
00443 *                    Scale X(1,1) to avoid overflow when updating
00444 *                    the right-hand side.
00445 *
00446                      IF( XNORM.GT.ONE ) THEN
00447                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00448                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00449                            SCALE = SCALE / XNORM
00450                         END IF
00451                      END IF
00452 *
00453 *                    Scale if necessary
00454 *
00455                      IF( SCALE.NE.ONE )
00456      $                  CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00457                      WORK( J+N ) = X( 1, 1 )
00458 *
00459 *                    Update right-hand side
00460 *
00461                      CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00462      $                           WORK( 1+N ), 1 )
00463 *
00464                   ELSE
00465 *
00466 *                    2-by-2 diagonal block
00467 *
00468                      CALL SLALN2( .FALSE., 2, 1, SMIN, ONE,
00469      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00470      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
00471      $                            SCALE, XNORM, IERR )
00472 *
00473 *                    Scale X(1,1) and X(2,1) to avoid overflow when
00474 *                    updating the right-hand side.
00475 *
00476                      IF( XNORM.GT.ONE ) THEN
00477                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00478                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00479                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00480                            X( 2, 1 ) = X( 2, 1 ) / XNORM
00481                            SCALE = SCALE / XNORM
00482                         END IF
00483                      END IF
00484 *
00485 *                    Scale if necessary
00486 *
00487                      IF( SCALE.NE.ONE )
00488      $                  CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00489                      WORK( J-1+N ) = X( 1, 1 )
00490                      WORK( J+N ) = X( 2, 1 )
00491 *
00492 *                    Update right-hand side
00493 *
00494                      CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00495      $                           WORK( 1+N ), 1 )
00496                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00497      $                           WORK( 1+N ), 1 )
00498                   END IF
00499    60          CONTINUE
00500 *
00501 *              Copy the vector x or Q*x to VR and normalize.
00502 *
00503                IF( .NOT.OVER ) THEN
00504                   CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
00505 *
00506                   II = ISAMAX( KI, VR( 1, IS ), 1 )
00507                   REMAX = ONE / ABS( VR( II, IS ) )
00508                   CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
00509 *
00510                   DO 70 K = KI + 1, N
00511                      VR( K, IS ) = ZERO
00512    70             CONTINUE
00513                ELSE
00514                   IF( KI.GT.1 )
00515      $               CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
00516      $                           WORK( 1+N ), 1, WORK( KI+N ),
00517      $                           VR( 1, KI ), 1 )
00518 *
00519                   II = ISAMAX( N, VR( 1, KI ), 1 )
00520                   REMAX = ONE / ABS( VR( II, KI ) )
00521                   CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
00522                END IF
00523 *
00524             ELSE
00525 *
00526 *              Complex right eigenvector.
00527 *
00528 *              Initial solve
00529 *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
00530 *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
00531 *
00532                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
00533                   WORK( KI-1+N ) = ONE
00534                   WORK( KI+N2 ) = WI / T( KI-1, KI )
00535                ELSE
00536                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
00537                   WORK( KI+N2 ) = ONE
00538                END IF
00539                WORK( KI+N ) = ZERO
00540                WORK( KI-1+N2 ) = ZERO
00541 *
00542 *              Form right-hand side
00543 *
00544                DO 80 K = 1, KI - 2
00545                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
00546                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
00547    80          CONTINUE
00548 *
00549 *              Solve upper quasi-triangular system:
00550 *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
00551 *
00552                JNXT = KI - 2
00553                DO 90 J = KI - 2, 1, -1
00554                   IF( J.GT.JNXT )
00555      $               GO TO 90
00556                   J1 = J
00557                   J2 = J
00558                   JNXT = J - 1
00559                   IF( J.GT.1 ) THEN
00560                      IF( T( J, J-1 ).NE.ZERO ) THEN
00561                         J1 = J - 1
00562                         JNXT = J - 2
00563                      END IF
00564                   END IF
00565 *
00566                   IF( J1.EQ.J2 ) THEN
00567 *
00568 *                    1-by-1 diagonal block
00569 *
00570                      CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00571      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
00572      $                            X, 2, SCALE, XNORM, IERR )
00573 *
00574 *                    Scale X(1,1) and X(1,2) to avoid overflow when
00575 *                    updating the right-hand side.
00576 *
00577                      IF( XNORM.GT.ONE ) THEN
00578                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00579                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00580                            X( 1, 2 ) = X( 1, 2 ) / XNORM
00581                            SCALE = SCALE / XNORM
00582                         END IF
00583                      END IF
00584 *
00585 *                    Scale if necessary
00586 *
00587                      IF( SCALE.NE.ONE ) THEN
00588                         CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00589                         CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00590                      END IF
00591                      WORK( J+N ) = X( 1, 1 )
00592                      WORK( J+N2 ) = X( 1, 2 )
00593 *
00594 *                    Update the right-hand side
00595 *
00596                      CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00597      $                           WORK( 1+N ), 1 )
00598                      CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
00599      $                           WORK( 1+N2 ), 1 )
00600 *
00601                   ELSE
00602 *
00603 *                    2-by-2 diagonal block
00604 *
00605                      CALL SLALN2( .FALSE., 2, 2, SMIN, ONE,
00606      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00607      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
00608      $                            XNORM, IERR )
00609 *
00610 *                    Scale X to avoid overflow when updating
00611 *                    the right-hand side.
00612 *
00613                      IF( XNORM.GT.ONE ) THEN
00614                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00615                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00616                            REC = ONE / XNORM
00617                            X( 1, 1 ) = X( 1, 1 )*REC
00618                            X( 1, 2 ) = X( 1, 2 )*REC
00619                            X( 2, 1 ) = X( 2, 1 )*REC
00620                            X( 2, 2 ) = X( 2, 2 )*REC
00621                            SCALE = SCALE*REC
00622                         END IF
00623                      END IF
00624 *
00625 *                    Scale if necessary
00626 *
00627                      IF( SCALE.NE.ONE ) THEN
00628                         CALL SSCAL( KI, SCALE, WORK( 1+N ), 1 )
00629                         CALL SSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00630                      END IF
00631                      WORK( J-1+N ) = X( 1, 1 )
00632                      WORK( J+N ) = X( 2, 1 )
00633                      WORK( J-1+N2 ) = X( 1, 2 )
00634                      WORK( J+N2 ) = X( 2, 2 )
00635 *
00636 *                    Update the right-hand side
00637 *
00638                      CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00639      $                           WORK( 1+N ), 1 )
00640                      CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00641      $                           WORK( 1+N ), 1 )
00642                      CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
00643      $                           WORK( 1+N2 ), 1 )
00644                      CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
00645      $                           WORK( 1+N2 ), 1 )
00646                   END IF
00647    90          CONTINUE
00648 *
00649 *              Copy the vector x or Q*x to VR and normalize.
00650 *
00651                IF( .NOT.OVER ) THEN
00652                   CALL SCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
00653                   CALL SCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
00654 *
00655                   EMAX = ZERO
00656                   DO 100 K = 1, KI
00657                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
00658      $                      ABS( VR( K, IS ) ) )
00659   100             CONTINUE
00660 *
00661                   REMAX = ONE / EMAX
00662                   CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
00663                   CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
00664 *
00665                   DO 110 K = KI + 1, N
00666                      VR( K, IS-1 ) = ZERO
00667                      VR( K, IS ) = ZERO
00668   110             CONTINUE
00669 *
00670                ELSE
00671 *
00672                   IF( KI.GT.2 ) THEN
00673                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00674      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
00675      $                           VR( 1, KI-1 ), 1 )
00676                      CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00677      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
00678      $                           VR( 1, KI ), 1 )
00679                   ELSE
00680                      CALL SSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
00681                      CALL SSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
00682                   END IF
00683 *
00684                   EMAX = ZERO
00685                   DO 120 K = 1, N
00686                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
00687      $                      ABS( VR( K, KI ) ) )
00688   120             CONTINUE
00689                   REMAX = ONE / EMAX
00690                   CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
00691                   CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
00692                END IF
00693             END IF
00694 *
00695             IS = IS - 1
00696             IF( IP.NE.0 )
00697      $         IS = IS - 1
00698   130       CONTINUE
00699             IF( IP.EQ.1 )
00700      $         IP = 0
00701             IF( IP.EQ.-1 )
00702      $         IP = 1
00703   140    CONTINUE
00704       END IF
00705 *
00706       IF( LEFTV ) THEN
00707 *
00708 *        Compute left eigenvectors.
00709 *
00710          IP = 0
00711          IS = 1
00712          DO 260 KI = 1, N
00713 *
00714             IF( IP.EQ.-1 )
00715      $         GO TO 250
00716             IF( KI.EQ.N )
00717      $         GO TO 150
00718             IF( T( KI+1, KI ).EQ.ZERO )
00719      $         GO TO 150
00720             IP = 1
00721 *
00722   150       CONTINUE
00723             IF( SOMEV ) THEN
00724                IF( .NOT.SELECT( KI ) )
00725      $            GO TO 250
00726             END IF
00727 *
00728 *           Compute the KI-th eigenvalue (WR,WI).
00729 *
00730             WR = T( KI, KI )
00731             WI = ZERO
00732             IF( IP.NE.0 )
00733      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
00734      $              SQRT( ABS( T( KI+1, KI ) ) )
00735             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00736 *
00737             IF( IP.EQ.0 ) THEN
00738 *
00739 *              Real left eigenvector.
00740 *
00741                WORK( KI+N ) = ONE
00742 *
00743 *              Form right-hand side
00744 *
00745                DO 160 K = KI + 1, N
00746                   WORK( K+N ) = -T( KI, K )
00747   160          CONTINUE
00748 *
00749 *              Solve the quasi-triangular system:
00750 *                 (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
00751 *
00752                VMAX = ONE
00753                VCRIT = BIGNUM
00754 *
00755                JNXT = KI + 1
00756                DO 170 J = KI + 1, N
00757                   IF( J.LT.JNXT )
00758      $               GO TO 170
00759                   J1 = J
00760                   J2 = J
00761                   JNXT = J + 1
00762                   IF( J.LT.N ) THEN
00763                      IF( T( J+1, J ).NE.ZERO ) THEN
00764                         J2 = J + 1
00765                         JNXT = J + 2
00766                      END IF
00767                   END IF
00768 *
00769                   IF( J1.EQ.J2 ) THEN
00770 *
00771 *                    1-by-1 diagonal block
00772 *
00773 *                    Scale if necessary to avoid overflow when forming
00774 *                    the right-hand side.
00775 *
00776                      IF( WORK( J ).GT.VCRIT ) THEN
00777                         REC = ONE / VMAX
00778                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00779                         VMAX = ONE
00780                         VCRIT = BIGNUM
00781                      END IF
00782 *
00783                      WORK( J+N ) = WORK( J+N ) -
00784      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
00785      $                             WORK( KI+1+N ), 1 )
00786 *
00787 *                    Solve (T(J,J)-WR)**T*X = WORK
00788 *
00789                      CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00790      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00791      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00792 *
00793 *                    Scale if necessary
00794 *
00795                      IF( SCALE.NE.ONE )
00796      $                  CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00797                      WORK( J+N ) = X( 1, 1 )
00798                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
00799                      VCRIT = BIGNUM / VMAX
00800 *
00801                   ELSE
00802 *
00803 *                    2-by-2 diagonal block
00804 *
00805 *                    Scale if necessary to avoid overflow when forming
00806 *                    the right-hand side.
00807 *
00808                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00809                      IF( BETA.GT.VCRIT ) THEN
00810                         REC = ONE / VMAX
00811                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00812                         VMAX = ONE
00813                         VCRIT = BIGNUM
00814                      END IF
00815 *
00816                      WORK( J+N ) = WORK( J+N ) -
00817      $                             SDOT( J-KI-1, T( KI+1, J ), 1,
00818      $                             WORK( KI+1+N ), 1 )
00819 *
00820                      WORK( J+1+N ) = WORK( J+1+N ) -
00821      $                               SDOT( J-KI-1, T( KI+1, J+1 ), 1,
00822      $                               WORK( KI+1+N ), 1 )
00823 *
00824 *                    Solve
00825 *                      [T(J,J)-WR   T(J,J+1)     ]**T* X = SCALE*( WORK1 )
00826 *                      [T(J+1,J)    T(J+1,J+1)-WR]               ( WORK2 )
00827 *
00828                      CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
00829      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00830      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00831 *
00832 *                    Scale if necessary
00833 *
00834                      IF( SCALE.NE.ONE )
00835      $                  CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00836                      WORK( J+N ) = X( 1, 1 )
00837                      WORK( J+1+N ) = X( 2, 1 )
00838 *
00839                      VMAX = MAX( ABS( WORK( J+N ) ),
00840      $                      ABS( WORK( J+1+N ) ), VMAX )
00841                      VCRIT = BIGNUM / VMAX
00842 *
00843                   END IF
00844   170          CONTINUE
00845 *
00846 *              Copy the vector x or Q*x to VL and normalize.
00847 *
00848                IF( .NOT.OVER ) THEN
00849                   CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
00850 *
00851                   II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
00852                   REMAX = ONE / ABS( VL( II, IS ) )
00853                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00854 *
00855                   DO 180 K = 1, KI - 1
00856                      VL( K, IS ) = ZERO
00857   180             CONTINUE
00858 *
00859                ELSE
00860 *
00861                   IF( KI.LT.N )
00862      $               CALL SGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
00863      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
00864      $                           VL( 1, KI ), 1 )
00865 *
00866                   II = ISAMAX( N, VL( 1, KI ), 1 )
00867                   REMAX = ONE / ABS( VL( II, KI ) )
00868                   CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
00869 *
00870                END IF
00871 *
00872             ELSE
00873 *
00874 *              Complex left eigenvector.
00875 *
00876 *               Initial solve:
00877 *                 ((T(KI,KI)    T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
00878 *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
00879 *
00880                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
00881                   WORK( KI+N ) = WI / T( KI, KI+1 )
00882                   WORK( KI+1+N2 ) = ONE
00883                ELSE
00884                   WORK( KI+N ) = ONE
00885                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
00886                END IF
00887                WORK( KI+1+N ) = ZERO
00888                WORK( KI+N2 ) = ZERO
00889 *
00890 *              Form right-hand side
00891 *
00892                DO 190 K = KI + 2, N
00893                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
00894                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
00895   190          CONTINUE
00896 *
00897 *              Solve complex quasi-triangular system:
00898 *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
00899 *
00900                VMAX = ONE
00901                VCRIT = BIGNUM
00902 *
00903                JNXT = KI + 2
00904                DO 200 J = KI + 2, N
00905                   IF( J.LT.JNXT )
00906      $               GO TO 200
00907                   J1 = J
00908                   J2 = J
00909                   JNXT = J + 1
00910                   IF( J.LT.N ) THEN
00911                      IF( T( J+1, J ).NE.ZERO ) THEN
00912                         J2 = J + 1
00913                         JNXT = J + 2
00914                      END IF
00915                   END IF
00916 *
00917                   IF( J1.EQ.J2 ) THEN
00918 *
00919 *                    1-by-1 diagonal block
00920 *
00921 *                    Scale if necessary to avoid overflow when
00922 *                    forming the right-hand side elements.
00923 *
00924                      IF( WORK( J ).GT.VCRIT ) THEN
00925                         REC = ONE / VMAX
00926                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00927                         CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00928                         VMAX = ONE
00929                         VCRIT = BIGNUM
00930                      END IF
00931 *
00932                      WORK( J+N ) = WORK( J+N ) -
00933      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
00934      $                             WORK( KI+2+N ), 1 )
00935                      WORK( J+N2 ) = WORK( J+N2 ) -
00936      $                              SDOT( J-KI-2, T( KI+2, J ), 1,
00937      $                              WORK( KI+2+N2 ), 1 )
00938 *
00939 *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
00940 *
00941                      CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00942      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00943      $                            -WI, X, 2, SCALE, XNORM, IERR )
00944 *
00945 *                    Scale if necessary
00946 *
00947                      IF( SCALE.NE.ONE ) THEN
00948                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00949                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
00950                      END IF
00951                      WORK( J+N ) = X( 1, 1 )
00952                      WORK( J+N2 ) = X( 1, 2 )
00953                      VMAX = MAX( ABS( WORK( J+N ) ),
00954      $                      ABS( WORK( J+N2 ) ), VMAX )
00955                      VCRIT = BIGNUM / VMAX
00956 *
00957                   ELSE
00958 *
00959 *                    2-by-2 diagonal block
00960 *
00961 *                    Scale if necessary to avoid overflow when forming
00962 *                    the right-hand side elements.
00963 *
00964                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00965                      IF( BETA.GT.VCRIT ) THEN
00966                         REC = ONE / VMAX
00967                         CALL SSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00968                         CALL SSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00969                         VMAX = ONE
00970                         VCRIT = BIGNUM
00971                      END IF
00972 *
00973                      WORK( J+N ) = WORK( J+N ) -
00974      $                             SDOT( J-KI-2, T( KI+2, J ), 1,
00975      $                             WORK( KI+2+N ), 1 )
00976 *
00977                      WORK( J+N2 ) = WORK( J+N2 ) -
00978      $                              SDOT( J-KI-2, T( KI+2, J ), 1,
00979      $                              WORK( KI+2+N2 ), 1 )
00980 *
00981                      WORK( J+1+N ) = WORK( J+1+N ) -
00982      $                               SDOT( J-KI-2, T( KI+2, J+1 ), 1,
00983      $                               WORK( KI+2+N ), 1 )
00984 *
00985                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
00986      $                                SDOT( J-KI-2, T( KI+2, J+1 ), 1,
00987      $                                WORK( KI+2+N2 ), 1 )
00988 *
00989 *                    Solve 2-by-2 complex linear equation
00990 *                      ([T(j,j)   T(j,j+1)  ]**T-(wr-i*wi)*I)*X = SCALE*B
00991 *                      ([T(j+1,j) T(j+1,j+1)]               )
00992 *
00993                      CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
00994      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00995      $                            -WI, X, 2, SCALE, XNORM, IERR )
00996 *
00997 *                    Scale if necessary
00998 *
00999                      IF( SCALE.NE.ONE ) THEN
01000                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
01001                         CALL SSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
01002                      END IF
01003                      WORK( J+N ) = X( 1, 1 )
01004                      WORK( J+N2 ) = X( 1, 2 )
01005                      WORK( J+1+N ) = X( 2, 1 )
01006                      WORK( J+1+N2 ) = X( 2, 2 )
01007                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
01008      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
01009                      VCRIT = BIGNUM / VMAX
01010 *
01011                   END IF
01012   200          CONTINUE
01013 *
01014 *              Copy the vector x or Q*x to VL and normalize.
01015 *
01016                IF( .NOT.OVER ) THEN
01017                   CALL SCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
01018                   CALL SCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
01019      $                        1 )
01020 *
01021                   EMAX = ZERO
01022                   DO 220 K = KI, N
01023                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
01024      $                      ABS( VL( K, IS+1 ) ) )
01025   220             CONTINUE
01026                   REMAX = ONE / EMAX
01027                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
01028                   CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
01029 *
01030                   DO 230 K = 1, KI - 1
01031                      VL( K, IS ) = ZERO
01032                      VL( K, IS+1 ) = ZERO
01033   230             CONTINUE
01034                ELSE
01035                   IF( KI.LT.N-1 ) THEN
01036                      CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
01037      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
01038      $                           VL( 1, KI ), 1 )
01039                      CALL SGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
01040      $                           LDVL, WORK( KI+2+N2 ), 1,
01041      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
01042                   ELSE
01043                      CALL SSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
01044                      CALL SSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
01045                   END IF
01046 *
01047                   EMAX = ZERO
01048                   DO 240 K = 1, N
01049                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
01050      $                      ABS( VL( K, KI+1 ) ) )
01051   240             CONTINUE
01052                   REMAX = ONE / EMAX
01053                   CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
01054                   CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
01055 *
01056                END IF
01057 *
01058             END IF
01059 *
01060             IS = IS + 1
01061             IF( IP.NE.0 )
01062      $         IS = IS + 1
01063   250       CONTINUE
01064             IF( IP.EQ.-1 )
01065      $         IP = 0
01066             IF( IP.EQ.1 )
01067      $         IP = -1
01068 *
01069   260    CONTINUE
01070 *
01071       END IF
01072 *
01073       RETURN
01074 *
01075 *     End of STREVC
01076 *
01077       END
 All Files Functions