LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
cgees.f
Go to the documentation of this file.
1 *> \brief <b> CGEES 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgees.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgees.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgees.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
22 * LDVS, WORK, LWORK, RWORK, 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 RWORK( * )
31 * COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
32 * ..
33 * .. Function Arguments ..
34 * LOGICAL SELECT
35 * EXTERNAL SELECT
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> CGEES computes for an N-by-N complex nonsymmetric matrix A, the
45 *> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
46 *> vectors Z. This gives the Schur factorization A = Z*T*(Z**H).
47 *>
48 *> Optionally, it also orders the eigenvalues on the diagonal of the
49 *> 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 complex matrix is in Schur form if it is upper triangular.
54 *> \endverbatim
55 *
56 * Arguments:
57 * ==========
58 *
59 *> \param[in] JOBVS
60 *> \verbatim
61 *> JOBVS is CHARACTER*1
62 *> = 'N': Schur vectors are not computed;
63 *> = 'V': Schur vectors are computed.
64 *> \endverbatim
65 *>
66 *> \param[in] SORT
67 *> \verbatim
68 *> SORT is CHARACTER*1
69 *> Specifies whether or not to order the eigenvalues on the
70 *> diagonal of the Schur form.
71 *> = 'N': Eigenvalues are not ordered:
72 *> = 'S': Eigenvalues are ordered (see SELECT).
73 *> \endverbatim
74 *>
75 *> \param[in] SELECT
76 *> \verbatim
77 *> SELECT is a LOGICAL FUNCTION of one COMPLEX argument
78 *> SELECT must be declared EXTERNAL in the calling subroutine.
79 *> If SORT = 'S', SELECT is used to select eigenvalues to order
80 *> to the top left of the Schur form.
81 *> IF SORT = 'N', SELECT is not referenced.
82 *> The eigenvalue W(j) is selected if SELECT(W(j)) is true.
83 *> \endverbatim
84 *>
85 *> \param[in] N
86 *> \verbatim
87 *> N is INTEGER
88 *> The order of the matrix A. N >= 0.
89 *> \endverbatim
90 *>
91 *> \param[in,out] A
92 *> \verbatim
93 *> A is COMPLEX array, dimension (LDA,N)
94 *> On entry, the N-by-N matrix A.
95 *> On exit, A has been overwritten by its Schur form T.
96 *> \endverbatim
97 *>
98 *> \param[in] LDA
99 *> \verbatim
100 *> LDA is INTEGER
101 *> The leading dimension of the array A. LDA >= max(1,N).
102 *> \endverbatim
103 *>
104 *> \param[out] SDIM
105 *> \verbatim
106 *> SDIM is INTEGER
107 *> If SORT = 'N', SDIM = 0.
108 *> If SORT = 'S', SDIM = number of eigenvalues for which
109 *> SELECT is true.
110 *> \endverbatim
111 *>
112 *> \param[out] W
113 *> \verbatim
114 *> W is COMPLEX array, dimension (N)
115 *> W contains the computed eigenvalues, in the same order that
116 *> they appear on the diagonal of the output Schur form T.
117 *> \endverbatim
118 *>
119 *> \param[out] VS
120 *> \verbatim
121 *> VS is COMPLEX array, dimension (LDVS,N)
122 *> If JOBVS = 'V', VS contains the unitary matrix Z of Schur
123 *> vectors.
124 *> If JOBVS = 'N', VS is not referenced.
125 *> \endverbatim
126 *>
127 *> \param[in] LDVS
128 *> \verbatim
129 *> LDVS is INTEGER
130 *> The leading dimension of the array VS. LDVS >= 1; if
131 *> JOBVS = 'V', LDVS >= N.
132 *> \endverbatim
133 *>
134 *> \param[out] WORK
135 *> \verbatim
136 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
137 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
138 *> \endverbatim
139 *>
140 *> \param[in] LWORK
141 *> \verbatim
142 *> LWORK is INTEGER
143 *> The dimension of the array WORK. LWORK >= max(1,2*N).
144 *> For good performance, LWORK must generally be larger.
145 *>
146 *> If LWORK = -1, then a workspace query is assumed; the routine
147 *> only calculates the optimal size of the WORK array, returns
148 *> this value as the first entry of the WORK array, and no error
149 *> message related to LWORK is issued by XERBLA.
150 *> \endverbatim
151 *>
152 *> \param[out] RWORK
153 *> \verbatim
154 *> RWORK is REAL array, dimension (N)
155 *> \endverbatim
156 *>
157 *> \param[out] BWORK
158 *> \verbatim
159 *> BWORK is LOGICAL array, dimension (N)
160 *> Not referenced if SORT = 'N'.
161 *> \endverbatim
162 *>
163 *> \param[out] INFO
164 *> \verbatim
165 *> INFO is INTEGER
166 *> = 0: successful exit
167 *> < 0: if INFO = -i, the i-th argument had an illegal value.
168 *> > 0: if INFO = i, and i is
169 *> <= N: the QR algorithm failed to compute all the
170 *> eigenvalues; elements 1:ILO-1 and i+1:N of W
171 *> contain those eigenvalues which have converged;
172 *> if JOBVS = 'V', VS contains the matrix which
173 *> reduces A to its partially converged Schur form.
174 *> = N+1: the eigenvalues could not be reordered because
175 *> some eigenvalues were too close to separate (the
176 *> problem is very ill-conditioned);
177 *> = N+2: after reordering, roundoff changed values of
178 *> some complex eigenvalues so that leading
179 *> eigenvalues in the Schur form no longer satisfy
180 *> SELECT = .TRUE.. This could also be caused by
181 *> underflow due to scaling.
182 *> \endverbatim
183 *
184 * Authors:
185 * ========
186 *
187 *> \author Univ. of Tennessee
188 *> \author Univ. of California Berkeley
189 *> \author Univ. of Colorado Denver
190 *> \author NAG Ltd.
191 *
192 *> \date November 2011
193 *
194 *> \ingroup complexGEeigen
195 *
196 * =====================================================================
197  SUBROUTINE cgees( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
198  \$ ldvs, work, lwork, rwork, bwork, info )
199 *
200 * -- LAPACK driver routine (version 3.4.0) --
201 * -- LAPACK is a software package provided by Univ. of Tennessee, --
202 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203 * November 2011
204 *
205 * .. Scalar Arguments ..
206  CHARACTER jobvs, sort
207  INTEGER info, lda, ldvs, lwork, n, sdim
208 * ..
209 * .. Array Arguments ..
210  LOGICAL bwork( * )
211  REAL rwork( * )
212  COMPLEX a( lda, * ), vs( ldvs, * ), w( * ), work( * )
213 * ..
214 * .. Function Arguments ..
215  LOGICAL select
216  EXTERNAL SELECT
217 * ..
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222  REAL zero, one
223  parameter( zero = 0.0e0, one = 1.0e0 )
224 * ..
225 * .. Local Scalars ..
226  LOGICAL lquery, scalea, wantst, wantvs
227  INTEGER hswork, i, ibal, icond, ierr, ieval, ihi, ilo,
228  \$ itau, iwrk, maxwrk, minwrk
229  REAL anrm, bignum, cscale, eps, s, sep, smlnum
230 * ..
231 * .. Local Arrays ..
232  REAL dum( 1 )
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL ccopy, cgebak, cgebal, cgehrd, chseqr, clacpy,
237 * ..
238 * .. External Functions ..
239  LOGICAL lsame
240  INTEGER ilaenv
241  REAL clange, slamch
242  EXTERNAL lsame, ilaenv, clange, slamch
243 * ..
244 * .. Intrinsic Functions ..
245  INTRINSIC max, sqrt
246 * ..
247 * .. Executable Statements ..
248 *
249 * Test the input arguments
250 *
251  info = 0
252  lquery = ( lwork.EQ.-1 )
253  wantvs = lsame( jobvs, 'V' )
254  wantst = lsame( sort, 'S' )
255  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
256  info = -1
257  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
258  info = -2
259  ELSE IF( n.LT.0 ) THEN
260  info = -4
261  ELSE IF( lda.LT.max( 1, n ) ) THEN
262  info = -6
263  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
264  info = -10
265  END IF
266 *
267 * Compute workspace
268 * (Note: Comments in the code beginning "Workspace:" describe the
269 * minimal amount of workspace needed at that point in the code,
270 * as well as the preferred amount for good performance.
271 * CWorkspace refers to complex workspace, and RWorkspace to real
272 * workspace. NB refers to the optimal block size for the
273 * immediately following subroutine, as returned by ILAENV.
274 * HSWORK refers to the workspace preferred by CHSEQR, as
275 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
276 * the worst case.)
277 *
278  IF( info.EQ.0 ) THEN
279  IF( n.EQ.0 ) THEN
280  minwrk = 1
281  maxwrk = 1
282  ELSE
283  maxwrk = n + n*ilaenv( 1, 'CGEHRD', ' ', n, 1, n, 0 )
284  minwrk = 2*n
285 *
286  CALL chseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
287  \$ work, -1, ieval )
288  hswork = work( 1 )
289 *
290  IF( .NOT.wantvs ) THEN
291  maxwrk = max( maxwrk, hswork )
292  ELSE
293  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'CUNGHR',
294  \$ ' ', n, 1, n, -1 ) )
295  maxwrk = max( maxwrk, hswork )
296  END IF
297  END IF
298  work( 1 ) = maxwrk
299 *
300  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
301  info = -12
302  END IF
303  END IF
304 *
305  IF( info.NE.0 ) THEN
306  CALL xerbla( 'CGEES ', -info )
307  return
308  ELSE IF( lquery ) THEN
309  return
310  END IF
311 *
312 * Quick return if possible
313 *
314  IF( n.EQ.0 ) THEN
315  sdim = 0
316  return
317  END IF
318 *
319 * Get machine constants
320 *
321  eps = slamch( 'P' )
322  smlnum = slamch( 'S' )
323  bignum = one / smlnum
324  CALL slabad( smlnum, bignum )
325  smlnum = sqrt( smlnum ) / eps
326  bignum = one / smlnum
327 *
328 * Scale A if max element outside range [SMLNUM,BIGNUM]
329 *
330  anrm = clange( 'M', n, n, a, lda, dum )
331  scalea = .false.
332  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
333  scalea = .true.
334  cscale = smlnum
335  ELSE IF( anrm.GT.bignum ) THEN
336  scalea = .true.
337  cscale = bignum
338  END IF
339  IF( scalea )
340  \$ CALL clascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
341 *
342 * Permute the matrix to make it more nearly triangular
343 * (CWorkspace: none)
344 * (RWorkspace: need N)
345 *
346  ibal = 1
347  CALL cgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
348 *
349 * Reduce to upper Hessenberg form
350 * (CWorkspace: need 2*N, prefer N+N*NB)
351 * (RWorkspace: none)
352 *
353  itau = 1
354  iwrk = n + itau
355  CALL cgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
356  \$ lwork-iwrk+1, ierr )
357 *
358  IF( wantvs ) THEN
359 *
360 * Copy Householder vectors to VS
361 *
362  CALL clacpy( 'L', n, n, a, lda, vs, ldvs )
363 *
364 * Generate unitary matrix in VS
365 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
366 * (RWorkspace: none)
367 *
368  CALL cunghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
369  \$ lwork-iwrk+1, ierr )
370  END IF
371 *
372  sdim = 0
373 *
374 * Perform QR iteration, accumulating Schur vectors in VS if desired
375 * (CWorkspace: need 1, prefer HSWORK (see comments) )
376 * (RWorkspace: none)
377 *
378  iwrk = itau
379  CALL chseqr( 'S', jobvs, n, ilo, ihi, a, lda, w, vs, ldvs,
380  \$ work( iwrk ), lwork-iwrk+1, ieval )
381  IF( ieval.GT.0 )
382  \$ info = ieval
383 *
384 * Sort eigenvalues if desired
385 *
386  IF( wantst .AND. info.EQ.0 ) THEN
387  IF( scalea )
388  \$ CALL clascl( 'G', 0, 0, cscale, anrm, n, 1, w, n, ierr )
389  DO 10 i = 1, n
390  bwork( i ) = SELECT( w( i ) )
391  10 continue
392 *
393 * Reorder eigenvalues and transform Schur vectors
394 * (CWorkspace: none)
395 * (RWorkspace: none)
396 *
397  CALL ctrsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, w, sdim,
398  \$ s, sep, work( iwrk ), lwork-iwrk+1, icond )
399  END IF
400 *
401  IF( wantvs ) THEN
402 *
403 * Undo balancing
404 * (CWorkspace: none)
405 * (RWorkspace: need N)
406 *
407  CALL cgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs, ldvs,
408  \$ ierr )
409  END IF
410 *
411  IF( scalea ) THEN
412 *
413 * Undo scaling for the Schur form of A
414 *
415  CALL clascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
416  CALL ccopy( n, a, lda+1, w, 1 )
417  END IF
418 *
419  work( 1 ) = maxwrk
420  return
421 *
422 * End of CGEES
423 *
424  END