LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
sgees.f
Go to the documentation of this file.
1 *> \brief <b> SGEES 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 SGEES + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgees.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgees.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgees.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGEES( 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 * REAL 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 *> SGEES 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 LOGICAL FUNCTION of two REAL 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 REAL 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 REAL array, dimension (N)
131 *> \endverbatim
132 *>
133 *> \param[out] WI
134 *> \verbatim
135 *> WI is REAL 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 REAL 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 REAL 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 realGEeigen
214 *
215 * =====================================================================
216  SUBROUTINE sgees( 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  REAL 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  REAL ZERO, ONE
242  parameter( zero = 0.0e0, one = 1.0e0 )
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  REAL ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
250 * ..
251 * .. Local Arrays ..
252  INTEGER IDUM( 1 )
253  REAL DUM( 1 )
254 * ..
255 * .. External Subroutines ..
256  EXTERNAL scopy, sgebak, sgebal, sgehrd, shseqr, slabad,
258 * ..
259 * .. External Functions ..
260  LOGICAL LSAME
261  INTEGER ILAENV
262  REAL SLAMCH, SLANGE
263  EXTERNAL lsame, ilaenv, slamch, slange
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 SHSEQR, 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, 'SGEHRD', ' ', n, 1, n, 0 )
304  minwrk = 3*n
305 *
306  CALL shseqr( '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  $ 'SORGHR', ' ', 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( 'SGEES ', -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 = slamch( 'P' )
342  smlnum = slamch( 'S' )
343  bignum = one / smlnum
344  CALL slabad( 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 = slange( '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 slascl( '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 sgebal( '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 sgehrd( 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 slacpy( '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 sorghr( 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 shseqr( '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 slascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
405  CALL slascl( '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 strsen( '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 sgebak( '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 slascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
435  CALL scopy( 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 slascl( '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 sswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
470  IF( n.GT.i+1 )
471  $ CALL sswap( n-i-1, a( i, i+2 ), lda,
472  $ a( i+1, i+2 ), lda )
473  IF( wantvs ) THEN
474  CALL sswap( 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 slascl( '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 SGEES
534 *
535  END
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
SGEBAK
Definition: sgebak.f:132
subroutine strsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
STRSEN
Definition: strsen.f:316
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:318
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:169
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine sgees(JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI, VS, LDVS, WORK, LWORK, BWORK, INFO)
SGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: sgees.f:218
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
SGEBAL
Definition: sgebal.f:162
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:141
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:128