![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTGSNA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTGSNA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsna.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsna.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsna.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTGSNA( 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 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( * ), 00033 * $ VL( LDVL, * ), VR( LDVR, * ), WORK( * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> DTGSNA estimates reciprocal condition numbers for specified 00043 *> eigenvalues and/or eigenvectors of a matrix pair (A, B) in 00044 *> generalized real Schur canonical form (or of any matrix pair 00045 *> (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where 00046 *> Z**T denotes the transpose of Z. 00047 *> 00048 *> (A, B) must be in generalized real Schur form (as returned by DGGES), 00049 *> i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal 00050 *> blocks. B is upper triangular. 00051 *> 00052 *> \endverbatim 00053 * 00054 * Arguments: 00055 * ========== 00056 * 00057 *> \param[in] JOB 00058 *> \verbatim 00059 *> JOB is CHARACTER*1 00060 *> Specifies whether condition numbers are required for 00061 *> eigenvalues (S) or eigenvectors (DIF): 00062 *> = 'E': for eigenvalues only (S); 00063 *> = 'V': for eigenvectors only (DIF); 00064 *> = 'B': for both eigenvalues and eigenvectors (S and DIF). 00065 *> \endverbatim 00066 *> 00067 *> \param[in] HOWMNY 00068 *> \verbatim 00069 *> HOWMNY is CHARACTER*1 00070 *> = 'A': compute condition numbers for all eigenpairs; 00071 *> = 'S': compute condition numbers for selected eigenpairs 00072 *> specified by the array SELECT. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] SELECT 00076 *> \verbatim 00077 *> SELECT is LOGICAL array, dimension (N) 00078 *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which 00079 *> condition numbers are required. To select condition numbers 00080 *> for the eigenpair corresponding to a real eigenvalue w(j), 00081 *> SELECT(j) must be set to .TRUE.. To select condition numbers 00082 *> corresponding to a complex conjugate pair of eigenvalues w(j) 00083 *> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be 00084 *> set to .TRUE.. 00085 *> If HOWMNY = 'A', SELECT is not referenced. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] N 00089 *> \verbatim 00090 *> N is INTEGER 00091 *> The order of the square matrix pair (A, B). N >= 0. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] A 00095 *> \verbatim 00096 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00097 *> The upper quasi-triangular matrix A in the pair (A,B). 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDA 00101 *> \verbatim 00102 *> LDA is INTEGER 00103 *> The leading dimension of the array A. LDA >= max(1,N). 00104 *> \endverbatim 00105 *> 00106 *> \param[in] B 00107 *> \verbatim 00108 *> B is DOUBLE PRECISION array, dimension (LDB,N) 00109 *> The upper triangular matrix B in the pair (A,B). 00110 *> \endverbatim 00111 *> 00112 *> \param[in] LDB 00113 *> \verbatim 00114 *> LDB is INTEGER 00115 *> The leading dimension of the array B. LDB >= max(1,N). 00116 *> \endverbatim 00117 *> 00118 *> \param[in] VL 00119 *> \verbatim 00120 *> VL is DOUBLE PRECISION array, dimension (LDVL,M) 00121 *> If JOB = 'E' or 'B', VL must contain left eigenvectors of 00122 *> (A, B), corresponding to the eigenpairs specified by HOWMNY 00123 *> and SELECT. The eigenvectors must be stored in consecutive 00124 *> columns of VL, as returned by DTGEVC. 00125 *> If JOB = 'V', VL is not referenced. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] LDVL 00129 *> \verbatim 00130 *> LDVL is INTEGER 00131 *> The leading dimension of the array VL. LDVL >= 1. 00132 *> If JOB = 'E' or 'B', LDVL >= N. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] VR 00136 *> \verbatim 00137 *> VR is DOUBLE PRECISION array, dimension (LDVR,M) 00138 *> If JOB = 'E' or 'B', VR must contain right eigenvectors of 00139 *> (A, B), corresponding to the eigenpairs specified by HOWMNY 00140 *> and SELECT. The eigenvectors must be stored in consecutive 00141 *> columns ov VR, as returned by DTGEVC. 00142 *> If JOB = 'V', VR is not referenced. 00143 *> \endverbatim 00144 *> 00145 *> \param[in] LDVR 00146 *> \verbatim 00147 *> LDVR is INTEGER 00148 *> The leading dimension of the array VR. LDVR >= 1. 00149 *> If JOB = 'E' or 'B', LDVR >= N. 00150 *> \endverbatim 00151 *> 00152 *> \param[out] S 00153 *> \verbatim 00154 *> S is DOUBLE PRECISION array, dimension (MM) 00155 *> If JOB = 'E' or 'B', the reciprocal condition numbers of the 00156 *> selected eigenvalues, stored in consecutive elements of the 00157 *> array. For a complex conjugate pair of eigenvalues two 00158 *> consecutive elements of S are set to the same value. Thus 00159 *> S(j), DIF(j), and the j-th columns of VL and VR all 00160 *> correspond to the same eigenpair (but not in general the 00161 *> j-th eigenpair, unless all eigenpairs are selected). 00162 *> If JOB = 'V', S is not referenced. 00163 *> \endverbatim 00164 *> 00165 *> \param[out] DIF 00166 *> \verbatim 00167 *> DIF is DOUBLE PRECISION array, dimension (MM) 00168 *> If JOB = 'V' or 'B', the estimated reciprocal condition 00169 *> numbers of the selected eigenvectors, stored in consecutive 00170 *> elements of the array. For a complex eigenvector two 00171 *> consecutive elements of DIF are set to the same value. If 00172 *> the eigenvalues cannot be reordered to compute DIF(j), DIF(j) 00173 *> is set to 0; this can only occur when the true value would be 00174 *> very small anyway. 00175 *> If JOB = 'E', DIF is not referenced. 00176 *> \endverbatim 00177 *> 00178 *> \param[in] MM 00179 *> \verbatim 00180 *> MM is INTEGER 00181 *> The number of elements in the arrays S and DIF. MM >= M. 00182 *> \endverbatim 00183 *> 00184 *> \param[out] M 00185 *> \verbatim 00186 *> M is INTEGER 00187 *> The number of elements of the arrays S and DIF used to store 00188 *> the specified condition numbers; for each selected real 00189 *> eigenvalue one element is used, and for each selected complex 00190 *> conjugate pair of eigenvalues, two elements are used. 00191 *> If HOWMNY = 'A', M is set to N. 00192 *> \endverbatim 00193 *> 00194 *> \param[out] WORK 00195 *> \verbatim 00196 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00197 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00198 *> \endverbatim 00199 *> 00200 *> \param[in] LWORK 00201 *> \verbatim 00202 *> LWORK is INTEGER 00203 *> The dimension of the array WORK. LWORK >= max(1,N). 00204 *> If JOB = 'V' or 'B' LWORK >= 2*N*(N+2)+16. 00205 *> 00206 *> If LWORK = -1, then a workspace query is assumed; the routine 00207 *> only calculates the optimal size of the WORK array, returns 00208 *> this value as the first entry of the WORK array, and no error 00209 *> message related to LWORK is issued by XERBLA. 00210 *> \endverbatim 00211 *> 00212 *> \param[out] IWORK 00213 *> \verbatim 00214 *> IWORK is INTEGER array, dimension (N + 6) 00215 *> If JOB = 'E', IWORK is not referenced. 00216 *> \endverbatim 00217 *> 00218 *> \param[out] INFO 00219 *> \verbatim 00220 *> INFO is INTEGER 00221 *> =0: Successful exit 00222 *> <0: If INFO = -i, the i-th argument had an illegal value 00223 *> \endverbatim 00224 * 00225 * Authors: 00226 * ======== 00227 * 00228 *> \author Univ. of Tennessee 00229 *> \author Univ. of California Berkeley 00230 *> \author Univ. of Colorado Denver 00231 *> \author NAG Ltd. 00232 * 00233 *> \date November 2011 00234 * 00235 *> \ingroup doubleOTHERcomputational 00236 * 00237 *> \par Further Details: 00238 * ===================== 00239 *> 00240 *> \verbatim 00241 *> 00242 *> The reciprocal of the condition number of a generalized eigenvalue 00243 *> w = (a, b) is defined as 00244 *> 00245 *> S(w) = (|u**TAv|**2 + |u**TBv|**2)**(1/2) / (norm(u)*norm(v)) 00246 *> 00247 *> where u and v are the left and right eigenvectors of (A, B) 00248 *> corresponding to w; |z| denotes the absolute value of the complex 00249 *> number, and norm(u) denotes the 2-norm of the vector u. 00250 *> The pair (a, b) corresponds to an eigenvalue w = a/b (= u**TAv/u**TBv) 00251 *> of the matrix pair (A, B). If both a and b equal zero, then (A B) is 00252 *> singular and S(I) = -1 is returned. 00253 *> 00254 *> An approximate error bound on the chordal distance between the i-th 00255 *> computed generalized eigenvalue w and the corresponding exact 00256 *> eigenvalue lambda is 00257 *> 00258 *> chord(w, lambda) <= EPS * norm(A, B) / S(I) 00259 *> 00260 *> where EPS is the machine precision. 00261 *> 00262 *> The reciprocal of the condition number DIF(i) of right eigenvector u 00263 *> and left eigenvector v corresponding to the generalized eigenvalue w 00264 *> is defined as follows: 00265 *> 00266 *> a) If the i-th eigenvalue w = (a,b) is real 00267 *> 00268 *> Suppose U and V are orthogonal transformations such that 00269 *> 00270 *> U**T*(A, B)*V = (S, T) = ( a * ) ( b * ) 1 00271 *> ( 0 S22 ),( 0 T22 ) n-1 00272 *> 1 n-1 1 n-1 00273 *> 00274 *> Then the reciprocal condition number DIF(i) is 00275 *> 00276 *> Difl((a, b), (S22, T22)) = sigma-min( Zl ), 00277 *> 00278 *> where sigma-min(Zl) denotes the smallest singular value of the 00279 *> 2(n-1)-by-2(n-1) matrix 00280 *> 00281 *> Zl = [ kron(a, In-1) -kron(1, S22) ] 00282 *> [ kron(b, In-1) -kron(1, T22) ] . 00283 *> 00284 *> Here In-1 is the identity matrix of size n-1. kron(X, Y) is the 00285 *> Kronecker product between the matrices X and Y. 00286 *> 00287 *> Note that if the default method for computing DIF(i) is wanted 00288 *> (see DLATDF), then the parameter DIFDRI (see below) should be 00289 *> changed from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). 00290 *> See DTGSYL for more details. 00291 *> 00292 *> b) If the i-th and (i+1)-th eigenvalues are complex conjugate pair, 00293 *> 00294 *> Suppose U and V are orthogonal transformations such that 00295 *> 00296 *> U**T*(A, B)*V = (S, T) = ( S11 * ) ( T11 * ) 2 00297 *> ( 0 S22 ),( 0 T22) n-2 00298 *> 2 n-2 2 n-2 00299 *> 00300 *> and (S11, T11) corresponds to the complex conjugate eigenvalue 00301 *> pair (w, conjg(w)). There exist unitary matrices U1 and V1 such 00302 *> that 00303 *> 00304 *> U1**T*S11*V1 = ( s11 s12 ) and U1**T*T11*V1 = ( t11 t12 ) 00305 *> ( 0 s22 ) ( 0 t22 ) 00306 *> 00307 *> where the generalized eigenvalues w = s11/t11 and 00308 *> conjg(w) = s22/t22. 00309 *> 00310 *> Then the reciprocal condition number DIF(i) is bounded by 00311 *> 00312 *> min( d1, max( 1, |real(s11)/real(s22)| )*d2 ) 00313 *> 00314 *> where, d1 = Difl((s11, t11), (s22, t22)) = sigma-min(Z1), where 00315 *> Z1 is the complex 2-by-2 matrix 00316 *> 00317 *> Z1 = [ s11 -s22 ] 00318 *> [ t11 -t22 ], 00319 *> 00320 *> This is done by computing (using real arithmetic) the 00321 *> roots of the characteristical polynomial det(Z1**T * Z1 - lambda I), 00322 *> where Z1**T denotes the transpose of Z1 and det(X) denotes 00323 *> the determinant of X. 00324 *> 00325 *> and d2 is an upper bound on Difl((S11, T11), (S22, T22)), i.e. an 00326 *> upper bound on sigma-min(Z2), where Z2 is (2n-2)-by-(2n-2) 00327 *> 00328 *> Z2 = [ kron(S11**T, In-2) -kron(I2, S22) ] 00329 *> [ kron(T11**T, In-2) -kron(I2, T22) ] 00330 *> 00331 *> Note that if the default method for computing DIF is wanted (see 00332 *> DLATDF), then the parameter DIFDRI (see below) should be changed 00333 *> from 3 to 4 (routine DLATDF(IJOB = 2 will be used)). See DTGSYL 00334 *> for more details. 00335 *> 00336 *> For each eigenvalue/vector specified by SELECT, DIF stores a 00337 *> Frobenius norm-based estimate of Difl. 00338 *> 00339 *> An approximate error bound for the i-th computed eigenvector VL(i) or 00340 *> VR(i) is given by 00341 *> 00342 *> EPS * norm(A, B) / DIF(i). 00343 *> 00344 *> See ref. [2-3] for more details and further references. 00345 *> \endverbatim 00346 * 00347 *> \par Contributors: 00348 * ================== 00349 *> 00350 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00351 *> Umea University, S-901 87 Umea, Sweden. 00352 * 00353 *> \par References: 00354 * ================ 00355 *> 00356 *> \verbatim 00357 *> 00358 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00359 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00360 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00361 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00362 *> 00363 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00364 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00365 *> Estimation: Theory, Algorithms and Software, 00366 *> Report UMINF - 94.04, Department of Computing Science, Umea 00367 *> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working 00368 *> Note 87. To appear in Numerical Algorithms, 1996. 00369 *> 00370 *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 00371 *> for Solving the Generalized Sylvester Equation and Estimating the 00372 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 00373 *> Department of Computing Science, Umea University, S-901 87 Umea, 00374 *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 00375 *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, 00376 *> No 1, 1996. 00377 *> \endverbatim 00378 *> 00379 * ===================================================================== 00380 SUBROUTINE DTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 00381 $ LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, 00382 $ IWORK, INFO ) 00383 * 00384 * -- LAPACK computational routine (version 3.4.0) -- 00385 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00386 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00387 * November 2011 00388 * 00389 * .. Scalar Arguments .. 00390 CHARACTER HOWMNY, JOB 00391 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N 00392 * .. 00393 * .. Array Arguments .. 00394 LOGICAL SELECT( * ) 00395 INTEGER IWORK( * ) 00396 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DIF( * ), S( * ), 00397 $ VL( LDVL, * ), VR( LDVR, * ), WORK( * ) 00398 * .. 00399 * 00400 * ===================================================================== 00401 * 00402 * .. Parameters .. 00403 INTEGER DIFDRI 00404 PARAMETER ( DIFDRI = 3 ) 00405 DOUBLE PRECISION ZERO, ONE, TWO, FOUR 00406 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, 00407 $ FOUR = 4.0D+0 ) 00408 * .. 00409 * .. Local Scalars .. 00410 LOGICAL LQUERY, PAIR, SOMCON, WANTBH, WANTDF, WANTS 00411 INTEGER I, IERR, IFST, ILST, IZ, K, KS, LWMIN, N1, N2 00412 DOUBLE PRECISION ALPHAI, ALPHAR, ALPRQT, BETA, C1, C2, COND, 00413 $ EPS, LNRM, RNRM, ROOT1, ROOT2, SCALE, SMLNUM, 00414 $ TMPII, TMPIR, TMPRI, TMPRR, UHAV, UHAVI, UHBV, 00415 $ UHBVI 00416 * .. 00417 * .. Local Arrays .. 00418 DOUBLE PRECISION DUMMY( 1 ), DUMMY1( 1 ) 00419 * .. 00420 * .. External Functions .. 00421 LOGICAL LSAME 00422 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2 00423 EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2 00424 * .. 00425 * .. External Subroutines .. 00426 EXTERNAL DGEMV, DLACPY, DLAG2, DTGEXC, DTGSYL, XERBLA 00427 * .. 00428 * .. Intrinsic Functions .. 00429 INTRINSIC MAX, MIN, SQRT 00430 * .. 00431 * .. Executable Statements .. 00432 * 00433 * Decode and test the input parameters 00434 * 00435 WANTBH = LSAME( JOB, 'B' ) 00436 WANTS = LSAME( JOB, 'E' ) .OR. WANTBH 00437 WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH 00438 * 00439 SOMCON = LSAME( HOWMNY, 'S' ) 00440 * 00441 INFO = 0 00442 LQUERY = ( LWORK.EQ.-1 ) 00443 * 00444 IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN 00445 INFO = -1 00446 ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN 00447 INFO = -2 00448 ELSE IF( N.LT.0 ) THEN 00449 INFO = -4 00450 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00451 INFO = -6 00452 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00453 INFO = -8 00454 ELSE IF( WANTS .AND. LDVL.LT.N ) THEN 00455 INFO = -10 00456 ELSE IF( WANTS .AND. LDVR.LT.N ) THEN 00457 INFO = -12 00458 ELSE 00459 * 00460 * Set M to the number of eigenpairs for which condition numbers 00461 * are required, and test MM. 00462 * 00463 IF( SOMCON ) THEN 00464 M = 0 00465 PAIR = .FALSE. 00466 DO 10 K = 1, N 00467 IF( PAIR ) THEN 00468 PAIR = .FALSE. 00469 ELSE 00470 IF( K.LT.N ) THEN 00471 IF( A( K+1, K ).EQ.ZERO ) THEN 00472 IF( SELECT( K ) ) 00473 $ M = M + 1 00474 ELSE 00475 PAIR = .TRUE. 00476 IF( SELECT( K ) .OR. SELECT( K+1 ) ) 00477 $ M = M + 2 00478 END IF 00479 ELSE 00480 IF( SELECT( N ) ) 00481 $ M = M + 1 00482 END IF 00483 END IF 00484 10 CONTINUE 00485 ELSE 00486 M = N 00487 END IF 00488 * 00489 IF( N.EQ.0 ) THEN 00490 LWMIN = 1 00491 ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN 00492 LWMIN = 2*N*( N + 2 ) + 16 00493 ELSE 00494 LWMIN = N 00495 END IF 00496 WORK( 1 ) = LWMIN 00497 * 00498 IF( MM.LT.M ) THEN 00499 INFO = -15 00500 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00501 INFO = -18 00502 END IF 00503 END IF 00504 * 00505 IF( INFO.NE.0 ) THEN 00506 CALL XERBLA( 'DTGSNA', -INFO ) 00507 RETURN 00508 ELSE IF( LQUERY ) THEN 00509 RETURN 00510 END IF 00511 * 00512 * Quick return if possible 00513 * 00514 IF( N.EQ.0 ) 00515 $ RETURN 00516 * 00517 * Get machine constants 00518 * 00519 EPS = DLAMCH( 'P' ) 00520 SMLNUM = DLAMCH( 'S' ) / EPS 00521 KS = 0 00522 PAIR = .FALSE. 00523 * 00524 DO 20 K = 1, N 00525 * 00526 * Determine whether A(k,k) begins a 1-by-1 or 2-by-2 block. 00527 * 00528 IF( PAIR ) THEN 00529 PAIR = .FALSE. 00530 GO TO 20 00531 ELSE 00532 IF( K.LT.N ) 00533 $ PAIR = A( K+1, K ).NE.ZERO 00534 END IF 00535 * 00536 * Determine whether condition numbers are required for the k-th 00537 * eigenpair. 00538 * 00539 IF( SOMCON ) THEN 00540 IF( PAIR ) THEN 00541 IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) ) 00542 $ GO TO 20 00543 ELSE 00544 IF( .NOT.SELECT( K ) ) 00545 $ GO TO 20 00546 END IF 00547 END IF 00548 * 00549 KS = KS + 1 00550 * 00551 IF( WANTS ) THEN 00552 * 00553 * Compute the reciprocal condition number of the k-th 00554 * eigenvalue. 00555 * 00556 IF( PAIR ) THEN 00557 * 00558 * Complex eigenvalue pair. 00559 * 00560 RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ), 00561 $ DNRM2( N, VR( 1, KS+1 ), 1 ) ) 00562 LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ), 00563 $ DNRM2( N, VL( 1, KS+1 ), 1 ) ) 00564 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO, 00565 $ WORK, 1 ) 00566 TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 ) 00567 TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 00568 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS+1 ), 1, 00569 $ ZERO, WORK, 1 ) 00570 TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 00571 TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 ) 00572 UHAV = TMPRR + TMPII 00573 UHAVI = TMPIR - TMPRI 00574 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO, 00575 $ WORK, 1 ) 00576 TMPRR = DDOT( N, WORK, 1, VL( 1, KS ), 1 ) 00577 TMPRI = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 00578 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS+1 ), 1, 00579 $ ZERO, WORK, 1 ) 00580 TMPII = DDOT( N, WORK, 1, VL( 1, KS+1 ), 1 ) 00581 TMPIR = DDOT( N, WORK, 1, VL( 1, KS ), 1 ) 00582 UHBV = TMPRR + TMPII 00583 UHBVI = TMPIR - TMPRI 00584 UHAV = DLAPY2( UHAV, UHAVI ) 00585 UHBV = DLAPY2( UHBV, UHBVI ) 00586 COND = DLAPY2( UHAV, UHBV ) 00587 S( KS ) = COND / ( RNRM*LNRM ) 00588 S( KS+1 ) = S( KS ) 00589 * 00590 ELSE 00591 * 00592 * Real eigenvalue. 00593 * 00594 RNRM = DNRM2( N, VR( 1, KS ), 1 ) 00595 LNRM = DNRM2( N, VL( 1, KS ), 1 ) 00596 CALL DGEMV( 'N', N, N, ONE, A, LDA, VR( 1, KS ), 1, ZERO, 00597 $ WORK, 1 ) 00598 UHAV = DDOT( N, WORK, 1, VL( 1, KS ), 1 ) 00599 CALL DGEMV( 'N', N, N, ONE, B, LDB, VR( 1, KS ), 1, ZERO, 00600 $ WORK, 1 ) 00601 UHBV = DDOT( N, WORK, 1, VL( 1, KS ), 1 ) 00602 COND = DLAPY2( UHAV, UHBV ) 00603 IF( COND.EQ.ZERO ) THEN 00604 S( KS ) = -ONE 00605 ELSE 00606 S( KS ) = COND / ( RNRM*LNRM ) 00607 END IF 00608 END IF 00609 END IF 00610 * 00611 IF( WANTDF ) THEN 00612 IF( N.EQ.1 ) THEN 00613 DIF( KS ) = DLAPY2( A( 1, 1 ), B( 1, 1 ) ) 00614 GO TO 20 00615 END IF 00616 * 00617 * Estimate the reciprocal condition number of the k-th 00618 * eigenvectors. 00619 IF( PAIR ) THEN 00620 * 00621 * Copy the 2-by 2 pencil beginning at (A(k,k), B(k, k)). 00622 * Compute the eigenvalue(s) at position K. 00623 * 00624 WORK( 1 ) = A( K, K ) 00625 WORK( 2 ) = A( K+1, K ) 00626 WORK( 3 ) = A( K, K+1 ) 00627 WORK( 4 ) = A( K+1, K+1 ) 00628 WORK( 5 ) = B( K, K ) 00629 WORK( 6 ) = B( K+1, K ) 00630 WORK( 7 ) = B( K, K+1 ) 00631 WORK( 8 ) = B( K+1, K+1 ) 00632 CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA, 00633 $ DUMMY1( 1 ), ALPHAR, DUMMY( 1 ), ALPHAI ) 00634 ALPRQT = ONE 00635 C1 = TWO*( ALPHAR*ALPHAR+ALPHAI*ALPHAI+BETA*BETA ) 00636 C2 = FOUR*BETA*BETA*ALPHAI*ALPHAI 00637 ROOT1 = C1 + SQRT( C1*C1-4.0D0*C2 ) 00638 ROOT2 = C2 / ROOT1 00639 ROOT1 = ROOT1 / TWO 00640 COND = MIN( SQRT( ROOT1 ), SQRT( ROOT2 ) ) 00641 END IF 00642 * 00643 * Copy the matrix (A, B) to the array WORK and swap the 00644 * diagonal block beginning at A(k,k) to the (1,1) position. 00645 * 00646 CALL DLACPY( 'Full', N, N, A, LDA, WORK, N ) 00647 CALL DLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N ) 00648 IFST = K 00649 ILST = 1 00650 * 00651 CALL DTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ), N, 00652 $ DUMMY, 1, DUMMY1, 1, IFST, ILST, 00653 $ WORK( N*N*2+1 ), LWORK-2*N*N, IERR ) 00654 * 00655 IF( IERR.GT.0 ) THEN 00656 * 00657 * Ill-conditioned problem - swap rejected. 00658 * 00659 DIF( KS ) = ZERO 00660 ELSE 00661 * 00662 * Reordering successful, solve generalized Sylvester 00663 * equation for R and L, 00664 * A22 * R - L * A11 = A12 00665 * B22 * R - L * B11 = B12, 00666 * and compute estimate of Difl((A11,B11), (A22, B22)). 00667 * 00668 N1 = 1 00669 IF( WORK( 2 ).NE.ZERO ) 00670 $ N1 = 2 00671 N2 = N - N1 00672 IF( N2.EQ.0 ) THEN 00673 DIF( KS ) = COND 00674 ELSE 00675 I = N*N + 1 00676 IZ = 2*N*N + 1 00677 CALL DTGSYL( 'N', DIFDRI, N2, N1, WORK( N*N1+N1+1 ), 00678 $ N, WORK, N, WORK( N1+1 ), N, 00679 $ WORK( N*N1+N1+I ), N, WORK( I ), N, 00680 $ WORK( N1+I ), N, SCALE, DIF( KS ), 00681 $ WORK( IZ+1 ), LWORK-2*N*N, IWORK, IERR ) 00682 * 00683 IF( PAIR ) 00684 $ DIF( KS ) = MIN( MAX( ONE, ALPRQT )*DIF( KS ), 00685 $ COND ) 00686 END IF 00687 END IF 00688 IF( PAIR ) 00689 $ DIF( KS+1 ) = DIF( KS ) 00690 END IF 00691 IF( PAIR ) 00692 $ KS = KS + 1 00693 * 00694 20 CONTINUE 00695 WORK( 1 ) = LWMIN 00696 RETURN 00697 * 00698 * End of DTGSNA 00699 * 00700 END