![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTRSNA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTRSNA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsna.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsna.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsna.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTRSNA( 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 * DOUBLE PRECISION RWORK( * ), S( * ), SEP( * ) 00032 * COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00033 * $ WORK( LDWORK, * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> ZTRSNA 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*16 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*16 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 *> ZHSEIN or ZTREVC. 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*16 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 *> ZHSEIN or ZTREVC. 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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 complex16OTHERcomputational 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 ZTRSNA( 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 DOUBLE PRECISION RWORK( * ), S( * ), SEP( * ) 00264 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00265 $ WORK( LDWORK, * ) 00266 * .. 00267 * 00268 * ===================================================================== 00269 * 00270 * .. Parameters .. 00271 DOUBLE PRECISION ZERO, ONE 00272 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D0+0 ) 00273 * .. 00274 * .. Local Scalars .. 00275 LOGICAL SOMCON, WANTBH, WANTS, WANTSP 00276 CHARACTER NORMIN 00277 INTEGER I, IERR, IX, J, K, KASE, KS 00278 DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM, 00279 $ XNORM 00280 COMPLEX*16 CDUM, PROD 00281 * .. 00282 * .. Local Arrays .. 00283 INTEGER ISAVE( 3 ) 00284 COMPLEX*16 DUMMY( 1 ) 00285 * .. 00286 * .. External Functions .. 00287 LOGICAL LSAME 00288 INTEGER IZAMAX 00289 DOUBLE PRECISION DLAMCH, DZNRM2 00290 COMPLEX*16 ZDOTC 00291 EXTERNAL LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC 00292 * .. 00293 * .. External Subroutines .. 00294 EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC 00295 * .. 00296 * .. Intrinsic Functions .. 00297 INTRINSIC ABS, DBLE, DIMAG, MAX 00298 * .. 00299 * .. Statement Functions .. 00300 DOUBLE PRECISION CABS1 00301 * .. 00302 * .. Statement Function definitions .. 00303 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 00304 * .. 00305 * .. Executable Statements .. 00306 * 00307 * Decode and test the input parameters 00308 * 00309 WANTBH = LSAME( JOB, 'B' ) 00310 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00311 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 00312 * 00313 SOMCON = LSAME( HOWMNY, 'S' ) 00314 * 00315 * Set M to the number of eigenpairs for which condition numbers are 00316 * to be computed. 00317 * 00318 IF( SOMCON ) THEN 00319 M = 0 00320 DO 10 J = 1, N 00321 IF( SELECT( J ) ) 00322 $ M = M + 1 00323 10 CONTINUE 00324 ELSE 00325 M = N 00326 END IF 00327 * 00328 INFO = 0 00329 IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN 00330 INFO = -1 00331 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00332 INFO = -2 00333 ELSE IF( N.LT.0 ) THEN 00334 INFO = -4 00335 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00336 INFO = -6 00337 ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN 00338 INFO = -8 00339 ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN 00340 INFO = -10 00341 ELSE IF( MM.LT.M ) THEN 00342 INFO = -13 00343 ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN 00344 INFO = -16 00345 END IF 00346 IF( INFO.NE.0 ) THEN 00347 CALL XERBLA( 'ZTRSNA', -INFO ) 00348 RETURN 00349 END IF 00350 * 00351 * Quick return if possible 00352 * 00353 IF( N.EQ.0 ) 00354 $ RETURN 00355 * 00356 IF( N.EQ.1 ) THEN 00357 IF( SOMCON ) THEN 00358 IF( .NOT.SELECT( 1 ) ) 00359 $ RETURN 00360 END IF 00361 IF( WANTS ) 00362 $ S( 1 ) = ONE 00363 IF( WANTSP ) 00364 $ SEP( 1 ) = ABS( T( 1, 1 ) ) 00365 RETURN 00366 END IF 00367 * 00368 * Get machine constants 00369 * 00370 EPS = DLAMCH( 'P' ) 00371 SMLNUM = DLAMCH( 'S' ) / EPS 00372 BIGNUM = ONE / SMLNUM 00373 CALL DLABAD( SMLNUM, BIGNUM ) 00374 * 00375 KS = 1 00376 DO 50 K = 1, N 00377 * 00378 IF( SOMCON ) THEN 00379 IF( .NOT.SELECT( K ) ) 00380 $ GO TO 50 00381 END IF 00382 * 00383 IF( WANTS ) THEN 00384 * 00385 * Compute the reciprocal condition number of the k-th 00386 * eigenvalue. 00387 * 00388 PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 ) 00389 RNRM = DZNRM2( N, VR( 1, KS ), 1 ) 00390 LNRM = DZNRM2( N, VL( 1, KS ), 1 ) 00391 S( KS ) = ABS( PROD ) / ( RNRM*LNRM ) 00392 * 00393 END IF 00394 * 00395 IF( WANTSP ) THEN 00396 * 00397 * Estimate the reciprocal condition number of the k-th 00398 * eigenvector. 00399 * 00400 * Copy the matrix T to the array WORK and swap the k-th 00401 * diagonal element to the (1,1) position. 00402 * 00403 CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK ) 00404 CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR ) 00405 * 00406 * Form C = T22 - lambda*I in WORK(2:N,2:N). 00407 * 00408 DO 20 I = 2, N 00409 WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 ) 00410 20 CONTINUE 00411 * 00412 * Estimate a lower bound for the 1-norm of inv(C**H). The 1st 00413 * and (N+1)th columns of WORK are used to store work vectors. 00414 * 00415 SEP( KS ) = ZERO 00416 EST = ZERO 00417 KASE = 0 00418 NORMIN = 'N' 00419 30 CONTINUE 00420 CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE ) 00421 * 00422 IF( KASE.NE.0 ) THEN 00423 IF( KASE.EQ.1 ) THEN 00424 * 00425 * Solve C**H*x = scale*b 00426 * 00427 CALL ZLATRS( 'Upper', 'Conjugate transpose', 00428 $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ), 00429 $ LDWORK, WORK, SCALE, RWORK, IERR ) 00430 ELSE 00431 * 00432 * Solve C*x = scale*b 00433 * 00434 CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit', 00435 $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK, 00436 $ SCALE, RWORK, IERR ) 00437 END IF 00438 NORMIN = 'Y' 00439 IF( SCALE.NE.ONE ) THEN 00440 * 00441 * Multiply by 1/SCALE if doing so will not cause 00442 * overflow. 00443 * 00444 IX = IZAMAX( N-1, WORK, 1 ) 00445 XNORM = CABS1( WORK( IX, 1 ) ) 00446 IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) 00447 $ GO TO 40 00448 CALL ZDRSCL( N, SCALE, WORK, 1 ) 00449 END IF 00450 GO TO 30 00451 END IF 00452 * 00453 SEP( KS ) = ONE / MAX( EST, SMLNUM ) 00454 END IF 00455 * 00456 40 CONTINUE 00457 KS = KS + 1 00458 50 CONTINUE 00459 RETURN 00460 * 00461 * End of ZTRSNA 00462 * 00463 END