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