LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ cggev3()

subroutine cggev3 ( 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 
)

CGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)

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

Purpose:
 CGGEV3 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.

          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
January 2015

Definition at line 218 of file cggev3.f.

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