LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ dggev3()

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

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

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

Purpose:
 DGGEV3 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

          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 223 of file dggev3.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, LWKOPT
252  DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
253  $ SMLNUM, TEMP
254 * ..
255 * .. Local Arrays ..
256  LOGICAL LDUMMA( 1 )
257 * ..
258 * .. External Subroutines ..
259  EXTERNAL dgeqrf, dggbak, dggbal, dgghd3, dhgeqz, dlabad,
261  $ xerbla
262 * ..
263 * .. External Functions ..
264  LOGICAL LSAME
265  DOUBLE PRECISION DLAMCH, DLANGE
266  EXTERNAL lsame, dlamch, dlange
267 * ..
268 * .. Intrinsic Functions ..
269  INTRINSIC abs, max, sqrt
270 * ..
271 * .. Executable Statements ..
272 *
273 * Decode the input arguments
274 *
275  IF( lsame( jobvl, 'N' ) ) THEN
276  ijobvl = 1
277  ilvl = .false.
278  ELSE IF( lsame( jobvl, 'V' ) ) THEN
279  ijobvl = 2
280  ilvl = .true.
281  ELSE
282  ijobvl = -1
283  ilvl = .false.
284  END IF
285 *
286  IF( lsame( jobvr, 'N' ) ) THEN
287  ijobvr = 1
288  ilvr = .false.
289  ELSE IF( lsame( jobvr, 'V' ) ) THEN
290  ijobvr = 2
291  ilvr = .true.
292  ELSE
293  ijobvr = -1
294  ilvr = .false.
295  END IF
296  ilv = ilvl .OR. ilvr
297 *
298 * Test the input arguments
299 *
300  info = 0
301  lquery = ( lwork.EQ.-1 )
302  IF( ijobvl.LE.0 ) THEN
303  info = -1
304  ELSE IF( ijobvr.LE.0 ) THEN
305  info = -2
306  ELSE IF( n.LT.0 ) THEN
307  info = -3
308  ELSE IF( lda.LT.max( 1, n ) ) THEN
309  info = -5
310  ELSE IF( ldb.LT.max( 1, n ) ) THEN
311  info = -7
312  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
313  info = -12
314  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
315  info = -14
316  ELSE IF( lwork.LT.max( 1, 8*n ) .AND. .NOT.lquery ) THEN
317  info = -16
318  END IF
319 *
320 * Compute workspace
321 *
322  IF( info.EQ.0 ) THEN
323  CALL dgeqrf( n, n, b, ldb, work, work, -1, ierr )
324  lwkopt = max(1, 8*n, 3*n+int( work( 1 ) ) )
325  CALL dormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work, -1,
326  $ ierr )
327  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
328  IF( ilvl ) THEN
329  CALL dorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
330  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
331  END IF
332  IF( ilv ) THEN
333  CALL dgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
334  $ ldvl, vr, ldvr, work, -1, ierr )
335  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
336  CALL dhgeqz( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
337  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
338  $ work, -1, ierr )
339  lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
340  ELSE
341  CALL dgghd3( 'N', 'N', n, 1, n, a, lda, b, ldb, vl, ldvl,
342  $ vr, ldvr, work, -1, ierr )
343  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
344  CALL dhgeqz( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
345  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
346  $ work, -1, ierr )
347  lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
348  END IF
349 
350  work( 1 ) = lwkopt
351  END IF
352 *
353  IF( info.NE.0 ) THEN
354  CALL xerbla( 'DGGEV3 ', -info )
355  RETURN
356  ELSE IF( lquery ) THEN
357  RETURN
358  END IF
359 *
360 * Quick return if possible
361 *
362  IF( n.EQ.0 )
363  $ RETURN
364 *
365 * Get machine constants
366 *
367  eps = dlamch( 'P' )
368  smlnum = dlamch( 'S' )
369  bignum = one / smlnum
370  CALL dlabad( smlnum, bignum )
371  smlnum = sqrt( smlnum ) / eps
372  bignum = one / smlnum
373 *
374 * Scale A if max element outside range [SMLNUM,BIGNUM]
375 *
376  anrm = dlange( 'M', n, n, a, lda, work )
377  ilascl = .false.
378  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
379  anrmto = smlnum
380  ilascl = .true.
381  ELSE IF( anrm.GT.bignum ) THEN
382  anrmto = bignum
383  ilascl = .true.
384  END IF
385  IF( ilascl )
386  $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
387 *
388 * Scale B if max element outside range [SMLNUM,BIGNUM]
389 *
390  bnrm = dlange( 'M', n, n, b, ldb, work )
391  ilbscl = .false.
392  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
393  bnrmto = smlnum
394  ilbscl = .true.
395  ELSE IF( bnrm.GT.bignum ) THEN
396  bnrmto = bignum
397  ilbscl = .true.
398  END IF
399  IF( ilbscl )
400  $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
401 *
402 * Permute the matrices A, B to isolate eigenvalues if possible
403 *
404  ileft = 1
405  iright = n + 1
406  iwrk = iright + n
407  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
408  $ work( iright ), work( iwrk ), ierr )
409 *
410 * Reduce B to triangular form (QR decomposition of B)
411 *
412  irows = ihi + 1 - ilo
413  IF( ilv ) THEN
414  icols = n + 1 - ilo
415  ELSE
416  icols = irows
417  END IF
418  itau = iwrk
419  iwrk = itau + irows
420  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
421  $ work( iwrk ), lwork+1-iwrk, ierr )
422 *
423 * Apply the orthogonal transformation to matrix A
424 *
425  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
426  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
427  $ lwork+1-iwrk, ierr )
428 *
429 * Initialize VL
430 *
431  IF( ilvl ) THEN
432  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
433  IF( irows.GT.1 ) THEN
434  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
435  $ vl( ilo+1, ilo ), ldvl )
436  END IF
437  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
438  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
439  END IF
440 *
441 * Initialize VR
442 *
443  IF( ilvr )
444  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
445 *
446 * Reduce to generalized Hessenberg form
447 *
448  IF( ilv ) THEN
449 *
450 * Eigenvectors requested -- work on whole matrix.
451 *
452  CALL dgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
453  $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
454  ELSE
455  CALL dgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
456  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
457  $ work( iwrk ), lwork+1-iwrk, ierr )
458  END IF
459 *
460 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
461 * Schur forms and Schur vectors)
462 *
463  iwrk = itau
464  IF( ilv ) THEN
465  chtemp = 'S'
466  ELSE
467  chtemp = 'E'
468  END IF
469  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
470  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
471  $ work( iwrk ), lwork+1-iwrk, 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 110
481  END IF
482 *
483 * Compute Eigenvectors
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 *
504  IF( ilvl ) THEN
505  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
506  $ work( iright ), n, vl, ldvl, ierr )
507  DO 50 jc = 1, n
508  IF( alphai( jc ).LT.zero )
509  $ GO TO 50
510  temp = zero
511  IF( alphai( jc ).EQ.zero ) THEN
512  DO 10 jr = 1, n
513  temp = max( temp, abs( vl( jr, jc ) ) )
514  10 CONTINUE
515  ELSE
516  DO 20 jr = 1, n
517  temp = max( temp, abs( vl( jr, jc ) )+
518  $ abs( vl( jr, jc+1 ) ) )
519  20 CONTINUE
520  END IF
521  IF( temp.LT.smlnum )
522  $ GO TO 50
523  temp = one / temp
524  IF( alphai( jc ).EQ.zero ) THEN
525  DO 30 jr = 1, n
526  vl( jr, jc ) = vl( jr, jc )*temp
527  30 CONTINUE
528  ELSE
529  DO 40 jr = 1, n
530  vl( jr, jc ) = vl( jr, jc )*temp
531  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
532  40 CONTINUE
533  END IF
534  50 CONTINUE
535  END IF
536  IF( ilvr ) THEN
537  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
538  $ work( iright ), n, vr, ldvr, ierr )
539  DO 100 jc = 1, n
540  IF( alphai( jc ).LT.zero )
541  $ GO TO 100
542  temp = zero
543  IF( alphai( jc ).EQ.zero ) THEN
544  DO 60 jr = 1, n
545  temp = max( temp, abs( vr( jr, jc ) ) )
546  60 CONTINUE
547  ELSE
548  DO 70 jr = 1, n
549  temp = max( temp, abs( vr( jr, jc ) )+
550  $ abs( vr( jr, jc+1 ) ) )
551  70 CONTINUE
552  END IF
553  IF( temp.LT.smlnum )
554  $ GO TO 100
555  temp = one / temp
556  IF( alphai( jc ).EQ.zero ) THEN
557  DO 80 jr = 1, n
558  vr( jr, jc ) = vr( jr, jc )*temp
559  80 CONTINUE
560  ELSE
561  DO 90 jr = 1, n
562  vr( jr, jc ) = vr( jr, jc )*temp
563  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
564  90 CONTINUE
565  END IF
566  100 CONTINUE
567  END IF
568 *
569 * End of eigenvector calculation
570 *
571  END IF
572 *
573 * Undo scaling if necessary
574 *
575  110 CONTINUE
576 *
577  IF( ilascl ) THEN
578  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
579  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
580  END IF
581 *
582  IF( ilbscl ) THEN
583  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
584  END IF
585 *
586  work( 1 ) = lwkopt
587  RETURN
588 *
589 * End of DGGEV3
590 *
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
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 dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
subroutine dgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DGGHD3
Definition: dgghd3.f:230
Here is the call graph for this function:
Here is the caller graph for this function: