LAPACK  3.8.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.
Date
April 2012
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 393 of file sggevx.f.

393 *
394 * -- LAPACK driver routine (version 3.7.0) --
395 * -- LAPACK is a software package provided by Univ. of Tennessee, --
396 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
397 * April 2012
398 *
399 * .. Scalar Arguments ..
400  CHARACTER balanc, jobvl, jobvr, sense
401  INTEGER ihi, ilo, info, lda, ldb, ldvl, ldvr, lwork, n
402  REAL abnrm, bbnrm
403 * ..
404 * .. Array Arguments ..
405  LOGICAL bwork( * )
406  INTEGER iwork( * )
407  REAL a( lda, * ), alphai( * ), alphar( * ),
408  $ b( ldb, * ), beta( * ), lscale( * ),
409  $ rconde( * ), rcondv( * ), rscale( * ),
410  $ vl( ldvl, * ), vr( ldvr, * ), work( * )
411 * ..
412 *
413 * =====================================================================
414 *
415 * .. Parameters ..
416  REAL zero, one
417  parameter( zero = 0.0e+0, one = 1.0e+0 )
418 * ..
419 * .. Local Scalars ..
420  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery, noscl,
421  $ pair, wantsb, wantse, wantsn, wantsv
422  CHARACTER chtemp
423  INTEGER i, icols, ierr, ijobvl, ijobvr, in, irows,
424  $ itau, iwrk, iwrk1, j, jc, jr, m, maxwrk,
425  $ minwrk, mm
426  REAL anrm, anrmto, bignum, bnrm, bnrmto, eps,
427  $ smlnum, temp
428 * ..
429 * .. Local Arrays ..
430  LOGICAL ldumma( 1 )
431 * ..
432 * .. External Subroutines ..
433  EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
435  $ stgsna, xerbla
436 * ..
437 * .. External Functions ..
438  LOGICAL lsame
439  INTEGER ilaenv
440  REAL slamch, slange
441  EXTERNAL lsame, ilaenv, slamch, slange
442 * ..
443 * .. Intrinsic Functions ..
444  INTRINSIC abs, max, sqrt
445 * ..
446 * .. Executable Statements ..
447 *
448 * Decode the input arguments
449 *
450  IF( lsame( jobvl, 'N' ) ) THEN
451  ijobvl = 1
452  ilvl = .false.
453  ELSE IF( lsame( jobvl, 'V' ) ) THEN
454  ijobvl = 2
455  ilvl = .true.
456  ELSE
457  ijobvl = -1
458  ilvl = .false.
459  END IF
460 *
461  IF( lsame( jobvr, 'N' ) ) THEN
462  ijobvr = 1
463  ilvr = .false.
464  ELSE IF( lsame( jobvr, 'V' ) ) THEN
465  ijobvr = 2
466  ilvr = .true.
467  ELSE
468  ijobvr = -1
469  ilvr = .false.
470  END IF
471  ilv = ilvl .OR. ilvr
472 *
473  noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
474  wantsn = lsame( sense, 'N' )
475  wantse = lsame( sense, 'E' )
476  wantsv = lsame( sense, 'V' )
477  wantsb = lsame( sense, 'B' )
478 *
479 * Test the input arguments
480 *
481  info = 0
482  lquery = ( lwork.EQ.-1 )
483  IF( .NOT.( noscl .OR. lsame( balanc, 'S' ) .OR.
484  $ lsame( balanc, 'B' ) ) ) THEN
485  info = -1
486  ELSE IF( ijobvl.LE.0 ) THEN
487  info = -2
488  ELSE IF( ijobvr.LE.0 ) THEN
489  info = -3
490  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
491  $ THEN
492  info = -4
493  ELSE IF( n.LT.0 ) THEN
494  info = -5
495  ELSE IF( lda.LT.max( 1, n ) ) THEN
496  info = -7
497  ELSE IF( ldb.LT.max( 1, n ) ) THEN
498  info = -9
499  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
500  info = -14
501  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
502  info = -16
503  END IF
504 *
505 * Compute workspace
506 * (Note: Comments in the code beginning "Workspace:" describe the
507 * minimal amount of workspace needed at that point in the code,
508 * as well as the preferred amount for good performance.
509 * NB refers to the optimal block size for the immediately
510 * following subroutine, as returned by ILAENV. The workspace is
511 * computed assuming ILO = 1 and IHI = N, the worst case.)
512 *
513  IF( info.EQ.0 ) THEN
514  IF( n.EQ.0 ) THEN
515  minwrk = 1
516  maxwrk = 1
517  ELSE
518  IF( noscl .AND. .NOT.ilv ) THEN
519  minwrk = 2*n
520  ELSE
521  minwrk = 6*n
522  END IF
523  IF( wantse ) THEN
524  minwrk = 10*n
525  ELSE IF( wantsv .OR. wantsb ) THEN
526  minwrk = 2*n*( n + 4 ) + 16
527  END IF
528  maxwrk = minwrk
529  maxwrk = max( maxwrk,
530  $ n + n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) )
531  maxwrk = max( maxwrk,
532  $ n + n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) )
533  IF( ilvl ) THEN
534  maxwrk = max( maxwrk, n +
535  $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, 0 ) )
536  END IF
537  END IF
538  work( 1 ) = maxwrk
539 *
540  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
541  info = -26
542  END IF
543  END IF
544 *
545  IF( info.NE.0 ) THEN
546  CALL xerbla( 'SGGEVX', -info )
547  RETURN
548  ELSE IF( lquery ) THEN
549  RETURN
550  END IF
551 *
552 * Quick return if possible
553 *
554  IF( n.EQ.0 )
555  $ RETURN
556 *
557 *
558 * Get machine constants
559 *
560  eps = slamch( 'P' )
561  smlnum = slamch( 'S' )
562  bignum = one / smlnum
563  CALL slabad( smlnum, bignum )
564  smlnum = sqrt( smlnum ) / eps
565  bignum = one / smlnum
566 *
567 * Scale A if max element outside range [SMLNUM,BIGNUM]
568 *
569  anrm = slange( 'M', n, n, a, lda, work )
570  ilascl = .false.
571  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
572  anrmto = smlnum
573  ilascl = .true.
574  ELSE IF( anrm.GT.bignum ) THEN
575  anrmto = bignum
576  ilascl = .true.
577  END IF
578  IF( ilascl )
579  $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
580 *
581 * Scale B if max element outside range [SMLNUM,BIGNUM]
582 *
583  bnrm = slange( 'M', n, n, b, ldb, work )
584  ilbscl = .false.
585  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
586  bnrmto = smlnum
587  ilbscl = .true.
588  ELSE IF( bnrm.GT.bignum ) THEN
589  bnrmto = bignum
590  ilbscl = .true.
591  END IF
592  IF( ilbscl )
593  $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
594 *
595 * Permute and/or balance the matrix pair (A,B)
596 * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
597 *
598  CALL sggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
599  $ work, ierr )
600 *
601 * Compute ABNRM and BBNRM
602 *
603  abnrm = slange( '1', n, n, a, lda, work( 1 ) )
604  IF( ilascl ) THEN
605  work( 1 ) = abnrm
606  CALL slascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
607  $ ierr )
608  abnrm = work( 1 )
609  END IF
610 *
611  bbnrm = slange( '1', n, n, b, ldb, work( 1 ) )
612  IF( ilbscl ) THEN
613  work( 1 ) = bbnrm
614  CALL slascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
615  $ ierr )
616  bbnrm = work( 1 )
617  END IF
618 *
619 * Reduce B to triangular form (QR decomposition of B)
620 * (Workspace: need N, prefer N*NB )
621 *
622  irows = ihi + 1 - ilo
623  IF( ilv .OR. .NOT.wantsn ) THEN
624  icols = n + 1 - ilo
625  ELSE
626  icols = irows
627  END IF
628  itau = 1
629  iwrk = itau + irows
630  CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
631  $ work( iwrk ), lwork+1-iwrk, ierr )
632 *
633 * Apply the orthogonal transformation to A
634 * (Workspace: need N, prefer N*NB)
635 *
636  CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
637  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
638  $ lwork+1-iwrk, ierr )
639 *
640 * Initialize VL and/or VR
641 * (Workspace: need N, prefer N*NB)
642 *
643  IF( ilvl ) THEN
644  CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
645  IF( irows.GT.1 ) THEN
646  CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
647  $ vl( ilo+1, ilo ), ldvl )
648  END IF
649  CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
650  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
651  END IF
652 *
653  IF( ilvr )
654  $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
655 *
656 * Reduce to generalized Hessenberg form
657 * (Workspace: none needed)
658 *
659  IF( ilv .OR. .NOT.wantsn ) THEN
660 *
661 * Eigenvectors requested -- work on whole matrix.
662 *
663  CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
664  $ ldvl, vr, ldvr, ierr )
665  ELSE
666  CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
667  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
668  END IF
669 *
670 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
671 * Schur forms and Schur vectors)
672 * (Workspace: need N)
673 *
674  IF( ilv .OR. .NOT.wantsn ) THEN
675  chtemp = 'S'
676  ELSE
677  chtemp = 'E'
678  END IF
679 *
680  CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
681  $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
682  $ lwork, ierr )
683  IF( ierr.NE.0 ) THEN
684  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
685  info = ierr
686  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
687  info = ierr - n
688  ELSE
689  info = n + 1
690  END IF
691  GO TO 130
692  END IF
693 *
694 * Compute Eigenvectors and estimate condition numbers if desired
695 * (Workspace: STGEVC: need 6*N
696 * STGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
697 * need N otherwise )
698 *
699  IF( ilv .OR. .NOT.wantsn ) THEN
700  IF( ilv ) THEN
701  IF( ilvl ) THEN
702  IF( ilvr ) THEN
703  chtemp = 'B'
704  ELSE
705  chtemp = 'L'
706  END IF
707  ELSE
708  chtemp = 'R'
709  END IF
710 *
711  CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
712  $ ldvl, vr, ldvr, n, in, work, ierr )
713  IF( ierr.NE.0 ) THEN
714  info = n + 2
715  GO TO 130
716  END IF
717  END IF
718 *
719  IF( .NOT.wantsn ) THEN
720 *
721 * compute eigenvectors (STGEVC) and estimate condition
722 * numbers (STGSNA). Note that the definition of the condition
723 * number is not invariant under transformation (u,v) to
724 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
725 * Schur form (S,T), Q and Z are orthogonal matrices. In order
726 * to avoid using extra 2*N*N workspace, we have to recalculate
727 * eigenvectors and estimate one condition numbers at a time.
728 *
729  pair = .false.
730  DO 20 i = 1, n
731 *
732  IF( pair ) THEN
733  pair = .false.
734  GO TO 20
735  END IF
736  mm = 1
737  IF( i.LT.n ) THEN
738  IF( a( i+1, i ).NE.zero ) THEN
739  pair = .true.
740  mm = 2
741  END IF
742  END IF
743 *
744  DO 10 j = 1, n
745  bwork( j ) = .false.
746  10 CONTINUE
747  IF( mm.EQ.1 ) THEN
748  bwork( i ) = .true.
749  ELSE IF( mm.EQ.2 ) THEN
750  bwork( i ) = .true.
751  bwork( i+1 ) = .true.
752  END IF
753 *
754  iwrk = mm*n + 1
755  iwrk1 = iwrk + mm*n
756 *
757 * Compute a pair of left and right eigenvectors.
758 * (compute workspace: need up to 4*N + 6*N)
759 *
760  IF( wantse .OR. wantsb ) THEN
761  CALL stgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
762  $ work( 1 ), n, work( iwrk ), n, mm, m,
763  $ work( iwrk1 ), ierr )
764  IF( ierr.NE.0 ) THEN
765  info = n + 2
766  GO TO 130
767  END IF
768  END IF
769 *
770  CALL stgsna( sense, 'S', bwork, n, a, lda, b, ldb,
771  $ work( 1 ), n, work( iwrk ), n, rconde( i ),
772  $ rcondv( i ), mm, m, work( iwrk1 ),
773  $ lwork-iwrk1+1, iwork, ierr )
774 *
775  20 CONTINUE
776  END IF
777  END IF
778 *
779 * Undo balancing on VL and VR and normalization
780 * (Workspace: none needed)
781 *
782  IF( ilvl ) THEN
783  CALL sggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
784  $ ldvl, ierr )
785 *
786  DO 70 jc = 1, n
787  IF( alphai( jc ).LT.zero )
788  $ GO TO 70
789  temp = zero
790  IF( alphai( jc ).EQ.zero ) THEN
791  DO 30 jr = 1, n
792  temp = max( temp, abs( vl( jr, jc ) ) )
793  30 CONTINUE
794  ELSE
795  DO 40 jr = 1, n
796  temp = max( temp, abs( vl( jr, jc ) )+
797  $ abs( vl( jr, jc+1 ) ) )
798  40 CONTINUE
799  END IF
800  IF( temp.LT.smlnum )
801  $ GO TO 70
802  temp = one / temp
803  IF( alphai( jc ).EQ.zero ) THEN
804  DO 50 jr = 1, n
805  vl( jr, jc ) = vl( jr, jc )*temp
806  50 CONTINUE
807  ELSE
808  DO 60 jr = 1, n
809  vl( jr, jc ) = vl( jr, jc )*temp
810  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
811  60 CONTINUE
812  END IF
813  70 CONTINUE
814  END IF
815  IF( ilvr ) THEN
816  CALL sggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
817  $ ldvr, ierr )
818  DO 120 jc = 1, n
819  IF( alphai( jc ).LT.zero )
820  $ GO TO 120
821  temp = zero
822  IF( alphai( jc ).EQ.zero ) THEN
823  DO 80 jr = 1, n
824  temp = max( temp, abs( vr( jr, jc ) ) )
825  80 CONTINUE
826  ELSE
827  DO 90 jr = 1, n
828  temp = max( temp, abs( vr( jr, jc ) )+
829  $ abs( vr( jr, jc+1 ) ) )
830  90 CONTINUE
831  END IF
832  IF( temp.LT.smlnum )
833  $ GO TO 120
834  temp = one / temp
835  IF( alphai( jc ).EQ.zero ) THEN
836  DO 100 jr = 1, n
837  vr( jr, jc ) = vr( jr, jc )*temp
838  100 CONTINUE
839  ELSE
840  DO 110 jr = 1, n
841  vr( jr, jc ) = vr( jr, jc )*temp
842  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
843  110 CONTINUE
844  END IF
845  120 CONTINUE
846  END IF
847 *
848 * Undo scaling if necessary
849 *
850  130 CONTINUE
851 *
852  IF( ilascl ) THEN
853  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
854  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
855  END IF
856 *
857  IF( ilbscl ) THEN
858  CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
859  END IF
860 *
861  work( 1 ) = maxwrk
862  RETURN
863 *
864 * End of SGGEVX
865 *
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:297
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:306
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:170
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:383
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:138
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
Definition: sggbak.f:149
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:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:112
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:145
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:130
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
Definition: sggbal.f:179
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
Definition: sgghrd.f:209
Here is the call graph for this function:
Here is the caller graph for this function: