NAG Library Manual, Mark 30
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_opt_handle_set_simplebounds (e04rhc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.0, 2024.
 */

/* Matrix completion problem (rank minimization) solved approximatelly
 *    by SDP via nuclear norm minimization formulated as:
 *       min    trace(X1) + trace(X2)
 *       s.t.   [ X1, Y; Y', X2 ] >=0
 *              0 <= Y_ij <= 1
 */

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

int main(void) {
  double const stol = 1e-5;
  double const time_limit = 180.0;

  Integer exit_status = 0;
  Integer dima, i, idblk, idx, idxend, idxobj, idxx, inform, j, lwork, m, maxs,
      n, nblk, nnz, nnzasum, nnzc, nnzu, nnzua, nnzuc, nvar, rank;
  double rdummy[1], rinfo[32], stats[32];
  double *a = 0, *bl = 0, *bu = 0, *c = 0, *s = 0, *work = 0, *x = 0, *y = 0;
  Integer *icola = 0, *idxc = 0, *irowa = 0, *nnza = 0;
  void *handle = 0;
  /* Nag Types */
  NagError fail;

#define Y(I, J) y[(J - 1) * m + I - 1]

  printf("nag_opt_handle_set_simplebounds (e04rhc) Example Program Results"
         "\n\n");
  fflush(stdout);

  /* Skip heading in data file. */
  scanf("%*[^\n]");
  /* Read in the problem size and ignore the rest of the line. */
  scanf("%" NAG_IFMT " %" NAG_IFMT "%*[^\n]", &m, &n);

  /* Allocate memory for matrix Y and read it in. */
  if (!(y = NAG_ALLOC(m * n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  for (i = 1; i <= m; i++)
    for (j = 1; j <= n; j++)
      scanf("%lf", &Y(i, j));
  scanf("%*[^\n]");

  /* Count the number of specified elements (i.e., nonnegative) */
  nnz = 0;
  for (i = 1; i <= m; i++)
    for (j = 1; j <= n; j++)
      if (Y(i, j) >= 0.0)
        nnz++;

  /* There are as many variables as missing entries in the Y matrix
   *  plus two full symmetric matrices m x m and n x n. */
  nvar = m * (m + 1) / 2 + n * (n + 1) / 2 + m * n - nnz;
  if (!(x = NAG_ALLOC(nvar, double)) || !(bl = NAG_ALLOC(nvar, double)) ||
      !(bu = NAG_ALLOC(nvar, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

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

  /* Create bounds for the missing entries in Y matrix to be between 0 and 1 */
  idxend = m * (m + 1) / 2 + n * (n + 1) / 2;
  for (idx = 0; idx < idxend; idx++) {
    bl[idx] = -1e+20;
    bu[idx] = 1e+20;
  }
  for (; idx < nvar; idx++) {
    bl[idx] = 0.0;
    bu[idx] = 1.0;
  }

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

  /* Allocate space for the objective - minimize trace of the matrix
   * constraint. There is no quadratic part in the objective. */
  nnzc = m + n;
  if (!(idxc = NAG_ALLOC(nnzc, Integer)) || !(c = NAG_ALLOC(nnzc, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Construct linear matrix inequality to request that
   *   [ X1, Y; Y', X2] is positive semidefinite. */

  /* How many nonzeros do we need? As many as number of variables
   * and the number of specified elements together. */
  nnzasum = m * (m + 1) / 2 + n * (n + 1) / 2 + m * n;
  if (!(nnza = NAG_ALLOC(nvar + 1, Integer)) ||
      !(irowa = NAG_ALLOC(nnzasum, Integer)) ||
      !(icola = NAG_ALLOC(nnzasum, Integer)) ||
      !(a = NAG_ALLOC(nnzasum, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  nnza[0] = nnz;
  for (i = 1; i <= nvar; i++)
    nnza[i] = 1;

  /* Copy Y to the upper block of A_0 with different sign (because -A_0)
   * (upper triangle) */
  idx = 0;
  for (i = 1; i <= m; i++)
    for (j = 1; j <= n; j++)
      if (Y(i, j) >= 0.0) {
        irowa[idx] = i;
        icola[idx] = m + j;
        a[idx] = -Y(i, j);
        idx++;
      }
  /* One matrix for each variable, A_i has just one nonzero - it binds
   * x_i with its position in the matrix constraint. Set also the objective. */
  /* 1,1 - block, X1 matrix (mxm) */
  idxobj = 0;
  idxx = 1;
  for (i = 1; i <= m; i++) {
    /* the next element is diagonal ==> part of the objective as a trace() */
    idxc[idxobj] = idxx;
    c[idxobj] = 1.0;
    idxobj++;
    for (j = i; j <= m; j++) {
      irowa[idx] = i;
      icola[idx] = j;
      a[idx] = 1.0;
      idx++;
      idxx++;
    }
  }
  /* 2,2 - block, X2 matrix (nxn) */
  for (i = 1; i <= n; i++) {
    /* the next element is diagonal ==> part of the objective as a trace() */
    idxc[idxobj] = idxx;
    c[idxobj] = 1.0;
    idxobj++;
    for (j = i; j <= n; j++) {
      irowa[idx] = m + i;
      icola[idx] = m + j;
      a[idx] = 1.0;
      idx++;
      idxx++;
    }
  }
  /* 1,2 - block, missing element in Y we are after */
  for (i = 1; i <= m; i++)
    for (j = 1; j <= n; j++)
      if (Y(i, j) < 0.0) {
        irowa[idx] = i;
        icola[idx] = m + j;
        a[idx] = 1.0;
        idx++;
      }

  /* nag_opt_handle_set_quadobj (e04rfc).
   * Add the sparse linear objective to the handle.*/
  nag_opt_handle_set_quadobj(handle, nnzc, idxc, c, 0, NULL, NULL, NULL,
                             NAGERR_DEFAULT);

  /* Just one matrix inequality of the dimension of the extended matrix. */
  nblk = 1;
  dima = m + n;
  idblk = 0;

  /* nag_opt_handle_set_linconstr (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);

  /* nag_opt_handle_opt_set (e04zmc).
   * Set optional arguments of the solver:
   * Completely turn off printing, allow timing and
   * turn on the monitor mode to stop every iteration. */
  nag_opt_handle_opt_set(handle, "Print File = -1", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Stats Time = Yes", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Monitor Frequency = 1", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Initial X = Automatic", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Dimacs = Check", NAGERR_DEFAULT);

  /* Pass the handle to the solver, we are not interested in
   * Lagrangian multipliers. */
  nnzu = 0;
  nnzuc = 0;
  nnzua = 0;
  while (1) {
    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 (inform == 0 || fail.code != NE_NOERROR) {
      /* Final exit, solver finished. */
      printf("Finished at iteration %2" NAG_IFMT
             ": objective %7.2f, avg.error %9.2e\n",
             (Integer)stats[0], rinfo[0],
             (rinfo[1] + rinfo[2] + rinfo[3]) / 3.0);
      fflush(stdout);
      break;
    } else {
      /* Monitor stop */
      printf("Monitor at iteration  %2" NAG_IFMT
             ": objective %7.2f, avg.error %9.2e\n",
             (Integer)stats[0], rinfo[0],
             (rinfo[1] + rinfo[2] + rinfo[3]) / 3.0);
      fflush(stdout);

      /* Check time limit and possibly stop the solver. */
      if (stats[7] > time_limit)
        inform = 0;
    }
  }

  if (fail.code == NE_NOERROR || fail.code == NW_NOT_CONVERGED) {
    /* Successful run, fill the missing elements in the matrix Y. */
    idx = m * (m + 1) / 2 + n * (n + 1) / 2;
    for (i = 1; i <= m; i++)
      for (j = 1; j <= n; j++)
        if (Y(i, j) < 0.0)
          Y(i, j) = x[idx++];

    /* nag_file_print_matrix_real_gen_comp (x04cbc).
     * Print the matrix. */
    nag_file_print_matrix_real_gen_comp(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, m, n, y, m, "%7.1f",
        "Completed Matrix", Nag_IntegerLabels, NULL, Nag_IntegerLabels, NULL,
        80, 0, NULL, NAGERR_DEFAULT);

    /* Compute rank of the matrix via SVD, use the fact that the order
     * of the singular values is descending. */
    lwork = 20 * (m > n ? m : n);
    maxs = m < n ? m : n;
    if (!(s = NAG_ALLOC(maxs, double)) ||
        !(work = NAG_ALLOC((lwork), double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    /* nag_lapackeig_dgesvd (f08kbc).
     * Compute the singular values */
    nag_lapackeig_dgesvd(Nag_ColMajor, Nag_NotU, Nag_NotVT, m, n, y, m, s,
                         rdummy, 1, rdummy, 1, work, NAGERR_DEFAULT);
    for (rank = 0; rank < maxs; rank++)
      if (s[rank] <= stol)
        break;
    printf("Rank is %" NAG_IFMT "\n", rank);
  } else if (fail.code == NE_USER_STOP) {
    printf("The given time limit was reached, run aborted.\n");
  } 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(bl);
  NAG_FREE(bu);
  NAG_FREE(c);
  NAG_FREE(s);
  NAG_FREE(work);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(icola);
  NAG_FREE(irowa);
  NAG_FREE(nnza);
  NAG_FREE(idxc);
  return exit_status;
}