LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> Download CGEESX + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeesx.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeesx.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeesx.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
20* VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
21* BWORK, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER JOBVS, SENSE, SORT
25* INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
26* REAL RCONDE, RCONDV
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*> CGEESX 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*> computes a reciprocal condition number for the average of the
51*> selected eigenvalues (RCONDE); and computes a reciprocal condition
52*> number for the right invariant subspace corresponding to the
53*> selected eigenvalues (RCONDV). The leading columns of Z form an
54*> orthonormal basis for this invariant subspace.
55*>
56*> For further explanation of the reciprocal condition numbers RCONDE
57*> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
58*> these quantities are called s and sep respectively).
59*>
60*> A complex matrix is in Schur form if it is upper triangular.
61*> \endverbatim
62*
63* Arguments:
64* ==========
65*
66*> \param[in] JOBVS
67*> \verbatim
68*> JOBVS is CHARACTER*1
69*> = 'N': Schur vectors are not computed;
70*> = 'V': Schur vectors are computed.
71*> \endverbatim
72*>
73*> \param[in] SORT
74*> \verbatim
75*> SORT is CHARACTER*1
76*> Specifies whether or not to order the eigenvalues on the
77*> diagonal of the Schur form.
78*> = 'N': Eigenvalues are not ordered;
79*> = 'S': Eigenvalues are ordered (see SELECT).
80*> \endverbatim
81*>
82*> \param[in] SELECT
83*> \verbatim
84*> SELECT is a LOGICAL FUNCTION of one COMPLEX argument
85*> SELECT must be declared EXTERNAL in the calling subroutine.
86*> If SORT = 'S', SELECT is used to select eigenvalues to order
87*> to the top left of the Schur form.
88*> If SORT = 'N', SELECT is not referenced.
89*> An eigenvalue W(j) is selected if SELECT(W(j)) is true.
90*> \endverbatim
91*>
92*> \param[in] SENSE
93*> \verbatim
94*> SENSE is CHARACTER*1
95*> Determines which reciprocal condition numbers are computed.
96*> = 'N': None are computed;
97*> = 'E': Computed for average of selected eigenvalues only;
98*> = 'V': Computed for selected right invariant subspace only;
99*> = 'B': Computed for both.
100*> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
101*> \endverbatim
102*>
103*> \param[in] N
104*> \verbatim
105*> N is INTEGER
106*> The order of the matrix A. N >= 0.
107*> \endverbatim
108*>
109*> \param[in,out] A
110*> \verbatim
111*> A is COMPLEX array, dimension (LDA, N)
112*> On entry, the N-by-N matrix A.
113*> On exit, A is overwritten by its Schur form T.
114*> \endverbatim
115*>
116*> \param[in] LDA
117*> \verbatim
118*> LDA is INTEGER
119*> The leading dimension of the array A. LDA >= max(1,N).
120*> \endverbatim
121*>
122*> \param[out] SDIM
123*> \verbatim
124*> SDIM is INTEGER
125*> If SORT = 'N', SDIM = 0.
126*> If SORT = 'S', SDIM = number of eigenvalues for which
127*> SELECT is true.
128*> \endverbatim
129*>
130*> \param[out] W
131*> \verbatim
132*> W is COMPLEX array, dimension (N)
133*> W contains the computed eigenvalues, in the same order
134*> that they appear on the diagonal of the output Schur form T.
135*> \endverbatim
136*>
137*> \param[out] VS
138*> \verbatim
139*> VS is COMPLEX array, dimension (LDVS,N)
140*> If JOBVS = 'V', VS contains the unitary matrix Z of Schur
141*> vectors.
142*> If JOBVS = 'N', VS is not referenced.
143*> \endverbatim
144*>
145*> \param[in] LDVS
146*> \verbatim
147*> LDVS is INTEGER
148*> The leading dimension of the array VS. LDVS >= 1, and if
149*> JOBVS = 'V', LDVS >= N.
150*> \endverbatim
151*>
152*> \param[out] RCONDE
153*> \verbatim
154*> RCONDE is REAL
155*> If SENSE = 'E' or 'B', RCONDE contains the reciprocal
156*> condition number for the average of the selected eigenvalues.
157*> Not referenced if SENSE = 'N' or 'V'.
158*> \endverbatim
159*>
160*> \param[out] RCONDV
161*> \verbatim
162*> RCONDV is REAL
163*> If SENSE = 'V' or 'B', RCONDV contains the reciprocal
164*> condition number for the selected right invariant subspace.
165*> Not referenced if SENSE = 'N' or 'E'.
166*> \endverbatim
167*>
168*> \param[out] WORK
169*> \verbatim
170*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
171*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
172*> \endverbatim
173*>
174*> \param[in] LWORK
175*> \verbatim
176*> LWORK is INTEGER
177*> The dimension of the array WORK. LWORK >= max(1,2*N).
178*> Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
179*> where SDIM is the number of selected eigenvalues computed by
180*> this routine. Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also
181*> that an error is only returned if LWORK < max(1,2*N), but if
182*> SENSE = 'E' or 'V' or 'B' this may not be large enough.
183*> For good performance, LWORK must generally be larger.
184*>
185*> If LWORK = -1, then a workspace query is assumed; the routine
186*> only calculates upper bound on the optimal size of the
187*> array WORK, returns this value as the first entry of the WORK
188*> array, and no error message related to LWORK is issued by
189*> XERBLA.
190*> \endverbatim
191*>
192*> \param[out] RWORK
193*> \verbatim
194*> RWORK is REAL array, dimension (N)
195*> \endverbatim
196*>
197*> \param[out] BWORK
198*> \verbatim
199*> BWORK is LOGICAL array, dimension (N)
200*> Not referenced if SORT = 'N'.
201*> \endverbatim
202*>
203*> \param[out] INFO
204*> \verbatim
205*> INFO is INTEGER
206*> = 0: successful exit
207*> < 0: if INFO = -i, the i-th argument had an illegal value.
208*> > 0: if INFO = i, and i is
209*> <= N: the QR algorithm failed to compute all the
210*> eigenvalues; elements 1:ILO-1 and i+1:N of W
211*> contain those eigenvalues which have converged; if
212*> JOBVS = 'V', VS contains the transformation which
213*> reduces A to its partially converged Schur form.
214*> = N+1: the eigenvalues could not be reordered because some
215*> eigenvalues were too close to separate (the problem
216*> is very ill-conditioned);
217*> = N+2: after reordering, roundoff changed values of some
218*> complex eigenvalues so that leading eigenvalues in
219*> the Schur form no longer satisfy SELECT=.TRUE. This
220*> could also be caused by underflow due to scaling.
221*> \endverbatim
222*
223* Authors:
224* ========
225*
226*> \author Univ. of Tennessee
227*> \author Univ. of California Berkeley
228*> \author Univ. of Colorado Denver
229*> \author NAG Ltd.
230*
231*> \ingroup geesx
232*
233* =====================================================================
234 SUBROUTINE cgeesx( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
235 $ W,
236 $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
237 $ BWORK, INFO )
238*
239* -- LAPACK driver routine --
240* -- LAPACK is a software package provided by Univ. of Tennessee, --
241* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242*
243* .. Scalar Arguments ..
244 CHARACTER JOBVS, SENSE, SORT
245 INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
246 REAL RCONDE, RCONDV
247* ..
248* .. Array Arguments ..
249 LOGICAL BWORK( * )
250 REAL RWORK( * )
251 COMPLEX A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
252* ..
253* .. Function Arguments ..
254 LOGICAL SELECT
255 EXTERNAL SELECT
256* ..
257*
258* =====================================================================
259*
260* .. Parameters ..
261 REAL ZERO, ONE
262 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
263* ..
264* .. Local Scalars ..
265 LOGICAL LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST,
266 $ WANTSV, WANTVS
267 INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
268 $ itau, iwrk, lwrk, maxwrk, minwrk
269 REAL ANRM, BIGNUM, CSCALE, EPS, SMLNUM
270* ..
271* .. Local Arrays ..
272 REAL DUM( 1 )
273* ..
274* .. External Subroutines ..
275 EXTERNAL ccopy, cgebak, cgebal, cgehrd, chseqr,
276 $ clacpy,
278* ..
279* .. External Functions ..
280 LOGICAL LSAME
281 INTEGER ILAENV
282 REAL CLANGE, SLAMCH, SROUNDUP_LWORK
283 EXTERNAL lsame, ilaenv, clange, slamch,
284 $ sroundup_lwork
285* ..
286* .. Intrinsic Functions ..
287 INTRINSIC max, sqrt
288* ..
289* .. Executable Statements ..
290*
291* Test the input arguments
292*
293 info = 0
294 wantvs = lsame( jobvs, 'V' )
295 wantst = lsame( sort, 'S' )
296 wantsn = lsame( sense, 'N' )
297 wantse = lsame( sense, 'E' )
298 wantsv = lsame( sense, 'V' )
299 wantsb = lsame( sense, 'B' )
300 lquery = ( lwork.EQ.-1 )
301*
302 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
303 info = -1
304 ELSE IF( ( .NOT.wantst ) .AND.
305 $ ( .NOT.lsame( sort, 'N' ) ) ) THEN
306 info = -2
307 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
308 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
309 info = -4
310 ELSE IF( n.LT.0 ) THEN
311 info = -5
312 ELSE IF( lda.LT.max( 1, n ) ) THEN
313 info = -7
314 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
315 info = -11
316 END IF
317*
318* Compute workspace
319* (Note: Comments in the code beginning "Workspace:" describe the
320* minimal amount of real workspace needed at that point in the
321* code, as well as the preferred amount for good performance.
322* CWorkspace refers to complex workspace, and RWorkspace to real
323* workspace. NB refers to the optimal block size for the
324* immediately following subroutine, as returned by ILAENV.
325* HSWORK refers to the workspace preferred by CHSEQR, as
326* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
327* the worst case.
328* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
329* depends on SDIM, which is computed by the routine CTRSEN later
330* in the code.)
331*
332 IF( info.EQ.0 ) THEN
333 IF( n.EQ.0 ) THEN
334 minwrk = 1
335 lwrk = 1
336 ELSE
337 maxwrk = n + n*ilaenv( 1, 'CGEHRD', ' ', n, 1, n, 0 )
338 minwrk = 2*n
339*
340 CALL chseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
341 $ work, -1, ieval )
342 hswork = int( work( 1 ) )
343*
344 IF( .NOT.wantvs ) THEN
345 maxwrk = max( maxwrk, hswork )
346 ELSE
347 maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1,
348 $ '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 ) = sroundup_lwork(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 smlnum = sqrt( smlnum ) / eps
383 bignum = one / smlnum
384*
385* Scale A if max element outside range [SMLNUM,BIGNUM]
386*
387 anrm = clange( 'M', n, n, a, lda, dum )
388 scalea = .false.
389 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
390 scalea = .true.
391 cscale = smlnum
392 ELSE IF( anrm.GT.bignum ) THEN
393 scalea = .true.
394 cscale = bignum
395 END IF
396 IF( scalea )
397 $ CALL clascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
398*
399*
400* Permute the matrix to make it more nearly triangular
401* (CWorkspace: none)
402* (RWorkspace: need N)
403*
404 ibal = 1
405 CALL cgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
406*
407* Reduce to upper Hessenberg form
408* (CWorkspace: need 2*N, prefer N+N*NB)
409* (RWorkspace: none)
410*
411 itau = 1
412 iwrk = n + itau
413 CALL cgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
414 $ lwork-iwrk+1, ierr )
415*
416 IF( wantvs ) THEN
417*
418* Copy Householder vectors to VS
419*
420 CALL clacpy( 'L', n, n, a, lda, vs, ldvs )
421*
422* Generate unitary matrix in VS
423* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
424* (RWorkspace: none)
425*
426 CALL cunghr( n, ilo, ihi, vs, ldvs, work( itau ),
427 $ 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,
459 $ sdim,
460 $ rconde, rcondv, work( iwrk ), lwork-iwrk+1,
461 $ icond )
462 IF( .NOT.wantsn )
463 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
464 IF( icond.EQ.-14 ) THEN
465*
466* Not enough complex workspace
467*
468 info = -15
469 END IF
470 END IF
471*
472 IF( wantvs ) THEN
473*
474* Undo balancing
475* (CWorkspace: none)
476* (RWorkspace: need N)
477*
478 CALL cgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs,
479 $ ldvs,
480 $ ierr )
481 END IF
482*
483 IF( scalea ) THEN
484*
485* Undo scaling for the Schur form of A
486*
487 CALL clascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
488 CALL ccopy( n, a, lda+1, w, 1 )
489 IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
490 dum( 1 ) = rcondv
491 CALL slascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1,
492 $ ierr )
493 rcondv = dum( 1 )
494 END IF
495 END IF
496*
497 work( 1 ) = sroundup_lwork(maxwrk)
498 RETURN
499*
500* End of CGEESX
501*
502 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
CGEBAK
Definition cgebak.f:129
subroutine cgebal(job, n, a, lda, ilo, ihi, scale, info)
CGEBAL
Definition cgebal.f:163
subroutine cgeesx(jobvs, sort, select, sense, n, a, lda, sdim, w, vs, ldvs, rconde, rcondv, work, lwork, rwork, bwork, info)
CGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition cgeesx.f:238
subroutine cgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
CGEHRD
Definition cgehrd.f:166
subroutine chseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
CHSEQR
Definition chseqr.f:297
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:142
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:142
subroutine ctrsen(job, compq, select, n, t, ldt, q, ldq, w, m, s, sep, work, lwork, info)
CTRSEN
Definition ctrsen.f:263
subroutine cunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
CUNGHR
Definition cunghr.f:125