/* nag_opt_sparse_convex_qp(e04nkc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 *
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static void qphess1(Integer ncolh, double x[], double hx[], Nag_Comm *comm);
static void qphess2(Integer ncolh, double x[], double hx[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

static void ex1(void);
static void ex2(void);


int main(void)
{
  /* Two examples are called: ex1() uses the
   * default settings to solve a problem while
   * ex2() solves another problem with some
   * of the optional parameters set by the user.
   */
  Vprintf("e04nkc Example Program Results.\n");
  ex1();
  ex2();
  return EXIT_SUCCESS;
}

static void ex1(void)
{
#define MAXN 10
#define MAXM 10
#define MAXBND MAXN+MAXM
#define MAXNNZ 50

  double a[MAXNNZ], bl[MAXBND], bu[MAXBND];
  double x[MAXBND];
  double sinf, obj;
  Integer i, icol, j, jcol;
  Integer iobj, ncolh;
  Integer m, n, nbnd, nnz;
  Integer ninf;
  Integer ha[MAXNNZ], ka[MAXN+1];
  static NagError fail;


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

  fail.print = TRUE;

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

  /* Read nnz, iobj, ncolh */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%ld", &nnz, &iobj, &ncolh);

  /* Read the matrix and set up ka */
  jcol = 1;
  ka[jcol-1] = 0;
  Vscanf(" %*[^\n]");
  for (i = 0; i < nnz; ++i)
    {
      /* a[i] stores the (ha[i], icol) element of matrix */
      Vscanf("%lf%ld%ld", &a[i], &ha[i], &icol);

      /* Check whether we have started a new column */
      if (icol == jcol+1)
        {
          ka[icol-1] = i;  /* Start of icol-th column in a */
          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; ++j)
            ka[j-1] = i;

          ka[icol-1] = i;
          jcol = icol;
        }
    }
  ka[n] = nnz;

  /* If the last columns are empty, set ka accordingly */
  if (n>icol)
    {
      for (j = icol; j<=n - 1; ++j)
        ka[j]=nnz;
    }

  /* Read the bounds */
  nbnd = n+m;
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bu[i]);

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

  /* Solve the problem */
  e04nkc(n, m, nnz, iobj, ncolh, qphess1, a, ha, ka, bl, bu,
         x, &ninf, &sinf, &obj, E04_DEFAULT, NAGCOMM_NULL, &fail);

  if (fail.code != NE_NOERROR)
    exit(EXIT_FAILURE);

}    /* ex1 */

static void qphess1(Integer ncolh,  double x[], double hx[], Nag_Comm *comm)
{
  /* Function to compute H*x. */
  hx[0] = 2.0*x[0];          /* 2.0  0.0  0.0  0.0  0.0  0.0  0.0  */
  hx[1] = 2.0*x[1];          /* 0.0  2.0  0.0  0.0  0.0  0.0  0.0  */
  hx[2] = 2.0*(x[2] + x[3]); /* 0.0  0.0  2.0  2.0  0.0  0.0  0.0  */
  hx[3] = hx[2];             /* 0.0  0.0  2.0  2.0  0.0  0.0  0.0  */
  hx[4] = 2.0*x[4];          /* 0.0  0.0  0.0  0.0  2.0  0.0  0.0  */
  hx[5] = 2.0*(x[5] + x[6]); /* 0.0  0.0  0.0  0.0  0.0  2.0  2.0  */
  hx[6] = hx[5];             /* 0.0  0.0  0.0  0.0  0.0  2.0  2.0  */

}  /* qphess1 */


/* Example 2 */

/* Declare a data structure for passing sparse Hessian data to qphess2 */

typedef struct
{
  double *hess;
  Integer *khess;
  Integer *hhess;
} HessianData;

static void ex2(void)
{
#define MAXHESSNNZ MAXNNZ
  char names[MAXBND][9];
  char *crnames[MAXBND];
  double a[MAXNNZ], bl[MAXBND], bu[MAXBND];
  double hess[MAXHESSNNZ];
  double x[MAXBND];
  double sinf, obj;
  Integer i, icol, j, jcol;
  Integer iobj, ncolh;
  Integer m, n, nbnd, nnz, nnz_hess;
  Integer ninf;
  Integer ha[MAXNNZ], ka[MAXN+1];
  Integer hhess[MAXHESSNNZ], khess[MAXN+1];
  HessianData hess_data;
  Nag_Comm comm;
  Nag_E04_Opt options;
  static NagError fail;

  Vprintf("\nExample 2: some optional parameters are set.\n");
  Vscanf(" %*[^\n]");

  fail.print = TRUE;

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

  /* Read nnz, iobj, ncolh */
  Vscanf(" %*[^\n]");
  Vscanf("%ld%ld%ld", &nnz, &iobj, &ncolh);

  /* Read the matrix and set up ka */
  jcol = 1;
  ka[jcol-1] = 0;
  Vscanf(" %*[^\n]");
  for (i = 0; i < nnz; ++i)
    {
      /* a[i] stores the (ha[i], icol) element of matrix */
      Vscanf("%lf%ld%ld", &a[i], &ha[i], &icol);

      /* Check whether we have started a new column */
      if (icol == jcol+1)
        {
          ka[icol-1] = i;  /* Start of icol-th column in a */
          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; ++j)
            ka[j-1] = i;

          ka[icol-1] = i;
          jcol = icol;
        }
    }
  ka[n] = nnz;

  /* If the last columns are empty, set ka accordingly */
  if (n>icol)
    {
      for (j = icol; j<=n - 1; ++j)
        ka[j]=nnz;
    }

  /* Read the bounds */
  nbnd = n+m;
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bu[i]);

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

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

  /* Read nnz_hess */
  Vscanf(" %*[^\n]");
  Vscanf("%ld", &nnz_hess);

  /* Read the hessian matrix and set up khess */
  jcol = 1;
  khess[jcol-1] = 0;
  Vscanf(" %*[^\n]");
  for (i = 0; i < nnz_hess; ++i)
    {
      /* hess[i] stores the (hhess[i], icol) element of matrix */
      Vscanf("%lf%ld%ld", &hess[i], &hhess[i], &icol);

      /* Check whether we have started a new column */
      if (icol == jcol+1)
        {
          khess[icol-1] = i;  /* Start of icol-th column in hess */
          jcol = icol;
        }
      else if (icol > jcol+1)
        {
          /* Index in hess 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 khess to i.
           */
          for (j = jcol+1; j < icol; ++j)
            khess[j-1] = i;

          khess[icol-1] = i;
        }
    }
  khess[ncolh] = nnz_hess;

  /* Initialize options structure */
  e04xxc(&options);
  options.crnames = crnames;

  /* Package up Hessian data for communication via comm */
  hess_data.hess = hess;
  hess_data.khess = khess;
  hess_data.hhess = hhess;

  comm.p = (Pointer)&hess_data;

  /* Solve the problem */
  e04nkc(n, m, nnz, iobj, ncolh, qphess2, a, ha, ka, bl, bu,
         x, &ninf, &sinf, &obj, &options, &comm, &fail);

  if (fail.code == NE_NOERROR)
    {
      Vprintf("\nPerturb the problem and re-solve with warm start.\n");
      for (i = 0; i < nnz_hess; ++i)
        hess[i] *= 1.001;

      options.start = Nag_Warm;
      options.print_level = Nag_Soln;
      e04nkc(n, m, nnz, iobj, ncolh, qphess2, a, ha, ka, bl, bu,
             x, &ninf, &sinf, &obj, &options, &comm, &fail);
    }

  /* Free memory NAG-allocated to members of options */
  e04xzc(&options, "", NAGERR_DEFAULT);

  if (fail.code != NE_NOERROR)
    exit(EXIT_FAILURE);

}  /* ex2 */

static void qphess2(Integer ncolh,  double x[], double hx[], Nag_Comm *comm)
{
  Integer i, j, jrow;
  HessianData *hd = (HessianData *)(comm->p);
  double *hess = hd->hess;
  Integer *hhess = hd->hhess;
  Integer *khess = hd->khess;


  for (i = 0; i < ncolh; ++i)
    hx[i] = 0.0;

  for (i = 0; i < ncolh; ++i)
    {
      /* For each column of Hessian */
      for (j = khess[i]; j < khess[i+1]; ++j)
        {
          /* For each non-zero in column */

          jrow = hhess[j] - 1;

          /* Using symmetry of hessian, add contribution
           * to hx of hess[j] as a diagonal or upper
           * triangular element of hessian.
           */
          hx[i] += hess[j]*x[jrow];

          /* If hess[j] not a diagonal element add its
           * contribution to hx as a lower triangular
           * element of hessian.
           */
          if (jrow != i)
            hx[jrow] += hess[j]*x[i];
        }
    }

}  /* qphess2 */