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