/* nag_opt_nlp_sparse (e04ugc) Example Program.
 *
 * Copyright 2000 Numerical Algorithms Group.
 *
 * NAG C Library
 *
 * Mark 6, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL confun(Integer ncnln, Integer njnln,
                            Integer nnzjac, const double x[], double conf[],
                            double conjac[], Nag_Comm *comm);

static void NAG_CALL objfun(Integer nonln,
                            const double x[], double *objf,
                            double objgrad[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

#define NAMES(I, J) names[(I)*9+J]
#define MAXNAMES 300

int main(int argc, char *argv[])
{
  FILE        *fpin, *fpout;
  char        *optionsfile = 0;
  char        *outfile = 0;
  Integer     exit_status = 0, *ha = 0, i, icol, iobj, j, jcol, *ka = 0, m, n,
              ncnln;
  Integer     ninf, njnln, nnz, nonln;
  Nag_E04_Opt options;
  char        **crnames = 0, *names = 0;
  double      *a = 0, *bl = 0, *bu = 0, obj, sinf, *xs = 0;
  NagError    fail;

  INIT_FAIL(fail);


  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  (void) nag_example_file_io(argc, argv, "-options", &optionsfile);
  (void) nag_example_file_io(argc, argv, "-nag_write", &outfile);
  if (!outfile)
    {
      outfile = NAG_ALLOC(7, char);
      strcpy(outfile, "stdout");
    }

  fprintf(fpout, "nag_opt_nlp_sparse (e04ugc) Example Program Results\n");

  /* Skip heading in data file*/
  fscanf(fpin, " %*[^\n]");
  /* Read the problem dimensions */
  fscanf(fpin, " %*[^\n]");
  fscanf(fpin, "%ld%ld", &n, &m);
  /* Read NCNLN, NONLN and NJNLN from data file. */
  fscanf(fpin, " %*[^\n]");
  fscanf(fpin, "%ld%ld%ld", &ncnln, &nonln, &njnln);
  /* Read NNZ, IOBJ */
  fscanf(fpin, " %*[^\n]");
  fscanf(fpin, "%ld%ld", &nnz, &iobj);

  if (!(a = NAG_ALLOC(nnz, double)) ||
      !(bl = NAG_ALLOC(n+m, double)) ||
      !(bu = NAG_ALLOC(n+m, double)) ||
      !(xs = NAG_ALLOC(n+m, double)) ||
      !(ha = NAG_ALLOC(nnz, Integer)) ||
      !(ka = NAG_ALLOC(n+1, Integer)) ||
      !(crnames = NAG_ALLOC(n+m, char *)) ||
      !(names = NAG_ALLOC((n+m)*9, char))
      )
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = 1;
      goto END;
    }

  /* Read the column and row names */
  fscanf(fpin, " %*[^\n]");
  fscanf(fpin, " %*[^']");
  for (i = 0; i < n+m; ++i)
    {
      fscanf(fpin, " '%8c'", &NAMES(i, 0));
      NAMES(i, 8) = '\0';
      crnames[i] = &NAMES(i, 0);
    }

  /* read the matrix and set up ka. */
  jcol = 1;
  ka[jcol - 1] = 0;
  fscanf(fpin, " %*[^\n]");
  for (i = 0; i < nnz; ++i)
    {
      /* a[i] stores (ha[i], icol) element of matrix */
      fscanf(fpin, "%lf%ld%ld", &a[i], &ha[i], &icol);
      if (icol < jcol)
        {
          /* Elements not ordered by increasing column index. */
          fprintf(fpout, "Element in column%5ld found after element in"
                  " column%5ld. Problem abandoned.\n", icol, jcol);
          exit_status = 1;
          goto END;
        }
      else if (icol == jcol + 1)
        {
          /* Index in a of the start of the icol-th column equals i. */
          ka[icol - 1] = i;
          jcol = icol;
        }
      else if (icol > jcol + 1)
        {
          /* Index in a of the start of the icol-th column equals i,
           *  but columns jcol+1,jcol+2,...,icol-1 are empty. Set the
           *  corresponding elements of ka to i.
           */
          for (j = jcol + 1; j <= icol - 1; ++j)
            ka[j - 1] = i;

          ka[icol - 1] = i;
          jcol = icol;
        }
    }
  ka[n] = nnz;
  if (n > icol)
    {
      /* Columns N,N-1,...,ICOL+1 are empty. Set the
       *  corresponding elements of ka accordingly. */
      for (j = icol; j <= n - 1; ++j)
        ka[j] = nnz;
    }

  /* Read the bounds */
  fscanf(fpin, " %*[^\n]");
  for (i = 0; i < n + m; ++i)
    fscanf(fpin, "%lf", &bl[i]);
  fscanf(fpin, " %*[^\n]");
  for (i = 0; i < n + m; ++i)
    fscanf(fpin, "%lf", &bu[i]);

  /* Read the initial estimate of x */
  fscanf(fpin, " %*[^\n]");
  for (i = 0; i < n; ++i)
    fscanf(fpin, "%lf", &xs[i]);

  /* Initialize the options structure */
  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options);
  strcpy(options.outfile, outfile);
  /* Read some option values from standard input */
  /* nag_opt_read (e04xyc).
   * Read options from a text file
   */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_read("e04ugc", optionsfile, &options, (Nag_Boolean) Nag_TRUE,
               options.outfile, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  /* Set some other options directly */
  options.major_iter_lim = 100;
  options.crnames = crnames;

  /* Solve the problem. */
  /* nag_opt_nlp_sparse (e04ugc), see above. */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_nlp_sparse(confun, objfun, n, m,
                     ncnln, nonln, njnln, iobj, nnz,
                     a, ha, ka, bl, bu, xs,
                     &ninf, &sinf, &obj, NAGCOMM_NULL,
                     &options, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_nlp_sparse (e04ugc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }

  /* We perturb the solution and solve the
   * same problem again using a warm start.
   */
  fprintf(fpout, "\n\n\nA run of the same example with a warm start:\n");
  fprintf(fpout, "--------------------------------------------\n");
  options.start = Nag_Warm;

  /* Modify some printing options */
  options.print_deriv = Nag_D_NoPrint;
  options.print_level = Nag_Iter;

  /* Perturb xs */
  for (i = 0; i < n+m; i++)
    xs[i] += 0.2;

  /* Reset multiplier estimates to 0.0 */
  if (ncnln > 0)
    {
      for (i = 0; i < ncnln; i++)
        options.lambda[n+i] = 0.0;
    }
  /* Solve the problem again. */
  /* nag_opt_nlp_sparse (e04ugc), see above. */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_nlp_sparse(confun, objfun, n, m,
                     ncnln, nonln, njnln, iobj, nnz,
                     a, ha, ka, bl, bu, xs,
                     &ninf, &sinf, &obj, NAGCOMM_NULL,
                     &options, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_nlp_sparse (e04ugc).\n%s\n",
              fail.message);
      exit_status = 1;
    }
  /* Free memory allocated by nag_opt_nlp_sparse (e04ugc) to pointers in options
   */
  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "all", &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (a) NAG_FREE(a);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  if (xs) NAG_FREE(xs);
  if (ha) NAG_FREE(ha);
  if (ka) NAG_FREE(ka);
  if (crnames) NAG_FREE(crnames);
  if (names) NAG_FREE(names);
  if (optionsfile) NAG_FREE(optionsfile);
  if (outfile) NAG_FREE(outfile);

  return exit_status;
}

/* Subroutine */
static void NAG_CALL confun(Integer ncnln, Integer njnln, Integer nnzjac,
                            const double x[], double conf[], double conjac[],
                            Nag_Comm *comm)
{
#define CONJAC(I) conjac[(I) -1]
#define CONF(I)   conf[(I) -1]
#define X(I)      x[(I) -1]
  /* Compute the nonlinear constraint functions and their Jacobian. */
  if (comm->flag == 0 || comm->flag == 2)
    {
      CONF(1) = sin(-X(1) - 0.25) * 1e3 + sin(-X(2) - 0.25) * 1e3;
      CONF(2) = sin(X(1) - 0.25) * 1e3 + sin(X(1) - X(2) - 0.25) * 1e3;
      CONF(3) = sin(X(2) - X(1) - 0.25) * 1e3 + sin(X(2) - 0.25) * 1e3;
    }
  if (comm->flag == 1 || comm->flag == 2)
    {
      /* Nonlinear Jacobian elements for column 1.0 */
      CONJAC(1) = cos(-X(1) - 0.25) * -1e3;
      CONJAC(2) = cos(X(1) - 0.25) * 1e3 + cos(X(1) - X(2) - 0.25) * 1e3;
      CONJAC(3) = cos(X(2) - X(1) - 0.25) * -1e3;
      /* Nonlinear Jacobian elements for column 2.0 */
      CONJAC(4) = cos(-X(2) - 0.25) * -1e3;
      CONJAC(5) = cos(X(1) - X(2) - 0.25) * -1e3;
      CONJAC(6) = cos(X(2) - X(1) - 0.25) * 1e3 + cos(X(2) - 0.25) * 1e3;
    }
}

static void NAG_CALL objfun(Integer nonln, const double x[], double *objf,
                            double objgrad[], Nag_Comm *comm)
{
#define OBJGRAD(I) objgrad[(I) -1]
#define X(I)       x[(I) -1]
  /* Compute the nonlinear part of the objective function and its grad */
  if (comm->flag == 0 || comm->flag == 2)
    *objf = X(3) * X(3) * X(3) * 1e-6 + X(4) * X(4) * X(4) * 2e-6 / 3.0;
  if (comm->flag == 1 || comm->flag == 2)
    {
      OBJGRAD(1) = 0.0;
      OBJGRAD(2) = 0.0;
      OBJGRAD(3) = X(3) * X(3) * 3e-6;
      OBJGRAD(4) = X(4) * X(4) * 2e-6;
    }
}