|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, 00002 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK, 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 ZLACN2 in place of ZLACON, 10 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 DOUBLE PRECISION RWORK( * ), S( * ), SEP( * ) 00019 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00020 $ WORK( LDWORK, * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * ZTRSNA estimates reciprocal condition numbers for specified 00027 * eigenvalues and/or right eigenvectors of a complex upper triangular 00028 * matrix T (or of any matrix Q*T*Q**H with Q unitary). 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * JOB (input) CHARACTER*1 00034 * Specifies whether condition numbers are required for 00035 * eigenvalues (S) or eigenvectors (SEP): 00036 * = 'E': for eigenvalues only (S); 00037 * = 'V': for eigenvectors only (SEP); 00038 * = 'B': for both eigenvalues and eigenvectors (S and SEP). 00039 * 00040 * HOWMNY (input) CHARACTER*1 00041 * = 'A': compute condition numbers for all eigenpairs; 00042 * = 'S': compute condition numbers for selected eigenpairs 00043 * specified by the array SELECT. 00044 * 00045 * SELECT (input) LOGICAL array, dimension (N) 00046 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which 00047 * condition numbers are required. To select condition numbers 00048 * for the j-th eigenpair, SELECT(j) must be set to .TRUE.. 00049 * If HOWMNY = 'A', SELECT is not referenced. 00050 * 00051 * N (input) INTEGER 00052 * The order of the matrix T. N >= 0. 00053 * 00054 * T (input) COMPLEX*16 array, dimension (LDT,N) 00055 * The upper triangular matrix T. 00056 * 00057 * LDT (input) INTEGER 00058 * The leading dimension of the array T. LDT >= max(1,N). 00059 * 00060 * VL (input) COMPLEX*16 array, dimension (LDVL,M) 00061 * If JOB = 'E' or 'B', VL must contain left eigenvectors of T 00062 * (or of any Q*T*Q**H with Q unitary), corresponding to the 00063 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors 00064 * must be stored in consecutive columns of VL, as returned by 00065 * ZHSEIN or ZTREVC. 00066 * If JOB = 'V', VL is not referenced. 00067 * 00068 * LDVL (input) INTEGER 00069 * The leading dimension of the array VL. 00070 * LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N. 00071 * 00072 * VR (input) COMPLEX*16 array, dimension (LDVR,M) 00073 * If JOB = 'E' or 'B', VR must contain right eigenvectors of T 00074 * (or of any Q*T*Q**H with Q unitary), corresponding to the 00075 * eigenpairs specified by HOWMNY and SELECT. The eigenvectors 00076 * must be stored in consecutive columns of VR, as returned by 00077 * ZHSEIN or ZTREVC. 00078 * If JOB = 'V', VR is not referenced. 00079 * 00080 * LDVR (input) INTEGER 00081 * The leading dimension of the array VR. 00082 * LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N. 00083 * 00084 * S (output) DOUBLE PRECISION array, dimension (MM) 00085 * If JOB = 'E' or 'B', the reciprocal condition numbers of the 00086 * selected eigenvalues, stored in consecutive elements of the 00087 * array. Thus S(j), SEP(j), and the j-th columns of VL and VR 00088 * all correspond to the same eigenpair (but not in general the 00089 * j-th eigenpair, unless all eigenpairs are selected). 00090 * If JOB = 'V', S is not referenced. 00091 * 00092 * SEP (output) DOUBLE PRECISION array, dimension (MM) 00093 * If JOB = 'V' or 'B', the estimated reciprocal condition 00094 * numbers of the selected eigenvectors, stored in consecutive 00095 * elements of the array. 00096 * If JOB = 'E', SEP is not referenced. 00097 * 00098 * MM (input) INTEGER 00099 * The number of elements in the arrays S (if JOB = 'E' or 'B') 00100 * and/or SEP (if JOB = 'V' or 'B'). MM >= M. 00101 * 00102 * M (output) INTEGER 00103 * The number of elements of the arrays S and/or SEP actually 00104 * used to store the estimated condition numbers. 00105 * If HOWMNY = 'A', M is set to N. 00106 * 00107 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N+6) 00108 * If JOB = 'E', WORK is not referenced. 00109 * 00110 * LDWORK (input) INTEGER 00111 * The leading dimension of the array WORK. 00112 * LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N. 00113 * 00114 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00115 * If JOB = 'E', RWORK is not referenced. 00116 * 00117 * INFO (output) INTEGER 00118 * = 0: successful exit 00119 * < 0: if INFO = -i, the i-th argument had an illegal value 00120 * 00121 * Further Details 00122 * =============== 00123 * 00124 * The reciprocal of the condition number of an eigenvalue lambda is 00125 * defined as 00126 * 00127 * S(lambda) = |v'*u| / (norm(u)*norm(v)) 00128 * 00129 * where u and v are the right and left eigenvectors of T corresponding 00130 * to lambda; v' denotes the conjugate transpose of v, and norm(u) 00131 * denotes the Euclidean norm. These reciprocal condition numbers always 00132 * lie between zero (very badly conditioned) and one (very well 00133 * conditioned). If n = 1, S(lambda) is defined to be 1. 00134 * 00135 * An approximate error bound for a computed eigenvalue W(i) is given by 00136 * 00137 * EPS * norm(T) / S(i) 00138 * 00139 * where EPS is the machine precision. 00140 * 00141 * The reciprocal of the condition number of the right eigenvector u 00142 * corresponding to lambda is defined as follows. Suppose 00143 * 00144 * T = ( lambda c ) 00145 * ( 0 T22 ) 00146 * 00147 * Then the reciprocal condition number is 00148 * 00149 * SEP( lambda, T22 ) = sigma-min( T22 - lambda*I ) 00150 * 00151 * where sigma-min denotes the smallest singular value. We approximate 00152 * the smallest singular value by the reciprocal of an estimate of the 00153 * one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is 00154 * defined to be abs(T(1,1)). 00155 * 00156 * An approximate error bound for a computed right eigenvector VR(i) 00157 * is given by 00158 * 00159 * EPS * norm(T) / SEP(i) 00160 * 00161 * ===================================================================== 00162 * 00163 * .. Parameters .. 00164 DOUBLE PRECISION ZERO, ONE 00165 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D0+0 ) 00166 * .. 00167 * .. Local Scalars .. 00168 LOGICAL SOMCON, WANTBH, WANTS, WANTSP 00169 CHARACTER NORMIN 00170 INTEGER I, IERR, IX, J, K, KASE, KS 00171 DOUBLE PRECISION BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM, 00172 $ XNORM 00173 COMPLEX*16 CDUM, PROD 00174 * .. 00175 * .. Local Arrays .. 00176 INTEGER ISAVE( 3 ) 00177 COMPLEX*16 DUMMY( 1 ) 00178 * .. 00179 * .. External Functions .. 00180 LOGICAL LSAME 00181 INTEGER IZAMAX 00182 DOUBLE PRECISION DLAMCH, DZNRM2 00183 COMPLEX*16 ZDOTC 00184 EXTERNAL LSAME, IZAMAX, DLAMCH, DZNRM2, ZDOTC 00185 * .. 00186 * .. External Subroutines .. 00187 EXTERNAL XERBLA, ZDRSCL, ZLACN2, ZLACPY, ZLATRS, ZTREXC 00188 * .. 00189 * .. Intrinsic Functions .. 00190 INTRINSIC ABS, DBLE, DIMAG, MAX 00191 * .. 00192 * .. Statement Functions .. 00193 DOUBLE PRECISION CABS1 00194 * .. 00195 * .. Statement Function definitions .. 00196 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 00197 * .. 00198 * .. Executable Statements .. 00199 * 00200 * Decode and test the input parameters 00201 * 00202 WANTBH = LSAME( JOB, 'B' ) 00203 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00204 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 00205 * 00206 SOMCON = LSAME( HOWMNY, 'S' ) 00207 * 00208 * Set M to the number of eigenpairs for which condition numbers are 00209 * to be computed. 00210 * 00211 IF( SOMCON ) THEN 00212 M = 0 00213 DO 10 J = 1, N 00214 IF( SELECT( J ) ) 00215 $ M = M + 1 00216 10 CONTINUE 00217 ELSE 00218 M = N 00219 END IF 00220 * 00221 INFO = 0 00222 IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN 00223 INFO = -1 00224 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00225 INFO = -2 00226 ELSE IF( N.LT.0 ) THEN 00227 INFO = -4 00228 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00229 INFO = -6 00230 ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN 00231 INFO = -8 00232 ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN 00233 INFO = -10 00234 ELSE IF( MM.LT.M ) THEN 00235 INFO = -13 00236 ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN 00237 INFO = -16 00238 END IF 00239 IF( INFO.NE.0 ) THEN 00240 CALL XERBLA( 'ZTRSNA', -INFO ) 00241 RETURN 00242 END IF 00243 * 00244 * Quick return if possible 00245 * 00246 IF( N.EQ.0 ) 00247 $ RETURN 00248 * 00249 IF( N.EQ.1 ) THEN 00250 IF( SOMCON ) THEN 00251 IF( .NOT.SELECT( 1 ) ) 00252 $ RETURN 00253 END IF 00254 IF( WANTS ) 00255 $ S( 1 ) = ONE 00256 IF( WANTSP ) 00257 $ SEP( 1 ) = ABS( T( 1, 1 ) ) 00258 RETURN 00259 END IF 00260 * 00261 * Get machine constants 00262 * 00263 EPS = DLAMCH( 'P' ) 00264 SMLNUM = DLAMCH( 'S' ) / EPS 00265 BIGNUM = ONE / SMLNUM 00266 CALL DLABAD( SMLNUM, BIGNUM ) 00267 * 00268 KS = 1 00269 DO 50 K = 1, N 00270 * 00271 IF( SOMCON ) THEN 00272 IF( .NOT.SELECT( K ) ) 00273 $ GO TO 50 00274 END IF 00275 * 00276 IF( WANTS ) THEN 00277 * 00278 * Compute the reciprocal condition number of the k-th 00279 * eigenvalue. 00280 * 00281 PROD = ZDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 ) 00282 RNRM = DZNRM2( N, VR( 1, KS ), 1 ) 00283 LNRM = DZNRM2( N, VL( 1, KS ), 1 ) 00284 S( KS ) = ABS( PROD ) / ( RNRM*LNRM ) 00285 * 00286 END IF 00287 * 00288 IF( WANTSP ) THEN 00289 * 00290 * Estimate the reciprocal condition number of the k-th 00291 * eigenvector. 00292 * 00293 * Copy the matrix T to the array WORK and swap the k-th 00294 * diagonal element to the (1,1) position. 00295 * 00296 CALL ZLACPY( 'Full', N, N, T, LDT, WORK, LDWORK ) 00297 CALL ZTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR ) 00298 * 00299 * Form C = T22 - lambda*I in WORK(2:N,2:N). 00300 * 00301 DO 20 I = 2, N 00302 WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 ) 00303 20 CONTINUE 00304 * 00305 * Estimate a lower bound for the 1-norm of inv(C'). The 1st 00306 * and (N+1)th columns of WORK are used to store work vectors. 00307 * 00308 SEP( KS ) = ZERO 00309 EST = ZERO 00310 KASE = 0 00311 NORMIN = 'N' 00312 30 CONTINUE 00313 CALL ZLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE ) 00314 * 00315 IF( KASE.NE.0 ) THEN 00316 IF( KASE.EQ.1 ) THEN 00317 * 00318 * Solve C'*x = scale*b 00319 * 00320 CALL ZLATRS( 'Upper', 'Conjugate transpose', 00321 $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ), 00322 $ LDWORK, WORK, SCALE, RWORK, IERR ) 00323 ELSE 00324 * 00325 * Solve C*x = scale*b 00326 * 00327 CALL ZLATRS( 'Upper', 'No transpose', 'Nonunit', 00328 $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK, 00329 $ SCALE, RWORK, IERR ) 00330 END IF 00331 NORMIN = 'Y' 00332 IF( SCALE.NE.ONE ) THEN 00333 * 00334 * Multiply by 1/SCALE if doing so will not cause 00335 * overflow. 00336 * 00337 IX = IZAMAX( N-1, WORK, 1 ) 00338 XNORM = CABS1( WORK( IX, 1 ) ) 00339 IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) 00340 $ GO TO 40 00341 CALL ZDRSCL( N, SCALE, WORK, 1 ) 00342 END IF 00343 GO TO 30 00344 END IF 00345 * 00346 SEP( KS ) = ONE / MAX( EST, SMLNUM ) 00347 END IF 00348 * 00349 40 CONTINUE 00350 KS = KS + 1 00351 50 CONTINUE 00352 RETURN 00353 * 00354 * End of ZTRSNA 00355 * 00356 END
1.7.3 
