LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgees.f
Go to the documentation of this file.
1*> \brief <b> DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGEES + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgees.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgees.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgees.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
22* VS, LDVS, WORK, LWORK, BWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBVS, SORT
26* INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
27* ..
28* .. Array Arguments ..
29* LOGICAL BWORK( * )
30* DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
31* $ WR( * )
32* ..
33* .. Function Arguments ..
34* LOGICAL SELECT
35* EXTERNAL SELECT
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> DGEES computes for an N-by-N real nonsymmetric matrix A, the
45*> eigenvalues, the real Schur form T, and, optionally, the matrix of
46*> Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
47*>
48*> Optionally, it also orders the eigenvalues on the diagonal of the
49*> real Schur form so that selected eigenvalues are at the top left.
50*> The leading columns of Z then form an orthonormal basis for the
51*> invariant subspace corresponding to the selected eigenvalues.
52*>
53*> A matrix is in real Schur form if it is upper quasi-triangular with
54*> 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
55*> form
56*> [ a b ]
57*> [ c a ]
58*>
59*> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
60*> \endverbatim
61*
62* Arguments:
63* ==========
64*
65*> \param[in] JOBVS
66*> \verbatim
67*> JOBVS is CHARACTER*1
68*> = 'N': Schur vectors are not computed;
69*> = 'V': Schur vectors are computed.
70*> \endverbatim
71*>
72*> \param[in] SORT
73*> \verbatim
74*> SORT is CHARACTER*1
75*> Specifies whether or not to order the eigenvalues on the
76*> diagonal of the Schur form.
77*> = 'N': Eigenvalues are not ordered;
78*> = 'S': Eigenvalues are ordered (see SELECT).
79*> \endverbatim
80*>
81*> \param[in] SELECT
82*> \verbatim
83*> SELECT is a LOGICAL FUNCTION of two DOUBLE PRECISION arguments
84*> SELECT must be declared EXTERNAL in the calling subroutine.
85*> If SORT = 'S', SELECT is used to select eigenvalues to sort
86*> to the top left of the Schur form.
87*> If SORT = 'N', SELECT is not referenced.
88*> An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
89*> SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
90*> conjugate pair of eigenvalues is selected, then both complex
91*> eigenvalues are selected.
92*> Note that a selected complex eigenvalue may no longer
93*> satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
94*> ordering may change the value of complex eigenvalues
95*> (especially if the eigenvalue is ill-conditioned); in this
96*> case INFO is set to N+2 (see INFO below).
97*> \endverbatim
98*>
99*> \param[in] N
100*> \verbatim
101*> N is INTEGER
102*> The order of the matrix A. N >= 0.
103*> \endverbatim
104*>
105*> \param[in,out] A
106*> \verbatim
107*> A is DOUBLE PRECISION array, dimension (LDA,N)
108*> On entry, the N-by-N matrix A.
109*> On exit, A has been overwritten by its real Schur form T.
110*> \endverbatim
111*>
112*> \param[in] LDA
113*> \verbatim
114*> LDA is INTEGER
115*> The leading dimension of the array A. LDA >= max(1,N).
116*> \endverbatim
117*>
118*> \param[out] SDIM
119*> \verbatim
120*> SDIM is INTEGER
121*> If SORT = 'N', SDIM = 0.
122*> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
123*> for which SELECT is true. (Complex conjugate
124*> pairs for which SELECT is true for either
125*> eigenvalue count as 2.)
126*> \endverbatim
127*>
128*> \param[out] WR
129*> \verbatim
130*> WR is DOUBLE PRECISION array, dimension (N)
131*> \endverbatim
132*>
133*> \param[out] WI
134*> \verbatim
135*> WI is DOUBLE PRECISION array, dimension (N)
136*> WR and WI contain the real and imaginary parts,
137*> respectively, of the computed eigenvalues in the same order
138*> that they appear on the diagonal of the output Schur form T.
139*> Complex conjugate pairs of eigenvalues will appear
140*> consecutively with the eigenvalue having the positive
141*> imaginary part first.
142*> \endverbatim
143*>
144*> \param[out] VS
145*> \verbatim
146*> VS is DOUBLE PRECISION array, dimension (LDVS,N)
147*> If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
148*> vectors.
149*> If JOBVS = 'N', VS is not referenced.
150*> \endverbatim
151*>
152*> \param[in] LDVS
153*> \verbatim
154*> LDVS is INTEGER
155*> The leading dimension of the array VS. LDVS >= 1; if
156*> JOBVS = 'V', LDVS >= N.
157*> \endverbatim
158*>
159*> \param[out] WORK
160*> \verbatim
161*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
162*> On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
163*> \endverbatim
164*>
165*> \param[in] LWORK
166*> \verbatim
167*> LWORK is INTEGER
168*> The dimension of the array WORK. LWORK >= max(1,3*N).
169*> For good performance, LWORK must generally be larger.
170*>
171*> If LWORK = -1, then a workspace query is assumed; the routine
172*> only calculates the optimal size of the WORK array, returns
173*> this value as the first entry of the WORK array, and no error
174*> message related to LWORK is issued by XERBLA.
175*> \endverbatim
176*>
177*> \param[out] BWORK
178*> \verbatim
179*> BWORK is LOGICAL array, dimension (N)
180*> Not referenced if SORT = 'N'.
181*> \endverbatim
182*>
183*> \param[out] INFO
184*> \verbatim
185*> INFO is INTEGER
186*> = 0: successful exit
187*> < 0: if INFO = -i, the i-th argument had an illegal value.
188*> > 0: if INFO = i, and i is
189*> <= N: the QR algorithm failed to compute all the
190*> eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
191*> contain those eigenvalues which have converged; if
192*> JOBVS = 'V', VS contains the matrix which reduces A
193*> to its partially converged Schur form.
194*> = N+1: the eigenvalues could not be reordered because some
195*> eigenvalues were too close to separate (the problem
196*> is very ill-conditioned);
197*> = N+2: after reordering, roundoff changed values of some
198*> complex eigenvalues so that leading eigenvalues in
199*> the Schur form no longer satisfy SELECT=.TRUE. This
200*> could also be caused by underflow due to scaling.
201*> \endverbatim
202*
203* Authors:
204* ========
205*
206*> \author Univ. of Tennessee
207*> \author Univ. of California Berkeley
208*> \author Univ. of Colorado Denver
209*> \author NAG Ltd.
210*
211*> \ingroup doubleGEeigen
212*
213* =====================================================================
214 SUBROUTINE dgees( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
215 $ VS, LDVS, WORK, LWORK, BWORK, INFO )
216*
217* -- LAPACK driver routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 CHARACTER JOBVS, SORT
223 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
224* ..
225* .. Array Arguments ..
226 LOGICAL BWORK( * )
227 DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
228 $ wr( * )
229* ..
230* .. Function Arguments ..
231 LOGICAL SELECT
232 EXTERNAL SELECT
233* ..
234*
235* =====================================================================
236*
237* .. Parameters ..
238 DOUBLE PRECISION ZERO, ONE
239 parameter( zero = 0.0d0, one = 1.0d0 )
240* ..
241* .. Local Scalars ..
242 LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
243 $ wantvs
244 INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
245 $ ihi, ilo, inxt, ip, itau, iwrk, maxwrk, minwrk
246 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
247* ..
248* .. Local Arrays ..
249 INTEGER IDUM( 1 )
250 DOUBLE PRECISION DUM( 1 )
251* ..
252* .. External Subroutines ..
253 EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
255* ..
256* .. External Functions ..
257 LOGICAL LSAME
258 INTEGER ILAENV
259 DOUBLE PRECISION DLAMCH, DLANGE
260 EXTERNAL lsame, ilaenv, dlamch, dlange
261* ..
262* .. Intrinsic Functions ..
263 INTRINSIC max, sqrt
264* ..
265* .. Executable Statements ..
266*
267* Test the input arguments
268*
269 info = 0
270 lquery = ( lwork.EQ.-1 )
271 wantvs = lsame( jobvs, 'V' )
272 wantst = lsame( sort, 'S' )
273 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
274 info = -1
275 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
276 info = -2
277 ELSE IF( n.LT.0 ) THEN
278 info = -4
279 ELSE IF( lda.LT.max( 1, n ) ) THEN
280 info = -6
281 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
282 info = -11
283 END IF
284*
285* Compute workspace
286* (Note: Comments in the code beginning "Workspace:" describe the
287* minimal amount of workspace needed at that point in the code,
288* as well as the preferred amount for good performance.
289* NB refers to the optimal block size for the immediately
290* following subroutine, as returned by ILAENV.
291* HSWORK refers to the workspace preferred by DHSEQR, as
292* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
293* the worst case.)
294*
295 IF( info.EQ.0 ) THEN
296 IF( n.EQ.0 ) THEN
297 minwrk = 1
298 maxwrk = 1
299 ELSE
300 maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
301 minwrk = 3*n
302*
303 CALL dhseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
304 $ work, -1, ieval )
305 hswork = int( work( 1 ) )
306*
307 IF( .NOT.wantvs ) THEN
308 maxwrk = max( maxwrk, n + hswork )
309 ELSE
310 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
311 $ 'DORGHR', ' ', n, 1, n, -1 ) )
312 maxwrk = max( maxwrk, n + hswork )
313 END IF
314 END IF
315 work( 1 ) = maxwrk
316*
317 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
318 info = -13
319 END IF
320 END IF
321*
322 IF( info.NE.0 ) THEN
323 CALL xerbla( 'DGEES ', -info )
324 RETURN
325 ELSE IF( lquery ) THEN
326 RETURN
327 END IF
328*
329* Quick return if possible
330*
331 IF( n.EQ.0 ) THEN
332 sdim = 0
333 RETURN
334 END IF
335*
336* Get machine constants
337*
338 eps = dlamch( 'P' )
339 smlnum = dlamch( 'S' )
340 bignum = one / smlnum
341 CALL dlabad( smlnum, bignum )
342 smlnum = sqrt( smlnum ) / eps
343 bignum = one / smlnum
344*
345* Scale A if max element outside range [SMLNUM,BIGNUM]
346*
347 anrm = dlange( 'M', n, n, a, lda, dum )
348 scalea = .false.
349 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
350 scalea = .true.
351 cscale = smlnum
352 ELSE IF( anrm.GT.bignum ) THEN
353 scalea = .true.
354 cscale = bignum
355 END IF
356 IF( scalea )
357 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
358*
359* Permute the matrix to make it more nearly triangular
360* (Workspace: need N)
361*
362 ibal = 1
363 CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
364*
365* Reduce to upper Hessenberg form
366* (Workspace: need 3*N, prefer 2*N+N*NB)
367*
368 itau = n + ibal
369 iwrk = n + itau
370 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
371 $ lwork-iwrk+1, ierr )
372*
373 IF( wantvs ) THEN
374*
375* Copy Householder vectors to VS
376*
377 CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
378*
379* Generate orthogonal matrix in VS
380* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
381*
382 CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
383 $ lwork-iwrk+1, ierr )
384 END IF
385*
386 sdim = 0
387*
388* Perform QR iteration, accumulating Schur vectors in VS if desired
389* (Workspace: need N+1, prefer N+HSWORK (see comments) )
390*
391 iwrk = itau
392 CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
393 $ work( iwrk ), lwork-iwrk+1, ieval )
394 IF( ieval.GT.0 )
395 $ info = ieval
396*
397* Sort eigenvalues if desired
398*
399 IF( wantst .AND. info.EQ.0 ) THEN
400 IF( scalea ) THEN
401 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
402 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
403 END IF
404 DO 10 i = 1, n
405 bwork( i ) = SELECT( wr( i ), wi( i ) )
406 10 CONTINUE
407*
408* Reorder eigenvalues and transform Schur vectors
409* (Workspace: none needed)
410*
411 CALL dtrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
412 $ sdim, s, sep, work( iwrk ), lwork-iwrk+1, idum, 1,
413 $ icond )
414 IF( icond.GT.0 )
415 $ info = n + icond
416 END IF
417*
418 IF( wantvs ) THEN
419*
420* Undo balancing
421* (Workspace: need N)
422*
423 CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
424 $ ierr )
425 END IF
426*
427 IF( scalea ) THEN
428*
429* Undo scaling for the Schur form of A
430*
431 CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
432 CALL dcopy( n, a, lda+1, wr, 1 )
433 IF( cscale.EQ.smlnum ) THEN
434*
435* If scaling back towards underflow, adjust WI if an
436* offdiagonal element of a 2-by-2 block in the Schur form
437* underflows.
438*
439 IF( ieval.GT.0 ) THEN
440 i1 = ieval + 1
441 i2 = ihi - 1
442 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi,
443 $ max( ilo-1, 1 ), ierr )
444 ELSE IF( wantst ) THEN
445 i1 = 1
446 i2 = n - 1
447 ELSE
448 i1 = ilo
449 i2 = ihi - 1
450 END IF
451 inxt = i1 - 1
452 DO 20 i = i1, i2
453 IF( i.LT.inxt )
454 $ GO TO 20
455 IF( wi( i ).EQ.zero ) THEN
456 inxt = i + 1
457 ELSE
458 IF( a( i+1, i ).EQ.zero ) THEN
459 wi( i ) = zero
460 wi( i+1 ) = zero
461 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
462 $ zero ) THEN
463 wi( i ) = zero
464 wi( i+1 ) = zero
465 IF( i.GT.1 )
466 $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
467 IF( n.GT.i+1 )
468 $ CALL dswap( n-i-1, a( i, i+2 ), lda,
469 $ a( i+1, i+2 ), lda )
470 IF( wantvs ) THEN
471 CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
472 END IF
473 a( i, i+1 ) = a( i+1, i )
474 a( i+1, i ) = zero
475 END IF
476 inxt = i + 2
477 END IF
478 20 CONTINUE
479 END IF
480*
481* Undo scaling for the imaginary part of the eigenvalues
482*
483 CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
484 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
485 END IF
486*
487 IF( wantst .AND. info.EQ.0 ) THEN
488*
489* Check if reordering successful
490*
491 lastsl = .true.
492 lst2sl = .true.
493 sdim = 0
494 ip = 0
495 DO 30 i = 1, n
496 cursl = SELECT( wr( i ), wi( i ) )
497 IF( wi( i ).EQ.zero ) THEN
498 IF( cursl )
499 $ sdim = sdim + 1
500 ip = 0
501 IF( cursl .AND. .NOT.lastsl )
502 $ info = n + 2
503 ELSE
504 IF( ip.EQ.1 ) THEN
505*
506* Last eigenvalue of conjugate pair
507*
508 cursl = cursl .OR. lastsl
509 lastsl = cursl
510 IF( cursl )
511 $ sdim = sdim + 2
512 ip = -1
513 IF( cursl .AND. .NOT.lst2sl )
514 $ info = n + 2
515 ELSE
516*
517* First eigenvalue of conjugate pair
518*
519 ip = 1
520 END IF
521 END IF
522 lst2sl = lastsl
523 lastsl = cursl
524 30 CONTINUE
525 END IF
526*
527 work( 1 ) = maxwrk
528 RETURN
529*
530* End of DGEES
531*
532 END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:160
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:130
subroutine dgees(JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI, VS, LDVS, WORK, LWORK, BWORK, INFO)
DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition: dgees.f:216
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:316
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
subroutine dtrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
DTRSEN
Definition: dtrsen.f:313