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

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

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

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

int main(void)
{
  /* Two examples are called: ex1() uses the
   * default settings to solve a problem while
   * ex2() solves the same problem with some
   * of the optional parameters set by the user,
   * perturbs the solution and then uses the
   * warm start facility.
   */
  Vprintf("e04ugc Example Program Results\n");
  if (ex1() != 0) return 1;
  if (ex2() != 0) return 1;
  return 0;
}

static int ex1(void)
{
  double  *a=0, *bl=0, *bu=0, *xs=0;
  Integer *ha=0, *ka=0;

  Integer exit_status=0;

  double obj, sinf;
  Integer m, n, nnz, ncnln, njnln, nonln;
  Integer ninf;
  Integer i, icol, j, jcol;
  Integer iobj;
  static NagError fail;


  Vprintf("\nExample 1: default options used.\n");
  Vscanf(" %*[^\n]"); /* Skip heading in data file. */
  Vscanf(" %*[^\n]");

  fail.print = TRUE;

  /* Read the problem dimensions */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld", &n, &m);

  /* Read NCNLN, NONLN and NJNLN from data file. */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%ld", &ncnln, &nonln, &njnln);

  /* Read NNZ, IOBJ */
  Vscanf(" %*[^\n]");
  Vscanf("%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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }
  /* read the matrix and set up ka. */
  jcol = 1;
  ka[jcol - 1] = 0;
  Vscanf(" %*[^\n]");
  for (i = 0; i < nnz; ++i)
    {
      /* a[i] stores (ha[i], icol) element of matrix */
      Vscanf("%lf%ld%ld", &a[i], &ha[i], &icol);
      if (icol < jcol)
        {
          /* Elements not ordered by increasing column index. */
          Vprintf("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 */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + m; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + m; ++i)
    Vscanf("%lf", &bu[i]);

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

  /* Solve the problem. */
  e04ugc  (confun, objfun, n, m,
           ncnln, nonln, njnln, iobj, nnz,
           a, ha, ka, bl, bu, xs,
           &ninf, &sinf, &obj, NAGCOMM_NULL,
           E04_DEFAULT, &fail);

 END:
  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);
  return exit_status;
}

/* Subroutine */
static void 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 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;
    }
}

#define MAXNAMES 300
static int ex2(void)
{
  char names   [MAXNAMES][9];
  char *crnames[MAXNAMES];

  double  *a=0 , *bl=0, *bu=0, *xs=0;
  Integer *ha=0, *ka=0;

  Integer exit_status=0;

  double obj, sinf;
  Integer m, n, nnz, ncnln, njnln, nonln;
  Integer ninf;
  Integer i, icol, j, jcol;
  Integer iobj;
  static NagError fail;
  Nag_E04_Opt options;


  Vprintf("\nExample 2: Use of the option structure.\n");
  Vscanf(" %*[^\n]");

  fail.print = TRUE;

  /* Read the problem dimensions */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld", &n, &m);

  /* Read NCNLN, NONLN and NJNLN from data file. */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%ld", &ncnln, &nonln, &njnln);

  /* Read NNZ, IOBJ */
  Vscanf(" %*[^\n]");
  Vscanf("%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)) )
    {
      Vprintf("Allocation failure\n");
      exit_status = 1;
      goto END;
    }

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


  /* read the matrix and set up ka. */
  jcol = 1;
  ka[jcol - 1] = 0;
  Vscanf(" %*[^\n]");
  for (i = 0; i < nnz; ++i)
    {
      /* a[i] stores (ha[i], icol) element of matrix */
      Vscanf("%lf%ld%ld", &a[i], &ha[i], &icol);
      if (icol < jcol)
        {
          /* Elements not ordered by increasing column index. */
          Vprintf("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 */
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + m; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]");
  for (i = 0; i < n + m; ++i)
    Vscanf("%lf", &bu[i]);

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

  /* Initialize the options structure */
  e04xxc(&options);

  /* Read some option values from standard input */
  e04xyc("e04ugc", "stdin", &options, (Boolean)TRUE,
         "stdout", NAGERR_DEFAULT);

  /* Set some other options directly */
  options.major_iter_lim = 100;
  options.crnames = crnames;

  /* Solve the problem. */
  e04ugc  (confun, objfun, n, m,
           ncnln, nonln, njnln, iobj, nnz,
           a, ha, ka, bl, bu, xs,
           &ninf, &sinf, &obj, NAGCOMM_NULL,
           &options, &fail);

  if (fail.code == NE_NOERROR)
    {
      /* We perturb the solution and solve the
       * same problem again using a warm start.
       */
      Vprintf("\n\n\nA run of the same example with a warm start:\n");
      Vprintf("--------------------------------------------\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. */
      e04ugc  (confun, objfun, n, m,
               ncnln, nonln, njnln, iobj, nnz,
               a, ha, ka, bl, bu, xs,
               &ninf, &sinf, &obj, NAGCOMM_NULL,
               &options, &fail);
    }
  exit_status = fail.code;

  /* Free memory allocated by e04ugc to pointers in options */
  e04xzc(&options, "all", NAGERR_DEFAULT);

END:
  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);
  return exit_status;
}