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