LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dggev()

subroutine dggev ( character  JOBVL,
character  JOBVR,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( * )  ALPHAR,
double precision, dimension( * )  ALPHAI,
double precision, dimension( * )  BETA,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

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

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

Purpose:
 DGGEV computes for a pair of N-by-N real 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 eigenvector v(j) corresponding to the eigenvalue lambda(j)
 of (A,B) satisfies

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

 The left eigenvector u(j) corresponding to the eigenvalue 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (N)
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
          the j-th eigenvalue is real; if positive, then the j-th and
          (j+1)-st eigenvalues are a complex conjugate pair, with
          ALPHAI(j+1) negative.

          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(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, ALPHAR and ALPHAI 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 DOUBLE PRECISION array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored one
          after another in the columns of VL, in the same order as
          their eigenvalues. If the j-th eigenvalue is real, then
          u(j) = VL(:,j), the j-th column of VL. If the j-th and
          (j+1)-th eigenvalues form a complex conjugate pair, then
          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
          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 DOUBLE PRECISION array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors v(j) are stored one
          after another in the columns of VR, in the same order as
          their eigenvalues. If the j-th eigenvalue is real, then
          v(j) = VR(:,j), the j-th column of VR. If the j-th and
          (j+1)-th eigenvalues form a complex conjugate pair, then
          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
          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 DOUBLE PRECISION 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,8*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]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 ALPHAR(j), ALPHAI(j), and BETA(j)
                should be correct for j=INFO+1,...,N.
          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
                =N+2: error return from DTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 228 of file dggev.f.

228 *
229 * -- LAPACK driver routine (version 3.7.0) --
230 * -- LAPACK is a software package provided by Univ. of Tennessee, --
231 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232 * April 2012
233 *
234 * .. Scalar Arguments ..
235  CHARACTER jobvl, jobvr
236  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
237 * ..
238 * .. Array Arguments ..
239  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
240  $ b( ldb, * ), beta( * ), vl( ldvl, * ),
241  $ vr( ldvr, * ), work( * )
242 * ..
243 *
244 * =====================================================================
245 *
246 * .. Parameters ..
247  DOUBLE PRECISION zero, one
248  parameter( zero = 0.0d+0, one = 1.0d+0 )
249 * ..
250 * .. Local Scalars ..
251  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery
252  CHARACTER chtemp
253  INTEGER icols, ierr, ihi, ijobvl, ijobvr, ileft, ilo,
254  $ in, iright, irows, itau, iwrk, jc, jr, maxwrk,
255  $ minwrk
256  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
257  $ smlnum, temp
258 * ..
259 * .. Local Arrays ..
260  LOGICAL ldumma( 1 )
261 * ..
262 * .. External Subroutines ..
263  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
265  $ xerbla
266 * ..
267 * .. External Functions ..
268  LOGICAL lsame
269  INTEGER ilaenv
270  DOUBLE PRECISION dlamch, dlange
271  EXTERNAL lsame, ilaenv, dlamch, dlange
272 * ..
273 * .. Intrinsic Functions ..
274  INTRINSIC abs, max, sqrt
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 = -12
319  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
320  info = -14
321  END IF
322 *
323 * Compute workspace
324 * (Note: Comments in the code beginning "Workspace:" describe the
325 * minimal amount of workspace needed at that point in the code,
326 * as well as the preferred amount for good performance.
327 * NB refers to the optimal block size for the immediately
328 * following subroutine, as returned by ILAENV. The workspace is
329 * computed assuming ILO = 1 and IHI = N, the worst case.)
330 *
331  IF( info.EQ.0 ) THEN
332  minwrk = max( 1, 8*n )
333  maxwrk = max( 1, n*( 7 +
334  $ ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 ) ) )
335  maxwrk = max( maxwrk, n*( 7 +
336  $ ilaenv( 1, 'DORMQR', ' ', n, 1, n, 0 ) ) )
337  IF( ilvl ) THEN
338  maxwrk = max( maxwrk, n*( 7 +
339  $ ilaenv( 1, 'DORGQR', ' ', n, 1, n, -1 ) ) )
340  END IF
341  work( 1 ) = maxwrk
342 *
343  IF( lwork.LT.minwrk .AND. .NOT.lquery )
344  $ info = -16
345  END IF
346 *
347  IF( info.NE.0 ) THEN
348  CALL xerbla( 'DGGEV ', -info )
349  RETURN
350  ELSE IF( lquery ) THEN
351  RETURN
352  END IF
353 *
354 * Quick return if possible
355 *
356  IF( n.EQ.0 )
357  $ RETURN
358 *
359 * Get machine constants
360 *
361  eps = dlamch( 'P' )
362  smlnum = dlamch( 'S' )
363  bignum = one / smlnum
364  CALL dlabad( smlnum, bignum )
365  smlnum = sqrt( smlnum ) / eps
366  bignum = one / smlnum
367 *
368 * Scale A if max element outside range [SMLNUM,BIGNUM]
369 *
370  anrm = dlange( 'M', n, n, a, lda, work )
371  ilascl = .false.
372  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
373  anrmto = smlnum
374  ilascl = .true.
375  ELSE IF( anrm.GT.bignum ) THEN
376  anrmto = bignum
377  ilascl = .true.
378  END IF
379  IF( ilascl )
380  $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
381 *
382 * Scale B if max element outside range [SMLNUM,BIGNUM]
383 *
384  bnrm = dlange( 'M', n, n, b, ldb, work )
385  ilbscl = .false.
386  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
387  bnrmto = smlnum
388  ilbscl = .true.
389  ELSE IF( bnrm.GT.bignum ) THEN
390  bnrmto = bignum
391  ilbscl = .true.
392  END IF
393  IF( ilbscl )
394  $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
395 *
396 * Permute the matrices A, B to isolate eigenvalues if possible
397 * (Workspace: need 6*N)
398 *
399  ileft = 1
400  iright = n + 1
401  iwrk = iright + n
402  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
403  $ work( iright ), work( iwrk ), ierr )
404 *
405 * Reduce B to triangular form (QR decomposition of B)
406 * (Workspace: need N, prefer N*NB)
407 *
408  irows = ihi + 1 - ilo
409  IF( ilv ) THEN
410  icols = n + 1 - ilo
411  ELSE
412  icols = irows
413  END IF
414  itau = iwrk
415  iwrk = itau + irows
416  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
417  $ work( iwrk ), lwork+1-iwrk, ierr )
418 *
419 * Apply the orthogonal transformation to matrix A
420 * (Workspace: need N, prefer N*NB)
421 *
422  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
423  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
424  $ lwork+1-iwrk, ierr )
425 *
426 * Initialize VL
427 * (Workspace: need N, prefer N*NB)
428 *
429  IF( ilvl ) THEN
430  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
431  IF( irows.GT.1 ) THEN
432  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
433  $ vl( ilo+1, ilo ), ldvl )
434  END IF
435  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
436  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
437  END IF
438 *
439 * Initialize VR
440 *
441  IF( ilvr )
442  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
443 *
444 * Reduce to generalized Hessenberg form
445 * (Workspace: none needed)
446 *
447  IF( ilv ) THEN
448 *
449 * Eigenvectors requested -- work on whole matrix.
450 *
451  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
452  $ ldvl, vr, ldvr, ierr )
453  ELSE
454  CALL dgghrd( '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 forms and Schur vectors)
460 * (Workspace: need N)
461 *
462  iwrk = itau
463  IF( ilv ) THEN
464  chtemp = 'S'
465  ELSE
466  chtemp = 'E'
467  END IF
468  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
469  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
470  $ work( iwrk ), lwork+1-iwrk, ierr )
471  IF( ierr.NE.0 ) THEN
472  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
473  info = ierr
474  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
475  info = ierr - n
476  ELSE
477  info = n + 1
478  END IF
479  GO TO 110
480  END IF
481 *
482 * Compute Eigenvectors
483 * (Workspace: need 6*N)
484 *
485  IF( ilv ) THEN
486  IF( ilvl ) THEN
487  IF( ilvr ) THEN
488  chtemp = 'B'
489  ELSE
490  chtemp = 'L'
491  END IF
492  ELSE
493  chtemp = 'R'
494  END IF
495  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
496  $ vr, ldvr, n, in, work( iwrk ), ierr )
497  IF( ierr.NE.0 ) THEN
498  info = n + 2
499  GO TO 110
500  END IF
501 *
502 * Undo balancing on VL and VR and normalization
503 * (Workspace: none needed)
504 *
505  IF( ilvl ) THEN
506  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
507  $ work( iright ), n, vl, ldvl, ierr )
508  DO 50 jc = 1, n
509  IF( alphai( jc ).LT.zero )
510  $ GO TO 50
511  temp = zero
512  IF( alphai( jc ).EQ.zero ) THEN
513  DO 10 jr = 1, n
514  temp = max( temp, abs( vl( jr, jc ) ) )
515  10 CONTINUE
516  ELSE
517  DO 20 jr = 1, n
518  temp = max( temp, abs( vl( jr, jc ) )+
519  $ abs( vl( jr, jc+1 ) ) )
520  20 CONTINUE
521  END IF
522  IF( temp.LT.smlnum )
523  $ GO TO 50
524  temp = one / temp
525  IF( alphai( jc ).EQ.zero ) THEN
526  DO 30 jr = 1, n
527  vl( jr, jc ) = vl( jr, jc )*temp
528  30 CONTINUE
529  ELSE
530  DO 40 jr = 1, n
531  vl( jr, jc ) = vl( jr, jc )*temp
532  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
533  40 CONTINUE
534  END IF
535  50 CONTINUE
536  END IF
537  IF( ilvr ) THEN
538  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
539  $ work( iright ), n, vr, ldvr, ierr )
540  DO 100 jc = 1, n
541  IF( alphai( jc ).LT.zero )
542  $ GO TO 100
543  temp = zero
544  IF( alphai( jc ).EQ.zero ) THEN
545  DO 60 jr = 1, n
546  temp = max( temp, abs( vr( jr, jc ) ) )
547  60 CONTINUE
548  ELSE
549  DO 70 jr = 1, n
550  temp = max( temp, abs( vr( jr, jc ) )+
551  $ abs( vr( jr, jc+1 ) ) )
552  70 CONTINUE
553  END IF
554  IF( temp.LT.smlnum )
555  $ GO TO 100
556  temp = one / temp
557  IF( alphai( jc ).EQ.zero ) THEN
558  DO 80 jr = 1, n
559  vr( jr, jc ) = vr( jr, jc )*temp
560  80 CONTINUE
561  ELSE
562  DO 90 jr = 1, n
563  vr( jr, jc ) = vr( jr, jc )*temp
564  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
565  90 CONTINUE
566  END IF
567  100 CONTINUE
568  END IF
569 *
570 * End of eigenvector calculation
571 *
572  END IF
573 *
574 * Undo scaling if necessary
575 *
576  110 CONTINUE
577 *
578  IF( ilascl ) THEN
579  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
580  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
581  END IF
582 *
583  IF( ilbscl ) THEN
584  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
585  END IF
586 *
587  work( 1 ) = maxwrk
588  RETURN
589 *
590 * End of DGGEV
591 *
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
Definition: dggbak.f:149
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
Definition: dgghrd.f:209
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
Definition: dggbal.f:179
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:169
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:130
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
Definition: dtgevc.f:297
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
Definition: dhgeqz.f:306
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:138
Here is the call graph for this function:
Here is the caller graph for this function: