/* nag_opt_miqp_mps_write (e04mwc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2017.
 */

#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 e04mw_name[9];

int main(void) {
  /* Scalars and arrays */

  /* Scalars */
  Integer exit_status = 0;
  double obj, objadd, sinf;
  Integer i, icol, iobj, j, jcol, lenc, m, n, ncolh, ne, ninf, nnzh, minmax,
    lintvar;
  Integer nname, ns;

  /* output file */
  char fname[] = "e04mwce.mps";

  /* Arrays */
  char prob[9];
  char **names;
  char pnames[5][9] = {"", "", "", "", ""};
  char (*crname)[9] = 0;
  double *acol = 0, *bl = 0, *bu = 0, *c = 0, *pi = 0, *rc = 0, *x = 0, *h = 0,
    *ruser = 0;
  Integer *helast = 0, *hs = 0, *inda = 0, *loca = 0, *iccolh = 0, *irowh = 0,
    *iuser = 0;

  /* Nag Types */
  Nag_E04State state;
  Nag_Start start;
  NagError fail;
  Nag_Comm comm;
  Nag_FileID fileid;

  INIT_FAIL(fail);

  printf("nag_opt_miqp_mps_write (e04mwc) Example Program Results\n\n");
  fflush(stdout);

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

  /* Read ne, iobj, ncolh, nnzh and nname from data file. */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &n, &m);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT
        "%*[^\n] ", &ne, &iobj, &ncolh, &nnzh, &nname);

  if (n >= 1 && m >= 1) {
    /* Allocate memory */
    if (!(names = NAG_ALLOC(n + m, char *)) ||
        !(acol = NAG_ALLOC(ne, double)) ||
        !(bl = NAG_ALLOC(m + n, double)) ||
        !(bu = NAG_ALLOC(m + n, double)) ||
        !(c = NAG_ALLOC(1, double)) ||
        !(pi = NAG_ALLOC(m, double)) ||
        !(rc = NAG_ALLOC(n + m, double)) ||
        !(x = NAG_ALLOC(n + m, double)) ||
        !(helast = NAG_ALLOC(n + m, Integer)) ||
        !(hs = NAG_ALLOC(n + m, Integer)) ||
        !(inda = NAG_ALLOC(ne, Integer)) ||
        !(loca = NAG_ALLOC(n + 1, Integer)) ||
        !(iccolh = NAG_ALLOC(ncolh+1, Integer)) ||
        !(irowh = NAG_ALLOC(nnzh, Integer)) ||
        !(h = NAG_ALLOC(nnzh, double)) ||
        !(crname = NAG_ALLOC(n + m, e04mw_name)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    printf("%s", "Either m or n invalid\n");
    exit_status = 1;
    return exit_status;
  }

  /* Read names fron file */
  for (i = 0; i < nname; ++i) {
    if (!(names[i] = NAG_ALLOC(9, char)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    scanf(" ' %8s '", names[i]);
  }
  scanf("%*[^\n] ");

  /* Read the matrix acol from data file. Set up loca. */
  jcol = 0;
  loca[jcol] = 1;
  for (i = 0; i < ne; ++i) {
    /* Element (inda[i], icol) is stored in acol[i]. */
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &acol[i], &inda[i], &icol);
    if (icol <= jcol) {
      /* Elements not ordered by increasing column index. */
      printf("%s%5" NAG_IFMT "%s%5" NAG_IFMT "%s%s\n", "Element in column",
             icol, " found after element in column", jcol + 1, ". Problem",
             " abandoned.");
    }
    while (jcol + 1 < icol) {
      /* Fill in loca[] for missing columns */
      jcol++;
      loca[jcol] = i + 1;
    }
  }
  while (jcol < n) {
    /* Fill in loca[] for remaining columns */
    jcol++;
    loca[jcol] = i + 1;
  }

  /* Read the matrix h, from data file. Set up iccolh */
  jcol = 0;
  iccolh[jcol] = 1;
  for (i = 0; i < nnzh; ++i) {
    /* Element (irowh[i], icol) is stored in h[i]. */
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &h[i], &irowh[i],
          &icol);
    if (icol <= jcol) {
      /* Elements not ordered by increasing column index. */
      printf("%s%5" NAG_IFMT "%s%5" NAG_IFMT "%s%s\n", "Element in column",
             icol, " found after element in column", jcol, ". Problem",
             " abandoned.");
    }
    while (jcol + 1 < icol) {
      /* Fill in iccolh[] for missing columns */
      jcol++;
      iccolh[jcol] = i + 1;
    }
  }
  while (jcol < ncolh) {
    /* Fill in iccolh[] for remaining columns */
    jcol++;
    iccolh[jcol] = i + 1;
  }

  /* Read lower and upper bound */
  for (i = 0; i < n + m; ++i) {
    scanf("%lf", &bl[i]);
  }
  scanf("%*[^\n] ");

  for (i = 0; i < n + m; ++i) {
    scanf("%lf", &bu[i]);
  }
  scanf("%*[^\n] ");

  /* Set cold start and hs*/
  start = Nag_Cold;
  for (i = 0; i < n; ++i) {
    hs[i] = 0.0;
  }

  /* Read the initial point x_0 */
  for (i = 0; i < n; ++i) {
    scanf("%lf", &x[i]);
  }
  scanf("%*[^\n] ");

  /* nag_opt_sparse_convex_qp_init (e04npc).
     Initialization function for nag_opt_sparse_convex_qp_solve (e04nqc) */
  nag_opt_sparse_convex_qp_init(&state, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Initialization of "
           "nag_opt_sparse_convex_qp_solve (e04nqc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

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

     nag_open_file (x04acc).
     Open unit number for reading, writing or appending, and
     associate unit with named file */
  nag_open_file("", 2, &fileid, NAGERR_DEFAULT);

  /* nag_opt_sparse_convex_qp_option_set_integer (e04ntc).
     Set a single option for nag_opt_sparse_convex_qp_solve (e04nqc)
     from an integer argument */
  nag_opt_sparse_convex_qp_option_set_integer("Print file", fileid, &state,
                                              &fail);
  if (fail.code != NE_NOERROR) {
    exit_status = 1;
    goto END;
  }

  /* We have no explicit objective vector so set lenc = 0; the
    objective vector is stored in row iobj of acol. */
  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 = 0; i < n + m; ++i) {
    helast[i] = 0;
  }

  /* Illustrate how to pass information to the user-supplied
     function qphx via the comm structure */
  comm.p = 0;

  if (!(iuser = NAG_ALLOC(ncolh + 1 + nnzh, Integer)) ||
    !(ruser = NAG_ALLOC(nnzh, double)))
  {
    printf("Allocation failure\n");
    exit_status = -3;
    goto END;
  }

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

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

  if (fail.code != NE_NOERROR) {
    printf("nag_opt_sparse_convex_qp_solve (e04nqc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print objective function and optimal x* */
  printf("Final objective value = %12.3e\n", obj);
  printf("Optimal X = ");

  for (i = 0; i < n; ++i) {
    printf("%9.2f%s", x[i], i % 7 == 6 || i == n - 1 ? "\n" : " ");
  }

  /* Set the rest of input for e04mwc. Due to iobj > 0, lenc = 0 and c and idxc
     will not be referenced (are optional arguments)
     c - optional
     idxc - optional
     intvar - optional */

  /* pnames */
  strcpy(pnames[0], "USRPNAME");
  strcpy(pnames[1], "OBJ.....");
  strcpy(pnames[2], "RHS.....");
  strcpy(pnames[3], "RANGE...");
  strcpy(pnames[4], "BOUND...");
  lintvar = 0;
  minmax = -1;

  /* transform char *names[] to (*crname)[9] */
  for (i = 0; i < n + m; i++)
    strcpy(crname[i], names[i]);

  /* open data file for writing */
  nag_open_file(fname, 1, &fileid, NAGERR_DEFAULT);

  /* nag_opt_miqp_mps_write (e04mwc).
    writes data for sparse lP, MILP, QP or MIQP problems to a
    file in MPS format. */
  nag_opt_miqp_mps_write(fileid, n, m, lenc, ne, ncolh, nnzh, lintvar,NULL,NULL,
          iobj, acol, inda, loca, bl, bu, pnames, crname, h, irowh,
          iccolh, minmax, NULL, &fail);

  if (fail.code != NE_NOERROR) {
    printf("nag_opt_miqp_mps_write (e04mwc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print name of the outfile */
  printf("\nMPS file was written: %s\n", fname);

  /* Close data file */
  nag_close_file(fileid, NAGERR_DEFAULT);

END:

  for (i = 0; i < n + m; i++) {
    NAG_FREE(names[i]);
  }
  NAG_FREE(acol);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(h);
  NAG_FREE(c);
  NAG_FREE(pi);
  NAG_FREE(rc);
  NAG_FREE(x);
  NAG_FREE(helast);
  NAG_FREE(hs);
  NAG_FREE(loca);
  NAG_FREE(iccolh);
  NAG_FREE(inda);
  NAG_FREE(irowh);
  NAG_FREE(iuser);
  NAG_FREE(ruser);
  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];
    }
  }
}