/* nag_opt_check_2nd_deriv (e04hdc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
  static void hess(Integer n, double xc[], double fhesl[],
                   double fhesd[], Nag_Comm *comm);
  static void funct(Integer n, double xc[], double *fc,
                    double gc[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{

#define NMAX 4

  Integer exit_status=0, i, j, k, n;
  NagError fail;
  Nag_Comm comm;
  double f, *g=0, *hesd=0, *hesl=0, *x=0;

  INIT_FAIL(fail);

#define X(I) x[(I)-1]
#define HESL(I) hesl[(I)-1]
#define HESD(I) hesd[(I)-1]
#define G(I) g[(I)-1]

  Vprintf("nag_opt_check_2nd_deriv (e04hdc) Example Program Results\n\n");

  /* Set up an arbitrary point at which to check the derivatives */
  n = NMAX;

  if (n>=1)
    {
      if ( !( hesd = NAG_ALLOC(n, double)) ||
           !( hesl = NAG_ALLOC(n*(n-1)/2, double)) ||
           !( g = NAG_ALLOC(n, double)) ||
           !( x = NAG_ALLOC(n, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      Vprintf("Invalid n.\n");
      exit_status = 1;
      return exit_status;
    }

  X(1) = 1.46;
  X(2) = -0.82;
  X(3) = 0.57;
  X(4) = 1.21;

  Vprintf("The test point is\n");
  for (j = 1; j <= n; ++j)
    Vprintf("%9.4f", X(j));
  Vprintf("\n");

  /* Check the 1st derivatives */
  /* nag_opt_check_deriv (e04hcc).
   * Derivative checker for use with nag_opt_bounds_deriv
   * (e04kbc)
   */
  nag_opt_check_deriv(n, funct, &X(1), &f, &G(1), &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_check_deriv (e04hcc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Check the 2nd derivatives */
  /* nag_opt_check_2nd_deriv (e04hdc).
   * Checks second derivatives of a user-defined function
   */
  nag_opt_check_2nd_deriv(n, funct, hess, &X(1), &G(1), &HESL(1), &HESD(1),
                          &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_check_2nd_deriv (e04hdc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  Vprintf("\n2nd derivatives are consistent with 1st derivatives.\n\n");
  Vprintf("At the test point, funct gives the function value, %12.4e\n",
          f);
  Vprintf("and the 1st derivatives\n");
  for (j = 1; j <= n; ++j)
    Vprintf("%12.3e%s", G(j), j%4?"":"\n");

  Vprintf("\nhess gives the lower triangle of the Hessian matrix\n");
  Vprintf("%12.3e\n", HESD(1));
  k = 1;
  for (i = 2; i <= n; ++i)
    {
      for (j = k; j <= k + i - 2; ++j)
        Vprintf("%12.3e", HESL(j));
      Vprintf("%12.3e\n", HESD(i));
      k = k + i - 1;
    }
 END:
  if (hesd) NAG_FREE(hesd);
  if (hesl) NAG_FREE(hesl);
  if (g) NAG_FREE(g);
  if (x) NAG_FREE(x);
  return exit_status;
}
static void funct(Integer n, double xc[], double *fc,
                  double gc[], Nag_Comm *comm)
{

  /* Routine to evaluate objective function and its 1st derivatives. */

#define GC(I) gc[(I)-1]
#define XC(I) xc[(I)-1]


  *fc = pow(XC(1)+10.0*XC(2), 2.0)
    + 5.0*pow(XC(3)-XC(4), 2.0)
    + pow(XC(2)-2.0*XC(3), 4.0)
    + 10.0*pow(XC(1)-XC(4), 4.0);

  GC(1) = 2.0*(XC(1)+10.0*XC(2)) +
    40.0*pow(XC(1)-XC(4),3.0);
  GC(2) = 20.0*(XC(1)+10.0*XC(2)) +
    4.0*pow(XC(2)-2.0*XC(3),3.0);
  GC(3) = 10.0*(XC(3)-XC(4)) -
    8.0*pow(XC(2)-2.0*XC(3),3.0);
  GC(4) = 10.0*(XC(4)-XC(3)) -
    40.0*pow(XC(1)-XC(4), 3.0);
}
static void hess(Integer n, double xc[], double fhesl[],
                 double fhesd[], Nag_Comm *comm)
{
  /* Routine to evaluate 2nd derivatives */

#define FHESD(I) fhesd[(I)-1]
#define FHESL(I) fhesl[(I)-1]
#define XC(I) xc[(I)-1]

  FHESD(1) = 2.0 + 120.0*pow(XC(1)-XC(4), 2.0);
  FHESD(2) = 200.0 + 12.0*pow(XC(2)-2.0*XC(3), 2.0);
  FHESD(3) = 10.0 + 48.0*pow(XC(2)-2.0*XC(3), 2.0);
  FHESD(4) = 10.0 + 120.0*pow(XC(1)-XC(4), 2.0);
  FHESL(1) = 20.0;
  FHESL(2) = 0.0;
  FHESL(3) = -24.0*pow(XC(2)-2.0*XC(3), 2.0);
  FHESL(4) = -120.0*pow(XC(1)-XC(4), 2.0);
  FHESL(5) = 0.0;
  FHESL(6) = -10.0;
}