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

/* Load a linear semidefinite programming problem from a sparse SDPA
 * file, formulate the problem via a handle, pass it to the solver
 * and print both primal and dual variables.
 */
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>
#include <nagx04.h>

int main(int argc, char *argv[])
{
  char fname_default[] = "e04rdce.opt";
  Integer filelst = 0;

  Integer exit_status = 0;
  Integer dima, idblk, idx, inform, j, k, maxnblk, maxnnz, maxnvar,
         nblk, nnz, nnzu, nnzuc, nnzua, nvar;
  char *fname = 0;
  char title[60];
  double rinfo[32], stats[32];
  double *a = 0, *cvec = 0, *ua = 0, *x = 0;
  Integer *blksizea = 0, *icola = 0, *irowa = 0, *nnza = 0;
  void *handle = 0;

  /* Nag Types */
  Nag_FileID infile;
  NagError fail;

  printf("nag_opt_sdp_read_sdpa (e04rdc) 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 SDPA file: %s\n", fname);
  fflush(stdout);

  /* nag_open_file (x04acc).
   * Open unit number for reading and associate unit with named file. */
  nag_open_file(fname, 0, &infile, NAGERR_DEFAULT);

  /* Go through the file for the first time to find the dimension
   * of the problem. Unless the file format is wrong, the function
   * should finish with fail.code = NE_INT_MAX (not enough space). */

  /* nag_opt_sdp_read_sdpa (e04rdc).
   * A reader of sparse SDPA data files for linear SDP problems. */
  INIT_FAIL(fail);
  nag_opt_sdp_read_sdpa(infile, 0, 0, 0, filelst, &nvar, &nblk, &nnz,
                        NULL, NULL, NULL, NULL, NULL, NULL, &fail);

  /* nag_close_file (x04adc).
   * Close file associated with given unit number. */
  nag_close_file(infile, NAGERR_DEFAULT);

  if (fail.code != NE_INT_MAX) {
    /* Possible problem with formatting, etc. Stop. */
    printf("Error from nag_opt_sdp_read_sdpa (e04rdc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  /* Allocate the right size of arrays for the data. */
  printf("Allocating space for the problem.\n");
  printf("nvar = %" NAG_IFMT "\n", nvar);
  printf("nblk = %" NAG_IFMT "\n", nblk);
  printf("nnz  = %" NAG_IFMT "\n", nnz);
  fflush(stdout);
  maxnvar = nvar;
  maxnblk = nblk;
  maxnnz = nnz;
  if (!(cvec = NAG_ALLOC(maxnvar, double)) ||
      !(nnza = NAG_ALLOC(maxnvar + 1, Integer)) ||
      !(irowa = NAG_ALLOC(maxnnz, Integer)) ||
      !(icola = NAG_ALLOC(maxnnz, Integer)) ||
      !(a = NAG_ALLOC(maxnnz, double)) ||
      !(blksizea = NAG_ALLOC(maxnblk, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Reopen the same file. */
  nag_open_file(fname, 0, &infile, NAGERR_DEFAULT);

  /* Read the problem data, there should be enough space this time.
   * nag_opt_sdp_read_sdpa (e04rdc).
   * A reader of sparse SDPA data files for linear SDP problems. */
  nag_opt_sdp_read_sdpa(infile, maxnvar, maxnblk, maxnnz, filelst, &nvar,
                        &nblk, &nnz, cvec, nnza, irowa, icola, a, blksizea,
                        NAGERR_DEFAULT);

  /* nag_close_file (x04adc).
   * Close file associated with given unit number. */
  nag_close_file(infile, NAGERR_DEFAULT);

  /* Problem was successfully decoded. */
  printf("Linear SDP problem was read, start formulating the problem\n");
  fflush(stdout);

  /* 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_linobj (e04rec).
   * Add a linear objective function to the formulation. */
  nag_opt_handle_set_linobj(handle, nvar, cvec, NAGERR_DEFAULT);

  dima = 0;
  for (k = 0; k < nblk; k++)
    dima += blksizea[k];
  idblk = 0;
  /* nag_opt_handle_set_linmatineq (e04rnc).
   * Add all linear matrix constraints to the formulation.*/
  nag_opt_handle_set_linmatineq(handle, nvar, dima, nnza, nnz, irowa, icola,
                                a, nblk, blksizea, &idblk, NAGERR_DEFAULT);

  printf("The problem formulation in a handle is completed.\n\n");
  fflush(stdout);

  /* nag_opt_handle_print (e04ryc).
   * Print overview of the handle. */
  nag_opt_handle_print(handle, 6, "Overview", NAGERR_DEFAULT);

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

  /* Compute memory needed for primal & dual variables. */

  /* There are no box constraints or linear constraints set
   * by e04rhc or by e04rjc, neither second order cone constraints.*/
  nnzu = 0;
  nnzuc = 0;

  /* Count size of the matrix multipliers, they are stored as packed
   * triangles respecting the block structure. */
  nnzua = 0;
  for (k = 0; k < nblk; k++)
    nnzua += blksizea[k] * (blksizea[k] + 1) / 2;

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

  /* Call the solver nag_opt_handle_solve_pennon (e04svc). */
  INIT_FAIL(fail);
  nag_opt_handle_solve_pennon(handle, nvar, x, nnzu, NULL, nnzuc, NULL,
                              nnzua, ua, 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 = 2;
    goto END;
  }

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

  /* Print packed lower triangles of the Lagrangian multipliers. */
  idx = 0;
  for (k = 0; k < nblk; k++) {
    sprintf(title, "Lagrangian multiplier for A_%" NAG_IFMT, k);
    nnz = blksizea[k] * (blksizea[k] + 1) / 2;
    /* nag_pack_real_mat_print (x04ccc).
     * Print real packed triangular matrix. */
    nag_pack_real_mat_print(Nag_ColMajor, Nag_Lower, Nag_NonUnitDiag,
                            blksizea[k], ua + idx, title, NULL,
                            NAGERR_DEFAULT);
    idx = idx + nnz;
  }

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(cvec);
  NAG_FREE(ua);
  NAG_FREE(x);
  NAG_FREE(blksizea);
  NAG_FREE(icola);
  NAG_FREE(irowa);
  NAG_FREE(nnza);
  return exit_status;
}