00001 SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 00002 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 00003 $ IWORK, 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 * .. Scalar Arguments .. 00011 CHARACTER HOWMNY, JOB 00012 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL SELECT( * ) 00016 INTEGER IWORK( * ) 00017 REAL DIF( * ), S( * ) 00018 COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 00019 $ VR( LDVR, * ), WORK( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CTGSNA estimates reciprocal condition numbers for specified 00026 * eigenvalues and/or eigenvectors of a matrix pair (A, B). 00027 * 00028 * (A, B) must be in generalized Schur canonical form, that is, A and 00029 * B are both upper triangular. 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * JOB (input) CHARACTER*1 00035 * Specifies whether condition numbers are required for 00036 * eigenvalues (S) or eigenvectors (DIF): 00037 * = 'E': for eigenvalues only (S); 00038 * = 'V': for eigenvectors only (DIF); 00039 * = 'B': for both eigenvalues and eigenvectors (S and DIF). 00040 * 00041 * HOWMNY (input) CHARACTER*1 00042 * = 'A': compute condition numbers for all eigenpairs; 00043 * = 'S': compute condition numbers for selected eigenpairs 00044 * specified by the array SELECT. 00045 * 00046 * SELECT (input) LOGICAL array, dimension (N) 00047 * If HOWMNY = 'S', SELECT specifies the eigenpairs for which 00048 * condition numbers are required. To select condition numbers 00049 * for the corresponding j-th eigenvalue and/or eigenvector, 00050 * SELECT(j) must be set to .TRUE.. 00051 * If HOWMNY = 'A', SELECT is not referenced. 00052 * 00053 * N (input) INTEGER 00054 * The order of the square matrix pair (A, B). N >= 0. 00055 * 00056 * A (input) COMPLEX array, dimension (LDA,N) 00057 * The upper triangular matrix A in the pair (A,B). 00058 * 00059 * LDA (input) INTEGER 00060 * The leading dimension of the array A. LDA >= max(1,N). 00061 * 00062 * B (input) COMPLEX array, dimension (LDB,N) 00063 * The upper triangular matrix B in the pair (A, B). 00064 * 00065 * LDB (input) INTEGER 00066 * The leading dimension of the array B. LDB >= max(1,N). 00067 * 00068 * VL (input) COMPLEX array, dimension (LDVL,M) 00069 * IF JOB = 'E' or 'B', VL must contain left eigenvectors of 00070 * (A, B), corresponding to the eigenpairs specified by HOWMNY 00071 * and SELECT. The eigenvectors must be stored in consecutive 00072 * columns of VL, as returned by CTGEVC. 00073 * If JOB = 'V', VL is not referenced. 00074 * 00075 * LDVL (input) INTEGER 00076 * The leading dimension of the array VL. LDVL >= 1; and 00077 * If JOB = 'E' or 'B', LDVL >= N. 00078 * 00079 * VR (input) COMPLEX array, dimension (LDVR,M) 00080 * IF JOB = 'E' or 'B', VR must contain right eigenvectors of 00081 * (A, B), corresponding to the eigenpairs specified by HOWMNY 00082 * and SELECT. The eigenvectors must be stored in consecutive 00083 * columns of VR, as returned by CTGEVC. 00084 * If JOB = 'V', VR is not referenced. 00085 * 00086 * LDVR (input) INTEGER 00087 * The leading dimension of the array VR. LDVR >= 1; 00088 * If JOB = 'E' or 'B', LDVR >= N. 00089 * 00090 * S (output) REAL array, dimension (MM) 00091 * If JOB = 'E' or 'B', the reciprocal condition numbers of the 00092 * selected eigenvalues, stored in consecutive elements of the 00093 * array. 00094 * If JOB = 'V', S is not referenced. 00095 * 00096 * DIF (output) REAL array, dimension (MM) 00097 * If JOB = 'V' or 'B', the estimated reciprocal condition 00098 * numbers of the selected eigenvectors, stored in consecutive 00099 * elements of the array. 00100 * If the eigenvalues cannot be reordered to compute DIF(j), 00101 * DIF(j) is set to 0; this can only occur when the true value 00102 * would be very small anyway. 00103 * For each eigenvalue/vector specified by SELECT, DIF stores 00104 * a Frobenius norm-based estimate of Difl. 00105 * If JOB = 'E', DIF is not referenced. 00106 * 00107 * MM (input) INTEGER 00108 * The number of elements in the arrays S and DIF. MM >= M. 00109 * 00110 * M (output) INTEGER 00111 * The number of elements of the arrays S and DIF used to store 00112 * the specified condition numbers; for each selected eigenvalue 00113 * one element is used. If HOWMNY = 'A', M is set to N. 00114 * 00115 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00116 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00117 * 00118 * LWORK (input) INTEGER 00119 * The dimension of the array WORK. LWORK >= max(1,N). 00120 * If JOB = 'V' or 'B', LWORK >= max(1,2*N*N). 00121 * 00122 * IWORK (workspace) INTEGER array, dimension (N+2) 00123 * If JOB = 'E', IWORK is not referenced. 00124 * 00125 * INFO (output) INTEGER 00126 * = 0: Successful exit 00127 * < 0: If INFO = -i, the i-th argument had an illegal value 00128 * 00129 * Further Details 00130 * =============== 00131 * 00132 * The reciprocal of the condition number of the i-th generalized 00133 * eigenvalue w = (a, b) is defined as 00134 * 00135 * S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v)) 00136 * 00137 * where u and v are the right and left eigenvectors of (A, B) 00138 * corresponding to w; |z| denotes the absolute value of the complex 00139 * number, and norm(u) denotes the 2-norm of the vector u. The pair 00140 * (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the 00141 * matrix pair (A, B). If both a and b equal zero, then (A,B) is 00142 * singular and S(I) = -1 is returned. 00143 * 00144 * An approximate error bound on the chordal distance between the i-th 00145 * computed generalized eigenvalue w and the corresponding exact 00146 * eigenvalue lambda is 00147 * 00148 * chord(w, lambda) <= EPS * norm(A, B) / S(I), 00149 * 00150 * where EPS is the machine precision. 00151 * 00152 * The reciprocal of the condition number of the right eigenvector u 00153 * and left eigenvector v corresponding to the generalized eigenvalue w 00154 * is defined as follows. Suppose 00155 * 00156 * (A, B) = ( a * ) ( b * ) 1 00157 * ( 0 A22 ),( 0 B22 ) n-1 00158 * 1 n-1 1 n-1 00159 * 00160 * Then the reciprocal condition number DIF(I) is 00161 * 00162 * Difl[(a, b), (A22, B22)] = sigma-min( Zl ) 00163 * 00164 * where sigma-min(Zl) denotes the smallest singular value of 00165 * 00166 * Zl = [ kron(a, In-1) -kron(1, A22) ] 00167 * [ kron(b, In-1) -kron(1, B22) ]. 00168 * 00169 * Here In-1 is the identity matrix of size n-1 and X' is the conjugate 00170 * transpose of X. kron(X, Y) is the Kronecker product between the 00171 * matrices X and Y. 00172 * 00173 * We approximate the smallest singular value of Zl with an upper 00174 * bound. This is done by CLATDF. 00175 * 00176 * An approximate error bound for a computed eigenvector VL(i) or 00177 * VR(i) is given by 00178 * 00179 * EPS * norm(A, B) / DIF(i). 00180 * 00181 * See ref. [2-3] for more details and further references. 00182 * 00183 * Based on contributions by 00184 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00185 * Umea University, S-901 87 Umea, Sweden. 00186 * 00187 * References 00188 * ========== 00189 * 00190 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00191 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00192 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00193 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00194 * 00195 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00196 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00197 * Estimation: Theory, Algorithms and Software, Report 00198 * UMINF - 94.04, Department of Computing Science, Umea University, 00199 * S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. 00200 * To appear in Numerical Algorithms, 1996. 00201 * 00202 * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 00203 * for Solving the Generalized Sylvester Equation and Estimating the 00204 * Separation between Regular Matrix Pairs, Report UMINF - 93.23, 00205 * Department of Computing Science, Umea University, S-901 87 Umea, 00206 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working 00207 * Note 75. 00208 * To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996. 00209 * 00210 * ===================================================================== 00211 * 00212 * .. Parameters .. 00213 REAL ZERO, ONE 00214 INTEGER IDIFJB 00215 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, IDIFJB = 3 ) 00216 * .. 00217 * .. Local Scalars .. 00218 LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS 00219 INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2 00220 REAL BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM 00221 COMPLEX YHAX, YHBX 00222 * .. 00223 * .. Local Arrays .. 00224 COMPLEX DUMMY( 1 ), DUMMY1( 1 ) 00225 * .. 00226 * .. External Functions .. 00227 LOGICAL LSAME 00228 REAL SCNRM2, SLAMCH, SLAPY2 00229 COMPLEX CDOTC 00230 EXTERNAL LSAME, SCNRM2, SLAMCH, SLAPY2, CDOTC 00231 * .. 00232 * .. External Subroutines .. 00233 EXTERNAL CGEMV, CLACPY, CTGEXC, CTGSYL, SLABAD, XERBLA 00234 * .. 00235 * .. Intrinsic Functions .. 00236 INTRINSIC ABS, CMPLX, MAX 00237 * .. 00238 * .. Executable Statements .. 00239 * 00240 * Decode and test the input parameters 00241 * 00242 WANTBH = LSAME( JOB, 'B' ) 00243 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00244 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH 00245 * 00246 SOMCON = LSAME( HOWMNY, 'S' ) 00247 * 00248 INFO = 0 00249 LQUERY = ( LWORK.EQ.-1 ) 00250 * 00251 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN 00252 INFO = -1 00253 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00254 INFO = -2 00255 ELSE IF( N.LT.0 ) THEN 00256 INFO = -4 00257 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00258 INFO = -6 00259 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00260 INFO = -8 00261 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN 00262 INFO = -10 00263 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN 00264 INFO = -12 00265 ELSE 00266 * 00267 * Set M to the number of eigenpairs for which condition numbers 00268 * are required, and test MM. 00269 * 00270 IF( SOMCON ) THEN 00271 M = 0 00272 DO 10 K = 1, N 00273 IF( SELECT( K ) ) 00274 $ M = M + 1 00275 10 CONTINUE 00276 ELSE 00277 M = N 00278 END IF 00279 * 00280 IF( N.EQ.0 ) THEN 00281 LWMIN = 1 00282 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN 00283 LWMIN = 2*N*N 00284 ELSE 00285 LWMIN = N 00286 END IF 00287 WORK( 1 ) = LWMIN 00288 * 00289 IF( MM.LT.M ) THEN 00290 INFO = -15 00291 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00292 INFO = -18 00293 END IF 00294 END IF 00295 * 00296 IF( INFO.NE.0 ) THEN 00297 CALL XERBLA( 'CTGSNA', -INFO ) 00298 RETURN 00299 ELSE IF( LQUERY ) THEN 00300 RETURN 00301 END IF 00302 * 00303 * Quick return if possible 00304 * 00305 IF( N.EQ.0 ) 00306 $ RETURN 00307 * 00308 * Get machine constants 00309 * 00310 EPS = SLAMCH( 'P' ) 00311 SMLNUM = SLAMCH( 'S' ) / EPS 00312 BIGNUM = ONE / SMLNUM 00313 CALL SLABAD( SMLNUM, BIGNUM ) 00314 KS = 0 00315 DO 20 K = 1, N 00316 * 00317 * Determine whether condition numbers are required for the k-th 00318 * eigenpair. 00319 * 00320 IF( SOMCON ) THEN 00321 IF( .NOT.SELECT( K ) ) 00322 $ GO TO 20 00323 END IF 00324 * 00325 KS = KS + 1 00326 * 00327 IF( WANTS ) THEN 00328 * 00329 * Compute the reciprocal condition number of the k-th 00330 * eigenvalue. 00331 * 00332 RNRM = SCNRM2( N, VR( 1, KS ), 1 ) 00333 LNRM = SCNRM2( N, VL( 1, KS ), 1 ) 00334 CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), A, LDA, 00335 $ VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 ) 00336 YHAX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 ) 00337 CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), B, LDB, 00338 $ VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 ) 00339 YHBX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 ) 00340 COND = SLAPY2( ABS( YHAX ), ABS( YHBX ) ) 00341 IF( COND.EQ.ZERO ) THEN 00342 S( KS ) = -ONE 00343 ELSE 00344 S( KS ) = COND / ( RNRM*LNRM ) 00345 END IF 00346 END IF 00347 * 00348 IF( WANTDF ) THEN 00349 IF( N.EQ.1 ) THEN 00350 DIF( KS ) = SLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) ) 00351 ELSE 00352 * 00353 * Estimate the reciprocal condition number of the k-th 00354 * eigenvectors. 00355 * 00356 * Copy the matrix (A, B) to the array WORK and move the 00357 * (k,k)th pair to the (1,1) position. 00358 * 00359 CALL CLACPY( 'Full', N, N, A, LDA, WORK, N ) 00360 CALL CLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 00361 IFST = K 00362 ILST = 1 00363 * 00364 CALL CTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), 00365 $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR ) 00366 * 00367 IF( IERR.GT.0 ) THEN 00368 * 00369 * Ill-conditioned problem - swap rejected. 00370 * 00371 DIF( KS ) = ZERO 00372 ELSE 00373 * 00374 * Reordering successful, solve generalized Sylvester 00375 * equation for R and L, 00376 * A22 * R - L * A11 = A12 00377 * B22 * R - L * B11 = B12, 00378 * and compute estimate of Difl[(A11,B11), (A22, B22)]. 00379 * 00380 N1 = 1 00381 N2 = N - N1 00382 I = N*N + 1 00383 CALL CTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ), 00384 $ N, WORK, N, WORK( N1+1 ), N, 00385 $ WORK( N*N1+N1+I ), N, WORK( I ), N, 00386 $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY, 00387 $ 1, IWORK, IERR ) 00388 END IF 00389 END IF 00390 END IF 00391 * 00392 20 CONTINUE 00393 WORK( 1 ) = LWMIN 00394 RETURN 00395 * 00396 * End of CTGSNA 00397 * 00398 END
1.7.2