/* nag_opt_sparse_convex_qp_solve (e04nqc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

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

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

int main(void)
{
  /* Scalars */
  double obj, objadd, sinf;
  Integer exit_status, i, icol, iobj, j, jcol, lenc, m, n, ncolh, ne, ninf;
  Integer nname, ns;
  Integer verbose_output;

  /* Arrays */
  char nag_enum_arg[40];
  char prob[9];
  char **names;
  double *acol = 0, *bl = 0, *bu = 0, *c = 0, *pi = 0, *rc = 0, *x = 0;
  Integer *helast = 0, *hs = 0, *inda = 0, *loca = 0;

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

  exit_status = 0;
  INIT_FAIL(fail);

  printf("nag_opt_sparse_convex_qp_solve (e04nqc) Example Program Results\n");
  fflush(stdout);

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

  /* Read ne, iobj, ncolh, start and nname from data file. */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &n, &m);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT " %39s %" NAG_IFMT "%*[^\n] ",
        &ne, &iobj, &ncolh, nag_enum_arg, &nname);
  /* nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  start = (Nag_Start) nag_enum_name_to_value(nag_enum_arg);
  if (n >= 1 && m >= 1) {
    /* Allocate memory */
    if (!(names = NAG_ALLOC(n + m, char *)) ||
        !(acol = NAG_ALLOC(ne, double)) ||
        !(bl = NAG_ALLOC(m + n, double)) ||
        !(bu = NAG_ALLOC(m + n, double)) ||
        !(c = NAG_ALLOC(1, double)) ||
        !(pi = NAG_ALLOC(m, double)) ||
        !(rc = NAG_ALLOC(n + m, double)) ||
        !(x = NAG_ALLOC(n + m, double)) ||
        !(helast = NAG_ALLOC(n + m, Integer)) ||
        !(hs = NAG_ALLOC(n + m, Integer)) ||
        !(inda = NAG_ALLOC(ne, Integer)) ||
        !(loca = NAG_ALLOC(n + 1, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("%s", "Either m or n invalid\n");
    exit_status = 1;
    return exit_status;
  }

  /* Read names from data file. */
  for (i = 1; i <= nname; ++i) {
    names[i - 1] = NAG_ALLOC(9, char);
    scanf(" ' %8s '", names[i - 1]);
  }
  scanf("%*[^\n] ");

  /* Read the matrix acol from data file. Set up LOCA. */
  jcol = 1;
  loca[jcol - 1] = 1;
  for (i = 1; i <= ne; ++i) {
    /* Element (inda[i-1], icol) is stored in acol[i-1]. */
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &acol[i - 1], &inda[i - 1],
          &icol);
    if (icol < jcol) {
      /* Elements not ordered by increasing column index. */
      printf("%s%5" NAG_IFMT "%s%5" NAG_IFMT "%s%s\n", "Element in column",
             icol, " found after element in column", jcol, ". Problem",
             " abandoned.");
    }
    else if (icol == jcol + 1) {
      /* Index in ACOL of the start of the ICOL-th column equals I. */
      loca[icol - 1] = i;
      jcol = icol;
    }
    else if (icol > jcol + 1) {
      /* Index in acol 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 loca to i. */
      for (j = jcol + 1; j <= icol - 1; ++j) {
        loca[j - 1] = i;
      }
      loca[icol - 1] = i;
      jcol = icol;
    }
  }
  loca[n] = ne + 1;

  if (n > icol) {
    /* Columns n,n-1,...,icol+1 are empty. Set the corresponding */
    /* elements of loca accordingly. */
    for (i = n; i >= icol + 1; --i) {
      loca[i - 1] = loca[i];
    }
  }

  /* Read bl, bu, hs and x from data file. */
  for (i = 1; i <= n + m; ++i) {
    scanf("%lf", &bl[i - 1]);
  }
  scanf("%*[^\n] ");

  for (i = 1; i <= n + m; ++i) {
    scanf("%lf", &bu[i - 1]);
  }
  scanf("%*[^\n] ");

  if (start == Nag_Cold) {
    for (i = 1; i <= n; ++i) {
      scanf("%" NAG_IFMT "", &hs[i - 1]);
    }
    scanf("%*[^\n] ");
  }
  else if (start == Nag_Warm) {
    for (i = 1; i <= n + m; ++i) {
      scanf("%" NAG_IFMT "", &hs[i - 1]);
    }
    scanf("%*[^\n] ");
  }
  for (i = 1; i <= n; ++i) {
    scanf("%lf", &x[i - 1]);
  }
  scanf("%*[^\n] ");

  printf("\nQP problem contains %3ld variables and %3"NAG_IFMT
         " linear constraints\n", n, m);

  /* nag_opt_sparse_convex_qp_init (e04npc).
   * Initialization function for
   * nag_opt_sparse_convex_qp_solve (e04nqc)
   */
  nag_opt_sparse_convex_qp_init(&state, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Initialization of "
           "nag_opt_sparse_convex_qp_solve (e04nqc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set this to 1 to cause e04nqc to produce intermediate
     progress output */
  verbose_output = 0;

  if (verbose_output)
    {
      /* By default nag_opt_sparse_convex_qp_solve (e04nqc) 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) {
        exit_status = 2;
        goto END;
      }

      /* nag_opt_sparse_convex_qp_option_set_integer (e04ntc).
       * Set a single option for nag_opt_sparse_convex_qp_solve (e04nqc)
       * from an integer argument
       */
      nag_opt_sparse_convex_qp_option_set_integer("Print file", fileid, &state,
                                                  &fail);
      if (fail.code != NE_NOERROR) {
        exit_status = 1;
        goto END;
      }
    }

  /* We have no explicit objective vector so set lenc = 0; the
   * objective vector is stored in row iobj of acol.
   */
  lenc = 0;
  objadd = 0.;
  strcpy(prob, "        ");

  /* Do not allow any elastic variables (i.e. they cannot be */
  /* infeasible). If we'd set optional argument "Elastic mode" to 0, */
  /* we wouldn't need to set the individual elements of array helast. */
  for (i = 1; i <= n + m; ++i) {
    helast[i - 1] = 0;
  }

  /* Illustrate how to pass information to the user-supplied
     function qphx via the comm structure */
  comm.p = 0;

  /* Solve the QP problem. */
  /* nag_opt_sparse_convex_qp_solve (e04nqc).
   * LP or QP problem (suitable for sparse problems)
   */
  nag_opt_sparse_convex_qp_solve(start, qphx, m, n, ne, nname, lenc, ncolh,
                                 iobj, objadd, prob, acol, inda, loca, bl, bu,
                                 c,
                                 (const char **) names, helast, hs, x, pi, rc,
                                 &ns, &ninf, &sinf, &obj, &state, &comm,
                                 &fail);

  if (fail.code != NE_NOERROR) {
    printf("nag_opt_sparse_convex_qp_solve (e04nqc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  if (fail.code != NE_NOERROR) {
    exit_status = 2;
    goto END;
  }

  printf("\n");
  printf("Final objective value = %12.3e\n", obj);
  printf("Optimal X = ");

  for (i = 1; i <= n; ++i) {
    printf("%9.2f%s", x[i - 1], i % 7 == 0 || i == n ? "\n" : " ");
  }

END:

  for (i = 0; i < n + m; i++) {
    NAG_FREE(names[i]);
  }
  NAG_FREE(names);
  NAG_FREE(acol);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(c);
  NAG_FREE(pi);
  NAG_FREE(rc);
  NAG_FREE(x);
  NAG_FREE(helast);
  NAG_FREE(hs);
  NAG_FREE(inda);
  NAG_FREE(loca);
  return exit_status;
}

static void NAG_CALL qphx(Integer ncolh, const double x[], double hx[],
                          Integer nstate, Nag_Comm *comm)
{
  /* Routine to compute H*x. (In this version of qphx, the Hessian
   * matrix H is not referenced explicitly.)
   */

  /* Parameter adjustments */
#define HX(I) hx[(I) -1]
#define X(I)  x[(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 qphx via the comm struct\n", comm->p);

  /* Function Body */
  HX(1) = X(1) * 2;
  HX(2) = X(2) * 2;
  HX(3) = (X(3) + X(4)) * 2;
  HX(4) = HX(3);
  HX(5) = X(5) * 2;
  HX(6) = (X(6) + X(7)) * 2;
  HX(7) = HX(6);
  return;
} /* qphx */