/* 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 <string.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nage04.h>
#include <nagx04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL qphess1(Integer n, Integer jthcol, double h[],
                             Integer tdh, double x[], double hx[],
                             Nag_Comm *comm);
static void NAG_CALL qphess2(Integer n, Integer jthcol, double h[],
                             Integer tdh, 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(int argc, char *argv[])
{
  FILE        *fpin, *fpout;
  char        *optionsfile = 0;
  char        *outfile = 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;
  NagError    fail;

  INIT_FAIL(fail);


  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);
  (void) nag_example_file_io(argc, argv, "-options", &optionsfile);
  (void) nag_example_file_io(argc, argv, "-nag_write", &outfile);
  if (!outfile)
    {
      outfile = NAG_ALLOC(7, char);
      strcpy(outfile, "stdout");
    }

  fprintf(fpout, "nag_opt_qp (e04nfc) Example Program Results\n");

  fscanf(fpin, " %*[^\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)))
        {
          fprintf(fpout, "Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdh = n;
    }
  else
    {
      fprintf(fpout, "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). */
  fscanf(fpin, " %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < n; ++i)
    fscanf(fpin, "%lf", &cvec[i]);

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

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

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

  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options); /* Initialise options structure */
  strcpy(options.outfile, outfile);
  /* 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
   */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_read("e04nfc", optionsfile, &options, print, options.outfile, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "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. */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *) 0, tdh,
             qphess1, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "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.
   */
  fprintf(fpout, "\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. */
  if (strcmp(outfile, "stdout"))
    fclose(fpout);
  nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh,
             qphess2, x, &objf, &options, NAGCOMM_NULL, &fail);
  if (strcmp(outfile, "stdout"))
    {
      fpout = fopen(outfile, "a");
    }
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "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)
    {
      fprintf(fpout, "Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  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);
  if (optionsfile) NAG_FREE(optionsfile);
  if (outfile) NAG_FREE(outfile);

  return exit_status;
}

static void NAG_CALL qphess1(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];
} /* qphess1 */

#undef H

static void NAG_CALL qphess2(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];
    }
} /* qphess2 */