Move LAPACK trunk into position.

1 files changed, 558 insertions, 0 deletions

diff --git a/SRC/sgges.f b/SRC/sgges.f
new file mode 100644
index 00000000..00f16a5e
--- /dev/null
+++ b/SRC/sgges.f
@@ -0,0 +1,558 @@
+      SUBROUTINE SGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
+     $                  SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
+     $                  LDVSR, WORK, LWORK, BWORK, INFO )
+*
+*  -- LAPACK driver routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          JOBVSL, JOBVSR, SORT
+      INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
+*     ..
+*     .. Array Arguments ..
+      LOGICAL            BWORK( * )
+      REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+     $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
+     $                   VSR( LDVSR, * ), WORK( * )
+*     ..
+*     .. Function Arguments ..
+      LOGICAL            SELCTG
+      EXTERNAL           SELCTG
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
+*  the generalized eigenvalues, the generalized real Schur form (S,T),
+*  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.The
+*  leading columns of VSL and VSR then form an orthonormal basis for the
+*  corresponding left and right eigenspaces (deflating subspaces).
+*
+*  (If only the generalized eigenvalues are needed, use the driver
+*  SGGEV instead, which is faster.)
+*
+*  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 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.
+*
+*
+*  Arguments
+*  =========
+*
+*  JOBVSL  (input) CHARACTER*1
+*          = 'N':  do not compute the left Schur vectors;
+*          = 'V':  compute the left Schur vectors.
+*
+*  JOBVSR  (input) CHARACTER*1
+*          = 'N':  do not compute the right Schur vectors;
+*          = 'V':  compute the right Schur vectors.
+*
+*  SORT    (input) 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);
+*
+*  SELCTG  (external procedure) 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 in the ill-conditioned case, a selected complex
+*          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
+*          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
+*          in this case.
+*
+*  N       (input) INTEGER
+*          The order of the matrices A, B, VSL, and VSR.  N >= 0.
+*
+*  A       (input/output) 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.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of A.  LDA >= max(1,N).
+*
+*  B       (input/output) 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.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of B.  LDB >= max(1,N).
+*
+*  SDIM    (output) 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.)
+*
+*  ALPHAR  (output) REAL array, dimension (N)
+*  ALPHAI  (output) REAL array, dimension (N)
+*  BETA    (output) 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).
+*
+*  VSL     (output) REAL array, dimension (LDVSL,N)
+*          If JOBVSL = 'V', VSL will contain the left Schur vectors.
+*          Not referenced if JOBVSL = 'N'.
+*
+*  LDVSL   (input) INTEGER
+*          The leading dimension of the matrix VSL. LDVSL >=1, and
+*          if JOBVSL = 'V', LDVSL >= N.
+*
+*  VSR     (output) REAL array, dimension (LDVSR,N)
+*          If JOBVSR = 'V', VSR will contain the right Schur vectors.
+*          Not referenced if JOBVSR = 'N'.
+*
+*  LDVSR   (input) INTEGER
+*          The leading dimension of the matrix VSR. LDVSR >= 1, and
+*          if JOBVSR = 'V', LDVSR >= N.
+*
+*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
+*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK.
+*          If N = 0, LWORK >= 1, else LWORK >= max(8*N,6*N+16).
+*          For good performance , LWORK must generally be larger.
+*
+*          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.
+*
+*  BWORK   (workspace) LOGICAL array, dimension (N)
+*          Not referenced if SORT = 'N'.
+*
+*  INFO    (output) 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.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
+     $                   LQUERY, LST2SL, WANTST
+      INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
+     $                   ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
+     $                   MINWRK
+      REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
+     $                   PVSR, SAFMAX, SAFMIN, SMLNUM
+*     ..
+*     .. Local Arrays ..
+      INTEGER            IDUM( 1 )
+      REAL               DIF( 2 )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHRD, SHGEQZ, SLABAD,
+     $                   SLACPY, SLASCL, SLASET, SORGQR, SORMQR, STGSEN,
+     $                   XERBLA
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAENV
+      REAL               SLAMCH, SLANGE
+      EXTERNAL           LSAME, ILAENV, SLAMCH, SLANGE
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, SQRT
+*     ..
+*     .. Executable Statements ..
+*
+*     Decode the input arguments
+*
+      IF( LSAME( JOBVSL, 'N' ) ) THEN
+         IJOBVL = 1
+         ILVSL = .FALSE.
+      ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
+         IJOBVL = 2
+         ILVSL = .TRUE.
+      ELSE
+         IJOBVL = -1
+         ILVSL = .FALSE.
+      END IF
+*
+      IF( LSAME( JOBVSR, 'N' ) ) THEN
+         IJOBVR = 1
+         ILVSR = .FALSE.
+      ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
+         IJOBVR = 2
+         ILVSR = .TRUE.
+      ELSE
+         IJOBVR = -1
+         ILVSR = .FALSE.
+      END IF
+*
+      WANTST = LSAME( SORT, 'S' )
+*
+*     Test the input arguments
+*
+      INFO = 0
+      LQUERY = ( LWORK.EQ.-1 )
+      IF( IJOBVL.LE.0 ) THEN
+         INFO = -1
+      ELSE IF( IJOBVR.LE.0 ) THEN
+         INFO = -2
+      ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
+         INFO = -3
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -5
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -7
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -9
+      ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
+         INFO = -15
+      ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
+         INFO = -17
+      END IF
+*
+*     Compute workspace
+*      (Note: Comments in the code beginning "Workspace:" describe the
+*       minimal amount of workspace needed at that point in the code,
+*       as well as the preferred amount for good performance.
+*       NB refers to the optimal block size for the immediately
+*       following subroutine, as returned by ILAENV.)
+*
+      IF( INFO.EQ.0 ) THEN
+         IF( N.GT.0 )THEN
+            MINWRK = MAX( 8*N, 6*N + 16 )
+            MAXWRK = MINWRK - N +
+     $               N*ILAENV( 1, 'SGEQRF', ' ', N, 1, N, 0 )
+            MAXWRK = MAX( MAXWRK, MINWRK - N +
+     $                    N*ILAENV( 1, 'SORMQR', ' ', N, 1, N, -1 ) )
+            IF( ILVSL ) THEN
+               MAXWRK = MAX( MAXWRK, MINWRK - N +
+     $                       N*ILAENV( 1, 'SORGQR', ' ', N, 1, N, -1 ) )
+            END IF
+         ELSE
+            MINWRK = 1
+            MAXWRK = 1
+         END IF
+         WORK( 1 ) = MAXWRK
+*
+         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
+     $      INFO = -19
+      END IF
+*
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SGGES ', -INFO )
+         RETURN
+      ELSE IF( LQUERY ) THEN
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 ) THEN
+         SDIM = 0
+         RETURN
+      END IF
+*
+*     Get machine constants
+*
+      EPS = SLAMCH( 'P' )
+      SAFMIN = SLAMCH( 'S' )
+      SAFMAX = ONE / SAFMIN
+      CALL SLABAD( SAFMIN, SAFMAX )
+      SMLNUM = SQRT( SAFMIN ) / EPS
+      BIGNUM = ONE / SMLNUM
+*
+*     Scale A if max element outside range [SMLNUM,BIGNUM]
+*
+      ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
+      ILASCL = .FALSE.
+      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
+         ANRMTO = SMLNUM
+         ILASCL = .TRUE.
+      ELSE IF( ANRM.GT.BIGNUM ) THEN
+         ANRMTO = BIGNUM
+         ILASCL = .TRUE.
+      END IF
+      IF( ILASCL )
+     $   CALL SLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
+*
+*     Scale B if max element outside range [SMLNUM,BIGNUM]
+*
+      BNRM = SLANGE( 'M', N, N, B, LDB, WORK )
+      ILBSCL = .FALSE.
+      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
+         BNRMTO = SMLNUM
+         ILBSCL = .TRUE.
+      ELSE IF( BNRM.GT.BIGNUM ) THEN
+         BNRMTO = BIGNUM
+         ILBSCL = .TRUE.
+      END IF
+      IF( ILBSCL )
+     $   CALL SLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
+*
+*     Permute the matrix to make it more nearly triangular
+*     (Workspace: need 6*N + 2*N space for storing balancing factors)
+*
+      ILEFT = 1
+      IRIGHT = N + 1
+      IWRK = IRIGHT + N
+      CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
+     $             WORK( IRIGHT ), WORK( IWRK ), IERR )
+*
+*     Reduce B to triangular form (QR decomposition of B)
+*     (Workspace: need N, prefer N*NB)
+*
+      IROWS = IHI + 1 - ILO
+      ICOLS = N + 1 - ILO
+      ITAU = IWRK
+      IWRK = ITAU + IROWS
+      CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
+     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
+*
+*     Apply the orthogonal transformation to matrix A
+*     (Workspace: need N, prefer N*NB)
+*
+      CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
+     $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
+     $             LWORK+1-IWRK, IERR )
+*
+*     Initialize VSL
+*     (Workspace: need N, prefer N*NB)
+*
+      IF( ILVSL ) THEN
+         CALL SLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
+         IF( IROWS.GT.1 ) THEN
+            CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
+     $                   VSL( ILO+1, ILO ), LDVSL )
+         END IF
+         CALL SORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
+     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
+      END IF
+*
+*     Initialize VSR
+*
+      IF( ILVSR )
+     $   CALL SLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
+*
+*     Reduce to generalized Hessenberg form
+*     (Workspace: none needed)
+*
+      CALL SGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
+     $             LDVSL, VSR, LDVSR, IERR )
+*
+*     Perform QZ algorithm, computing Schur vectors if desired
+*     (Workspace: need N)
+*
+      IWRK = ITAU
+      CALL SHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
+     $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
+     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
+      IF( IERR.NE.0 ) THEN
+         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
+            INFO = IERR
+         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
+            INFO = IERR - N
+         ELSE
+            INFO = N + 1
+         END IF
+         GO TO 40
+      END IF
+*
+*     Sort eigenvalues ALPHA/BETA if desired
+*     (Workspace: need 4*N+16 )
+*
+      SDIM = 0
+      IF( WANTST ) THEN
+*
+*        Undo scaling on eigenvalues before SELCTGing
+*
+         IF( ILASCL ) THEN
+            CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
+     $                   IERR )
+            CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
+     $                   IERR )
+         END IF
+         IF( ILBSCL )
+     $      CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
+*
+*        Select eigenvalues
+*
+         DO 10 I = 1, N
+            BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
+   10    CONTINUE
+*
+         CALL STGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
+     $                ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
+     $                PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
+     $                IERR )
+         IF( IERR.EQ.1 )
+     $      INFO = N + 3
+*
+      END IF
+*
+*     Apply back-permutation to VSL and VSR
+*     (Workspace: none needed)
+*
+      IF( ILVSL )
+     $   CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
+     $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
+*
+      IF( ILVSR )
+     $   CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
+     $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
+*
+*     Check if unscaling would cause over/underflow, if so, rescale 
+*     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of 
+*     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
+*
+      IF( ILASCL )THEN
+         DO 50 I = 1, N 
+            IF( ALPHAI( I ).NE.ZERO ) THEN 
+               IF( ( ALPHAR( I )/SAFMAX ).GT.( ANRMTO/ANRM ) .OR.
+     $             ( SAFMIN/ALPHAR( I ) ).GT.( ANRM/ANRMTO ) ) THEN
+                  WORK( 1 ) = ABS( A( I, I )/ALPHAR( I ) )
+                  BETA( I ) = BETA( I )*WORK( 1 )
+                  ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
+                  ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
+               ELSE IF( ( ALPHAI( I )/SAFMAX ).GT.( ANRMTO/ANRM ) .OR.
+     $             ( SAFMIN/ALPHAI( I ) ).GT.( ANRM/ANRMTO ) ) THEN
+                  WORK( 1 ) = ABS( A( I, I+1 )/ALPHAI( I ) )
+                  BETA( I ) = BETA( I )*WORK( 1 )
+                  ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
+                  ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
+               END IF
+            END IF
+   50    CONTINUE
+      END IF 
+*
+      IF( ILBSCL )THEN 
+         DO 60 I = 1, N
+            IF( ALPHAI( I ).NE.ZERO ) THEN
+                IF( ( BETA( I )/SAFMAX ).GT.( BNRMTO/BNRM ) .OR.
+     $              ( SAFMIN/BETA( I ) ).GT.( BNRM/BNRMTO ) ) THEN
+                   WORK( 1 ) = ABS(B( I, I )/BETA( I ))
+                   BETA( I ) = BETA( I )*WORK( 1 )
+                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
+                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
+                END IF 
+             END IF
+   60    CONTINUE 
+      END IF 
+*
+*     Undo scaling
+*
+      IF( ILASCL ) THEN
+         CALL SLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
+         CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
+         CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
+      END IF
+*
+      IF( ILBSCL ) THEN
+         CALL SLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
+         CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
+      END IF
+*
+      IF( WANTST ) THEN
+*
+*        Check if reordering is correct
+*
+         LASTSL = .TRUE.
+         LST2SL = .TRUE.
+         SDIM = 0
+         IP = 0
+         DO 30 I = 1, N
+            CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
+            IF( ALPHAI( I ).EQ.ZERO ) THEN
+               IF( CURSL )
+     $            SDIM = SDIM + 1
+               IP = 0
+               IF( CURSL .AND. .NOT.LASTSL )
+     $            INFO = N + 2
+            ELSE
+               IF( IP.EQ.1 ) THEN
+*
+*                 Last eigenvalue of conjugate pair
+*
+                  CURSL = CURSL .OR. LASTSL
+                  LASTSL = CURSL
+                  IF( CURSL )
+     $               SDIM = SDIM + 2
+                  IP = -1
+                  IF( CURSL .AND. .NOT.LST2SL )
+     $               INFO = N + 2
+               ELSE
+*
+*                 First eigenvalue of conjugate pair
+*
+                  IP = 1
+               END IF
+            END IF
+            LST2SL = LASTSL
+            LASTSL = CURSL
+   30    CONTINUE
+*
+      END IF
+*
+   40 CONTINUE
+*
+      WORK( 1 ) = MAXWRK
+*
+      RETURN
+*
+*     End of SGGES
+*
+      END