LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sggevx()

subroutine sggevx ( character  BALANC,
character  JOBVL,
character  JOBVR,
character  SENSE,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( * )  ALPHAR,
real, dimension( * )  ALPHAI,
real, dimension( * )  BETA,
real, dimension( ldvl, * )  VL,
integer  LDVL,
real, dimension( ldvr, * )  VR,
integer  LDVR,
integer  ILO,
integer  IHI,
real, dimension( * )  LSCALE,
real, dimension( * )  RSCALE,
real  ABNRM,
real  BBNRM,
real, dimension( * )  RCONDE,
real, dimension( * )  RCONDV,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

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

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

Purpose:
 SGGEVX 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 REAL 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 REAL 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 REAL array, dimension (N)
[out]ALPHAI
          ALPHAI is REAL array, dimension (N)
[out]BETA
          BETA is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL
          The one-norm of the balanced matrix A.
[out]BBNRM
          BBNRM is REAL
          The one-norm of the balanced matrix B.
[out]RCONDE
          RCONDE is REAL 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 REAL 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 REAL 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', 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 SHGEQZ.
                =N+2: error return from STGEVC.
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 sggevx.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  REAL ABNRM, BBNRM
400 * ..
401 * .. Array Arguments ..
402  LOGICAL BWORK( * )
403  INTEGER IWORK( * )
404  REAL 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  REAL ZERO, ONE
414  parameter( zero = 0.0e+0, one = 1.0e+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  REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
424  $ SMLNUM, TEMP
425 * ..
426 * .. Local Arrays ..
427  LOGICAL LDUMMA( 1 )
428 * ..
429 * .. External Subroutines ..
430  EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
432  $ stgsna, xerbla
433 * ..
434 * .. External Functions ..
435  LOGICAL LSAME
436  INTEGER ILAENV
437  REAL SLAMCH, SLANGE
438  EXTERNAL lsame, ilaenv, slamch, slange
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.( noscl .OR. lsame( balanc, 'S' ) .OR.
481  $ lsame( balanc, 'B' ) ) ) THEN
482  info = -1
483  ELSE IF( ijobvl.LE.0 ) THEN
484  info = -2
485  ELSE IF( ijobvr.LE.0 ) THEN
486  info = -3
487  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
488  $ THEN
489  info = -4
490  ELSE IF( n.LT.0 ) THEN
491  info = -5
492  ELSE IF( lda.LT.max( 1, n ) ) THEN
493  info = -7
494  ELSE IF( ldb.LT.max( 1, n ) ) THEN
495  info = -9
496  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
497  info = -14
498  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
499  info = -16
500  END IF
501 *
502 * Compute workspace
503 * (Note: Comments in the code beginning "Workspace:" describe the
504 * minimal amount of workspace needed at that point in the code,
505 * as well as the preferred amount for good performance.
506 * NB refers to the optimal block size for the immediately
507 * following subroutine, as returned by ILAENV. The workspace is
508 * computed assuming ILO = 1 and IHI = N, the worst case.)
509 *
510  IF( info.EQ.0 ) THEN
511  IF( n.EQ.0 ) THEN
512  minwrk = 1
513  maxwrk = 1
514  ELSE
515  IF( noscl .AND. .NOT.ilv ) THEN
516  minwrk = 2*n
517  ELSE
518  minwrk = 6*n
519  END IF
520  IF( wantse ) THEN
521  minwrk = 10*n
522  ELSE IF( wantsv .OR. wantsb ) THEN
523  minwrk = 2*n*( n + 4 ) + 16
524  END IF
525  maxwrk = minwrk
526  maxwrk = max( maxwrk,
527  $ n + n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) )
528  maxwrk = max( maxwrk,
529  $ n + n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) )
530  IF( ilvl ) THEN
531  maxwrk = max( maxwrk, n +
532  $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, 0 ) )
533  END IF
534  END IF
535  work( 1 ) = maxwrk
536 *
537  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
538  info = -26
539  END IF
540  END IF
541 *
542  IF( info.NE.0 ) THEN
543  CALL xerbla( 'SGGEVX', -info )
544  RETURN
545  ELSE IF( lquery ) THEN
546  RETURN
547  END IF
548 *
549 * Quick return if possible
550 *
551  IF( n.EQ.0 )
552  $ RETURN
553 *
554 *
555 * Get machine constants
556 *
557  eps = slamch( 'P' )
558  smlnum = slamch( 'S' )
559  bignum = one / smlnum
560  CALL slabad( smlnum, bignum )
561  smlnum = sqrt( smlnum ) / eps
562  bignum = one / smlnum
563 *
564 * Scale A if max element outside range [SMLNUM,BIGNUM]
565 *
566  anrm = slange( 'M', n, n, a, lda, work )
567  ilascl = .false.
568  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
569  anrmto = smlnum
570  ilascl = .true.
571  ELSE IF( anrm.GT.bignum ) THEN
572  anrmto = bignum
573  ilascl = .true.
574  END IF
575  IF( ilascl )
576  $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
577 *
578 * Scale B if max element outside range [SMLNUM,BIGNUM]
579 *
580  bnrm = slange( 'M', n, n, b, ldb, work )
581  ilbscl = .false.
582  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
583  bnrmto = smlnum
584  ilbscl = .true.
585  ELSE IF( bnrm.GT.bignum ) THEN
586  bnrmto = bignum
587  ilbscl = .true.
588  END IF
589  IF( ilbscl )
590  $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
591 *
592 * Permute and/or balance the matrix pair (A,B)
593 * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
594 *
595  CALL sggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
596  $ work, ierr )
597 *
598 * Compute ABNRM and BBNRM
599 *
600  abnrm = slange( '1', n, n, a, lda, work( 1 ) )
601  IF( ilascl ) THEN
602  work( 1 ) = abnrm
603  CALL slascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
604  $ ierr )
605  abnrm = work( 1 )
606  END IF
607 *
608  bbnrm = slange( '1', n, n, b, ldb, work( 1 ) )
609  IF( ilbscl ) THEN
610  work( 1 ) = bbnrm
611  CALL slascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
612  $ ierr )
613  bbnrm = work( 1 )
614  END IF
615 *
616 * Reduce B to triangular form (QR decomposition of B)
617 * (Workspace: need N, prefer N*NB )
618 *
619  irows = ihi + 1 - ilo
620  IF( ilv .OR. .NOT.wantsn ) THEN
621  icols = n + 1 - ilo
622  ELSE
623  icols = irows
624  END IF
625  itau = 1
626  iwrk = itau + irows
627  CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
628  $ work( iwrk ), lwork+1-iwrk, ierr )
629 *
630 * Apply the orthogonal transformation to A
631 * (Workspace: need N, prefer N*NB)
632 *
633  CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
634  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
635  $ lwork+1-iwrk, ierr )
636 *
637 * Initialize VL and/or VR
638 * (Workspace: need N, prefer N*NB)
639 *
640  IF( ilvl ) THEN
641  CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
642  IF( irows.GT.1 ) THEN
643  CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
644  $ vl( ilo+1, ilo ), ldvl )
645  END IF
646  CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
647  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
648  END IF
649 *
650  IF( ilvr )
651  $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
652 *
653 * Reduce to generalized Hessenberg form
654 * (Workspace: none needed)
655 *
656  IF( ilv .OR. .NOT.wantsn ) THEN
657 *
658 * Eigenvectors requested -- work on whole matrix.
659 *
660  CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
661  $ ldvl, vr, ldvr, ierr )
662  ELSE
663  CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
664  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
665  END IF
666 *
667 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
668 * Schur forms and Schur vectors)
669 * (Workspace: need N)
670 *
671  IF( ilv .OR. .NOT.wantsn ) THEN
672  chtemp = 'S'
673  ELSE
674  chtemp = 'E'
675  END IF
676 *
677  CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
678  $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
679  $ lwork, ierr )
680  IF( ierr.NE.0 ) THEN
681  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
682  info = ierr
683  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
684  info = ierr - n
685  ELSE
686  info = n + 1
687  END IF
688  GO TO 130
689  END IF
690 *
691 * Compute Eigenvectors and estimate condition numbers if desired
692 * (Workspace: STGEVC: need 6*N
693 * STGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
694 * need N otherwise )
695 *
696  IF( ilv .OR. .NOT.wantsn ) THEN
697  IF( ilv ) THEN
698  IF( ilvl ) THEN
699  IF( ilvr ) THEN
700  chtemp = 'B'
701  ELSE
702  chtemp = 'L'
703  END IF
704  ELSE
705  chtemp = 'R'
706  END IF
707 *
708  CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
709  $ ldvl, vr, ldvr, n, in, work, ierr )
710  IF( ierr.NE.0 ) THEN
711  info = n + 2
712  GO TO 130
713  END IF
714  END IF
715 *
716  IF( .NOT.wantsn ) THEN
717 *
718 * compute eigenvectors (STGEVC) and estimate condition
719 * numbers (STGSNA). Note that the definition of the condition
720 * number is not invariant under transformation (u,v) to
721 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
722 * Schur form (S,T), Q and Z are orthogonal matrices. In order
723 * to avoid using extra 2*N*N workspace, we have to recalculate
724 * eigenvectors and estimate one condition numbers at a time.
725 *
726  pair = .false.
727  DO 20 i = 1, n
728 *
729  IF( pair ) THEN
730  pair = .false.
731  GO TO 20
732  END IF
733  mm = 1
734  IF( i.LT.n ) THEN
735  IF( a( i+1, i ).NE.zero ) THEN
736  pair = .true.
737  mm = 2
738  END IF
739  END IF
740 *
741  DO 10 j = 1, n
742  bwork( j ) = .false.
743  10 CONTINUE
744  IF( mm.EQ.1 ) THEN
745  bwork( i ) = .true.
746  ELSE IF( mm.EQ.2 ) THEN
747  bwork( i ) = .true.
748  bwork( i+1 ) = .true.
749  END IF
750 *
751  iwrk = mm*n + 1
752  iwrk1 = iwrk + mm*n
753 *
754 * Compute a pair of left and right eigenvectors.
755 * (compute workspace: need up to 4*N + 6*N)
756 *
757  IF( wantse .OR. wantsb ) THEN
758  CALL stgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
759  $ work( 1 ), n, work( iwrk ), n, mm, m,
760  $ work( iwrk1 ), ierr )
761  IF( ierr.NE.0 ) THEN
762  info = n + 2
763  GO TO 130
764  END IF
765  END IF
766 *
767  CALL stgsna( sense, 'S', bwork, n, a, lda, b, ldb,
768  $ work( 1 ), n, work( iwrk ), n, rconde( i ),
769  $ rcondv( i ), mm, m, work( iwrk1 ),
770  $ lwork-iwrk1+1, iwork, ierr )
771 *
772  20 CONTINUE
773  END IF
774  END IF
775 *
776 * Undo balancing on VL and VR and normalization
777 * (Workspace: none needed)
778 *
779  IF( ilvl ) THEN
780  CALL sggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
781  $ ldvl, ierr )
782 *
783  DO 70 jc = 1, n
784  IF( alphai( jc ).LT.zero )
785  $ GO TO 70
786  temp = zero
787  IF( alphai( jc ).EQ.zero ) THEN
788  DO 30 jr = 1, n
789  temp = max( temp, abs( vl( jr, jc ) ) )
790  30 CONTINUE
791  ELSE
792  DO 40 jr = 1, n
793  temp = max( temp, abs( vl( jr, jc ) )+
794  $ abs( vl( jr, jc+1 ) ) )
795  40 CONTINUE
796  END IF
797  IF( temp.LT.smlnum )
798  $ GO TO 70
799  temp = one / temp
800  IF( alphai( jc ).EQ.zero ) THEN
801  DO 50 jr = 1, n
802  vl( jr, jc ) = vl( jr, jc )*temp
803  50 CONTINUE
804  ELSE
805  DO 60 jr = 1, n
806  vl( jr, jc ) = vl( jr, jc )*temp
807  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
808  60 CONTINUE
809  END IF
810  70 CONTINUE
811  END IF
812  IF( ilvr ) THEN
813  CALL sggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
814  $ ldvr, ierr )
815  DO 120 jc = 1, n
816  IF( alphai( jc ).LT.zero )
817  $ GO TO 120
818  temp = zero
819  IF( alphai( jc ).EQ.zero ) THEN
820  DO 80 jr = 1, n
821  temp = max( temp, abs( vr( jr, jc ) ) )
822  80 CONTINUE
823  ELSE
824  DO 90 jr = 1, n
825  temp = max( temp, abs( vr( jr, jc ) )+
826  $ abs( vr( jr, jc+1 ) ) )
827  90 CONTINUE
828  END IF
829  IF( temp.LT.smlnum )
830  $ GO TO 120
831  temp = one / temp
832  IF( alphai( jc ).EQ.zero ) THEN
833  DO 100 jr = 1, n
834  vr( jr, jc ) = vr( jr, jc )*temp
835  100 CONTINUE
836  ELSE
837  DO 110 jr = 1, n
838  vr( jr, jc ) = vr( jr, jc )*temp
839  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
840  110 CONTINUE
841  END IF
842  120 CONTINUE
843  END IF
844 *
845 * Undo scaling if necessary
846 *
847  130 CONTINUE
848 *
849  IF( ilascl ) THEN
850  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
851  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
852  END IF
853 *
854  IF( ilbscl ) THEN
855  CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
856  END IF
857 *
858  work( 1 ) = maxwrk
859  RETURN
860 *
861 * End of SGGEVX
862 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
Definition: sggbak.f:147
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
Definition: sggbal.f:177
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:295
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:145
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
Definition: shgeqz.f:304
subroutine stgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
STGSNA
Definition: stgsna.f:381
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:128
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:168
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
Definition: sgghrd.f:207
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: