SRC/dggesx.f

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