LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ sggesx()

subroutine sggesx ( character  JOBVSL,
character  JOBVSR,
character  SORT,
external  SELCTG,
character  SENSE,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
integer  SDIM,
real, dimension( * )  ALPHAR,
real, dimension( * )  ALPHAI,
real, dimension( * )  BETA,
real, dimension( ldvsl, * )  VSL,
integer  LDVSL,
real, dimension( ldvsr, * )  VSR,
integer  LDVSR,
real, dimension( 2 )  RCONDE,
real, dimension( 2 )  RCONDV,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

SGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

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

Purpose:
 SGGESX computes for a pair of N-by-N real nonsymmetric matrices
 (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
 optionally, the left and/or right matrices of Schur vectors (VSL and
 VSR).  This gives the generalized Schur factorization

      (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )

 Optionally, it also orders the eigenvalues so that a selected cluster
 of eigenvalues appears in the leading diagonal blocks of the upper
 quasi-triangular matrix S and the upper triangular matrix T; computes
 a reciprocal condition number for the average of the selected
 eigenvalues (RCONDE); and computes a reciprocal condition number for
 the right and left deflating subspaces corresponding to the selected
 eigenvalues (RCONDV). The leading columns of VSL and VSR then form
 an orthonormal basis for the corresponding left and right eigenspaces
 (deflating subspaces).

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
 or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
 usually represented as the pair (alpha,beta), as there is a
 reasonable interpretation for beta=0 or for both being zero.

 A pair of matrices (S,T) is in generalized real Schur form if T is
 upper triangular with non-negative diagonal and S is block upper
 triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
 to real generalized eigenvalues, while 2-by-2 blocks of S will be
 "standardized" by making the corresponding elements of T have the
 form:
         [  a  0  ]
         [  0  b  ]

 and the pair of corresponding 2-by-2 blocks in S and T will have a
 complex conjugate pair of generalized eigenvalues.
Parameters
[in]JOBVSL
          JOBVSL is CHARACTER*1
          = 'N':  do not compute the left Schur vectors;
          = 'V':  compute the left Schur vectors.
[in]JOBVSR
          JOBVSR is CHARACTER*1
          = 'N':  do not compute the right Schur vectors;
          = 'V':  compute the right Schur vectors.
[in]SORT
          SORT is CHARACTER*1
          Specifies whether or not to order the eigenvalues on the
          diagonal of the generalized Schur form.
          = 'N':  Eigenvalues are not ordered;
          = 'S':  Eigenvalues are ordered (see SELCTG).
[in]SELCTG
          SELCTG is a LOGICAL FUNCTION of three REAL arguments
          SELCTG must be declared EXTERNAL in the calling subroutine.
          If SORT = 'N', SELCTG is not referenced.
          If SORT = 'S', SELCTG is used to select eigenvalues to sort
          to the top left of the Schur form.
          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
          one of a complex conjugate pair of eigenvalues is selected,
          then both complex eigenvalues are selected.
          Note that a selected complex eigenvalue may no longer satisfy
          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
          since ordering may change the value of complex eigenvalues
          (especially if the eigenvalue is ill-conditioned), in this
          case INFO is set to N+3.
[in]SENSE
          SENSE is CHARACTER*1
          Determines which reciprocal condition numbers are computed.
          = 'N' : None are computed;
          = 'E' : Computed for average of selected eigenvalues only;
          = 'V' : Computed for selected deflating subspaces only;
          = 'B' : Computed for both.
          If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
[in]N
          N is INTEGER
          The order of the matrices A, B, VSL, and VSR.  N >= 0.
[in,out]A
          A is REAL array, dimension (LDA, N)
          On entry, the first of the pair of matrices.
          On exit, A has been overwritten by its generalized Schur
          form S.
[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 second of the pair of matrices.
          On exit, B has been overwritten by its generalized Schur
          form T.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]SDIM
          SDIM is INTEGER
          If SORT = 'N', SDIM = 0.
          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
          for which SELCTG is true.  (Complex conjugate pairs for which
          SELCTG is true for either eigenvalue count as 2.)
[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.  ALPHAR(j) + ALPHAI(j)*i
          and BETA(j),j=1,...,N  are the diagonals of the complex Schur
          form (S,T) that would result if the 2-by-2 diagonal blocks of
          the real Schur form of (A,B) were further reduced to
          triangular form using 2-by-2 complex unitary transformations.
          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.
          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]VSL
          VSL is REAL array, dimension (LDVSL,N)
          If JOBVSL = 'V', VSL will contain the left Schur vectors.
          Not referenced if JOBVSL = 'N'.
[in]LDVSL
          LDVSL is INTEGER
          The leading dimension of the matrix VSL. LDVSL >=1, and
          if JOBVSL = 'V', LDVSL >= N.
[out]VSR
          VSR is REAL array, dimension (LDVSR,N)
          If JOBVSR = 'V', VSR will contain the right Schur vectors.
          Not referenced if JOBVSR = 'N'.
[in]LDVSR
          LDVSR is INTEGER
          The leading dimension of the matrix VSR. LDVSR >= 1, and
          if JOBVSR = 'V', LDVSR >= N.
[out]RCONDE
          RCONDE is REAL array, dimension ( 2 )
          If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
          reciprocal condition numbers for the average of the selected
          eigenvalues.
          Not referenced if SENSE = 'N' or 'V'.
[out]RCONDV
          RCONDV is REAL array, dimension ( 2 )
          If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
          reciprocal condition numbers for the selected deflating
          subspaces.
          Not referenced if SENSE = 'N' or 'E'.
[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.
          If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
          LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
          LWORK >= max( 8*N, 6*N+16 ).
          Note that 2*SDIM*(N-SDIM) <= N*N/2.
          Note also that an error is only returned if
          LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
          this may not be large enough.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the bound on the optimal size of the WORK
          array and the minimum size of the IWORK array, returns these
          values as the first entries of the WORK and IWORK arrays, and
          no error message related to LWORK or LIWORK is issued by
          XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.
          If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
          LIWORK >= N+6.

          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the bound on the optimal size of the
          WORK array and the minimum size of the IWORK array, returns
          these values as the first entries of the WORK and IWORK
          arrays, and no error message related to LWORK or LIWORK is
          issued by XERBLA.
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
          Not referenced if SORT = 'N'.
[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.  (A,B) are not in Schur
                form, 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: after reordering, roundoff changed values of
                      some complex eigenvalues so that leading
                      eigenvalues in the Generalized Schur form no
                      longer satisfy SELCTG=.TRUE.  This could also
                      be caused due to scaling.
                =N+3: reordering failed in STGSEN.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017
Further Details:
  An approximate (asymptotic) bound on the average absolute error of
  the selected eigenvalues is

       EPS * norm((A, B)) / RCONDE( 1 ).

  An approximate (asymptotic) bound on the maximum angular error in
  the computed deflating subspaces is

       EPS * norm((A, B)) / RCONDV( 2 ).

  See LAPACK User's Guide, section 4.11 for more information.

Definition at line 367 of file sggesx.f.

367 *
368 * -- LAPACK driver routine (version 3.7.1) --
369 * -- LAPACK is a software package provided by Univ. of Tennessee, --
370 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
371 * June 2017
372 *
373 * .. Scalar Arguments ..
374  CHARACTER jobvsl, jobvsr, sense, sort
375  INTEGER info, lda, ldb, ldvsl, ldvsr, liwork, lwork, n,
376  $ sdim
377 * ..
378 * .. Array Arguments ..
379  LOGICAL bwork( * )
380  INTEGER iwork( * )
381  REAL a( lda, * ), alphai( * ), alphar( * ),
382  $ b( ldb, * ), beta( * ), rconde( 2 ),
383  $ rcondv( 2 ), vsl( ldvsl, * ), vsr( ldvsr, * ),
384  $ work( * )
385 * ..
386 * .. Function Arguments ..
387  LOGICAL selctg
388  EXTERNAL selctg
389 * ..
390 *
391 * =====================================================================
392 *
393 * .. Parameters ..
394  REAL zero, one
395  parameter( zero = 0.0e+0, one = 1.0e+0 )
396 * ..
397 * .. Local Scalars ..
398  LOGICAL cursl, ilascl, ilbscl, ilvsl, ilvsr, lastsl,
399  $ lquery, lst2sl, wantsb, wantse, wantsn, wantst,
400  $ wantsv
401  INTEGER i, icols, ierr, ihi, ijob, ijobvl, ijobvr,
402  $ ileft, ilo, ip, iright, irows, itau, iwrk,
403  $ liwmin, lwrk, maxwrk, minwrk
404  REAL anrm, anrmto, bignum, bnrm, bnrmto, eps, pl,
405  $ pr, safmax, safmin, smlnum
406 * ..
407 * .. Local Arrays ..
408  REAL dif( 2 )
409 * ..
410 * .. External Subroutines ..
411  EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
413  $ xerbla
414 * ..
415 * .. External Functions ..
416  LOGICAL lsame
417  INTEGER ilaenv
418  REAL slamch, slange
419  EXTERNAL lsame, ilaenv, slamch, slange
420 * ..
421 * .. Intrinsic Functions ..
422  INTRINSIC abs, max, sqrt
423 * ..
424 * .. Executable Statements ..
425 *
426 * Decode the input arguments
427 *
428  IF( lsame( jobvsl, 'N' ) ) THEN
429  ijobvl = 1
430  ilvsl = .false.
431  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
432  ijobvl = 2
433  ilvsl = .true.
434  ELSE
435  ijobvl = -1
436  ilvsl = .false.
437  END IF
438 *
439  IF( lsame( jobvsr, 'N' ) ) THEN
440  ijobvr = 1
441  ilvsr = .false.
442  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
443  ijobvr = 2
444  ilvsr = .true.
445  ELSE
446  ijobvr = -1
447  ilvsr = .false.
448  END IF
449 *
450  wantst = lsame( sort, 'S' )
451  wantsn = lsame( sense, 'N' )
452  wantse = lsame( sense, 'E' )
453  wantsv = lsame( sense, 'V' )
454  wantsb = lsame( sense, 'B' )
455  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
456  IF( wantsn ) THEN
457  ijob = 0
458  ELSE IF( wantse ) THEN
459  ijob = 1
460  ELSE IF( wantsv ) THEN
461  ijob = 2
462  ELSE IF( wantsb ) THEN
463  ijob = 4
464  END IF
465 *
466 * Test the input arguments
467 *
468  info = 0
469  IF( ijobvl.LE.0 ) THEN
470  info = -1
471  ELSE IF( ijobvr.LE.0 ) THEN
472  info = -2
473  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
474  info = -3
475  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
476  $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
477  info = -5
478  ELSE IF( n.LT.0 ) THEN
479  info = -6
480  ELSE IF( lda.LT.max( 1, n ) ) THEN
481  info = -8
482  ELSE IF( ldb.LT.max( 1, n ) ) THEN
483  info = -10
484  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
485  info = -16
486  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
487  info = -18
488  END IF
489 *
490 * Compute workspace
491 * (Note: Comments in the code beginning "Workspace:" describe the
492 * minimal amount of workspace needed at that point in the code,
493 * as well as the preferred amount for good performance.
494 * NB refers to the optimal block size for the immediately
495 * following subroutine, as returned by ILAENV.)
496 *
497  IF( info.EQ.0 ) THEN
498  IF( n.GT.0) THEN
499  minwrk = max( 8*n, 6*n + 16 )
500  maxwrk = minwrk - n +
501  $ n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 )
502  maxwrk = max( maxwrk, minwrk - n +
503  $ n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, -1 ) )
504  IF( ilvsl ) THEN
505  maxwrk = max( maxwrk, minwrk - n +
506  $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, -1 ) )
507  END IF
508  lwrk = maxwrk
509  IF( ijob.GE.1 )
510  $ lwrk = max( lwrk, n*n/2 )
511  ELSE
512  minwrk = 1
513  maxwrk = 1
514  lwrk = 1
515  END IF
516  work( 1 ) = lwrk
517  IF( wantsn .OR. n.EQ.0 ) THEN
518  liwmin = 1
519  ELSE
520  liwmin = n + 6
521  END IF
522  iwork( 1 ) = liwmin
523 *
524  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
525  info = -22
526  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
527  info = -24
528  END IF
529  END IF
530 *
531  IF( info.NE.0 ) THEN
532  CALL xerbla( 'SGGESX', -info )
533  RETURN
534  ELSE IF (lquery) THEN
535  RETURN
536  END IF
537 *
538 * Quick return if possible
539 *
540  IF( n.EQ.0 ) THEN
541  sdim = 0
542  RETURN
543  END IF
544 *
545 * Get machine constants
546 *
547  eps = slamch( 'P' )
548  safmin = slamch( 'S' )
549  safmax = one / safmin
550  CALL slabad( safmin, safmax )
551  smlnum = sqrt( safmin ) / eps
552  bignum = one / smlnum
553 *
554 * Scale A if max element outside range [SMLNUM,BIGNUM]
555 *
556  anrm = slange( 'M', n, n, a, lda, work )
557  ilascl = .false.
558  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
559  anrmto = smlnum
560  ilascl = .true.
561  ELSE IF( anrm.GT.bignum ) THEN
562  anrmto = bignum
563  ilascl = .true.
564  END IF
565  IF( ilascl )
566  $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
567 *
568 * Scale B if max element outside range [SMLNUM,BIGNUM]
569 *
570  bnrm = slange( 'M', n, n, b, ldb, work )
571  ilbscl = .false.
572  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
573  bnrmto = smlnum
574  ilbscl = .true.
575  ELSE IF( bnrm.GT.bignum ) THEN
576  bnrmto = bignum
577  ilbscl = .true.
578  END IF
579  IF( ilbscl )
580  $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
581 *
582 * Permute the matrix to make it more nearly triangular
583 * (Workspace: need 6*N + 2*N for permutation parameters)
584 *
585  ileft = 1
586  iright = n + 1
587  iwrk = iright + n
588  CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
589  $ work( iright ), work( iwrk ), ierr )
590 *
591 * Reduce B to triangular form (QR decomposition of B)
592 * (Workspace: need N, prefer N*NB)
593 *
594  irows = ihi + 1 - ilo
595  icols = n + 1 - ilo
596  itau = iwrk
597  iwrk = itau + irows
598  CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
599  $ work( iwrk ), lwork+1-iwrk, ierr )
600 *
601 * Apply the orthogonal transformation to matrix A
602 * (Workspace: need N, prefer N*NB)
603 *
604  CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
605  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
606  $ lwork+1-iwrk, ierr )
607 *
608 * Initialize VSL
609 * (Workspace: need N, prefer N*NB)
610 *
611  IF( ilvsl ) THEN
612  CALL slaset( 'Full', n, n, zero, one, vsl, ldvsl )
613  IF( irows.GT.1 ) THEN
614  CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
615  $ vsl( ilo+1, ilo ), ldvsl )
616  END IF
617  CALL sorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
618  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
619  END IF
620 *
621 * Initialize VSR
622 *
623  IF( ilvsr )
624  $ CALL slaset( 'Full', n, n, zero, one, vsr, ldvsr )
625 *
626 * Reduce to generalized Hessenberg form
627 * (Workspace: none needed)
628 *
629  CALL sgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
630  $ ldvsl, vsr, ldvsr, ierr )
631 *
632  sdim = 0
633 *
634 * Perform QZ algorithm, computing Schur vectors if desired
635 * (Workspace: need N)
636 *
637  iwrk = itau
638  CALL shgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
639  $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
640  $ work( iwrk ), lwork+1-iwrk, ierr )
641  IF( ierr.NE.0 ) THEN
642  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
643  info = ierr
644  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
645  info = ierr - n
646  ELSE
647  info = n + 1
648  END IF
649  GO TO 50
650  END IF
651 *
652 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
653 * condition number(s)
654 * (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
655 * otherwise, need 8*(N+1) )
656 *
657  IF( wantst ) THEN
658 *
659 * Undo scaling on eigenvalues before SELCTGing
660 *
661  IF( ilascl ) THEN
662  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
663  $ ierr )
664  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
665  $ ierr )
666  END IF
667  IF( ilbscl )
668  $ CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
669 *
670 * Select eigenvalues
671 *
672  DO 10 i = 1, n
673  bwork( i ) = selctg( alphar( i ), alphai( i ), beta( i ) )
674  10 CONTINUE
675 *
676 * Reorder eigenvalues, transform Generalized Schur vectors, and
677 * compute reciprocal condition numbers
678 *
679  CALL stgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
680  $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
681  $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
682  $ iwork, liwork, ierr )
683 *
684  IF( ijob.GE.1 )
685  $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
686  IF( ierr.EQ.-22 ) THEN
687 *
688 * not enough real workspace
689 *
690  info = -22
691  ELSE
692  IF( ijob.EQ.1 .OR. ijob.EQ.4 ) THEN
693  rconde( 1 ) = pl
694  rconde( 2 ) = pr
695  END IF
696  IF( ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
697  rcondv( 1 ) = dif( 1 )
698  rcondv( 2 ) = dif( 2 )
699  END IF
700  IF( ierr.EQ.1 )
701  $ info = n + 3
702  END IF
703 *
704  END IF
705 *
706 * Apply permutation to VSL and VSR
707 * (Workspace: none needed)
708 *
709  IF( ilvsl )
710  $ CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
711  $ work( iright ), n, vsl, ldvsl, ierr )
712 *
713  IF( ilvsr )
714  $ CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
715  $ work( iright ), n, vsr, ldvsr, ierr )
716 *
717 * Check if unscaling would cause over/underflow, if so, rescale
718 * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
719 * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
720 *
721  IF( ilascl ) THEN
722  DO 20 i = 1, n
723  IF( alphai( i ).NE.zero ) THEN
724  IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
725  $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) )
726  $ THEN
727  work( 1 ) = abs( a( i, i ) / alphar( i ) )
728  beta( i ) = beta( i )*work( 1 )
729  alphar( i ) = alphar( i )*work( 1 )
730  alphai( i ) = alphai( i )*work( 1 )
731  ELSE IF( ( alphai( i ) / safmax ).GT.( anrmto / anrm )
732  $ .OR. ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
733  $ THEN
734  work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
735  beta( i ) = beta( i )*work( 1 )
736  alphar( i ) = alphar( i )*work( 1 )
737  alphai( i ) = alphai( i )*work( 1 )
738  END IF
739  END IF
740  20 CONTINUE
741  END IF
742 *
743  IF( ilbscl ) THEN
744  DO 25 i = 1, n
745  IF( alphai( i ).NE.zero ) THEN
746  IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
747  $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) ) THEN
748  work( 1 ) = abs( b( i, i ) / beta( i ) )
749  beta( i ) = beta( i )*work( 1 )
750  alphar( i ) = alphar( i )*work( 1 )
751  alphai( i ) = alphai( i )*work( 1 )
752  END IF
753  END IF
754  25 CONTINUE
755  END IF
756 *
757 * Undo scaling
758 *
759  IF( ilascl ) THEN
760  CALL slascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
761  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
762  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
763  END IF
764 *
765  IF( ilbscl ) THEN
766  CALL slascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
767  CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
768  END IF
769 *
770  IF( wantst ) THEN
771 *
772 * Check if reordering is correct
773 *
774  lastsl = .true.
775  lst2sl = .true.
776  sdim = 0
777  ip = 0
778  DO 40 i = 1, n
779  cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
780  IF( alphai( i ).EQ.zero ) THEN
781  IF( cursl )
782  $ sdim = sdim + 1
783  ip = 0
784  IF( cursl .AND. .NOT.lastsl )
785  $ info = n + 2
786  ELSE
787  IF( ip.EQ.1 ) THEN
788 *
789 * Last eigenvalue of conjugate pair
790 *
791  cursl = cursl .OR. lastsl
792  lastsl = cursl
793  IF( cursl )
794  $ sdim = sdim + 2
795  ip = -1
796  IF( cursl .AND. .NOT.lst2sl )
797  $ info = n + 2
798  ELSE
799 *
800 * First eigenvalue of conjugate pair
801 *
802  ip = 1
803  END IF
804  END IF
805  lst2sl = lastsl
806  lastsl = cursl
807  40 CONTINUE
808 *
809  END IF
810 *
811  50 CONTINUE
812 *
813  work( 1 ) = maxwrk
814  iwork( 1 ) = liwmin
815 *
816  RETURN
817 *
818 * End of SGGESX
819 *
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 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 stgsen(IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK, LWORK, IWORK, LIWORK, INFO)
STGSEN
Definition: stgsen.f:453
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: