/* nag_nag_opt_sparse_nlp_jacobian (e04vjc) Example Program.
 *
 * Copyright 2004 Numerical Algorithms Group.
 *
 * Mark 8, 2004.
 */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>
#include <nagx04.h>

static void usrfun(Integer *status, Integer n, const double x[],
                   Integer needf, Integer nf, double f[], Integer needg,
                   Integer leng, double g[], Nag_Comm *comm);
int main(void)
{

  /* Scalars */
  double objadd, sinf;
  Integer exit_status, i, lena, leng, n, nea, neg, nf, nfname, ninf, ns, nxname;
  Integer objrow, start_int;

  /* Arrays */
  char **fnames=0, prob[9], **xnames=0;
  double *a=0, *f=0, *flow=0, *fmul=0, *fupp=0, *x=0, *xlow=0, *xmul=0, *xupp=0;
  Integer *fstate=0, *iafun=0, *igfun=0, *javar=0, *jgvar=0, *xstate=0;

  /* Nag Types*/
  Nag_E04State state;
  NagError fail;
  Nag_Comm comm;
  Nag_Start start;
  Nag_FileID fileid;

  exit_status = 0;
  INIT_FAIL(fail);
  Vprintf("nag_opt_sparse_nlp_jacobian (e04vjc) Example Program Results\n");

  /* Skip heading in data file */
  Vscanf("%*[^\n] ");

  Vscanf("%ld%ld%*[^\n] ",   &n,  &nf);
  if (n > 0 && nf > 0)
    {
      nfname = 1;
      nxname = 1;
      lena = 300;
      leng = 300;
      /* Allocate memory */
      if ( !(fnames = NAG_ALLOC(nfname, char *)) ||
           !(xnames = NAG_ALLOC(nxname, char *)) ||
           !(a = NAG_ALLOC(lena, double)) ||
           !(f = NAG_ALLOC(nf, double)) ||
           !(flow = NAG_ALLOC(nf, double)) ||
           !(fmul = NAG_ALLOC(nf, double)) ||
           !(fupp = NAG_ALLOC(nf, double)) ||
           !(x = NAG_ALLOC(n, double)) ||
           !(xlow = NAG_ALLOC(n, double)) ||
           !(xmul = NAG_ALLOC(n, double)) ||
           !(xupp = NAG_ALLOC(n, double)) ||
           !(fstate = NAG_ALLOC(nf, Integer)) ||
           !(iafun = NAG_ALLOC(lena, Integer)) ||
           !(igfun = NAG_ALLOC(leng, Integer)) ||
           !(javar = NAG_ALLOC(lena, Integer)) ||
           !(jgvar = NAG_ALLOC(leng, Integer)) ||
           !(xstate = NAG_ALLOC(n, Integer)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      Vprintf("Invalid n or nf\n");
      exit_status = 1;
      goto END;
    }

  /* Call nag_opt_sparse_nlp_init (e04vgc) to initialise e04vjc. */
  /* nag_opt_sparse_nlp_init (e04vgc).
   * Initialization function for nag_opt_sparse_nlp_solve
   * (e04vhc)
   */
  nag_opt_sparse_nlp_init(&state, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Initialisation of nag_opt_sparse_nlp_init (e04vgc) failed.\n");
      exit_status = 1;
      goto END;
    }

  /* Read the bounds on the variables. */
  for (i = 1; i <= n; ++i)
    {
      Vscanf("%lf%lf%*[^\n] ",   &xlow[i - 1],  &xupp[i - 1]);
    }

  for (i = 1; i <= n; ++i)
    {
      x[i - 1] = 0.;
    }

  /* Determine the Jacobian structure. */
  /* nag_opt_sparse_nlp_jacobian (e04vjc).
   * Determine the pattern of nonzeros in the Jacobian matrix
   * for nag_opt_sparse_nlp_solve (e04vhc)
   */
  nag_opt_sparse_nlp_jacobian(nf, n, usrfun, iafun, javar, a, lena, &nea, igfun,
                              jgvar, leng, &neg, x, xlow, xupp, &state, &comm,
                              &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("nag_opt_sparse_nlp_jacobian (e04vjc) failed to determine the"
              " Jacobian structure\n");
      exit_status = 1;
      goto END;
    }

  /* Print the Jacobian structure. */

  Vprintf("\n");
  Vprintf("NEA (the number of non-zero entries in A) = %3ld\n",   nea);

  Vprintf("  I     IAFUN(I)   JAVAR(I)          A(I)\n");
  Vprintf("----    --------   --------   -----------\n");

  for (i = 1; i <= nea; ++i)
    {
      Vprintf("%3ld%10ld%10ld%18.4e\n", i, iafun[i - 1],
              javar[i - 1], a[i - 1]);
    }

  Vprintf("\n");
  Vprintf("NEG (the number of non-zero entries in G) = %3ld\n",   neg);
  Vprintf("  I     IGFUN(I)   JGVAR(I)\n");
  Vprintf("----    --------   --------\n");

  for (i = 1; i <= neg; ++i)
    {
      Vprintf("%3ld%10ld%10ld\n", i, igfun[i - 1],
              jgvar[i - 1]);
    }

  /* Now that we have the determined the structure of the
   * Jacobian, set up the information necessary to solve
   * the optimization problem.
   */
  start_int = 0;
  if (start_int == 0)
    {
      start = Nag_Cold;
    }
  else
    {
      start = Nag_Warm;
    }
  strcpy(prob, "e04vjce");
  objadd = 0.0;
  for (i = 1; i <= n; ++i)
    {
      x[i - 1] = 0.;
      xstate[i - 1] = 0;
      xmul[i - 1] = 0.;
    }
  for (i = 1; i <= nf; ++i)
    {
      f[i - 1] = 0.;
      fstate[i - 1] = 0;
      fmul[i - 1] = 0.;
    }

  /* The row containing the objective function. */
  Vscanf("%ld%*[^\n] ",   &objrow);

  /* Read the bounds on the functions. */
  for (i = 1; i <= nf; ++i)
    {
      Vscanf("%lf%lf%*[^\n] ",   &flow[i - 1],  &fupp[i - 1]);
    }

  /* By default nag_opt_sparse_nlp_solve (e04vhc) does not print monitoring
   * information. Call nag_open_file (x04acc) to set the print file fileid
   */
  /* nag_open_file (x04acc).
   * Open unit number for reading, writing or appending, and
   * associate unit with named file
   */
  nag_open_file("", 2, &fileid, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Fileid could not be obtained.\n");
      exit_status = 1;
      goto END;
    }
  /* nag_opt_sparse_nlp_option_set_integer (e04vmc).
   * Set a single option for nag_opt_sparse_nlp_solve (e04vhc)
   * from an integer argument
   */
  nag_opt_sparse_nlp_option_set_integer("Print file", fileid, &state, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("nag_opt_sparse_nlp_option_set_integer (e04vmc) failed to set"
              " Print File\n");
      exit_status = 1;
      goto END;
    }

  /* Tell nag_opt_sparse_nlp_solve (e04vhc) that we supply no derivatives in
   *  usrfun. */
  /* nag_opt_sparse_nlp_option_set_string (e04vlc).
   * Set a single option for nag_opt_sparse_nlp_solve (e04vhc)
   * from a character string
   */
  nag_opt_sparse_nlp_option_set_string("Derivative option 0", &state, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("nag_opt_sparse_nlp_option_set_string (e04vlc) failed to set"
              " Derivative option\n");
      exit_status = 1;
      goto END;
    }
  for (i = 1; i <= nfname; ++i)
    {
      fnames[i -1] = NAG_ALLOC(9, char);
      strcpy(fnames[i-1],"");
    }

  for (i = 1; i <= nxname; ++i)
    {
      xnames[i-1] = NAG_ALLOC(9, char);
      strcpy(xnames[i-1],"");
    }

  /* Solve the problem. */
  /* nag_opt_sparse_nlp_solve (e04vhc).
   * General sparse nonlinear optimizer
   */
  nag_opt_sparse_nlp_solve(start, nf, n, nxname, nfname, objadd, objrow, prob,
                           usrfun, iafun, javar, a, lena, nea, igfun, jgvar,
                           leng, neg, xlow, xupp, xnames, flow, fupp, fnames, x,
                           xstate, xmul, f, fstate, fmul, &ns, &ninf, &sinf,
                           &state, &comm, &fail);

  if (fail.code == NE_NOERROR || fail.code == NW_NOT_FEASIBLE)
    {
      Vprintf("Final objective value = %11.1f\n",   f[objrow - 1]);
      Vprintf("Optimal X = ");

      for (i = 1; i <= n; ++i)
        {
          Vprintf("%9.2f%s", x[i - 1], i%7 == 0 || i == n ?"\n":" ");
        }
    }
  else
    {
      Vprintf ("Error message from nag_opt_sparse_nlp_solve (e04vhc) %s\n",
               fail.message);
    }

 END:
  for (i=0; i < nxname; i++)
    {
      NAG_FREE(xnames[i]);
    }
  for (i=0; i < nfname; i++)
    {
      NAG_FREE(fnames[i]);
    }
  if (fnames) NAG_FREE(fnames);
  if (xnames) NAG_FREE(xnames);
  if (a) NAG_FREE(a);
  if (f) NAG_FREE(f);
  if (flow) NAG_FREE(flow);
  if (fmul) NAG_FREE(fmul);
  if (fupp) NAG_FREE(fupp);
  if (x) NAG_FREE(x);
  if (xlow) NAG_FREE(xlow);
  if (xmul) NAG_FREE(xmul);
  if (xupp) NAG_FREE(xupp);
  if (fstate) NAG_FREE(fstate);
  if (iafun) NAG_FREE(iafun);
  if (igfun) NAG_FREE(igfun);
  if (javar) NAG_FREE(javar);
  if (jgvar) NAG_FREE(jgvar);
  if (xstate) NAG_FREE(xstate);



  return exit_status;
}

static void usrfun(Integer *status, Integer n, const double x[],
                   Integer needf, Integer nf, double f[], Integer needg,
                   Integer leng, double g[], Nag_Comm *comm)
{
  /* Parameter adjustments */
#define X(I) x[(I)-1]
#define F(I) f[(I)-1]
#define G(I) g[(I)-1]

  /* Function Body */
  if (needf > 0)
    {
      F(1) = sin(-X(1) - .25) * 1e3 + sin(-X(2) - .25) * 1e3 - X(3);
      F(2) = sin(X(1) - .25) * 1e3 + sin(X(1) - X(2) - .25) * 1e3 - X(4);
      F(3) = sin(X(2) - X(1) - .25) * 1e3 + sin(X(2) - .25) * 1e3;
      F(4) = -X(1) + X(2);
      F(5) = X(1) - X(2);
      F(6) = X(3) * (X(3) * X(3)) * 1e-6 + X(4) * (X(4) * X(4)) * 2e-6 / 3.
        + X(3) * 3 + X(4) * 2;
    }

  return;
} /* usrfun */