LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
212 *
213 *> \ingroup doubleGEeigen
214 *
215 * =====================================================================
216  SUBROUTINE dgees( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
217  $ vs, ldvs, work, lwork, bwork, info )
218 *
219 * -- LAPACK driver routine (version 3.4.0) --
220 * -- LAPACK is a software package provided by Univ. of Tennessee, --
221 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222 * November 2011
223 *
224 * .. Scalar Arguments ..
225  CHARACTER jobvs, sort
226  INTEGER info, lda, ldvs, lwork, n, sdim
227 * ..
228 * .. Array Arguments ..
229  LOGICAL bwork( * )
230  DOUBLE PRECISION a( lda, * ), vs( ldvs, * ), wi( * ), work( * ),
231  $ wr( * )
232 * ..
233 * .. Function Arguments ..
234  LOGICAL select
235  EXTERNAL SELECT
236 * ..
237 *
238 * =====================================================================
239 *
240 * .. Parameters ..
241  DOUBLE PRECISION zero, one
242  parameter( zero = 0.0d0, one = 1.0d0 )
243 * ..
244 * .. Local Scalars ..
245  LOGICAL cursl, lastsl, lquery, lst2sl, scalea, wantst,
246  $ wantvs
247  INTEGER hswork, i, i1, i2, ibal, icond, ierr, ieval,
248  $ ihi, ilo, inxt, ip, itau, iwrk, maxwrk, minwrk
249  DOUBLE PRECISION anrm, bignum, cscale, eps, s, sep, smlnum
250 * ..
251 * .. Local Arrays ..
252  INTEGER idum( 1 )
253  DOUBLE PRECISION dum( 1 )
254 * ..
255 * .. External Subroutines ..
256  EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
258 * ..
259 * .. External Functions ..
260  LOGICAL lsame
261  INTEGER ilaenv
262  DOUBLE PRECISION dlamch, dlange
263  EXTERNAL lsame, ilaenv, dlamch, dlange
264 * ..
265 * .. Intrinsic Functions ..
266  INTRINSIC max, sqrt
267 * ..
268 * .. Executable Statements ..
269 *
270 * Test the input arguments
271 *
272  info = 0
273  lquery = ( lwork.EQ.-1 )
274  wantvs = lsame( jobvs, 'V' )
275  wantst = lsame( sort, 'S' )
276  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
277  info = -1
278  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
279  info = -2
280  ELSE IF( n.LT.0 ) THEN
281  info = -4
282  ELSE IF( lda.LT.max( 1, n ) ) THEN
283  info = -6
284  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
285  info = -11
286  END IF
287 *
288 * Compute workspace
289 * (Note: Comments in the code beginning "Workspace:" describe the
290 * minimal amount of workspace needed at that point in the code,
291 * as well as the preferred amount for good performance.
292 * NB refers to the optimal block size for the immediately
293 * following subroutine, as returned by ILAENV.
294 * HSWORK refers to the workspace preferred by DHSEQR, as
295 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
296 * the worst case.)
297 *
298  IF( info.EQ.0 ) THEN
299  IF( n.EQ.0 ) THEN
300  minwrk = 1
301  maxwrk = 1
302  ELSE
303  maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
304  minwrk = 3*n
305 *
306  CALL dhseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
307  $ work, -1, ieval )
308  hswork = work( 1 )
309 *
310  IF( .NOT.wantvs ) THEN
311  maxwrk = max( maxwrk, n + hswork )
312  ELSE
313  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
314  $ 'DORGHR', ' ', n, 1, n, -1 ) )
315  maxwrk = max( maxwrk, n + hswork )
316  END IF
317  END IF
318  work( 1 ) = maxwrk
319 *
320  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
321  info = -13
322  END IF
323  END IF
324 *
325  IF( info.NE.0 ) THEN
326  CALL xerbla( 'DGEES ', -info )
327  return
328  ELSE IF( lquery ) THEN
329  return
330  END IF
331 *
332 * Quick return if possible
333 *
334  IF( n.EQ.0 ) THEN
335  sdim = 0
336  return
337  END IF
338 *
339 * Get machine constants
340 *
341  eps = dlamch( 'P' )
342  smlnum = dlamch( 'S' )
343  bignum = one / smlnum
344  CALL dlabad( smlnum, bignum )
345  smlnum = sqrt( smlnum ) / eps
346  bignum = one / smlnum
347 *
348 * Scale A if max element outside range [SMLNUM,BIGNUM]
349 *
350  anrm = dlange( 'M', n, n, a, lda, dum )
351  scalea = .false.
352  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
353  scalea = .true.
354  cscale = smlnum
355  ELSE IF( anrm.GT.bignum ) THEN
356  scalea = .true.
357  cscale = bignum
358  END IF
359  IF( scalea )
360  $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
361 *
362 * Permute the matrix to make it more nearly triangular
363 * (Workspace: need N)
364 *
365  ibal = 1
366  CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
367 *
368 * Reduce to upper Hessenberg form
369 * (Workspace: need 3*N, prefer 2*N+N*NB)
370 *
371  itau = n + ibal
372  iwrk = n + itau
373  CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
374  $ lwork-iwrk+1, ierr )
375 *
376  IF( wantvs ) THEN
377 *
378 * Copy Householder vectors to VS
379 *
380  CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
381 *
382 * Generate orthogonal matrix in VS
383 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
384 *
385  CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
386  $ lwork-iwrk+1, ierr )
387  END IF
388 *
389  sdim = 0
390 *
391 * Perform QR iteration, accumulating Schur vectors in VS if desired
392 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
393 *
394  iwrk = itau
395  CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
396  $ work( iwrk ), lwork-iwrk+1, ieval )
397  IF( ieval.GT.0 )
398  $ info = ieval
399 *
400 * Sort eigenvalues if desired
401 *
402  IF( wantst .AND. info.EQ.0 ) THEN
403  IF( scalea ) THEN
404  CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
405  CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
406  END IF
407  DO 10 i = 1, n
408  bwork( i ) = SELECT( wr( i ), wi( i ) )
409  10 continue
410 *
411 * Reorder eigenvalues and transform Schur vectors
412 * (Workspace: none needed)
413 *
414  CALL dtrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
415  $ sdim, s, sep, work( iwrk ), lwork-iwrk+1, idum, 1,
416  $ icond )
417  IF( icond.GT.0 )
418  $ info = n + icond
419  END IF
420 *
421  IF( wantvs ) THEN
422 *
423 * Undo balancing
424 * (Workspace: need N)
425 *
426  CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
427  $ ierr )
428  END IF
429 *
430  IF( scalea ) THEN
431 *
432 * Undo scaling for the Schur form of A
433 *
434  CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
435  CALL dcopy( n, a, lda+1, wr, 1 )
436  IF( cscale.EQ.smlnum ) THEN
437 *
438 * If scaling back towards underflow, adjust WI if an
439 * offdiagonal element of a 2-by-2 block in the Schur form
440 * underflows.
441 *
442  IF( ieval.GT.0 ) THEN
443  i1 = ieval + 1
444  i2 = ihi - 1
445  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi,
446  $ max( ilo-1, 1 ), ierr )
447  ELSE IF( wantst ) THEN
448  i1 = 1
449  i2 = n - 1
450  ELSE
451  i1 = ilo
452  i2 = ihi - 1
453  END IF
454  inxt = i1 - 1
455  DO 20 i = i1, i2
456  IF( i.LT.inxt )
457  $ go to 20
458  IF( wi( i ).EQ.zero ) THEN
459  inxt = i + 1
460  ELSE
461  IF( a( i+1, i ).EQ.zero ) THEN
462  wi( i ) = zero
463  wi( i+1 ) = zero
464  ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
465  $ zero ) THEN
466  wi( i ) = zero
467  wi( i+1 ) = zero
468  IF( i.GT.1 )
469  $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
470  IF( n.GT.i+1 )
471  $ CALL dswap( n-i-1, a( i, i+2 ), lda,
472  $ a( i+1, i+2 ), lda )
473  IF( wantvs ) THEN
474  CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
475  END IF
476  a( i, i+1 ) = a( i+1, i )
477  a( i+1, i ) = zero
478  END IF
479  inxt = i + 2
480  END IF
481  20 continue
482  END IF
483 *
484 * Undo scaling for the imaginary part of the eigenvalues
485 *
486  CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
487  $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
488  END IF
489 *
490  IF( wantst .AND. info.EQ.0 ) THEN
491 *
492 * Check if reordering successful
493 *
494  lastsl = .true.
495  lst2sl = .true.
496  sdim = 0
497  ip = 0
498  DO 30 i = 1, n
499  cursl = SELECT( wr( i ), wi( i ) )
500  IF( wi( i ).EQ.zero ) THEN
501  IF( cursl )
502  $ sdim = sdim + 1
503  ip = 0
504  IF( cursl .AND. .NOT.lastsl )
505  $ info = n + 2
506  ELSE
507  IF( ip.EQ.1 ) THEN
508 *
509 * Last eigenvalue of conjugate pair
510 *
511  cursl = cursl .OR. lastsl
512  lastsl = cursl
513  IF( cursl )
514  $ sdim = sdim + 2
515  ip = -1
516  IF( cursl .AND. .NOT.lst2sl )
517  $ info = n + 2
518  ELSE
519 *
520 * First eigenvalue of conjugate pair
521 *
522  ip = 1
523  END IF
524  END IF
525  lst2sl = lastsl
526  lastsl = cursl
527  30 continue
528  END IF
529 *
530  work( 1 ) = maxwrk
531  return
532 *
533 * End of DGEES
534 *
535  END