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

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