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