/* nag_opt_handle_solve_dfls (e04ffc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */

#include <stdio.h>
#include <nag.h>
#include <assert.h>

typedef struct pdata
{
  int    ny, nz;
  double *y, *z;
} pdata;

static void free_pdata(pdata pd);

#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL objfun(Integer nvar, const double x[],
                            Integer nres, double rx[],
                            Integer *inform, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  const int    defbounds = 1;
  const double infbnd = 1.0e20;

  pdata        pd;
  Integer      nvar, nres, isparse, nnzrd;
  double       x[4] = { 0.25, 0.39, 0.415, 0.39 };
  double       rinfo[100], stats[100];
  double       *rx, *lx, *ux;
  void         *handle;
  Integer      exit_status = 0;

  /* Nag Types */
  Nag_Comm     comm;
  NagError     fail;

  printf("nag_opt_handle_solve_dfls (e04ffc) Example Program Results\n\n");
  fflush(stdout);
  
  /* Fill the problem data structure */
  nvar = 4;
  nres = 11;
  pd.ny = nres;
  pd.nz = nres;
  pd.y = NAG_ALLOC(pd.ny,double); assert(pd.y);
  pd.z = NAG_ALLOC(pd.nz,double); assert(pd.z);
  pd.y[0]  = 4.0   ; pd.z[0] = 0.1957;
  pd.y[1]  = 2.0   ; pd.z[1] = 0.1947;
  pd.y[2]  = 1.0   ; pd.z[2] = 0.1735;
  pd.y[3]  = 0.5   ; pd.z[3] = 0.16;
  pd.y[4]  = 0.25  ; pd.z[4] = 0.0844;
  pd.y[5]  = 0.167 ; pd.z[5] = 0.0627;
  pd.y[6]  = 0.125 ; pd.z[6] = 0.0456;
  pd.y[7]  = 0.1   ; pd.z[7] = 0.0342;
  pd.y[8]  = 0.0833; pd.z[8] = 0.0323;
  pd.y[9]  = 0.0714; pd.z[9] = 0.0235;
  pd.y[10] = 0.0625; pd.z[10] = 0.0246;

  /* nag_opt_handle_init (e04rac).
   * Initialize the handle
   */
  nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);

  /* nag_opt_handle_set_nlnls (e04rmc)
   * Define residuals structure, isparse=0 means the residual structure is
   * dense => irowrd and icolrd arguments can be NULL
   */
  isparse = 0;
  nnzrd = 1;
  nag_opt_handle_set_nlnls(handle, nres, isparse, nnzrd, NULL, NULL,
                           NAGERR_DEFAULT);

  /* nag_opt_handle_opt_set (e04zmc)
   * Set options
   */
  /* Relax the main convergence criteria a bit */
  nag_opt_handle_opt_set(handle, "DFLS Trust Region Tolerance = 5.0e-6",
                         NAGERR_DEFAULT);
  /* Print the solution */
  nag_opt_handle_opt_set(handle, "Print Solution = X", NAGERR_DEFAULT);

  /* Optionally define bounds for the second and the fourth variable */
  if (defbounds)
    {
      lx = NAG_ALLOC(nvar, double); assert(lx);
      ux = NAG_ALLOC(nvar, double); assert(ux);
      lx[0] = -infbnd; ux[0] = infbnd;
      lx[1] = 0.2; ux[1] = 1.0;
      lx[2] = -infbnd; ux[2] = infbnd;
      lx[3] = 0.3; ux[3] = infbnd;
      /* nag_opt_handle_set_simplebounds (e04rhc) */
      nag_opt_handle_set_simplebounds(handle, nvar, lx, ux, NAGERR_DEFAULT);
    }

  /* nag_opt_handle_solve_dfls (e04ffc)
   * Call the solver
   */
  rx = NAG_ALLOC(nres, double); assert(rx);
  comm.p = &pd;
  INIT_FAIL(fail);
  nag_opt_handle_solve_dfls(handle, objfun, NULL, nvar, x, nres, rx,
                                 rinfo, stats, &comm, &fail);
  if (fail.code != NE_NOERROR){
    printf("Error from nag_opt_handle_solve_dfls (e04ffc).\n%s\n",
           fail.message);
    exit_status = 1;
  }

  /* Clean data */
  if (handle)
    /* nag_opt_handle_free (e04rzc).
     * Destroy the problem handle and deallocate all the memory used 
     */
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);
  free_pdata(pd);
  if (rx)
    NAG_FREE(rx);
  if (lx)
    NAG_FREE(lx);
  if (ux)
    NAG_FREE(ux);

  return exit_status;
}

static void NAG_CALL objfun(Integer nvar, const double x[],
                            Integer nres, double rx[],
                            Integer *inform, Nag_Comm *comm)
{
  pdata  *pd;
  int    i;
  double r1, r2;

  /* Interrupt the solver if the comm structure is not correctly initialized */
  if (!comm || !(comm->p))
    {
      *inform = -1;
      return;
    }

  /* Extract the problem data from the comm structure */
  pd = (pdata *) comm->p;

  /* Interrupt the solver if the data does not correspond to the problem */
  if (nvar != 4 || nres != 11 || pd->ny != nres || pd->nz != nres)
    {
      *inform = -1;
      return;
    }

  /* Fill the residuals values */
  for (i = 0; i < nres; i++)
    {
      r1 = pd->y[i] * (pd->y[i]+x[1]);
      r2 = pd->y[i] * (pd->y[i]+x[2]) + x[3];
      rx[i] = pd->z[i] - x[0]*r1/r2;
    }
}

static void free_pdata(pdata pd)
{
  if (pd.y)
    NAG_FREE(pd.y);
  if (pd.z)
    NAG_FREE(pd.z);
}