/* nag_opt_check_2nd_deriv(e04hdc) Example Program. 
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 */
#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nage04.h>

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);

main()
{
  double hesd[4];
  double hesl[6],  f;
  double g[4];
  double x[4];
  
  Integer n;
  Integer i, j, k;
  
  Nag_Comm comm;
  
#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("e04hdc Example Program Results\n\n");
  
  /*  Set up an arbitrary point at which to check the derivatives */
  n = 4;
  X(1) = 1.46;
  X(2) = -.82;
  X(3) = .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 */
  e04hcc(n,  funct,  &X(1),  &f,  &G(1),  &comm, NAGERR_DEFAULT);
  
  /*  Check the 2nd derivatives */
  e04hdc(n,  funct,  hess,  &X(1),  &G(1),  &HESL(1),  &HESD(1),  
	 &comm, NAGERR_DEFAULT);
  
  Vprintf("\n2nd derivatives are consistent with 1st derivatives.\n\n");
  Vprintf("%s%12.4e\n",
	  "At the test point, funct gives the function value, ", 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;
    }
  exit(EXIT_SUCCESS);
}

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;
}