/* nag_opt_handle_set_simplebounds (e04rhc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

/* 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 <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>
#include <nagf08.h>
#include <nagx04.h>
#include <math.h>

int
main(void)
{
  double const stol = 1e-5;
  double const time_limit = 120.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 == 1)
        {
          /* 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;
        }
      else
        {
          /* 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;
        }
    }

  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_gen_real_mat_print_comp (x04cbc).
       * Print the matrix. */
      nag_gen_real_mat_print_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_dgesvd (f08kbc).
       * Compute the singular values */
      nag_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;
}