00001 SUBROUTINE CTRSNA( 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 CLACN2 in place of CLACON, 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 REAL RWORK( * ), S( * ), SEP( * ) 00019 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), 00020 $ WORK( LDWORK, * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * CTRSNA 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 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 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 * CHSEIN or CTREVC. 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 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 * CHSEIN or CTREVC. 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) REAL 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) REAL 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 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) REAL 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 REAL ZERO, ONE 00165 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0+0 ) 00166 * .. 00167 * .. Local Scalars .. 00168 LOGICAL SOMCON, WANTBH, WANTS, WANTSP 00169 CHARACTER NORMIN 00170 INTEGER I, IERR, IX, J, K, KASE, KS 00171 REAL BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM, 00172 $ XNORM 00173 COMPLEX CDUM, PROD 00174 * .. 00175 * .. Local Arrays .. 00176 INTEGER ISAVE( 3 ) 00177 COMPLEX DUMMY( 1 ) 00178 * .. 00179 * .. External Functions .. 00180 LOGICAL LSAME 00181 INTEGER ICAMAX 00182 REAL SCNRM2, SLAMCH 00183 COMPLEX CDOTC 00184 EXTERNAL LSAME, ICAMAX, SCNRM2, SLAMCH, CDOTC 00185 * .. 00186 * .. External Subroutines .. 00187 EXTERNAL CLACN2, CLACPY, CLATRS, CSRSCL, CTREXC, SLABAD, 00188 $ XERBLA 00189 * .. 00190 * .. Intrinsic Functions .. 00191 INTRINSIC ABS, AIMAG, MAX, REAL 00192 * .. 00193 * .. Statement Functions .. 00194 REAL CABS1 00195 * .. 00196 * .. Statement Function definitions .. 00197 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00198 * .. 00199 * .. Executable Statements .. 00200 * 00201 * Decode and test the input parameters 00202 * 00203 WANTBH = LSAME( JOB, 'B' ) 00204 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00205 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 00206 * 00207 SOMCON = LSAME( HOWMNY, 'S' ) 00208 * 00209 * Set M to the number of eigenpairs for which condition numbers are 00210 * to be computed. 00211 * 00212 IF( SOMCON ) THEN 00213 M = 0 00214 DO 10 J = 1, N 00215 IF( SELECT( J ) ) 00216 $ M = M + 1 00217 10 CONTINUE 00218 ELSE 00219 M = N 00220 END IF 00221 * 00222 INFO = 0 00223 IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN 00224 INFO = -1 00225 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00226 INFO = -2 00227 ELSE IF( N.LT.0 ) THEN 00228 INFO = -4 00229 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00230 INFO = -6 00231 ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN 00232 INFO = -8 00233 ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN 00234 INFO = -10 00235 ELSE IF( MM.LT.M ) THEN 00236 INFO = -13 00237 ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN 00238 INFO = -16 00239 END IF 00240 IF( INFO.NE.0 ) THEN 00241 CALL XERBLA( 'CTRSNA', -INFO ) 00242 RETURN 00243 END IF 00244 * 00245 * Quick return if possible 00246 * 00247 IF( N.EQ.0 ) 00248 $ RETURN 00249 * 00250 IF( N.EQ.1 ) THEN 00251 IF( SOMCON ) THEN 00252 IF( .NOT.SELECT( 1 ) ) 00253 $ RETURN 00254 END IF 00255 IF( WANTS ) 00256 $ S( 1 ) = ONE 00257 IF( WANTSP ) 00258 $ SEP( 1 ) = ABS( T( 1, 1 ) ) 00259 RETURN 00260 END IF 00261 * 00262 * Get machine constants 00263 * 00264 EPS = SLAMCH( 'P' ) 00265 SMLNUM = SLAMCH( 'S' ) / EPS 00266 BIGNUM = ONE / SMLNUM 00267 CALL SLABAD( SMLNUM, BIGNUM ) 00268 * 00269 KS = 1 00270 DO 50 K = 1, N 00271 * 00272 IF( SOMCON ) THEN 00273 IF( .NOT.SELECT( K ) ) 00274 $ GO TO 50 00275 END IF 00276 * 00277 IF( WANTS ) THEN 00278 * 00279 * Compute the reciprocal condition number of the k-th 00280 * eigenvalue. 00281 * 00282 PROD = CDOTC( N, VR( 1, KS ), 1, VL( 1, KS ), 1 ) 00283 RNRM = SCNRM2( N, VR( 1, KS ), 1 ) 00284 LNRM = SCNRM2( N, VL( 1, KS ), 1 ) 00285 S( KS ) = ABS( PROD ) / ( RNRM*LNRM ) 00286 * 00287 END IF 00288 * 00289 IF( WANTSP ) THEN 00290 * 00291 * Estimate the reciprocal condition number of the k-th 00292 * eigenvector. 00293 * 00294 * Copy the matrix T to the array WORK and swap the k-th 00295 * diagonal element to the (1,1) position. 00296 * 00297 CALL CLACPY( 'Full', N, N, T, LDT, WORK, LDWORK ) 00298 CALL CTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, K, 1, IERR ) 00299 * 00300 * Form C = T22 - lambda*I in WORK(2:N,2:N). 00301 * 00302 DO 20 I = 2, N 00303 WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 ) 00304 20 CONTINUE 00305 * 00306 * Estimate a lower bound for the 1-norm of inv(C'). The 1st 00307 * and (N+1)th columns of WORK are used to store work vectors. 00308 * 00309 SEP( KS ) = ZERO 00310 EST = ZERO 00311 KASE = 0 00312 NORMIN = 'N' 00313 30 CONTINUE 00314 CALL CLACN2( N-1, WORK( 1, N+1 ), WORK, EST, KASE, ISAVE ) 00315 * 00316 IF( KASE.NE.0 ) THEN 00317 IF( KASE.EQ.1 ) THEN 00318 * 00319 * Solve C'*x = scale*b 00320 * 00321 CALL CLATRS( 'Upper', 'Conjugate transpose', 00322 $ 'Nonunit', NORMIN, N-1, WORK( 2, 2 ), 00323 $ LDWORK, WORK, SCALE, RWORK, IERR ) 00324 ELSE 00325 * 00326 * Solve C*x = scale*b 00327 * 00328 CALL CLATRS( 'Upper', 'No transpose', 'Nonunit', 00329 $ NORMIN, N-1, WORK( 2, 2 ), LDWORK, WORK, 00330 $ SCALE, RWORK, IERR ) 00331 END IF 00332 NORMIN = 'Y' 00333 IF( SCALE.NE.ONE ) THEN 00334 * 00335 * Multiply by 1/SCALE if doing so will not cause 00336 * overflow. 00337 * 00338 IX = ICAMAX( N-1, WORK, 1 ) 00339 XNORM = CABS1( WORK( IX, 1 ) ) 00340 IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO ) 00341 $ GO TO 40 00342 CALL CSRSCL( N, SCALE, WORK, 1 ) 00343 END IF 00344 GO TO 30 00345 END IF 00346 * 00347 SEP( KS ) = ONE / MAX( EST, SMLNUM ) 00348 END IF 00349 * 00350 40 CONTINUE 00351 KS = KS + 1 00352 50 CONTINUE 00353 RETURN 00354 * 00355 * End of CTRSNA 00356 * 00357 END
1.7.2