      SUBROUTINE DS3GRD(ALPHA, B, D, ETA0, FX, G, IRC, P, W, X)
C
C  ***  COMPUTE FINITE DIFFERENCE GRADIENT BY STWEART*S SCHEME  ***
C
C     ***  PARAMETERS  ***
C
      INTEGER IRC, P
      DOUBLE PRECISION ALPHA(P), B(2,P), D(P), ETA0, FX, G(P), W(6),
     1                 X(P)
C
C.......................................................................
C
C     ***  PURPOSE  ***
C
C        THIS SUBROUTINE USES AN EMBELLISHED FORM OF THE FINITE-DIFFER-
C     ENCE SCHEME PROPOSED BY STEWART (REF. 1) TO APPROXIMATE THE
C     GRADIENT OF THE FUNCTION F(X), WHOSE VALUES ARE SUPPLIED BY
C     REVERSE COMMUNICATION.
C
C     ***  PARAMETER DESCRIPTION  ***
C
C  ALPHA IN  (APPROXIMATE) DIAGONAL ELEMENTS OF THE HESSIAN OF F(X).
C      B IN  ARRAY OF SIMPLE LOWER AND UPPER BOUNDS ON X.  X MUST
C             SATISFY B(1,I) .LE. X(I) .LE. B(2,I), I = 1(1)P.
C             FOR ALL I WITH B(1,I) .GE. B(2,I), DS3GRD SIMPLY
C             SETS G(I) TO 0.
C      D IN  SCALE VECTOR SUCH THAT D(I)*X(I), I = 1,...,P, ARE IN
C             COMPARABLE UNITS.
C   ETA0 IN  ESTIMATED BOUND ON RELATIVE ERROR IN THE FUNCTION VALUE...
C             (TRUE VALUE) = (COMPUTED VALUE)*(1+E),   WHERE
C             ABS(E) .LE. ETA0.
C     FX I/O ON INPUT,  FX  MUST BE THE COMPUTED VALUE OF F(X).  ON
C             OUTPUT WITH IRC = 0, FX HAS BEEN RESTORED TO ITS ORIGINAL
C             VALUE, THE ONE IT HAD WHEN DS3GRD WAS LAST CALLED WITH
C             IRC = 0.
C      G I/O ON INPUT WITH IRC = 0, G SHOULD CONTAIN AN APPROXIMATION
C             TO THE GRADIENT OF F NEAR X, E.G., THE GRADIENT AT THE
C             PREVIOUS ITERATE.  WHEN DS3GRD RETURNS WITH IRC = 0, G IS
C             THE DESIRED FINITE-DIFFERENCE APPROXIMATION TO THE
C             GRADIENT AT X.
C    IRC I/O INPUT/RETURN CODE... BEFORE THE VERY FIRST CALL ON DS3GRD,
C             THE CALLER MUST SET IRC TO 0.  WHENEVER DS3GRD RETURNS A
C             NONZERO VALUE (OF AT MOST P) FOR IRC, IT HAS PERTURBED
C             SOME COMPONENT OF X... THE CALLER SHOULD EVALUATE F(X)
C             AND CALL DS3GRD AGAIN WITH FX = F(X).  IF B PREVENTS
C             ESTIMATING G(I) I.E., IF THERE IS AN I WITH
C             B(1,I) .LT. B(2,I) BUT WITH B(1,I) SO CLOSE TO B(2,I)
C             THAT THE FINITE-DIFFERENCING STEPS CANNOT BE CHOSEN,
C             THEN DS3GRD RETURNS WITH IRC .GT. P.
C      P IN  THE NUMBER OF VARIABLES (COMPONENTS OF X) ON WHICH F
C             DEPENDS.
C      X I/O ON INPUT WITH IRC = 0, X IS THE POINT AT WHICH THE
C             GRADIENT OF F IS DESIRED.  ON OUTPUT WITH IRC NONZERO, X
C             IS THE POINT AT WHICH F SHOULD BE EVALUATED.  ON OUTPUT
C             WITH IRC = 0, X HAS BEEN RESTORED TO ITS ORIGINAL VALUE
C             (THE ONE IT HAD WHEN DS3GRD WAS LAST CALLED WITH IRC = 0)
C             AND G CONTAINS THE DESIRED GRADIENT APPROXIMATION.
C      W I/O WORK VECTOR OF LENGTH 6 IN WHICH DS3GRD SAVES CERTAIN
C             QUANTITIES WHILE THE CALLER IS EVALUATING F(X) AT A
C             PERTURBED X.
C
C     ***  APPLICATION AND USAGE RESTRICTIONS  ***
C
C        THIS ROUTINE IS INTENDED FOR USE WITH QUASI-NEWTON ROUTINES
C     FOR UNCONSTRAINED MINIMIZATION (IN WHICH CASE  ALPHA  COMES FROM
C     THE DIAGONAL OF THE QUASI-NEWTON HESSIAN APPROXIMATION).
C
C     ***  ALGORITHM NOTES  ***
C
C        THIS CODE DEPARTS FROM THE SCHEME PROPOSED BY STEWART (REF. 1)
C     IN ITS GUARDING AGAINST OVERLY LARGE OR SMALL STEP SIZES AND ITS
C     HANDLING OF SPECIAL CASES (SUCH AS ZERO COMPONENTS OF ALPHA OR G).
C
C     ***  REFERENCES  ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C        METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C        J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C
C     ***  HISTORY  ***
C
C     DESIGNED AND CODED BY DAVID M. GAY (SUMMER 1977/SUMMER 1980).
C
C     ***  GENERAL  ***
C
C        THIS ROUTINE WAS PREPARED IN CONNECTION WITH WORK SUPPORTED BY
C     THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS76-00324 AND
C     MCS-7906671.
C
C.......................................................................
C
C     *****  EXTERNAL FUNCTION  *****
C
      DOUBLE PRECISION DR7MDC
      EXTERNAL DR7MDC
C DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS.
C
C     ***** INTRINSIC FUNCTIONS *****
C/+
      DOUBLE PRECISION DSQRT
C/
C     ***** LOCAL VARIABLES *****
C
      LOGICAL HIT
      INTEGER FH, FX0, HSAVE, I, XISAVE
      DOUBLE PRECISION AAI, AFX, AFXETA, AGI, ALPHAI, AXI, AXIBAR,
     1                 DISCON, ETA, GI, H, HMIN, XI, XIH
      DOUBLE PRECISION C2000, FOUR, HMAX0, HMIN0, H0, MACHEP, ONE, P002,
     1                 THREE, TWO, ZERO
C
C/6
C     DATA C2000/2.0D+3/, FOUR/4.0D+0/, HMAX0/0.02D+0/, HMIN0/5.0D+1/,
C    1     ONE/1.0D+0/, P002/0.002D+0/, THREE/3.0D+0/,
C    2     TWO/2.0D+0/, ZERO/0.0D+0/
C/7
      PARAMETER (C2000=2.0D+3, FOUR=4.0D+0, HMAX0=0.02D+0, HMIN0=5.0D+1,
     1     ONE=1.0D+0, P002=0.002D+0, THREE=3.0D+0,
     2     TWO=2.0D+0, ZERO=0.0D+0)
C/
C/6
C     DATA FH/3/, FX0/4/, HSAVE/5/, XISAVE/6/
C/7
      PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6)
C/
C
C---------------------------------  BODY  ------------------------------
C
      IF (IRC) 80, 10, 210
C
C     ***  FRESH START -- GET MACHINE-DEPENDENT CONSTANTS  ***
C
C     STORE MACHEP IN W(1) AND H0 IN W(2), WHERE MACHEP IS THE UNIT
C     ROUNDOFF (THE SMALLEST POSITIVE NUMBER SUCH THAT
C     1 + MACHEP .GT. 1  AND  1 - MACHEP .LT. 1),  AND  H0 IS THE
C     SQUARE-ROOT OF MACHEP.
C
 10   W(1) = DR7MDC(3)
      W(2) = DSQRT(W(1))
C
      W(FX0) = FX
C
C     ***  INCREMENT  I  AND START COMPUTING  G(I)  ***
C
 20   I = IABS(IRC) + 1
      IF (I .GT. P) GO TO 220
         IRC = I
         IF (B(1,I) .LT. B(2,I)) GO TO 30
            G(I) = ZERO
            GO TO 20
 30      AFX = DABS(W(FX0))
         MACHEP = W(1)
         H0 = W(2)
         HMIN = HMIN0 * MACHEP
         XI = X(I)
         W(XISAVE) = XI
         AXI = DABS(XI)
         AXIBAR = DMAX1(AXI, ONE/D(I))
         GI = G(I)
         AGI = DABS(GI)
         ETA = DABS(ETA0)
         IF (AFX .GT. ZERO) ETA = DMAX1(ETA, AGI*AXI*MACHEP/AFX)
         ALPHAI = ALPHA(I)
         IF (ALPHAI .EQ. ZERO) GO TO 130
         IF (GI .EQ. ZERO .OR. FX .EQ. ZERO) GO TO 140
         AFXETA = AFX*ETA
         AAI = DABS(ALPHAI)
C
C        *** COMPUTE H = STEWART*S FORWARD-DIFFERENCE STEP SIZE.
C
         IF (GI**2 .LE. AFXETA*AAI) GO TO 40
              H = TWO*DSQRT(AFXETA/AAI)
              H = H*(ONE - AAI*H/(THREE*AAI*H + FOUR*AGI))
              GO TO 50
C40      H = TWO*(AFXETA*AGI/(AAI**2))**(ONE/THREE)
 40      H = TWO * (AFXETA*AGI)**(ONE/THREE) * AAI**(-TWO/THREE)
         H = H*(ONE - TWO*AGI/(THREE*AAI*H + FOUR*AGI))
C
C        ***  ENSURE THAT  H  IS NOT INSIGNIFICANTLY SMALL  ***
C
 50      H = DMAX1(H, HMIN*AXIBAR)
C
C        *** USE FORWARD DIFFERENCE IF BOUND ON TRUNCATION ERROR IS AT
C        *** MOST 10**-3.
C
         IF (AAI*H .LE. P002*AGI) GO TO 120
C
C        *** COMPUTE H = STEWART*S STEP FOR CENTRAL DIFFERENCE.
C
         DISCON = C2000*AFXETA
         H = DISCON/(AGI + DSQRT(GI**2 + AAI*DISCON))
C
C        ***  ENSURE THAT  H  IS NEITHER TOO SMALL NOR TOO BIG  ***
C
         H = DMAX1(H, HMIN*AXIBAR)
         IF (H .GE. HMAX0*AXIBAR) H = AXIBAR * H0**(TWO/THREE)
C
C        ***  COMPUTE CENTRAL DIFFERENCE  ***
C
         XIH = XI + H
         IF (XI - H .LT. B(1,I)) GO TO 60
         IRC = -I
         IF (XIH .LE. B(2,I)) GO TO 200
            H = -H
            XIH = XI + H
            IF (XI + TWO*H .LT. B(1,I)) GO TO 190
            GO TO 70
 60      IF (XI + TWO*H .GT. B(2,I)) GO TO 190
C        *** MUST DO OFF-SIDE CENTRAL DIFFERENCE ***
 70      IRC = -(I + P)
         GO TO 200
C
 80      I = -IRC
         IF (I .LE. P) GO TO 100
         I = I - P
         IF (I .GT. P) GO TO 90
         W(FH) = FX
         H = TWO * W(HSAVE)
         XIH = W(XISAVE) + H
         IRC = IRC - P
         GO TO 200
C
C    *** FINISH OFF-SIDE CENTRAL DIFFERENCE ***
C
 90      I = I - P
         G(I) = (FOUR*W(FH) - FX - THREE*W(FX0)) / W(HSAVE)
         IRC = I
         X(I) = W(XISAVE)
         GO TO 20
C
 100     H = -W(HSAVE)
         IF (H .GT. ZERO) GO TO 110
         W(FH) = FX
         XIH = W(XISAVE) + H
         GO TO 200
C
 110     G(I) = (W(FH) - FX) / (TWO * H)
         X(I) = W(XISAVE)
         GO TO 20
C
C     ***  COMPUTE FORWARD DIFFERENCES IN VARIOUS CASES  ***
C
 120     IF (H .GE. HMAX0*AXIBAR) H = H0 * AXIBAR
         IF (ALPHAI*GI .LT. ZERO) H = -H
         GO TO 150
 130     H = AXIBAR
         GO TO 150
 140     H = H0 * AXIBAR
C
 150     HIT = .FALSE.
 160     XIH = XI + H
         IF (H .GT. ZERO) GO TO 170
            IF (XIH .GE. B(1,I)) GO TO 200
            GO TO 180
 170     IF (XIH .LE. B(2,I)) GO TO 200
 180        IF (HIT) GO TO 190
            HIT = .TRUE.
            H = -H
            GO TO 160
C
C        *** ERROR RETURN...
 190     IRC = I + P
         GO TO 230
C
C        *** RETURN FOR NEW FUNCTION VALUE...
 200     X(I) = XIH
         W(HSAVE) = H
         GO TO 999
C
C     ***  COMPUTE ACTUAL FORWARD DIFFERENCE  ***
C
 210     G(IRC) = (FX - W(FX0)) / W(HSAVE)
         X(IRC) = W(XISAVE)
         GO TO 20
C
C  ***  RESTORE FX AND INDICATE THAT G HAS BEEN COMPUTED  ***
C
 220  IRC = 0
 230  FX = W(FX0)
C
 999  RETURN
C  ***  LAST LINE OF DS3GRD FOLLOWS  ***
      END
