/* nag_opt_estimate_deriv(e04xac) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 * Mark 7 revised, 2001.
 *
 */

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

static void ex1(void);
static void ex2(void);
#ifdef __cplusplus
extern "C" {
#endif
static void objfun(Integer n, double x[], double *objf,
                   double g[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

static void objfun(Integer n, double x[], double *objf,
                   double g[], Nag_Comm *comm)
{
  double a, asq, b, bsq, c, csq, d, dsq;

  a = x[0] + 10.0*x[1];
  b = x[2] - x[3];
  c = x[1] - 2.0*x[2];
  d = x[0] - x[3];
  asq = a*a;
  bsq = b*b;
  csq = c*c;
  dsq = d*d;
  *objf = asq + 5.0*bsq + csq*csq + 10.0*dsq*dsq;
  if (comm->flag == 2)
    {
      g[0] =   2.0*a + 40.0*d*dsq;
      g[1] =  20.0*a + 4.0*c*csq;
      g[2] =  10.0*b -  8.0*c*csq;
      g[3] = -10.0*b - 40.0*d*dsq;

    }
}
/* objfun */

int main(void)
{
  Vprintf("e04xac Example Program Results\n");
  ex1();
  ex2();
  return EXIT_SUCCESS;
}

static void ex1(void)
{
#define MAXN 4

  /* Local variables */
  Integer n, tdh;

  double objf;
  double x[MAXN];
  double g[MAXN], h[MAXN][MAXN];

  Nag_DerivInfo deriv_info[MAXN];
  static NagError fail;

  n = MAXN;
  tdh = MAXN;
  Vprintf("\nExample 1: default options\n");

  x[0] = 3.0;
  x[1] = -1.0;
  x[2] = 0.0;
  x[3] = 1.0;
  fail.print = TRUE;

  /* Pass null pointers for the h_central and h_forward parameters
   * as we do not need these values.
   */
  e04xac(n, x, objfun, &objf, g, (double*)0, (double*)0,
         &h[0][0], tdh, deriv_info, E04_DEFAULT, NAGCOMM_NULL, &fail);

} /* ex1 */

static void ex2(void)
{
  /* Local variables */
  Integer i, j;
  Integer n, tdh;

  double objf;
  double x[MAXN];
  double h_central[MAXN];
  double h_forward[MAXN];
  double g[MAXN], h[MAXN][MAXN], hess_diag[MAXN];


  Nag_DerivInfo deriv_info[MAXN];
  Nag_E04_Opt options;
  static NagError fail;

  n = MAXN;
  tdh = MAXN;
  x[0] = 3.0;
  x[1] = -1.0;
  x[2] = 0.0;
  x[3] = 1.0;
  fail.print = TRUE;

  Vprintf("\nExample 2: some options are set\n");

  e04xxc(&options);
  options.list = FALSE;
  options.print_deriv = Nag_D_NoPrint;

  options.deriv_want = Nag_Grad_HessDiag;

  Vprintf("\nEstimate gradient and Hessian diagonals given function only\n");

  /* Note: it is acceptable to pass an array of length n (hess_diag)
   * as the Hessian parameter in this case.
   */
  e04xac(n, x, objfun, &objf, g, h_forward, h_central,
         hess_diag, tdh, deriv_info, &options, NAGCOMM_NULL, &fail);

  Vprintf("\nFunction value: %12.4e\n", objf);
  Vprintf("Estimated gradient vector\n");
  for (i = 0; i < n; ++i)
    Vprintf("%12.4e  ", g[i]);
  Vprintf("\nEstimated Hessian matrix diagonal\n");
  for (i = 0; i < n; ++i)
    Vprintf("%12.4e  ", hess_diag[i]);
  Vprintf("\n");

  options.deriv_want = Nag_HessFull;

  Vprintf("\nEstimate full Hessian given function and gradients\n");
  e04xac(n, x, objfun, &objf, g, h_forward, h_central,
         &h[0][0], tdh, deriv_info, &options, NAGCOMM_NULL, &fail);

  Vprintf("\nFunction value: %12.4e\n", objf);
  Vprintf("Computed gradient vector\n");
  for (i = 0; i < n; ++i)
    Vprintf("%12.4e  ", g[i]);
  Vprintf("\nEstimated Hessian matrix\n");
  for (i = 0; i < n; ++i)
    {
      for (j = 0; j < n; ++j)
        Vprintf("%12.4e  ", h[i][j]);
      Vprintf("\n");
    }
  exit(EXIT_SUCCESS);
} /* ex2 */