/* nag_opt_handle_set_linconstr (e04rjc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 */

/* Read in LP/QP problem stored in a MPS file, formulated it
 * as a handle and pass it to the solver.
 */
#include <nag.h>
#include <stdio.h>

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

int main(int argc, char *argv[]) {
  char fname_default[] = "e04rjce.opt";
  Integer mpslst = 1;

  Integer exit_status = 0;
  Integer idlc, idx, idx_dest, inform, iobj, j, lintvar, m, maxlintvar, maxm,
      maxn, maxncolh, maxnnz, maxnnzh, minmax, n, ncolh, nname, nnz, nnzc, nnzh,
      nnzu, nnzua, nnzuc;
  char *fname = 0;
  char pnames[5][9] = {"", "", "", "", ""};
  double rinfo[32], stats[32];
  double *a = 0, *bl = 0, *bu = 0, *c = 0, *h = 0, *x = 0;
  Integer *iccola = 0, *iccolh = 0, *icola = 0, *icolh = 0, *idxc = 0,
          *irowa = 0, *irowh = 0;
  char(*crname)[9] = 0;
  void *handle = 0;
  /* Nag Types */
  Nag_FileID fileid;
  NagError fail;

  printf("nag_opt_handle_set_linconstr (e04rjc) Example Program Results\n\n");

  /* Use the first command line argument as the filename or
   * choose default filename stored in 'fname_default'. */
  if (argc > 1)
    fname = argv[1];
  else
    fname = fname_default;
  printf("Reading MPS file: %s\n", fname);
  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 or QP problem.
   * Query call to estimate the size of the problem. */
  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 = -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)) ||
      !(crname = NAG_ALLOC(maxn + maxm, e04mx_name)) ||
      !(x = NAG_ALLOC(maxn, double)) ||
      !(icolh = NAG_ALLOC(maxnnzh, Integer)) ||
      !(icola = NAG_ALLOC(maxnnz, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Reopen the same file. */
  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, NULL, NAGERR_DEFAULT);

  nag_file_close(fileid, NAGERR_DEFAULT);

  printf("MPS/QPS file read.\n");
  fflush(stdout);

  /* Data has been read. Set up the problem to the solver. */

  /* nag_opt_handle_init (e04rac).
   * Initialize an empty problem handle with n variables. */
  nag_opt_handle_init(&handle, n, NAGERR_DEFAULT);

  /* iccola[] was returned as 1-based, i.e., iccola[0]=1.
   * Change it to 0-based to simplify C operations. */
  for (j = 0; j <= n; j++)
    iccola[j]--;

  /* Move the linear objective from a to c. */
  if (iobj > 0) {
    /* Shift bounds. */
    for (j = n + iobj - 1; j < n + m - 1; j++) {
      bl[j] = bl[j + 1];
      bu[j] = bu[j + 1];
    }
    m--;

    /* Count how many nonzeros will be in c. */
    nnzc = 0;
    for (idx = 0; idx < nnz; idx++)
      if (irowa[idx] == iobj)
        nnzc++;

    if (!(idxc = NAG_ALLOC(nnzc, Integer)) || !(c = NAG_ALLOC(nnzc, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    /* Extract row iobj form column compressed matrix a. */
    idx = 0;
    idx_dest = 0;
    nnzc = 0;
    for (j = 0; j < n; j++) {
      for (; idx < iccola[j + 1]; idx++) {
        if (irowa[idx] < iobj) {
          a[idx_dest] = a[idx];
          irowa[idx_dest] = irowa[idx];
          idx_dest++;
        } else if (irowa[idx] == iobj) {
          idxc[nnzc] = j + 1;
          c[nnzc] = a[idx];
          nnzc++;
        } else {
          a[idx_dest] = a[idx];
          irowa[idx_dest] = irowa[idx] - 1;
          idx_dest++;
        }
      }
      iccola[j + 1] = idx_dest;
    }
    nnz = idx_dest;
  } else
    /* There is no linear part of the objective function. */
    nnzc = 0;

  /* Convert (decompress) iccola[] to icola[]. */
  for (j = 0; j < n; j++)
    for (idx = iccola[j]; idx < iccola[j + 1]; idx++)
      icola[idx] = j + 1;

  /* Add objective function to the problem formulation. */
  if (nnzh == 0)
    /* nag_opt_handle_set_quadobj (e04rfc).
     * Define the objective as a (sparse) linear function (no nonzeros in H).*/
    nag_opt_handle_set_quadobj(handle, nnzc, idxc, c, 0, NULL, NULL, NULL,
                               NAGERR_DEFAULT);
  else {
    /* The objective is a quadratic function.
     * Firstly, transform (decompress) iccolh[] to icolh[], both are 1-based.
     * Secondly, e04mxc() returned a lower triangle but here
     * the upper triangle is needed ==> swap rows and columns. */
    for (j = 0; j < ncolh; j++)
      for (idx = iccolh[j]; idx < iccolh[j + 1]; idx++)
        icolh[idx - 1] = j + 1;
    /* nag_opt_handle_set_quadobj (e04rfc).
     * Add the quadratic objective to the handle.*/
    nag_opt_handle_set_quadobj(handle, nnzc, idxc, c, nnzh, icolh, irowh, h,
                               NAGERR_DEFAULT);
  }

  /* nag_opt_handle_set_simplebounds (e04rhc).
   * Define bounds on the variables. */
  nag_opt_handle_set_simplebounds(handle, n, bl, bu, NAGERR_DEFAULT);

  idlc = 0;
  /* nag_opt_handle_set_linconstr (e04rjc).
   * Add a block of linear constraints to the problem formulation.
   * Bounds of the constraints are stored after bounds on the variables. */
  nag_opt_handle_set_linconstr(handle, m, bl + n, bu + n, nnz, irowa, icola, a,
                               &idlc, NAGERR_DEFAULT);

  printf("The problem was set-up\n");
  fflush(stdout);

  /* Set optional arguments of the solver by calling
   * nag_opt_handle_opt_set (e04zmc). */
  nag_opt_handle_opt_set(handle, "Print Options = No", NAGERR_DEFAULT);

  /* Set up a starting point and call the solver,
   * ignore Lagrangian multipliers U/UA. */
  for (j = 0; j < n; j++)
    x[j] = 0.0;
  nnzu = 0;
  nnzuc = 0;
  nnzua = 0;

  /* Call the solver nag_opt_handle_solve_pennon (e04svc). */
  INIT_FAIL(fail);
  nag_opt_handle_solve_pennon(handle, n, x, nnzu, NULL, nnzuc, NULL, nnzua,
                              NULL, rinfo, stats, &inform, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_handle_solve_pennon (e04svc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nOptimal solution:\n");
  for (j = 0; j < n; j++)
    printf("  %f\n", x[j]);
  fflush(stdout);

END:

  /* nag_opt_handle_free (e04rzc).
   * Destroy the problem handle and deallocate all the memory. */
  if (handle)
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);

  NAG_FREE(a);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(idxc);
  NAG_FREE(c);
  NAG_FREE(h);
  NAG_FREE(x);
  NAG_FREE(iccola);
  NAG_FREE(iccolh);
  NAG_FREE(icola);
  NAG_FREE(icolh);
  NAG_FREE(irowa);
  NAG_FREE(irowh);
  NAG_FREE(crname);
  return exit_status;
}