home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / lapack / dgeesx.f < prev    next >
Text File  |  1996-09-28  |  18KB  |  493 lines
  1.       SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
  2.      $                   WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
  3.      $                   IWORK, LIWORK, BWORK, INFO )
  4. *
  5. *  -- LAPACK driver routine (version 2.0) --
  6. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  7. *     Courant Institute, Argonne National Lab, and Rice University
  8. *     March 31, 1993
  9. *
  10. *     .. Scalar Arguments ..
  11.       CHARACTER          JOBVS, SENSE, SORT
  12.       INTEGER            INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
  13.       DOUBLE PRECISION   RCONDE, RCONDV
  14. *     ..
  15. *     .. Array Arguments ..
  16.       LOGICAL            BWORK( * )
  17.       INTEGER            IWORK( * )
  18.       DOUBLE PRECISION   A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
  19.      $                   WR( * )
  20. *     ..
  21. *     .. Function Arguments ..
  22.       LOGICAL            SELECT
  23.       EXTERNAL           SELECT
  24. *     ..
  25. *
  26. *  Purpose
  27. *  =======
  28. *
  29. *  DGEESX computes for an N-by-N real nonsymmetric matrix A, the
  30. *  eigenvalues, the real Schur form T, and, optionally, the matrix of
  31. *  Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
  32. *
  33. *  Optionally, it also orders the eigenvalues on the diagonal of the
  34. *  real Schur form so that selected eigenvalues are at the top left;
  35. *  computes a reciprocal condition number for the average of the
  36. *  selected eigenvalues (RCONDE); and computes a reciprocal condition
  37. *  number for the right invariant subspace corresponding to the
  38. *  selected eigenvalues (RCONDV).  The leading columns of Z form an
  39. *  orthonormal basis for this invariant subspace.
  40. *
  41. *  For further explanation of the reciprocal condition numbers RCONDE
  42. *  and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
  43. *  these quantities are called s and sep respectively).
  44. *
  45. *  A real matrix is in real Schur form if it is upper quasi-triangular
  46. *  with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
  47. *  the form
  48. *            [  a  b  ]
  49. *            [  c  a  ]
  50. *
  51. *  where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
  52. *
  53. *  Arguments
  54. *  =========
  55. *
  56. *  JOBVS   (input) CHARACTER*1
  57. *          = 'N': Schur vectors are not computed;
  58. *          = 'V': Schur vectors are computed.
  59. *
  60. *  SORT    (input) CHARACTER*1
  61. *          Specifies whether or not to order the eigenvalues on the
  62. *          diagonal of the Schur form.
  63. *          = 'N': Eigenvalues are not ordered;
  64. *          = 'S': Eigenvalues are ordered (see SELECT).
  65. *
  66. *  SELECT  (input) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
  67. *          SELECT must be declared EXTERNAL in the calling subroutine.
  68. *          If SORT = 'S', SELECT is used to select eigenvalues to sort
  69. *          to the top left of the Schur form.
  70. *          If SORT = 'N', SELECT is not referenced.
  71. *          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
  72. *          SELECT(WR(j),WI(j)) is true; i.e., if either one of a
  73. *          complex conjugate pair of eigenvalues is selected, then both
  74. *          are.  Note that a selected complex eigenvalue may no longer
  75. *          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
  76. *          ordering may change the value of complex eigenvalues
  77. *          (especially if the eigenvalue is ill-conditioned); in this
  78. *          case INFO may be set to N+3 (see INFO below).
  79. *
  80. *  SENSE   (input) CHARACTER*1
  81. *          Determines which reciprocal condition numbers are computed.
  82. *          = 'N': None are computed;
  83. *          = 'E': Computed for average of selected eigenvalues only;
  84. *          = 'V': Computed for selected right invariant subspace only;
  85. *          = 'B': Computed for both.
  86. *          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
  87. *
  88. *  N       (input) INTEGER
  89. *          The order of the matrix A. N >= 0.
  90. *
  91. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
  92. *          On entry, the N-by-N matrix A.
  93. *          On exit, A is overwritten by its real Schur form T.
  94. *
  95. *  LDA     (input) INTEGER
  96. *          The leading dimension of the array A.  LDA >= max(1,N).
  97. *
  98. *  SDIM    (output) INTEGER
  99. *          If SORT = 'N', SDIM = 0.
  100. *          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
  101. *                         for which SELECT is true. (Complex conjugate
  102. *                         pairs for which SELECT is true for either
  103. *                         eigenvalue count as 2.)
  104. *
  105. *  WR      (output) DOUBLE PRECISION array, dimension (N)
  106. *  WI      (output) DOUBLE PRECISION array, dimension (N)
  107. *          WR and WI contain the real and imaginary parts, respectively,
  108. *          of the computed eigenvalues, in the same order that they
  109. *          appear on the diagonal of the output Schur form T.  Complex
  110. *          conjugate pairs of eigenvalues appear consecutively with the
  111. *          eigenvalue having the positive imaginary part first.
  112. *
  113. *  VS      (output) DOUBLE PRECISION array, dimension (LDVS,N)
  114. *          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
  115. *          vectors.
  116. *          If JOBVS = 'N', VS is not referenced.
  117. *
  118. *  LDVS    (input) INTEGER
  119. *          The leading dimension of the array VS.  LDVS >= 1, and if
  120. *          JOBVS = 'V', LDVS >= N.
  121. *
  122. *  RCONDE  (output) DOUBLE PRECISION
  123. *          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
  124. *          condition number for the average of the selected eigenvalues.
  125. *          Not referenced if SENSE = 'N' or 'V'.
  126. *
  127. *  RCONDV  (output) DOUBLE PRECISION
  128. *          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
  129. *          condition number for the selected right invariant subspace.
  130. *          Not referenced if SENSE = 'N' or 'E'.
  131. *
  132. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  133. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  134. *
  135. *  LWORK   (input) INTEGER
  136. *          The dimension of the array WORK.  LWORK >= max(1,3*N).
  137. *          Also, if SENSE = 'E' or 'V' or 'B',
  138. *          LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
  139. *          selected eigenvalues computed by this routine.  Note that
  140. *          N+2*SDIM*(N-SDIM) <= N+N*N/2.
  141. *          For good performance, LWORK must generally be larger.
  142. *
  143. *  IWORK   (workspace) INTEGER array, dimension (LIWORK)
  144. *          Not referenced if SENSE = 'N' or 'E'.
  145. *
  146. *  LIWORK  (input) INTEGER
  147. *          The dimension of the array IWORK.
  148. *          LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
  149. *
  150. *  BWORK   (workspace) LOGICAL array, dimension (N)
  151. *          Not referenced if SORT = 'N'.
  152. *
  153. *  INFO    (output) INTEGER
  154. *          = 0: successful exit
  155. *          < 0: if INFO = -i, the i-th argument had an illegal value.
  156. *          > 0: if INFO = i, and i is
  157. *             <= N: the QR algorithm failed to compute all the
  158. *                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
  159. *                   contain those eigenvalues which have converged; if
  160. *                   JOBVS = 'V', VS contains the transformation which
  161. *                   reduces A to its partially converged Schur form.
  162. *             = N+1: the eigenvalues could not be reordered because some
  163. *                   eigenvalues were too close to separate (the problem
  164. *                   is very ill-conditioned);
  165. *             = N+2: after reordering, roundoff changed values of some
  166. *                   complex eigenvalues so that leading eigenvalues in
  167. *                   the Schur form no longer satisfy SELECT=.TRUE.  This
  168. *                   could also be caused by underflow due to scaling.
  169. *
  170. *  =====================================================================
  171. *
  172. *     .. Parameters ..
  173.       DOUBLE PRECISION   ZERO, ONE
  174.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  175. *     ..
  176. *     .. Local Scalars ..
  177.       LOGICAL            CURSL, LASTSL, LST2SL, SCALEA, WANTSB, WANTSE,
  178.      $                   WANTSN, WANTST, WANTSV, WANTVS
  179.       INTEGER            HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
  180.      $                   IHI, ILO, INXT, IP, ITAU, IWRK, K, MAXB,
  181.      $                   MAXWRK, MINWRK
  182.       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SMLNUM
  183. *     ..
  184. *     .. Local Arrays ..
  185.       DOUBLE PRECISION   DUM( 1 )
  186. *     ..
  187. *     .. External Subroutines ..
  188.       EXTERNAL           DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD,
  189.      $                   DLACPY, DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
  190. *     ..
  191. *     .. External Functions ..
  192.       LOGICAL            LSAME
  193.       INTEGER            ILAENV
  194.       DOUBLE PRECISION   DLAMCH, DLANGE
  195.       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
  196. *     ..
  197. *     .. Intrinsic Functions ..
  198.       INTRINSIC          MAX, MIN, SQRT
  199. *     ..
  200. *     .. Executable Statements ..
  201. *
  202. *     Test the input arguments
  203. *
  204.       INFO = 0
  205.       WANTVS = LSAME( JOBVS, 'V' )
  206.       WANTST = LSAME( SORT, 'S' )
  207.       WANTSN = LSAME( SENSE, 'N' )
  208.       WANTSE = LSAME( SENSE, 'E' )
  209.       WANTSV = LSAME( SENSE, 'V' )
  210.       WANTSB = LSAME( SENSE, 'B' )
  211.       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
  212.          INFO = -1
  213.       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
  214.          INFO = -2
  215.       ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
  216.      $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
  217.          INFO = -4
  218.       ELSE IF( N.LT.0 ) THEN
  219.          INFO = -5
  220.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  221.          INFO = -7
  222.       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
  223.          INFO = -12
  224.       END IF
  225. *
  226. *     Compute workspace
  227. *      (Note: Comments in the code beginning "RWorkspace:" describe the
  228. *       minimal amount of real workspace needed at that point in the
  229. *       code, as well as the preferred amount for good performance.
  230. *       IWorkspace refers to integer workspace.
  231. *       NB refers to the optimal block size for the immediately
  232. *       following subroutine, as returned by ILAENV.
  233. *       HSWORK refers to the workspace preferred by DHSEQR, as
  234. *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
  235. *       the worst case.
  236. *       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
  237. *       depends on SDIM, which is computed by the routine DTRSEN later
  238. *       in the code.)
  239. *
  240.       MINWRK = 1
  241.       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
  242.          MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
  243.          MINWRK = MAX( 1, 3*N )
  244.          IF( .NOT.WANTVS ) THEN
  245.             MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SN', N, 1, N, -1 ), 2 )
  246.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SN', N, 1,
  247.      $          N, -1 ) ) )
  248.             HSWORK = MAX( K*( K+2 ), 2*N )
  249.             MAXWRK = MAX( MAXWRK, N+HSWORK, 1 )
  250.          ELSE
  251.             MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
  252.      $               ILAENV( 1, 'DORGHR', ' ', N, 1, N, -1 ) )
  253.             MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SV', N, 1, N, -1 ), 2 )
  254.             K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SV', N, 1,
  255.      $          N, -1 ) ) )
  256.             HSWORK = MAX( K*( K+2 ), 2*N )
  257.             MAXWRK = MAX( MAXWRK, N+HSWORK, 1 )
  258.          END IF
  259.          WORK( 1 ) = MAXWRK
  260.       END IF
  261.       IF( LWORK.LT.MINWRK ) THEN
  262.          INFO = -16
  263.       END IF
  264.       IF( INFO.NE.0 ) THEN
  265.          CALL XERBLA( 'DGEESX', -INFO )
  266.          RETURN
  267.       END IF
  268. *
  269. *     Quick return if possible
  270. *
  271.       IF( N.EQ.0 ) THEN
  272.          SDIM = 0
  273.          RETURN
  274.       END IF
  275. *
  276. *     Get machine constants
  277. *
  278.       EPS = DLAMCH( 'P' )
  279.       SMLNUM = DLAMCH( 'S' )
  280.       BIGNUM = ONE / SMLNUM
  281.       CALL DLABAD( SMLNUM, BIGNUM )
  282.       SMLNUM = SQRT( SMLNUM ) / EPS
  283.       BIGNUM = ONE / SMLNUM
  284. *
  285. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  286. *
  287.       ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
  288.       SCALEA = .FALSE.
  289.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  290.          SCALEA = .TRUE.
  291.          CSCALE = SMLNUM
  292.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  293.          SCALEA = .TRUE.
  294.          CSCALE = BIGNUM
  295.       END IF
  296.       IF( SCALEA )
  297.      $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
  298. *
  299. *     Permute the matrix to make it more nearly triangular
  300. *     (RWorkspace: need N)
  301. *
  302.       IBAL = 1
  303.       CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
  304. *
  305. *     Reduce to upper Hessenberg form
  306. *     (RWorkspace: need 3*N, prefer 2*N+N*NB)
  307. *
  308.       ITAU = N + IBAL
  309.       IWRK = N + ITAU
  310.       CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
  311.      $             LWORK-IWRK+1, IERR )
  312. *
  313.       IF( WANTVS ) THEN
  314. *
  315. *        Copy Householder vectors to VS
  316. *
  317.          CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
  318. *
  319. *        Generate orthogonal matrix in VS
  320. *        (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
  321. *
  322.          CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
  323.      $                LWORK-IWRK+1, IERR )
  324.       END IF
  325. *
  326.       SDIM = 0
  327. *
  328. *     Perform QR iteration, accumulating Schur vectors in VS if desired
  329. *     (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
  330. *
  331.       IWRK = ITAU
  332.       CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
  333.      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
  334.       IF( IEVAL.GT.0 )
  335.      $   INFO = IEVAL
  336. *
  337. *     Sort eigenvalues if desired
  338. *
  339.       IF( WANTST .AND. INFO.EQ.0 ) THEN
  340.          IF( SCALEA ) THEN
  341.             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
  342.             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
  343.          END IF
  344.          DO 10 I = 1, N
  345.             BWORK( I ) = SELECT( WR( I ), WI( I ) )
  346.    10    CONTINUE
  347. *
  348. *        Reorder eigenvalues, transform Schur vectors, and compute
  349. *        reciprocal condition numbers
  350. *        (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
  351. *                     otherwise, need N )
  352. *        (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
  353. *                     otherwise, need 0 )
  354. *
  355.          CALL DTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
  356.      $                SDIM, RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
  357.      $                IWORK, LIWORK, ICOND )
  358.          IF( .NOT.WANTSN )
  359.      $      MAXWRK = MAX( MAXWRK, N+2*SDIM*( N-SDIM ) )
  360.          IF( ICOND.EQ.-15 ) THEN
  361. *
  362. *           Not enough real workspace
  363. *
  364.             INFO = -16
  365.          ELSE IF( ICOND.EQ.-17 ) THEN
  366. *
  367. *           Not enough integer workspace
  368. *
  369.             INFO = -18
  370.          ELSE IF( ICOND.GT.0 ) THEN
  371. *
  372. *           DTRSEN failed to reorder or to restore standard Schur form
  373. *
  374.             INFO = ICOND + N
  375.          END IF
  376.       END IF
  377. *
  378.       IF( WANTVS ) THEN
  379. *
  380. *        Undo balancing
  381. *        (RWorkspace: need N)
  382. *
  383.          CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
  384.      $                IERR )
  385.       END IF
  386. *
  387.       IF( SCALEA ) THEN
  388. *
  389. *        Undo scaling for the Schur form of A
  390. *
  391.          CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
  392.          CALL DCOPY( N, A, LDA+1, WR, 1 )
  393.          IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
  394.             DUM( 1 ) = RCONDV
  395.             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
  396.             RCONDV = DUM( 1 )
  397.          END IF
  398.          IF( CSCALE.EQ.SMLNUM ) THEN
  399. *
  400. *           If scaling back towards underflow, adjust WI if an
  401. *           offdiagonal element of a 2-by-2 block in the Schur form
  402. *           underflows.
  403. *
  404.             IF( IEVAL.GT.0 ) THEN
  405.                I1 = IEVAL + 1
  406.                I2 = IHI - 1
  407.                CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
  408.      $                      IERR )
  409.             ELSE IF( WANTST ) THEN
  410.                I1 = 1
  411.                I2 = N - 1
  412.             ELSE
  413.                I1 = ILO
  414.                I2 = IHI - 1
  415.             END IF
  416.             INXT = I1 - 1
  417.             DO 20 I = I1, I2
  418.                IF( I.LT.INXT )
  419.      $            GO TO 20
  420.                IF( WI( I ).EQ.ZERO ) THEN
  421.                   INXT = I + 1
  422.                ELSE
  423.                   IF( A( I+1, I ).EQ.ZERO ) THEN
  424.                      WI( I ) = ZERO
  425.                      WI( I+1 ) = ZERO
  426.                   ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
  427.      $                     ZERO ) THEN
  428.                      WI( I ) = ZERO
  429.                      WI( I+1 ) = ZERO
  430.                      IF( I.GT.1 )
  431.      $                  CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
  432.                      IF( N.GT.I+1 )
  433.      $                  CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
  434.      $                              A( I+1, I+2 ), LDA )
  435.                      CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
  436.                      A( I, I+1 ) = A( I+1, I )
  437.                      A( I+1, I ) = ZERO
  438.                   END IF
  439.                   INXT = I + 2
  440.                END IF
  441.    20       CONTINUE
  442.          END IF
  443.          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
  444.      $                WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
  445.       END IF
  446. *
  447.       IF( WANTST .AND. INFO.EQ.0 ) THEN
  448. *
  449. *        Check if reordering successful
  450. *
  451.          LASTSL = .TRUE.
  452.          LST2SL = .TRUE.
  453.          SDIM = 0
  454.          IP = 0
  455.          DO 30 I = 1, N
  456.             CURSL = SELECT( WR( I ), WI( I ) )
  457.             IF( WI( I ).EQ.ZERO ) THEN
  458.                IF( CURSL )
  459.      $            SDIM = SDIM + 1
  460.                IP = 0
  461.                IF( CURSL .AND. .NOT.LASTSL )
  462.      $            INFO = N + 2
  463.             ELSE
  464.                IF( IP.EQ.1 ) THEN
  465. *
  466. *                 Last eigenvalue of conjugate pair
  467. *
  468.                   CURSL = CURSL .OR. LASTSL
  469.                   LASTSL = CURSL
  470.                   IF( CURSL )
  471.      $               SDIM = SDIM + 2
  472.                   IP = -1
  473.                   IF( CURSL .AND. .NOT.LST2SL )
  474.      $               INFO = N + 2
  475.                ELSE
  476. *
  477. *                 First eigenvalue of conjugate pair
  478. *
  479.                   IP = 1
  480.                END IF
  481.             END IF
  482.             LST2SL = LASTSL
  483.             LASTSL = CURSL
  484.    30    CONTINUE
  485.       END IF
  486. *
  487.       WORK( 1 ) = MAXWRK
  488.       RETURN
  489. *
  490. *     End of DGEESX
  491. *
  492.       END
  493.