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