![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CTRSNA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CTRSNA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrsna.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrsna.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrsna.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00022 * LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, 00023 * INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER HOWMNY, JOB 00027 * INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N 00028 * .. 00029 * .. Array Arguments .. 00030 * LOGICAL SELECT( * ) 00031 * REAL RWORK( * ), S( * ), SEP( * ) 00032 * COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00033 * $ WORK( LDWORK, * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> CTRSNA estimates reciprocal condition numbers for specified 00043 *> eigenvalues and/or right eigenvectors of a complex upper triangular 00044 *> matrix T (or of any matrix Q*T*Q**H with Q unitary). 00045 *> \endverbatim 00046 * 00047 * Arguments: 00048 * ========== 00049 * 00050 *> \param[in] JOB 00051 *> \verbatim 00052 *> JOB is CHARACTER*1 00053 *> Specifies whether condition numbers are required for 00054 *> eigenvalues (S) or eigenvectors (SEP): 00055 *> = 'E': for eigenvalues only (S); 00056 *> = 'V': for eigenvectors only (SEP); 00057 *> = 'B': for both eigenvalues and eigenvectors (S and SEP). 00058 *> \endverbatim 00059 *> 00060 *> \param[in] HOWMNY 00061 *> \verbatim 00062 *> HOWMNY is CHARACTER*1 00063 *> = 'A': compute condition numbers for all eigenpairs; 00064 *> = 'S': compute condition numbers for selected eigenpairs 00065 *> specified by the array SELECT. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] SELECT 00069 *> \verbatim 00070 *> SELECT is LOGICAL array, dimension (N) 00071 *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which 00072 *> condition numbers are required. To select condition numbers 00073 *> for the j-th eigenpair, SELECT(j) must be set to .TRUE.. 00074 *> If HOWMNY = 'A', SELECT is not referenced. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] N 00078 *> \verbatim 00079 *> N is INTEGER 00080 *> The order of the matrix T. N >= 0. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] T 00084 *> \verbatim 00085 *> T is COMPLEX array, dimension (LDT,N) 00086 *> The upper triangular matrix T. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDT 00090 *> \verbatim 00091 *> LDT is INTEGER 00092 *> The leading dimension of the array T. LDT >= max(1,N). 00093 *> \endverbatim 00094 *> 00095 *> \param[in] VL 00096 *> \verbatim 00097 *> VL is COMPLEX array, dimension (LDVL,M) 00098 *> If JOB = 'E' or 'B', VL must contain left eigenvectors of T 00099 *> (or of any Q*T*Q**H with Q unitary), corresponding to the 00100 *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors 00101 *> must be stored in consecutive columns of VL, as returned by 00102 *> CHSEIN or CTREVC. 00103 *> If JOB = 'V', VL is not referenced. 00104 *> \endverbatim 00105 *> 00106 *> \param[in] LDVL 00107 *> \verbatim 00108 *> LDVL is INTEGER 00109 *> The leading dimension of the array VL. 00110 *> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] VR 00114 *> \verbatim 00115 *> VR is COMPLEX array, dimension (LDVR,M) 00116 *> If JOB = 'E' or 'B', VR must contain right eigenvectors of T 00117 *> (or of any Q*T*Q**H with Q unitary), corresponding to the 00118 *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors 00119 *> must be stored in consecutive columns of VR, as returned by 00120 *> CHSEIN or CTREVC. 00121 *> If JOB = 'V', VR is not referenced. 00122 *> \endverbatim 00123 *> 00124 *> \param[in] LDVR 00125 *> \verbatim 00126 *> LDVR is INTEGER 00127 *> The leading dimension of the array VR. 00128 *> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. 00129 *> \endverbatim 00130 *> 00131 *> \param[out] S 00132 *> \verbatim 00133 *> S is REAL array, dimension (MM) 00134 *> If JOB = 'E' or 'B', the reciprocal condition numbers of the 00135 *> selected eigenvalues, stored in consecutive elements of the 00136 *> array. Thus S(j), SEP(j), and the j-th columns of VL and VR 00137 *> all correspond to the same eigenpair (but not in general the 00138 *> j-th eigenpair, unless all eigenpairs are selected). 00139 *> If JOB = 'V', S is not referenced. 00140 *> \endverbatim 00141 *> 00142 *> \param[out] SEP 00143 *> \verbatim 00144 *> SEP is REAL array, dimension (MM) 00145 *> If JOB = 'V' or 'B', the estimated reciprocal condition 00146 *> numbers of the selected eigenvectors, stored in consecutive 00147 *> elements of the array. 00148 *> If JOB = 'E', SEP is not referenced. 00149 *> \endverbatim 00150 *> 00151 *> \param[in] MM 00152 *> \verbatim 00153 *> MM is INTEGER 00154 *> The number of elements in the arrays S (if JOB = 'E' or 'B') 00155 *> and/or SEP (if JOB = 'V' or 'B'). MM >= M. 00156 *> \endverbatim 00157 *> 00158 *> \param[out] M 00159 *> \verbatim 00160 *> M is INTEGER 00161 *> The number of elements of the arrays S and/or SEP actually 00162 *> used to store the estimated condition numbers. 00163 *> If HOWMNY = 'A', M is set to N. 00164 *> \endverbatim 00165 *> 00166 *> \param[out] WORK 00167 *> \verbatim 00168 *> WORK is COMPLEX array, dimension (LDWORK,N+6) 00169 *> If JOB = 'E', WORK is not referenced. 00170 *> \endverbatim 00171 *> 00172 *> \param[in] LDWORK 00173 *> \verbatim 00174 *> LDWORK is INTEGER 00175 *> The leading dimension of the array WORK. 00176 *> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. 00177 *> \endverbatim 00178 *> 00179 *> \param[out] RWORK 00180 *> \verbatim 00181 *> RWORK is REAL array, dimension (N) 00182 *> If JOB = 'E', RWORK is not referenced. 00183 *> \endverbatim 00184 *> 00185 *> \param[out] INFO 00186 *> \verbatim 00187 *> INFO is INTEGER 00188 *> = 0: successful exit 00189 *> < 0: if INFO = -i, the i-th argument had an illegal value 00190 *> \endverbatim 00191 * 00192 * Authors: 00193 * ======== 00194 * 00195 *> \author Univ. of Tennessee 00196 *> \author Univ. of California Berkeley 00197 *> \author Univ. of Colorado Denver 00198 *> \author NAG Ltd. 00199 * 00200 *> \date November 2011 00201 * 00202 *> \ingroup complexOTHERcomputational 00203 * 00204 *> \par Further Details: 00205 * ===================== 00206 *> 00207 *> \verbatim 00208 *> 00209 *> The reciprocal of the condition number of an eigenvalue lambda is 00210 *> defined as 00211 *> 00212 *> S(lambda) = |v**H*u| / (norm(u)*norm(v)) 00213 *> 00214 *> where u and v are the right and left eigenvectors of T corresponding 00215 *> to lambda; v**H denotes the conjugate transpose of v, and norm(u) 00216 *> denotes the Euclidean norm. These reciprocal condition numbers always 00217 *> lie between zero (very badly conditioned) and one (very well 00218 *> conditioned). If n = 1, S(lambda) is defined to be 1. 00219 *> 00220 *> An approximate error bound for a computed eigenvalue W(i) is given by 00221 *> 00222 *> EPS * norm(T) / S(i) 00223 *> 00224 *> where EPS is the machine precision. 00225 *> 00226 *> The reciprocal of the condition number of the right eigenvector u 00227 *> corresponding to lambda is defined as follows. Suppose 00228 *> 00229 *> T = ( lambda c ) 00230 *> ( 0 T22 ) 00231 *> 00232 *> Then the reciprocal condition number is 00233 *> 00234 *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I ) 00235 *> 00236 *> where sigma-min denotes the smallest singular value. We approximate 00237 *> the smallest singular value by the reciprocal of an estimate of the 00238 *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is 00239 *> defined to be abs(T(1,1)). 00240 *> 00241 *> An approximate error bound for a computed right eigenvector VR(i) 00242 *> is given by 00243 *> 00244 *> EPS * norm(T) / SEP(i) 00245 *> \endverbatim 00246 *> 00247 * ===================================================================== 00248 SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00249 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, 00250 $ INFO ) 00251 * 00252 * -- LAPACK computational routine (version 3.4.0) -- 00253 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00254 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00255 * November 2011 00256 * 00257 * .. Scalar Arguments .. 00258 CHARACTER HOWMNY, JOB 00259 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N 00260 * .. 00261 * .. Array Arguments .. 00262 LOGICAL SELECT( * ) 00263 REAL RWORK( * ), S( * ), SEP( * ) 00264 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00265 $ WORK( LDWORK, * ) 00266 * .. 00267 * 00268 * ===================================================================== 00269 * 00270 * .. Parameters .. 00271 REAL ZERO, ONE 00272 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0+0 ) 00273 * .. 00274 * .. Local Scalars .. 00275 LOGICAL SOMCON, WANTBH, WANTS, WANTSP 00276 CHARACTER NORMIN 00277 INTEGER I, IERR, IX, J, K, KASE, KS 00278 REAL BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM, 00279 $ XNORM 00280 COMPLEX CDUM, PROD 00281 * .. 00282 * .. Local Arrays .. 00283 INTEGER ISAVE( 3 ) 00284 COMPLEX DUMMY( 1 ) 00285 * .. 00286 * .. External Functions .. 00287 LOGICAL LSAME 00288 INTEGER ICAMAX 00289 REAL SCNRM2, SLAMCH 00290 COMPLEX CDOTC 00291 EXTERNAL LSAME, ICAMAX, SCNRM2, SLAMCH, CDOTC 00292 * .. 00293 * .. External Subroutines .. 00294 EXTERNAL CLACN2, CLACPY, CLATRS, CSRSCL, CTREXC, SLABAD, 00295 $ XERBLA 00296 * .. 00297 * .. Intrinsic Functions .. 00298 INTRINSIC ABS, AIMAG, MAX, REAL 00299 * .. 00300 * .. Statement Functions .. 00301 REAL CABS1 00302 * .. 00303 * .. Statement Function definitions .. 00304 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00305 * .. 00306 * .. Executable Statements .. 00307 * 00308 * Decode and test the input parameters 00309 * 00310 WANTBH = LSAME( JOB, 'B' ) 00311 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00312 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 00313 * 00314 SOMCON = LSAME( HOWMNY, 'S' ) 00315 * 00316 * Set M to the number of eigenpairs for which condition numbers are 00317 * to be computed. 00318 * 00319 IF( SOMCON ) THEN 00320 M = 0 00321 DO 10 J = 1, N 00322 IF( SELECT( J ) ) 00323 $ M = M + 1 00324 10 CONTINUE 00325 ELSE 00326 M = N 00327 END IF 00328 * 00329 INFO = 0 00330 IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN 00331 INFO = -1 00332 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00333 INFO = -2 00334 ELSE IF( N.LT.0 ) THEN 00335 INFO = -4 00336 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00337 INFO = -6 00338 ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN 00339 INFO = -8 00340 ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN 00341 INFO = -10 00342 ELSE IF( MM.LT.M ) THEN 00343 INFO = -13 00344 ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN 00345 INFO = -16 00346 END IF 00347 IF( INFO.NE.0 ) THEN 00348 CALL XERBLA( 'CTRSNA', -INFO ) 00349 RETURN 00350 END IF 00351 * 00352 * Quick return if possible 00353 * 00354 IF( N.EQ.0 ) 00355 $ RETURN 00356 * 00357 IF( N.EQ.1 ) THEN 00358 IF( SOMCON ) THEN 00359 IF( .NOT.SELECT( 1 ) ) 00360 $ RETURN 00361 END IF 00362 IF( WANTS ) 00363 $ S( 1 ) = ONE 00364 IF( WANTSP ) 00365 $ SEP( 1 ) = ABS( T( 1, 1 ) ) 00366 RETURN 00367 END IF 00368 * 00369 * Get machine constants 00370 * 00371 EPS = SLAMCH( 'P' ) 00372 SMLNUM = SLAMCH( 'S' ) / EPS 00373 BIGNUM = ONE / SMLNUM 00374 CALL SLABAD( SMLNUM, BIGNUM ) 00375 * 00376 KS = 1 00377 DO 50 K = 1, N 00378 * 00379 IF( SOMCON ) THEN 00380 IF( .NOT.SELECT( K ) ) 00381 $ GO TO 50 00382 END IF 00383 * 00384 IF( WANTS ) THEN 00385 * 00386 * Compute the reciprocal condition number of the k-th 00387 * eigenvalue. 00388 * 00389 PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 ) 00390 RNRM = SCNRM2( N, VR( 1, KS ), 1 ) 00391 LNRM = SCNRM2( N, VL( 1, KS ), 1 ) 00392 S( KS ) = ABS( PROD ) / ( RNRM*LNRM ) 00393 * 00394 END IF 00395 * 00396 IF( WANTSP ) THEN 00397 * 00398 * Estimate the reciprocal condition number of the k-th 00399 * eigenvector. 00400 * 00401 * Copy the matrix T to the array WORK and swap the k-th 00402 * diagonal element to the (1,1) position. 00403 * 00404 CALL CLACPY( 'Full', N, N, T, LDT, WORK, LDWORK ) 00405 CALL CTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR ) 00406 * 00407 * Form C = T22 - lambda*I in WORK(2:N,2:N). 00408 * 00409 DO 20 I = 2, N 00410 WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 ) 00411 20 CONTINUE 00412 * 00413 * Estimate a lower bound for the 1-norm of inv(C**H). The 1st 00414 * and (N+1)th columns of WORK are used to store work vectors. 00415 * 00416 SEP( KS ) = ZERO 00417 EST = ZERO 00418 KASE = 0 00419 NORMIN = 'N' 00420 30 CONTINUE 00421 CALL CLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE ) 00422 * 00423 IF( KASE.NE.0 ) THEN 00424 IF( KASE.EQ.1 ) THEN 00425 * 00426 * Solve C**H*x = scale*b 00427 * 00428 CALL CLATRS( 'Upper', 'Conjugate transpose', 00429 $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ), 00430 $ LDWORK, WORK, SCALE, RWORK, IERR ) 00431 ELSE 00432 * 00433 * Solve C*x = scale*b 00434 * 00435 CALL CLATRS( 'Upper', 'No transpose', 'Nonunit', 00436 $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK, 00437 $ SCALE, RWORK, IERR ) 00438 END IF 00439 NORMIN = 'Y' 00440 IF( SCALE.NE.ONE ) THEN 00441 * 00442 * Multiply by 1/SCALE if doing so will not cause 00443 * overflow. 00444 * 00445 IX = ICAMAX( N-1, WORK, 1 ) 00446 XNORM = CABS1( WORK( IX, 1 ) ) 00447 IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) 00448 $ GO TO 40 00449 CALL CSRSCL( N, SCALE, WORK, 1 ) 00450 END IF 00451 GO TO 30 00452 END IF 00453 * 00454 SEP( KS ) = ONE / MAX( EST, SMLNUM ) 00455 END IF 00456 * 00457 40 CONTINUE 00458 KS = KS + 1 00459 50 CONTINUE 00460 RETURN 00461 * 00462 * End of CTRSNA 00463 * 00464 END