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