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