LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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:146
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: