LAPACK  3.8.0
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.
Date
January 2015

Definition at line 228 of file dggev3.f.

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