00001 SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, 00002 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, 00003 $ PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * January 2007 00009 * 00010 * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. 00011 * 00012 * .. Scalar Arguments .. 00013 LOGICAL WANTQ, WANTZ 00014 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, 00015 $ M, N 00016 DOUBLE PRECISION PL, PR 00017 * .. 00018 * .. Array Arguments .. 00019 LOGICAL SELECT( * ) 00020 INTEGER IWORK( * ) 00021 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 00022 $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), 00023 $ WORK( * ), Z( LDZ, * ) 00024 * .. 00025 * 00026 * Purpose 00027 * ======= 00028 * 00029 * DTGSEN reorders the generalized real Schur decomposition of a real 00030 * matrix pair (A, B) (in terms of an orthonormal equivalence trans- 00031 * formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues 00032 * appears in the leading diagonal blocks of the upper quasi-triangular 00033 * matrix A and the upper triangular B. The leading columns of Q and 00034 * Z form orthonormal bases of the corresponding left and right eigen- 00035 * spaces (deflating subspaces). (A, B) must be in generalized real 00036 * Schur canonical form (as returned by DGGES), i.e. A is block upper 00037 * triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper 00038 * triangular. 00039 * 00040 * DTGSEN also computes the generalized eigenvalues 00041 * 00042 * w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) 00043 * 00044 * of the reordered matrix pair (A, B). 00045 * 00046 * Optionally, DTGSEN computes the estimates of reciprocal condition 00047 * numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), 00048 * (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) 00049 * between the matrix pairs (A11, B11) and (A22,B22) that correspond to 00050 * the selected cluster and the eigenvalues outside the cluster, resp., 00051 * and norms of "projections" onto left and right eigenspaces w.r.t. 00052 * the selected cluster in the (1,1)-block. 00053 * 00054 * Arguments 00055 * ========= 00056 * 00057 * IJOB (input) INTEGER 00058 * Specifies whether condition numbers are required for the 00059 * cluster of eigenvalues (PL and PR) or the deflating subspaces 00060 * (Difu and Difl): 00061 * =0: Only reorder w.r.t. SELECT. No extras. 00062 * =1: Reciprocal of norms of "projections" onto left and right 00063 * eigenspaces w.r.t. the selected cluster (PL and PR). 00064 * =2: Upper bounds on Difu and Difl. F-norm-based estimate 00065 * (DIF(1:2)). 00066 * =3: Estimate of Difu and Difl. 1-norm-based estimate 00067 * (DIF(1:2)). 00068 * About 5 times as expensive as IJOB = 2. 00069 * =4: Compute PL, PR and DIF (i.e. 0, 1 and 2 above): Economic 00070 * version to get it all. 00071 * =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above) 00072 * 00073 * WANTQ (input) LOGICAL 00074 * .TRUE. : update the left transformation matrix Q; 00075 * .FALSE.: do not update Q. 00076 * 00077 * WANTZ (input) LOGICAL 00078 * .TRUE. : update the right transformation matrix Z; 00079 * .FALSE.: do not update Z. 00080 * 00081 * SELECT (input) LOGICAL array, dimension (N) 00082 * SELECT specifies the eigenvalues in the selected cluster. 00083 * To select a real eigenvalue w(j), SELECT(j) must be set to 00084 * .TRUE.. To select a complex conjugate pair of eigenvalues 00085 * w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, 00086 * either SELECT(j) or SELECT(j+1) or both must be set to 00087 * .TRUE.; a complex conjugate pair of eigenvalues must be 00088 * either both included in the cluster or both excluded. 00089 * 00090 * N (input) INTEGER 00091 * The order of the matrices A and B. N >= 0. 00092 * 00093 * A (input/output) DOUBLE PRECISION array, dimension(LDA,N) 00094 * On entry, the upper quasi-triangular matrix A, with (A, B) in 00095 * generalized real Schur canonical form. 00096 * On exit, A is overwritten by the reordered matrix A. 00097 * 00098 * LDA (input) INTEGER 00099 * The leading dimension of the array A. LDA >= max(1,N). 00100 * 00101 * B (input/output) DOUBLE PRECISION array, dimension(LDB,N) 00102 * On entry, the upper triangular matrix B, with (A, B) in 00103 * generalized real Schur canonical form. 00104 * On exit, B is overwritten by the reordered matrix B. 00105 * 00106 * LDB (input) INTEGER 00107 * The leading dimension of the array B. LDB >= max(1,N). 00108 * 00109 * ALPHAR (output) DOUBLE PRECISION array, dimension (N) 00110 * ALPHAI (output) DOUBLE PRECISION array, dimension (N) 00111 * BETA (output) DOUBLE PRECISION array, dimension (N) 00112 * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will 00113 * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i 00114 * and BETA(j),j=1,...,N are the diagonals of the complex Schur 00115 * form (S,T) that would result if the 2-by-2 diagonal blocks of 00116 * the real generalized Schur form of (A,B) were further reduced 00117 * to triangular form using complex unitary transformations. 00118 * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if 00119 * positive, then the j-th and (j+1)-st eigenvalues are a 00120 * complex conjugate pair, with ALPHAI(j+1) negative. 00121 * 00122 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) 00123 * On entry, if WANTQ = .TRUE., Q is an N-by-N matrix. 00124 * On exit, Q has been postmultiplied by the left orthogonal 00125 * transformation matrix which reorder (A, B); The leading M 00126 * columns of Q form orthonormal bases for the specified pair of 00127 * left eigenspaces (deflating subspaces). 00128 * If WANTQ = .FALSE., Q is not referenced. 00129 * 00130 * LDQ (input) INTEGER 00131 * The leading dimension of the array Q. LDQ >= 1; 00132 * and if WANTQ = .TRUE., LDQ >= N. 00133 * 00134 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) 00135 * On entry, if WANTZ = .TRUE., Z is an N-by-N matrix. 00136 * On exit, Z has been postmultiplied by the left orthogonal 00137 * transformation matrix which reorder (A, B); The leading M 00138 * columns of Z form orthonormal bases for the specified pair of 00139 * left eigenspaces (deflating subspaces). 00140 * If WANTZ = .FALSE., Z is not referenced. 00141 * 00142 * LDZ (input) INTEGER 00143 * The leading dimension of the array Z. LDZ >= 1; 00144 * If WANTZ = .TRUE., LDZ >= N. 00145 * 00146 * M (output) INTEGER 00147 * The dimension of the specified pair of left and right eigen- 00148 * spaces (deflating subspaces). 0 <= M <= N. 00149 * 00150 * PL (output) DOUBLE PRECISION 00151 * PR (output) DOUBLE PRECISION 00152 * If IJOB = 1, 4 or 5, PL, PR are lower bounds on the 00153 * reciprocal of the norm of "projections" onto left and right 00154 * eigenspaces with respect to the selected cluster. 00155 * 0 < PL, PR <= 1. 00156 * If M = 0 or M = N, PL = PR = 1. 00157 * If IJOB = 0, 2 or 3, PL and PR are not referenced. 00158 * 00159 * DIF (output) DOUBLE PRECISION array, dimension (2). 00160 * If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl. 00161 * If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on 00162 * Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based 00163 * estimates of Difu and Difl. 00164 * If M = 0 or N, DIF(1:2) = F-norm([A, B]). 00165 * If IJOB = 0 or 1, DIF is not referenced. 00166 * 00167 * WORK (workspace/output) DOUBLE PRECISION array, 00168 * dimension (MAX(1,LWORK)) 00169 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00170 * 00171 * LWORK (input) INTEGER 00172 * The dimension of the array WORK. LWORK >= 4*N+16. 00173 * If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). 00174 * If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)). 00175 * 00176 * If LWORK = -1, then a workspace query is assumed; the routine 00177 * only calculates the optimal size of the WORK array, returns 00178 * this value as the first entry of the WORK array, and no error 00179 * message related to LWORK is issued by XERBLA. 00180 * 00181 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00182 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00183 * 00184 * LIWORK (input) INTEGER 00185 * The dimension of the array IWORK. LIWORK >= 1. 00186 * If IJOB = 1, 2 or 4, LIWORK >= N+6. 00187 * If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6). 00188 * 00189 * If LIWORK = -1, then a workspace query is assumed; the 00190 * routine only calculates the optimal size of the IWORK array, 00191 * returns this value as the first entry of the IWORK array, and 00192 * no error message related to LIWORK is issued by XERBLA. 00193 * 00194 * INFO (output) INTEGER 00195 * =0: Successful exit. 00196 * <0: If INFO = -i, the i-th argument had an illegal value. 00197 * =1: Reordering of (A, B) failed because the transformed 00198 * matrix pair (A, B) would be too far from generalized 00199 * Schur form; the problem is very ill-conditioned. 00200 * (A, B) may have been partially reordered. 00201 * If requested, 0 is returned in DIF(*), PL and PR. 00202 * 00203 * Further Details 00204 * =============== 00205 * 00206 * DTGSEN first collects the selected eigenvalues by computing 00207 * orthogonal U and W that move them to the top left corner of (A, B). 00208 * In other words, the selected eigenvalues are the eigenvalues of 00209 * (A11, B11) in: 00210 * 00211 * U'*(A, B)*W = (A11 A12) (B11 B12) n1 00212 * ( 0 A22),( 0 B22) n2 00213 * n1 n2 n1 n2 00214 * 00215 * where N = n1+n2 and U' means the transpose of U. The first n1 columns 00216 * of U and W span the specified pair of left and right eigenspaces 00217 * (deflating subspaces) of (A, B). 00218 * 00219 * If (A, B) has been obtained from the generalized real Schur 00220 * decomposition of a matrix pair (C, D) = Q*(A, B)*Z', then the 00221 * reordered generalized real Schur form of (C, D) is given by 00222 * 00223 * (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)', 00224 * 00225 * and the first n1 columns of Q*U and Z*W span the corresponding 00226 * deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.). 00227 * 00228 * Note that if the selected eigenvalue is sufficiently ill-conditioned, 00229 * then its value may differ significantly from its value before 00230 * reordering. 00231 * 00232 * The reciprocal condition numbers of the left and right eigenspaces 00233 * spanned by the first n1 columns of U and W (or Q*U and Z*W) may 00234 * be returned in DIF(1:2), corresponding to Difu and Difl, resp. 00235 * 00236 * The Difu and Difl are defined as: 00237 * 00238 * Difu[(A11, B11), (A22, B22)] = sigma-min( Zu ) 00239 * and 00240 * Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)], 00241 * 00242 * where sigma-min(Zu) is the smallest singular value of the 00243 * (2*n1*n2)-by-(2*n1*n2) matrix 00244 * 00245 * Zu = [ kron(In2, A11) -kron(A22', In1) ] 00246 * [ kron(In2, B11) -kron(B22', In1) ]. 00247 * 00248 * Here, Inx is the identity matrix of size nx and A22' is the 00249 * transpose of A22. kron(X, Y) is the Kronecker product between 00250 * the matrices X and Y. 00251 * 00252 * When DIF(2) is small, small changes in (A, B) can cause large changes 00253 * in the deflating subspace. An approximate (asymptotic) bound on the 00254 * maximum angular error in the computed deflating subspaces is 00255 * 00256 * EPS * norm((A, B)) / DIF(2), 00257 * 00258 * where EPS is the machine precision. 00259 * 00260 * The reciprocal norm of the projectors on the left and right 00261 * eigenspaces associated with (A11, B11) may be returned in PL and PR. 00262 * They are computed as follows. First we compute L and R so that 00263 * P*(A, B)*Q is block diagonal, where 00264 * 00265 * P = ( I -L ) n1 Q = ( I R ) n1 00266 * ( 0 I ) n2 and ( 0 I ) n2 00267 * n1 n2 n1 n2 00268 * 00269 * and (L, R) is the solution to the generalized Sylvester equation 00270 * 00271 * A11*R - L*A22 = -A12 00272 * B11*R - L*B22 = -B12 00273 * 00274 * Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2). 00275 * An approximate (asymptotic) bound on the average absolute error of 00276 * the selected eigenvalues is 00277 * 00278 * EPS * norm((A, B)) / PL. 00279 * 00280 * There are also global error bounds which valid for perturbations up 00281 * to a certain restriction: A lower bound (x) on the smallest 00282 * F-norm(E,F) for which an eigenvalue of (A11, B11) may move and 00283 * coalesce with an eigenvalue of (A22, B22) under perturbation (E,F), 00284 * (i.e. (A + E, B + F), is 00285 * 00286 * x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)). 00287 * 00288 * An approximate bound on x can be computed from DIF(1:2), PL and PR. 00289 * 00290 * If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed 00291 * (L', R') and unperturbed (L, R) left and right deflating subspaces 00292 * associated with the selected cluster in the (1,1)-blocks can be 00293 * bounded as 00294 * 00295 * max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2)) 00296 * max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2)) 00297 * 00298 * See LAPACK User's Guide section 4.11 or the following references 00299 * for more information. 00300 * 00301 * Note that if the default method for computing the Frobenius-norm- 00302 * based estimate DIF is not wanted (see DLATDF), then the parameter 00303 * IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF 00304 * (IJOB = 2 will be used)). See DTGSYL for more details. 00305 * 00306 * Based on contributions by 00307 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00308 * Umea University, S-901 87 Umea, Sweden. 00309 * 00310 * References 00311 * ========== 00312 * 00313 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00314 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00315 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00316 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00317 * 00318 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00319 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00320 * Estimation: Theory, Algorithms and Software, 00321 * Report UMINF - 94.04, Department of Computing Science, Umea 00322 * University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working 00323 * Note 87. To appear in Numerical Algorithms, 1996. 00324 * 00325 * [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 00326 * for Solving the Generalized Sylvester Equation and Estimating the 00327 * Separation between Regular Matrix Pairs, Report UMINF - 93.23, 00328 * Department of Computing Science, Umea University, S-901 87 Umea, 00329 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working 00330 * Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 00331 * 1996. 00332 * 00333 * ===================================================================== 00334 * 00335 * .. Parameters .. 00336 INTEGER IDIFJB 00337 PARAMETER ( IDIFJB = 3 ) 00338 DOUBLE PRECISION ZERO, ONE 00339 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00340 * .. 00341 * .. Local Scalars .. 00342 LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2, 00343 $ WANTP 00344 INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN, 00345 $ MN2, N1, N2 00346 DOUBLE PRECISION DSCALE, DSUM, EPS, RDSCAL, SMLNUM 00347 * .. 00348 * .. Local Arrays .. 00349 INTEGER ISAVE( 3 ) 00350 * .. 00351 * .. External Subroutines .. 00352 EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL, 00353 $ XERBLA 00354 * .. 00355 * .. External Functions .. 00356 DOUBLE PRECISION DLAMCH 00357 EXTERNAL DLAMCH 00358 * .. 00359 * .. Intrinsic Functions .. 00360 INTRINSIC MAX, SIGN, SQRT 00361 * .. 00362 * .. Executable Statements .. 00363 * 00364 * Decode and test the input parameters 00365 * 00366 INFO = 0 00367 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00368 * 00369 IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN 00370 INFO = -1 00371 ELSE IF( N.LT.0 ) THEN 00372 INFO = -5 00373 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00374 INFO = -7 00375 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00376 INFO = -9 00377 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00378 INFO = -14 00379 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00380 INFO = -16 00381 END IF 00382 * 00383 IF( INFO.NE.0 ) THEN 00384 CALL XERBLA( 'DTGSEN', -INFO ) 00385 RETURN 00386 END IF 00387 * 00388 * Get machine constants 00389 * 00390 EPS = DLAMCH( 'P' ) 00391 SMLNUM = DLAMCH( 'S' ) / EPS 00392 IERR = 0 00393 * 00394 WANTP = IJOB.EQ.1 .OR. IJOB.GE.4 00395 WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4 00396 WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5 00397 WANTD = WANTD1 .OR. WANTD2 00398 * 00399 * Set M to the dimension of the specified pair of deflating 00400 * subspaces. 00401 * 00402 M = 0 00403 PAIR = .FALSE. 00404 DO 10 K = 1, N 00405 IF( PAIR ) THEN 00406 PAIR = .FALSE. 00407 ELSE 00408 IF( K.LT.N ) THEN 00409 IF( A( K+1, K ).EQ.ZERO ) THEN 00410 IF( SELECT( K ) ) 00411 $ M = M + 1 00412 ELSE 00413 PAIR = .TRUE. 00414 IF( SELECT( K ) .OR. SELECT( K+1 ) ) 00415 $ M = M + 2 00416 END IF 00417 ELSE 00418 IF( SELECT( N ) ) 00419 $ M = M + 1 00420 END IF 00421 END IF 00422 10 CONTINUE 00423 * 00424 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 00425 LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) ) 00426 LIWMIN = MAX( 1, N+6 ) 00427 ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN 00428 LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) ) 00429 LIWMIN = MAX( 1, 2*M*( N-M ), N+6 ) 00430 ELSE 00431 LWMIN = MAX( 1, 4*N+16 ) 00432 LIWMIN = 1 00433 END IF 00434 * 00435 WORK( 1 ) = LWMIN 00436 IWORK( 1 ) = LIWMIN 00437 * 00438 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00439 INFO = -22 00440 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00441 INFO = -24 00442 END IF 00443 * 00444 IF( INFO.NE.0 ) THEN 00445 CALL XERBLA( 'DTGSEN', -INFO ) 00446 RETURN 00447 ELSE IF( LQUERY ) THEN 00448 RETURN 00449 END IF 00450 * 00451 * Quick return if possible. 00452 * 00453 IF( M.EQ.N .OR. M.EQ.0 ) THEN 00454 IF( WANTP ) THEN 00455 PL = ONE 00456 PR = ONE 00457 END IF 00458 IF( WANTD ) THEN 00459 DSCALE = ZERO 00460 DSUM = ONE 00461 DO 20 I = 1, N 00462 CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM ) 00463 CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM ) 00464 20 CONTINUE 00465 DIF( 1 ) = DSCALE*SQRT( DSUM ) 00466 DIF( 2 ) = DIF( 1 ) 00467 END IF 00468 GO TO 60 00469 END IF 00470 * 00471 * Collect the selected blocks at the top-left corner of (A, B). 00472 * 00473 KS = 0 00474 PAIR = .FALSE. 00475 DO 30 K = 1, N 00476 IF( PAIR ) THEN 00477 PAIR = .FALSE. 00478 ELSE 00479 * 00480 SWAP = SELECT( K ) 00481 IF( K.LT.N ) THEN 00482 IF( A( K+1, K ).NE.ZERO ) THEN 00483 PAIR = .TRUE. 00484 SWAP = SWAP .OR. SELECT( K+1 ) 00485 END IF 00486 END IF 00487 * 00488 IF( SWAP ) THEN 00489 KS = KS + 1 00490 * 00491 * Swap the K-th block to position KS. 00492 * Perform the reordering of diagonal blocks in (A, B) 00493 * by orthogonal transformation matrices and update 00494 * Q and Z accordingly (if requested): 00495 * 00496 KK = K 00497 IF( K.NE.KS ) 00498 $ CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00499 $ Z, LDZ, KK, KS, WORK, LWORK, IERR ) 00500 * 00501 IF( IERR.GT.0 ) THEN 00502 * 00503 * Swap is rejected: exit. 00504 * 00505 INFO = 1 00506 IF( WANTP ) THEN 00507 PL = ZERO 00508 PR = ZERO 00509 END IF 00510 IF( WANTD ) THEN 00511 DIF( 1 ) = ZERO 00512 DIF( 2 ) = ZERO 00513 END IF 00514 GO TO 60 00515 END IF 00516 * 00517 IF( PAIR ) 00518 $ KS = KS + 1 00519 END IF 00520 END IF 00521 30 CONTINUE 00522 IF( WANTP ) THEN 00523 * 00524 * Solve generalized Sylvester equation for R and L 00525 * and compute PL and PR. 00526 * 00527 N1 = M 00528 N2 = N - M 00529 I = N1 + 1 00530 IJB = 0 00531 CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 ) 00532 CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ), 00533 $ N1 ) 00534 CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 00535 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1, 00536 $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 00537 $ LWORK-2*N1*N2, IWORK, IERR ) 00538 * 00539 * Estimate the reciprocal of norms of "projections" onto left 00540 * and right eigenspaces. 00541 * 00542 RDSCAL = ZERO 00543 DSUM = ONE 00544 CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM ) 00545 PL = RDSCAL*SQRT( DSUM ) 00546 IF( PL.EQ.ZERO ) THEN 00547 PL = ONE 00548 ELSE 00549 PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) ) 00550 END IF 00551 RDSCAL = ZERO 00552 DSUM = ONE 00553 CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM ) 00554 PR = RDSCAL*SQRT( DSUM ) 00555 IF( PR.EQ.ZERO ) THEN 00556 PR = ONE 00557 ELSE 00558 PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) ) 00559 END IF 00560 END IF 00561 * 00562 IF( WANTD ) THEN 00563 * 00564 * Compute estimates of Difu and Difl. 00565 * 00566 IF( WANTD1 ) THEN 00567 N1 = M 00568 N2 = N - M 00569 I = N1 + 1 00570 IJB = IDIFJB 00571 * 00572 * Frobenius norm-based Difu-estimate. 00573 * 00574 CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 00575 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), 00576 $ N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ), 00577 $ LWORK-2*N1*N2, IWORK, IERR ) 00578 * 00579 * Frobenius norm-based Difl-estimate. 00580 * 00581 CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK, 00582 $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ), 00583 $ N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ), 00584 $ LWORK-2*N1*N2, IWORK, IERR ) 00585 ELSE 00586 * 00587 * 00588 * Compute 1-norm-based estimates of Difu and Difl using 00589 * reversed communication with DLACN2. In each step a 00590 * generalized Sylvester equation or a transposed variant 00591 * is solved. 00592 * 00593 KASE = 0 00594 N1 = M 00595 N2 = N - M 00596 I = N1 + 1 00597 IJB = 0 00598 MN2 = 2*N1*N2 00599 * 00600 * 1-norm-based estimate of Difu. 00601 * 00602 40 CONTINUE 00603 CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ), 00604 $ KASE, ISAVE ) 00605 IF( KASE.NE.0 ) THEN 00606 IF( KASE.EQ.1 ) THEN 00607 * 00608 * Solve generalized Sylvester equation. 00609 * 00610 CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, 00611 $ WORK, N1, B, LDB, B( I, I ), LDB, 00612 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 00613 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00614 $ IERR ) 00615 ELSE 00616 * 00617 * Solve the transposed variant. 00618 * 00619 CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA, 00620 $ WORK, N1, B, LDB, B( I, I ), LDB, 00621 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 00622 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00623 $ IERR ) 00624 END IF 00625 GO TO 40 00626 END IF 00627 DIF( 1 ) = DSCALE / DIF( 1 ) 00628 * 00629 * 1-norm-based estimate of Difl. 00630 * 00631 50 CONTINUE 00632 CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ), 00633 $ KASE, ISAVE ) 00634 IF( KASE.NE.0 ) THEN 00635 IF( KASE.EQ.1 ) THEN 00636 * 00637 * Solve generalized Sylvester equation. 00638 * 00639 CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, 00640 $ WORK, N2, B( I, I ), LDB, B, LDB, 00641 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 00642 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00643 $ IERR ) 00644 ELSE 00645 * 00646 * Solve the transposed variant. 00647 * 00648 CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA, 00649 $ WORK, N2, B( I, I ), LDB, B, LDB, 00650 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 00651 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00652 $ IERR ) 00653 END IF 00654 GO TO 50 00655 END IF 00656 DIF( 2 ) = DSCALE / DIF( 2 ) 00657 * 00658 END IF 00659 END IF 00660 * 00661 60 CONTINUE 00662 * 00663 * Compute generalized eigenvalues of reordered pair (A, B) and 00664 * normalize the generalized Schur form. 00665 * 00666 PAIR = .FALSE. 00667 DO 80 K = 1, N 00668 IF( PAIR ) THEN 00669 PAIR = .FALSE. 00670 ELSE 00671 * 00672 IF( K.LT.N ) THEN 00673 IF( A( K+1, K ).NE.ZERO ) THEN 00674 PAIR = .TRUE. 00675 END IF 00676 END IF 00677 * 00678 IF( PAIR ) THEN 00679 * 00680 * Compute the eigenvalue(s) at position K. 00681 * 00682 WORK( 1 ) = A( K, K ) 00683 WORK( 2 ) = A( K+1, K ) 00684 WORK( 3 ) = A( K, K+1 ) 00685 WORK( 4 ) = A( K+1, K+1 ) 00686 WORK( 5 ) = B( K, K ) 00687 WORK( 6 ) = B( K+1, K ) 00688 WORK( 7 ) = B( K, K+1 ) 00689 WORK( 8 ) = B( K+1, K+1 ) 00690 CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ), 00691 $ BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ), 00692 $ ALPHAI( K ) ) 00693 ALPHAI( K+1 ) = -ALPHAI( K ) 00694 * 00695 ELSE 00696 * 00697 IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN 00698 * 00699 * If B(K,K) is negative, make it positive 00700 * 00701 DO 70 I = 1, N 00702 A( K, I ) = -A( K, I ) 00703 B( K, I ) = -B( K, I ) 00704 IF( WANTQ ) Q( I, K ) = -Q( I, K ) 00705 70 CONTINUE 00706 END IF 00707 * 00708 ALPHAR( K ) = A( K, K ) 00709 ALPHAI( K ) = ZERO 00710 BETA( K ) = B( K, K ) 00711 * 00712 END IF 00713 END IF 00714 80 CONTINUE 00715 * 00716 WORK( 1 ) = LWMIN 00717 IWORK( 1 ) = LIWMIN 00718 * 00719 RETURN 00720 * 00721 * End of DTGSEN 00722 * 00723 END
1.7.2