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

◆ dggevx()

subroutine dggevx ( character  balanc,
character  jobvl,
character  jobvr,
character  sense,
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,
integer  ilo,
integer  ihi,
double precision, dimension( * )  lscale,
double precision, dimension( * )  rscale,
double precision  abnrm,
double precision  bbnrm,
double precision, dimension( * )  rconde,
double precision, dimension( * )  rcondv,
double precision, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
logical, dimension( * )  bwork,
integer  info 
)

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

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

Purpose:
 DGGEVX 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.

 Optionally also, it computes a balancing transformation to improve
 the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
 LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
 the eigenvalues (RCONDE), and reciprocal condition numbers for the
 right eigenvectors (RCONDV).

 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]BALANC
          BALANC is CHARACTER*1
          Specifies the balance option to be performed.
          = 'N':  do not diagonally scale or permute;
          = 'P':  permute only;
          = 'S':  scale only;
          = 'B':  both permute and scale.
          Computed reciprocal condition numbers will be for the
          matrices after permuting and/or balancing. Permuting does
          not change condition numbers (in exact arithmetic), but
          balancing does.
[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]SENSE
          SENSE is CHARACTER*1
          Determines which reciprocal condition numbers are computed.
          = 'N': none are computed;
          = 'E': computed for eigenvalues only;
          = 'V': computed for eigenvectors only;
          = 'B': computed for eigenvalues and 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. If JOBVL='V' or JOBVR='V'
          or both, then A contains the first part of the real Schur
          form of the "balanced" versions of the input A and B.
[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. If JOBVL='V' or JOBVR='V'
          or both, then B contains the second part of the real Schur
          form of the "balanced" versions of the input A and B.
[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 will be scaled so the largest component have
          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 will be scaled so the largest component have
          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]ILO
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are integer values such that on exit
          A(i,j) = 0 and B(i,j) = 0 if i > j and
          j = 1,...,ILO-1 or i = IHI+1,...,N.
          If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
[out]LSCALE
          LSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the left side of A and B.  If PL(j) is the index of the
          row interchanged with row j, and DL(j) is the scaling
          factor applied to row j, then
            LSCALE(j) = PL(j)  for j = 1,...,ILO-1
                      = DL(j)  for j = ILO,...,IHI
                      = PL(j)  for j = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]RSCALE
          RSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the right side of A and B.  If PR(j) is the index of the
          column interchanged with column j, and DR(j) is the scaling
          factor applied to column j, then
            RSCALE(j) = PR(j)  for j = 1,...,ILO-1
                      = DR(j)  for j = ILO,...,IHI
                      = PR(j)  for j = IHI+1,...,N
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]ABNRM
          ABNRM is DOUBLE PRECISION
          The one-norm of the balanced matrix A.
[out]BBNRM
          BBNRM is DOUBLE PRECISION
          The one-norm of the balanced matrix B.
[out]RCONDE
          RCONDE is DOUBLE PRECISION array, dimension (N)
          If SENSE = 'E' or 'B', the reciprocal condition numbers of
          the eigenvalues, stored in consecutive elements of the array.
          For a complex conjugate pair of eigenvalues two consecutive
          elements of RCONDE are set to the same value. Thus RCONDE(j),
          RCONDV(j), and the j-th columns of VL and VR all correspond
          to the j-th eigenpair.
          If SENSE = 'N or 'V', RCONDE is not referenced.
[out]RCONDV
          RCONDV is DOUBLE PRECISION array, dimension (N)
          If SENSE = 'V' or 'B', the estimated reciprocal condition
          numbers of the eigenvectors, stored in consecutive elements
          of the array. For a complex eigenvector two consecutive
          elements of RCONDV are set to the same value. If the
          eigenvalues cannot be reordered to compute RCONDV(j),
          RCONDV(j) is set to 0; this can only occur when the true
          value would be very small anyway.
          If SENSE = 'N' or 'E', RCONDV is not referenced.
[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,2*N).
          If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
          LWORK >= max(1,6*N).
          If SENSE = 'E' or 'B', LWORK >= max(1,10*N).
          If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.

          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]IWORK
          IWORK is INTEGER array, dimension (N+6)
          If SENSE = 'E', IWORK is not referenced.
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
          If SENSE = 'N', BWORK is not referenced.
[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.
Further Details:
  Balancing a matrix pair (A,B) includes, first, permuting rows and
  columns to isolate eigenvalues, second, applying diagonal similarity
  transformation to the rows and columns to make the rows and columns
  as close in norm as possible. The computed reciprocal condition
  numbers correspond to the balanced matrix. Permuting rows and columns
  will not change the condition numbers (in exact arithmetic) but
  diagonal scaling will.  For further explanation of balancing, see
  section 4.11.1.2 of LAPACK Users' Guide.

  An approximate error bound on the chordal distance between the i-th
  computed generalized eigenvalue w and the corresponding exact
  eigenvalue lambda is

       chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)

  An approximate error bound for the angle between the i-th computed
  eigenvector VL(i) or VR(i) is given by

       EPS * norm(ABNRM, BBNRM) / DIF(i).

  For further explanation of the reciprocal condition numbers RCONDE
  and RCONDV, see section 4.11 of LAPACK User's Guide.

Definition at line 387 of file dggevx.f.

391*
392* -- LAPACK driver routine --
393* -- LAPACK is a software package provided by Univ. of Tennessee, --
394* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
395*
396* .. Scalar Arguments ..
397 CHARACTER BALANC, JOBVL, JOBVR, SENSE
398 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
399 DOUBLE PRECISION ABNRM, BBNRM
400* ..
401* .. Array Arguments ..
402 LOGICAL BWORK( * )
403 INTEGER IWORK( * )
404 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
405 $ B( LDB, * ), BETA( * ), LSCALE( * ),
406 $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
407 $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
408* ..
409*
410* =====================================================================
411*
412* .. Parameters ..
413 DOUBLE PRECISION ZERO, ONE
414 parameter( zero = 0.0d+0, one = 1.0d+0 )
415* ..
416* .. Local Scalars ..
417 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
418 $ PAIR, WANTSB, WANTSE, WANTSN, WANTSV
419 CHARACTER CHTEMP
420 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
421 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
422 $ MINWRK, MM
423 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
424 $ SMLNUM, TEMP
425* ..
426* .. Local Arrays ..
427 LOGICAL LDUMMA( 1 )
428* ..
429* .. External Subroutines ..
430 EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
432 $ xerbla
433* ..
434* .. External Functions ..
435 LOGICAL LSAME
436 INTEGER ILAENV
437 DOUBLE PRECISION DLAMCH, DLANGE
438 EXTERNAL lsame, ilaenv, dlamch, dlange
439* ..
440* .. Intrinsic Functions ..
441 INTRINSIC abs, max, sqrt
442* ..
443* .. Executable Statements ..
444*
445* Decode the input arguments
446*
447 IF( lsame( jobvl, 'N' ) ) THEN
448 ijobvl = 1
449 ilvl = .false.
450 ELSE IF( lsame( jobvl, 'V' ) ) THEN
451 ijobvl = 2
452 ilvl = .true.
453 ELSE
454 ijobvl = -1
455 ilvl = .false.
456 END IF
457*
458 IF( lsame( jobvr, 'N' ) ) THEN
459 ijobvr = 1
460 ilvr = .false.
461 ELSE IF( lsame( jobvr, 'V' ) ) THEN
462 ijobvr = 2
463 ilvr = .true.
464 ELSE
465 ijobvr = -1
466 ilvr = .false.
467 END IF
468 ilv = ilvl .OR. ilvr
469*
470 noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
471 wantsn = lsame( sense, 'N' )
472 wantse = lsame( sense, 'E' )
473 wantsv = lsame( sense, 'V' )
474 wantsb = lsame( sense, 'B' )
475*
476* Test the input arguments
477*
478 info = 0
479 lquery = ( lwork.EQ.-1 )
480 IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc,
481 $ 'S' ) .OR. lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) )
482 $ THEN
483 info = -1
484 ELSE IF( ijobvl.LE.0 ) THEN
485 info = -2
486 ELSE IF( ijobvr.LE.0 ) THEN
487 info = -3
488 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
489 $ THEN
490 info = -4
491 ELSE IF( n.LT.0 ) THEN
492 info = -5
493 ELSE IF( lda.LT.max( 1, n ) ) THEN
494 info = -7
495 ELSE IF( ldb.LT.max( 1, n ) ) THEN
496 info = -9
497 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
498 info = -14
499 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
500 info = -16
501 END IF
502*
503* Compute workspace
504* (Note: Comments in the code beginning "Workspace:" describe the
505* minimal amount of workspace needed at that point in the code,
506* as well as the preferred amount for good performance.
507* NB refers to the optimal block size for the immediately
508* following subroutine, as returned by ILAENV. The workspace is
509* computed assuming ILO = 1 and IHI = N, the worst case.)
510*
511 IF( info.EQ.0 ) THEN
512 IF( n.EQ.0 ) THEN
513 minwrk = 1
514 maxwrk = 1
515 ELSE
516 IF( noscl .AND. .NOT.ilv ) THEN
517 minwrk = 2*n
518 ELSE
519 minwrk = 6*n
520 END IF
521 IF( wantse .OR. wantsb ) THEN
522 minwrk = 10*n
523 END IF
524 IF( wantsv .OR. wantsb ) THEN
525 minwrk = max( minwrk, 2*n*( n + 4 ) + 16 )
526 END IF
527 maxwrk = minwrk
528 maxwrk = max( maxwrk,
529 $ n + n*ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 ) )
530 maxwrk = max( maxwrk,
531 $ n + n*ilaenv( 1, 'DORMQR', ' ', n, 1, n, 0 ) )
532 IF( ilvl ) THEN
533 maxwrk = max( maxwrk, n +
534 $ n*ilaenv( 1, 'DORGQR', ' ', n, 1, n, 0 ) )
535 END IF
536 END IF
537 work( 1 ) = maxwrk
538*
539 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
540 info = -26
541 END IF
542 END IF
543*
544 IF( info.NE.0 ) THEN
545 CALL xerbla( 'DGGEVX', -info )
546 RETURN
547 ELSE IF( lquery ) THEN
548 RETURN
549 END IF
550*
551* Quick return if possible
552*
553 IF( n.EQ.0 )
554 $ RETURN
555*
556*
557* Get machine constants
558*
559 eps = dlamch( 'P' )
560 smlnum = dlamch( 'S' )
561 bignum = one / smlnum
562 smlnum = sqrt( smlnum ) / eps
563 bignum = one / smlnum
564*
565* Scale A if max element outside range [SMLNUM,BIGNUM]
566*
567 anrm = dlange( 'M', n, n, a, lda, work )
568 ilascl = .false.
569 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
570 anrmto = smlnum
571 ilascl = .true.
572 ELSE IF( anrm.GT.bignum ) THEN
573 anrmto = bignum
574 ilascl = .true.
575 END IF
576 IF( ilascl )
577 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
578*
579* Scale B if max element outside range [SMLNUM,BIGNUM]
580*
581 bnrm = dlange( 'M', n, n, b, ldb, work )
582 ilbscl = .false.
583 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
584 bnrmto = smlnum
585 ilbscl = .true.
586 ELSE IF( bnrm.GT.bignum ) THEN
587 bnrmto = bignum
588 ilbscl = .true.
589 END IF
590 IF( ilbscl )
591 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
592*
593* Permute and/or balance the matrix pair (A,B)
594* (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
595*
596 CALL dggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
597 $ work, ierr )
598*
599* Compute ABNRM and BBNRM
600*
601 abnrm = dlange( '1', n, n, a, lda, work( 1 ) )
602 IF( ilascl ) THEN
603 work( 1 ) = abnrm
604 CALL dlascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
605 $ ierr )
606 abnrm = work( 1 )
607 END IF
608*
609 bbnrm = dlange( '1', n, n, b, ldb, work( 1 ) )
610 IF( ilbscl ) THEN
611 work( 1 ) = bbnrm
612 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
613 $ ierr )
614 bbnrm = work( 1 )
615 END IF
616*
617* Reduce B to triangular form (QR decomposition of B)
618* (Workspace: need N, prefer N*NB )
619*
620 irows = ihi + 1 - ilo
621 IF( ilv .OR. .NOT.wantsn ) THEN
622 icols = n + 1 - ilo
623 ELSE
624 icols = irows
625 END IF
626 itau = 1
627 iwrk = itau + irows
628 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
629 $ work( iwrk ), lwork+1-iwrk, ierr )
630*
631* Apply the orthogonal transformation to A
632* (Workspace: need N, prefer N*NB)
633*
634 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
635 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
636 $ lwork+1-iwrk, ierr )
637*
638* Initialize VL and/or VR
639* (Workspace: need N, prefer N*NB)
640*
641 IF( ilvl ) THEN
642 CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
643 IF( irows.GT.1 ) THEN
644 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
645 $ vl( ilo+1, ilo ), ldvl )
646 END IF
647 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
648 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
649 END IF
650*
651 IF( ilvr )
652 $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
653*
654* Reduce to generalized Hessenberg form
655* (Workspace: none needed)
656*
657 IF( ilv .OR. .NOT.wantsn ) THEN
658*
659* Eigenvectors requested -- work on whole matrix.
660*
661 CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
662 $ ldvl, vr, ldvr, ierr )
663 ELSE
664 CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
665 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
666 END IF
667*
668* Perform QZ algorithm (Compute eigenvalues, and optionally, the
669* Schur forms and Schur vectors)
670* (Workspace: need N)
671*
672 IF( ilv .OR. .NOT.wantsn ) THEN
673 chtemp = 'S'
674 ELSE
675 chtemp = 'E'
676 END IF
677*
678 CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
679 $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
680 $ lwork, ierr )
681 IF( ierr.NE.0 ) THEN
682 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
683 info = ierr
684 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
685 info = ierr - n
686 ELSE
687 info = n + 1
688 END IF
689 GO TO 130
690 END IF
691*
692* Compute Eigenvectors and estimate condition numbers if desired
693* (Workspace: DTGEVC: need 6*N
694* DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
695* need N otherwise )
696*
697 IF( ilv .OR. .NOT.wantsn ) THEN
698 IF( ilv ) THEN
699 IF( ilvl ) THEN
700 IF( ilvr ) THEN
701 chtemp = 'B'
702 ELSE
703 chtemp = 'L'
704 END IF
705 ELSE
706 chtemp = 'R'
707 END IF
708*
709 CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
710 $ ldvl, vr, ldvr, n, in, work, ierr )
711 IF( ierr.NE.0 ) THEN
712 info = n + 2
713 GO TO 130
714 END IF
715 END IF
716*
717 IF( .NOT.wantsn ) THEN
718*
719* compute eigenvectors (DTGEVC) and estimate condition
720* numbers (DTGSNA). Note that the definition of the condition
721* number is not invariant under transformation (u,v) to
722* (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
723* Schur form (S,T), Q and Z are orthogonal matrices. In order
724* to avoid using extra 2*N*N workspace, we have to recalculate
725* eigenvectors and estimate one condition numbers at a time.
726*
727 pair = .false.
728 DO 20 i = 1, n
729*
730 IF( pair ) THEN
731 pair = .false.
732 GO TO 20
733 END IF
734 mm = 1
735 IF( i.LT.n ) THEN
736 IF( a( i+1, i ).NE.zero ) THEN
737 pair = .true.
738 mm = 2
739 END IF
740 END IF
741*
742 DO 10 j = 1, n
743 bwork( j ) = .false.
744 10 CONTINUE
745 IF( mm.EQ.1 ) THEN
746 bwork( i ) = .true.
747 ELSE IF( mm.EQ.2 ) THEN
748 bwork( i ) = .true.
749 bwork( i+1 ) = .true.
750 END IF
751*
752 iwrk = mm*n + 1
753 iwrk1 = iwrk + mm*n
754*
755* Compute a pair of left and right eigenvectors.
756* (compute workspace: need up to 4*N + 6*N)
757*
758 IF( wantse .OR. wantsb ) THEN
759 CALL dtgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
760 $ work( 1 ), n, work( iwrk ), n, mm, m,
761 $ work( iwrk1 ), ierr )
762 IF( ierr.NE.0 ) THEN
763 info = n + 2
764 GO TO 130
765 END IF
766 END IF
767*
768 CALL dtgsna( sense, 'S', bwork, n, a, lda, b, ldb,
769 $ work( 1 ), n, work( iwrk ), n, rconde( i ),
770 $ rcondv( i ), mm, m, work( iwrk1 ),
771 $ lwork-iwrk1+1, iwork, ierr )
772*
773 20 CONTINUE
774 END IF
775 END IF
776*
777* Undo balancing on VL and VR and normalization
778* (Workspace: none needed)
779*
780 IF( ilvl ) THEN
781 CALL dggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
782 $ ldvl, ierr )
783*
784 DO 70 jc = 1, n
785 IF( alphai( jc ).LT.zero )
786 $ GO TO 70
787 temp = zero
788 IF( alphai( jc ).EQ.zero ) THEN
789 DO 30 jr = 1, n
790 temp = max( temp, abs( vl( jr, jc ) ) )
791 30 CONTINUE
792 ELSE
793 DO 40 jr = 1, n
794 temp = max( temp, abs( vl( jr, jc ) )+
795 $ abs( vl( jr, jc+1 ) ) )
796 40 CONTINUE
797 END IF
798 IF( temp.LT.smlnum )
799 $ GO TO 70
800 temp = one / temp
801 IF( alphai( jc ).EQ.zero ) THEN
802 DO 50 jr = 1, n
803 vl( jr, jc ) = vl( jr, jc )*temp
804 50 CONTINUE
805 ELSE
806 DO 60 jr = 1, n
807 vl( jr, jc ) = vl( jr, jc )*temp
808 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
809 60 CONTINUE
810 END IF
811 70 CONTINUE
812 END IF
813 IF( ilvr ) THEN
814 CALL dggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
815 $ ldvr, ierr )
816 DO 120 jc = 1, n
817 IF( alphai( jc ).LT.zero )
818 $ GO TO 120
819 temp = zero
820 IF( alphai( jc ).EQ.zero ) THEN
821 DO 80 jr = 1, n
822 temp = max( temp, abs( vr( jr, jc ) ) )
823 80 CONTINUE
824 ELSE
825 DO 90 jr = 1, n
826 temp = max( temp, abs( vr( jr, jc ) )+
827 $ abs( vr( jr, jc+1 ) ) )
828 90 CONTINUE
829 END IF
830 IF( temp.LT.smlnum )
831 $ GO TO 120
832 temp = one / temp
833 IF( alphai( jc ).EQ.zero ) THEN
834 DO 100 jr = 1, n
835 vr( jr, jc ) = vr( jr, jc )*temp
836 100 CONTINUE
837 ELSE
838 DO 110 jr = 1, n
839 vr( jr, jc ) = vr( jr, jc )*temp
840 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
841 110 CONTINUE
842 END IF
843 120 CONTINUE
844 END IF
845*
846* Undo scaling if necessary
847*
848 130 CONTINUE
849*
850 IF( ilascl ) THEN
851 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
852 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
853 END IF
854*
855 IF( ilbscl ) THEN
856 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
857 END IF
858*
859 work( 1 ) = maxwrk
860 RETURN
861*
862* End of DGGEVX
863*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
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
subroutine dgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
DGGHRD
Definition dgghrd.f:207
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dtgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
DTGEVC
Definition dtgevc.f:295
subroutine dtgsna(job, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, ldvr, s, dif, mm, m, work, lwork, iwork, info)
DTGSNA
Definition dtgsna.f:381
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
Here is the call graph for this function:
Here is the caller graph for this function: