' **********************************************************************
' Numerical Hessian and Gradient
' **********************************************************************

SUB NumHessGrad (Func AS FUNCTION(X() AS DOUBLE) AS DOUBLE,_
                 X() AS DOUBLE, Eta AS DOUBLE,_ 
                 G() AS DOUBLE, H() AS DOUBLE)
' ----------------------------------------------------------------------
' Compute the hessian H and the gradient G of function Func at point X
' using finite difference approximations.
' Eta is the relative increment used to compute derivatives.
' ----------------------------------------------------------------------

  CONST EtaMin = 6.06D-6  ' Min. relative increment (~ MachEp^(1/3))

  DIM AS INTEGER L, N, I, J
  DIM AS DOUBLE  Temp1, Temp2, F, F2plus

  L = LBOUND(X): N = UBOUND(X)

  DIM AS DOUBLE Delta(L TO N)   ' Increment
  DIM AS DOUBLE Xminus(L TO N)  ' X - Delta
  DIM AS DOUBLE Xplus(L TO N)   ' X + Delta
  DIM AS DOUBLE Fminus(L TO N)  ' F(X - Delta)
  DIM AS DOUBLE Fplus(L TO N)   ' F(X + Delta)

  IF Eta < EtaMin THEN Eta = EtaMin

  F = Func(X())
  
  FOR I = L TO N
    IF X(I) <> 0 THEN 
      Delta(I) = Eta * ABS(X(I)) 
    ELSE 
      Delta(I) = Eta
    END IF  
    Xplus(I) = X(I) + Delta(I)
    Xminus(I) = X(I) - Delta(I)
  NEXT I

  FOR I = L TO N
    Temp1 = X(I)
    X(I) = Xminus(I)
    Fminus(I) = Func(X())
    X(I) = Xplus(I)
    Fplus(I) = Func(X())
    X(I) = Temp1
  NEXT I

  FOR I = L TO N
    G(I) = (Fplus(I) - Fminus(I)) / (2 * Delta(I))
    H(I, I) = (Fplus(I) + Fminus(I) - 2 * F) / (Delta(I) * Delta(I))
  NEXT I

  FOR I = L TO N - 1
    Temp1 = X(I)
    X(I) = Xplus(I)
    FOR J = I + L TO N
      Temp2 = X(J)
      X(J) = Xplus(J)
      F2plus = Func(X())
      H(I, J) = (F2plus - Fplus(I) - Fplus(J) + F) / (Delta(I) * Delta(J))
      H(J, I) = H(I, J)
      X(J) = Temp2
    NEXT J
    X(I) = Temp1
  NEXT I
END SUB

