LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ 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, dlabad,
432  $ dtgsna, 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  CALL dlabad( smlnum, bignum )
563  smlnum = sqrt( smlnum ) / eps
564  bignum = one / smlnum
565 *
566 * Scale A if max element outside range [SMLNUM,BIGNUM]
567 *
568  anrm = dlange( 'M', n, n, a, lda, work )
569  ilascl = .false.
570  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
571  anrmto = smlnum
572  ilascl = .true.
573  ELSE IF( anrm.GT.bignum ) THEN
574  anrmto = bignum
575  ilascl = .true.
576  END IF
577  IF( ilascl )
578  $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
579 *
580 * Scale B if max element outside range [SMLNUM,BIGNUM]
581 *
582  bnrm = dlange( 'M', n, n, b, ldb, work )
583  ilbscl = .false.
584  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
585  bnrmto = smlnum
586  ilbscl = .true.
587  ELSE IF( bnrm.GT.bignum ) THEN
588  bnrmto = bignum
589  ilbscl = .true.
590  END IF
591  IF( ilbscl )
592  $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
593 *
594 * Permute and/or balance the matrix pair (A,B)
595 * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
596 *
597  CALL dggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
598  $ work, ierr )
599 *
600 * Compute ABNRM and BBNRM
601 *
602  abnrm = dlange( '1', n, n, a, lda, work( 1 ) )
603  IF( ilascl ) THEN
604  work( 1 ) = abnrm
605  CALL dlascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
606  $ ierr )
607  abnrm = work( 1 )
608  END IF
609 *
610  bbnrm = dlange( '1', n, n, b, ldb, work( 1 ) )
611  IF( ilbscl ) THEN
612  work( 1 ) = bbnrm
613  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
614  $ ierr )
615  bbnrm = work( 1 )
616  END IF
617 *
618 * Reduce B to triangular form (QR decomposition of B)
619 * (Workspace: need N, prefer N*NB )
620 *
621  irows = ihi + 1 - ilo
622  IF( ilv .OR. .NOT.wantsn ) THEN
623  icols = n + 1 - ilo
624  ELSE
625  icols = irows
626  END IF
627  itau = 1
628  iwrk = itau + irows
629  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
630  $ work( iwrk ), lwork+1-iwrk, ierr )
631 *
632 * Apply the orthogonal transformation to A
633 * (Workspace: need N, prefer N*NB)
634 *
635  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
636  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
637  $ lwork+1-iwrk, ierr )
638 *
639 * Initialize VL and/or VR
640 * (Workspace: need N, prefer N*NB)
641 *
642  IF( ilvl ) THEN
643  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
644  IF( irows.GT.1 ) THEN
645  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
646  $ vl( ilo+1, ilo ), ldvl )
647  END IF
648  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
649  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
650  END IF
651 *
652  IF( ilvr )
653  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
654 *
655 * Reduce to generalized Hessenberg form
656 * (Workspace: none needed)
657 *
658  IF( ilv .OR. .NOT.wantsn ) THEN
659 *
660 * Eigenvectors requested -- work on whole matrix.
661 *
662  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
663  $ ldvl, vr, ldvr, ierr )
664  ELSE
665  CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
666  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
667  END IF
668 *
669 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
670 * Schur forms and Schur vectors)
671 * (Workspace: need N)
672 *
673  IF( ilv .OR. .NOT.wantsn ) THEN
674  chtemp = 'S'
675  ELSE
676  chtemp = 'E'
677  END IF
678 *
679  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
680  $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
681  $ lwork, ierr )
682  IF( ierr.NE.0 ) THEN
683  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
684  info = ierr
685  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
686  info = ierr - n
687  ELSE
688  info = n + 1
689  END IF
690  GO TO 130
691  END IF
692 *
693 * Compute Eigenvectors and estimate condition numbers if desired
694 * (Workspace: DTGEVC: need 6*N
695 * DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
696 * need N otherwise )
697 *
698  IF( ilv .OR. .NOT.wantsn ) THEN
699  IF( ilv ) THEN
700  IF( ilvl ) THEN
701  IF( ilvr ) THEN
702  chtemp = 'B'
703  ELSE
704  chtemp = 'L'
705  END IF
706  ELSE
707  chtemp = 'R'
708  END IF
709 *
710  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
711  $ ldvl, vr, ldvr, n, in, work, ierr )
712  IF( ierr.NE.0 ) THEN
713  info = n + 2
714  GO TO 130
715  END IF
716  END IF
717 *
718  IF( .NOT.wantsn ) THEN
719 *
720 * compute eigenvectors (DTGEVC) and estimate condition
721 * numbers (DTGSNA). Note that the definition of the condition
722 * number is not invariant under transformation (u,v) to
723 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
724 * Schur form (S,T), Q and Z are orthogonal matrices. In order
725 * to avoid using extra 2*N*N workspace, we have to recalculate
726 * eigenvectors and estimate one condition numbers at a time.
727 *
728  pair = .false.
729  DO 20 i = 1, n
730 *
731  IF( pair ) THEN
732  pair = .false.
733  GO TO 20
734  END IF
735  mm = 1
736  IF( i.LT.n ) THEN
737  IF( a( i+1, i ).NE.zero ) THEN
738  pair = .true.
739  mm = 2
740  END IF
741  END IF
742 *
743  DO 10 j = 1, n
744  bwork( j ) = .false.
745  10 CONTINUE
746  IF( mm.EQ.1 ) THEN
747  bwork( i ) = .true.
748  ELSE IF( mm.EQ.2 ) THEN
749  bwork( i ) = .true.
750  bwork( i+1 ) = .true.
751  END IF
752 *
753  iwrk = mm*n + 1
754  iwrk1 = iwrk + mm*n
755 *
756 * Compute a pair of left and right eigenvectors.
757 * (compute workspace: need up to 4*N + 6*N)
758 *
759  IF( wantse .OR. wantsb ) THEN
760  CALL dtgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
761  $ work( 1 ), n, work( iwrk ), n, mm, m,
762  $ work( iwrk1 ), ierr )
763  IF( ierr.NE.0 ) THEN
764  info = n + 2
765  GO TO 130
766  END IF
767  END IF
768 *
769  CALL dtgsna( sense, 'S', bwork, n, a, lda, b, ldb,
770  $ work( 1 ), n, work( iwrk ), n, rconde( i ),
771  $ rcondv( i ), mm, m, work( iwrk1 ),
772  $ lwork-iwrk1+1, iwork, ierr )
773 *
774  20 CONTINUE
775  END IF
776  END IF
777 *
778 * Undo balancing on VL and VR and normalization
779 * (Workspace: none needed)
780 *
781  IF( ilvl ) THEN
782  CALL dggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
783  $ ldvl, ierr )
784 *
785  DO 70 jc = 1, n
786  IF( alphai( jc ).LT.zero )
787  $ GO TO 70
788  temp = zero
789  IF( alphai( jc ).EQ.zero ) THEN
790  DO 30 jr = 1, n
791  temp = max( temp, abs( vl( jr, jc ) ) )
792  30 CONTINUE
793  ELSE
794  DO 40 jr = 1, n
795  temp = max( temp, abs( vl( jr, jc ) )+
796  $ abs( vl( jr, jc+1 ) ) )
797  40 CONTINUE
798  END IF
799  IF( temp.LT.smlnum )
800  $ GO TO 70
801  temp = one / temp
802  IF( alphai( jc ).EQ.zero ) THEN
803  DO 50 jr = 1, n
804  vl( jr, jc ) = vl( jr, jc )*temp
805  50 CONTINUE
806  ELSE
807  DO 60 jr = 1, n
808  vl( jr, jc ) = vl( jr, jc )*temp
809  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
810  60 CONTINUE
811  END IF
812  70 CONTINUE
813  END IF
814  IF( ilvr ) THEN
815  CALL dggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
816  $ ldvr, ierr )
817  DO 120 jc = 1, n
818  IF( alphai( jc ).LT.zero )
819  $ GO TO 120
820  temp = zero
821  IF( alphai( jc ).EQ.zero ) THEN
822  DO 80 jr = 1, n
823  temp = max( temp, abs( vr( jr, jc ) ) )
824  80 CONTINUE
825  ELSE
826  DO 90 jr = 1, n
827  temp = max( temp, abs( vr( jr, jc ) )+
828  $ abs( vr( jr, jc+1 ) ) )
829  90 CONTINUE
830  END IF
831  IF( temp.LT.smlnum )
832  $ GO TO 120
833  temp = one / temp
834  IF( alphai( jc ).EQ.zero ) THEN
835  DO 100 jr = 1, n
836  vr( jr, jc ) = vr( jr, jc )*temp
837  100 CONTINUE
838  ELSE
839  DO 110 jr = 1, n
840  vr( jr, jc ) = vr( jr, jc )*temp
841  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
842  110 CONTINUE
843  END IF
844  120 CONTINUE
845  END IF
846 *
847 * Undo scaling if necessary
848 *
849  130 CONTINUE
850 *
851  IF( ilascl ) THEN
852  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
853  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
854  END IF
855 *
856  IF( ilbscl ) THEN
857  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
858  END IF
859 *
860  work( 1 ) = maxwrk
861  RETURN
862 *
863 * End of DGGEVX
864 *
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
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
Here is the call graph for this function:
Here is the caller graph for this function: