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

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

#ifdef E04NKC
  static void qphess1(Integer ncolh, double x[], double hx[], Nag_Comm *comm);
#else
  static void qphess1(Integer ncolh, const double x[], double hx[], Integer nstate, Nag_Comm *comm);
#endif

int main(void)
{

  Vprintf("e04nkc Example Program Results.\n");
#define MAXNNZ 50

  Integer exit_status=0, i, icol, iobj, j, jcol, *ka=0, m, n, nbnd;
  Integer lenc, ncolh, ninf, nname, ns;
  Integer   nnz;
  NagError fail;
  double  *x, obj, sinf, objadd ;
  double a[]={0.02, 0.02, 0.03, 1.00, 0.70, 0.02, 0.15, -200.00, 0.06,
		0.75, 0.03, 0.04, 0.05, 0.04, 1.00, -2000.00, 0.02,
		1.00, 0.01, 0.08, 0.08, 0.80, -2000.00, 1.00, 0.12,
		0.02, 0.02, 0.75, 0.04, -2000.00, 0.01, 0.80, 0.02,
		1.00, 0.02, 0.06, 0.02, -2000.00, 1.00, 0.01, 0.01,
		0.97, 0.01, 400.00, 0.97, 0.03, 1.00, 400.00};
  Integer ha[]={7, 5, 3, 1, 6, 4, 2, 8, 7, 6, 5, 4, 3, 2, 1, 8, 2, 1,
		  4, 3, 7, 6, 8, 1, 7, 3, 4, 6, 2, 8, 5, 6, 7, 1, 2, 3,
		  4, 8, 1, 2, 3, 6, 7, 8, 7, 2, 1, 8};
  Integer icolarray[48]={1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,
		       3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6,
		       6, 6, 6, 6, 7, 7, 7, 7};
  char *names[] = {"...X1...", "...X2...", "...X3...", "...X4...", "...X5...",
		   "...X6...", "...X7...", "..ROW1..", "..ROW2..", "..ROW3..",
		   "..ROW4..", "..ROW5..", "..ROW6..", "..ROW7..", "..COST.."};
  double bl[] = {0.0, 0.0, 4.0e+02, 1.0e+02, 0.0, 0.0, 0.0, 2.0e+03, 
		    -1.0e+25, -1.0e+25, -1.0e+25, -1.0e+25, 1.5e+03, 2.5e+02, -1.0e+25};
  double bu[] = {2.0e+02, 2.5e+03, 8.0e+02, 7.0e+02, 1.5e+03, 1.0e+25, 1.0e+25, 2.0e+03, 
		   6.0e+01, 1.0e+02, 4.0e+01, 3.0e+01, 1.0e+25, 3.0e+02, 1.0e+25};
  Nag_E04State state;
  Nag_Start start;
  Nag_Comm comm;
  Nag_FileID fileid;
  char prob[9];
  Integer *helast=0, *hs=0, *inda=0, *loca=0;
  double *acol=0, *c=0, *pi=0, *rc=0;

  INIT_FAIL(fail);

  /* Problem dimensions */
  n = 7;
  m = 8;

  /* nnz, iobj, ncolh */
  nnz = 48;
  iobj = 8;
  ncolh = 7;

  if (n>=1 &&  m>=1 && nnz>=1 && nnz<=n*m)
    {
      nbnd = n + m;
      if ( !( x = NAG_ALLOC(nbnd, double)) ||
           !(c = NAG_ALLOC(1, double)) ||
           !(pi = NAG_ALLOC(m, double)) ||
           !(rc = NAG_ALLOC(n+m, double)) ||
           !( ka = NAG_ALLOC(n+1, Integer)) ||
           !(helast = NAG_ALLOC(n+m, Integer)) ||
           !(hs = NAG_ALLOC(n+m, Integer)) 
	   ) 
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      Vprintf("Invalid n or m or nnz.\n");
      exit_status = 1;
      return exit_status;
    }

  /* Read the matrix and set up ka */
  jcol = 1;
  ka[jcol-1] = 0;
  for (i = 0; i < nnz; ++i)
    {
      /* a[i] stores the (ha[i], icol) element of matrix */
      icol = icolarray[i];
 
      /* 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;
    }

  /* The initial estimate of x */
  for (i = 0; i < n; ++i)
    x[i] = 0.0;

  /* Solve the problem */
#ifdef E04NKC
  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)
           {
             Vprintf("Error from e04nkc.\n%s\n", fail.message);
             exit_status = 1;
             goto END;
           }
#else
  /* Call e04npc to initialise e04nqc. */
  e04npc(&state, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Initialisation of e04nqc failed.\n");
      exit_status = 1;
      goto END;
    }
  /* By default e04nqc does not print monitoring */
  /* information. Call x04acc to set the print file fileid */
  x04acc("", 2, &fileid, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Fileid could not be obtained.\n");
      exit_status = 1;
      goto END;
    }

  e04ntc("Print file", fileid, &state, &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Files stream could not be set.\n");
      exit_status = 1;
      goto END;
    }

  /* We have no explicit objective vector so set lenc = 0; the 
  * objective vector is stored in row iobj of a.
  */
  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;
      hs[i-1] = 0.0;
    }

  /* Solve the QP problem. */
  start = Nag_Cold;
  nname = 15;
  for (i=0; i<n+1; i++)
    ka[i]++;
  e04nqc(start, qphess1, m, n, nnz, nname, lenc, ncolh, iobj,
	 objadd, prob, a, ha, ka, bl, bu, c, names, helast, 
	 hs, x, pi, rc, &ns, &ninf, &sinf, &obj, &state, &comm, &fail);
        
  Vprintf("\n");
  Vprintf("On exit from e04nqc, fail.message = %s\n", fail.message);
  if (fail.code == NE_NOERROR)
    {
      Vprintf("Final objective value = %11.3e\n",   obj);
      Vprintf("Optimal X = ");

      for (i = 1; i <= n; ++i)
        {
          Vprintf("%9.2f%s", x[i - 1], i%7 == 0 || i == n ?"\n":" ");
        }
    }
#endif
 END:
  if (c) NAG_FREE(c);
  if (pi) NAG_FREE(pi);
  if (rc) NAG_FREE(rc);
  if (helast) NAG_FREE(helast);
  if (x) NAG_FREE(x);
  if (ka) NAG_FREE(ka);
  if (hs) NAG_FREE(hs);
  return exit_status;
}

#ifdef E04NKC
  static void qphess1(Integer ncolh, double x[], double hx[], Nag_Comm *comm)
#else
static void qphess1(Integer ncolh,  const double x[], double hx[], Integer nstate, Nag_Comm *comm)
#endif
{
#define HX(I) hx[(I)-1]
#define X(I) x[(I)-1]

  /* 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);


}  /* qphess1 */