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