LAPACK: SRC/dggesx.f Source File

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