/* nag_opt_handle_solve_pennon (e04svc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.3, 2024.
 */

/* Compute Lovasz theta number of the given graph G on the input
 * via semidefinite programming as
 *   min {lambda_max(H) | H is nv x nv symmetric matrix where
 *                        h_ij=1 if ij is not an edge or if i==j}
 */

#include <nag.h>
#include <stdio.h>

int main(void) {
  Integer exit_status = 0;
  Integer i, idblk, idx, inform, ivalue, j, maxe, ne, nnzasum, nnzu, nnzua,
      nnzuc, nv, nvar;
  double rvalue, c[1], rinfo[32], stats[32], *a = 0, *x = 0;
  Integer idxc[1], *icola = 0, *irowa = 0, *nnza = 0, *va = 0, *vb = 0;
  char cvalue[41];
  void *handle = 0;
  /* Nag Types */
  NagError fail;
  Nag_VariableType optype;

  printf("nag_opt_handle_solve_pennon (e04svc) Example Program Results\n\n");
  fflush(stdout);

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

  /* Read in the number of vertices and edges of the graph. */
  scanf("%" NAG_IFMT "%*[^\n]", &nv);
  scanf("%" NAG_IFMT "%*[^\n]", &ne);

  if (!(va = NAG_ALLOC(ne, Integer)) || !(vb = NAG_ALLOC(ne, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in edges of the graph, accept only 1<=i<j<=nv. */
  maxe = ne;
  ne = 0;
  for (idx = 0; idx < maxe; idx++) {
    scanf("%" NAG_IFMT " %" NAG_IFMT "%*[^\n]", &i, &j);
    if (i >= 1 && i < j && j <= nv) {
      va[ne] = i;
      vb[ne] = j;
      ne++;
    }
  }

  /* Number of variables (same as edges in the graph plus one). */
  nvar = ne + 1;

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

  idxc[0] = 1;
  c[0] = 1.0;
  /* nag_opt_handle_set_quadobj (e04rfc).
   * Add the quadratic objective to the handle. */
  nag_opt_handle_set_quadobj(handle, 1, idxc, c, 0, NULL, NULL, NULL,
                             NAGERR_DEFAULT);

  /* Generate matrix constraint as:
   *   sum_{ij is edge in G} x_ij*E_ij + t*I - J >=0
   * where J is the all-ones matrix. Its dimension is the same
   * as the number of vertices. */

  /* Total number of nonzeros: */
  nnzasum = ne + nv + (nv + 1) * nv / 2;
  if (!(nnza = NAG_ALLOC(nvar + 1, Integer)) ||
      !(irowa = NAG_ALLOC(nnzasum, Integer)) ||
      !(icola = NAG_ALLOC(nnzasum, Integer)) ||
      !(a = NAG_ALLOC(nnzasum, double)) || !(x = NAG_ALLOC(nvar, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* A_0 is all ones matrix. */
  idx = 0;
  nnza[0] = (nv + 1) * nv / 2;
  for (i = 1; i <= nv; i++)
    for (j = i; j <= nv; j++) {
      irowa[idx] = i;
      icola[idx] = j;
      a[idx] = 1.0;
      idx++;
    }
  /* A_1 is the identity. */
  nnza[1] = nv;
  for (i = 1; i <= nv; i++) {
    irowa[idx] = i;
    icola[idx] = i;
    a[idx] = 1.0;
    idx++;
  }
  /* A_2, A_3, ..., A_{ne+1} match the E_ij matrices. */
  for (i = 0; i < ne; i++) {
    nnza[2 + i] = 1;
    irowa[idx] = va[i];
    icola[idx] = vb[i];
    a[idx] = 1.0;
    idx++;
  }
  idblk = 0;

  /* nag_opt_handle_set_linconstr (e04rnc).
   * Add the linear matrix constraint to the problem formulation. */
  nag_opt_handle_set_linmatineq(handle, nvar, nv, nnza, nnzasum, irowa, icola,
                                a, 1, NULL, &idblk, NAGERR_DEFAULT);

  /* nag_opt_handle_opt_set (e04zmc).
   * Set optional arguments of the solver */
  nag_opt_handle_opt_set(handle, "Initial X = Automatic", NAGERR_DEFAULT);

  /* Pass the handle to the solver, we are not interested in
   * Lagrangian multipliers. */
  nnzu = 0;
  nnzuc = 0;
  nnzua = 0;
  INIT_FAIL(fail);
  /* nag_opt_handle_solve_pennon (e04svc). */
  nag_opt_handle_solve_pennon(handle, nvar, x, nnzu, NULL, nnzuc, NULL, nnzua,
                              NULL, rinfo, stats, &inform, &fail);

  if (fail.code == NE_NOERROR || fail.code == NW_NOT_CONVERGED) {
    /* Retrieve some of the settings by calling
     * nag_opt_handle_opt_get (e04znc). */
    nag_opt_handle_opt_get(handle, "Hessian Density", &ivalue, &rvalue, cvalue,
                           40, &optype, NAGERR_DEFAULT);
    printf("The solver chose to use %s hessian", cvalue);
    nag_opt_handle_opt_get(handle, "Linesearch Mode", &ivalue, &rvalue, cvalue,
                           40, &optype, NAGERR_DEFAULT);
    printf(" and %s as linesearch.\n", cvalue);
    printf("Lovasz theta number of the given graph is %7.2f.\n", rinfo[0]);
  } else {
    printf("Error from nag_opt_handle_solve_pennon (e04svc).\n%s\n",
           fail.message);
    exit_status = 1;
  }

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(x);
  NAG_FREE(icola);
  NAG_FREE(irowa);
  NAG_FREE(nnza);
  NAG_FREE(va);
  NAG_FREE(vb);
  return exit_status;
}