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