/* nag_opt_miqp_mps_read (e04mxc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.3, 2024.
 */
#include <nag.h>
#include <stdio.h>
#include <string.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, verbose_output;
  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_file_open (x04acc).
     Open unit number for reading and associate unit with named file. */
  nag_file_open(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_file_close (x04adc).
     Close file associated with given unit number. */
  nag_file_close(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_file_open(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_file_close(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_qpconvex2_sparse_init (e04npc).
     Initialization function for nag_opt_qpconvex2_sparse_solve (e04nqc). */
  nag_opt_qpconvex2_sparse_init(&state, NAGERR_DEFAULT);

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

  /* Set this to 1 to cause intermediate progress output to be produced. */
  verbose_output = 0;

  if (verbose_output == 1) {
    /* By default nag_opt_qpconvex2_sparse_solve (e04nqc) does not print
     * monitoring information. Call nag_file_open (x04acc) to set the print
     * file fileid and then call
     * nag_opt_qpconvex2_sparse_option_integer_set (e04ntc) to register that
     * setting with the solver.
     */
    nag_file_open("", 2, &fileid, NAGERR_DEFAULT);
    nag_opt_qpconvex2_sparse_option_integer_set("Print file", fileid, &state,
                                                NAGERR_DEFAULT);
  } else {
    printf("\nProblem contains %3" NAG_IFMT " variables and %3" NAG_IFMT
           " linear constraints\n",
           n, m);
  }
  /* nag_opt_qpconvex2_sparse_solve (e04nqc).
     LP or QP problem (suitable for sparse problems). */
  nag_opt_qpconvex2_sparse_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);

  if (verbose_output == 0) {
    printf("\n");
    printf("Final objective value = %12.3e\n", obj);
    printf("Optimal x = ");

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

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];
    }
  }
}