SUBROUTINE  F7DHB(B, D, G, IRT, IV, LIV, LV, P, V, X)
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING
C  ***  AT V(IV(FDH)) = V(-IV(H)).  HONOR SIMPLE BOUNDS IN B.
C
C  ***  IF IV(COVREQ) .GE. 0 THEN  F7DHB USES GRADIENT DIFFERENCES,
C  ***  OTHERWISE FUNCTION DIFFERENCES.  STORAGE IN V IS AS IN G7LIT.
C
C IRT VALUES...
C     1 = COMPUTE FUNCTION VALUE, I.E., V(F).
C     2 = COMPUTE G.
C     3 = DONE.
C
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER IRT, LIV, LV, P
      INTEGER IV(LIV)
      REAL B(2,P), D(P), G(P), V(LV), X(P)
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL OFFSID
      INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2,
     1        NEWM1, PP1O2, STPI, STPM, STP0
      REAL DEL, DEL0, T, XM, XM1
      REAL HALF, HLIM, ONE, TWO, ZERO
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL V7CPY,  V7SCP
C
C V7CPY.... COPY ONE VECTOR TO ANOTHER.
C  V7SCP... COPY SCALAR TO ALL COMPONENTS OF A VECTOR.
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H, KAGQT, MODE,
     1        NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE
C
C/6
C     DATA HALF/0.5E+0/, HLIM/0.1E+0/, ONE/1.E+0/, TWO/2.E+0/,
C    1     ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, HLIM=0.1E+0, ONE=1.E+0, TWO=2.E+0,
     1           ZERO=0.E+0)
C/
C
C/6
C     DATA COVREQ/15/, DELTA/52/, DELTA0/44/, DLTFDC/42/, F/10/,
C    1     FDH/74/, FX/53/, H/56/, KAGQT/33/, MODE/35/, NFGCAL/7/,
C    2     SAVEI/63/, SWITCH/12/, TOOBIG/2/, W/65/, XMSAVE/51/
C/7
      PARAMETER (COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42, F=10,
     1           FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7,
     2           SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51)
C/
C
C+++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++++
C
      IRT = 4
      KIND = IV(COVREQ)
      M = IV(MODE)
      IF (M .GT. 0) GO TO 10
         HES = IABS(IV(H))
         IV(H) = -HES
         IV(FDH) = 0
         IV(KAGQT) = -1
         V(FX) = V(F)
C        *** SUPPLY ZEROS IN CASE B(1,I) = B(2,I) FOR SOME I ***
         CALL  V7SCP(P*(P+1)/2, V(HES), ZERO)
 10   IF (M .GT. P) GO TO 999
      IF (KIND .LT. 0) GO TO 120
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND
C  ***  GRADIENT VALUES.
C
      GSAVE1 = IV(W) + P
      IF (M .GT. 0) GO TO 20
C        ***  FIRST CALL ON  F7DHB.  SET GSAVE = G, TAKE FIRST STEP  ***
         CALL V7CPY(P, V(GSAVE1), G)
         IV(SWITCH) = IV(NFGCAL)
         GO TO 80
C
 20   DEL = V(DELTA)
      X(M) = V(XMSAVE)
      IF (IV(TOOBIG) .EQ. 0) GO TO 30
C
C     ***  HANDLE OVERSIZE V(DELTA)  ***
C
         DEL0 = V(DELTA0) * AMAX1(ONE/D(M),  ABS(X(M)))
         DEL = HALF * DEL
         IF ( ABS(DEL/DEL0) .LE. HLIM) GO TO 140
C
 30   HES = -IV(H)
C
C  ***  SET  G = (G - GSAVE)/DEL  ***
C
      DEL = ONE / DEL
      DO 40 I = 1, P
         G(I) = DEL * (G(I) - V(GSAVE1))
         GSAVE1 = GSAVE1 + 1
 40      CONTINUE
C
C  ***  ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX  ***
C
      K = HES + M*(M-1)/2
      L = K + M - 2
      IF (M .EQ. 1) GO TO 60
C
C  ***  SET  H(I,M) = 0.5 * (H(I,M) + G(I))  FOR I = 1 TO M-1  ***
C
      MM1 = M - 1
      DO 50 I = 1, MM1
         IF (B(1,I) .LT. B(2,I)) V(K) = HALF * (V(K) + G(I))
         K = K + 1
 50      CONTINUE
C
C  ***  ADD  H(I,M) = G(I)  FOR I = M TO P  ***
C
 60   L = L + 1
      DO 70 I = M, P
         IF (B(1,I) .LT. B(2,I)) V(L) = G(I)
         L = L + I
 70      CONTINUE
C
 80   M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 340
      IF (B(1,M) .GE. B(2,M)) GO TO 80
C
C  ***  CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE  ***
C
      DEL = V(DELTA0) * AMAX1(ONE/D(M),  ABS(X(M)))
      XM = X(M)
      IF (XM .LT. ZERO) GO TO 90
         XM1 = XM + DEL
         IF (XM1 .LE. B(2,M)) GO TO 110
           XM1 = XM - DEL
           IF (XM1 .GE. B(1,M)) GO TO 100
           GO TO 280
 90    XM1 = XM - DEL
       IF (XM1 .GE. B(1,M)) GO TO 100
       XM1 = XM + DEL
       IF (XM1 .LE. B(2,M)) GO TO 110
       GO TO 280
C
 100  DEL = -DEL
 110  V(XMSAVE) = XM
      X(M) = XM1
      V(DELTA) = DEL
      IRT = 2
      GO TO 999
C
C  ***  COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY.
C
 120  STP0 = IV(W) + P - 1
      MM1 = M - 1
      MM1O2 = M*MM1/2
      HES = -IV(H)
      IF (M .GT. 0) GO TO 130
C        ***  FIRST CALL ON  F7DHB.  ***
         IV(SAVEI) = 0
         GO TO 240
C
 130  IF (IV(TOOBIG) .EQ. 0) GO TO 150
C        ***  PUNT IN THE EVENT OF AN OVERSIZE STEP  ***
 140     IV(FDH) = -2
         GO TO 350
 150  I = IV(SAVEI)
      IF (I .GT. 0) GO TO 190
C
C  ***  SAVE F(X + STP(M)*E(M)) IN H(P,M)  ***
C
      PP1O2 = P * (P-1) / 2
      HPM = HES + PP1O2 + MM1
      V(HPM) = V(F)
C
C  ***  START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H.  ***
C
      NEWM1 = 1
      GO TO 260
 160  HMI = HES + MM1O2
      IF (MM1 .EQ. 0) GO TO 180
      HPI = HES + PP1O2
      DO 170 I = 1, MM1
         T = ZERO
         IF (B(1,I) .LT. B(2,I)) T = V(FX) - (V(F) + V(HPI))
         V(HMI) = T
         HMI = HMI + 1
         HPI = HPI + 1
 170     CONTINUE
 180  V(HMI) = V(F) - TWO*V(FX)
      IF (OFFSID) V(HMI) = V(FX) - TWO*V(F)
C
C  ***  COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H.  ***
C
      I = 0
      GO TO 200
C
 190  X(I) = V(DELTA)
C
C  ***  FINISH COMPUTING H(M,I)  ***
C
      STPI = STP0 + I
      HMI = HES + MM1O2 + I - 1
      STPM = STP0 + M
      V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM))
 200  I = I + 1
      IF (I .GT. M) GO TO 230
         IF (B(1,I) .LT. B(2,I)) GO TO 210
         GO TO 200
C
 210  IV(SAVEI) = I
      STPI = STP0 + I
      V(DELTA) = X(I)
      X(I) = X(I) + V(STPI)
      IRT = 1
      IF (I .LT. M) GO TO 999
      NEWM1 = 2
      GO TO 260
 220  X(M) = V(XMSAVE) - DEL
      IF (OFFSID) X(M) = V(XMSAVE) + TWO*DEL
      GO TO 999
C
 230  IV(SAVEI) = 0
      X(M) = V(XMSAVE)
C
 240  M = M + 1
      IV(MODE) = M
      IF (M .GT. P) GO TO 330
      IF (B(1,M) .LT. B(2,M)) GO TO 250
      GO TO 240
C
C  ***  PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H.
C  ***  COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN
C  ***  F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR.
C
 250  V(XMSAVE) = X(M)
      NEWM1 = 3
 260  XM = V(XMSAVE)
      DEL = V(DLTFDC) * AMAX1(ONE/D(M),  ABS(XM))
      XM1 = XM + DEL
      OFFSID = .FALSE.
      IF (XM1 .LE. B(2,M)) GO TO 270
         OFFSID = .TRUE.
         XM1 = XM - DEL
         IF (XM - TWO*DEL .GE. B(1,M)) GO TO 300
         GO TO 280
 270   IF (XM-DEL .GE. B(1,M)) GO TO 290
       OFFSID = .TRUE.
       IF (XM + TWO*DEL .LE. B(2,M)) GO TO 310
C
 280  IV(FDH) = -2
      GO TO 350
C
 290  IF (XM .GE. ZERO) GO TO 310
      XM1 = XM - DEL
 300  DEL = -DEL
 310  GO TO (160, 220, 320), NEWM1
 320  X(M) = XM1
      STPM = STP0 + M
      V(STPM) = DEL
      IRT = 1
      GO TO 999
C
C  ***  HANDLE SPECIAL CASE OF B(1,P) = B(2,P) -- CLEAR SCRATCH VALUES
C  ***  FROM LAST ROW OF FDH...
C
 330  IF (B(1,P) .LT. B(2,P)) GO TO 340
         I = HES + P*(P-1)/2
         CALL  V7SCP(P, V(I), ZERO)
C
C  ***  RESTORE V(F), ETC.  ***
C
 340  IV(FDH) = HES
 350  V(F) = V(FX)
      IRT = 3
      IF (KIND .LT. 0) GO TO 999
         IV(NFGCAL) = IV(SWITCH)
         GSAVE1 = IV(W) + P
         CALL V7CPY(P, G, V(GSAVE1))
         GO TO 999
C
 999  RETURN
C  ***  LAST LINE OF  F7DHB FOLLOWS  ***
      END

.