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