|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZTGSEN( 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 ZLACN2 in place of ZLACON, 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 DOUBLE PRECISION PL, PR 00017 * .. 00018 * .. Array Arguments .. 00019 LOGICAL SELECT( * ) 00020 INTEGER IWORK( * ) 00021 DOUBLE PRECISION DIF( * ) 00022 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), 00023 $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * ) 00024 * .. 00025 * 00026 * Purpose 00027 * ======= 00028 * 00029 * ZTGSEN 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 * ZTGSEN 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*16 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*16 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*16 array, dimension (N) 00105 * BETA (output) COMPLEX*16 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*16 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*16 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) DOUBLE PRECISION 00140 * PR (output) DOUBLE PRECISION 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) DOUBLE PRECISION 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 ZLACN2. 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*16 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 * ZTGSEN 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 ZLATDF), then the parameter 00292 * IDIFJB (see below) should be changed from 3 to 4 (routine ZLATDF 00293 * (IJOB = 2 will be used)). See ZTGSYL 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 DOUBLE PRECISION ZERO, ONE 00328 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+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 DOUBLE PRECISION DSCALE, DSUM, RDSCAL, SAFMIN 00335 COMPLEX*16 TEMP1, TEMP2 00336 * .. 00337 * .. Local Arrays .. 00338 INTEGER ISAVE( 3 ) 00339 * .. 00340 * .. External Subroutines .. 00341 EXTERNAL XERBLA, ZLACN2, ZLACPY, ZLASSQ, ZSCAL, ZTGEXC, 00342 $ ZTGSYL 00343 * .. 00344 * .. Intrinsic Functions .. 00345 INTRINSIC ABS, DCMPLX, DCONJG, MAX, SQRT 00346 * .. 00347 * .. External Functions .. 00348 DOUBLE PRECISION DLAMCH 00349 EXTERNAL DLAMCH 00350 * .. 00351 * .. Executable Statements .. 00352 * 00353 * Decode and test the input parameters 00354 * 00355 INFO = 0 00356 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00357 * 00358 IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN 00359 INFO = -1 00360 ELSE IF( N.LT.0 ) THEN 00361 INFO = -5 00362 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00363 INFO = -7 00364 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00365 INFO = -9 00366 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00367 INFO = -13 00368 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00369 INFO = -15 00370 END IF 00371 * 00372 IF( INFO.NE.0 ) THEN 00373 CALL XERBLA( 'ZTGSEN', -INFO ) 00374 RETURN 00375 END IF 00376 * 00377 IERR = 0 00378 * 00379 WANTP = IJOB.EQ.1 .OR. IJOB.GE.4 00380 WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4 00381 WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5 00382 WANTD = WANTD1 .OR. WANTD2 00383 * 00384 * Set M to the dimension of the specified pair of deflating 00385 * subspaces. 00386 * 00387 M = 0 00388 DO 10 K = 1, N 00389 ALPHA( K ) = A( K, K ) 00390 BETA( K ) = B( K, K ) 00391 IF( K.LT.N ) THEN 00392 IF( SELECT( K ) ) 00393 $ M = M + 1 00394 ELSE 00395 IF( SELECT( N ) ) 00396 $ M = M + 1 00397 END IF 00398 10 CONTINUE 00399 * 00400 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 00401 LWMIN = MAX( 1, 2*M*( N-M ) ) 00402 LIWMIN = MAX( 1, N+2 ) 00403 ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN 00404 LWMIN = MAX( 1, 4*M*( N-M ) ) 00405 LIWMIN = MAX( 1, 2*M*( N-M ), N+2 ) 00406 ELSE 00407 LWMIN = 1 00408 LIWMIN = 1 00409 END IF 00410 * 00411 WORK( 1 ) = LWMIN 00412 IWORK( 1 ) = LIWMIN 00413 * 00414 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00415 INFO = -21 00416 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00417 INFO = -23 00418 END IF 00419 * 00420 IF( INFO.NE.0 ) THEN 00421 CALL XERBLA( 'ZTGSEN', -INFO ) 00422 RETURN 00423 ELSE IF( LQUERY ) THEN 00424 RETURN 00425 END IF 00426 * 00427 * Quick return if possible. 00428 * 00429 IF( M.EQ.N .OR. M.EQ.0 ) THEN 00430 IF( WANTP ) THEN 00431 PL = ONE 00432 PR = ONE 00433 END IF 00434 IF( WANTD ) THEN 00435 DSCALE = ZERO 00436 DSUM = ONE 00437 DO 20 I = 1, N 00438 CALL ZLASSQ( N, A( 1, I ), 1, DSCALE, DSUM ) 00439 CALL ZLASSQ( N, B( 1, I ), 1, DSCALE, DSUM ) 00440 20 CONTINUE 00441 DIF( 1 ) = DSCALE*SQRT( DSUM ) 00442 DIF( 2 ) = DIF( 1 ) 00443 END IF 00444 GO TO 70 00445 END IF 00446 * 00447 * Get machine constant 00448 * 00449 SAFMIN = DLAMCH( 'S' ) 00450 * 00451 * Collect the selected blocks at the top-left corner of (A, B). 00452 * 00453 KS = 0 00454 DO 30 K = 1, N 00455 SWAP = SELECT( K ) 00456 IF( SWAP ) THEN 00457 KS = KS + 1 00458 * 00459 * Swap the K-th block to position KS. Compute unitary Q 00460 * and Z that will swap adjacent diagonal blocks in (A, B). 00461 * 00462 IF( K.NE.KS ) 00463 $ CALL ZTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00464 $ LDZ, K, KS, IERR ) 00465 * 00466 IF( IERR.GT.0 ) THEN 00467 * 00468 * Swap is rejected: exit. 00469 * 00470 INFO = 1 00471 IF( WANTP ) THEN 00472 PL = ZERO 00473 PR = ZERO 00474 END IF 00475 IF( WANTD ) THEN 00476 DIF( 1 ) = ZERO 00477 DIF( 2 ) = ZERO 00478 END IF 00479 GO TO 70 00480 END IF 00481 END IF 00482 30 CONTINUE 00483 IF( WANTP ) THEN 00484 * 00485 * Solve generalized Sylvester equation for R and L: 00486 * A11 * R - L * A22 = A12 00487 * B11 * R - L * B22 = B12 00488 * 00489 N1 = M 00490 N2 = N - M 00491 I = N1 + 1 00492 CALL ZLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 ) 00493 CALL ZLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ), 00494 $ N1 ) 00495 IJB = 0 00496 CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 00497 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1, 00498 $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 00499 $ LWORK-2*N1*N2, IWORK, IERR ) 00500 * 00501 * Estimate the reciprocal of norms of "projections" onto 00502 * left and right eigenspaces 00503 * 00504 RDSCAL = ZERO 00505 DSUM = ONE 00506 CALL ZLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM ) 00507 PL = RDSCAL*SQRT( DSUM ) 00508 IF( PL.EQ.ZERO ) THEN 00509 PL = ONE 00510 ELSE 00511 PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) ) 00512 END IF 00513 RDSCAL = ZERO 00514 DSUM = ONE 00515 CALL ZLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM ) 00516 PR = RDSCAL*SQRT( DSUM ) 00517 IF( PR.EQ.ZERO ) THEN 00518 PR = ONE 00519 ELSE 00520 PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) ) 00521 END IF 00522 END IF 00523 IF( WANTD ) THEN 00524 * 00525 * Compute estimates Difu and Difl. 00526 * 00527 IF( WANTD1 ) THEN 00528 N1 = M 00529 N2 = N - M 00530 I = N1 + 1 00531 IJB = IDIFJB 00532 * 00533 * Frobenius norm-based Difu estimate. 00534 * 00535 CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 00536 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), 00537 $ N1, DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 00538 $ LWORK-2*N1*N2, IWORK, IERR ) 00539 * 00540 * Frobenius norm-based Difl estimate. 00541 * 00542 CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK, 00543 $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ), 00544 $ N2, DSCALE, DIF( 2 ), WORK( N1*N2*2+1 ), 00545 $ LWORK-2*N1*N2, IWORK, IERR ) 00546 ELSE 00547 * 00548 * Compute 1-norm-based estimates of Difu and Difl using 00549 * reversed communication with ZLACN2. In each step a 00550 * generalized Sylvester equation or a transposed variant 00551 * is solved. 00552 * 00553 KASE = 0 00554 N1 = M 00555 N2 = N - M 00556 I = N1 + 1 00557 IJB = 0 00558 MN2 = 2*N1*N2 00559 * 00560 * 1-norm-based estimate of Difu. 00561 * 00562 40 CONTINUE 00563 CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 1 ), KASE, 00564 $ ISAVE ) 00565 IF( KASE.NE.0 ) THEN 00566 IF( KASE.EQ.1 ) THEN 00567 * 00568 * Solve generalized Sylvester equation 00569 * 00570 CALL ZTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, 00571 $ WORK, N1, B, LDB, B( I, I ), LDB, 00572 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 00573 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 00574 $ IERR ) 00575 ELSE 00576 * 00577 * Solve the transposed variant. 00578 * 00579 CALL ZTGSYL( 'C', IJB, N1, N2, A, LDA, A( I, I ), LDA, 00580 $ WORK, N1, B, LDB, B( I, I ), LDB, 00581 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 00582 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 00583 $ IERR ) 00584 END IF 00585 GO TO 40 00586 END IF 00587 DIF( 1 ) = DSCALE / DIF( 1 ) 00588 * 00589 * 1-norm-based estimate of Difl. 00590 * 00591 50 CONTINUE 00592 CALL ZLACN2( MN2, WORK( MN2+1 ), WORK, DIF( 2 ), KASE, 00593 $ ISAVE ) 00594 IF( KASE.NE.0 ) THEN 00595 IF( KASE.EQ.1 ) THEN 00596 * 00597 * Solve generalized Sylvester equation 00598 * 00599 CALL ZTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, 00600 $ WORK, N2, B( I, I ), LDB, B, LDB, 00601 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 00602 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 00603 $ IERR ) 00604 ELSE 00605 * 00606 * Solve the transposed variant. 00607 * 00608 CALL ZTGSYL( 'C', IJB, N2, N1, A( I, I ), LDA, A, LDA, 00609 $ WORK, N2, B, LDB, B( I, I ), LDB, 00610 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 00611 $ WORK( N1*N2*2+1 ), LWORK-2*N1*N2, IWORK, 00612 $ IERR ) 00613 END IF 00614 GO TO 50 00615 END IF 00616 DIF( 2 ) = DSCALE / DIF( 2 ) 00617 END IF 00618 END IF 00619 * 00620 * If B(K,K) is complex, make it real and positive (normalization 00621 * of the generalized Schur form) and Store the generalized 00622 * eigenvalues of reordered pair (A, B) 00623 * 00624 DO 60 K = 1, N 00625 DSCALE = ABS( B( K, K ) ) 00626 IF( DSCALE.GT.SAFMIN ) THEN 00627 TEMP1 = DCONJG( B( K, K ) / DSCALE ) 00628 TEMP2 = B( K, K ) / DSCALE 00629 B( K, K ) = DSCALE 00630 CALL ZSCAL( N-K, TEMP1, B( K, K+1 ), LDB ) 00631 CALL ZSCAL( N-K+1, TEMP1, A( K, K ), LDA ) 00632 IF( WANTQ ) 00633 $ CALL ZSCAL( N, TEMP2, Q( 1, K ), 1 ) 00634 ELSE 00635 B( K, K ) = DCMPLX( ZERO, ZERO ) 00636 END IF 00637 * 00638 ALPHA( K ) = A( K, K ) 00639 BETA( K ) = B( K, K ) 00640 * 00641 60 CONTINUE 00642 * 00643 70 CONTINUE 00644 * 00645 WORK( 1 ) = LWMIN 00646 IWORK( 1 ) = LIWMIN 00647 * 00648 RETURN 00649 * 00650 * End of ZTGSEN 00651 * 00652 END
1.7.3 
