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