![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTGSEN 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTGSEN + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsen.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsen.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsen.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTGSEN( 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 * DOUBLE PRECISION PL, PR 00030 * .. 00031 * .. Array Arguments .. 00032 * LOGICAL SELECT( * ) 00033 * INTEGER IWORK( * ) 00034 * DOUBLE PRECISION 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 *> DTGSEN 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 DGGES), 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 *> DTGSEN 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, DTGSEN 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 00156 *> \endverbatim 00157 *> 00158 *> \param[out] ALPHAI 00159 *> \verbatim 00160 *> ALPHAI is DOUBLE PRECISION array, dimension (N) 00161 *> \endverbatim 00162 *> 00163 *> \param[out] BETA 00164 *> \verbatim 00165 *> BETA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 00224 *> \endverbatim 00225 00226 *> \param[out] PR 00227 *> \verbatim 00228 *> PR is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, 00252 *> dimension (MAX(1,LWORK)) 00253 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00254 *> \endverbatim 00255 *> 00256 *> \param[in] LWORK 00257 *> \verbatim 00258 *> LWORK is INTEGER 00259 *> The dimension of the array WORK. LWORK >= 4*N+16. 00260 *> If IJOB = 1, 2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)). 00261 *> If IJOB = 3 or 5, LWORK >= MAX(4*N+16, 4*M*(N-M)). 00262 *> 00263 *> If LWORK = -1, then a workspace query is assumed; the routine 00264 *> only calculates the optimal size of the WORK array, returns 00265 *> this value as the first entry of the WORK array, and no error 00266 *> message related to LWORK is issued by XERBLA. 00267 *> \endverbatim 00268 *> 00269 *> \param[out] IWORK 00270 *> \verbatim 00271 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00272 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00273 *> \endverbatim 00274 *> 00275 *> \param[in] LIWORK 00276 *> \verbatim 00277 *> LIWORK is INTEGER 00278 *> The dimension of the array IWORK. LIWORK >= 1. 00279 *> If IJOB = 1, 2 or 4, LIWORK >= N+6. 00280 *> If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M), N+6). 00281 *> 00282 *> If LIWORK = -1, then a workspace query is assumed; the 00283 *> routine only calculates the optimal size of the IWORK array, 00284 *> returns this value as the first entry of the IWORK array, and 00285 *> no error message related to LIWORK is issued by XERBLA. 00286 *> \endverbatim 00287 *> 00288 *> \param[out] INFO 00289 *> \verbatim 00290 *> INFO is INTEGER 00291 *> =0: Successful exit. 00292 *> <0: If INFO = -i, the i-th argument had an illegal value. 00293 *> =1: Reordering of (A, B) failed because the transformed 00294 *> matrix pair (A, B) would be too far from generalized 00295 *> Schur form; the problem is very ill-conditioned. 00296 *> (A, B) may have been partially reordered. 00297 *> If requested, 0 is returned in DIF(*), PL and PR. 00298 *> \endverbatim 00299 * 00300 * Authors: 00301 * ======== 00302 * 00303 *> \author Univ. of Tennessee 00304 *> \author Univ. of California Berkeley 00305 *> \author Univ. of Colorado Denver 00306 *> \author NAG Ltd. 00307 * 00308 *> \date November 2011 00309 * 00310 *> \ingroup doubleOTHERcomputational 00311 * 00312 *> \par Further Details: 00313 * ===================== 00314 *> 00315 *> \verbatim 00316 *> 00317 *> DTGSEN first collects the selected eigenvalues by computing 00318 *> orthogonal U and W that move them to the top left corner of (A, B). 00319 *> In other words, the selected eigenvalues are the eigenvalues of 00320 *> (A11, B11) in: 00321 *> 00322 *> U**T*(A, B)*W = (A11 A12) (B11 B12) n1 00323 *> ( 0 A22),( 0 B22) n2 00324 *> n1 n2 n1 n2 00325 *> 00326 *> where N = n1+n2 and U**T means the transpose of U. The first n1 columns 00327 *> of U and W span the specified pair of left and right eigenspaces 00328 *> (deflating subspaces) of (A, B). 00329 *> 00330 *> If (A, B) has been obtained from the generalized real Schur 00331 *> decomposition of a matrix pair (C, D) = Q*(A, B)*Z**T, then the 00332 *> reordered generalized real Schur form of (C, D) is given by 00333 *> 00334 *> (C, D) = (Q*U)*(U**T*(A, B)*W)*(Z*W)**T, 00335 *> 00336 *> and the first n1 columns of Q*U and Z*W span the corresponding 00337 *> deflating subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.). 00338 *> 00339 *> Note that if the selected eigenvalue is sufficiently ill-conditioned, 00340 *> then its value may differ significantly from its value before 00341 *> reordering. 00342 *> 00343 *> The reciprocal condition numbers of the left and right eigenspaces 00344 *> spanned by the first n1 columns of U and W (or Q*U and Z*W) may 00345 *> be returned in DIF(1:2), corresponding to Difu and Difl, resp. 00346 *> 00347 *> The Difu and Difl are defined as: 00348 *> 00349 *> Difu[(A11, B11), (A22, B22)] = sigma-min( Zu ) 00350 *> and 00351 *> Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)], 00352 *> 00353 *> where sigma-min(Zu) is the smallest singular value of the 00354 *> (2*n1*n2)-by-(2*n1*n2) matrix 00355 *> 00356 *> Zu = [ kron(In2, A11) -kron(A22**T, In1) ] 00357 *> [ kron(In2, B11) -kron(B22**T, In1) ]. 00358 *> 00359 *> Here, Inx is the identity matrix of size nx and A22**T is the 00360 *> transpose of A22. kron(X, Y) is the Kronecker product between 00361 *> the matrices X and Y. 00362 *> 00363 *> When DIF(2) is small, small changes in (A, B) can cause large changes 00364 *> in the deflating subspace. An approximate (asymptotic) bound on the 00365 *> maximum angular error in the computed deflating subspaces is 00366 *> 00367 *> EPS * norm((A, B)) / DIF(2), 00368 *> 00369 *> where EPS is the machine precision. 00370 *> 00371 *> The reciprocal norm of the projectors on the left and right 00372 *> eigenspaces associated with (A11, B11) may be returned in PL and PR. 00373 *> They are computed as follows. First we compute L and R so that 00374 *> P*(A, B)*Q is block diagonal, where 00375 *> 00376 *> P = ( I -L ) n1 Q = ( I R ) n1 00377 *> ( 0 I ) n2 and ( 0 I ) n2 00378 *> n1 n2 n1 n2 00379 *> 00380 *> and (L, R) is the solution to the generalized Sylvester equation 00381 *> 00382 *> A11*R - L*A22 = -A12 00383 *> B11*R - L*B22 = -B12 00384 *> 00385 *> Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2). 00386 *> An approximate (asymptotic) bound on the average absolute error of 00387 *> the selected eigenvalues is 00388 *> 00389 *> EPS * norm((A, B)) / PL. 00390 *> 00391 *> There are also global error bounds which valid for perturbations up 00392 *> to a certain restriction: A lower bound (x) on the smallest 00393 *> F-norm(E,F) for which an eigenvalue of (A11, B11) may move and 00394 *> coalesce with an eigenvalue of (A22, B22) under perturbation (E,F), 00395 *> (i.e. (A + E, B + F), is 00396 *> 00397 *> x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)). 00398 *> 00399 *> An approximate bound on x can be computed from DIF(1:2), PL and PR. 00400 *> 00401 *> If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed 00402 *> (L', R') and unperturbed (L, R) left and right deflating subspaces 00403 *> associated with the selected cluster in the (1,1)-blocks can be 00404 *> bounded as 00405 *> 00406 *> max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2)) 00407 *> max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2)) 00408 *> 00409 *> See LAPACK User's Guide section 4.11 or the following references 00410 *> for more information. 00411 *> 00412 *> Note that if the default method for computing the Frobenius-norm- 00413 *> based estimate DIF is not wanted (see DLATDF), then the parameter 00414 *> IDIFJB (see below) should be changed from 3 to 4 (routine DLATDF 00415 *> (IJOB = 2 will be used)). See DTGSYL for more details. 00416 *> \endverbatim 00417 * 00418 *> \par Contributors: 00419 * ================== 00420 *> 00421 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00422 *> Umea University, S-901 87 Umea, Sweden. 00423 * 00424 *> \par References: 00425 * ================ 00426 *> 00427 *> \verbatim 00428 *> 00429 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00430 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00431 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00432 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00433 *> 00434 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00435 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00436 *> Estimation: Theory, Algorithms and Software, 00437 *> Report UMINF - 94.04, Department of Computing Science, Umea 00438 *> University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working 00439 *> Note 87. To appear in Numerical Algorithms, 1996. 00440 *> 00441 *> [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 00442 *> for Solving the Generalized Sylvester Equation and Estimating the 00443 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23, 00444 *> Department of Computing Science, Umea University, S-901 87 Umea, 00445 *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working 00446 *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 00447 *> 1996. 00448 *> \endverbatim 00449 *> 00450 * ===================================================================== 00451 SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, 00452 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, 00453 $ PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO ) 00454 * 00455 * -- LAPACK computational routine (version 3.4.0) -- 00456 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00457 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00458 * November 2011 00459 * 00460 * .. Scalar Arguments .. 00461 LOGICAL WANTQ, WANTZ 00462 INTEGER IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, 00463 $ M, N 00464 DOUBLE PRECISION PL, PR 00465 * .. 00466 * .. Array Arguments .. 00467 LOGICAL SELECT( * ) 00468 INTEGER IWORK( * ) 00469 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 00470 $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), 00471 $ WORK( * ), Z( LDZ, * ) 00472 * .. 00473 * 00474 * ===================================================================== 00475 * 00476 * .. Parameters .. 00477 INTEGER IDIFJB 00478 PARAMETER ( IDIFJB = 3 ) 00479 DOUBLE PRECISION ZERO, ONE 00480 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00481 * .. 00482 * .. Local Scalars .. 00483 LOGICAL LQUERY, PAIR, SWAP, WANTD, WANTD1, WANTD2, 00484 $ WANTP 00485 INTEGER I, IERR, IJB, K, KASE, KK, KS, LIWMIN, LWMIN, 00486 $ MN2, N1, N2 00487 DOUBLE PRECISION DSCALE, DSUM, EPS, RDSCAL, SMLNUM 00488 * .. 00489 * .. Local Arrays .. 00490 INTEGER ISAVE( 3 ) 00491 * .. 00492 * .. External Subroutines .. 00493 EXTERNAL DLACN2, DLACPY, DLAG2, DLASSQ, DTGEXC, DTGSYL, 00494 $ XERBLA 00495 * .. 00496 * .. External Functions .. 00497 DOUBLE PRECISION DLAMCH 00498 EXTERNAL DLAMCH 00499 * .. 00500 * .. Intrinsic Functions .. 00501 INTRINSIC MAX, SIGN, SQRT 00502 * .. 00503 * .. Executable Statements .. 00504 * 00505 * Decode and test the input parameters 00506 * 00507 INFO = 0 00508 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00509 * 00510 IF( IJOB.LT.0 .OR. IJOB.GT.5 ) THEN 00511 INFO = -1 00512 ELSE IF( N.LT.0 ) THEN 00513 INFO = -5 00514 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00515 INFO = -7 00516 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00517 INFO = -9 00518 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00519 INFO = -14 00520 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00521 INFO = -16 00522 END IF 00523 * 00524 IF( INFO.NE.0 ) THEN 00525 CALL XERBLA( 'DTGSEN', -INFO ) 00526 RETURN 00527 END IF 00528 * 00529 * Get machine constants 00530 * 00531 EPS = DLAMCH( 'P' ) 00532 SMLNUM = DLAMCH( 'S' ) / EPS 00533 IERR = 0 00534 * 00535 WANTP = IJOB.EQ.1 .OR. IJOB.GE.4 00536 WANTD1 = IJOB.EQ.2 .OR. IJOB.EQ.4 00537 WANTD2 = IJOB.EQ.3 .OR. IJOB.EQ.5 00538 WANTD = WANTD1 .OR. WANTD2 00539 * 00540 * Set M to the dimension of the specified pair of deflating 00541 * subspaces. 00542 * 00543 M = 0 00544 PAIR = .FALSE. 00545 DO 10 K = 1, N 00546 IF( PAIR ) THEN 00547 PAIR = .FALSE. 00548 ELSE 00549 IF( K.LT.N ) THEN 00550 IF( A( K+1, K ).EQ.ZERO ) THEN 00551 IF( SELECT( K ) ) 00552 $ M = M + 1 00553 ELSE 00554 PAIR = .TRUE. 00555 IF( SELECT( K ) .OR. SELECT( K+1 ) ) 00556 $ M = M + 2 00557 END IF 00558 ELSE 00559 IF( SELECT( N ) ) 00560 $ M = M + 1 00561 END IF 00562 END IF 00563 10 CONTINUE 00564 * 00565 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 00566 LWMIN = MAX( 1, 4*N+16, 2*M*( N-M ) ) 00567 LIWMIN = MAX( 1, N+6 ) 00568 ELSE IF( IJOB.EQ.3 .OR. IJOB.EQ.5 ) THEN 00569 LWMIN = MAX( 1, 4*N+16, 4*M*( N-M ) ) 00570 LIWMIN = MAX( 1, 2*M*( N-M ), N+6 ) 00571 ELSE 00572 LWMIN = MAX( 1, 4*N+16 ) 00573 LIWMIN = 1 00574 END IF 00575 * 00576 WORK( 1 ) = LWMIN 00577 IWORK( 1 ) = LIWMIN 00578 * 00579 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00580 INFO = -22 00581 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00582 INFO = -24 00583 END IF 00584 * 00585 IF( INFO.NE.0 ) THEN 00586 CALL XERBLA( 'DTGSEN', -INFO ) 00587 RETURN 00588 ELSE IF( LQUERY ) THEN 00589 RETURN 00590 END IF 00591 * 00592 * Quick return if possible. 00593 * 00594 IF( M.EQ.N .OR. M.EQ.0 ) THEN 00595 IF( WANTP ) THEN 00596 PL = ONE 00597 PR = ONE 00598 END IF 00599 IF( WANTD ) THEN 00600 DSCALE = ZERO 00601 DSUM = ONE 00602 DO 20 I = 1, N 00603 CALL DLASSQ( N, A( 1, I ), 1, DSCALE, DSUM ) 00604 CALL DLASSQ( N, B( 1, I ), 1, DSCALE, DSUM ) 00605 20 CONTINUE 00606 DIF( 1 ) = DSCALE*SQRT( DSUM ) 00607 DIF( 2 ) = DIF( 1 ) 00608 END IF 00609 GO TO 60 00610 END IF 00611 * 00612 * Collect the selected blocks at the top-left corner of (A, B). 00613 * 00614 KS = 0 00615 PAIR = .FALSE. 00616 DO 30 K = 1, N 00617 IF( PAIR ) THEN 00618 PAIR = .FALSE. 00619 ELSE 00620 * 00621 SWAP = SELECT( K ) 00622 IF( K.LT.N ) THEN 00623 IF( A( K+1, K ).NE.ZERO ) THEN 00624 PAIR = .TRUE. 00625 SWAP = SWAP .OR. SELECT( K+1 ) 00626 END IF 00627 END IF 00628 * 00629 IF( SWAP ) THEN 00630 KS = KS + 1 00631 * 00632 * Swap the K-th block to position KS. 00633 * Perform the reordering of diagonal blocks in (A, B) 00634 * by orthogonal transformation matrices and update 00635 * Q and Z accordingly (if requested): 00636 * 00637 KK = K 00638 IF( K.NE.KS ) 00639 $ CALL DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00640 $ Z, LDZ, KK, KS, WORK, LWORK, IERR ) 00641 * 00642 IF( IERR.GT.0 ) THEN 00643 * 00644 * Swap is rejected: exit. 00645 * 00646 INFO = 1 00647 IF( WANTP ) THEN 00648 PL = ZERO 00649 PR = ZERO 00650 END IF 00651 IF( WANTD ) THEN 00652 DIF( 1 ) = ZERO 00653 DIF( 2 ) = ZERO 00654 END IF 00655 GO TO 60 00656 END IF 00657 * 00658 IF( PAIR ) 00659 $ KS = KS + 1 00660 END IF 00661 END IF 00662 30 CONTINUE 00663 IF( WANTP ) THEN 00664 * 00665 * Solve generalized Sylvester equation for R and L 00666 * and compute PL and PR. 00667 * 00668 N1 = M 00669 N2 = N - M 00670 I = N1 + 1 00671 IJB = 0 00672 CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 ) 00673 CALL DLACPY( 'Full', N1, N2, B( 1, I ), LDB, WORK( N1*N2+1 ), 00674 $ N1 ) 00675 CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 00676 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), N1, 00677 $ DSCALE, DIF( 1 ), WORK( N1*N2*2+1 ), 00678 $ LWORK-2*N1*N2, IWORK, IERR ) 00679 * 00680 * Estimate the reciprocal of norms of "projections" onto left 00681 * and right eigenspaces. 00682 * 00683 RDSCAL = ZERO 00684 DSUM = ONE 00685 CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM ) 00686 PL = RDSCAL*SQRT( DSUM ) 00687 IF( PL.EQ.ZERO ) THEN 00688 PL = ONE 00689 ELSE 00690 PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) ) 00691 END IF 00692 RDSCAL = ZERO 00693 DSUM = ONE 00694 CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM ) 00695 PR = RDSCAL*SQRT( DSUM ) 00696 IF( PR.EQ.ZERO ) THEN 00697 PR = ONE 00698 ELSE 00699 PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) ) 00700 END IF 00701 END IF 00702 * 00703 IF( WANTD ) THEN 00704 * 00705 * Compute estimates of Difu and Difl. 00706 * 00707 IF( WANTD1 ) THEN 00708 N1 = M 00709 N2 = N - M 00710 I = N1 + 1 00711 IJB = IDIFJB 00712 * 00713 * Frobenius norm-based Difu-estimate. 00714 * 00715 CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, WORK, 00716 $ N1, B, LDB, B( I, I ), LDB, WORK( N1*N2+1 ), 00717 $ N1, DSCALE, DIF( 1 ), WORK( 2*N1*N2+1 ), 00718 $ LWORK-2*N1*N2, IWORK, IERR ) 00719 * 00720 * Frobenius norm-based Difl-estimate. 00721 * 00722 CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, WORK, 00723 $ N2, B( I, I ), LDB, B, LDB, WORK( N1*N2+1 ), 00724 $ N2, DSCALE, DIF( 2 ), WORK( 2*N1*N2+1 ), 00725 $ LWORK-2*N1*N2, IWORK, IERR ) 00726 ELSE 00727 * 00728 * 00729 * Compute 1-norm-based estimates of Difu and Difl using 00730 * reversed communication with DLACN2. In each step a 00731 * generalized Sylvester equation or a transposed variant 00732 * is solved. 00733 * 00734 KASE = 0 00735 N1 = M 00736 N2 = N - M 00737 I = N1 + 1 00738 IJB = 0 00739 MN2 = 2*N1*N2 00740 * 00741 * 1-norm-based estimate of Difu. 00742 * 00743 40 CONTINUE 00744 CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 1 ), 00745 $ KASE, ISAVE ) 00746 IF( KASE.NE.0 ) THEN 00747 IF( KASE.EQ.1 ) THEN 00748 * 00749 * Solve generalized Sylvester equation. 00750 * 00751 CALL DTGSYL( 'N', IJB, N1, N2, A, LDA, A( I, I ), LDA, 00752 $ WORK, N1, B, LDB, B( I, I ), LDB, 00753 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 00754 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00755 $ IERR ) 00756 ELSE 00757 * 00758 * Solve the transposed variant. 00759 * 00760 CALL DTGSYL( 'T', IJB, N1, N2, A, LDA, A( I, I ), LDA, 00761 $ WORK, N1, B, LDB, B( I, I ), LDB, 00762 $ WORK( N1*N2+1 ), N1, DSCALE, DIF( 1 ), 00763 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00764 $ IERR ) 00765 END IF 00766 GO TO 40 00767 END IF 00768 DIF( 1 ) = DSCALE / DIF( 1 ) 00769 * 00770 * 1-norm-based estimate of Difl. 00771 * 00772 50 CONTINUE 00773 CALL DLACN2( MN2, WORK( MN2+1 ), WORK, IWORK, DIF( 2 ), 00774 $ KASE, ISAVE ) 00775 IF( KASE.NE.0 ) THEN 00776 IF( KASE.EQ.1 ) THEN 00777 * 00778 * Solve generalized Sylvester equation. 00779 * 00780 CALL DTGSYL( 'N', IJB, N2, N1, A( I, I ), LDA, A, LDA, 00781 $ WORK, N2, B( I, I ), LDB, B, LDB, 00782 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 00783 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00784 $ IERR ) 00785 ELSE 00786 * 00787 * Solve the transposed variant. 00788 * 00789 CALL DTGSYL( 'T', IJB, N2, N1, A( I, I ), LDA, A, LDA, 00790 $ WORK, N2, B( I, I ), LDB, B, LDB, 00791 $ WORK( N1*N2+1 ), N2, DSCALE, DIF( 2 ), 00792 $ WORK( 2*N1*N2+1 ), LWORK-2*N1*N2, IWORK, 00793 $ IERR ) 00794 END IF 00795 GO TO 50 00796 END IF 00797 DIF( 2 ) = DSCALE / DIF( 2 ) 00798 * 00799 END IF 00800 END IF 00801 * 00802 60 CONTINUE 00803 * 00804 * Compute generalized eigenvalues of reordered pair (A, B) and 00805 * normalize the generalized Schur form. 00806 * 00807 PAIR = .FALSE. 00808 DO 80 K = 1, N 00809 IF( PAIR ) THEN 00810 PAIR = .FALSE. 00811 ELSE 00812 * 00813 IF( K.LT.N ) THEN 00814 IF( A( K+1, K ).NE.ZERO ) THEN 00815 PAIR = .TRUE. 00816 END IF 00817 END IF 00818 * 00819 IF( PAIR ) THEN 00820 * 00821 * Compute the eigenvalue(s) at position K. 00822 * 00823 WORK( 1 ) = A( K, K ) 00824 WORK( 2 ) = A( K+1, K ) 00825 WORK( 3 ) = A( K, K+1 ) 00826 WORK( 4 ) = A( K+1, K+1 ) 00827 WORK( 5 ) = B( K, K ) 00828 WORK( 6 ) = B( K+1, K ) 00829 WORK( 7 ) = B( K, K+1 ) 00830 WORK( 8 ) = B( K+1, K+1 ) 00831 CALL DLAG2( WORK, 2, WORK( 5 ), 2, SMLNUM*EPS, BETA( K ), 00832 $ BETA( K+1 ), ALPHAR( K ), ALPHAR( K+1 ), 00833 $ ALPHAI( K ) ) 00834 ALPHAI( K+1 ) = -ALPHAI( K ) 00835 * 00836 ELSE 00837 * 00838 IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN 00839 * 00840 * If B(K,K) is negative, make it positive 00841 * 00842 DO 70 I = 1, N 00843 A( K, I ) = -A( K, I ) 00844 B( K, I ) = -B( K, I ) 00845 IF( WANTQ ) Q( I, K ) = -Q( I, K ) 00846 70 CONTINUE 00847 END IF 00848 * 00849 ALPHAR( K ) = A( K, K ) 00850 ALPHAI( K ) = ZERO 00851 BETA( K ) = B( K, K ) 00852 * 00853 END IF 00854 END IF 00855 80 CONTINUE 00856 * 00857 WORK( 1 ) = LWMIN 00858 IWORK( 1 ) = LIWMIN 00859 * 00860 RETURN 00861 * 00862 * End of DTGSEN 00863 * 00864 END