![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZGEESX 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 ZGEESX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeesx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeesx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeesx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, 00022 * VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, 00023 * BWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBVS, SENSE, SORT 00027 * INTEGER INFO, LDA, LDVS, LWORK, N, SDIM 00028 * DOUBLE PRECISION RCONDE, RCONDV 00029 * .. 00030 * .. Array Arguments .. 00031 * LOGICAL BWORK( * ) 00032 * DOUBLE PRECISION RWORK( * ) 00033 * COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) 00034 * .. 00035 * .. Function Arguments .. 00036 * LOGICAL SELECT 00037 * EXTERNAL SELECT 00038 * .. 00039 * 00040 * 00041 *> \par Purpose: 00042 * ============= 00043 *> 00044 *> \verbatim 00045 *> 00046 *> ZGEESX computes for an N-by-N complex nonsymmetric matrix A, the 00047 *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur 00048 *> vectors Z. This gives the Schur factorization A = Z*T*(Z**H). 00049 *> 00050 *> Optionally, it also orders the eigenvalues on the diagonal of the 00051 *> Schur form so that selected eigenvalues are at the top left; 00052 *> computes a reciprocal condition number for the average of the 00053 *> selected eigenvalues (RCONDE); and computes a reciprocal condition 00054 *> number for the right invariant subspace corresponding to the 00055 *> selected eigenvalues (RCONDV). The leading columns of Z form an 00056 *> orthonormal basis for this invariant subspace. 00057 *> 00058 *> For further explanation of the reciprocal condition numbers RCONDE 00059 *> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where 00060 *> these quantities are called s and sep respectively). 00061 *> 00062 *> A complex matrix is in Schur form if it is upper triangular. 00063 *> \endverbatim 00064 * 00065 * Arguments: 00066 * ========== 00067 * 00068 *> \param[in] JOBVS 00069 *> \verbatim 00070 *> JOBVS is CHARACTER*1 00071 *> = 'N': Schur vectors are not computed; 00072 *> = 'V': Schur vectors are computed. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] SORT 00076 *> \verbatim 00077 *> SORT is CHARACTER*1 00078 *> Specifies whether or not to order the eigenvalues on the 00079 *> diagonal of the Schur form. 00080 *> = 'N': Eigenvalues are not ordered; 00081 *> = 'S': Eigenvalues are ordered (see SELECT). 00082 *> \endverbatim 00083 *> 00084 *> \param[in] SELECT 00085 *> \verbatim 00086 *> SELECT is procedure) LOGICAL FUNCTION of one COMPLEX*16 argument 00087 *> SELECT must be declared EXTERNAL in the calling subroutine. 00088 *> If SORT = 'S', SELECT is used to select eigenvalues to order 00089 *> to the top left of the Schur form. 00090 *> If SORT = 'N', SELECT is not referenced. 00091 *> An eigenvalue W(j) is selected if SELECT(W(j)) is true. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] SENSE 00095 *> \verbatim 00096 *> SENSE is CHARACTER*1 00097 *> Determines which reciprocal condition numbers are computed. 00098 *> = 'N': None are computed; 00099 *> = 'E': Computed for average of selected eigenvalues only; 00100 *> = 'V': Computed for selected right invariant subspace only; 00101 *> = 'B': Computed for both. 00102 *> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] N 00106 *> \verbatim 00107 *> N is INTEGER 00108 *> The order of the matrix A. N >= 0. 00109 *> \endverbatim 00110 *> 00111 *> \param[in,out] A 00112 *> \verbatim 00113 *> A is COMPLEX*16 array, dimension (LDA, N) 00114 *> On entry, the N-by-N matrix A. 00115 *> On exit, A is overwritten by its Schur form T. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] LDA 00119 *> \verbatim 00120 *> LDA is INTEGER 00121 *> The leading dimension of the array A. LDA >= max(1,N). 00122 *> \endverbatim 00123 *> 00124 *> \param[out] SDIM 00125 *> \verbatim 00126 *> SDIM is INTEGER 00127 *> If SORT = 'N', SDIM = 0. 00128 *> If SORT = 'S', SDIM = number of eigenvalues for which 00129 *> SELECT is true. 00130 *> \endverbatim 00131 *> 00132 *> \param[out] W 00133 *> \verbatim 00134 *> W is COMPLEX*16 array, dimension (N) 00135 *> W contains the computed eigenvalues, in the same order 00136 *> that they appear on the diagonal of the output Schur form T. 00137 *> \endverbatim 00138 *> 00139 *> \param[out] VS 00140 *> \verbatim 00141 *> VS is COMPLEX*16 array, dimension (LDVS,N) 00142 *> If JOBVS = 'V', VS contains the unitary matrix Z of Schur 00143 *> vectors. 00144 *> If JOBVS = 'N', VS is not referenced. 00145 *> \endverbatim 00146 *> 00147 *> \param[in] LDVS 00148 *> \verbatim 00149 *> LDVS is INTEGER 00150 *> The leading dimension of the array VS. LDVS >= 1, and if 00151 *> JOBVS = 'V', LDVS >= N. 00152 *> \endverbatim 00153 *> 00154 *> \param[out] RCONDE 00155 *> \verbatim 00156 *> RCONDE is DOUBLE PRECISION 00157 *> If SENSE = 'E' or 'B', RCONDE contains the reciprocal 00158 *> condition number for the average of the selected eigenvalues. 00159 *> Not referenced if SENSE = 'N' or 'V'. 00160 *> \endverbatim 00161 *> 00162 *> \param[out] RCONDV 00163 *> \verbatim 00164 *> RCONDV is DOUBLE PRECISION 00165 *> If SENSE = 'V' or 'B', RCONDV contains the reciprocal 00166 *> condition number for the selected right invariant subspace. 00167 *> Not referenced if SENSE = 'N' or 'E'. 00168 *> \endverbatim 00169 *> 00170 *> \param[out] WORK 00171 *> \verbatim 00172 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 00173 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00174 *> \endverbatim 00175 *> 00176 *> \param[in] LWORK 00177 *> \verbatim 00178 *> LWORK is INTEGER 00179 *> The dimension of the array WORK. LWORK >= max(1,2*N). 00180 *> Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM), 00181 *> where SDIM is the number of selected eigenvalues computed by 00182 *> this routine. Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also 00183 *> that an error is only returned if LWORK < max(1,2*N), but if 00184 *> SENSE = 'E' or 'V' or 'B' this may not be large enough. 00185 *> For good performance, LWORK must generally be larger. 00186 *> 00187 *> If LWORK = -1, then a workspace query is assumed; the routine 00188 *> only calculates upper bound on the optimal size of the 00189 *> array WORK, returns this value as the first entry of the WORK 00190 *> array, and no error message related to LWORK is issued by 00191 *> XERBLA. 00192 *> \endverbatim 00193 *> 00194 *> \param[out] RWORK 00195 *> \verbatim 00196 *> RWORK is DOUBLE PRECISION array, dimension (N) 00197 *> \endverbatim 00198 *> 00199 *> \param[out] BWORK 00200 *> \verbatim 00201 *> BWORK is LOGICAL array, dimension (N) 00202 *> Not referenced if SORT = 'N'. 00203 *> \endverbatim 00204 *> 00205 *> \param[out] INFO 00206 *> \verbatim 00207 *> INFO is INTEGER 00208 *> = 0: successful exit 00209 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00210 *> > 0: if INFO = i, and i is 00211 *> <= N: the QR algorithm failed to compute all the 00212 *> eigenvalues; elements 1:ILO-1 and i+1:N of W 00213 *> contain those eigenvalues which have converged; if 00214 *> JOBVS = 'V', VS contains the transformation which 00215 *> reduces A to its partially converged Schur form. 00216 *> = N+1: the eigenvalues could not be reordered because some 00217 *> eigenvalues were too close to separate (the problem 00218 *> is very ill-conditioned); 00219 *> = N+2: after reordering, roundoff changed values of some 00220 *> complex eigenvalues so that leading eigenvalues in 00221 *> the Schur form no longer satisfy SELECT=.TRUE. This 00222 *> could also be caused by underflow due to scaling. 00223 *> \endverbatim 00224 * 00225 * Authors: 00226 * ======== 00227 * 00228 *> \author Univ. of Tennessee 00229 *> \author Univ. of California Berkeley 00230 *> \author Univ. of Colorado Denver 00231 *> \author NAG Ltd. 00232 * 00233 *> \date November 2011 00234 * 00235 *> \ingroup complex16GEeigen 00236 * 00237 * ===================================================================== 00238 SUBROUTINE ZGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, 00239 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, 00240 $ BWORK, INFO ) 00241 * 00242 * -- LAPACK driver routine (version 3.4.0) -- 00243 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00244 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00245 * November 2011 00246 * 00247 * .. Scalar Arguments .. 00248 CHARACTER JOBVS, SENSE, SORT 00249 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM 00250 DOUBLE PRECISION RCONDE, RCONDV 00251 * .. 00252 * .. Array Arguments .. 00253 LOGICAL BWORK( * ) 00254 DOUBLE PRECISION RWORK( * ) 00255 COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) 00256 * .. 00257 * .. Function Arguments .. 00258 LOGICAL SELECT 00259 EXTERNAL SELECT 00260 * .. 00261 * 00262 * ===================================================================== 00263 * 00264 * .. Parameters .. 00265 DOUBLE PRECISION ZERO, ONE 00266 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00267 * .. 00268 * .. Local Scalars .. 00269 LOGICAL LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST, 00270 $ WANTSV, WANTVS 00271 INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO, 00272 $ ITAU, IWRK, LWRK, MAXWRK, MINWRK 00273 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM 00274 * .. 00275 * .. Local Arrays .. 00276 DOUBLE PRECISION DUM( 1 ) 00277 * .. 00278 * .. External Subroutines .. 00279 EXTERNAL DLABAD, DLASCL, XERBLA, ZCOPY, ZGEBAK, ZGEBAL, 00280 $ ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR 00281 * .. 00282 * .. External Functions .. 00283 LOGICAL LSAME 00284 INTEGER ILAENV 00285 DOUBLE PRECISION DLAMCH, ZLANGE 00286 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE 00287 * .. 00288 * .. Intrinsic Functions .. 00289 INTRINSIC MAX, SQRT 00290 * .. 00291 * .. Executable Statements .. 00292 * 00293 * Test the input arguments 00294 * 00295 INFO = 0 00296 WANTVS = LSAME( JOBVS, 'V' ) 00297 WANTST = LSAME( SORT, 'S' ) 00298 WANTSN = LSAME( SENSE, 'N' ) 00299 WANTSE = LSAME( SENSE, 'E' ) 00300 WANTSV = LSAME( SENSE, 'V' ) 00301 WANTSB = LSAME( SENSE, 'B' ) 00302 LQUERY = ( LWORK.EQ.-1 ) 00303 * 00304 IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN 00305 INFO = -1 00306 ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN 00307 INFO = -2 00308 ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR. 00309 $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN 00310 INFO = -4 00311 ELSE IF( N.LT.0 ) THEN 00312 INFO = -5 00313 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00314 INFO = -7 00315 ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN 00316 INFO = -11 00317 END IF 00318 * 00319 * Compute workspace 00320 * (Note: Comments in the code beginning "Workspace:" describe the 00321 * minimal amount of real workspace needed at that point in the 00322 * code, as well as the preferred amount for good performance. 00323 * CWorkspace refers to complex workspace, and RWorkspace to real 00324 * workspace. NB refers to the optimal block size for the 00325 * immediately following subroutine, as returned by ILAENV. 00326 * HSWORK refers to the workspace preferred by ZHSEQR, as 00327 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 00328 * the worst case. 00329 * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed 00330 * depends on SDIM, which is computed by the routine ZTRSEN later 00331 * in the code.) 00332 * 00333 IF( INFO.EQ.0 ) THEN 00334 IF( N.EQ.0 ) THEN 00335 MINWRK = 1 00336 LWRK = 1 00337 ELSE 00338 MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 ) 00339 MINWRK = 2*N 00340 * 00341 CALL ZHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS, 00342 $ WORK, -1, IEVAL ) 00343 HSWORK = WORK( 1 ) 00344 * 00345 IF( .NOT.WANTVS ) THEN 00346 MAXWRK = MAX( MAXWRK, HSWORK ) 00347 ELSE 00348 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR', 00349 $ ' ', N, 1, N, -1 ) ) 00350 MAXWRK = MAX( MAXWRK, HSWORK ) 00351 END IF 00352 LWRK = MAXWRK 00353 IF( .NOT.WANTSN ) 00354 $ LWRK = MAX( LWRK, ( N*N )/2 ) 00355 END IF 00356 WORK( 1 ) = LWRK 00357 * 00358 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00359 INFO = -15 00360 END IF 00361 END IF 00362 * 00363 IF( INFO.NE.0 ) THEN 00364 CALL XERBLA( 'ZGEESX', -INFO ) 00365 RETURN 00366 ELSE IF( LQUERY ) THEN 00367 RETURN 00368 END IF 00369 * 00370 * Quick return if possible 00371 * 00372 IF( N.EQ.0 ) THEN 00373 SDIM = 0 00374 RETURN 00375 END IF 00376 * 00377 * Get machine constants 00378 * 00379 EPS = DLAMCH( 'P' ) 00380 SMLNUM = DLAMCH( 'S' ) 00381 BIGNUM = ONE / SMLNUM 00382 CALL DLABAD( SMLNUM, BIGNUM ) 00383 SMLNUM = SQRT( SMLNUM ) / EPS 00384 BIGNUM = ONE / SMLNUM 00385 * 00386 * Scale A if max element outside range [SMLNUM,BIGNUM] 00387 * 00388 ANRM = ZLANGE( 'M', N, N, A, LDA, DUM ) 00389 SCALEA = .FALSE. 00390 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00391 SCALEA = .TRUE. 00392 CSCALE = SMLNUM 00393 ELSE IF( ANRM.GT.BIGNUM ) THEN 00394 SCALEA = .TRUE. 00395 CSCALE = BIGNUM 00396 END IF 00397 IF( SCALEA ) 00398 $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 00399 * 00400 * 00401 * Permute the matrix to make it more nearly triangular 00402 * (CWorkspace: none) 00403 * (RWorkspace: need N) 00404 * 00405 IBAL = 1 00406 CALL ZGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR ) 00407 * 00408 * Reduce to upper Hessenberg form 00409 * (CWorkspace: need 2*N, prefer N+N*NB) 00410 * (RWorkspace: none) 00411 * 00412 ITAU = 1 00413 IWRK = N + ITAU 00414 CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 00415 $ LWORK-IWRK+1, IERR ) 00416 * 00417 IF( WANTVS ) THEN 00418 * 00419 * Copy Householder vectors to VS 00420 * 00421 CALL ZLACPY( 'L', N, N, A, LDA, VS, LDVS ) 00422 * 00423 * Generate unitary matrix in VS 00424 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) 00425 * (RWorkspace: none) 00426 * 00427 CALL ZUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ), 00428 $ LWORK-IWRK+1, IERR ) 00429 END IF 00430 * 00431 SDIM = 0 00432 * 00433 * Perform QR iteration, accumulating Schur vectors in VS if desired 00434 * (CWorkspace: need 1, prefer HSWORK (see comments) ) 00435 * (RWorkspace: none) 00436 * 00437 IWRK = ITAU 00438 CALL ZHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS, 00439 $ WORK( IWRK ), LWORK-IWRK+1, IEVAL ) 00440 IF( IEVAL.GT.0 ) 00441 $ INFO = IEVAL 00442 * 00443 * Sort eigenvalues if desired 00444 * 00445 IF( WANTST .AND. INFO.EQ.0 ) THEN 00446 IF( SCALEA ) 00447 $ CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR ) 00448 DO 10 I = 1, N 00449 BWORK( I ) = SELECT( W( I ) ) 00450 10 CONTINUE 00451 * 00452 * Reorder eigenvalues, transform Schur vectors, and compute 00453 * reciprocal condition numbers 00454 * (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM) 00455 * otherwise, need none ) 00456 * (RWorkspace: none) 00457 * 00458 CALL ZTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM, 00459 $ RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1, 00460 $ ICOND ) 00461 IF( .NOT.WANTSN ) 00462 $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) ) 00463 IF( ICOND.EQ.-14 ) THEN 00464 * 00465 * Not enough complex workspace 00466 * 00467 INFO = -15 00468 END IF 00469 END IF 00470 * 00471 IF( WANTVS ) THEN 00472 * 00473 * Undo balancing 00474 * (CWorkspace: none) 00475 * (RWorkspace: need N) 00476 * 00477 CALL ZGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS, 00478 $ IERR ) 00479 END IF 00480 * 00481 IF( SCALEA ) THEN 00482 * 00483 * Undo scaling for the Schur form of A 00484 * 00485 CALL ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR ) 00486 CALL ZCOPY( N, A, LDA+1, W, 1 ) 00487 IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN 00488 DUM( 1 ) = RCONDV 00489 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR ) 00490 RCONDV = DUM( 1 ) 00491 END IF 00492 END IF 00493 * 00494 WORK( 1 ) = MAXWRK 00495 RETURN 00496 * 00497 * End of ZGEESX 00498 * 00499 END