|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE STRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00002 $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, 00003 $ INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH. 00011 * 00012 * .. Scalar Arguments .. 00013 CHARACTER HOWMNY, JOB 00014 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N 00015 * .. 00016 * .. Array Arguments .. 00017 LOGICAL SELECT( * ) 00018 INTEGER IWORK( * ) 00019 REAL S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ), 00020 $ VR( LDVR, * ), WORK( LDWORK, * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * STRSNA estimates reciprocal condition numbers for specified 00027 * eigenvalues and/or right eigenvectors of a real upper 00028 * quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q 00029 * orthogonal). 00030 * 00031 * T must be in Schur canonical form (as returned by SHSEQR), that is, 00032 * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 00033 * 2-by-2 diagonal block has its diagonal elements equal and its 00034 * off-diagonal elements of opposite sign. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * JOB (input) CHARACTER*1 00040 * Specifies whether condition numbers are required for 00041 * eigenvalues (S) or eigenvectors (SEP): 00042 * = 'E': for eigenvalues only (S); 00043 * = 'V': for eigenvectors only (SEP); 00044 * = 'B': for both eigenvalues and eigenvectors (S and SEP). 00045 * 00046 * HOWMNY (input) CHARACTER*1 00047 * = 'A': compute condition numbers for all eigenpairs; 00048 * = 'S': compute condition numbers for selected eigenpairs 00049 * specified by the array SELECT. 00050 * 00051 * SELECT (input) LOGICAL array, dimension (N) 00052 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which 00053 * condition numbers are required. To select condition numbers 00054 * for the eigenpair corresponding to a real eigenvalue w(j), 00055 * SELECT(j) must be set to .TRUE.. To select condition numbers 00056 * corresponding to a complex conjugate pair of eigenvalues w(j) 00057 * and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be 00058 * set to .TRUE.. 00059 * If HOWMNY = 'A', SELECT is not referenced. 00060 * 00061 * N (input) INTEGER 00062 * The order of the matrix T. N >= 0. 00063 * 00064 * T (input) REAL array, dimension (LDT,N) 00065 * The upper quasi-triangular matrix T, in Schur canonical form. 00066 * 00067 * LDT (input) INTEGER 00068 * The leading dimension of the array T. LDT >= max(1,N). 00069 * 00070 * VL (input) REAL array, dimension (LDVL,M) 00071 * If JOB = 'E' or 'B', VL must contain left eigenvectors of T 00072 * (or of any Q*T*Q**T with Q orthogonal), corresponding to the 00073 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors 00074 * must be stored in consecutive columns of VL, as returned by 00075 * SHSEIN or STREVC. 00076 * If JOB = 'V', VL is not referenced. 00077 * 00078 * LDVL (input) INTEGER 00079 * The leading dimension of the array VL. 00080 * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. 00081 * 00082 * VR (input) REAL array, dimension (LDVR,M) 00083 * If JOB = 'E' or 'B', VR must contain right eigenvectors of T 00084 * (or of any Q*T*Q**T with Q orthogonal), corresponding to the 00085 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors 00086 * must be stored in consecutive columns of VR, as returned by 00087 * SHSEIN or STREVC. 00088 * If JOB = 'V', VR is not referenced. 00089 * 00090 * LDVR (input) INTEGER 00091 * The leading dimension of the array VR. 00092 * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. 00093 * 00094 * S (output) REAL array, dimension (MM) 00095 * If JOB = 'E' or 'B', the reciprocal condition numbers of the 00096 * selected eigenvalues, stored in consecutive elements of the 00097 * array. For a complex conjugate pair of eigenvalues two 00098 * consecutive elements of S are set to the same value. Thus 00099 * S(j), SEP(j), and the j-th columns of VL and VR all 00100 * correspond to the same eigenpair (but not in general the 00101 * j-th eigenpair, unless all eigenpairs are selected). 00102 * If JOB = 'V', S is not referenced. 00103 * 00104 * SEP (output) REAL array, dimension (MM) 00105 * If JOB = 'V' or 'B', the estimated reciprocal condition 00106 * numbers of the selected eigenvectors, stored in consecutive 00107 * elements of the array. For a complex eigenvector two 00108 * consecutive elements of SEP are set to the same value. If 00109 * the eigenvalues cannot be reordered to compute SEP(j), SEP(j) 00110 * is set to 0; this can only occur when the true value would be 00111 * very small anyway. 00112 * If JOB = 'E', SEP is not referenced. 00113 * 00114 * MM (input) INTEGER 00115 * The number of elements in the arrays S (if JOB = 'E' or 'B') 00116 * and/or SEP (if JOB = 'V' or 'B'). MM >= M. 00117 * 00118 * M (output) INTEGER 00119 * The number of elements of the arrays S and/or SEP actually 00120 * used to store the estimated condition numbers. 00121 * If HOWMNY = 'A', M is set to N. 00122 * 00123 * WORK (workspace) REAL array, dimension (LDWORK,N+6) 00124 * If JOB = 'E', WORK is not referenced. 00125 * 00126 * LDWORK (input) INTEGER 00127 * The leading dimension of the array WORK. 00128 * LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. 00129 * 00130 * IWORK (workspace) INTEGER array, dimension (2*(N-1)) 00131 * If JOB = 'E', IWORK is not referenced. 00132 * 00133 * INFO (output) INTEGER 00134 * = 0: successful exit 00135 * < 0: if INFO = -i, the i-th argument had an illegal value 00136 * 00137 * Further Details 00138 * =============== 00139 * 00140 * The reciprocal of the condition number of an eigenvalue lambda is 00141 * defined as 00142 * 00143 * S(lambda) = |v'*u| / (norm(u)*norm(v)) 00144 * 00145 * where u and v are the right and left eigenvectors of T corresponding 00146 * to lambda; v' denotes the conjugate-transpose of v, and norm(u) 00147 * denotes the Euclidean norm. These reciprocal condition numbers always 00148 * lie between zero (very badly conditioned) and one (very well 00149 * conditioned). If n = 1, S(lambda) is defined to be 1. 00150 * 00151 * An approximate error bound for a computed eigenvalue W(i) is given by 00152 * 00153 * EPS * norm(T) / S(i) 00154 * 00155 * where EPS is the machine precision. 00156 * 00157 * The reciprocal of the condition number of the right eigenvector u 00158 * corresponding to lambda is defined as follows. Suppose 00159 * 00160 * T = ( lambda c ) 00161 * ( 0 T22 ) 00162 * 00163 * Then the reciprocal condition number is 00164 * 00165 * SEP( lambda, T22 ) = sigma-min( T22 - lambda*I ) 00166 * 00167 * where sigma-min denotes the smallest singular value. We approximate 00168 * the smallest singular value by the reciprocal of an estimate of the 00169 * one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is 00170 * defined to be abs(T(1,1)). 00171 * 00172 * An approximate error bound for a computed right eigenvector VR(i) 00173 * is given by 00174 * 00175 * EPS * norm(T) / SEP(i) 00176 * 00177 * ===================================================================== 00178 * 00179 * .. Parameters .. 00180 REAL ZERO, ONE, TWO 00181 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 00182 * .. 00183 * .. Local Scalars .. 00184 LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP 00185 INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN 00186 REAL BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM, 00187 $ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN 00188 * .. 00189 * .. Local Arrays .. 00190 INTEGER ISAVE( 3 ) 00191 REAL DUMMY( 1 ) 00192 * .. 00193 * .. External Functions .. 00194 LOGICAL LSAME 00195 REAL SDOT, SLAMCH, SLAPY2, SNRM2 00196 EXTERNAL LSAME, SDOT, SLAMCH, SLAPY2, SNRM2 00197 * .. 00198 * .. External Subroutines .. 00199 EXTERNAL SLABAD, SLACN2, SLACPY, SLAQTR, STREXC, XERBLA 00200 * .. 00201 * .. Intrinsic Functions .. 00202 INTRINSIC ABS, MAX, SQRT 00203 * .. 00204 * .. Executable Statements .. 00205 * 00206 * Decode and test the input parameters 00207 * 00208 WANTBH = LSAME( JOB, 'B' ) 00209 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00210 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 00211 * 00212 SOMCON = LSAME( HOWMNY, 'S' ) 00213 * 00214 INFO = 0 00215 IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN 00216 INFO = -1 00217 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00218 INFO = -2 00219 ELSE IF( N.LT.0 ) THEN 00220 INFO = -4 00221 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00222 INFO = -6 00223 ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN 00224 INFO = -8 00225 ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN 00226 INFO = -10 00227 ELSE 00228 * 00229 * Set M to the number of eigenpairs for which condition numbers 00230 * are required, and test MM. 00231 * 00232 IF( SOMCON ) THEN 00233 M = 0 00234 PAIR = .FALSE. 00235 DO 10 K = 1, N 00236 IF( PAIR ) THEN 00237 PAIR = .FALSE. 00238 ELSE 00239 IF( K.LT.N ) THEN 00240 IF( T( K+1, K ).EQ.ZERO ) THEN 00241 IF( SELECT( K ) ) 00242 $ M = M + 1 00243 ELSE 00244 PAIR = .TRUE. 00245 IF( SELECT( K ) .OR. SELECT( K+1 ) ) 00246 $ M = M + 2 00247 END IF 00248 ELSE 00249 IF( SELECT( N ) ) 00250 $ M = M + 1 00251 END IF 00252 END IF 00253 10 CONTINUE 00254 ELSE 00255 M = N 00256 END IF 00257 * 00258 IF( MM.LT.M ) THEN 00259 INFO = -13 00260 ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN 00261 INFO = -16 00262 END IF 00263 END IF 00264 IF( INFO.NE.0 ) THEN 00265 CALL XERBLA( 'STRSNA', -INFO ) 00266 RETURN 00267 END IF 00268 * 00269 * Quick return if possible 00270 * 00271 IF( N.EQ.0 ) 00272 $ RETURN 00273 * 00274 IF( N.EQ.1 ) THEN 00275 IF( SOMCON ) THEN 00276 IF( .NOT.SELECT( 1 ) ) 00277 $ RETURN 00278 END IF 00279 IF( WANTS ) 00280 $ S( 1 ) = ONE 00281 IF( WANTSP ) 00282 $ SEP( 1 ) = ABS( T( 1, 1 ) ) 00283 RETURN 00284 END IF 00285 * 00286 * Get machine constants 00287 * 00288 EPS = SLAMCH( 'P' ) 00289 SMLNUM = SLAMCH( 'S' ) / EPS 00290 BIGNUM = ONE / SMLNUM 00291 CALL SLABAD( SMLNUM, BIGNUM ) 00292 * 00293 KS = 0 00294 PAIR = .FALSE. 00295 DO 60 K = 1, N 00296 * 00297 * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block. 00298 * 00299 IF( PAIR ) THEN 00300 PAIR = .FALSE. 00301 GO TO 60 00302 ELSE 00303 IF( K.LT.N ) 00304 $ PAIR = T( K+1, K ).NE.ZERO 00305 END IF 00306 * 00307 * Determine whether condition numbers are required for the k-th 00308 * eigenpair. 00309 * 00310 IF( SOMCON ) THEN 00311 IF( PAIR ) THEN 00312 IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) ) 00313 $ GO TO 60 00314 ELSE 00315 IF( .NOT.SELECT( K ) ) 00316 $ GO TO 60 00317 END IF 00318 END IF 00319 * 00320 KS = KS + 1 00321 * 00322 IF( WANTS ) THEN 00323 * 00324 * Compute the reciprocal condition number of the k-th 00325 * eigenvalue. 00326 * 00327 IF( .NOT.PAIR ) THEN 00328 * 00329 * Real eigenvalue. 00330 * 00331 PROD = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 ) 00332 RNRM = SNRM2( N, VR( 1, KS ), 1 ) 00333 LNRM = SNRM2( N, VL( 1, KS ), 1 ) 00334 S( KS ) = ABS( PROD ) / ( RNRM*LNRM ) 00335 ELSE 00336 * 00337 * Complex eigenvalue. 00338 * 00339 PROD1 = SDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 ) 00340 PROD1 = PROD1 + SDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ), 00341 $ 1 ) 00342 PROD2 = SDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 ) 00343 PROD2 = PROD2 - SDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ), 00344 $ 1 ) 00345 RNRM = SLAPY2( SNRM2( N, VR( 1, KS ), 1 ), 00346 $ SNRM2( N, VR( 1, KS+1 ), 1 ) ) 00347 LNRM = SLAPY2( SNRM2( N, VL( 1, KS ), 1 ), 00348 $ SNRM2( N, VL( 1, KS+1 ), 1 ) ) 00349 COND = SLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM ) 00350 S( KS ) = COND 00351 S( KS+1 ) = COND 00352 END IF 00353 END IF 00354 * 00355 IF( WANTSP ) THEN 00356 * 00357 * Estimate the reciprocal condition number of the k-th 00358 * eigenvector. 00359 * 00360 * Copy the matrix T to the array WORK and swap the diagonal 00361 * block beginning at T(k,k) to the (1,1) position. 00362 * 00363 CALL SLACPY( 'Full', N, N, T, LDT, WORK, LDWORK ) 00364 IFST = K 00365 ILST = 1 00366 CALL STREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST, 00367 $ WORK( 1, N+1 ), IERR ) 00368 * 00369 IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN 00370 * 00371 * Could not swap because blocks not well separated 00372 * 00373 SCALE = ONE 00374 EST = BIGNUM 00375 ELSE 00376 * 00377 * Reordering successful 00378 * 00379 IF( WORK( 2, 1 ).EQ.ZERO ) THEN 00380 * 00381 * Form C = T22 - lambda*I in WORK(2:N,2:N). 00382 * 00383 DO 20 I = 2, N 00384 WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 ) 00385 20 CONTINUE 00386 N2 = 1 00387 NN = N - 1 00388 ELSE 00389 * 00390 * Triangularize the 2 by 2 block by unitary 00391 * transformation U = [ cs i*ss ] 00392 * [ i*ss cs ]. 00393 * such that the (1,1) position of WORK is complex 00394 * eigenvalue lambda with positive imaginary part. (2,2) 00395 * position of WORK is the complex eigenvalue lambda 00396 * with negative imaginary part. 00397 * 00398 MU = SQRT( ABS( WORK( 1, 2 ) ) )* 00399 $ SQRT( ABS( WORK( 2, 1 ) ) ) 00400 DELTA = SLAPY2( MU, WORK( 2, 1 ) ) 00401 CS = MU / DELTA 00402 SN = -WORK( 2, 1 ) / DELTA 00403 * 00404 * Form 00405 * 00406 * C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ] 00407 * [ mu ] 00408 * [ .. ] 00409 * [ .. ] 00410 * [ mu ] 00411 * where C' is conjugate transpose of complex matrix C, 00412 * and RWORK is stored starting in the N+1-st column of 00413 * WORK. 00414 * 00415 DO 30 J = 3, N 00416 WORK( 2, J ) = CS*WORK( 2, J ) 00417 WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 ) 00418 30 CONTINUE 00419 WORK( 2, 2 ) = ZERO 00420 * 00421 WORK( 1, N+1 ) = TWO*MU 00422 DO 40 I = 2, N - 1 00423 WORK( I, N+1 ) = SN*WORK( 1, I+1 ) 00424 40 CONTINUE 00425 N2 = 2 00426 NN = 2*( N-1 ) 00427 END IF 00428 * 00429 * Estimate norm(inv(C')) 00430 * 00431 EST = ZERO 00432 KASE = 0 00433 50 CONTINUE 00434 CALL SLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK, 00435 $ EST, KASE, ISAVE ) 00436 IF( KASE.NE.0 ) THEN 00437 IF( KASE.EQ.1 ) THEN 00438 IF( N2.EQ.1 ) THEN 00439 * 00440 * Real eigenvalue: solve C'*x = scale*c. 00441 * 00442 CALL SLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ), 00443 $ LDWORK, DUMMY, DUMM, SCALE, 00444 $ WORK( 1, N+4 ), WORK( 1, N+6 ), 00445 $ IERR ) 00446 ELSE 00447 * 00448 * Complex eigenvalue: solve 00449 * C'*(p+iq) = scale*(c+id) in real arithmetic. 00450 * 00451 CALL SLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ), 00452 $ LDWORK, WORK( 1, N+1 ), MU, SCALE, 00453 $ WORK( 1, N+4 ), WORK( 1, N+6 ), 00454 $ IERR ) 00455 END IF 00456 ELSE 00457 IF( N2.EQ.1 ) THEN 00458 * 00459 * Real eigenvalue: solve C*x = scale*c. 00460 * 00461 CALL SLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ), 00462 $ LDWORK, DUMMY, DUMM, SCALE, 00463 $ WORK( 1, N+4 ), WORK( 1, N+6 ), 00464 $ IERR ) 00465 ELSE 00466 * 00467 * Complex eigenvalue: solve 00468 * C*(p+iq) = scale*(c+id) in real arithmetic. 00469 * 00470 CALL SLAQTR( .FALSE., .FALSE., N-1, 00471 $ WORK( 2, 2 ), LDWORK, 00472 $ WORK( 1, N+1 ), MU, SCALE, 00473 $ WORK( 1, N+4 ), WORK( 1, N+6 ), 00474 $ IERR ) 00475 * 00476 END IF 00477 END IF 00478 * 00479 GO TO 50 00480 END IF 00481 END IF 00482 * 00483 SEP( KS ) = SCALE / MAX( EST, SMLNUM ) 00484 IF( PAIR ) 00485 $ SEP( KS+1 ) = SEP( KS ) 00486 END IF 00487 * 00488 IF( PAIR ) 00489 $ KS = KS + 1 00490 * 00491 60 CONTINUE 00492 RETURN 00493 * 00494 * End of STRSNA 00495 * 00496 END
1.7.3 
