SRC/dgges3.f

   1 *> \brief <b> DGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGGES3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgges3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgges3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgges3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
  22 *                          SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
  23 *                          LDVSR, WORK, LWORK, BWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBVSL, JOBVSR, SORT
  27 *       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       LOGICAL            BWORK( * )
  31 *       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
  32 *      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
  33 *      $                   VSR( LDVSR, * ), WORK( * )
  34 *       ..
  35 *       .. Function Arguments ..
  36 *       LOGICAL            SELCTG
  37 *       EXTERNAL           SELCTG
  38 *       ..
  39 *
  40 *
  41 *> \par Purpose:
  42 *  =============
  43 *>
  44 *> \verbatim
  45 *>
  46 *> DGGES3 computes for a pair of N-by-N real nonsymmetric matrices (A,B),
  47 *> the generalized eigenvalues, the generalized real Schur form (S,T),
  48 *> optionally, the left and/or right matrices of Schur vectors (VSL and
  49 *> VSR). This gives the generalized Schur factorization
  50 *>
  51 *>          (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
  52 *>
  53 *> Optionally, it also orders the eigenvalues so that a selected cluster
  54 *> of eigenvalues appears in the leading diagonal blocks of the upper
  55 *> quasi-triangular matrix S and the upper triangular matrix T.The
  56 *> leading columns of VSL and VSR then form an orthonormal basis for the
  57 *> corresponding left and right eigenspaces (deflating subspaces).
  58 *>
  59 *> (If only the generalized eigenvalues are needed, use the driver
  60 *> DGGEV instead, which is faster.)
  61 *>
  62 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
  63 *> or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
  64 *> usually represented as the pair (alpha,beta), as there is a
  65 *> reasonable interpretation for beta=0 or both being zero.
  66 *>
  67 *> A pair of matrices (S,T) is in generalized real Schur form if T is
  68 *> upper triangular with non-negative diagonal and S is block upper
  69 *> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
  70 *> to real generalized eigenvalues, while 2-by-2 blocks of S will be
  71 *> "standardized" by making the corresponding elements of T have the
  72 *> form:
  73 *>         [  a  0  ]
  74 *>         [  0  b  ]
  75 *>
  76 *> and the pair of corresponding 2-by-2 blocks in S and T will have a
  77 *> complex conjugate pair of generalized eigenvalues.
  78 *>
  79 *> \endverbatim
  80 *
  81 *  Arguments:
  82 *  ==========
  83 *
  84 *> \param[in] JOBVSL
  85 *> \verbatim
  86 *>          JOBVSL is CHARACTER*1
  87 *>          = 'N':  do not compute the left Schur vectors;
  88 *>          = 'V':  compute the left Schur vectors.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] JOBVSR
  92 *> \verbatim
  93 *>          JOBVSR is CHARACTER*1
  94 *>          = 'N':  do not compute the right Schur vectors;
  95 *>          = 'V':  compute the right Schur vectors.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] SORT
  99 *> \verbatim
 100 *>          SORT is CHARACTER*1
 101 *>          Specifies whether or not to order the eigenvalues on the
 102 *>          diagonal of the generalized Schur form.
 103 *>          = 'N':  Eigenvalues are not ordered;
 104 *>          = 'S':  Eigenvalues are ordered (see SELCTG);
 105 *> \endverbatim
 106 *>
 107 *> \param[in] SELCTG
 108 *> \verbatim
 109 *>          SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments
 110 *>          SELCTG must be declared EXTERNAL in the calling subroutine.
 111 *>          If SORT = 'N', SELCTG is not referenced.
 112 *>          If SORT = 'S', SELCTG is used to select eigenvalues to sort
 113 *>          to the top left of the Schur form.
 114 *>          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
 115 *>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
 116 *>          one of a complex conjugate pair of eigenvalues is selected,
 117 *>          then both complex eigenvalues are selected.
 118 *>
 119 *>          Note that in the ill-conditioned case, a selected complex
 120 *>          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
 121 *>          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
 122 *>          in this case.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] N
 126 *> \verbatim
 127 *>          N is INTEGER
 128 *>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
 129 *> \endverbatim
 130 *>
 131 *> \param[in,out] A
 132 *> \verbatim
 133 *>          A is DOUBLE PRECISION array, dimension (LDA, N)
 134 *>          On entry, the first of the pair of matrices.
 135 *>          On exit, A has been overwritten by its generalized Schur
 136 *>          form S.
 137 *> \endverbatim
 138 *>
 139 *> \param[in] LDA
 140 *> \verbatim
 141 *>          LDA is INTEGER
 142 *>          The leading dimension of A.  LDA >= max(1,N).
 143 *> \endverbatim
 144 *>
 145 *> \param[in,out] B
 146 *> \verbatim
 147 *>          B is DOUBLE PRECISION array, dimension (LDB, N)
 148 *>          On entry, the second of the pair of matrices.
 149 *>          On exit, B has been overwritten by its generalized Schur
 150 *>          form T.
 151 *> \endverbatim
 152 *>
 153 *> \param[in] LDB
 154 *> \verbatim
 155 *>          LDB is INTEGER
 156 *>          The leading dimension of B.  LDB >= max(1,N).
 157 *> \endverbatim
 158 *>
 159 *> \param[out] SDIM
 160 *> \verbatim
 161 *>          SDIM is INTEGER
 162 *>          If SORT = 'N', SDIM = 0.
 163 *>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
 164 *>          for which SELCTG is true.  (Complex conjugate pairs for which
 165 *>          SELCTG is true for either eigenvalue count as 2.)
 166 *> \endverbatim
 167 *>
 168 *> \param[out] ALPHAR
 169 *> \verbatim
 170 *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
 171 *> \endverbatim
 172 *>
 173 *> \param[out] ALPHAI
 174 *> \verbatim
 175 *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
 176 *> \endverbatim
 177 *>
 178 *> \param[out] BETA
 179 *> \verbatim
 180 *>          BETA is DOUBLE PRECISION array, dimension (N)
 181 *>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
 182 *>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
 183 *>          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
 184 *>          form (S,T) that would result if the 2-by-2 diagonal blocks of
 185 *>          the real Schur form of (A,B) were further reduced to
 186 *>          triangular form using 2-by-2 complex unitary transformations.
 187 *>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
 188 *>          positive, then the j-th and (j+1)-st eigenvalues are a
 189 *>          complex conjugate pair, with ALPHAI(j+1) negative.
 190 *>
 191 *>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
 192 *>          may easily over- or underflow, and BETA(j) may even be zero.
 193 *>          Thus, the user should avoid naively computing the ratio.
 194 *>          However, ALPHAR and ALPHAI will be always less than and
 195 *>          usually comparable with norm(A) in magnitude, and BETA always
 196 *>          less than and usually comparable with norm(B).
 197 *> \endverbatim
 198 *>
 199 *> \param[out] VSL
 200 *> \verbatim
 201 *>          VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
 202 *>          If JOBVSL = 'V', VSL will contain the left Schur vectors.
 203 *>          Not referenced if JOBVSL = 'N'.
 204 *> \endverbatim
 205 *>
 206 *> \param[in] LDVSL
 207 *> \verbatim
 208 *>          LDVSL is INTEGER
 209 *>          The leading dimension of the matrix VSL. LDVSL >=1, and
 210 *>          if JOBVSL = 'V', LDVSL >= N.
 211 *> \endverbatim
 212 *>
 213 *> \param[out] VSR
 214 *> \verbatim
 215 *>          VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
 216 *>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
 217 *>          Not referenced if JOBVSR = 'N'.
 218 *> \endverbatim
 219 *>
 220 *> \param[in] LDVSR
 221 *> \verbatim
 222 *>          LDVSR is INTEGER
 223 *>          The leading dimension of the matrix VSR. LDVSR >= 1, and
 224 *>          if JOBVSR = 'V', LDVSR >= N.
 225 *> \endverbatim
 226 *>
 227 *> \param[out] WORK
 228 *> \verbatim
 229 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 230 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 231 *> \endverbatim
 232 *>
 233 *> \param[in] LWORK
 234 *> \verbatim
 235 *>          LWORK is INTEGER
 236 *>          The dimension of the array WORK.
 237 *>
 238 *>          If LWORK = -1, then a workspace query is assumed; the routine
 239 *>          only calculates the optimal size of the WORK array, returns
 240 *>          this value as the first entry of the WORK array, and no error
 241 *>          message related to LWORK is issued by XERBLA.
 242 *> \endverbatim
 243 *>
 244 *> \param[out] BWORK
 245 *> \verbatim
 246 *>          BWORK is LOGICAL array, dimension (N)
 247 *>          Not referenced if SORT = 'N'.
 248 *> \endverbatim
 249 *>
 250 *> \param[out] INFO
 251 *> \verbatim
 252 *>          INFO is INTEGER
 253 *>          = 0:  successful exit
 254 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 255 *>          = 1,...,N:
 256 *>                The QZ iteration failed.  (A,B) are not in Schur
 257 *>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
 258 *>                be correct for j=INFO+1,...,N.
 259 *>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
 260 *>                =N+2: after reordering, roundoff changed values of
 261 *>                      some complex eigenvalues so that leading
 262 *>                      eigenvalues in the Generalized Schur form no
 263 *>                      longer satisfy SELCTG=.TRUE.  This could also
 264 *>                      be caused due to scaling.
 265 *>                =N+3: reordering failed in DTGSEN.
 266 *> \endverbatim
 267 *
 268 *  Authors:
 269 *  ========
 270 *
 271 *> \author Univ. of Tennessee
 272 *> \author Univ. of California Berkeley
 273 *> \author Univ. of Colorado Denver
 274 *> \author NAG Ltd.
 275 *
 276 *> \date January 2015
 277 *
 278 *> \ingroup doubleGEeigen
 279 *
 280 *  =====================================================================
 281       SUBROUTINE DGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
 282      $                   LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
 283      $                   VSR, LDVSR, WORK, LWORK, BWORK, INFO )
 284 *
 285 *  -- LAPACK driver routine (version 3.6.0) --
 286 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 287 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 288 *     January 2015
 289 *
 290 *     .. Scalar Arguments ..
 291       CHARACTER          JOBVSL, JOBVSR, SORT
 292       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
 293 *     ..
 294 *     .. Array Arguments ..
 295       LOGICAL            BWORK( * )
 296       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
 297      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
 298      $                   VSR( LDVSR, * ), WORK( * )
 299 *     ..
 300 *     .. Function Arguments ..
 301       LOGICAL            SELCTG
 302       EXTERNAL           SELCTG
 303 *     ..
 304 *
 305 *  =====================================================================
 306 *
 307 *     .. Parameters ..
 308       DOUBLE PRECISION   ZERO, ONE
 309       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 310 *     ..
 311 *     .. Local Scalars ..
 312       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
 313      $                   LQUERY, LST2SL, WANTST
 314       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
 315      $                   ILO, IP, IRIGHT, IROWS, ITAU, IWRK, LWKOPT
 316       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
 317      $                   PVSR, SAFMAX, SAFMIN, SMLNUM
 318 *     ..
 319 *     .. Local Arrays ..
 320       INTEGER            IDUM( 1 )
 321       DOUBLE PRECISION   DIF( 2 )
 322 *     ..
 323 *     .. External Subroutines ..
 324       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHD3, DHGEQZ, DLABAD,
 325      $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
 326      $                   XERBLA
 327 *     ..
 328 *     .. External Functions ..
 329       LOGICAL            LSAME
 330       DOUBLE PRECISION   DLAMCH, DLANGE
 331       EXTERNAL           LSAME, DLAMCH, DLANGE
 332 *     ..
 333 *     .. Intrinsic Functions ..
 334       INTRINSIC          ABS, MAX, SQRT
 335 *     ..
 336 *     .. Executable Statements ..
 337 *
 338 *     Decode the input arguments
 339 *
 340       IF( LSAME( JOBVSL, 'N' ) ) THEN
 341          IJOBVL = 1
 342          ILVSL = .FALSE.
 343       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
 344          IJOBVL = 2
 345          ILVSL = .TRUE.
 346       ELSE
 347          IJOBVL = -1
 348          ILVSL = .FALSE.
 349       END IF
 350 *
 351       IF( LSAME( JOBVSR, 'N' ) ) THEN
 352          IJOBVR = 1
 353          ILVSR = .FALSE.
 354       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
 355          IJOBVR = 2
 356          ILVSR = .TRUE.
 357       ELSE
 358          IJOBVR = -1
 359          ILVSR = .FALSE.
 360       END IF
 361 *
 362       WANTST = LSAME( SORT, 'S' )
 363 *
 364 *     Test the input arguments
 365 *
 366       INFO = 0
 367       LQUERY = ( LWORK.EQ.-1 )
 368       IF( IJOBVL.LE.0 ) THEN
 369          INFO = -1
 370       ELSE IF( IJOBVR.LE.0 ) THEN
 371          INFO = -2
 372       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
 373          INFO = -3
 374       ELSE IF( N.LT.0 ) THEN
 375          INFO = -5
 376       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 377          INFO = -7
 378       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 379          INFO = -9
 380       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
 381          INFO = -15
 382       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
 383          INFO = -17
 384       ELSE IF( LWORK.LT.6*N+16 .AND. .NOT.LQUERY ) THEN
 385          INFO = -19
 386       END IF
 387 *
 388 *     Compute workspace
 389 *
 390       IF( INFO.EQ.0 ) THEN
 391          CALL DGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
 392          LWKOPT = MAX( 6*N+16, 3*N+INT( WORK ( 1 ) ) )
 393          CALL DORMQR( 'L', 'T', N, N, N, B, LDB, WORK, A, LDA, WORK,
 394      $                -1, IERR )
 395          LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
 396          IF( ILVSL ) THEN
 397             CALL DORGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1, IERR )
 398             LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
 399          END IF
 400          CALL DGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL,
 401      $                LDVSL, VSR, LDVSR, WORK, -1, IERR )
 402          LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
 403          CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB,
 404      $                ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
 405      $                WORK, -1, IERR )
 406          LWKOPT = MAX( LWKOPT, 2*N+INT( WORK ( 1 ) ) )
 407          IF( WANTST ) THEN
 408             CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
 409      $                   ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
 410      $                   SDIM, PVSL, PVSR, DIF, WORK, -1, IDUM, 1,
 411      $                   IERR )
 412             LWKOPT = MAX( LWKOPT, 2*N+INT( WORK ( 1 ) ) )
 413          END IF
 414          WORK( 1 ) = LWKOPT
 415       END IF
 416 *
 417       IF( INFO.NE.0 ) THEN
 418          CALL XERBLA( 'DGGES3 ', -INFO )
 419          RETURN
 420       ELSE IF( LQUERY ) THEN
 421          RETURN
 422       END IF
 423 *
 424 *     Quick return if possible
 425 *
 426       IF( N.EQ.0 ) THEN
 427          SDIM = 0
 428          RETURN
 429       END IF
 430 *
 431 *     Get machine constants
 432 *
 433       EPS = DLAMCH( 'P' )
 434       SAFMIN = DLAMCH( 'S' )
 435       SAFMAX = ONE / SAFMIN
 436       CALL DLABAD( SAFMIN, SAFMAX )
 437       SMLNUM = SQRT( SAFMIN ) / EPS
 438       BIGNUM = ONE / SMLNUM
 439 *
 440 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 441 *
 442       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
 443       ILASCL = .FALSE.
 444       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 445          ANRMTO = SMLNUM
 446          ILASCL = .TRUE.
 447       ELSE IF( ANRM.GT.BIGNUM ) THEN
 448          ANRMTO = BIGNUM
 449          ILASCL = .TRUE.
 450       END IF
 451       IF( ILASCL )
 452      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 453 *
 454 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 455 *
 456       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
 457       ILBSCL = .FALSE.
 458       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 459          BNRMTO = SMLNUM
 460          ILBSCL = .TRUE.
 461       ELSE IF( BNRM.GT.BIGNUM ) THEN
 462          BNRMTO = BIGNUM
 463          ILBSCL = .TRUE.
 464       END IF
 465       IF( ILBSCL )
 466      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 467 *
 468 *     Permute the matrix to make it more nearly triangular
 469 *
 470       ILEFT = 1
 471       IRIGHT = N + 1
 472       IWRK = IRIGHT + N
 473       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
 474      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
 475 *
 476 *     Reduce B to triangular form (QR decomposition of B)
 477 *
 478       IROWS = IHI + 1 - ILO
 479       ICOLS = N + 1 - ILO
 480       ITAU = IWRK
 481       IWRK = ITAU + IROWS
 482       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 483      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 484 *
 485 *     Apply the orthogonal transformation to matrix A
 486 *
 487       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 488      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 489      $             LWORK+1-IWRK, IERR )
 490 *
 491 *     Initialize VSL
 492 *
 493       IF( ILVSL ) THEN
 494          CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
 495          IF( IROWS.GT.1 ) THEN
 496             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 497      $                   VSL( ILO+1, ILO ), LDVSL )
 498          END IF
 499          CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
 500      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 501       END IF
 502 *
 503 *     Initialize VSR
 504 *
 505       IF( ILVSR )
 506      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
 507 *
 508 *     Reduce to generalized Hessenberg form
 509 *
 510       CALL DGGHD3( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
 511      $             LDVSL, VSR, LDVSR, WORK( IWRK ), LWORK+1-IWRK,
 512      $             IERR )
 513 *
 514 *     Perform QZ algorithm, computing Schur vectors if desired
 515 *
 516       IWRK = ITAU
 517       CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
 518      $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
 519      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 520       IF( IERR.NE.0 ) THEN
 521          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
 522             INFO = IERR
 523          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
 524             INFO = IERR - N
 525          ELSE
 526             INFO = N + 1
 527          END IF
 528          GO TO 50
 529       END IF
 530 *
 531 *     Sort eigenvalues ALPHA/BETA if desired
 532 *
 533       SDIM = 0
 534       IF( WANTST ) THEN
 535 *
 536 *        Undo scaling on eigenvalues before SELCTGing
 537 *
 538          IF( ILASCL ) THEN
 539             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
 540      $                   IERR )
 541             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
 542      $                   IERR )
 543          END IF
 544          IF( ILBSCL )
 545      $      CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 546 *
 547 *        Select eigenvalues
 548 *
 549          DO 10 I = 1, N
 550             BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
 551    10    CONTINUE
 552 *
 553          CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
 554      $                ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
 555      $                PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
 556      $                IERR )
 557          IF( IERR.EQ.1 )
 558      $      INFO = N + 3
 559 *
 560       END IF
 561 *
 562 *     Apply back-permutation to VSL and VSR
 563 *
 564       IF( ILVSL )
 565      $   CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
 566      $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
 567 *
 568       IF( ILVSR )
 569      $   CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
 570      $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
 571 *
 572 *     Check if unscaling would cause over/underflow, if so, rescale
 573 *     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
 574 *     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
 575 *
 576       IF( ILASCL ) THEN
 577          DO 20 I = 1, N
 578             IF( ALPHAI( I ).NE.ZERO ) THEN
 579                IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
 580      $             ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
 581                   WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
 582                   BETA( I ) = BETA( I )*WORK( 1 )
 583                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
 584                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
 585                ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
 586      $                  ( ANRMTO / ANRM ) .OR.
 587      $                  ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
 588      $                   THEN
 589                   WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
 590                   BETA( I ) = BETA( I )*WORK( 1 )
 591                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
 592                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
 593                END IF
 594             END IF
 595    20    CONTINUE
 596       END IF
 597 *
 598       IF( ILBSCL ) THEN
 599          DO 30 I = 1, N
 600             IF( ALPHAI( I ).NE.ZERO ) THEN
 601                IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
 602      $             ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
 603                   WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
 604                   BETA( I ) = BETA( I )*WORK( 1 )
 605                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
 606                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
 607                END IF
 608             END IF
 609    30    CONTINUE
 610       END IF
 611 *
 612 *     Undo scaling
 613 *
 614       IF( ILASCL ) THEN
 615          CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
 616          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
 617          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
 618       END IF
 619 *
 620       IF( ILBSCL ) THEN
 621          CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
 622          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 623       END IF
 624 *
 625       IF( WANTST ) THEN
 626 *
 627 *        Check if reordering is correct
 628 *
 629          LASTSL = .TRUE.
 630          LST2SL = .TRUE.
 631          SDIM = 0
 632          IP = 0
 633          DO 40 I = 1, N
 634             CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
 635             IF( ALPHAI( I ).EQ.ZERO ) THEN
 636                IF( CURSL )
 637      $            SDIM = SDIM + 1
 638                IP = 0
 639                IF( CURSL .AND. .NOT.LASTSL )
 640      $            INFO = N + 2
 641             ELSE
 642                IF( IP.EQ.1 ) THEN
 643 *
 644 *                 Last eigenvalue of conjugate pair
 645 *
 646                   CURSL = CURSL .OR. LASTSL
 647                   LASTSL = CURSL
 648                   IF( CURSL )
 649      $               SDIM = SDIM + 2
 650                   IP = -1
 651                   IF( CURSL .AND. .NOT.LST2SL )
 652      $               INFO = N + 2
 653                ELSE
 654 *
 655 *                 First eigenvalue of conjugate pair
 656 *
 657                   IP = 1
 658                END IF
 659             END IF
 660             LST2SL = LASTSL
 661             LASTSL = CURSL
 662    40    CONTINUE
 663 *
 664       END IF
 665 *
 666    50 CONTINUE
 667 *
 668       WORK( 1 ) = LWKOPT
 669 *
 670       RETURN
 671 *
 672 *     End of DGGES3
 673 *
 674       END