LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dtrevc.f
Go to the documentation of this file.
00001 *> \brief \b DTREVC
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DTREVC + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DTREVC( 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 *       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00031 *      $                   WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DTREVC 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 DHSEQR.
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DHSEQR).
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 DOUBLE PRECISION 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 DHSEQR).
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 DOUBLE PRECISION 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 doubleOTHERcomputational
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 DTREVC( 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       DOUBLE PRECISION   T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00237      $                   WORK( * )
00238 *     ..
00239 *
00240 *  =====================================================================
00241 *
00242 *     .. Parameters ..
00243       DOUBLE PRECISION   ZERO, ONE
00244       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+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       DOUBLE PRECISION   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            IDAMAX
00256       DOUBLE PRECISION   DDOT, DLAMCH
00257       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
00258 *     ..
00259 *     .. External Subroutines ..
00260       EXTERNAL           DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
00261 *     ..
00262 *     .. Intrinsic Functions ..
00263       INTRINSIC          ABS, MAX, SQRT
00264 *     ..
00265 *     .. Local Arrays ..
00266       DOUBLE PRECISION   X( 2, 2 )
00267 *     ..
00268 *     .. Executable Statements ..
00269 *
00270 *     Decode and test the input parameters
00271 *
00272       BOTHV = LSAME( SIDE, 'B' )
00273       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00274       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00275 *
00276       ALLV = LSAME( HOWMNY, 'A' )
00277       OVER = LSAME( HOWMNY, 'B' )
00278       SOMEV = LSAME( HOWMNY, 'S' )
00279 *
00280       INFO = 0
00281       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00282          INFO = -1
00283       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
00284          INFO = -2
00285       ELSE IF( N.LT.0 ) THEN
00286          INFO = -4
00287       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00288          INFO = -6
00289       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00290          INFO = -8
00291       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00292          INFO = -10
00293       ELSE
00294 *
00295 *        Set M to the number of columns required to store the selected
00296 *        eigenvectors, standardize the array SELECT if necessary, and
00297 *        test MM.
00298 *
00299          IF( SOMEV ) THEN
00300             M = 0
00301             PAIR = .FALSE.
00302             DO 10 J = 1, N
00303                IF( PAIR ) THEN
00304                   PAIR = .FALSE.
00305                   SELECT( J ) = .FALSE.
00306                ELSE
00307                   IF( J.LT.N ) THEN
00308                      IF( T( J+1, J ).EQ.ZERO ) THEN
00309                         IF( SELECT( J ) )
00310      $                     M = M + 1
00311                      ELSE
00312                         PAIR = .TRUE.
00313                         IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
00314                            SELECT( J ) = .TRUE.
00315                            M = M + 2
00316                         END IF
00317                      END IF
00318                   ELSE
00319                      IF( SELECT( N ) )
00320      $                  M = M + 1
00321                   END IF
00322                END IF
00323    10       CONTINUE
00324          ELSE
00325             M = N
00326          END IF
00327 *
00328          IF( MM.LT.M ) THEN
00329             INFO = -11
00330          END IF
00331       END IF
00332       IF( INFO.NE.0 ) THEN
00333          CALL XERBLA( 'DTREVC', -INFO )
00334          RETURN
00335       END IF
00336 *
00337 *     Quick return if possible.
00338 *
00339       IF( N.EQ.0 )
00340      $   RETURN
00341 *
00342 *     Set the constants to control overflow.
00343 *
00344       UNFL = DLAMCH( 'Safe minimum' )
00345       OVFL = ONE / UNFL
00346       CALL DLABAD( UNFL, OVFL )
00347       ULP = DLAMCH( 'Precision' )
00348       SMLNUM = UNFL*( N / ULP )
00349       BIGNUM = ( ONE-ULP ) / SMLNUM
00350 *
00351 *     Compute 1-norm of each column of strictly upper triangular
00352 *     part of T to control overflow in triangular solver.
00353 *
00354       WORK( 1 ) = ZERO
00355       DO 30 J = 2, N
00356          WORK( J ) = ZERO
00357          DO 20 I = 1, J - 1
00358             WORK( J ) = WORK( J ) + ABS( T( I, J ) )
00359    20    CONTINUE
00360    30 CONTINUE
00361 *
00362 *     Index IP is used to specify the real or complex eigenvalue:
00363 *       IP = 0, real eigenvalue,
00364 *            1, first of conjugate complex pair: (wr,wi)
00365 *           -1, second of conjugate complex pair: (wr,wi)
00366 *
00367       N2 = 2*N
00368 *
00369       IF( RIGHTV ) THEN
00370 *
00371 *        Compute right eigenvectors.
00372 *
00373          IP = 0
00374          IS = M
00375          DO 140 KI = N, 1, -1
00376 *
00377             IF( IP.EQ.1 )
00378      $         GO TO 130
00379             IF( KI.EQ.1 )
00380      $         GO TO 40
00381             IF( T( KI, KI-1 ).EQ.ZERO )
00382      $         GO TO 40
00383             IP = -1
00384 *
00385    40       CONTINUE
00386             IF( SOMEV ) THEN
00387                IF( IP.EQ.0 ) THEN
00388                   IF( .NOT.SELECT( KI ) )
00389      $               GO TO 130
00390                ELSE
00391                   IF( .NOT.SELECT( KI-1 ) )
00392      $               GO TO 130
00393                END IF
00394             END IF
00395 *
00396 *           Compute the KI-th eigenvalue (WR,WI).
00397 *
00398             WR = T( KI, KI )
00399             WI = ZERO
00400             IF( IP.NE.0 )
00401      $         WI = SQRT( ABS( T( KI, KI-1 ) ) )*
00402      $              SQRT( ABS( T( KI-1, KI ) ) )
00403             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00404 *
00405             IF( IP.EQ.0 ) THEN
00406 *
00407 *              Real right eigenvector
00408 *
00409                WORK( KI+N ) = ONE
00410 *
00411 *              Form right-hand side
00412 *
00413                DO 50 K = 1, KI - 1
00414                   WORK( K+N ) = -T( K, KI )
00415    50          CONTINUE
00416 *
00417 *              Solve the upper quasi-triangular system:
00418 *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
00419 *
00420                JNXT = KI - 1
00421                DO 60 J = KI - 1, 1, -1
00422                   IF( J.GT.JNXT )
00423      $               GO TO 60
00424                   J1 = J
00425                   J2 = J
00426                   JNXT = J - 1
00427                   IF( J.GT.1 ) THEN
00428                      IF( T( J, J-1 ).NE.ZERO ) THEN
00429                         J1 = J - 1
00430                         JNXT = J - 2
00431                      END IF
00432                   END IF
00433 *
00434                   IF( J1.EQ.J2 ) THEN
00435 *
00436 *                    1-by-1 diagonal block
00437 *
00438                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00439      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00440      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00441 *
00442 *                    Scale X(1,1) to avoid overflow when updating
00443 *                    the right-hand side.
00444 *
00445                      IF( XNORM.GT.ONE ) THEN
00446                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00447                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00448                            SCALE = SCALE / XNORM
00449                         END IF
00450                      END IF
00451 *
00452 *                    Scale if necessary
00453 *
00454                      IF( SCALE.NE.ONE )
00455      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00456                      WORK( J+N ) = X( 1, 1 )
00457 *
00458 *                    Update right-hand side
00459 *
00460                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00461      $                           WORK( 1+N ), 1 )
00462 *
00463                   ELSE
00464 *
00465 *                    2-by-2 diagonal block
00466 *
00467                      CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
00468      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00469      $                            WORK( J-1+N ), N, WR, ZERO, X, 2,
00470      $                            SCALE, XNORM, IERR )
00471 *
00472 *                    Scale X(1,1) and X(2,1) to avoid overflow when
00473 *                    updating the right-hand side.
00474 *
00475                      IF( XNORM.GT.ONE ) THEN
00476                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00477                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00478                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00479                            X( 2, 1 ) = X( 2, 1 ) / XNORM
00480                            SCALE = SCALE / XNORM
00481                         END IF
00482                      END IF
00483 *
00484 *                    Scale if necessary
00485 *
00486                      IF( SCALE.NE.ONE )
00487      $                  CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00488                      WORK( J-1+N ) = X( 1, 1 )
00489                      WORK( J+N ) = X( 2, 1 )
00490 *
00491 *                    Update right-hand side
00492 *
00493                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00494      $                           WORK( 1+N ), 1 )
00495                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00496      $                           WORK( 1+N ), 1 )
00497                   END IF
00498    60          CONTINUE
00499 *
00500 *              Copy the vector x or Q*x to VR and normalize.
00501 *
00502                IF( .NOT.OVER ) THEN
00503                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
00504 *
00505                   II = IDAMAX( KI, VR( 1, IS ), 1 )
00506                   REMAX = ONE / ABS( VR( II, IS ) )
00507                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
00508 *
00509                   DO 70 K = KI + 1, N
00510                      VR( K, IS ) = ZERO
00511    70             CONTINUE
00512                ELSE
00513                   IF( KI.GT.1 )
00514      $               CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
00515      $                           WORK( 1+N ), 1, WORK( KI+N ),
00516      $                           VR( 1, KI ), 1 )
00517 *
00518                   II = IDAMAX( N, VR( 1, KI ), 1 )
00519                   REMAX = ONE / ABS( VR( II, KI ) )
00520                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
00521                END IF
00522 *
00523             ELSE
00524 *
00525 *              Complex right eigenvector.
00526 *
00527 *              Initial solve
00528 *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
00529 *                [ (T(KI,KI-1)   T(KI,KI)   )               ]
00530 *
00531                IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
00532                   WORK( KI-1+N ) = ONE
00533                   WORK( KI+N2 ) = WI / T( KI-1, KI )
00534                ELSE
00535                   WORK( KI-1+N ) = -WI / T( KI, KI-1 )
00536                   WORK( KI+N2 ) = ONE
00537                END IF
00538                WORK( KI+N ) = ZERO
00539                WORK( KI-1+N2 ) = ZERO
00540 *
00541 *              Form right-hand side
00542 *
00543                DO 80 K = 1, KI - 2
00544                   WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
00545                   WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
00546    80          CONTINUE
00547 *
00548 *              Solve upper quasi-triangular system:
00549 *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
00550 *
00551                JNXT = KI - 2
00552                DO 90 J = KI - 2, 1, -1
00553                   IF( J.GT.JNXT )
00554      $               GO TO 90
00555                   J1 = J
00556                   J2 = J
00557                   JNXT = J - 1
00558                   IF( J.GT.1 ) THEN
00559                      IF( T( J, J-1 ).NE.ZERO ) THEN
00560                         J1 = J - 1
00561                         JNXT = J - 2
00562                      END IF
00563                   END IF
00564 *
00565                   IF( J1.EQ.J2 ) THEN
00566 *
00567 *                    1-by-1 diagonal block
00568 *
00569                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00570      $                            LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
00571      $                            X, 2, SCALE, XNORM, IERR )
00572 *
00573 *                    Scale X(1,1) and X(1,2) to avoid overflow when
00574 *                    updating the right-hand side.
00575 *
00576                      IF( XNORM.GT.ONE ) THEN
00577                         IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
00578                            X( 1, 1 ) = X( 1, 1 ) / XNORM
00579                            X( 1, 2 ) = X( 1, 2 ) / XNORM
00580                            SCALE = SCALE / XNORM
00581                         END IF
00582                      END IF
00583 *
00584 *                    Scale if necessary
00585 *
00586                      IF( SCALE.NE.ONE ) THEN
00587                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00588                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00589                      END IF
00590                      WORK( J+N ) = X( 1, 1 )
00591                      WORK( J+N2 ) = X( 1, 2 )
00592 *
00593 *                    Update the right-hand side
00594 *
00595                      CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
00596      $                           WORK( 1+N ), 1 )
00597                      CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
00598      $                           WORK( 1+N2 ), 1 )
00599 *
00600                   ELSE
00601 *
00602 *                    2-by-2 diagonal block
00603 *
00604                      CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
00605      $                            T( J-1, J-1 ), LDT, ONE, ONE,
00606      $                            WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
00607      $                            XNORM, IERR )
00608 *
00609 *                    Scale X to avoid overflow when updating
00610 *                    the right-hand side.
00611 *
00612                      IF( XNORM.GT.ONE ) THEN
00613                         BETA = MAX( WORK( J-1 ), WORK( J ) )
00614                         IF( BETA.GT.BIGNUM / XNORM ) THEN
00615                            REC = ONE / XNORM
00616                            X( 1, 1 ) = X( 1, 1 )*REC
00617                            X( 1, 2 ) = X( 1, 2 )*REC
00618                            X( 2, 1 ) = X( 2, 1 )*REC
00619                            X( 2, 2 ) = X( 2, 2 )*REC
00620                            SCALE = SCALE*REC
00621                         END IF
00622                      END IF
00623 *
00624 *                    Scale if necessary
00625 *
00626                      IF( SCALE.NE.ONE ) THEN
00627                         CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
00628                         CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
00629                      END IF
00630                      WORK( J-1+N ) = X( 1, 1 )
00631                      WORK( J+N ) = X( 2, 1 )
00632                      WORK( J-1+N2 ) = X( 1, 2 )
00633                      WORK( J+N2 ) = X( 2, 2 )
00634 *
00635 *                    Update the right-hand side
00636 *
00637                      CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
00638      $                           WORK( 1+N ), 1 )
00639                      CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
00640      $                           WORK( 1+N ), 1 )
00641                      CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
00642      $                           WORK( 1+N2 ), 1 )
00643                      CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
00644      $                           WORK( 1+N2 ), 1 )
00645                   END IF
00646    90          CONTINUE
00647 *
00648 *              Copy the vector x or Q*x to VR and normalize.
00649 *
00650                IF( .NOT.OVER ) THEN
00651                   CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
00652                   CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
00653 *
00654                   EMAX = ZERO
00655                   DO 100 K = 1, KI
00656                      EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
00657      $                      ABS( VR( K, IS ) ) )
00658   100             CONTINUE
00659 *
00660                   REMAX = ONE / EMAX
00661                   CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
00662                   CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
00663 *
00664                   DO 110 K = KI + 1, N
00665                      VR( K, IS-1 ) = ZERO
00666                      VR( K, IS ) = ZERO
00667   110             CONTINUE
00668 *
00669                ELSE
00670 *
00671                   IF( KI.GT.2 ) THEN
00672                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00673      $                           WORK( 1+N ), 1, WORK( KI-1+N ),
00674      $                           VR( 1, KI-1 ), 1 )
00675                      CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
00676      $                           WORK( 1+N2 ), 1, WORK( KI+N2 ),
00677      $                           VR( 1, KI ), 1 )
00678                   ELSE
00679                      CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
00680                      CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
00681                   END IF
00682 *
00683                   EMAX = ZERO
00684                   DO 120 K = 1, N
00685                      EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
00686      $                      ABS( VR( K, KI ) ) )
00687   120             CONTINUE
00688                   REMAX = ONE / EMAX
00689                   CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
00690                   CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
00691                END IF
00692             END IF
00693 *
00694             IS = IS - 1
00695             IF( IP.NE.0 )
00696      $         IS = IS - 1
00697   130       CONTINUE
00698             IF( IP.EQ.1 )
00699      $         IP = 0
00700             IF( IP.EQ.-1 )
00701      $         IP = 1
00702   140    CONTINUE
00703       END IF
00704 *
00705       IF( LEFTV ) THEN
00706 *
00707 *        Compute left eigenvectors.
00708 *
00709          IP = 0
00710          IS = 1
00711          DO 260 KI = 1, N
00712 *
00713             IF( IP.EQ.-1 )
00714      $         GO TO 250
00715             IF( KI.EQ.N )
00716      $         GO TO 150
00717             IF( T( KI+1, KI ).EQ.ZERO )
00718      $         GO TO 150
00719             IP = 1
00720 *
00721   150       CONTINUE
00722             IF( SOMEV ) THEN
00723                IF( .NOT.SELECT( KI ) )
00724      $            GO TO 250
00725             END IF
00726 *
00727 *           Compute the KI-th eigenvalue (WR,WI).
00728 *
00729             WR = T( KI, KI )
00730             WI = ZERO
00731             IF( IP.NE.0 )
00732      $         WI = SQRT( ABS( T( KI, KI+1 ) ) )*
00733      $              SQRT( ABS( T( KI+1, KI ) ) )
00734             SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
00735 *
00736             IF( IP.EQ.0 ) THEN
00737 *
00738 *              Real left eigenvector.
00739 *
00740                WORK( KI+N ) = ONE
00741 *
00742 *              Form right-hand side
00743 *
00744                DO 160 K = KI + 1, N
00745                   WORK( K+N ) = -T( KI, K )
00746   160          CONTINUE
00747 *
00748 *              Solve the quasi-triangular system:
00749 *                 (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
00750 *
00751                VMAX = ONE
00752                VCRIT = BIGNUM
00753 *
00754                JNXT = KI + 1
00755                DO 170 J = KI + 1, N
00756                   IF( J.LT.JNXT )
00757      $               GO TO 170
00758                   J1 = J
00759                   J2 = J
00760                   JNXT = J + 1
00761                   IF( J.LT.N ) THEN
00762                      IF( T( J+1, J ).NE.ZERO ) THEN
00763                         J2 = J + 1
00764                         JNXT = J + 2
00765                      END IF
00766                   END IF
00767 *
00768                   IF( J1.EQ.J2 ) THEN
00769 *
00770 *                    1-by-1 diagonal block
00771 *
00772 *                    Scale if necessary to avoid overflow when forming
00773 *                    the right-hand side.
00774 *
00775                      IF( WORK( J ).GT.VCRIT ) THEN
00776                         REC = ONE / VMAX
00777                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00778                         VMAX = ONE
00779                         VCRIT = BIGNUM
00780                      END IF
00781 *
00782                      WORK( J+N ) = WORK( J+N ) -
00783      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
00784      $                             WORK( KI+1+N ), 1 )
00785 *
00786 *                    Solve (T(J,J)-WR)**T*X = WORK
00787 *
00788                      CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
00789      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00790      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00791 *
00792 *                    Scale if necessary
00793 *
00794                      IF( SCALE.NE.ONE )
00795      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00796                      WORK( J+N ) = X( 1, 1 )
00797                      VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
00798                      VCRIT = BIGNUM / VMAX
00799 *
00800                   ELSE
00801 *
00802 *                    2-by-2 diagonal block
00803 *
00804 *                    Scale if necessary to avoid overflow when forming
00805 *                    the right-hand side.
00806 *
00807                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00808                      IF( BETA.GT.VCRIT ) THEN
00809                         REC = ONE / VMAX
00810                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00811                         VMAX = ONE
00812                         VCRIT = BIGNUM
00813                      END IF
00814 *
00815                      WORK( J+N ) = WORK( J+N ) -
00816      $                             DDOT( J-KI-1, T( KI+1, J ), 1,
00817      $                             WORK( KI+1+N ), 1 )
00818 *
00819                      WORK( J+1+N ) = WORK( J+1+N ) -
00820      $                               DDOT( J-KI-1, T( KI+1, J+1 ), 1,
00821      $                               WORK( KI+1+N ), 1 )
00822 *
00823 *                    Solve
00824 *                      [T(J,J)-WR   T(J,J+1)     ]**T * X = SCALE*( WORK1 )
00825 *                      [T(J+1,J)    T(J+1,J+1)-WR]                ( WORK2 )
00826 *
00827                      CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
00828      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00829      $                            ZERO, X, 2, SCALE, XNORM, IERR )
00830 *
00831 *                    Scale if necessary
00832 *
00833                      IF( SCALE.NE.ONE )
00834      $                  CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00835                      WORK( J+N ) = X( 1, 1 )
00836                      WORK( J+1+N ) = X( 2, 1 )
00837 *
00838                      VMAX = MAX( ABS( WORK( J+N ) ),
00839      $                      ABS( WORK( J+1+N ) ), VMAX )
00840                      VCRIT = BIGNUM / VMAX
00841 *
00842                   END IF
00843   170          CONTINUE
00844 *
00845 *              Copy the vector x or Q*x to VL and normalize.
00846 *
00847                IF( .NOT.OVER ) THEN
00848                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
00849 *
00850                   II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
00851                   REMAX = ONE / ABS( VL( II, IS ) )
00852                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00853 *
00854                   DO 180 K = 1, KI - 1
00855                      VL( K, IS ) = ZERO
00856   180             CONTINUE
00857 *
00858                ELSE
00859 *
00860                   IF( KI.LT.N )
00861      $               CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
00862      $                           WORK( KI+1+N ), 1, WORK( KI+N ),
00863      $                           VL( 1, KI ), 1 )
00864 *
00865                   II = IDAMAX( N, VL( 1, KI ), 1 )
00866                   REMAX = ONE / ABS( VL( II, KI ) )
00867                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
00868 *
00869                END IF
00870 *
00871             ELSE
00872 *
00873 *              Complex left eigenvector.
00874 *
00875 *               Initial solve:
00876 *                 ((T(KI,KI)    T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
00877 *                 ((T(KI+1,KI) T(KI+1,KI+1))                )
00878 *
00879                IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
00880                   WORK( KI+N ) = WI / T( KI, KI+1 )
00881                   WORK( KI+1+N2 ) = ONE
00882                ELSE
00883                   WORK( KI+N ) = ONE
00884                   WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
00885                END IF
00886                WORK( KI+1+N ) = ZERO
00887                WORK( KI+N2 ) = ZERO
00888 *
00889 *              Form right-hand side
00890 *
00891                DO 190 K = KI + 2, N
00892                   WORK( K+N ) = -WORK( KI+N )*T( KI, K )
00893                   WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
00894   190          CONTINUE
00895 *
00896 *              Solve complex quasi-triangular system:
00897 *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
00898 *
00899                VMAX = ONE
00900                VCRIT = BIGNUM
00901 *
00902                JNXT = KI + 2
00903                DO 200 J = KI + 2, N
00904                   IF( J.LT.JNXT )
00905      $               GO TO 200
00906                   J1 = J
00907                   J2 = J
00908                   JNXT = J + 1
00909                   IF( J.LT.N ) THEN
00910                      IF( T( J+1, J ).NE.ZERO ) THEN
00911                         J2 = J + 1
00912                         JNXT = J + 2
00913                      END IF
00914                   END IF
00915 *
00916                   IF( J1.EQ.J2 ) THEN
00917 *
00918 *                    1-by-1 diagonal block
00919 *
00920 *                    Scale if necessary to avoid overflow when
00921 *                    forming the right-hand side elements.
00922 *
00923                      IF( WORK( J ).GT.VCRIT ) THEN
00924                         REC = ONE / VMAX
00925                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00926                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00927                         VMAX = ONE
00928                         VCRIT = BIGNUM
00929                      END IF
00930 *
00931                      WORK( J+N ) = WORK( J+N ) -
00932      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
00933      $                             WORK( KI+2+N ), 1 )
00934                      WORK( J+N2 ) = WORK( J+N2 ) -
00935      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
00936      $                              WORK( KI+2+N2 ), 1 )
00937 *
00938 *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
00939 *
00940                      CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
00941      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00942      $                            -WI, X, 2, SCALE, XNORM, IERR )
00943 *
00944 *                    Scale if necessary
00945 *
00946                      IF( SCALE.NE.ONE ) THEN
00947                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
00948                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
00949                      END IF
00950                      WORK( J+N ) = X( 1, 1 )
00951                      WORK( J+N2 ) = X( 1, 2 )
00952                      VMAX = MAX( ABS( WORK( J+N ) ),
00953      $                      ABS( WORK( J+N2 ) ), VMAX )
00954                      VCRIT = BIGNUM / VMAX
00955 *
00956                   ELSE
00957 *
00958 *                    2-by-2 diagonal block
00959 *
00960 *                    Scale if necessary to avoid overflow when forming
00961 *                    the right-hand side elements.
00962 *
00963                      BETA = MAX( WORK( J ), WORK( J+1 ) )
00964                      IF( BETA.GT.VCRIT ) THEN
00965                         REC = ONE / VMAX
00966                         CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
00967                         CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
00968                         VMAX = ONE
00969                         VCRIT = BIGNUM
00970                      END IF
00971 *
00972                      WORK( J+N ) = WORK( J+N ) -
00973      $                             DDOT( J-KI-2, T( KI+2, J ), 1,
00974      $                             WORK( KI+2+N ), 1 )
00975 *
00976                      WORK( J+N2 ) = WORK( J+N2 ) -
00977      $                              DDOT( J-KI-2, T( KI+2, J ), 1,
00978      $                              WORK( KI+2+N2 ), 1 )
00979 *
00980                      WORK( J+1+N ) = WORK( J+1+N ) -
00981      $                               DDOT( J-KI-2, T( KI+2, J+1 ), 1,
00982      $                               WORK( KI+2+N ), 1 )
00983 *
00984                      WORK( J+1+N2 ) = WORK( J+1+N2 ) -
00985      $                                DDOT( J-KI-2, T( KI+2, J+1 ), 1,
00986      $                                WORK( KI+2+N2 ), 1 )
00987 *
00988 *                    Solve 2-by-2 complex linear equation
00989 *                      ([T(j,j)   T(j,j+1)  ]**T-(wr-i*wi)*I)*X = SCALE*B
00990 *                      ([T(j+1,j) T(j+1,j+1)]               )
00991 *
00992                      CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
00993      $                            LDT, ONE, ONE, WORK( J+N ), N, WR,
00994      $                            -WI, X, 2, SCALE, XNORM, IERR )
00995 *
00996 *                    Scale if necessary
00997 *
00998                      IF( SCALE.NE.ONE ) THEN
00999                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
01000                         CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
01001                      END IF
01002                      WORK( J+N ) = X( 1, 1 )
01003                      WORK( J+N2 ) = X( 1, 2 )
01004                      WORK( J+1+N ) = X( 2, 1 )
01005                      WORK( J+1+N2 ) = X( 2, 2 )
01006                      VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
01007      $                      ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
01008                      VCRIT = BIGNUM / VMAX
01009 *
01010                   END IF
01011   200          CONTINUE
01012 *
01013 *              Copy the vector x or Q*x to VL and normalize.
01014 *
01015                IF( .NOT.OVER ) THEN
01016                   CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
01017                   CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
01018      $                        1 )
01019 *
01020                   EMAX = ZERO
01021                   DO 220 K = KI, N
01022                      EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
01023      $                      ABS( VL( K, IS+1 ) ) )
01024   220             CONTINUE
01025                   REMAX = ONE / EMAX
01026                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
01027                   CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
01028 *
01029                   DO 230 K = 1, KI - 1
01030                      VL( K, IS ) = ZERO
01031                      VL( K, IS+1 ) = ZERO
01032   230             CONTINUE
01033                ELSE
01034                   IF( KI.LT.N-1 ) THEN
01035                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
01036      $                           LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
01037      $                           VL( 1, KI ), 1 )
01038                      CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
01039      $                           LDVL, WORK( KI+2+N2 ), 1,
01040      $                           WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
01041                   ELSE
01042                      CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
01043                      CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
01044                   END IF
01045 *
01046                   EMAX = ZERO
01047                   DO 240 K = 1, N
01048                      EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
01049      $                      ABS( VL( K, KI+1 ) ) )
01050   240             CONTINUE
01051                   REMAX = ONE / EMAX
01052                   CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
01053                   CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
01054 *
01055                END IF
01056 *
01057             END IF
01058 *
01059             IS = IS + 1
01060             IF( IP.NE.0 )
01061      $         IS = IS + 1
01062   250       CONTINUE
01063             IF( IP.EQ.-1 )
01064      $         IP = 0
01065             IF( IP.EQ.1 )
01066      $         IP = -1
01067 *
01068   260    CONTINUE
01069 *
01070       END IF
01071 *
01072       RETURN
01073 *
01074 *     End of DTREVC
01075 *
01076       END
 All Files Functions