LAPACK: SRC/ctrevc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
00002      $                   LDVR, MM, M, WORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          HOWMNY, SIDE
00011       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       LOGICAL            SELECT( * )
00015       REAL               RWORK( * )
00016       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00017      $                   WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CTREVC computes some or all of the right and/or left eigenvectors of
00024 *  a complex upper triangular matrix T.
00025 *  Matrices of this type are produced by the Schur factorization of
00026 *  a complex general matrix:  A = Q*T*Q**H, as computed by CHSEQR.
00027 *  
00028 *  The right eigenvector x and the left eigenvector y of T corresponding
00029 *  to an eigenvalue w are defined by:
00030 *  
00031 *               T*x = w*x,     (y**H)*T = w*(y**H)
00032 *  
00033 *  where y**H denotes the conjugate transpose of the vector y.
00034 *  The eigenvalues are not input to this routine, but are read directly
00035 *  from the diagonal of T.
00036 *  
00037 *  This routine returns the matrices X and/or Y of right and left
00038 *  eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
00039 *  input matrix.  If Q is the unitary factor that reduces a matrix A to
00040 *  Schur form T, then Q*X and Q*Y are the matrices of right and left
00041 *  eigenvectors of A.
00042 *
00043 *  Arguments
00044 *  =========
00045 *
00046 *  SIDE    (input) CHARACTER*1
00047 *          = 'R':  compute right eigenvectors only;
00048 *          = 'L':  compute left eigenvectors only;
00049 *          = 'B':  compute both right and left eigenvectors.
00050 *
00051 *  HOWMNY  (input) CHARACTER*1
00052 *          = 'A':  compute all right and/or left eigenvectors;
00053 *          = 'B':  compute all right and/or left eigenvectors,
00054 *                  backtransformed using the matrices supplied in
00055 *                  VR and/or VL;
00056 *          = 'S':  compute selected right and/or left eigenvectors,
00057 *                  as indicated by the logical array SELECT.
00058 *
00059 *  SELECT  (input) LOGICAL array, dimension (N)
00060 *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
00061 *          computed.
00062 *          The eigenvector corresponding to the j-th eigenvalue is
00063 *          computed if SELECT(j) = .TRUE..
00064 *          Not referenced if HOWMNY = 'A' or 'B'.
00065 *
00066 *  N       (input) INTEGER
00067 *          The order of the matrix T. N >= 0.
00068 *
00069 *  T       (input/output) COMPLEX array, dimension (LDT,N)
00070 *          The upper triangular matrix T.  T is modified, but restored
00071 *          on exit.
00072 *
00073 *  LDT     (input) INTEGER
00074 *          The leading dimension of the array T. LDT >= max(1,N).
00075 *
00076 *  VL      (input/output) COMPLEX array, dimension (LDVL,MM)
00077 *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
00078 *          contain an N-by-N matrix Q (usually the unitary matrix Q of
00079 *          Schur vectors returned by CHSEQR).
00080 *          On exit, if SIDE = 'L' or 'B', VL contains:
00081 *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
00082 *          if HOWMNY = 'B', the matrix Q*Y;
00083 *          if HOWMNY = 'S', the left eigenvectors of T specified by
00084 *                           SELECT, stored consecutively in the columns
00085 *                           of VL, in the same order as their
00086 *                           eigenvalues.
00087 *          Not referenced if SIDE = 'R'.
00088 *
00089 *  LDVL    (input) INTEGER
00090 *          The leading dimension of the array VL.  LDVL >= 1, and if
00091 *          SIDE = 'L' or 'B', LDVL >= N.
00092 *
00093 *  VR      (input/output) COMPLEX array, dimension (LDVR,MM)
00094 *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
00095 *          contain an N-by-N matrix Q (usually the unitary matrix Q of
00096 *          Schur vectors returned by CHSEQR).
00097 *          On exit, if SIDE = 'R' or 'B', VR contains:
00098 *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
00099 *          if HOWMNY = 'B', the matrix Q*X;
00100 *          if HOWMNY = 'S', the right eigenvectors of T specified by
00101 *                           SELECT, stored consecutively in the columns
00102 *                           of VR, in the same order as their
00103 *                           eigenvalues.
00104 *          Not referenced if SIDE = 'L'.
00105 *
00106 *  LDVR    (input) INTEGER
00107 *          The leading dimension of the array VR.  LDVR >= 1, and if
00108 *          SIDE = 'R' or 'B'; LDVR >= N.
00109 *
00110 *  MM      (input) INTEGER
00111 *          The number of columns in the arrays VL and/or VR. MM >= M.
00112 *
00113 *  M       (output) INTEGER
00114 *          The number of columns in the arrays VL and/or VR actually
00115 *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
00116 *          is set to N.  Each selected eigenvector occupies one
00117 *          column.
00118 *
00119 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00120 *
00121 *  RWORK   (workspace) REAL array, dimension (N)
00122 *
00123 *  INFO    (output) INTEGER
00124 *          = 0:  successful exit
00125 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00126 *
00127 *  Further Details
00128 *  ===============
00129 *
00130 *  The algorithm used in this program is basically backward (forward)
00131 *  substitution, with scaling to make the the code robust against
00132 *  possible overflow.
00133 *
00134 *  Each eigenvector is normalized so that the element of largest
00135 *  magnitude has magnitude 1; here the magnitude of a complex number
00136 *  (x,y) is taken to be |x| + |y|.
00137 *
00138 *  =====================================================================
00139 *
00140 *     .. Parameters ..
00141       REAL               ZERO, ONE
00142       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00143       COMPLEX            CMZERO, CMONE
00144       PARAMETER          ( CMZERO = ( 0.0E+0, 0.0E+0 ),
00145      $                   CMONE = ( 1.0E+0, 0.0E+0 ) )
00146 *     ..
00147 *     .. Local Scalars ..
00148       LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
00149       INTEGER            I, II, IS, J, K, KI
00150       REAL               OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
00151       COMPLEX            CDUM
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       INTEGER            ICAMAX
00156       REAL               SCASUM, SLAMCH
00157       EXTERNAL           LSAME, ICAMAX, SCASUM, SLAMCH
00158 *     ..
00159 *     .. External Subroutines ..
00160       EXTERNAL           CCOPY, CGEMV, CLATRS, CSSCAL, SLABAD, XERBLA
00161 *     ..
00162 *     .. Intrinsic Functions ..
00163       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL
00164 *     ..
00165 *     .. Statement Functions ..
00166       REAL               CABS1
00167 *     ..
00168 *     .. Statement Function definitions ..
00169       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00170 *     ..
00171 *     .. Executable Statements ..
00172 *
00173 *     Decode and test the input parameters
00174 *
00175       BOTHV = LSAME( SIDE, 'B' )
00176       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00177       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00178 *
00179       ALLV = LSAME( HOWMNY, 'A' )
00180       OVER = LSAME( HOWMNY, 'B' )
00181       SOMEV = LSAME( HOWMNY, 'S' )
00182 *
00183 *     Set M to the number of columns required to store the selected
00184 *     eigenvectors.
00185 *
00186       IF( SOMEV ) THEN
00187          M = 0
00188          DO 10 J = 1, N
00189             IF( SELECT( J ) )
00190      $         M = M + 1
00191    10    CONTINUE
00192       ELSE
00193          M = N
00194       END IF
00195 *
00196       INFO = 0
00197       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00198          INFO = -1
00199       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
00200          INFO = -2
00201       ELSE IF( N.LT.0 ) THEN
00202          INFO = -4
00203       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
00204          INFO = -6
00205       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00206          INFO = -8
00207       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00208          INFO = -10
00209       ELSE IF( MM.LT.M ) THEN
00210          INFO = -11
00211       END IF
00212       IF( INFO.NE.0 ) THEN
00213          CALL XERBLA( 'CTREVC', -INFO )
00214          RETURN
00215       END IF
00216 *
00217 *     Quick return if possible.
00218 *
00219       IF( N.EQ.0 )
00220      $   RETURN
00221 *
00222 *     Set the constants to control overflow.
00223 *
00224       UNFL = SLAMCH( 'Safe minimum' )
00225       OVFL = ONE / UNFL
00226       CALL SLABAD( UNFL, OVFL )
00227       ULP = SLAMCH( 'Precision' )
00228       SMLNUM = UNFL*( N / ULP )
00229 *
00230 *     Store the diagonal elements of T in working array WORK.
00231 *
00232       DO 20 I = 1, N
00233          WORK( I+N ) = T( I, I )
00234    20 CONTINUE
00235 *
00236 *     Compute 1-norm of each column of strictly upper triangular
00237 *     part of T to control overflow in triangular solver.
00238 *
00239       RWORK( 1 ) = ZERO
00240       DO 30 J = 2, N
00241          RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 )
00242    30 CONTINUE
00243 *
00244       IF( RIGHTV ) THEN
00245 *
00246 *        Compute right eigenvectors.
00247 *
00248          IS = M
00249          DO 80 KI = N, 1, -1
00250 *
00251             IF( SOMEV ) THEN
00252                IF( .NOT.SELECT( KI ) )
00253      $            GO TO 80
00254             END IF
00255             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
00256 *
00257             WORK( 1 ) = CMONE
00258 *
00259 *           Form right-hand side.
00260 *
00261             DO 40 K = 1, KI - 1
00262                WORK( K ) = -T( K, KI )
00263    40       CONTINUE
00264 *
00265 *           Solve the triangular system:
00266 *              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
00267 *
00268             DO 50 K = 1, KI - 1
00269                T( K, K ) = T( K, K ) - T( KI, KI )
00270                IF( CABS1( T( K, K ) ).LT.SMIN )
00271      $            T( K, K ) = SMIN
00272    50       CONTINUE
00273 *
00274             IF( KI.GT.1 ) THEN
00275                CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
00276      $                      KI-1, T, LDT, WORK( 1 ), SCALE, RWORK,
00277      $                      INFO )
00278                WORK( KI ) = SCALE
00279             END IF
00280 *
00281 *           Copy the vector x or Q*x to VR and normalize.
00282 *
00283             IF( .NOT.OVER ) THEN
00284                CALL CCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
00285 *
00286                II = ICAMAX( KI, VR( 1, IS ), 1 )
00287                REMAX = ONE / CABS1( VR( II, IS ) )
00288                CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 )
00289 *
00290                DO 60 K = KI + 1, N
00291                   VR( K, IS ) = CMZERO
00292    60          CONTINUE
00293             ELSE
00294                IF( KI.GT.1 )
00295      $            CALL CGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
00296      $                        1, CMPLX( SCALE ), VR( 1, KI ), 1 )
00297 *
00298                II = ICAMAX( N, VR( 1, KI ), 1 )
00299                REMAX = ONE / CABS1( VR( II, KI ) )
00300                CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 )
00301             END IF
00302 *
00303 *           Set back the original diagonal elements of T.
00304 *
00305             DO 70 K = 1, KI - 1
00306                T( K, K ) = WORK( K+N )
00307    70       CONTINUE
00308 *
00309             IS = IS - 1
00310    80    CONTINUE
00311       END IF
00312 *
00313       IF( LEFTV ) THEN
00314 *
00315 *        Compute left eigenvectors.
00316 *
00317          IS = 1
00318          DO 130 KI = 1, N
00319 *
00320             IF( SOMEV ) THEN
00321                IF( .NOT.SELECT( KI ) )
00322      $            GO TO 130
00323             END IF
00324             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
00325 *
00326             WORK( N ) = CMONE
00327 *
00328 *           Form right-hand side.
00329 *
00330             DO 90 K = KI + 1, N
00331                WORK( K ) = -CONJG( T( KI, K ) )
00332    90       CONTINUE
00333 *
00334 *           Solve the triangular system:
00335 *              (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK.
00336 *
00337             DO 100 K = KI + 1, N
00338                T( K, K ) = T( K, K ) - T( KI, KI )
00339                IF( CABS1( T( K, K ) ).LT.SMIN )
00340      $            T( K, K ) = SMIN
00341   100       CONTINUE
00342 *
00343             IF( KI.LT.N ) THEN
00344                CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
00345      $                      'Y', N-KI, T( KI+1, KI+1 ), LDT,
00346      $                      WORK( KI+1 ), SCALE, RWORK, INFO )
00347                WORK( KI ) = SCALE
00348             END IF
00349 *
00350 *           Copy the vector x or Q*x to VL and normalize.
00351 *
00352             IF( .NOT.OVER ) THEN
00353                CALL CCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
00354 *
00355                II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
00356                REMAX = ONE / CABS1( VL( II, IS ) )
00357                CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
00358 *
00359                DO 110 K = 1, KI - 1
00360                   VL( K, IS ) = CMZERO
00361   110          CONTINUE
00362             ELSE
00363                IF( KI.LT.N )
00364      $            CALL CGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
00365      $                        WORK( KI+1 ), 1, CMPLX( SCALE ),
00366      $                        VL( 1, KI ), 1 )
00367 *
00368                II = ICAMAX( N, VL( 1, KI ), 1 )
00369                REMAX = ONE / CABS1( VL( II, KI ) )
00370                CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 )
00371             END IF
00372 *
00373 *           Set back the original diagonal elements of T.
00374 *
00375             DO 120 K = KI + 1, N
00376                T( K, K ) = WORK( K+N )
00377   120       CONTINUE
00378 *
00379             IS = IS + 1
00380   130    CONTINUE
00381       END IF
00382 *
00383       RETURN
00384 *
00385 *     End of CTREVC
00386 *
00387       END