/* nag_opt_miqp_mps_read (e04mxc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */
#include <stdio.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>
#include <nagx04.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

/* Make a typedef for convenience when allocating crname. */
typedef char e04mx_name[9];

int main(void)
{
  Integer exit_status = 0;
  double obj, objadd, sinf;
  Integer i, iobj, j, lenc, lintvar, m, maxlintvar, maxm, maxn, maxncolh,
    maxnnz, maxnnzh, minmax, n, ncolh, ninf, nname, nnz, nnzh, ns;
  Integer mpslst = 1;
  double *a = 0, *bl = 0, *bu = 0, *h = 0, *pi = 0, *rc = 0, *ruser = 0, *x = 0;
  Integer *helast = 0, *hs = 0, *iccola = 0, *iccolh = 0, *intvar = 0,
    *irowa = 0, *irowh = 0, *iuser = 0;
  char pnames[5][9] = {"", "", "", "", ""};
  char (*crname)[9] = 0;
  char **names = 0;
  char fname[] = "e04mxce.opt";

  /* Nag Types */
  Nag_Boolean readints = Nag_FALSE;
  Nag_E04State state;
  NagError fail;
  Nag_Comm comm;
  Nag_FileID fileid;

  INIT_FAIL(fail);

  printf("nag_opt_miqp_mps_read (e04mxc) Example Program Results\n");
  fflush(stdout);

  /* nag_open_file (x04acc).
     Open unit number for reading and associate unit with named file. */
  nag_open_file(fname, 0, &fileid, NAGERR_DEFAULT);

  /* nag_opt_miqp_mps_read (e04mxc).
     Reads MPS data file defining LP, QP, MILP or MIQP problem.
     Query call. */
  nag_opt_miqp_mps_read(fileid, 0, 0, 0, 0, 0, 0, mpslst, &n, &m, &nnz, &ncolh,
                        &nnzh, &lintvar, NULL, NULL, NULL, NULL, NULL, NULL,
                        NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
                        NAGERR_DEFAULT);

  /* nag_close_file (x04adc).
     Close file associated with given unit number. */
  nag_close_file(fileid, NAGERR_DEFAULT);

  maxm = m;
  maxn = n;
  maxnnz = nnz;
  maxnnzh = nnzh;
  maxncolh = ncolh;
  maxlintvar = (readints && lintvar > 0) ? lintvar : -1;

  if (!(irowa = NAG_ALLOC(maxnnz, Integer)) ||
      !(iccola = NAG_ALLOC(maxn + 1, Integer)) ||
      !(a = NAG_ALLOC(maxnnz, double)) ||
      !(bl = NAG_ALLOC(maxn + maxm, double)) ||
      !(bu = NAG_ALLOC(maxn + maxm, double)) ||
      !(irowh = NAG_ALLOC(maxnnzh, Integer)) ||
      !(iccolh = NAG_ALLOC(maxncolh + 1, Integer)) ||
      !(h = NAG_ALLOC(maxnnzh, double)) ||
      (maxlintvar > 0 && !(intvar = NAG_ALLOC(maxlintvar, Integer))) ||
      !(crname = NAG_ALLOC(maxn + maxm, e04mx_name)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  nag_open_file(fname, 0, &fileid, NAGERR_DEFAULT);

  /* Full call to the reader. */
  nag_opt_miqp_mps_read(fileid, maxn, maxm, maxnnz, maxncolh, maxnnzh,
                        maxlintvar, mpslst, &n, &m, &nnz, &ncolh, &nnzh,
                        &lintvar, &iobj, a, irowa, iccola, bl, bu, pnames,
                        &nname, crname, h, irowh, iccolh, &minmax,
                        (maxlintvar > 0) ? intvar : NULL, NAGERR_DEFAULT);

  nag_close_file(fileid, NAGERR_DEFAULT);

  /* Data has been read. Set up and run the solver.
     We have no explicit objective vector so set lenc = 0; the
     objective vector is stored in row iobj of a. */
  lenc = 0;
  objadd = 0.0;

  if (!(helast = NAG_ALLOC(n + m, Integer)) ||
      !(x = NAG_ALLOC(n + m, double)) ||
      !(pi = NAG_ALLOC(m, double)) ||
      !(rc = NAG_ALLOC(n + m, double)) ||
      !(hs = NAG_ALLOC(n + m, Integer)) ||
      !(iuser = NAG_ALLOC(ncolh + 1 + nnzh, Integer)) ||
      !(ruser = NAG_ALLOC(nnzh, double)) ||
      !(names = NAG_ALLOC(n + m, char *)))
    {
      printf("Allocation failure\n");
      exit_status = -3;
      goto END;
    }

  /* Transform char (*crname)[9] to char *names[] to be compatible with the
     solver. */
  for (i = 0; i < n + m; i++)
    names[i] = crname[i];

  for (i = 0; i < n + m; i++)
    helast[i] = 0;
  for (i = 0; i < n + m; i++)
    hs[i] = 0;
  for (i = 0; i < n + m; i++)
    x[i] = MIN(MAX(0.0, bl[i]), bu[i]);

  if (ncolh > 0)
    {
      /* Store the nonzeros of H in ruser for use by qphx. */
      for (i = 0; i < nnzh; i++)
        ruser[i] = h[i];
      /* Store iccolh and irowh in iuser for use by qphx. */
      for (i = 0; i < ncolh + 1; i++)
        iuser[i] = iccolh[i];
      for (i = ncolh + 1, j = 0; i < nnzh+ncolh+1; i++, j++)
        iuser[i] = irowh[j];
      comm.iuser = iuser;
      comm.user = ruser;
    }

  /* nag_opt_sparse_convex_qp_init (e04npc).
     Initialization function for nag_opt_sparse_convex_qp_solve (e04nqc). */
  nag_opt_sparse_convex_qp_init(&state, NAGERR_DEFAULT);

  /* Use nag_opt_sparse_convex_qp_option_set_string (e04nsc) to change
     the direction of optimization. Minimization is assumed by default. */
  if (minmax == 1)
    nag_opt_sparse_convex_qp_option_set_string("Maximize", &state,
                                               NAGERR_DEFAULT);
  else if (minmax == 0)
    nag_opt_sparse_convex_qp_option_set_string("Feasible Point", &state,
                                               NAGERR_DEFAULT);

  /* 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 and then call
     nag_opt_sparse_convex_qp_option_set_integer (e04ntc) to register that
     setting with the solver. */
  nag_open_file("", 2, &fileid, NAGERR_DEFAULT);
  nag_opt_sparse_convex_qp_option_set_integer("Print file", fileid, &state,
                                              NAGERR_DEFAULT);

  /* nag_opt_sparse_convex_qp_solve (e04nqc).
     LP or QP problem (suitable for sparse problems). */
  nag_opt_sparse_convex_qp_solve(Nag_Cold, qphx, m, n, nnz, nname, lenc,
                                 ncolh, iobj, objadd, pnames[0], a, irowa,
                                 iccola, bl, bu, NULL, (const char**)names,
                                 helast, hs, x, pi, rc, &ns, &ninf, &sinf,
                                 &obj, &state, &comm, &fail);

 END:
  NAG_FREE(a);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(h);
  NAG_FREE(pi);
  NAG_FREE(rc);
  NAG_FREE(ruser);
  NAG_FREE(x);
  NAG_FREE(helast);
  NAG_FREE(hs);
  NAG_FREE(iccola);
  NAG_FREE(iccolh);
  NAG_FREE(intvar);
  NAG_FREE(irowa);
  NAG_FREE(irowh);
  NAG_FREE(iuser);
  NAG_FREE(crname);
  NAG_FREE(names);
  return exit_status;
}

static void NAG_CALL qphx(Integer ncolh, const double x[], double hx[],
                          Integer nstate, Nag_Comm *comm)
{
  /* Function to compute H*x.
     Note: comm->iuser and comm->user contain the following data:
     comm->user[0:nnzh-1] = h[0:nnzh-1]
     comm->iuser[0:ncolh] = iccolh[0:ncolh]
     comm->iuser[ncolh+1:nnzh+ncolh] = irowh[0:nnzh-1] */

  Integer i, end, icol, idx, irow, start;

  for (i = 0; i < ncolh; i++)
    hx[i] = 0.0;
  for (icol = 0; icol < ncolh; icol++)
    {
      start = comm->iuser[icol];
      end = comm->iuser[icol + 1] - 1;
      for (idx = start-1; idx < end; idx++)
        {
          irow = comm->iuser[ncolh + 1 + idx] - 1;
          hx[irow] = hx[irow] + x[icol] * comm->user[idx];
          if (irow != icol)
            hx[icol] = hx[icol] + x[irow] * comm->user[idx];
        }
    }
}