SRC/dgges.f

   1 *> \brief <b> DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGGES + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgges.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgges.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgges.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGGES( 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 *> DGGES 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 *>          If N = 0, LWORK >= 1, else LWORK >= 8*N+16.
 238 *>          For good performance , LWORK must generally be larger.
 239 *>
 240 *>          If LWORK = -1, then a workspace query is assumed; the routine
 241 *>          only calculates the optimal size of the WORK array, returns
 242 *>          this value as the first entry of the WORK array, and no error
 243 *>          message related to LWORK is issued by XERBLA.
 244 *> \endverbatim
 245 *>
 246 *> \param[out] BWORK
 247 *> \verbatim
 248 *>          BWORK is LOGICAL array, dimension (N)
 249 *>          Not referenced if SORT = 'N'.
 250 *> \endverbatim
 251 *>
 252 *> \param[out] INFO
 253 *> \verbatim
 254 *>          INFO is INTEGER
 255 *>          = 0:  successful exit
 256 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 257 *>          = 1,...,N:
 258 *>                The QZ iteration failed.  (A,B) are not in Schur
 259 *>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
 260 *>                be correct for j=INFO+1,...,N.
 261 *>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
 262 *>                =N+2: after reordering, roundoff changed values of
 263 *>                      some complex eigenvalues so that leading
 264 *>                      eigenvalues in the Generalized Schur form no
 265 *>                      longer satisfy SELCTG=.TRUE.  This could also
 266 *>                      be caused due to scaling.
 267 *>                =N+3: reordering failed in DTGSEN.
 268 *> \endverbatim
 269 *
 270 *  Authors:
 271 *  ========
 272 *
 273 *> \author Univ. of Tennessee
 274 *> \author Univ. of California Berkeley
 275 *> \author Univ. of Colorado Denver
 276 *> \author NAG Ltd.
 277 *
 278 *> \date November 2011
 279 *
 280 *> \ingroup doubleGEeigen
 281 *
 282 *  =====================================================================
 283       SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
 284      $                  SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
 285      $                  LDVSR, WORK, LWORK, BWORK, INFO )
 286 *
 287 *  -- LAPACK driver routine (version 3.4.0) --
 288 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 289 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 290 *     November 2011
 291 *
 292 *     .. Scalar Arguments ..
 293       CHARACTER          JOBVSL, JOBVSR, SORT
 294       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
 295 *     ..
 296 *     .. Array Arguments ..
 297       LOGICAL            BWORK( * )
 298       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
 299      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
 300      $                   VSR( LDVSR, * ), WORK( * )
 301 *     ..
 302 *     .. Function Arguments ..
 303       LOGICAL            SELCTG
 304       EXTERNAL           SELCTG
 305 *     ..
 306 *
 307 *  =====================================================================
 308 *
 309 *     .. Parameters ..
 310       DOUBLE PRECISION   ZERO, ONE
 311       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 312 *     ..
 313 *     .. Local Scalars ..
 314       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
 315      $                   LQUERY, LST2SL, WANTST
 316       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
 317      $                   ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
 318      $                   MINWRK
 319       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
 320      $                   PVSR, SAFMAX, SAFMIN, SMLNUM
 321 *     ..
 322 *     .. Local Arrays ..
 323       INTEGER            IDUM( 1 )
 324       DOUBLE PRECISION   DIF( 2 )
 325 *     ..
 326 *     .. External Subroutines ..
 327       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
 328      $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
 329      $                   XERBLA
 330 *     ..
 331 *     .. External Functions ..
 332       LOGICAL            LSAME
 333       INTEGER            ILAENV
 334       DOUBLE PRECISION   DLAMCH, DLANGE
 335       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
 336 *     ..
 337 *     .. Intrinsic Functions ..
 338       INTRINSIC          ABS, MAX, SQRT
 339 *     ..
 340 *     .. Executable Statements ..
 341 *
 342 *     Decode the input arguments
 343 *
 344       IF( LSAME( JOBVSL, 'N' ) ) THEN
 345          IJOBVL = 1
 346          ILVSL = .FALSE.
 347       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
 348          IJOBVL = 2
 349          ILVSL = .TRUE.
 350       ELSE
 351          IJOBVL = -1
 352          ILVSL = .FALSE.
 353       END IF
 354 *
 355       IF( LSAME( JOBVSR, 'N' ) ) THEN
 356          IJOBVR = 1
 357          ILVSR = .FALSE.
 358       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
 359          IJOBVR = 2
 360          ILVSR = .TRUE.
 361       ELSE
 362          IJOBVR = -1
 363          ILVSR = .FALSE.
 364       END IF
 365 *
 366       WANTST = LSAME( SORT, 'S' )
 367 *
 368 *     Test the input arguments
 369 *
 370       INFO = 0
 371       LQUERY = ( LWORK.EQ.-1 )
 372       IF( IJOBVL.LE.0 ) THEN
 373          INFO = -1
 374       ELSE IF( IJOBVR.LE.0 ) THEN
 375          INFO = -2
 376       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
 377          INFO = -3
 378       ELSE IF( N.LT.0 ) THEN
 379          INFO = -5
 380       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 381          INFO = -7
 382       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 383          INFO = -9
 384       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
 385          INFO = -15
 386       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
 387          INFO = -17
 388       END IF
 389 *
 390 *     Compute workspace
 391 *      (Note: Comments in the code beginning "Workspace:" describe the
 392 *       minimal amount of workspace needed at that point in the code,
 393 *       as well as the preferred amount for good performance.
 394 *       NB refers to the optimal block size for the immediately
 395 *       following subroutine, as returned by ILAENV.)
 396 *
 397       IF( INFO.EQ.0 ) THEN
 398          IF( N.GT.0 )THEN
 399             MINWRK = MAX( 8*N, 6*N + 16 )
 400             MAXWRK = MINWRK - N +
 401      $               N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
 402             MAXWRK = MAX( MAXWRK, MINWRK - N +
 403      $                    N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
 404             IF( ILVSL ) THEN
 405                MAXWRK = MAX( MAXWRK, MINWRK - N +
 406      $                       N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
 407             END IF
 408          ELSE
 409             MINWRK = 1
 410             MAXWRK = 1
 411          END IF
 412          WORK( 1 ) = MAXWRK
 413 *
 414          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
 415      $      INFO = -19
 416       END IF
 417 *
 418       IF( INFO.NE.0 ) THEN
 419          CALL XERBLA( 'DGGES ', -INFO )
 420          RETURN
 421       ELSE IF( LQUERY ) THEN
 422          RETURN
 423       END IF
 424 *
 425 *     Quick return if possible
 426 *
 427       IF( N.EQ.0 ) THEN
 428          SDIM = 0
 429          RETURN
 430       END IF
 431 *
 432 *     Get machine constants
 433 *
 434       EPS = DLAMCH( 'P' )
 435       SAFMIN = DLAMCH( 'S' )
 436       SAFMAX = ONE / SAFMIN
 437       CALL DLABAD( SAFMIN, SAFMAX )
 438       SMLNUM = SQRT( SAFMIN ) / EPS
 439       BIGNUM = ONE / SMLNUM
 440 *
 441 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 442 *
 443       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
 444       ILASCL = .FALSE.
 445       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 446          ANRMTO = SMLNUM
 447          ILASCL = .TRUE.
 448       ELSE IF( ANRM.GT.BIGNUM ) THEN
 449          ANRMTO = BIGNUM
 450          ILASCL = .TRUE.
 451       END IF
 452       IF( ILASCL )
 453      $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
 454 *
 455 *     Scale B if max element outside range [SMLNUM,BIGNUM]
 456 *
 457       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
 458       ILBSCL = .FALSE.
 459       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
 460          BNRMTO = SMLNUM
 461          ILBSCL = .TRUE.
 462       ELSE IF( BNRM.GT.BIGNUM ) THEN
 463          BNRMTO = BIGNUM
 464          ILBSCL = .TRUE.
 465       END IF
 466       IF( ILBSCL )
 467      $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
 468 *
 469 *     Permute the matrix to make it more nearly triangular
 470 *     (Workspace: need 6*N + 2*N space for storing balancing factors)
 471 *
 472       ILEFT = 1
 473       IRIGHT = N + 1
 474       IWRK = IRIGHT + N
 475       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
 476      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
 477 *
 478 *     Reduce B to triangular form (QR decomposition of B)
 479 *     (Workspace: need N, prefer N*NB)
 480 *
 481       IROWS = IHI + 1 - ILO
 482       ICOLS = N + 1 - ILO
 483       ITAU = IWRK
 484       IWRK = ITAU + IROWS
 485       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
 486      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 487 *
 488 *     Apply the orthogonal transformation to matrix A
 489 *     (Workspace: need N, prefer N*NB)
 490 *
 491       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
 492      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
 493      $             LWORK+1-IWRK, IERR )
 494 *
 495 *     Initialize VSL
 496 *     (Workspace: need N, prefer N*NB)
 497 *
 498       IF( ILVSL ) THEN
 499          CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
 500          IF( IROWS.GT.1 ) THEN
 501             CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
 502      $                   VSL( ILO+1, ILO ), LDVSL )
 503          END IF
 504          CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
 505      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
 506       END IF
 507 *
 508 *     Initialize VSR
 509 *
 510       IF( ILVSR )
 511      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
 512 *
 513 *     Reduce to generalized Hessenberg form
 514 *     (Workspace: none needed)
 515 *
 516       CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
 517      $             LDVSL, VSR, LDVSR, IERR )
 518 *
 519 *     Perform QZ algorithm, computing Schur vectors if desired
 520 *     (Workspace: need N)
 521 *
 522       IWRK = ITAU
 523       CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
 524      $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
 525      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
 526       IF( IERR.NE.0 ) THEN
 527          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
 528             INFO = IERR
 529          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
 530             INFO = IERR - N
 531          ELSE
 532             INFO = N + 1
 533          END IF
 534          GO TO 50
 535       END IF
 536 *
 537 *     Sort eigenvalues ALPHA/BETA if desired
 538 *     (Workspace: need 4*N+16 )
 539 *
 540       SDIM = 0
 541       IF( WANTST ) THEN
 542 *
 543 *        Undo scaling on eigenvalues before SELCTGing
 544 *
 545          IF( ILASCL ) THEN
 546             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
 547      $                   IERR )
 548             CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
 549      $                   IERR )
 550          END IF
 551          IF( ILBSCL )
 552      $      CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 553 *
 554 *        Select eigenvalues
 555 *
 556          DO 10 I = 1, N
 557             BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
 558    10    CONTINUE
 559 *
 560          CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
 561      $                ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
 562      $                PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
 563      $                IERR )
 564          IF( IERR.EQ.1 )
 565      $      INFO = N + 3
 566 *
 567       END IF
 568 *
 569 *     Apply back-permutation to VSL and VSR
 570 *     (Workspace: none needed)
 571 *
 572       IF( ILVSL )
 573      $   CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
 574      $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
 575 *
 576       IF( ILVSR )
 577      $   CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
 578      $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
 579 *
 580 *     Check if unscaling would cause over/underflow, if so, rescale
 581 *     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
 582 *     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
 583 *
 584       IF( ILASCL ) THEN
 585          DO 20 I = 1, N
 586             IF( ALPHAI( I ).NE.ZERO ) THEN
 587                IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
 588      $             ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
 589                   WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
 590                   BETA( I ) = BETA( I )*WORK( 1 )
 591                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
 592                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
 593                ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
 594      $                  ( ANRMTO / ANRM ) .OR.
 595      $                  ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
 596      $                   THEN
 597                   WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
 598                   BETA( I ) = BETA( I )*WORK( 1 )
 599                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
 600                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
 601                END IF
 602             END IF
 603    20    CONTINUE
 604       END IF
 605 *
 606       IF( ILBSCL ) THEN
 607          DO 30 I = 1, N
 608             IF( ALPHAI( I ).NE.ZERO ) THEN
 609                IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
 610      $             ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
 611                   WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
 612                   BETA( I ) = BETA( I )*WORK( 1 )
 613                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
 614                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
 615                END IF
 616             END IF
 617    30    CONTINUE
 618       END IF
 619 *
 620 *     Undo scaling
 621 *
 622       IF( ILASCL ) THEN
 623          CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
 624          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
 625          CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
 626       END IF
 627 *
 628       IF( ILBSCL ) THEN
 629          CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
 630          CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
 631       END IF
 632 *
 633       IF( WANTST ) THEN
 634 *
 635 *        Check if reordering is correct
 636 *
 637          LASTSL = .TRUE.
 638          LST2SL = .TRUE.
 639          SDIM = 0
 640          IP = 0
 641          DO 40 I = 1, N
 642             CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
 643             IF( ALPHAI( I ).EQ.ZERO ) THEN
 644                IF( CURSL )
 645      $            SDIM = SDIM + 1
 646                IP = 0
 647                IF( CURSL .AND. .NOT.LASTSL )
 648      $            INFO = N + 2
 649             ELSE
 650                IF( IP.EQ.1 ) THEN
 651 *
 652 *                 Last eigenvalue of conjugate pair
 653 *
 654                   CURSL = CURSL .OR. LASTSL
 655                   LASTSL = CURSL
 656                   IF( CURSL )
 657      $               SDIM = SDIM + 2
 658                   IP = -1
 659                   IF( CURSL .AND. .NOT.LST2SL )
 660      $               INFO = N + 2
 661                ELSE
 662 *
 663 *                 First eigenvalue of conjugate pair
 664 *
 665                   IP = 1
 666                END IF
 667             END IF
 668             LST2SL = LASTSL
 669             LASTSL = CURSL
 670    40    CONTINUE
 671 *
 672       END IF
 673 *
 674    50 CONTINUE
 675 *
 676       WORK( 1 ) = MAXWRK
 677 *
 678       RETURN
 679 *
 680 *     End of DGGES
 681 *
 682       END