LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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 SHGEQZ,
                =N+2: error return from STGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 219 of file cggev.f.

219 *
220 * -- LAPACK driver routine (version 3.7.0) --
221 * -- LAPACK is a software package provided by Univ. of Tennessee, --
222 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223 * April 2012
224 *
225 * .. Scalar Arguments ..
226  CHARACTER jobvl, jobvr
227  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
228 * ..
229 * .. Array Arguments ..
230  REAL rwork( * )
231  COMPLEX a( lda, * ), alpha( * ), b( ldb, * ),
232  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
233  $ work( * )
234 * ..
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239  REAL zero, one
240  parameter( zero = 0.0e0, one = 1.0e0 )
241  COMPLEX czero, cone
242  parameter( czero = ( 0.0e0, 0.0e0 ),
243  $ cone = ( 1.0e0, 0.0e0 ) )
244 * ..
245 * .. Local Scalars ..
246  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery
247  CHARACTER chtemp
248  INTEGER icols, ierr, ihi, ijobvl, ijobvr, ileft, ilo,
249  $ in, iright, irows, irwrk, itau, iwrk, jc, jr,
250  $ lwkmin, lwkopt
251  REAL anrm, anrmto, bignum, bnrm, bnrmto, eps,
252  $ smlnum, temp
253  COMPLEX x
254 * ..
255 * .. Local Arrays ..
256  LOGICAL ldumma( 1 )
257 * ..
258 * .. External Subroutines ..
259  EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
261  $ xerbla
262 * ..
263 * .. External Functions ..
264  LOGICAL lsame
265  INTEGER ilaenv
266  REAL clange, slamch
267  EXTERNAL lsame, ilaenv, clange, slamch
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC abs, aimag, max, REAL, sqrt
271 * ..
272 * .. Statement Functions ..
273  REAL abs1
274 * ..
275 * .. Statement Function definitions ..
276  abs1( x ) = abs( REAL( X ) ) + abs( aimag( x ) )
277 * ..
278 * .. Executable Statements ..
279 *
280 * Decode the input arguments
281 *
282  IF( lsame( jobvl, 'N' ) ) THEN
283  ijobvl = 1
284  ilvl = .false.
285  ELSE IF( lsame( jobvl, 'V' ) ) THEN
286  ijobvl = 2
287  ilvl = .true.
288  ELSE
289  ijobvl = -1
290  ilvl = .false.
291  END IF
292 *
293  IF( lsame( jobvr, 'N' ) ) THEN
294  ijobvr = 1
295  ilvr = .false.
296  ELSE IF( lsame( jobvr, 'V' ) ) THEN
297  ijobvr = 2
298  ilvr = .true.
299  ELSE
300  ijobvr = -1
301  ilvr = .false.
302  END IF
303  ilv = ilvl .OR. ilvr
304 *
305 * Test the input arguments
306 *
307  info = 0
308  lquery = ( lwork.EQ.-1 )
309  IF( ijobvl.LE.0 ) THEN
310  info = -1
311  ELSE IF( ijobvr.LE.0 ) THEN
312  info = -2
313  ELSE IF( n.LT.0 ) THEN
314  info = -3
315  ELSE IF( lda.LT.max( 1, n ) ) THEN
316  info = -5
317  ELSE IF( ldb.LT.max( 1, n ) ) THEN
318  info = -7
319  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
320  info = -11
321  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
322  info = -13
323  END IF
324 *
325 * Compute workspace
326 * (Note: Comments in the code beginning "Workspace:" describe the
327 * minimal amount of workspace needed at that point in the code,
328 * as well as the preferred amount for good performance.
329 * NB refers to the optimal block size for the immediately
330 * following subroutine, as returned by ILAENV. The workspace is
331 * computed assuming ILO = 1 and IHI = N, the worst case.)
332 *
333  IF( info.EQ.0 ) THEN
334  lwkmin = max( 1, 2*n )
335  lwkopt = max( 1, n + n*ilaenv( 1, 'CGEQRF', ' ', n, 1, n, 0 ) )
336  lwkopt = max( lwkopt, n +
337  $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
338  IF( ilvl ) THEN
339  lwkopt = max( lwkopt, n +
340  $ n*ilaenv( 1, 'CUNGQR', ' ', n, 1, n, -1 ) )
341  END IF
342  work( 1 ) = lwkopt
343 *
344  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
345  $ info = -15
346  END IF
347 *
348  IF( info.NE.0 ) THEN
349  CALL xerbla( 'CGGEV ', -info )
350  RETURN
351  ELSE IF( lquery ) THEN
352  RETURN
353  END IF
354 *
355 * Quick return if possible
356 *
357  IF( n.EQ.0 )
358  $ RETURN
359 *
360 * Get machine constants
361 *
362  eps = slamch( 'E' )*slamch( 'B' )
363  smlnum = slamch( 'S' )
364  bignum = one / smlnum
365  CALL slabad( smlnum, bignum )
366  smlnum = sqrt( smlnum ) / eps
367  bignum = one / smlnum
368 *
369 * Scale A if max element outside range [SMLNUM,BIGNUM]
370 *
371  anrm = clange( 'M', n, n, a, lda, rwork )
372  ilascl = .false.
373  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
374  anrmto = smlnum
375  ilascl = .true.
376  ELSE IF( anrm.GT.bignum ) THEN
377  anrmto = bignum
378  ilascl = .true.
379  END IF
380  IF( ilascl )
381  $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
382 *
383 * Scale B if max element outside range [SMLNUM,BIGNUM]
384 *
385  bnrm = clange( 'M', n, n, b, ldb, rwork )
386  ilbscl = .false.
387  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
388  bnrmto = smlnum
389  ilbscl = .true.
390  ELSE IF( bnrm.GT.bignum ) THEN
391  bnrmto = bignum
392  ilbscl = .true.
393  END IF
394  IF( ilbscl )
395  $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
396 *
397 * Permute the matrices A, B to isolate eigenvalues if possible
398 * (Real Workspace: need 6*N)
399 *
400  ileft = 1
401  iright = n + 1
402  irwrk = iright + n
403  CALL cggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
404  $ rwork( iright ), rwork( irwrk ), ierr )
405 *
406 * Reduce B to triangular form (QR decomposition of B)
407 * (Complex Workspace: need N, prefer N*NB)
408 *
409  irows = ihi + 1 - ilo
410  IF( ilv ) THEN
411  icols = n + 1 - ilo
412  ELSE
413  icols = irows
414  END IF
415  itau = 1
416  iwrk = itau + irows
417  CALL cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
418  $ work( iwrk ), lwork+1-iwrk, ierr )
419 *
420 * Apply the orthogonal transformation to matrix A
421 * (Complex Workspace: need N, prefer N*NB)
422 *
423  CALL cunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
424  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
425  $ lwork+1-iwrk, ierr )
426 *
427 * Initialize VL
428 * (Complex Workspace: need N, prefer N*NB)
429 *
430  IF( ilvl ) THEN
431  CALL claset( 'Full', n, n, czero, cone, vl, ldvl )
432  IF( irows.GT.1 ) THEN
433  CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
434  $ vl( ilo+1, ilo ), ldvl )
435  END IF
436  CALL cungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
437  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
438  END IF
439 *
440 * Initialize VR
441 *
442  IF( ilvr )
443  $ CALL claset( 'Full', n, n, czero, cone, vr, ldvr )
444 *
445 * Reduce to generalized Hessenberg form
446 *
447  IF( ilv ) THEN
448 *
449 * Eigenvectors requested -- work on whole matrix.
450 *
451  CALL cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
452  $ ldvl, vr, ldvr, ierr )
453  ELSE
454  CALL cgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
455  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
456  END IF
457 *
458 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
459 * Schur form and Schur vectors)
460 * (Complex Workspace: need N)
461 * (Real Workspace: need N)
462 *
463  iwrk = itau
464  IF( ilv ) THEN
465  chtemp = 'S'
466  ELSE
467  chtemp = 'E'
468  END IF
469  CALL chgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
470  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
471  $ lwork+1-iwrk, rwork( irwrk ), ierr )
472  IF( ierr.NE.0 ) THEN
473  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
474  info = ierr
475  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
476  info = ierr - n
477  ELSE
478  info = n + 1
479  END IF
480  GO TO 70
481  END IF
482 *
483 * Compute Eigenvectors
484 * (Real Workspace: need 2*N)
485 * (Complex Workspace: need 2*N)
486 *
487  IF( ilv ) THEN
488  IF( ilvl ) THEN
489  IF( ilvr ) THEN
490  chtemp = 'B'
491  ELSE
492  chtemp = 'L'
493  END IF
494  ELSE
495  chtemp = 'R'
496  END IF
497 *
498  CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
499  $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
500  $ ierr )
501  IF( ierr.NE.0 ) THEN
502  info = n + 2
503  GO TO 70
504  END IF
505 *
506 * Undo balancing on VL and VR and normalization
507 * (Workspace: none needed)
508 *
509  IF( ilvl ) THEN
510  CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
511  $ rwork( iright ), n, vl, ldvl, ierr )
512  DO 30 jc = 1, n
513  temp = zero
514  DO 10 jr = 1, n
515  temp = max( temp, abs1( vl( jr, jc ) ) )
516  10 CONTINUE
517  IF( temp.LT.smlnum )
518  $ GO TO 30
519  temp = one / temp
520  DO 20 jr = 1, n
521  vl( jr, jc ) = vl( jr, jc )*temp
522  20 CONTINUE
523  30 CONTINUE
524  END IF
525  IF( ilvr ) THEN
526  CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
527  $ rwork( iright ), n, vr, ldvr, ierr )
528  DO 60 jc = 1, n
529  temp = zero
530  DO 40 jr = 1, n
531  temp = max( temp, abs1( vr( jr, jc ) ) )
532  40 CONTINUE
533  IF( temp.LT.smlnum )
534  $ GO TO 60
535  temp = one / temp
536  DO 50 jr = 1, n
537  vr( jr, jc ) = vr( jr, jc )*temp
538  50 CONTINUE
539  60 CONTINUE
540  END IF
541  END IF
542 *
543 * Undo scaling if necessary
544 *
545  70 CONTINUE
546 *
547  IF( ilascl )
548  $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
549 *
550  IF( ilbscl )
551  $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
552 *
553  work( 1 ) = lwkopt
554  RETURN
555 *
556 * End of CGGEV
557 *
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC
Definition: ctgevc.f:221
subroutine cggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
CGGBAL
Definition: cggbal.f:179
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:108
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:130
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:145
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:286
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
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:117
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:138
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine cgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
CGGHRD
Definition: cgghrd.f:206
subroutine cggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
CGGBAK
Definition: cggbak.f:150
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
Here is the call graph for this function:
Here is the caller graph for this function: