/* nag_opt_qp (e04nfc) Example Program.
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

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

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

int main(void)
{
  Integer exit_status_ex1=0;
  Integer exit_status_ex2=0;

  /* 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("nag_opt_qp (e04nfc) Example Program Results.\n");

  exit_status_ex1 = ex1();
  exit_status_ex2 = ex2();

  return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1;
}

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

static int ex1(void)
{
  Integer exit_status=0, i, j, n, nbnd, nclin, tda, tdh;
  NagError fail;
  double *a=0, bigbnd, *bl=0, *bu=0, *cvec=0, *h=0, objf, *x=0;

  INIT_FAIL(fail);

  Vprintf("\nExample 1: default options used.\n");

  /* Define the problem. This example is due to Bunch and Kaufman,
   * `A computational method for the indefinite quadratic programming
   * problem ', Linear Algebra and its Applications, 34, 341-370 (1980).
   *
   * h  = the QP Hessian matrix.
   * a  = the general 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.
   *
   * Set the actual problem dimensions.
   * n      = the number of variables.
   * nclin  = the number of general linear constraints (may be 0).
   */
  n = 8;
  nclin = 7;
  if (n>0 && nclin>=0)
    {
      nbnd = nclin + n;
      if ( !( a = NAG_ALLOC(nclin*n, double)) ||
           !( h = NAG_ALLOC(n*n, double)) ||
           !( x = NAG_ALLOC(n, double)) ||
           !( cvec = NAG_ALLOC(n, double)) ||
           !( bl = NAG_ALLOC(nbnd, double)) ||
           !( bu = NAG_ALLOC(nbnd, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdh = n;
    }
  else
    {
      Vprintf("Invalid n or nclin.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Define the value used to denote ``infinite'' bounds. */
  bigbnd = 1.0e20;

  for (i = 0; i < nclin; ++i)
    for (j = 0; j < n; ++j)
      A(i,j) = 0.0;

  for (i = 0; i < nclin; ++i)
    {
      A(i,i) = -1.0;
      A(i,i+1) = 1.0;
      bl[n + i] = -1.0 - 0.05*(double)i;
      bu[n + i] = bigbnd;
    }

  for (j = 0; j < n; ++j)
    {
      bl[j] = -(double)(j+1) - 0.1*(double)(j);
      bu[j] = (double)(j+1);
      cvec[j] = (double)(7 - j);
    }

  for (i = 0; i < n; ++i)
    {
      for (j = i+1; j < n; ++j)
        H(i,j) = (double)(ABS(i-j));
      H(i,i) = 1.69;
    }

  /* Set the initial estimate of the solution. */
  x[0] = -1.0;
  x[1] = -2.0;
  x[2] = -3.0;
  x[3] = -4.0;
  x[4] = -5.0;
  x[5] = -6.0;
  x[6] = -7.0;
  x[7] = -8.0;

  /* Solve the QP problem. */
  /* nag_opt_qp (e04nfc).
   * Quadratic programming
   */
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh,
             NULLFN, x, &objf, E04_DEFAULT, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  Vprintf("Re-solve problem with the Hessian defined by function qphess1.\n");

  /* Set a new initial estimate of the solution. */
  x[0] = -1.0;
  x[1] = 12.0;
  x[2] = -3.0;
  x[3] = 14.0;
  x[4] = -5.0;
  x[5] = 16.0;
  x[6] = -7.0;
  x[7] = 18.0;

  /* Solve the QP problem. */
  /* nag_opt_qp (e04nfc), see above. */
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *)0, tdh,
             qphess1, x, &objf, E04_DEFAULT, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (a) NAG_FREE(a);
  if (h) NAG_FREE(h);
  if (x) NAG_FREE(x);
  if (cvec) NAG_FREE(cvec);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  return exit_status;
}      /* ex1 */
static void qphess1(Integer n, Integer jthcol, double h[], Integer tdh,
                    double x[], double hx[], Nag_Comm *comm)
{
  /* qphess1 computes the vector  Hx = (H)*x  for some matrix H
   * that defines the Hessian of the required QP problem.
   *
   * In this version qphess the Hessian matrix is implicit.
   * The array h[] is not accessed. There is no special coding
   * for the case jthcol > 0
   */
  Integer i, j;
  double sum;

  for (i = 0; i < n; ++i)
    {
      sum = 1.69*x[i];
      for (j = 0; j < n; ++j)
        sum += x[j]*(double)ABS(i - j);
      hx[i] = sum;
    }
}                               /* qphess1 */

static int ex2(void)
{
  Nag_Boolean print;
  Integer exit_status=0, i, j, n, nbnd, nclin, tda, tdh;
  NagError fail;
  Nag_E04_Opt options;
  double *a=0, *bl=0, *bu=0, *cvec=0, *h=0, objf, *x=0;

  INIT_FAIL(fail);


  Vprintf("\nExample 2: some optional parameters are set.\n");
  Vscanf(" %*[^\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)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdh = n;
    }
  else
    {
      Vprintf("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). */
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < n; ++i)
    Vscanf("%lf",&cvec[i]);

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

  /* Read the bounds. */
  nbnd = n + nclin;
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nbnd; ++i)
    Vscanf("%lf", &bl[i]);
  Vscanf(" %*[^\n]"); /* Skip heading in data file */
  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]);

  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options); /* Initialise 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", "stdin", &options, print, "stdout", &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("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 qphess2.
   */
  /* nag_opt_qp (e04nfc), see above. */
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *)0, tdh,
             qphess2, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("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 qphess3().
   * Only the final solution from the results is printed.
   */
  Vprintf("\nA run of the same example with a warm start:\n");

  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,
             qphess3, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("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)
    {
      Vprintf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (x) NAG_FREE(x);
  if (cvec) NAG_FREE(cvec);
  if (a) NAG_FREE(a);
  if (h) NAG_FREE(h);
  if (bl) NAG_FREE(bl);
  if (bu) NAG_FREE(bu);
  return exit_status;

} /* ex2 */
static void qphess2(Integer n, Integer jthcol, double h[], Integer tdh,
                    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.
   */

  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];
} /* qphess2 */

#undef H
static void qphess3(Integer n, Integer jthcol, double h[], Integer tdh,
                    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 (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];
    }
} /* qphess3 */