LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cggev()

subroutine cggev ( character  jobvl,
character  jobvr,
integer  n,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( ldb, * )  b,
integer  ldb,
complex, dimension( * )  alpha,
complex, dimension( * )  beta,
complex, dimension( ldvl, * )  vl,
integer  ldvl,
complex, dimension( ldvr, * )  vr,
integer  ldvr,
complex, dimension( * )  work,
integer  lwork,
real, dimension( * )  rwork,
integer  info 
)

CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

Download CGGEV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CGGEV computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, and optionally, the left and/or
 right generalized eigenvectors.

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar
 lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
 singular. It is usually represented as the pair (alpha,beta), as
 there is a reasonable interpretation for beta=0, and even for both
 being zero.

 The right generalized eigenvector v(j) corresponding to the
 generalized eigenvalue lambda(j) of (A,B) satisfies

              A * v(j) = lambda(j) * B * v(j).

 The left generalized eigenvector u(j) corresponding to the
 generalized eigenvalues lambda(j) of (A,B) satisfies

              u(j)**H * A = lambda(j) * u(j)**H * B

 where u(j)**H is the conjugate-transpose of u(j).
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors.
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is COMPLEX array, dimension (LDA, N)
          On entry, the matrix A in the pair (A,B).
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX array, dimension (LDB, N)
          On entry, the matrix B in the pair (A,B).
          On exit, B has been overwritten.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX array, dimension (N)
[out]BETA
          BETA is COMPLEX array, dimension (N)
          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VL
          VL is COMPLEX array, dimension (LDVL,N)
          If JOBVL = 'V', the left generalized eigenvectors u(j) are
          stored one after another in the columns of VL, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX array, dimension (LDVR,N)
          If JOBVR = 'V', the right generalized eigenvectors v(j) are
          stored one after another in the columns of VR, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[out]WORK
          WORK is COMPLEX array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,2*N).
          For good performance, LWORK must generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]RWORK
          RWORK is REAL array, dimension (8*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          =1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHA(j) and BETA(j) should be
                correct for j=INFO+1,...,N.
          > N:  =N+1: other then QZ iteration failed in CHGEQZ,
                =N+2: error return from CTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 215 of file cggev.f.

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* ..
259* .. External Functions ..
260 LOGICAL LSAME
261 INTEGER ILAENV
262 REAL CLANGE, SLAMCH, 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, 0 ) )
332 lwkopt = max( lwkopt, n +
333 $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
334 IF( ilvl ) THEN
335 lwkopt = max( lwkopt, n +
336 $ n*ilaenv( 1, 'CUNGQR', ' ', n, 1, n, -1 ) )
337 END IF
338 work( 1 ) = sroundup_lwork(lwkopt)
339*
340 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
341 $ info = -15
342 END IF
343*
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'CGGEV ', -info )
346 RETURN
347 ELSE IF( lquery ) THEN
348 RETURN
349 END IF
350*
351* Quick return if possible
352*
353 IF( n.EQ.0 )
354 $ RETURN
355*
356* Get machine constants
357*
358 eps = slamch( 'E' )*slamch( 'B' )
359 smlnum = slamch( 'S' )
360 bignum = one / smlnum
361 smlnum = sqrt( smlnum ) / eps
362 bignum = one / smlnum
363*
364* Scale A if max element outside range [SMLNUM,BIGNUM]
365*
366 anrm = clange( 'M', n, n, a, lda, rwork )
367 ilascl = .false.
368 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
369 anrmto = smlnum
370 ilascl = .true.
371 ELSE IF( anrm.GT.bignum ) THEN
372 anrmto = bignum
373 ilascl = .true.
374 END IF
375 IF( ilascl )
376 $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
377*
378* Scale B if max element outside range [SMLNUM,BIGNUM]
379*
380 bnrm = clange( 'M', n, n, b, ldb, rwork )
381 ilbscl = .false.
382 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
383 bnrmto = smlnum
384 ilbscl = .true.
385 ELSE IF( bnrm.GT.bignum ) THEN
386 bnrmto = bignum
387 ilbscl = .true.
388 END IF
389 IF( ilbscl )
390 $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
391*
392* Permute the matrices A, B to isolate eigenvalues if possible
393* (Real Workspace: need 6*N)
394*
395 ileft = 1
396 iright = n + 1
397 irwrk = iright + n
398 CALL cggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
399 $ rwork( iright ), rwork( irwrk ), ierr )
400*
401* Reduce B to triangular form (QR decomposition of B)
402* (Complex Workspace: need N, prefer N*NB)
403*
404 irows = ihi + 1 - ilo
405 IF( ilv ) THEN
406 icols = n + 1 - ilo
407 ELSE
408 icols = irows
409 END IF
410 itau = 1
411 iwrk = itau + irows
412 CALL cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
413 $ work( iwrk ), lwork+1-iwrk, ierr )
414*
415* Apply the orthogonal transformation to matrix A
416* (Complex Workspace: need N, prefer N*NB)
417*
418 CALL cunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
419 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
420 $ lwork+1-iwrk, ierr )
421*
422* Initialize VL
423* (Complex Workspace: need N, prefer N*NB)
424*
425 IF( ilvl ) THEN
426 CALL claset( 'Full', n, n, czero, cone, vl, ldvl )
427 IF( irows.GT.1 ) THEN
428 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
429 $ vl( ilo+1, ilo ), ldvl )
430 END IF
431 CALL cungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
432 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
433 END IF
434*
435* Initialize VR
436*
437 IF( ilvr )
438 $ CALL claset( 'Full', n, n, czero, cone, vr, ldvr )
439*
440* Reduce to generalized Hessenberg form
441*
442 IF( ilv ) THEN
443*
444* Eigenvectors requested -- work on whole matrix.
445*
446 CALL cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
447 $ ldvl, vr, ldvr, ierr )
448 ELSE
449 CALL cgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
450 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
451 END IF
452*
453* Perform QZ algorithm (Compute eigenvalues, and optionally, the
454* Schur form and Schur vectors)
455* (Complex Workspace: need N)
456* (Real Workspace: need N)
457*
458 iwrk = itau
459 IF( ilv ) THEN
460 chtemp = 'S'
461 ELSE
462 chtemp = 'E'
463 END IF
464 CALL chgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
465 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
466 $ lwork+1-iwrk, rwork( irwrk ), ierr )
467 IF( ierr.NE.0 ) THEN
468 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
469 info = ierr
470 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
471 info = ierr - n
472 ELSE
473 info = n + 1
474 END IF
475 GO TO 70
476 END IF
477*
478* Compute Eigenvectors
479* (Real Workspace: need 2*N)
480* (Complex Workspace: need 2*N)
481*
482 IF( ilv ) THEN
483 IF( ilvl ) THEN
484 IF( ilvr ) THEN
485 chtemp = 'B'
486 ELSE
487 chtemp = 'L'
488 END IF
489 ELSE
490 chtemp = 'R'
491 END IF
492*
493 CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
494 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
495 $ ierr )
496 IF( ierr.NE.0 ) THEN
497 info = n + 2
498 GO TO 70
499 END IF
500*
501* Undo balancing on VL and VR and normalization
502* (Workspace: none needed)
503*
504 IF( ilvl ) THEN
505 CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
506 $ rwork( iright ), n, vl, ldvl, ierr )
507 DO 30 jc = 1, n
508 temp = zero
509 DO 10 jr = 1, n
510 temp = max( temp, abs1( vl( jr, jc ) ) )
511 10 CONTINUE
512 IF( temp.LT.smlnum )
513 $ GO TO 30
514 temp = one / temp
515 DO 20 jr = 1, n
516 vl( jr, jc ) = vl( jr, jc )*temp
517 20 CONTINUE
518 30 CONTINUE
519 END IF
520 IF( ilvr ) THEN
521 CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
522 $ rwork( iright ), n, vr, ldvr, ierr )
523 DO 60 jc = 1, n
524 temp = zero
525 DO 40 jr = 1, n
526 temp = max( temp, abs1( vr( jr, jc ) ) )
527 40 CONTINUE
528 IF( temp.LT.smlnum )
529 $ GO TO 60
530 temp = one / temp
531 DO 50 jr = 1, n
532 vr( jr, jc ) = vr( jr, jc )*temp
533 50 CONTINUE
534 60 CONTINUE
535 END IF
536 END IF
537*
538* Undo scaling if necessary
539*
540 70 CONTINUE
541*
542 IF( ilascl )
543 $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
544*
545 IF( ilbscl )
546 $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
547*
548 work( 1 ) = sroundup_lwork(lwkopt)
549 RETURN
550*
551* End of CGGEV
552*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
subroutine cggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
CGGBAK
Definition cggbak.f:148
subroutine cggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
CGGBAL
Definition cggbal.f:177
subroutine cgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
CGGHRD
Definition cgghrd.f:204
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
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 cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
Here is the call graph for this function:
Here is the caller graph for this function: