LAPACK 3.11.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 $ 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*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
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 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 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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: