![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CTGSNA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CTGSNA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsna.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsna.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsna.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 00022 * LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 00023 * IWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER HOWMNY, JOB 00027 * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 00028 * .. 00029 * .. Array Arguments .. 00030 * LOGICAL SELECT( * ) 00031 * INTEGER IWORK( * ) 00032 * REAL DIF( * ), S( * ) 00033 * COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 00034 * $ VR( LDVR, * ), WORK( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> CTGSNA estimates reciprocal condition numbers for specified 00044 *> eigenvalues and/or eigenvectors of a matrix pair (A, B). 00045 *> 00046 *> (A, B) must be in generalized Schur canonical form, that is, A and 00047 *> B are both upper triangular. 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 00057 *> eigenvalues (S) or eigenvectors (DIF): 00058 *> = 'E': for eigenvalues only (S); 00059 *> = 'V': for eigenvectors only (DIF); 00060 *> = 'B': for both eigenvalues and eigenvectors (S and DIF). 00061 *> \endverbatim 00062 *> 00063 *> \param[in] HOWMNY 00064 *> \verbatim 00065 *> HOWMNY is CHARACTER*1 00066 *> = 'A': compute condition numbers for all eigenpairs; 00067 *> = 'S': compute condition numbers for selected eigenpairs 00068 *> specified by the array SELECT. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] SELECT 00072 *> \verbatim 00073 *> SELECT is LOGICAL array, dimension (N) 00074 *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which 00075 *> condition numbers are required. To select condition numbers 00076 *> for the corresponding j-th eigenvalue and/or eigenvector, 00077 *> SELECT(j) must be set to .TRUE.. 00078 *> If HOWMNY = 'A', SELECT is not referenced. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] N 00082 *> \verbatim 00083 *> N is INTEGER 00084 *> The order of the square matrix pair (A, B). N >= 0. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] A 00088 *> \verbatim 00089 *> A is COMPLEX array, dimension (LDA,N) 00090 *> The upper triangular matrix A in the pair (A,B). 00091 *> \endverbatim 00092 *> 00093 *> \param[in] LDA 00094 *> \verbatim 00095 *> LDA is INTEGER 00096 *> The leading dimension of the array A. LDA >= max(1,N). 00097 *> \endverbatim 00098 *> 00099 *> \param[in] B 00100 *> \verbatim 00101 *> B is COMPLEX array, dimension (LDB,N) 00102 *> The upper triangular matrix B in the pair (A, B). 00103 *> \endverbatim 00104 *> 00105 *> \param[in] LDB 00106 *> \verbatim 00107 *> LDB is INTEGER 00108 *> The leading dimension of the array B. LDB >= max(1,N). 00109 *> \endverbatim 00110 *> 00111 *> \param[in] VL 00112 *> \verbatim 00113 *> VL is COMPLEX array, dimension (LDVL,M) 00114 *> IF JOB = 'E' or 'B', VL must contain left eigenvectors of 00115 *> (A, B), corresponding to the eigenpairs specified by HOWMNY 00116 *> and SELECT. The eigenvectors must be stored in consecutive 00117 *> columns of VL, as returned by CTGEVC. 00118 *> If JOB = 'V', VL is not referenced. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] LDVL 00122 *> \verbatim 00123 *> LDVL is INTEGER 00124 *> The leading dimension of the array VL. LDVL >= 1; and 00125 *> If JOB = 'E' or 'B', LDVL >= N. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] VR 00129 *> \verbatim 00130 *> VR is COMPLEX array, dimension (LDVR,M) 00131 *> IF JOB = 'E' or 'B', VR must contain right eigenvectors of 00132 *> (A, B), corresponding to the eigenpairs specified by HOWMNY 00133 *> and SELECT. The eigenvectors must be stored in consecutive 00134 *> columns of VR, as returned by CTGEVC. 00135 *> If JOB = 'V', VR is not referenced. 00136 *> \endverbatim 00137 *> 00138 *> \param[in] LDVR 00139 *> \verbatim 00140 *> LDVR is INTEGER 00141 *> The leading dimension of the array VR. LDVR >= 1; 00142 *> If JOB = 'E' or 'B', LDVR >= N. 00143 *> \endverbatim 00144 *> 00145 *> \param[out] S 00146 *> \verbatim 00147 *> S is REAL array, dimension (MM) 00148 *> If JOB = 'E' or 'B', the reciprocal condition numbers of the 00149 *> selected eigenvalues, stored in consecutive elements of the 00150 *> array. 00151 *> If JOB = 'V', S is not referenced. 00152 *> \endverbatim 00153 *> 00154 *> \param[out] DIF 00155 *> \verbatim 00156 *> DIF is REAL array, dimension (MM) 00157 *> If JOB = 'V' or 'B', the estimated reciprocal condition 00158 *> numbers of the selected eigenvectors, stored in consecutive 00159 *> elements of the array. 00160 *> If the eigenvalues cannot be reordered to compute DIF(j), 00161 *> DIF(j) is set to 0; this can only occur when the true value 00162 *> would be very small anyway. 00163 *> For each eigenvalue/vector specified by SELECT, DIF stores 00164 *> a Frobenius norm-based estimate of Difl. 00165 *> If JOB = 'E', DIF is not referenced. 00166 *> \endverbatim 00167 *> 00168 *> \param[in] MM 00169 *> \verbatim 00170 *> MM is INTEGER 00171 *> The number of elements in the arrays S and DIF. MM >= M. 00172 *> \endverbatim 00173 *> 00174 *> \param[out] M 00175 *> \verbatim 00176 *> M is INTEGER 00177 *> The number of elements of the arrays S and DIF used to store 00178 *> the specified condition numbers; for each selected eigenvalue 00179 *> one element is used. If HOWMNY = 'A', M is set to N. 00180 *> \endverbatim 00181 *> 00182 *> \param[out] WORK 00183 *> \verbatim 00184 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00185 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00186 *> \endverbatim 00187 *> 00188 *> \param[in] LWORK 00189 *> \verbatim 00190 *> LWORK is INTEGER 00191 *> The dimension of the array WORK. LWORK >= max(1,N). 00192 *> If JOB = 'V' or 'B', LWORK >= max(1,2*N*N). 00193 *> \endverbatim 00194 *> 00195 *> \param[out] IWORK 00196 *> \verbatim 00197 *> IWORK is INTEGER array, dimension (N+2) 00198 *> If JOB = 'E', IWORK is not referenced. 00199 *> \endverbatim 00200 *> 00201 *> \param[out] INFO 00202 *> \verbatim 00203 *> INFO is INTEGER 00204 *> = 0: Successful exit 00205 *> < 0: If INFO = -i, the i-th argument had an illegal value 00206 *> \endverbatim 00207 * 00208 * Authors: 00209 * ======== 00210 * 00211 *> \author Univ. of Tennessee 00212 *> \author Univ. of California Berkeley 00213 *> \author Univ. of Colorado Denver 00214 *> \author NAG Ltd. 00215 * 00216 *> \date November 2011 00217 * 00218 *> \ingroup complexOTHERcomputational 00219 * 00220 *> \par Further Details: 00221 * ===================== 00222 *> 00223 *> \verbatim 00224 *> 00225 *> The reciprocal of the condition number of the i-th generalized 00226 *> eigenvalue w = (a, b) is defined as 00227 *> 00228 *> S(I) = (|v**HAu|**2 + |v**HBu|**2)**(1/2) / (norm(u)*norm(v)) 00229 *> 00230 *> where u and v are the right and left eigenvectors of (A, B) 00231 *> corresponding to w; |z| denotes the absolute value of the complex 00232 *> number, and norm(u) denotes the 2-norm of the vector u. The pair 00233 *> (a, b) corresponds to an eigenvalue w = a/b (= v**HAu/v**HBu) of the 00234 *> matrix pair (A, B). If both a and b equal zero, then (A,B) is 00235 *> singular and S(I) = -1 is returned. 00236 *> 00237 *> An approximate error bound on the chordal distance between the i-th 00238 *> computed generalized eigenvalue w and the corresponding exact 00239 *> eigenvalue lambda is 00240 *> 00241 *> chord(w, lambda) <= EPS * norm(A, B) / S(I), 00242 *> 00243 *> where EPS is the machine precision. 00244 *> 00245 *> The reciprocal of the condition number of the right eigenvector u 00246 *> and left eigenvector v corresponding to the generalized eigenvalue w 00247 *> is defined as follows. Suppose 00248 *> 00249 *> (A, B) = ( a * ) ( b * ) 1 00250 *> ( 0 A22 ),( 0 B22 ) n-1 00251 *> 1 n-1 1 n-1 00252 *> 00253 *> Then the reciprocal condition number DIF(I) is 00254 *> 00255 *> Difl[(a, b), (A22, B22)] = sigma-min( Zl ) 00256 *> 00257 *> where sigma-min(Zl) denotes the smallest singular value of 00258 *> 00259 *> Zl = [ kron(a, In-1) -kron(1, A22) ] 00260 *> [ kron(b, In-1) -kron(1, B22) ]. 00261 *> 00262 *> Here In-1 is the identity matrix of size n-1 and X**H is the conjugate 00263 *> transpose of X. kron(X, Y) is the Kronecker product between the 00264 *> matrices X and Y. 00265 *> 00266 *> We approximate the smallest singular value of Zl with an upper 00267 *> bound. This is done by CLATDF. 00268 *> 00269 *> An approximate error bound for a computed eigenvector VL(i) or 00270 *> VR(i) is given by 00271 *> 00272 *> EPS * norm(A, B) / DIF(i). 00273 *> 00274 *> See ref. [2-3] for more details and further references. 00275 *> \endverbatim 00276 * 00277 *> \par Contributors: 00278 * ================== 00279 *> 00280 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00281 *> Umea University, S-901 87 Umea, Sweden. 00282 * 00283 *> \par References: 00284 * ================ 00285 *> 00286 *> \verbatim 00287 *> 00288 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00289 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00290 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00291 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00292 *> 00293 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00294 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00295 *> Estimation: Theory, Algorithms and Software, Report 00296 *> UMINF - 94.04, Department of Computing Science, Umea University, 00297 *> S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. 00298 *> To appear in Numerical Algorithms, 1996. 00299 *> 00300 *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 00301 *> for Solving the Generalized Sylvester Equation and Estimating the 00302 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 00303 *> Department of Computing Science, Umea University, S-901 87 Umea, 00304 *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 00305 *> Note 75. 00306 *> To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996. 00307 *> \endverbatim 00308 *> 00309 * ===================================================================== 00310 SUBROUTINE CTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 00311 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 00312 $ IWORK, INFO ) 00313 * 00314 * -- LAPACK computational routine (version 3.4.0) -- 00315 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00316 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00317 * November 2011 00318 * 00319 * .. Scalar Arguments .. 00320 CHARACTER HOWMNY, JOB 00321 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 00322 * .. 00323 * .. Array Arguments .. 00324 LOGICAL SELECT( * ) 00325 INTEGER IWORK( * ) 00326 REAL DIF( * ), S( * ) 00327 COMPLEX A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 00328 $ VR( LDVR, * ), WORK( * ) 00329 * .. 00330 * 00331 * ===================================================================== 00332 * 00333 * .. Parameters .. 00334 REAL ZERO, ONE 00335 INTEGER IDIFJB 00336 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, IDIFJB = 3 ) 00337 * .. 00338 * .. Local Scalars .. 00339 LOGICAL LQUERY, SOMCON, WANTBH, WANTDF, WANTS 00340 INTEGER I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2 00341 REAL BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM 00342 COMPLEX YHAX, YHBX 00343 * .. 00344 * .. Local Arrays .. 00345 COMPLEX DUMMY( 1 ), DUMMY1( 1 ) 00346 * .. 00347 * .. External Functions .. 00348 LOGICAL LSAME 00349 REAL SCNRM2, SLAMCH, SLAPY2 00350 COMPLEX CDOTC 00351 EXTERNAL LSAME, SCNRM2, SLAMCH, SLAPY2, CDOTC 00352 * .. 00353 * .. External Subroutines .. 00354 EXTERNAL CGEMV, CLACPY, CTGEXC, CTGSYL, SLABAD, XERBLA 00355 * .. 00356 * .. Intrinsic Functions .. 00357 INTRINSIC ABS, CMPLX, MAX 00358 * .. 00359 * .. Executable Statements .. 00360 * 00361 * Decode and test the input parameters 00362 * 00363 WANTBH = LSAME( JOB, 'B' ) 00364 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00365 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH 00366 * 00367 SOMCON = LSAME( HOWMNY, 'S' ) 00368 * 00369 INFO = 0 00370 LQUERY = ( LWORK.EQ.-1 ) 00371 * 00372 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN 00373 INFO = -1 00374 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00375 INFO = -2 00376 ELSE IF( N.LT.0 ) THEN 00377 INFO = -4 00378 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00379 INFO = -6 00380 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00381 INFO = -8 00382 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN 00383 INFO = -10 00384 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN 00385 INFO = -12 00386 ELSE 00387 * 00388 * Set M to the number of eigenpairs for which condition numbers 00389 * are required, and test MM. 00390 * 00391 IF( SOMCON ) THEN 00392 M = 0 00393 DO 10 K = 1, N 00394 IF( SELECT( K ) ) 00395 $ M = M + 1 00396 10 CONTINUE 00397 ELSE 00398 M = N 00399 END IF 00400 * 00401 IF( N.EQ.0 ) THEN 00402 LWMIN = 1 00403 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN 00404 LWMIN = 2*N*N 00405 ELSE 00406 LWMIN = N 00407 END IF 00408 WORK( 1 ) = LWMIN 00409 * 00410 IF( MM.LT.M ) THEN 00411 INFO = -15 00412 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00413 INFO = -18 00414 END IF 00415 END IF 00416 * 00417 IF( INFO.NE.0 ) THEN 00418 CALL XERBLA( 'CTGSNA', -INFO ) 00419 RETURN 00420 ELSE IF( LQUERY ) THEN 00421 RETURN 00422 END IF 00423 * 00424 * Quick return if possible 00425 * 00426 IF( N.EQ.0 ) 00427 $ RETURN 00428 * 00429 * Get machine constants 00430 * 00431 EPS = SLAMCH( 'P' ) 00432 SMLNUM = SLAMCH( 'S' ) / EPS 00433 BIGNUM = ONE / SMLNUM 00434 CALL SLABAD( SMLNUM, BIGNUM ) 00435 KS = 0 00436 DO 20 K = 1, N 00437 * 00438 * Determine whether condition numbers are required for the k-th 00439 * eigenpair. 00440 * 00441 IF( SOMCON ) THEN 00442 IF( .NOT.SELECT( K ) ) 00443 $ GO TO 20 00444 END IF 00445 * 00446 KS = KS + 1 00447 * 00448 IF( WANTS ) THEN 00449 * 00450 * Compute the reciprocal condition number of the k-th 00451 * eigenvalue. 00452 * 00453 RNRM = SCNRM2( N, VR( 1, KS ), 1 ) 00454 LNRM = SCNRM2( N, VL( 1, KS ), 1 ) 00455 CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), A, LDA, 00456 $ VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 ) 00457 YHAX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 ) 00458 CALL CGEMV( 'N', N, N, CMPLX( ONE, ZERO ), B, LDB, 00459 $ VR( 1, KS ), 1, CMPLX( ZERO, ZERO ), WORK, 1 ) 00460 YHBX = CDOTC( N, WORK, 1, VL( 1, KS ), 1 ) 00461 COND = SLAPY2( ABS( YHAX ), ABS( YHBX ) ) 00462 IF( COND.EQ.ZERO ) THEN 00463 S( KS ) = -ONE 00464 ELSE 00465 S( KS ) = COND / ( RNRM*LNRM ) 00466 END IF 00467 END IF 00468 * 00469 IF( WANTDF ) THEN 00470 IF( N.EQ.1 ) THEN 00471 DIF( KS ) = SLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) ) 00472 ELSE 00473 * 00474 * Estimate the reciprocal condition number of the k-th 00475 * eigenvectors. 00476 * 00477 * Copy the matrix (A, B) to the array WORK and move the 00478 * (k,k)th pair to the (1,1) position. 00479 * 00480 CALL CLACPY( 'Full', N, N, A, LDA, WORK, N ) 00481 CALL CLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 00482 IFST = K 00483 ILST = 1 00484 * 00485 CALL CTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), 00486 $ N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR ) 00487 * 00488 IF( IERR.GT.0 ) THEN 00489 * 00490 * Ill-conditioned problem - swap rejected. 00491 * 00492 DIF( KS ) = ZERO 00493 ELSE 00494 * 00495 * Reordering successful, solve generalized Sylvester 00496 * equation for R and L, 00497 * A22 * R - L * A11 = A12 00498 * B22 * R - L * B11 = B12, 00499 * and compute estimate of Difl[(A11,B11), (A22, B22)]. 00500 * 00501 N1 = 1 00502 N2 = N - N1 00503 I = N*N + 1 00504 CALL CTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ), 00505 $ N, WORK, N, WORK( N1+1 ), N, 00506 $ WORK( N*N1+N1+I ), N, WORK( I ), N, 00507 $ WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY, 00508 $ 1, IWORK, IERR ) 00509 END IF 00510 END IF 00511 END IF 00512 * 00513 20 CONTINUE 00514 WORK( 1 ) = LWMIN 00515 RETURN 00516 * 00517 * End of CTGSNA 00518 * 00519 END