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