LAPACK: SRC/dgges.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
00002      $                  SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
00003      $                  LDVSR, WORK, LWORK, BWORK, INFO )
00004 *
00005 *  -- LAPACK driver routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          JOBVSL, JOBVSR, SORT
00012       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            BWORK( * )
00016       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00017      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
00018      $                   VSR( LDVSR, * ), WORK( * )
00019 *     ..
00020 *     .. Function Arguments ..
00021       LOGICAL            SELCTG
00022       EXTERNAL           SELCTG
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *  DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
00029 *  the generalized eigenvalues, the generalized real Schur form (S,T),
00030 *  optionally, the left and/or right matrices of Schur vectors (VSL and
00031 *  VSR). This gives the generalized Schur factorization
00032 *
00033 *           (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
00034 *
00035 *  Optionally, it also orders the eigenvalues so that a selected cluster
00036 *  of eigenvalues appears in the leading diagonal blocks of the upper
00037 *  quasi-triangular matrix S and the upper triangular matrix T.The
00038 *  leading columns of VSL and VSR then form an orthonormal basis for the
00039 *  corresponding left and right eigenspaces (deflating subspaces).
00040 *
00041 *  (If only the generalized eigenvalues are needed, use the driver
00042 *  DGGEV instead, which is faster.)
00043 *
00044 *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
00045 *  or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
00046 *  usually represented as the pair (alpha,beta), as there is a
00047 *  reasonable interpretation for beta=0 or both being zero.
00048 *
00049 *  A pair of matrices (S,T) is in generalized real Schur form if T is
00050 *  upper triangular with non-negative diagonal and S is block upper
00051 *  triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
00052 *  to real generalized eigenvalues, while 2-by-2 blocks of S will be
00053 *  "standardized" by making the corresponding elements of T have the
00054 *  form:
00055 *          [  a  0  ]
00056 *          [  0  b  ]
00057 *
00058 *  and the pair of corresponding 2-by-2 blocks in S and T will have a
00059 *  complex conjugate pair of generalized eigenvalues.
00060 *
00061 *
00062 *  Arguments
00063 *  =========
00064 *
00065 *  JOBVSL  (input) CHARACTER*1
00066 *          = 'N':  do not compute the left Schur vectors;
00067 *          = 'V':  compute the left Schur vectors.
00068 *
00069 *  JOBVSR  (input) CHARACTER*1
00070 *          = 'N':  do not compute the right Schur vectors;
00071 *          = 'V':  compute the right Schur vectors.
00072 *
00073 *  SORT    (input) CHARACTER*1
00074 *          Specifies whether or not to order the eigenvalues on the
00075 *          diagonal of the generalized Schur form.
00076 *          = 'N':  Eigenvalues are not ordered;
00077 *          = 'S':  Eigenvalues are ordered (see SELCTG);
00078 *
00079 *  SELCTG  (external procedure) LOGICAL FUNCTION of three DOUBLE PRECISION arguments
00080 *          SELCTG must be declared EXTERNAL in the calling subroutine.
00081 *          If SORT = 'N', SELCTG is not referenced.
00082 *          If SORT = 'S', SELCTG is used to select eigenvalues to sort
00083 *          to the top left of the Schur form.
00084 *          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
00085 *          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
00086 *          one of a complex conjugate pair of eigenvalues is selected,
00087 *          then both complex eigenvalues are selected.
00088 *
00089 *          Note that in the ill-conditioned case, a selected complex
00090 *          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
00091 *          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
00092 *          in this case.
00093 *
00094 *  N       (input) INTEGER
00095 *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
00096 *
00097 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
00098 *          On entry, the first of the pair of matrices.
00099 *          On exit, A has been overwritten by its generalized Schur
00100 *          form S.
00101 *
00102 *  LDA     (input) INTEGER
00103 *          The leading dimension of A.  LDA >= max(1,N).
00104 *
00105 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
00106 *          On entry, the second of the pair of matrices.
00107 *          On exit, B has been overwritten by its generalized Schur
00108 *          form T.
00109 *
00110 *  LDB     (input) INTEGER
00111 *          The leading dimension of B.  LDB >= max(1,N).
00112 *
00113 *  SDIM    (output) INTEGER
00114 *          If SORT = 'N', SDIM = 0.
00115 *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
00116 *          for which SELCTG is true.  (Complex conjugate pairs for which
00117 *          SELCTG is true for either eigenvalue count as 2.)
00118 *
00119 *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
00120 *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
00121 *  BETA    (output) DOUBLE PRECISION array, dimension (N)
00122 *          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
00123 *          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
00124 *          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
00125 *          form (S,T) that would result if the 2-by-2 diagonal blocks of
00126 *          the real Schur form of (A,B) were further reduced to
00127 *          triangular form using 2-by-2 complex unitary transformations.
00128 *          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
00129 *          positive, then the j-th and (j+1)-st eigenvalues are a
00130 *          complex conjugate pair, with ALPHAI(j+1) negative.
00131 *
00132 *          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
00133 *          may easily over- or underflow, and BETA(j) may even be zero.
00134 *          Thus, the user should avoid naively computing the ratio.
00135 *          However, ALPHAR and ALPHAI will be always less than and
00136 *          usually comparable with norm(A) in magnitude, and BETA always
00137 *          less than and usually comparable with norm(B).
00138 *
00139 *  VSL     (output) DOUBLE PRECISION array, dimension (LDVSL,N)
00140 *          If JOBVSL = 'V', VSL will contain the left Schur vectors.
00141 *          Not referenced if JOBVSL = 'N'.
00142 *
00143 *  LDVSL   (input) INTEGER
00144 *          The leading dimension of the matrix VSL. LDVSL >=1, and
00145 *          if JOBVSL = 'V', LDVSL >= N.
00146 *
00147 *  VSR     (output) DOUBLE PRECISION array, dimension (LDVSR,N)
00148 *          If JOBVSR = 'V', VSR will contain the right Schur vectors.
00149 *          Not referenced if JOBVSR = 'N'.
00150 *
00151 *  LDVSR   (input) INTEGER
00152 *          The leading dimension of the matrix VSR. LDVSR >= 1, and
00153 *          if JOBVSR = 'V', LDVSR >= N.
00154 *
00155 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00156 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00157 *
00158 *  LWORK   (input) INTEGER
00159 *          The dimension of the array WORK.
00160 *          If N = 0, LWORK >= 1, else LWORK >= 8*N+16.
00161 *          For good performance , LWORK must generally be larger.
00162 *
00163 *          If LWORK = -1, then a workspace query is assumed; the routine
00164 *          only calculates the optimal size of the WORK array, returns
00165 *          this value as the first entry of the WORK array, and no error
00166 *          message related to LWORK is issued by XERBLA.
00167 *
00168 *  BWORK   (workspace) LOGICAL array, dimension (N)
00169 *          Not referenced if SORT = 'N'.
00170 *
00171 *  INFO    (output) INTEGER
00172 *          = 0:  successful exit
00173 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00174 *          = 1,...,N:
00175 *                The QZ iteration failed.  (A,B) are not in Schur
00176 *                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
00177 *                be correct for j=INFO+1,...,N.
00178 *          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
00179 *                =N+2: after reordering, roundoff changed values of
00180 *                      some complex eigenvalues so that leading
00181 *                      eigenvalues in the Generalized Schur form no
00182 *                      longer satisfy SELCTG=.TRUE.  This could also
00183 *                      be caused due to scaling.
00184 *                =N+3: reordering failed in DTGSEN.
00185 *
00186 *  =====================================================================
00187 *
00188 *     .. Parameters ..
00189       DOUBLE PRECISION   ZERO, ONE
00190       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00191 *     ..
00192 *     .. Local Scalars ..
00193       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
00194      $                   LQUERY, LST2SL, WANTST
00195       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
00196      $                   ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
00197      $                   MINWRK
00198       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
00199      $                   PVSR, SAFMAX, SAFMIN, SMLNUM
00200 *     ..
00201 *     .. Local Arrays ..
00202       INTEGER            IDUM( 1 )
00203       DOUBLE PRECISION   DIF( 2 )
00204 *     ..
00205 *     .. External Subroutines ..
00206       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
00207      $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
00208      $                   XERBLA
00209 *     ..
00210 *     .. External Functions ..
00211       LOGICAL            LSAME
00212       INTEGER            ILAENV
00213       DOUBLE PRECISION   DLAMCH, DLANGE
00214       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
00215 *     ..
00216 *     .. Intrinsic Functions ..
00217       INTRINSIC          ABS, MAX, SQRT
00218 *     ..
00219 *     .. Executable Statements ..
00220 *
00221 *     Decode the input arguments
00222 *
00223       IF( LSAME( JOBVSL, 'N' ) ) THEN
00224          IJOBVL = 1
00225          ILVSL = .FALSE.
00226       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00227          IJOBVL = 2
00228          ILVSL = .TRUE.
00229       ELSE
00230          IJOBVL = -1
00231          ILVSL = .FALSE.
00232       END IF
00233 *
00234       IF( LSAME( JOBVSR, 'N' ) ) THEN
00235          IJOBVR = 1
00236          ILVSR = .FALSE.
00237       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00238          IJOBVR = 2
00239          ILVSR = .TRUE.
00240       ELSE
00241          IJOBVR = -1
00242          ILVSR = .FALSE.
00243       END IF
00244 *
00245       WANTST = LSAME( SORT, 'S' )
00246 *
00247 *     Test the input arguments
00248 *
00249       INFO = 0
00250       LQUERY = ( LWORK.EQ.-1 )
00251       IF( IJOBVL.LE.0 ) THEN
00252          INFO = -1
00253       ELSE IF( IJOBVR.LE.0 ) THEN
00254          INFO = -2
00255       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
00256          INFO = -3
00257       ELSE IF( N.LT.0 ) THEN
00258          INFO = -5
00259       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00260          INFO = -7
00261       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00262          INFO = -9
00263       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00264          INFO = -15
00265       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00266          INFO = -17
00267       END IF
00268 *
00269 *     Compute workspace
00270 *      (Note: Comments in the code beginning "Workspace:" describe the
00271 *       minimal amount of workspace needed at that point in the code,
00272 *       as well as the preferred amount for good performance.
00273 *       NB refers to the optimal block size for the immediately
00274 *       following subroutine, as returned by ILAENV.)
00275 *
00276       IF( INFO.EQ.0 ) THEN
00277          IF( N.GT.0 )THEN
00278             MINWRK = MAX( 8*N, 6*N + 16 )
00279             MAXWRK = MINWRK - N +
00280      $               N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
00281             MAXWRK = MAX( MAXWRK, MINWRK - N +
00282      $                    N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
00283             IF( ILVSL ) THEN
00284                MAXWRK = MAX( MAXWRK, MINWRK - N +
00285      $                       N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
00286             END IF
00287          ELSE
00288             MINWRK = 1
00289             MAXWRK = 1
00290          END IF
00291          WORK( 1 ) = MAXWRK
00292 *
00293          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
00294      $      INFO = -19
00295       END IF
00296 *
00297       IF( INFO.NE.0 ) THEN
00298          CALL XERBLA( 'DGGES ', -INFO )
00299          RETURN
00300       ELSE IF( LQUERY ) THEN
00301          RETURN
00302       END IF
00303 *
00304 *     Quick return if possible
00305 *
00306       IF( N.EQ.0 ) THEN
00307          SDIM = 0
00308          RETURN
00309       END IF
00310 *
00311 *     Get machine constants
00312 *
00313       EPS = DLAMCH( 'P' )
00314       SAFMIN = DLAMCH( 'S' )
00315       SAFMAX = ONE / SAFMIN
00316       CALL DLABAD( SAFMIN, SAFMAX )
00317       SMLNUM = SQRT( SAFMIN ) / EPS
00318       BIGNUM = ONE / SMLNUM
00319 *
00320 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00321 *
00322       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
00323       ILASCL = .FALSE.
00324       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00325          ANRMTO = SMLNUM
00326          ILASCL = .TRUE.
00327       ELSE IF( ANRM.GT.BIGNUM ) THEN
00328          ANRMTO = BIGNUM
00329          ILASCL = .TRUE.
00330       END IF
00331       IF( ILASCL )
00332      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
00333 *
00334 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00335 *
00336       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
00337       ILBSCL = .FALSE.
00338       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00339          BNRMTO = SMLNUM
00340          ILBSCL = .TRUE.
00341       ELSE IF( BNRM.GT.BIGNUM ) THEN
00342          BNRMTO = BIGNUM
00343          ILBSCL = .TRUE.
00344       END IF
00345       IF( ILBSCL )
00346      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
00347 *
00348 *     Permute the matrix to make it more nearly triangular
00349 *     (Workspace: need 6*N + 2*N space for storing balancing factors)
00350 *
00351       ILEFT = 1
00352       IRIGHT = N + 1
00353       IWRK = IRIGHT + N
00354       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
00355      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
00356 *
00357 *     Reduce B to triangular form (QR decomposition of B)
00358 *     (Workspace: need N, prefer N*NB)
00359 *
00360       IROWS = IHI + 1 - ILO
00361       ICOLS = N + 1 - ILO
00362       ITAU = IWRK
00363       IWRK = ITAU + IROWS
00364       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00365      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00366 *
00367 *     Apply the orthogonal transformation to matrix A
00368 *     (Workspace: need N, prefer N*NB)
00369 *
00370       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00371      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
00372      $             LWORK+1-IWRK, IERR )
00373 *
00374 *     Initialize VSL
00375 *     (Workspace: need N, prefer N*NB)
00376 *
00377       IF( ILVSL ) THEN
00378          CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
00379          IF( IROWS.GT.1 ) THEN
00380             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00381      $                   VSL( ILO+1, ILO ), LDVSL )
00382          END IF
00383          CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00384      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
00385       END IF
00386 *
00387 *     Initialize VSR
00388 *
00389       IF( ILVSR )
00390      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
00391 *
00392 *     Reduce to generalized Hessenberg form
00393 *     (Workspace: none needed)
00394 *
00395       CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00396      $             LDVSL, VSR, LDVSR, IERR )
00397 *
00398 *     Perform QZ algorithm, computing Schur vectors if desired
00399 *     (Workspace: need N)
00400 *
00401       IWRK = ITAU
00402       CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00403      $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
00404      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00405       IF( IERR.NE.0 ) THEN
00406          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
00407             INFO = IERR
00408          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
00409             INFO = IERR - N
00410          ELSE
00411             INFO = N + 1
00412          END IF
00413          GO TO 50
00414       END IF
00415 *
00416 *     Sort eigenvalues ALPHA/BETA if desired
00417 *     (Workspace: need 4*N+16 )
00418 *
00419       SDIM = 0
00420       IF( WANTST ) THEN
00421 *
00422 *        Undo scaling on eigenvalues before SELCTGing
00423 *
00424          IF( ILASCL ) THEN
00425             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
00426      $                   IERR )
00427             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
00428      $                   IERR )
00429          END IF
00430          IF( ILBSCL )
00431      $      CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00432 *
00433 *        Select eigenvalues
00434 *
00435          DO 10 I = 1, N
00436             BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
00437    10    CONTINUE
00438 *
00439          CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
00440      $                ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
00441      $                PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
00442      $                IERR )
00443          IF( IERR.EQ.1 )
00444      $      INFO = N + 3
00445 *
00446       END IF
00447 *
00448 *     Apply back-permutation to VSL and VSR
00449 *     (Workspace: none needed)
00450 *
00451       IF( ILVSL )
00452      $   CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
00453      $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
00454 *
00455       IF( ILVSR )
00456      $   CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
00457      $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
00458 *
00459 *     Check if unscaling would cause over/underflow, if so, rescale
00460 *     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
00461 *     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
00462 *
00463       IF( ILASCL ) THEN
00464          DO 20 I = 1, N
00465             IF( ALPHAI( I ).NE.ZERO ) THEN
00466                IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
00467      $             ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
00468                   WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
00469                   BETA( I ) = BETA( I )*WORK( 1 )
00470                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
00471                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
00472                ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
00473      $                  ( ANRMTO / ANRM ) .OR.
00474      $                  ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
00475      $                   THEN
00476                   WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
00477                   BETA( I ) = BETA( I )*WORK( 1 )
00478                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
00479                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
00480                END IF
00481             END IF
00482    20    CONTINUE
00483       END IF
00484 *
00485       IF( ILBSCL ) THEN
00486          DO 30 I = 1, N
00487             IF( ALPHAI( I ).NE.ZERO ) THEN
00488                IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
00489      $             ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
00490                   WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
00491                   BETA( I ) = BETA( I )*WORK( 1 )
00492                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
00493                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
00494                END IF
00495             END IF
00496    30    CONTINUE
00497       END IF
00498 *
00499 *     Undo scaling
00500 *
00501       IF( ILASCL ) THEN
00502          CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
00503          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
00504          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
00505       END IF
00506 *
00507       IF( ILBSCL ) THEN
00508          CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
00509          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00510       END IF
00511 *
00512       IF( WANTST ) THEN
00513 *
00514 *        Check if reordering is correct
00515 *
00516          LASTSL = .TRUE.
00517          LST2SL = .TRUE.
00518          SDIM = 0
00519          IP = 0
00520          DO 40 I = 1, N
00521             CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
00522             IF( ALPHAI( I ).EQ.ZERO ) THEN
00523                IF( CURSL )
00524      $            SDIM = SDIM + 1
00525                IP = 0
00526                IF( CURSL .AND. .NOT.LASTSL )
00527      $            INFO = N + 2
00528             ELSE
00529                IF( IP.EQ.1 ) THEN
00530 *
00531 *                 Last eigenvalue of conjugate pair
00532 *
00533                   CURSL = CURSL .OR. LASTSL
00534                   LASTSL = CURSL
00535                   IF( CURSL )
00536      $               SDIM = SDIM + 2
00537                   IP = -1
00538                   IF( CURSL .AND. .NOT.LST2SL )
00539      $               INFO = N + 2
00540                ELSE
00541 *
00542 *                 First eigenvalue of conjugate pair
00543 *
00544                   IP = 1
00545                END IF
00546             END IF
00547             LST2SL = LASTSL
00548             LASTSL = CURSL
00549    40    CONTINUE
00550 *
00551       END IF
00552 *
00553    50 CONTINUE
00554 *
00555       WORK( 1 ) = MAXWRK
00556 *
00557       RETURN
00558 *
00559 *     End of DGGES
00560 *
00561       END