/* nag_opt_qp (e04nfc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 *
 */

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

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL qphess1(Integer n, Integer jthcol, const double h[],
                               Integer tdh, const double x[], double hx[],
                               Nag_Comm *comm);
  static void NAG_CALL qphess2(Integer n, Integer jthcol, const double h[],
                               Integer tdh, const double x[], double hx[],
                               Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

#define A(I, J) a[(I) *tda + J]
#define H(I, J) h[(I) *tdh + J]

int main(void)
{
  const char *optionsfile = "e04nfce.opt";
  static double ruser[2] = { -1.0, -1.0 };
  Nag_Boolean print;
  Integer exit_status = 0, i, j, n, nbnd, nclin, tda, tdh;
  Nag_E04_Opt options;
  double *a = 0, *bl = 0, *bu = 0, *cvec = 0, *h = 0, objf, *x = 0;
  Nag_Comm comm;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_opt_qp (e04nfc) Example Program Results\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

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

  /* Set the actual problem dimensions.
   * n = the number of variables.
   * nclin = the number of general linear constraints (may be 0).
   */
  n = 7;
  nclin = 7;
  if (n > 0 && nclin >= 0) {
    nbnd = n + nclin;
    if (!(x = NAG_ALLOC(n, double)) ||
        !(cvec = NAG_ALLOC(n, double)) ||
        !(a = NAG_ALLOC(nclin * n, double)) ||
        !(h = NAG_ALLOC(n * n, double)) ||
        !(bl = NAG_ALLOC(nbnd, double)) || !(bu = NAG_ALLOC(nbnd, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    tda = n;
    tdh = n;
  }
  else {
    printf("Invalid n or nclin.\n");
    exit_status = 1;
    return exit_status;
  }
  /* cvec   = the coefficients of the explicit linear term of f(x).
   * a      = the linear constraint matrix.
   * bl     = the lower bounds on x and A*x.
   * bu     = the upper bounds on x and A*x.
   * x      = the initial estimate of the solution.
   */

  /* Read the coefficients of the explicit linear term of f(x). */
  scanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < n; ++i)
    scanf("%lf", &cvec[i]);

  /* Read the linear constraint matrix A. */
  scanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nclin; ++i)
    for (j = 0; j < n; ++j)
      scanf("%lf", &A(i, j));

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

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

  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options); /* Initialize options structure */
  /* Set one option directly
   * Bounds  >=   inf_bound will be treated as plus  infinity.
   * Bounds  <=  -inf_bound will be treated as minus infinity.
   */
  options.inf_bound = 1.0e21;

  /* Read remaining option values from file */
  print = Nag_TRUE;
  /* nag_opt_read (e04xyc).
   * Read options from a text file
   */
  nag_opt_read("e04nfc", optionsfile, &options, print, "stdout", &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_read (e04xyc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Solve the problem from a cold start.
   * The Hessian is defined implicitly by function qphess1.
   */
  /* nag_opt_qp (e04nfc), see above. */
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *) 0, tdh,
             qphess1, x, &objf, &options, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* The following is for illustrative purposes only. We do a warm
   * start with the final working set of the previous run.
   * This time we store the Hessian explicitly in h[][], and use
   * the corresponding function qphess2().
   * Only the final solution from the results is printed.
   */
  printf("\nA run of the same example with a warm start:\n");
  fflush(stdout);

  options.start = Nag_Warm;
  options.print_level = Nag_Soln;

  for (i = 0; i < n; ++i) {
    for (j = 0; j < n; ++j)
      H(i, j) = 0.0;
    if (i <= 4)
      H(i, i) = 2.0;
    else
      H(i, i) = -2.0;
  }
  H(2, 3) = 2.0;
  H(3, 2) = 2.0;
  H(5, 6) = -2.0;
  H(6, 5) = -2.0;

  /* Solve the problem again. */
  /* nag_opt_qp (e04nfc), see above. */
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh,
             qphess2, x, &objf, &options, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
    exit_status = 1;
  }
  /* Free memory allocated by nag_opt_qp (e04nfc) 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) {
    printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

END:
  NAG_FREE(x);
  NAG_FREE(cvec);
  NAG_FREE(a);
  NAG_FREE(h);
  NAG_FREE(bl);
  NAG_FREE(bu);

  return exit_status;
}

static void NAG_CALL qphess1(Integer n, Integer jthcol, const double h[],
                             Integer tdh, const double x[], double hx[],
                             Nag_Comm *comm)
{
  /* In this version of qphess the Hessian matrix is implicit.
   * The array h[] is not accessed. There is no special coding
   * for the case jthcol > 0.
   */

  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback qphess1, first invocation.)\n");
    fflush(stdout);
    comm->user[0] = 0.0;
  }

  hx[0] = 2.0 * x[0];
  hx[1] = 2.0 * x[1];
  hx[2] = 2.0 * (x[2] + x[3]);
  hx[3] = hx[2];
  hx[4] = 2.0 * x[4];
  hx[5] = -2.0 * (x[5] + x[6]);
  hx[6] = hx[5];
} /* qphess1 */

#undef H

static void NAG_CALL qphess2(Integer n, Integer jthcol, const double h[],
                             Integer tdh, const double x[], double hx[],
                             Nag_Comm *comm)
{
  /* In this version of qphess, the matrix H is stored in h[]
   * as a full two-dimensional array.
   */

#define H(I, J) h[(I) *tdh + (J)]

  Integer i, j;

  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback qphess2, first invocation.)\n");
    fflush(stdout);
    comm->user[1] = 0.0;
  }

  if (jthcol != 0) {
    /* Special case -- extract one column of  H. */
    j = jthcol - 1;
    for (i = 0; i < n; ++i)
      hx[i] = H(i, j);
  }
  else {
    /* Normal Case. */
    for (i = 0; i < n; ++i)
      hx[i] = 0.0;

    for (i = 0; i < n; ++i)
      for (j = 0; j < n; ++j)
        hx[i] += H(i, j) * x[j];
  }
} /* qphess2 */