 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, 00002 $ SEP, WORK, LWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER COMPQ, JOB 00013 INTEGER INFO, LDQ, LDT, LWORK, M, N 00014 REAL S, SEP 00015 * .. 00016 * .. Array Arguments .. 00017 LOGICAL SELECT( * ) 00018 COMPLEX Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CTRSEN reorders the Schur factorization of a complex matrix 00025 * A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in 00026 * the leading positions on the diagonal of the upper triangular matrix 00027 * T, and the leading columns of Q form an orthonormal basis of the 00028 * corresponding right invariant subspace. 00029 * 00030 * Optionally the routine computes the reciprocal condition numbers of 00031 * the cluster of eigenvalues and/or the invariant subspace. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * JOB (input) CHARACTER*1 00037 * Specifies whether condition numbers are required for the 00038 * cluster of eigenvalues (S) or the invariant subspace (SEP): 00039 * = 'N': none; 00040 * = 'E': for eigenvalues only (S); 00041 * = 'V': for invariant subspace only (SEP); 00042 * = 'B': for both eigenvalues and invariant subspace (S and 00043 * SEP). 00044 * 00045 * COMPQ (input) CHARACTER*1 00046 * = 'V': update the matrix Q of Schur vectors; 00047 * = 'N': do not update Q. 00048 * 00049 * SELECT (input) LOGICAL array, dimension (N) 00050 * SELECT specifies the eigenvalues in the selected cluster. To 00051 * select the j-th eigenvalue, SELECT(j) must be set to .TRUE.. 00052 * 00053 * N (input) INTEGER 00054 * The order of the matrix T. N >= 0. 00055 * 00056 * T (input/output) COMPLEX array, dimension (LDT,N) 00057 * On entry, the upper triangular matrix T. 00058 * On exit, T is overwritten by the reordered matrix T, with the 00059 * selected eigenvalues as the leading diagonal elements. 00060 * 00061 * LDT (input) INTEGER 00062 * The leading dimension of the array T. LDT >= max(1,N). 00063 * 00064 * Q (input/output) COMPLEX array, dimension (LDQ,N) 00065 * On entry, if COMPQ = 'V', the matrix Q of Schur vectors. 00066 * On exit, if COMPQ = 'V', Q has been postmultiplied by the 00067 * unitary transformation matrix which reorders T; the leading M 00068 * columns of Q form an orthonormal basis for the specified 00069 * invariant subspace. 00070 * If COMPQ = 'N', Q is not referenced. 00071 * 00072 * LDQ (input) INTEGER 00073 * The leading dimension of the array Q. 00074 * LDQ >= 1; and if COMPQ = 'V', LDQ >= N. 00075 * 00076 * W (output) COMPLEX array, dimension (N) 00077 * The reordered eigenvalues of T, in the same order as they 00078 * appear on the diagonal of T. 00079 * 00080 * M (output) INTEGER 00081 * The dimension of the specified invariant subspace. 00082 * 0 <= M <= N. 00083 * 00084 * S (output) REAL 00085 * If JOB = 'E' or 'B', S is a lower bound on the reciprocal 00086 * condition number for the selected cluster of eigenvalues. 00087 * S cannot underestimate the true reciprocal condition number 00088 * by more than a factor of sqrt(N). If M = 0 or N, S = 1. 00089 * If JOB = 'N' or 'V', S is not referenced. 00090 * 00091 * SEP (output) REAL 00092 * If JOB = 'V' or 'B', SEP is the estimated reciprocal 00093 * condition number of the specified invariant subspace. If 00094 * M = 0 or N, SEP = norm(T). 00095 * If JOB = 'N' or 'E', SEP is not referenced. 00096 * 00097 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00098 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00099 * 00100 * LWORK (input) INTEGER 00101 * The dimension of the array WORK. 00102 * If JOB = 'N', LWORK >= 1; 00103 * if JOB = 'E', LWORK = max(1,M*(N-M)); 00104 * if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). 00105 * 00106 * If LWORK = -1, then a workspace query is assumed; the routine 00107 * only calculates the optimal size of the WORK array, returns 00108 * this value as the first entry of the WORK array, and no error 00109 * message related to LWORK is issued by XERBLA. 00110 * 00111 * INFO (output) INTEGER 00112 * = 0: successful exit 00113 * < 0: if INFO = -i, the i-th argument had an illegal value 00114 * 00115 * Further Details 00116 * =============== 00117 * 00118 * CTRSEN first collects the selected eigenvalues by computing a unitary 00119 * transformation Z to move them to the top left corner of T. In other 00120 * words, the selected eigenvalues are the eigenvalues of T11 in: 00121 * 00122 * Z**H * T * Z = ( T11 T12 ) n1 00123 * ( 0 T22 ) n2 00124 * n1 n2 00125 * 00126 * where N = n1+n2. The first 00127 * n1 columns of Z span the specified invariant subspace of T. 00128 * 00129 * If T has been obtained from the Schur factorization of a matrix 00130 * A = Q*T*Q**H, then the reordered Schur factorization of A is given by 00131 * A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the 00132 * corresponding invariant subspace of A. 00133 * 00134 * The reciprocal condition number of the average of the eigenvalues of 00135 * T11 may be returned in S. S lies between 0 (very badly conditioned) 00136 * and 1 (very well conditioned). It is computed as follows. First we 00137 * compute R so that 00138 * 00139 * P = ( I R ) n1 00140 * ( 0 0 ) n2 00141 * n1 n2 00142 * 00143 * is the projector on the invariant subspace associated with T11. 00144 * R is the solution of the Sylvester equation: 00145 * 00146 * T11*R - R*T22 = T12. 00147 * 00148 * Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote 00149 * the two-norm of M. Then S is computed as the lower bound 00150 * 00151 * (1 + F-norm(R)**2)**(-1/2) 00152 * 00153 * on the reciprocal of 2-norm(P), the true reciprocal condition number. 00154 * S cannot underestimate 1 / 2-norm(P) by more than a factor of 00155 * sqrt(N). 00156 * 00157 * An approximate error bound for the computed average of the 00158 * eigenvalues of T11 is 00159 * 00160 * EPS * norm(T) / S 00161 * 00162 * where EPS is the machine precision. 00163 * 00164 * The reciprocal condition number of the right invariant subspace 00165 * spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. 00166 * SEP is defined as the separation of T11 and T22: 00167 * 00168 * sep( T11, T22 ) = sigma-min( C ) 00169 * 00170 * where sigma-min(C) is the smallest singular value of the 00171 * n1*n2-by-n1*n2 matrix 00172 * 00173 * C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) 00174 * 00175 * I(m) is an m by m identity matrix, and kprod denotes the Kronecker 00176 * product. We estimate sigma-min(C) by the reciprocal of an estimate of 00177 * the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) 00178 * cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). 00179 * 00180 * When SEP is small, small changes in T can cause large changes in 00181 * the invariant subspace. An approximate bound on the maximum angular 00182 * error in the computed right invariant subspace is 00183 * 00184 * EPS * norm(T) / SEP 00185 * 00186 * ===================================================================== 00187 * 00188 * .. Parameters .. 00189 REAL ZERO, ONE 00190 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00191 * .. 00192 * .. Local Scalars .. 00193 LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP 00194 INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN 00195 REAL EST, RNORM, SCALE 00196 * .. 00197 * .. Local Arrays .. 00198 INTEGER ISAVE( 3 ) 00199 REAL RWORK( 1 ) 00200 * .. 00201 * .. External Functions .. 00202 LOGICAL LSAME 00203 REAL CLANGE 00204 EXTERNAL LSAME, CLANGE 00205 * .. 00206 * .. External Subroutines .. 00207 EXTERNAL CLACN2, CLACPY, CTREXC, CTRSYL, XERBLA 00208 * .. 00209 * .. Intrinsic Functions .. 00210 INTRINSIC MAX, SQRT 00211 * .. 00212 * .. Executable Statements .. 00213 * 00214 * Decode and test the input parameters. 00215 * 00216 WANTBH = LSAME( JOB, 'B' ) 00217 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00218 WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH 00219 WANTQ = LSAME( COMPQ, 'V' ) 00220 * 00221 * Set M to the number of selected eigenvalues. 00222 * 00223 M = 0 00224 DO 10 K = 1, N 00225 IF( SELECT( K ) ) 00226 $ M = M + 1 00227 10 CONTINUE 00228 * 00229 N1 = M 00230 N2 = N - M 00231 NN = N1*N2 00232 * 00233 INFO = 0 00234 LQUERY = ( LWORK.EQ.-1 ) 00235 * 00236 IF( WANTSP ) THEN 00237 LWMIN = MAX( 1, 2*NN ) 00238 ELSE IF( LSAME( JOB, 'N' ) ) THEN 00239 LWMIN = 1 00240 ELSE IF( LSAME( JOB, 'E' ) ) THEN 00241 LWMIN = MAX( 1, NN ) 00242 END IF 00243 * 00244 IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP ) 00245 $ THEN 00246 INFO = -1 00247 ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN 00248 INFO = -2 00249 ELSE IF( N.LT.0 ) THEN 00250 INFO = -4 00251 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00252 INFO = -6 00253 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00254 INFO = -8 00255 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00256 INFO = -14 00257 END IF 00258 * 00259 IF( INFO.EQ.0 ) THEN 00260 WORK( 1 ) = LWMIN 00261 END IF 00262 * 00263 IF( INFO.NE.0 ) THEN 00264 CALL XERBLA( 'CTRSEN', -INFO ) 00265 RETURN 00266 ELSE IF( LQUERY ) THEN 00267 RETURN 00268 END IF 00269 * 00270 * Quick return if possible 00271 * 00272 IF( M.EQ.N .OR. M.EQ.0 ) THEN 00273 IF( WANTS ) 00274 $ S = ONE 00275 IF( WANTSP ) 00276 $ SEP = CLANGE( '1', N, N, T, LDT, RWORK ) 00277 GO TO 40 00278 END IF 00279 * 00280 * Collect the selected eigenvalues at the top left corner of T. 00281 * 00282 KS = 0 00283 DO 20 K = 1, N 00284 IF( SELECT( K ) ) THEN 00285 KS = KS + 1 00286 * 00287 * Swap the K-th eigenvalue to position KS. 00288 * 00289 IF( K.NE.KS ) 00290 $ CALL CTREXC( COMPQ, N, T, LDT, Q, LDQ, K, KS, IERR ) 00291 END IF 00292 20 CONTINUE 00293 * 00294 IF( WANTS ) THEN 00295 * 00296 * Solve the Sylvester equation for R: 00297 * 00298 * T11*R - R*T22 = scale*T12 00299 * 00300 CALL CLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 ) 00301 CALL CTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ), 00302 $ LDT, WORK, N1, SCALE, IERR ) 00303 * 00304 * Estimate the reciprocal of the condition number of the cluster 00305 * of eigenvalues. 00306 * 00307 RNORM = CLANGE( 'F', N1, N2, WORK, N1, RWORK ) 00308 IF( RNORM.EQ.ZERO ) THEN 00309 S = ONE 00310 ELSE 00311 S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )* 00312 $ SQRT( RNORM ) ) 00313 END IF 00314 END IF 00315 * 00316 IF( WANTSP ) THEN 00317 * 00318 * Estimate sep(T11,T22). 00319 * 00320 EST = ZERO 00321 KASE = 0 00322 30 CONTINUE 00323 CALL CLACN2( NN, WORK( NN+1 ), WORK, EST, KASE, ISAVE ) 00324 IF( KASE.NE.0 ) THEN 00325 IF( KASE.EQ.1 ) THEN 00326 * 00327 * Solve T11*R - R*T22 = scale*X. 00328 * 00329 CALL CTRSYL( 'N', 'N', -1, N1, N2, T, LDT, 00330 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE, 00331 $ IERR ) 00332 ELSE 00333 * 00334 * Solve T11**H*R - R*T22**H = scale*X. 00335 * 00336 CALL CTRSYL( 'C', 'C', -1, N1, N2, T, LDT, 00337 $ T( N1+1, N1+1 ), LDT, WORK, N1, SCALE, 00338 $ IERR ) 00339 END IF 00340 GO TO 30 00341 END IF 00342 * 00343 SEP = SCALE / EST 00344 END IF 00345 * 00346 40 CONTINUE 00347 * 00348 * Copy reordered eigenvalues to W. 00349 * 00350 DO 50 K = 1, N 00351 W( K ) = T( K, K ) 00352 50 CONTINUE 00353 * 00354 WORK( 1 ) = LWMIN 00355 * 00356 RETURN 00357 * 00358 * End of CTRSEN 00359 * 00360 END
1.7.3 
