LAPACK  3.9.1
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.

Definition at line 224 of file dggev.f.

226 *
227 * -- LAPACK driver routine --
228 * -- LAPACK is a software package provided by Univ. of Tennessee, --
229 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230 *
231 * .. Scalar Arguments ..
232  CHARACTER JOBVL, JOBVR
233  INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
234 * ..
235 * .. Array Arguments ..
236  DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
237  $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
238  $ VR( LDVR, * ), WORK( * )
239 * ..
240 *
241 * =====================================================================
242 *
243 * .. Parameters ..
244  DOUBLE PRECISION ZERO, ONE
245  parameter( zero = 0.0d+0, one = 1.0d+0 )
246 * ..
247 * .. Local Scalars ..
248  LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
249  CHARACTER CHTEMP
250  INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
251  $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
252  $ MINWRK
253  DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
254  $ SMLNUM, TEMP
255 * ..
256 * .. Local Arrays ..
257  LOGICAL LDUMMA( 1 )
258 * ..
259 * .. External Subroutines ..
260  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
262  $ xerbla
263 * ..
264 * .. External Functions ..
265  LOGICAL LSAME
266  INTEGER ILAENV
267  DOUBLE PRECISION DLAMCH, DLANGE
268  EXTERNAL lsame, ilaenv, dlamch, dlange
269 * ..
270 * .. Intrinsic Functions ..
271  INTRINSIC abs, max, sqrt
272 * ..
273 * .. Executable Statements ..
274 *
275 * Decode the input arguments
276 *
277  IF( lsame( jobvl, 'N' ) ) THEN
278  ijobvl = 1
279  ilvl = .false.
280  ELSE IF( lsame( jobvl, 'V' ) ) THEN
281  ijobvl = 2
282  ilvl = .true.
283  ELSE
284  ijobvl = -1
285  ilvl = .false.
286  END IF
287 *
288  IF( lsame( jobvr, 'N' ) ) THEN
289  ijobvr = 1
290  ilvr = .false.
291  ELSE IF( lsame( jobvr, 'V' ) ) THEN
292  ijobvr = 2
293  ilvr = .true.
294  ELSE
295  ijobvr = -1
296  ilvr = .false.
297  END IF
298  ilv = ilvl .OR. ilvr
299 *
300 * Test the input arguments
301 *
302  info = 0
303  lquery = ( lwork.EQ.-1 )
304  IF( ijobvl.LE.0 ) THEN
305  info = -1
306  ELSE IF( ijobvr.LE.0 ) THEN
307  info = -2
308  ELSE IF( n.LT.0 ) THEN
309  info = -3
310  ELSE IF( lda.LT.max( 1, n ) ) THEN
311  info = -5
312  ELSE IF( ldb.LT.max( 1, n ) ) THEN
313  info = -7
314  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
315  info = -12
316  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
317  info = -14
318  END IF
319 *
320 * Compute workspace
321 * (Note: Comments in the code beginning "Workspace:" describe the
322 * minimal amount of workspace needed at that point in the code,
323 * as well as the preferred amount for good performance.
324 * NB refers to the optimal block size for the immediately
325 * following subroutine, as returned by ILAENV. The workspace is
326 * computed assuming ILO = 1 and IHI = N, the worst case.)
327 *
328  IF( info.EQ.0 ) THEN
329  minwrk = max( 1, 8*n )
330  maxwrk = max( 1, n*( 7 +
331  $ ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 ) ) )
332  maxwrk = max( maxwrk, n*( 7 +
333  $ ilaenv( 1, 'DORMQR', ' ', n, 1, n, 0 ) ) )
334  IF( ilvl ) THEN
335  maxwrk = max( maxwrk, n*( 7 +
336  $ ilaenv( 1, 'DORGQR', ' ', n, 1, n, -1 ) ) )
337  END IF
338  work( 1 ) = maxwrk
339 *
340  IF( lwork.LT.minwrk .AND. .NOT.lquery )
341  $ info = -16
342  END IF
343 *
344  IF( info.NE.0 ) THEN
345  CALL xerbla( 'DGGEV ', -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 = dlamch( 'P' )
359  smlnum = dlamch( 'S' )
360  bignum = one / smlnum
361  CALL dlabad( smlnum, bignum )
362  smlnum = sqrt( smlnum ) / eps
363  bignum = one / smlnum
364 *
365 * Scale A if max element outside range [SMLNUM,BIGNUM]
366 *
367  anrm = dlange( 'M', n, n, a, lda, work )
368  ilascl = .false.
369  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
370  anrmto = smlnum
371  ilascl = .true.
372  ELSE IF( anrm.GT.bignum ) THEN
373  anrmto = bignum
374  ilascl = .true.
375  END IF
376  IF( ilascl )
377  $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
378 *
379 * Scale B if max element outside range [SMLNUM,BIGNUM]
380 *
381  bnrm = dlange( 'M', n, n, b, ldb, work )
382  ilbscl = .false.
383  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
384  bnrmto = smlnum
385  ilbscl = .true.
386  ELSE IF( bnrm.GT.bignum ) THEN
387  bnrmto = bignum
388  ilbscl = .true.
389  END IF
390  IF( ilbscl )
391  $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
392 *
393 * Permute the matrices A, B to isolate eigenvalues if possible
394 * (Workspace: need 6*N)
395 *
396  ileft = 1
397  iright = n + 1
398  iwrk = iright + n
399  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
400  $ work( iright ), work( iwrk ), ierr )
401 *
402 * Reduce B to triangular form (QR decomposition of B)
403 * (Workspace: need N, prefer N*NB)
404 *
405  irows = ihi + 1 - ilo
406  IF( ilv ) THEN
407  icols = n + 1 - ilo
408  ELSE
409  icols = irows
410  END IF
411  itau = iwrk
412  iwrk = itau + irows
413  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
414  $ work( iwrk ), lwork+1-iwrk, ierr )
415 *
416 * Apply the orthogonal transformation to matrix A
417 * (Workspace: need N, prefer N*NB)
418 *
419  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
420  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
421  $ lwork+1-iwrk, ierr )
422 *
423 * Initialize VL
424 * (Workspace: need N, prefer N*NB)
425 *
426  IF( ilvl ) THEN
427  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
428  IF( irows.GT.1 ) THEN
429  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
430  $ vl( ilo+1, ilo ), ldvl )
431  END IF
432  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
433  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
434  END IF
435 *
436 * Initialize VR
437 *
438  IF( ilvr )
439  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
440 *
441 * Reduce to generalized Hessenberg form
442 * (Workspace: none needed)
443 *
444  IF( ilv ) THEN
445 *
446 * Eigenvectors requested -- work on whole matrix.
447 *
448  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
449  $ ldvl, vr, ldvr, ierr )
450  ELSE
451  CALL dgghrd( '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 forms and Schur vectors)
457 * (Workspace: need N)
458 *
459  iwrk = itau
460  IF( ilv ) THEN
461  chtemp = 'S'
462  ELSE
463  chtemp = 'E'
464  END IF
465  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
466  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
467  $ work( iwrk ), lwork+1-iwrk, ierr )
468  IF( ierr.NE.0 ) THEN
469  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
470  info = ierr
471  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
472  info = ierr - n
473  ELSE
474  info = n + 1
475  END IF
476  GO TO 110
477  END IF
478 *
479 * Compute Eigenvectors
480 * (Workspace: need 6*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  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
493  $ vr, ldvr, n, in, work( iwrk ), ierr )
494  IF( ierr.NE.0 ) THEN
495  info = n + 2
496  GO TO 110
497  END IF
498 *
499 * Undo balancing on VL and VR and normalization
500 * (Workspace: none needed)
501 *
502  IF( ilvl ) THEN
503  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
504  $ work( iright ), n, vl, ldvl, ierr )
505  DO 50 jc = 1, n
506  IF( alphai( jc ).LT.zero )
507  $ GO TO 50
508  temp = zero
509  IF( alphai( jc ).EQ.zero ) THEN
510  DO 10 jr = 1, n
511  temp = max( temp, abs( vl( jr, jc ) ) )
512  10 CONTINUE
513  ELSE
514  DO 20 jr = 1, n
515  temp = max( temp, abs( vl( jr, jc ) )+
516  $ abs( vl( jr, jc+1 ) ) )
517  20 CONTINUE
518  END IF
519  IF( temp.LT.smlnum )
520  $ GO TO 50
521  temp = one / temp
522  IF( alphai( jc ).EQ.zero ) THEN
523  DO 30 jr = 1, n
524  vl( jr, jc ) = vl( jr, jc )*temp
525  30 CONTINUE
526  ELSE
527  DO 40 jr = 1, n
528  vl( jr, jc ) = vl( jr, jc )*temp
529  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
530  40 CONTINUE
531  END IF
532  50 CONTINUE
533  END IF
534  IF( ilvr ) THEN
535  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
536  $ work( iright ), n, vr, ldvr, ierr )
537  DO 100 jc = 1, n
538  IF( alphai( jc ).LT.zero )
539  $ GO TO 100
540  temp = zero
541  IF( alphai( jc ).EQ.zero ) THEN
542  DO 60 jr = 1, n
543  temp = max( temp, abs( vr( jr, jc ) ) )
544  60 CONTINUE
545  ELSE
546  DO 70 jr = 1, n
547  temp = max( temp, abs( vr( jr, jc ) )+
548  $ abs( vr( jr, jc+1 ) ) )
549  70 CONTINUE
550  END IF
551  IF( temp.LT.smlnum )
552  $ GO TO 100
553  temp = one / temp
554  IF( alphai( jc ).EQ.zero ) THEN
555  DO 80 jr = 1, n
556  vr( jr, jc ) = vr( jr, jc )*temp
557  80 CONTINUE
558  ELSE
559  DO 90 jr = 1, n
560  vr( jr, jc ) = vr( jr, jc )*temp
561  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
562  90 CONTINUE
563  END IF
564  100 CONTINUE
565  END IF
566 *
567 * End of eigenvector calculation
568 *
569  END IF
570 *
571 * Undo scaling if necessary
572 *
573  110 CONTINUE
574 *
575  IF( ilascl ) THEN
576  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
577  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
578  END IF
579 *
580  IF( ilbscl ) THEN
581  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
582  END IF
583 *
584  work( 1 ) = maxwrk
585  RETURN
586 *
587 * End of DGGEV
588 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
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:110
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 dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
Definition: dggbak.f:147
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
Definition: dggbal.f:177
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:114
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:304
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:145
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
Definition: dtgevc.f:295
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:128
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
Definition: dgghrd.f:207
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: