LAPACK: SRC/ztrevc.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZTREVC( 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       DOUBLE PRECISION   RWORK( * )
00016       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
00017      $                   WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  ZTREVC 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 ZHSEQR.
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*16 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*16 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 ZHSEQR).
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*16 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 ZHSEQR).
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*16 array, dimension (2*N)
00120 *
00121 *  RWORK   (workspace) DOUBLE PRECISION 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       DOUBLE PRECISION   ZERO, ONE
00142       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00143       COMPLEX*16         CMZERO, CMONE
00144       PARAMETER          ( CMZERO = ( 0.0D+0, 0.0D+0 ),
00145      $                   CMONE = ( 1.0D+0, 0.0D+0 ) )
00146 *     ..
00147 *     .. Local Scalars ..
00148       LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
00149       INTEGER            I, II, IS, J, K, KI
00150       DOUBLE PRECISION   OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
00151       COMPLEX*16         CDUM
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       INTEGER            IZAMAX
00156       DOUBLE PRECISION   DLAMCH, DZASUM
00157       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZASUM
00158 *     ..
00159 *     .. External Subroutines ..
00160       EXTERNAL           XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
00161 *     ..
00162 *     .. Intrinsic Functions ..
00163       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
00164 *     ..
00165 *     .. Statement Functions ..
00166       DOUBLE PRECISION   CABS1
00167 *     ..
00168 *     .. Statement Function definitions ..
00169       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( 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( 'ZTREVC', -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 = DLAMCH( 'Safe minimum' )
00225       OVFL = ONE / UNFL
00226       CALL DLABAD( UNFL, OVFL )
00227       ULP = DLAMCH( '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 ) = DZASUM( 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 ZLATRS( '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 ZCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
00285 *
00286                II = IZAMAX( KI, VR( 1, IS ), 1 )
00287                REMAX = ONE / CABS1( VR( II, IS ) )
00288                CALL ZDSCAL( 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 ZGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
00296      $                        1, DCMPLX( SCALE ), VR( 1, KI ), 1 )
00297 *
00298                II = IZAMAX( N, VR( 1, KI ), 1 )
00299                REMAX = ONE / CABS1( VR( II, KI ) )
00300                CALL ZDSCAL( 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 ) = -DCONJG( 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 ZLATRS( '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 ZCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
00354 *
00355                II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
00356                REMAX = ONE / CABS1( VL( II, IS ) )
00357                CALL ZDSCAL( 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 ZGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
00365      $                        WORK( KI+1 ), 1, DCMPLX( SCALE ),
00366      $                        VL( 1, KI ), 1 )
00367 *
00368                II = IZAMAX( N, VL( 1, KI ), 1 )
00369                REMAX = ONE / CABS1( VL( II, KI ) )
00370                CALL ZDSCAL( 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 ZTREVC
00386 *
00387       END