/* nag_opt_sparse_nlp_jacobian (e04vjc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL usrfun(Integer *status, Integer n, const double x[],
                              Integer needf, Integer nf, double f[],
                              Integer needg, Integer leng, double g[],
                              Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  double objadd, sinf;
  Integer exit_status = 0;
  Integer i, lena, leng, n, nea, neg, nf, nfname, ninf, ns, nxname, objrow;

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

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

  /* By default e04vhc does not print monitoring information.
     Define SHOW_MONITORING_INFO to turn it on - see further below. */
#ifdef SHOW_MONITORING_INFORMATION
  Nag_FileID fileid;
#endif

  exit_status = 0;
  INIT_FAIL(fail);

  printf("nag_opt_sparse_nlp_jacobian (e04vjc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\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)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("Invalid n or nf\n");
    exit_status = 1;
    goto END;
  }

  /* 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) {
    printf("Initialization of nag_opt_sparse_nlp_init (e04vgc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

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

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

  /* Illustrate how to pass information to the user-supplied
     function usrfun via the comm structure */
  comm.p = 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) {
    printf("nag_opt_sparse_nlp_jacobian (e04vjc) failed to determine the"
           " Jacobian structure\n");
    exit_status = 1;
    goto END;
  }

  /* Print the Jacobian structure. */

  printf("\n");
  printf("NEA (the number of nonzero entries in A) = %3" NAG_IFMT "\n", nea);

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

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

  printf("\n");
  printf("NEG (the number of nonzero entries in G) = %3" NAG_IFMT "\n", neg);
  printf("  I     IGFUN(I)   JGVAR(I)\n");
  printf("----    --------   --------\n");

  for (i = 1; i <= neg; ++i) {
    printf("%3" NAG_IFMT "%10" NAG_IFMT "%10" NAG_IFMT "\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 = Nag_Cold;
  strcpy(prob, "        ");
  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. */
  scanf("%" NAG_IFMT "%*[^\n] ", &objrow);

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

#ifdef SHOW_MONITORING_INFO
  /* 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) {
    exit_status = 2;
    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) {
    exit_status = 1;
    goto END;
  }
#endif

  /* 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) {
    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
   */
  fflush(stdout);
  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, (const char **) xnames,
                           flow, fupp, (const char **) fnames, x, xstate,
                           xmul, f, fstate, fmul, &ns, &ninf, &sinf, &state,
                           &comm, &fail);
  if (n > 0 && nf > 0) {
    for (i = 0; i < nxname; i++)
      NAG_FREE(xnames[i]);
    for (i = 0; i < nfname; i++)
      NAG_FREE(fnames[i]);
  }
  if (fail.code == NE_NOERROR || fail.code == NW_NOT_FEASIBLE) {
    printf("\n");
    printf("Final objective value = %11.1f\n", f[objrow - 1]);
    printf("Optimal X = ");

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

  if (fail.code != NE_NOERROR)
    exit_status = 2;

END:
  NAG_FREE(fnames);
  NAG_FREE(xnames);
  NAG_FREE(a);
  NAG_FREE(f);
  NAG_FREE(flow);
  NAG_FREE(fmul);
  NAG_FREE(fupp);
  NAG_FREE(x);
  NAG_FREE(xlow);
  NAG_FREE(xmul);
  NAG_FREE(xupp);
  NAG_FREE(fstate);
  NAG_FREE(iafun);
  NAG_FREE(igfun);
  NAG_FREE(javar);
  NAG_FREE(jgvar);
  NAG_FREE(xstate);

  return exit_status;
}

static void NAG_CALL 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]

  /* Check whether information came from the main program
     via the comm structure. Even if it was, we ignore it
     in this example. */
  if (comm->p)
    printf("Pointer %p was passed to usrfun via the comm struct\n", comm->p);

  /* 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 */