![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGGESX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggesx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggesx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggesx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, 00022 * B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, 00023 * LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, 00024 * IWORK, LIWORK, BWORK, INFO ) 00025 * 00026 * .. Scalar Arguments .. 00027 * CHARACTER JOBVSL, JOBVSR, SENSE, SORT 00028 * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N, 00029 * $ SDIM 00030 * .. 00031 * .. Array Arguments .. 00032 * LOGICAL BWORK( * ) 00033 * INTEGER IWORK( * ) 00034 * REAL RCONDE( 2 ), RCONDV( 2 ), RWORK( * ) 00035 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00036 * $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), 00037 * $ WORK( * ) 00038 * .. 00039 * .. Function Arguments .. 00040 * LOGICAL SELCTG 00041 * EXTERNAL SELCTG 00042 * .. 00043 * 00044 * 00045 *> \par Purpose: 00046 * ============= 00047 *> 00048 *> \verbatim 00049 *> 00050 *> CGGESX computes for a pair of N-by-N complex nonsymmetric matrices 00051 *> (A,B), the generalized eigenvalues, the complex Schur form (S,T), 00052 *> and, optionally, the left and/or right matrices of Schur vectors (VSL 00053 *> and VSR). This gives the generalized Schur factorization 00054 *> 00055 *> (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H ) 00056 *> 00057 *> where (VSR)**H is the conjugate-transpose of VSR. 00058 *> 00059 *> Optionally, it also orders the eigenvalues so that a selected cluster 00060 *> of eigenvalues appears in the leading diagonal blocks of the upper 00061 *> triangular matrix S and the upper triangular matrix T; computes 00062 *> a reciprocal condition number for the average of the selected 00063 *> eigenvalues (RCONDE); and computes a reciprocal condition number for 00064 *> the right and left deflating subspaces corresponding to the selected 00065 *> eigenvalues (RCONDV). The leading columns of VSL and VSR then form 00066 *> an orthonormal basis for the corresponding left and right eigenspaces 00067 *> (deflating subspaces). 00068 *> 00069 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w 00070 *> or a ratio alpha/beta = w, such that A - w*B is singular. It is 00071 *> usually represented as the pair (alpha,beta), as there is a 00072 *> reasonable interpretation for beta=0 or for both being zero. 00073 *> 00074 *> A pair of matrices (S,T) is in generalized complex Schur form if T is 00075 *> upper triangular with non-negative diagonal and S is upper 00076 *> triangular. 00077 *> \endverbatim 00078 * 00079 * Arguments: 00080 * ========== 00081 * 00082 *> \param[in] JOBVSL 00083 *> \verbatim 00084 *> JOBVSL is CHARACTER*1 00085 *> = 'N': do not compute the left Schur vectors; 00086 *> = 'V': compute the left Schur vectors. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] JOBVSR 00090 *> \verbatim 00091 *> JOBVSR is CHARACTER*1 00092 *> = 'N': do not compute the right Schur vectors; 00093 *> = 'V': compute the right Schur vectors. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] SORT 00097 *> \verbatim 00098 *> SORT is CHARACTER*1 00099 *> Specifies whether or not to order the eigenvalues on the 00100 *> diagonal of the generalized Schur form. 00101 *> = 'N': Eigenvalues are not ordered; 00102 *> = 'S': Eigenvalues are ordered (see SELCTG). 00103 *> \endverbatim 00104 *> 00105 *> \param[in] SELCTG 00106 *> \verbatim 00107 *> SELCTG is procedure) LOGICAL FUNCTION of two COMPLEX arguments 00108 *> SELCTG must be declared EXTERNAL in the calling subroutine. 00109 *> If SORT = 'N', SELCTG is not referenced. 00110 *> If SORT = 'S', SELCTG is used to select eigenvalues to sort 00111 *> to the top left of the Schur form. 00112 *> Note that a selected complex eigenvalue may no longer satisfy 00113 *> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since 00114 *> ordering may change the value of complex eigenvalues 00115 *> (especially if the eigenvalue is ill-conditioned), in this 00116 *> case INFO is set to N+3 see INFO below). 00117 *> \endverbatim 00118 *> 00119 *> \param[in] SENSE 00120 *> \verbatim 00121 *> SENSE is CHARACTER*1 00122 *> Determines which reciprocal condition numbers are computed. 00123 *> = 'N' : None are computed; 00124 *> = 'E' : Computed for average of selected eigenvalues only; 00125 *> = 'V' : Computed for selected deflating subspaces only; 00126 *> = 'B' : Computed for both. 00127 *> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'. 00128 *> \endverbatim 00129 *> 00130 *> \param[in] N 00131 *> \verbatim 00132 *> N is INTEGER 00133 *> The order of the matrices A, B, VSL, and VSR. N >= 0. 00134 *> \endverbatim 00135 *> 00136 *> \param[in,out] A 00137 *> \verbatim 00138 *> A is COMPLEX array, dimension (LDA, N) 00139 *> On entry, the first of the pair of matrices. 00140 *> On exit, A has been overwritten by its generalized Schur 00141 *> form S. 00142 *> \endverbatim 00143 *> 00144 *> \param[in] LDA 00145 *> \verbatim 00146 *> LDA is INTEGER 00147 *> The leading dimension of A. LDA >= max(1,N). 00148 *> \endverbatim 00149 *> 00150 *> \param[in,out] B 00151 *> \verbatim 00152 *> B is COMPLEX array, dimension (LDB, N) 00153 *> On entry, the second of the pair of matrices. 00154 *> On exit, B has been overwritten by its generalized Schur 00155 *> form T. 00156 *> \endverbatim 00157 *> 00158 *> \param[in] LDB 00159 *> \verbatim 00160 *> LDB is INTEGER 00161 *> The leading dimension of B. LDB >= max(1,N). 00162 *> \endverbatim 00163 *> 00164 *> \param[out] SDIM 00165 *> \verbatim 00166 *> SDIM is INTEGER 00167 *> If SORT = 'N', SDIM = 0. 00168 *> If SORT = 'S', SDIM = number of eigenvalues (after sorting) 00169 *> for which SELCTG is true. 00170 *> \endverbatim 00171 *> 00172 *> \param[out] ALPHA 00173 *> \verbatim 00174 *> ALPHA is COMPLEX array, dimension (N) 00175 *> \endverbatim 00176 *> 00177 *> \param[out] BETA 00178 *> \verbatim 00179 *> BETA is COMPLEX array, dimension (N) 00180 *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the 00181 *> generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are 00182 *> the diagonals of the complex Schur form (S,T). BETA(j) will 00183 *> be non-negative real. 00184 *> 00185 *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or 00186 *> underflow, and BETA(j) may even be zero. Thus, the user 00187 *> should avoid naively computing the ratio alpha/beta. 00188 *> However, ALPHA will be always less than and usually 00189 *> comparable with norm(A) in magnitude, and BETA always less 00190 *> than and usually comparable with norm(B). 00191 *> \endverbatim 00192 *> 00193 *> \param[out] VSL 00194 *> \verbatim 00195 *> VSL is COMPLEX array, dimension (LDVSL,N) 00196 *> If JOBVSL = 'V', VSL will contain the left Schur vectors. 00197 *> Not referenced if JOBVSL = 'N'. 00198 *> \endverbatim 00199 *> 00200 *> \param[in] LDVSL 00201 *> \verbatim 00202 *> LDVSL is INTEGER 00203 *> The leading dimension of the matrix VSL. LDVSL >=1, and 00204 *> if JOBVSL = 'V', LDVSL >= N. 00205 *> \endverbatim 00206 *> 00207 *> \param[out] VSR 00208 *> \verbatim 00209 *> VSR is COMPLEX array, dimension (LDVSR,N) 00210 *> If JOBVSR = 'V', VSR will contain the right Schur vectors. 00211 *> Not referenced if JOBVSR = 'N'. 00212 *> \endverbatim 00213 *> 00214 *> \param[in] LDVSR 00215 *> \verbatim 00216 *> LDVSR is INTEGER 00217 *> The leading dimension of the matrix VSR. LDVSR >= 1, and 00218 *> if JOBVSR = 'V', LDVSR >= N. 00219 *> \endverbatim 00220 *> 00221 *> \param[out] RCONDE 00222 *> \verbatim 00223 *> RCONDE is REAL array, dimension ( 2 ) 00224 *> If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the 00225 *> reciprocal condition numbers for the average of the selected 00226 *> eigenvalues. 00227 *> Not referenced if SENSE = 'N' or 'V'. 00228 *> \endverbatim 00229 *> 00230 *> \param[out] RCONDV 00231 *> \verbatim 00232 *> RCONDV is REAL array, dimension ( 2 ) 00233 *> If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the 00234 *> reciprocal condition number for the selected deflating 00235 *> subspaces. 00236 *> Not referenced if SENSE = 'N' or 'E'. 00237 *> \endverbatim 00238 *> 00239 *> \param[out] WORK 00240 *> \verbatim 00241 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00242 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00243 *> \endverbatim 00244 *> 00245 *> \param[in] LWORK 00246 *> \verbatim 00247 *> LWORK is INTEGER 00248 *> The dimension of the array WORK. 00249 *> If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B', 00250 *> LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else 00251 *> LWORK >= MAX(1,2*N). Note that 2*SDIM*(N-SDIM) <= N*N/2. 00252 *> Note also that an error is only returned if 00253 *> LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may 00254 *> not be large enough. 00255 *> 00256 *> If LWORK = -1, then a workspace query is assumed; the routine 00257 *> only calculates the bound on the optimal size of the WORK 00258 *> array and the minimum size of the IWORK array, returns these 00259 *> values as the first entries of the WORK and IWORK arrays, and 00260 *> no error message related to LWORK or LIWORK is issued by 00261 *> XERBLA. 00262 *> \endverbatim 00263 *> 00264 *> \param[out] RWORK 00265 *> \verbatim 00266 *> RWORK is REAL array, dimension ( 8*N ) 00267 *> Real workspace. 00268 *> \endverbatim 00269 *> 00270 *> \param[out] IWORK 00271 *> \verbatim 00272 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00273 *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 00274 *> \endverbatim 00275 *> 00276 *> \param[in] LIWORK 00277 *> \verbatim 00278 *> LIWORK is INTEGER 00279 *> The dimension of the array WORK. 00280 *> If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise 00281 *> LIWORK >= N+2. 00282 *> 00283 *> If LIWORK = -1, then a workspace query is assumed; the 00284 *> routine only calculates the bound on the optimal size of the 00285 *> WORK array and the minimum size of the IWORK array, returns 00286 *> these values as the first entries of the WORK and IWORK 00287 *> arrays, and no error message related to LWORK or LIWORK is 00288 *> issued by XERBLA. 00289 *> \endverbatim 00290 *> 00291 *> \param[out] BWORK 00292 *> \verbatim 00293 *> BWORK is LOGICAL array, dimension (N) 00294 *> Not referenced if SORT = 'N'. 00295 *> \endverbatim 00296 *> 00297 *> \param[out] INFO 00298 *> \verbatim 00299 *> INFO is INTEGER 00300 *> = 0: successful exit 00301 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00302 *> = 1,...,N: 00303 *> The QZ iteration failed. (A,B) are not in Schur 00304 *> form, but ALPHA(j) and BETA(j) should be correct for 00305 *> j=INFO+1,...,N. 00306 *> > N: =N+1: other than QZ iteration failed in CHGEQZ 00307 *> =N+2: after reordering, roundoff changed values of 00308 *> some complex eigenvalues so that leading 00309 *> eigenvalues in the Generalized Schur form no 00310 *> longer satisfy SELCTG=.TRUE. This could also 00311 *> be caused due to scaling. 00312 *> =N+3: reordering failed in CTGSEN. 00313 *> \endverbatim 00314 * 00315 * Authors: 00316 * ======== 00317 * 00318 *> \author Univ. of Tennessee 00319 *> \author Univ. of California Berkeley 00320 *> \author Univ. of Colorado Denver 00321 *> \author NAG Ltd. 00322 * 00323 *> \date November 2011 00324 * 00325 *> \ingroup complexGEeigen 00326 * 00327 * ===================================================================== 00328 SUBROUTINE CGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, 00329 $ B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, 00330 $ LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, 00331 $ IWORK, LIWORK, BWORK, INFO ) 00332 * 00333 * -- LAPACK driver routine (version 3.4.0) -- 00334 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00335 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00336 * November 2011 00337 * 00338 * .. Scalar Arguments .. 00339 CHARACTER JOBVSL, JOBVSR, SENSE, SORT 00340 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N, 00341 $ SDIM 00342 * .. 00343 * .. Array Arguments .. 00344 LOGICAL BWORK( * ) 00345 INTEGER IWORK( * ) 00346 REAL RCONDE( 2 ), RCONDV( 2 ), RWORK( * ) 00347 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00348 $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), 00349 $ WORK( * ) 00350 * .. 00351 * .. Function Arguments .. 00352 LOGICAL SELCTG 00353 EXTERNAL SELCTG 00354 * .. 00355 * 00356 * ===================================================================== 00357 * 00358 * .. Parameters .. 00359 REAL ZERO, ONE 00360 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00361 COMPLEX CZERO, CONE 00362 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00363 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00364 * .. 00365 * .. Local Scalars .. 00366 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL, 00367 $ LQUERY, WANTSB, WANTSE, WANTSN, WANTST, WANTSV 00368 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR, 00369 $ ILEFT, ILO, IRIGHT, IROWS, IRWRK, ITAU, IWRK, 00370 $ LIWMIN, LWRK, MAXWRK, MINWRK 00371 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL, 00372 $ PR, SMLNUM 00373 * .. 00374 * .. Local Arrays .. 00375 REAL DIF( 2 ) 00376 * .. 00377 * .. External Subroutines .. 00378 EXTERNAL CGEQRF, CGGBAK, CGGBAL, CGGHRD, CHGEQZ, CLACPY, 00379 $ CLASCL, CLASET, CTGSEN, CUNGQR, CUNMQR, SLABAD, 00380 $ XERBLA 00381 * .. 00382 * .. External Functions .. 00383 LOGICAL LSAME 00384 INTEGER ILAENV 00385 REAL CLANGE, SLAMCH 00386 EXTERNAL LSAME, ILAENV, CLANGE, SLAMCH 00387 * .. 00388 * .. Intrinsic Functions .. 00389 INTRINSIC MAX, SQRT 00390 * .. 00391 * .. Executable Statements .. 00392 * 00393 * Decode the input arguments 00394 * 00395 IF( LSAME( JOBVSL, 'N' ) ) THEN 00396 IJOBVL = 1 00397 ILVSL = .FALSE. 00398 ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN 00399 IJOBVL = 2 00400 ILVSL = .TRUE. 00401 ELSE 00402 IJOBVL = -1 00403 ILVSL = .FALSE. 00404 END IF 00405 * 00406 IF( LSAME( JOBVSR, 'N' ) ) THEN 00407 IJOBVR = 1 00408 ILVSR = .FALSE. 00409 ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN 00410 IJOBVR = 2 00411 ILVSR = .TRUE. 00412 ELSE 00413 IJOBVR = -1 00414 ILVSR = .FALSE. 00415 END IF 00416 * 00417 WANTST = LSAME( SORT, 'S' ) 00418 WANTSN = LSAME( SENSE, 'N' ) 00419 WANTSE = LSAME( SENSE, 'E' ) 00420 WANTSV = LSAME( SENSE, 'V' ) 00421 WANTSB = LSAME( SENSE, 'B' ) 00422 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00423 IF( WANTSN ) THEN 00424 IJOB = 0 00425 ELSE IF( WANTSE ) THEN 00426 IJOB = 1 00427 ELSE IF( WANTSV ) THEN 00428 IJOB = 2 00429 ELSE IF( WANTSB ) THEN 00430 IJOB = 4 00431 END IF 00432 * 00433 * Test the input arguments 00434 * 00435 INFO = 0 00436 IF( IJOBVL.LE.0 ) THEN 00437 INFO = -1 00438 ELSE IF( IJOBVR.LE.0 ) THEN 00439 INFO = -2 00440 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 00441 INFO = -3 00442 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 00443 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 00444 INFO = -5 00445 ELSE IF( N.LT.0 ) THEN 00446 INFO = -6 00447 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00448 INFO = -8 00449 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00450 INFO = -10 00451 ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN 00452 INFO = -15 00453 ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN 00454 INFO = -17 00455 END IF 00456 * 00457 * Compute workspace 00458 * (Note: Comments in the code beginning "Workspace:" describe the 00459 * minimal amount of workspace needed at that point in the code, 00460 * as well as the preferred amount for good performance. 00461 * NB refers to the optimal block size for the immediately 00462 * following subroutine, as returned by ILAENV.) 00463 * 00464 IF( INFO.EQ.0 ) THEN 00465 IF( N.GT.0) THEN 00466 MINWRK = 2*N 00467 MAXWRK = N*(1 + ILAENV( 1, 'CGEQRF', ' ', N, 1, N, 0 ) ) 00468 MAXWRK = MAX( MAXWRK, N*( 1 + 00469 $ ILAENV( 1, 'CUNMQR', ' ', N, 1, N, -1 ) ) ) 00470 IF( ILVSL ) THEN 00471 MAXWRK = MAX( MAXWRK, N*( 1 + 00472 $ ILAENV( 1, 'CUNGQR', ' ', N, 1, N, -1 ) ) ) 00473 END IF 00474 LWRK = MAXWRK 00475 IF( IJOB.GE.1 ) 00476 $ LWRK = MAX( LWRK, N*N/2 ) 00477 ELSE 00478 MINWRK = 1 00479 MAXWRK = 1 00480 LWRK = 1 00481 END IF 00482 WORK( 1 ) = LWRK 00483 IF( WANTSN .OR. N.EQ.0 ) THEN 00484 LIWMIN = 1 00485 ELSE 00486 LIWMIN = N + 2 00487 END IF 00488 IWORK( 1 ) = LIWMIN 00489 * 00490 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00491 INFO = -21 00492 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY) THEN 00493 INFO = -24 00494 END IF 00495 END IF 00496 * 00497 IF( INFO.NE.0 ) THEN 00498 CALL XERBLA( 'CGGESX', -INFO ) 00499 RETURN 00500 ELSE IF (LQUERY) THEN 00501 RETURN 00502 END IF 00503 * 00504 * Quick return if possible 00505 * 00506 IF( N.EQ.0 ) THEN 00507 SDIM = 0 00508 RETURN 00509 END IF 00510 * 00511 * Get machine constants 00512 * 00513 EPS = SLAMCH( 'P' ) 00514 SMLNUM = SLAMCH( 'S' ) 00515 BIGNUM = ONE / SMLNUM 00516 CALL SLABAD( SMLNUM, BIGNUM ) 00517 SMLNUM = SQRT( SMLNUM ) / EPS 00518 BIGNUM = ONE / SMLNUM 00519 * 00520 * Scale A if max element outside range [SMLNUM,BIGNUM] 00521 * 00522 ANRM = CLANGE( 'M', N, N, A, LDA, RWORK ) 00523 ILASCL = .FALSE. 00524 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00525 ANRMTO = SMLNUM 00526 ILASCL = .TRUE. 00527 ELSE IF( ANRM.GT.BIGNUM ) THEN 00528 ANRMTO = BIGNUM 00529 ILASCL = .TRUE. 00530 END IF 00531 IF( ILASCL ) 00532 $ CALL CLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR ) 00533 * 00534 * Scale B if max element outside range [SMLNUM,BIGNUM] 00535 * 00536 BNRM = CLANGE( 'M', N, N, B, LDB, RWORK ) 00537 ILBSCL = .FALSE. 00538 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00539 BNRMTO = SMLNUM 00540 ILBSCL = .TRUE. 00541 ELSE IF( BNRM.GT.BIGNUM ) THEN 00542 BNRMTO = BIGNUM 00543 ILBSCL = .TRUE. 00544 END IF 00545 IF( ILBSCL ) 00546 $ CALL CLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR ) 00547 * 00548 * Permute the matrix to make it more nearly triangular 00549 * (Real Workspace: need 6*N) 00550 * 00551 ILEFT = 1 00552 IRIGHT = N + 1 00553 IRWRK = IRIGHT + N 00554 CALL CGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ), 00555 $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR ) 00556 * 00557 * Reduce B to triangular form (QR decomposition of B) 00558 * (Complex Workspace: need N, prefer N*NB) 00559 * 00560 IROWS = IHI + 1 - ILO 00561 ICOLS = N + 1 - ILO 00562 ITAU = 1 00563 IWRK = ITAU + IROWS 00564 CALL CGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 00565 $ WORK( IWRK ), LWORK+1-IWRK, IERR ) 00566 * 00567 * Apply the unitary transformation to matrix A 00568 * (Complex Workspace: need N, prefer N*NB) 00569 * 00570 CALL CUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 00571 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ), 00572 $ LWORK+1-IWRK, IERR ) 00573 * 00574 * Initialize VSL 00575 * (Complex Workspace: need N, prefer N*NB) 00576 * 00577 IF( ILVSL ) THEN 00578 CALL CLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL ) 00579 IF( IROWS.GT.1 ) THEN 00580 CALL CLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 00581 $ VSL( ILO+1, ILO ), LDVSL ) 00582 END IF 00583 CALL CUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL, 00584 $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR ) 00585 END IF 00586 * 00587 * Initialize VSR 00588 * 00589 IF( ILVSR ) 00590 $ CALL CLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR ) 00591 * 00592 * Reduce to generalized Hessenberg form 00593 * (Workspace: none needed) 00594 * 00595 CALL CGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL, 00596 $ LDVSL, VSR, LDVSR, IERR ) 00597 * 00598 SDIM = 0 00599 * 00600 * Perform QZ algorithm, computing Schur vectors if desired 00601 * (Complex Workspace: need N) 00602 * (Real Workspace: need N) 00603 * 00604 IWRK = ITAU 00605 CALL CHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, 00606 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWRK ), 00607 $ LWORK+1-IWRK, RWORK( IRWRK ), IERR ) 00608 IF( IERR.NE.0 ) THEN 00609 IF( IERR.GT.0 .AND. IERR.LE.N ) THEN 00610 INFO = IERR 00611 ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN 00612 INFO = IERR - N 00613 ELSE 00614 INFO = N + 1 00615 END IF 00616 GO TO 40 00617 END IF 00618 * 00619 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of 00620 * condition number(s) 00621 * 00622 IF( WANTST ) THEN 00623 * 00624 * Undo scaling on eigenvalues before SELCTGing 00625 * 00626 IF( ILASCL ) 00627 $ CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 00628 IF( ILBSCL ) 00629 $ CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 00630 * 00631 * Select eigenvalues 00632 * 00633 DO 10 I = 1, N 00634 BWORK( I ) = SELCTG( ALPHA( I ), BETA( I ) ) 00635 10 CONTINUE 00636 * 00637 * Reorder eigenvalues, transform Generalized Schur vectors, and 00638 * compute reciprocal condition numbers 00639 * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM)) 00640 * otherwise, need 1 ) 00641 * 00642 CALL CTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, 00643 $ ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PL, PR, 00644 $ DIF, WORK( IWRK ), LWORK-IWRK+1, IWORK, LIWORK, 00645 $ IERR ) 00646 * 00647 IF( IJOB.GE.1 ) 00648 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 00649 IF( IERR.EQ.-21 ) THEN 00650 * 00651 * not enough complex workspace 00652 * 00653 INFO = -21 00654 ELSE 00655 IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN 00656 RCONDE( 1 ) = PL 00657 RCONDE( 2 ) = PR 00658 END IF 00659 IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN 00660 RCONDV( 1 ) = DIF( 1 ) 00661 RCONDV( 2 ) = DIF( 2 ) 00662 END IF 00663 IF( IERR.EQ.1 ) 00664 $ INFO = N + 3 00665 END IF 00666 * 00667 END IF 00668 * 00669 * Apply permutation to VSL and VSR 00670 * (Workspace: none needed) 00671 * 00672 IF( ILVSL ) 00673 $ CALL CGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ), 00674 $ RWORK( IRIGHT ), N, VSL, LDVSL, IERR ) 00675 * 00676 IF( ILVSR ) 00677 $ CALL CGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ), 00678 $ RWORK( IRIGHT ), N, VSR, LDVSR, IERR ) 00679 * 00680 * Undo scaling 00681 * 00682 IF( ILASCL ) THEN 00683 CALL CLASCL( 'U', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR ) 00684 CALL CLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR ) 00685 END IF 00686 * 00687 IF( ILBSCL ) THEN 00688 CALL CLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR ) 00689 CALL CLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR ) 00690 END IF 00691 * 00692 IF( WANTST ) THEN 00693 * 00694 * Check if reordering is correct 00695 * 00696 LASTSL = .TRUE. 00697 SDIM = 0 00698 DO 30 I = 1, N 00699 CURSL = SELCTG( ALPHA( I ), BETA( I ) ) 00700 IF( CURSL ) 00701 $ SDIM = SDIM + 1 00702 IF( CURSL .AND. .NOT.LASTSL ) 00703 $ INFO = N + 2 00704 LASTSL = CURSL 00705 30 CONTINUE 00706 * 00707 END IF 00708 * 00709 40 CONTINUE 00710 * 00711 WORK( 1 ) = MAXWRK 00712 IWORK( 1 ) = LIWMIN 00713 * 00714 RETURN 00715 * 00716 * End of CGGESX 00717 * 00718 END