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