/* nag_opt_handle_set_linmatineq (e04rnc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

/* Compute E-optimal experiment design via semidefinite programming,
 * this can be done as follows
 *   max {lambda_min(A) | A = sum x_i*v_i*v_i^T, x_i>=0, sum x_i = 1}
 * where v_i are given vectors.
 */

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

int main(void) {
  const double big = 1e+20;
  const double tol = 0.00001;

  Integer exit_status = 0;
  Integer dima, i, idlc, idblk, idx, inform, j, k, m, nblk, nnzasum, nnzb, nnzu,
      nnzua, nnzuc, nvar, p;
  double blb[1], bub[1], c[1], rinfo[32], stats[32];
  Integer idxc[1];
  double *a = 0, *b = 0, *bl = 0, *bu = 0, *v = 0, *x = 0;
  Integer *icola = 0, *icolb = 0, *irowa = 0, *irowb = 0, *nnza = 0;
  void *handle = 0;
  /* Nag Types */
  NagError fail;

#define V(I, J) v[((I)-1) * p + (J)-1]

  printf("nag_opt_handle_set_linmatineq (e04rnc) Example Program Results\n\n");
  fflush(stdout);

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

  /* Read in the number of vectors and their size. */
  scanf("%" NAG_IFMT "%*[^\n]", &m);
  scanf("%" NAG_IFMT "%*[^\n]", &p);

  /* Allocate memory for the vectors and read them in. */
  if (!(v = NAG_ALLOC(m * p, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  for (i = 1; i <= m; i++)
    for (j = 1; j <= p; j++)
      scanf("%lf", &V(i, j));
  scanf("%*[^\n]");

  /* Variables of the problem will be x_1, ..., x_m (weights of the vectors)
   * and t (artificial variable for minimum eigenvalue). */
  nvar = m + 1;

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

  /* nag_opt_handle_set_quadobj (e04rfc).
   * Add the objective function to the handle: max t. */
  idxc[0] = m + 1;
  c[0] = 1.0;
  nag_opt_handle_set_quadobj(handle, 1, idxc, c, 0, NULL, NULL, NULL,
                             NAGERR_DEFAULT);

  if (!(bl = NAG_ALLOC(nvar, double)) || !(bu = NAG_ALLOC(nvar, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* nag_opt_handle_set_simplebounds (e04rhc).
   * Add simple bounds on variables to the problem formulation. */
  for (i = 0; i < m; i++)
    bl[i] = 0.0;
  bl[m] = -big;
  for (i = 0; i <= m; i++)
    bu[i] = big;
  nag_opt_handle_set_simplebounds(handle, nvar, bl, bu, NAGERR_DEFAULT);

  /* Create linear constraint: sum x_i = 1. */
  nnzb = m;
  if (!(irowb = NAG_ALLOC(nnzb, Integer)) ||
      !(icolb = NAG_ALLOC(nnzb, Integer)) || !(b = NAG_ALLOC(nnzb, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  for (i = 0; i < m; i++) {
    /* irowb, icolb use 1-based indices */
    irowb[i] = 1;
    icolb[i] = i + 1;
    b[i] = 1.0;
  }
  blb[0] = 1.0;
  bub[0] = 1.0;
  idlc = 0;

  /* nag_opt_handle_set_linconstr (e04rjc).
   * Add the linear constraint to the problem formulation. */
  nag_opt_handle_set_linconstr(handle, 1, blb, bub, nnzb, irowb, icolb, b,
                               &idlc, NAGERR_DEFAULT);

  /* Generate matrix constraint as:
   *   sum_{i=1}^m x_i*v_i*v_i^T - t*I >=0 */

  nblk = 1;
  idblk = 0;
  dima = p;
  /* Total number of nonzeros */
  nnzasum = p + m * (p + 1) * p / 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 empty */
  nnza[0] = 0;
  /* A_1, A_2, ..., A_m are v_i*v_i^T */
  idx = 0;
  for (k = 1; k <= m; k++) {
    nnza[k] = (p + 1) * p / 2;
    for (i = 1; i <= p; i++)
      for (j = i; j <= p; j++) {
        irowa[idx] = i;
        icola[idx] = j;
        a[idx] = V(k, i) * V(k, j);
        idx++;
      }
  }
  /* A_{m+1} is the -identity */
  nnza[m + 1] = p;
  for (i = 1; i <= p; i++) {
    irowa[idx] = i;
    icola[idx] = i;
    a[idx] = -1.0;
    idx++;
  }

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

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

  /* Pass the handle to the solver, we are not interested in
   * Lagrangian multipliers.
   * nag_opt_handle_solve_pennon (e04svc). */
  nnzu = 0;
  nnzuc = 0;
  nnzua = 0;
  INIT_FAIL(fail);
  nag_opt_handle_solve_pennon(handle, nvar, 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;
  }

  /* Print results */
  printf("\n  Weight        Row of design matrix");
  for (i = 1; i <= m; i++)
    if (x[i - 1] > tol) {
      printf("\n %7.2f      ", x[i - 1]);
      for (j = 1; j <= p; j++)
        printf("%7.2f ", V(i, j));
    }
  printf("\n only those rows with a weight > %7.1e are shown\n", tol);

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(b);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(v);
  NAG_FREE(x);
  NAG_FREE(icola);
  NAG_FREE(icolb);
  NAG_FREE(irowa);
  NAG_FREE(irowb);
  NAG_FREE(nnza);
  return exit_status;
}