/* nag_opt_bounds_2nd_deriv (e04lbc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

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

#define NMAX 4

#ifdef __cplusplus
extern "C" {
#endif
  static void funct(Integer n,  double xc[], double *fc, double gc[],
                    Nag_Comm *comm);
  static void hess(Integer n,  double xc[], double fhesl[],
                   double fhesd[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
static int ex1(void);
static int ex2(void);
static void funct(Integer n, double xc[], double *fc, double gc[],
                  Nag_Comm *comm)
{
  /* Function to evaluate objective function and its 1st derivatives. */
  double term1, term1_sq;
  double term2, term2_sq;
  double term3, term3_sq, term3_cu;
  double term4, term4_sq, term4_cu;

  term1 = xc[0] + 10.0*xc[1];
  term1_sq = term1*term1;

  term2 = xc[2] - xc[3];
  term2_sq = term2*term2;

  term3 = xc[1] - 2.0*xc[2];
  term3_sq = term3*term3;
  term3_cu = term3*term3_sq;

  term4 = xc[0] - xc[3];
  term4_sq = term4*term4;
  term4_cu = term4_sq*term4;

  *fc = term1_sq + 5.0*term2_sq
    + term3_sq*term3_sq + 10.0*term4_sq*term4_sq;

  gc[0] =   2.0*term1 + 40.0*term4_cu;
  gc[1] =  20.0*term1 +  4.0*term3_cu;
  gc[2] =  10.0*term2 -  8.0*term3_cu;
  gc[3] = -10.0*term2 - 40.0*term4_cu;
}
/* funct */
static void hess(Integer n, double xc[], double fhesl[],
                 double fhesd[], Nag_Comm *comm)
{
  /* Routine to evaluate 2nd derivatives */
  double term3_sq;
  double term4_sq;

  term3_sq = (xc[1] - 2.0*xc[2])*(xc[1] - 2.0*xc[2]);
  term4_sq = (xc[0] - xc[3])*(xc[0] - xc[3]);

  fhesd[0] = 2.0 + 120.0*term4_sq;
  fhesd[1] = 200.0 + 12.0*term3_sq;
  fhesd[2] = 10.0 + 48.0*term3_sq;
  fhesd[3] = 10.0 + 120.0*term4_sq;

  fhesl[0] = 20.0;
  fhesl[1] = 0.0;
  fhesl[2] = -24.0*term3_sq;
  fhesl[3] = -120.0*term4_sq;
  fhesl[4] = 0.0;
  fhesl[5] = -10.0;
}
/* hess */

int main(void)
{

  Integer exit_status_ex1=0;
  Integer exit_status_ex2=0;


  /* Two examples are called, ex1() which uses the
   * default settings to solve the problem and
   * ex2() which solves the same problem with
   * some optional parameters set by the user.
   */

  Vprintf("nag_opt_bounds_2nd_deriv (e04lbc) Example Program Results.\n");

  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();
  return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1;
}

static int ex1(void)
{
  Integer exit_status=0;
  Integer n = NMAX;
  double *x=0, *bl=0, *bu=0, *g=0;
  double f;
  NagError fail;

  INIT_FAIL(fail);


  /* Function Body */

  Vprintf("\nnag_opt_bounds_2nd_deriv (e04lbc) example 1: no option"
          " setting.\n");

  if (n>=1)
    {
      if ( !( x = NAG_ALLOC(n, double)) ||
           !( bl = NAG_ALLOC(n, double)) ||
           !( bu = NAG_ALLOC(n, double)) ||
           !( g = 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[0] = 1.46;
  x[1] = -0.82;
  x[2] = 0.57;
  x[3] = 1.21;

  bl[0] = 1.0;
  bu[0] = 3.0;
  bl[1] = -2.0;
  bu[1] = 0.0;

  /* x[2] is not bounded, so we set bl[2] to a large negative
   *  number and bu[2] to a large positive number
   */
  bl[2] = -1e6;
  bu[2] = 1e6;
  bl[3] = 1.0;
  bu[3] = 3.0;

  /* Set up starting point */
  x[0] = 3.0;
  x[1] = -1.0;
  x[2] = 0.0;
  x[3] = 1.0;

  /* nag_opt_bounds_2nd_deriv (e04lbc).
   * Solves bound constrained problems (first and second
   * derivatives required)
   */
  nag_opt_bounds_2nd_deriv(n, funct, hess, Nag_Bounds, bl, bu, x, &f, g,
                           E04_DEFAULT, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error/Warning from nag_opt_bounds_2nd_deriv (e04lbc).\n%s\n",
              fail.message);
      if (fail.code != NW_COND_MIN)
        exit_status = 1;
      goto END;
    }
 END:
  if (x) NAG_FREE(x);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (g) NAG_FREE(g);
  return exit_status;
}
/* ex1 */

static int ex2(void)
{

  Nag_Boolean print;
  Integer exit_status=0, n=NMAX;
  NagError fail;
  Nag_Comm comm;
  Nag_E04_Opt options;
  double *bl=0, *bu=0, f, *g=0, *x=0;

  INIT_FAIL(fail);

  /* Function Body */
  Vprintf("\n\nnag_opt_bounds_2nd_deriv (e04lbc) example 2: using option"
          " setting.\n");

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

  bl[0] = 1.0;
  bu[0] = 3.0;
  bl[1] = -2.0;
  bu[1] = 0.0;

  /* x[2] is not bounded, so we set bl[2] to a large negative
   *  number and bu[2] to a large positive number
   */
  bl[2] = -1e6;
  bu[2] = 1e6;
  bl[3] = 1.0;
  bu[3] = 3.0;

  /* Set up starting point */
  x[0] = 3.0;
  x[1] = -1.0;
  x[2] = 0.0;
  x[3] = 1.0;

  print = Nag_TRUE;
  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options);
  /* nag_opt_read (e04xyc).
   * Read options from a text file
   */
  nag_opt_read("e04lbc", "stdin", &options, print, "stdout", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from e0xyc.\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  /* nag_opt_bounds_2nd_deriv (e04lbc), see above. */
  nag_opt_bounds_2nd_deriv(n, funct, hess, Nag_Bounds, bl, bu, x, &f, g,
                           &options, &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from e04lbc/Warning.\n%s\n", fail.message);
      if (fail.code != NW_COND_MIN)
        exit_status = 1;
    }

  /* Free memory allocated by nag_opt_bounds_deriv (e04kbc) to pointers hesd,
   * hesl and state.
   */
  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "all", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_bounds_2nd_deriv (e04lbc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (x) NAG_FREE(x);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (g) NAG_FREE(g);
  return exit_status;
} /* ex2 */