LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cggev.f
Go to the documentation of this file.
1*> \brief <b> CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors 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 CGGEV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggev.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggev.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggev.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
22* VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBVL, JOBVR
26* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27* ..
28* .. Array Arguments ..
29* REAL RWORK( * )
30* COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
31* $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
32* $ WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> CGGEV computes for a pair of N-by-N complex nonsymmetric matrices
42*> (A,B), the generalized eigenvalues, and optionally, the left and/or
43*> right generalized eigenvectors.
44*>
45*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
46*> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
47*> singular. It is usually represented as the pair (alpha,beta), as
48*> there is a reasonable interpretation for beta=0, and even for both
49*> being zero.
50*>
51*> The right generalized eigenvector v(j) corresponding to the
52*> generalized eigenvalue lambda(j) of (A,B) satisfies
53*>
54*> A * v(j) = lambda(j) * B * v(j).
55*>
56*> The left generalized eigenvector u(j) corresponding to the
57*> generalized eigenvalues lambda(j) of (A,B) satisfies
58*>
59*> u(j)**H * A = lambda(j) * u(j)**H * B
60*>
61*> where u(j)**H is the conjugate-transpose of u(j).
62*> \endverbatim
63*
64* Arguments:
65* ==========
66*
67*> \param[in] JOBVL
68*> \verbatim
69*> JOBVL is CHARACTER*1
70*> = 'N': do not compute the left generalized eigenvectors;
71*> = 'V': compute the left generalized eigenvectors.
72*> \endverbatim
73*>
74*> \param[in] JOBVR
75*> \verbatim
76*> JOBVR is CHARACTER*1
77*> = 'N': do not compute the right generalized eigenvectors;
78*> = 'V': compute the right generalized eigenvectors.
79*> \endverbatim
80*>
81*> \param[in] N
82*> \verbatim
83*> N is INTEGER
84*> The order of the matrices A, B, VL, and VR. N >= 0.
85*> \endverbatim
86*>
87*> \param[in,out] A
88*> \verbatim
89*> A is COMPLEX array, dimension (LDA, N)
90*> On entry, the matrix A in the pair (A,B).
91*> On exit, A has been overwritten.
92*> \endverbatim
93*>
94*> \param[in] LDA
95*> \verbatim
96*> LDA is INTEGER
97*> The leading dimension of A. LDA >= max(1,N).
98*> \endverbatim
99*>
100*> \param[in,out] B
101*> \verbatim
102*> B is COMPLEX array, dimension (LDB, N)
103*> On entry, the matrix B in the pair (A,B).
104*> On exit, B has been overwritten.
105*> \endverbatim
106*>
107*> \param[in] LDB
108*> \verbatim
109*> LDB is INTEGER
110*> The leading dimension of B. LDB >= max(1,N).
111*> \endverbatim
112*>
113*> \param[out] ALPHA
114*> \verbatim
115*> ALPHA is COMPLEX array, dimension (N)
116*> \endverbatim
117*>
118*> \param[out] BETA
119*> \verbatim
120*> BETA is COMPLEX array, dimension (N)
121*> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
122*> generalized eigenvalues.
123*>
124*> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
125*> underflow, and BETA(j) may even be zero. Thus, the user
126*> should avoid naively computing the ratio alpha/beta.
127*> However, ALPHA will be always less than and usually
128*> comparable with norm(A) in magnitude, and BETA always less
129*> than and usually comparable with norm(B).
130*> \endverbatim
131*>
132*> \param[out] VL
133*> \verbatim
134*> VL is COMPLEX array, dimension (LDVL,N)
135*> If JOBVL = 'V', the left generalized eigenvectors u(j) are
136*> stored one after another in the columns of VL, in the same
137*> order as their eigenvalues.
138*> Each eigenvector is scaled so the largest component has
139*> abs(real part) + abs(imag. part) = 1.
140*> Not referenced if JOBVL = 'N'.
141*> \endverbatim
142*>
143*> \param[in] LDVL
144*> \verbatim
145*> LDVL is INTEGER
146*> The leading dimension of the matrix VL. LDVL >= 1, and
147*> if JOBVL = 'V', LDVL >= N.
148*> \endverbatim
149*>
150*> \param[out] VR
151*> \verbatim
152*> VR is COMPLEX array, dimension (LDVR,N)
153*> If JOBVR = 'V', the right generalized eigenvectors v(j) are
154*> stored one after another in the columns of VR, in the same
155*> order as their eigenvalues.
156*> Each eigenvector is scaled so the largest component has
157*> abs(real part) + abs(imag. part) = 1.
158*> Not referenced if JOBVR = 'N'.
159*> \endverbatim
160*>
161*> \param[in] LDVR
162*> \verbatim
163*> LDVR is INTEGER
164*> The leading dimension of the matrix VR. LDVR >= 1, and
165*> if JOBVR = 'V', LDVR >= N.
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*> For good performance, LWORK must generally be larger.
179*>
180*> If LWORK = -1, then a workspace query is assumed; the routine
181*> only calculates the optimal size of the WORK array, returns
182*> this value as the first entry of the WORK array, and no error
183*> message related to LWORK is issued by XERBLA.
184*> \endverbatim
185*>
186*> \param[out] RWORK
187*> \verbatim
188*> RWORK is REAL array, dimension (8*N)
189*> \endverbatim
190*>
191*> \param[out] INFO
192*> \verbatim
193*> INFO is INTEGER
194*> = 0: successful exit
195*> < 0: if INFO = -i, the i-th argument had an illegal value.
196*> =1,...,N:
197*> The QZ iteration failed. No eigenvectors have been
198*> calculated, but ALPHA(j) and BETA(j) should be
199*> correct for j=INFO+1,...,N.
200*> > N: =N+1: other then QZ iteration failed in CHGEQZ,
201*> =N+2: error return from CTGEVC.
202*> \endverbatim
203*
204* Authors:
205* ========
206*
207*> \author Univ. of Tennessee
208*> \author Univ. of California Berkeley
209*> \author Univ. of Colorado Denver
210*> \author NAG Ltd.
211*
212*> \ingroup complexGEeigen
213*
214* =====================================================================
215 SUBROUTINE cggev( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
216 $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
217*
218* -- LAPACK driver routine --
219* -- LAPACK is a software package provided by Univ. of Tennessee, --
220* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221*
222* .. Scalar Arguments ..
223 CHARACTER JOBVL, JOBVR
224 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
225* ..
226* .. Array Arguments ..
227 REAL RWORK( * )
228 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
229 $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
230 $ work( * )
231* ..
232*
233* =====================================================================
234*
235* .. Parameters ..
236 REAL ZERO, ONE
237 parameter( zero = 0.0e0, one = 1.0e0 )
238 COMPLEX CZERO, CONE
239 parameter( czero = ( 0.0e0, 0.0e0 ),
240 $ cone = ( 1.0e0, 0.0e0 ) )
241* ..
242* .. Local Scalars ..
243 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
244 CHARACTER CHTEMP
245 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
246 $ in, iright, irows, irwrk, itau, iwrk, jc, jr,
247 $ lwkmin, lwkopt
248 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
249 $ smlnum, temp
250 COMPLEX X
251* ..
252* .. Local Arrays ..
253 LOGICAL LDUMMA( 1 )
254* ..
255* .. External Subroutines ..
256 EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
258 $ xerbla
259* ..
260* .. External Functions ..
261 LOGICAL LSAME
262 INTEGER ILAENV
263 REAL CLANGE, SLAMCH
264 EXTERNAL lsame, ilaenv, clange, slamch
265* ..
266* .. Intrinsic Functions ..
267 INTRINSIC abs, aimag, max, real, sqrt
268* ..
269* .. Statement Functions ..
270 REAL ABS1
271* ..
272* .. Statement Function definitions ..
273 abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
274* ..
275* .. Executable Statements ..
276*
277* Decode the input arguments
278*
279 IF( lsame( jobvl, 'N' ) ) THEN
280 ijobvl = 1
281 ilvl = .false.
282 ELSE IF( lsame( jobvl, 'V' ) ) THEN
283 ijobvl = 2
284 ilvl = .true.
285 ELSE
286 ijobvl = -1
287 ilvl = .false.
288 END IF
289*
290 IF( lsame( jobvr, 'N' ) ) THEN
291 ijobvr = 1
292 ilvr = .false.
293 ELSE IF( lsame( jobvr, 'V' ) ) THEN
294 ijobvr = 2
295 ilvr = .true.
296 ELSE
297 ijobvr = -1
298 ilvr = .false.
299 END IF
300 ilv = ilvl .OR. ilvr
301*
302* Test the input arguments
303*
304 info = 0
305 lquery = ( lwork.EQ.-1 )
306 IF( ijobvl.LE.0 ) THEN
307 info = -1
308 ELSE IF( ijobvr.LE.0 ) THEN
309 info = -2
310 ELSE IF( n.LT.0 ) THEN
311 info = -3
312 ELSE IF( lda.LT.max( 1, n ) ) THEN
313 info = -5
314 ELSE IF( ldb.LT.max( 1, n ) ) THEN
315 info = -7
316 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
317 info = -11
318 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
319 info = -13
320 END IF
321*
322* Compute workspace
323* (Note: Comments in the code beginning "Workspace:" describe the
324* minimal amount of workspace needed at that point in the code,
325* as well as the preferred amount for good performance.
326* NB refers to the optimal block size for the immediately
327* following subroutine, as returned by ILAENV. The workspace is
328* computed assuming ILO = 1 and IHI = N, the worst case.)
329*
330 IF( info.EQ.0 ) THEN
331 lwkmin = max( 1, 2*n )
332 lwkopt = max( 1, n + n*ilaenv( 1, 'CGEQRF', ' ', n, 1, n, 0 ) )
333 lwkopt = max( lwkopt, n +
334 $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
335 IF( ilvl ) THEN
336 lwkopt = max( lwkopt, n +
337 $ n*ilaenv( 1, 'CUNGQR', ' ', n, 1, n, -1 ) )
338 END IF
339 work( 1 ) = lwkopt
340*
341 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
342 $ info = -15
343 END IF
344*
345 IF( info.NE.0 ) THEN
346 CALL xerbla( 'CGGEV ', -info )
347 RETURN
348 ELSE IF( lquery ) THEN
349 RETURN
350 END IF
351*
352* Quick return if possible
353*
354 IF( n.EQ.0 )
355 $ RETURN
356*
357* Get machine constants
358*
359 eps = slamch( 'E' )*slamch( 'B' )
360 smlnum = slamch( 'S' )
361 bignum = one / smlnum
362 CALL slabad( smlnum, bignum )
363 smlnum = sqrt( smlnum ) / eps
364 bignum = one / smlnum
365*
366* Scale A if max element outside range [SMLNUM,BIGNUM]
367*
368 anrm = clange( 'M', n, n, a, lda, rwork )
369 ilascl = .false.
370 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
371 anrmto = smlnum
372 ilascl = .true.
373 ELSE IF( anrm.GT.bignum ) THEN
374 anrmto = bignum
375 ilascl = .true.
376 END IF
377 IF( ilascl )
378 $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
379*
380* Scale B if max element outside range [SMLNUM,BIGNUM]
381*
382 bnrm = clange( 'M', n, n, b, ldb, rwork )
383 ilbscl = .false.
384 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
385 bnrmto = smlnum
386 ilbscl = .true.
387 ELSE IF( bnrm.GT.bignum ) THEN
388 bnrmto = bignum
389 ilbscl = .true.
390 END IF
391 IF( ilbscl )
392 $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
393*
394* Permute the matrices A, B to isolate eigenvalues if possible
395* (Real Workspace: need 6*N)
396*
397 ileft = 1
398 iright = n + 1
399 irwrk = iright + n
400 CALL cggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
401 $ rwork( iright ), rwork( irwrk ), ierr )
402*
403* Reduce B to triangular form (QR decomposition of B)
404* (Complex Workspace: need N, prefer N*NB)
405*
406 irows = ihi + 1 - ilo
407 IF( ilv ) THEN
408 icols = n + 1 - ilo
409 ELSE
410 icols = irows
411 END IF
412 itau = 1
413 iwrk = itau + irows
414 CALL cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
415 $ work( iwrk ), lwork+1-iwrk, ierr )
416*
417* Apply the orthogonal transformation to matrix A
418* (Complex Workspace: need N, prefer N*NB)
419*
420 CALL cunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
421 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
422 $ lwork+1-iwrk, ierr )
423*
424* Initialize VL
425* (Complex Workspace: need N, prefer N*NB)
426*
427 IF( ilvl ) THEN
428 CALL claset( 'Full', n, n, czero, cone, vl, ldvl )
429 IF( irows.GT.1 ) THEN
430 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
431 $ vl( ilo+1, ilo ), ldvl )
432 END IF
433 CALL cungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
434 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
435 END IF
436*
437* Initialize VR
438*
439 IF( ilvr )
440 $ CALL claset( 'Full', n, n, czero, cone, vr, ldvr )
441*
442* Reduce to generalized Hessenberg form
443*
444 IF( ilv ) THEN
445*
446* Eigenvectors requested -- work on whole matrix.
447*
448 CALL cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
449 $ ldvl, vr, ldvr, ierr )
450 ELSE
451 CALL cgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
452 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
453 END IF
454*
455* Perform QZ algorithm (Compute eigenvalues, and optionally, the
456* Schur form and Schur vectors)
457* (Complex Workspace: need N)
458* (Real Workspace: need N)
459*
460 iwrk = itau
461 IF( ilv ) THEN
462 chtemp = 'S'
463 ELSE
464 chtemp = 'E'
465 END IF
466 CALL chgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
467 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
468 $ lwork+1-iwrk, rwork( irwrk ), ierr )
469 IF( ierr.NE.0 ) THEN
470 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
471 info = ierr
472 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
473 info = ierr - n
474 ELSE
475 info = n + 1
476 END IF
477 GO TO 70
478 END IF
479*
480* Compute Eigenvectors
481* (Real Workspace: need 2*N)
482* (Complex Workspace: need 2*N)
483*
484 IF( ilv ) THEN
485 IF( ilvl ) THEN
486 IF( ilvr ) THEN
487 chtemp = 'B'
488 ELSE
489 chtemp = 'L'
490 END IF
491 ELSE
492 chtemp = 'R'
493 END IF
494*
495 CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
496 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
497 $ ierr )
498 IF( ierr.NE.0 ) THEN
499 info = n + 2
500 GO TO 70
501 END IF
502*
503* Undo balancing on VL and VR and normalization
504* (Workspace: none needed)
505*
506 IF( ilvl ) THEN
507 CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
508 $ rwork( iright ), n, vl, ldvl, ierr )
509 DO 30 jc = 1, n
510 temp = zero
511 DO 10 jr = 1, n
512 temp = max( temp, abs1( vl( jr, jc ) ) )
513 10 CONTINUE
514 IF( temp.LT.smlnum )
515 $ GO TO 30
516 temp = one / temp
517 DO 20 jr = 1, n
518 vl( jr, jc ) = vl( jr, jc )*temp
519 20 CONTINUE
520 30 CONTINUE
521 END IF
522 IF( ilvr ) THEN
523 CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
524 $ rwork( iright ), n, vr, ldvr, ierr )
525 DO 60 jc = 1, n
526 temp = zero
527 DO 40 jr = 1, n
528 temp = max( temp, abs1( vr( jr, jc ) ) )
529 40 CONTINUE
530 IF( temp.LT.smlnum )
531 $ GO TO 60
532 temp = one / temp
533 DO 50 jr = 1, n
534 vr( jr, jc ) = vr( jr, jc )*temp
535 50 CONTINUE
536 60 CONTINUE
537 END IF
538 END IF
539*
540* Undo scaling if necessary
541*
542 70 CONTINUE
543*
544 IF( ilascl )
545 $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
546*
547 IF( ilbscl )
548 $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
549*
550 work( 1 ) = lwkopt
551 RETURN
552*
553* End of CGGEV
554*
555 END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
CGGBAL
Definition: cggbal.f:177
subroutine cggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
CGGBAK
Definition: cggbak.f:148
subroutine chgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
CHGEQZ
Definition: chgeqz.f:284
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:146
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC
Definition: ctgevc.f:219
subroutine cggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: cggev.f:217
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
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:143
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:168
subroutine cgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
CGGHRD
Definition: cgghrd.f:204
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:128