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