/* nag_opt_handle_print (e04ryc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

/* Demonstrate the life-cycle of a handle of a typical BMI-SDP problem
 * by printing the evolution of the HANDLE in certain stages to
 * the standard output.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>
#include <nagx04.h>

int main(void)
{
  Integer exit_status = 0;
  Integer dima, idblk, inform, nblk, nnzasum, nvar;
  double a[6], bl[2], bu[2], cvec[2], rinfo[32], stats[32], x[2];
  Integer icola[6], irowa[6], nnza[3], qi[1], qj[1];
  void *h = 0;
  /* Nag Types */
  Nag_FileID fileid;
  NagError fail;

  printf("nag_opt_handle_print (e04ryc) Example Program Results\n\n");
  fflush(stdout);

  /* Get Nag_FileID associated with the standard output by
   * nag_open_file (x04aac). */
  nag_open_file(NULL, 1, &fileid, NAGERR_DEFAULT);

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

  /* Anything can be defined at this phase. */
  printf("Freshly created handle\n");
  fflush(stdout);
  /* nag_opt_handle_print (e04ryc).
   * Print information currently stored in the handle. */
  nag_opt_handle_print(h, fileid, "Overview", NAGERR_DEFAULT);

  /* nag_opt_handle_set_linobj (e04rec).
   * Define linear objective (min y). */
  cvec[0] = 0.0;
  cvec[1] = 1.0;
  nag_opt_handle_set_linobj(h, nvar, cvec, NAGERR_DEFAULT);

  /* nag_opt_handle_set_simplebounds (e04rhc).
   * Add simple bounds (x>=0, -3<=y<=3) to the problem formulation. */
  bl[0] = 0.0;
  bu[0] = 1.0e20;
  bl[1] = -3.0;
  bu[1] = 3.0;
  nag_opt_handle_set_simplebounds(h, nvar, bl, bu, NAGERR_DEFAULT);

  /* The simple bounds and the objective are set and cannot be changed. */
  printf("\nHandle after definition of simple bounds and the objective\n");
  fflush(stdout);
  nag_opt_handle_print(h, fileid, "Overview,Objective,Simple Bounds",
                       NAGERR_DEFAULT);

  /* Definition of the first (linear) matrix constraint
   *    ( 1   x-1   y )
   *    (x-1  3/4   0 )  >= 0
   *    ( y    0   16 )
   * only upper triangles, thus we have matrices
   *       ( 1  -1  0 )        ( 0  1  0 )       ( 0  0  1 )
   * A0 = -(   3/4  0 ),  A1 = (    0  0 ), A2 = (    0  0 )
   *       (       16 )        (       0 )       (       0 )
   * Note: don't forget the minus at A0 term! */
  dima = 3;
  nnzasum = 6;
  nblk = 1;
  idblk = 0;
  /* A0 */
  irowa[0] = 1;
  icola[0] = 1;
  a[0] = -1.0;
  irowa[1] = 1;
  icola[1] = 2;
  a[1] = 1.0;
  irowa[2] = 2;
  icola[2] = 2;
  a[2] = -0.75;
  irowa[3] = 3;
  icola[3] = 3;
  a[3] = -16.0;
  nnza[0] = 4;
  /* A1 */
  irowa[4] = 1;
  icola[4] = 2;
  a[4] = 1.0;
  nnza[1] = 1;
  /* A2 */
  irowa[5] = 1;
  icola[5] = 3;
  a[5] = 1.0;
  nnza[2] = 1;
  /* nag_opt_handle_set_linmatineq (e04rnc).
   * Add a linear matrix inequality constraint to the problem formulation. */
  nag_opt_handle_set_linmatineq(h, nvar, dima, nnza, nnzasum, irowa, icola,
                                a, nblk, NULL, &idblk, NAGERR_DEFAULT);

  /* It is possible to add or extend existing matrix constraints. */
  printf("\nHandle after definition of the 1st matrix constraint\n");
  fflush(stdout);
  nag_opt_handle_print(h, fileid, "Overview,Matrix Constraints",
                       NAGERR_DEFAULT);

  /* Definition of the absolute term and linear part of BMI
   *    (  x   -xy )
   *    (-xy    1  ) >= 0
   * thus
   *       ( 0  0 )        ( 1  0 )
   * A0 = -(    1 ),  A1 = (    0 ), A2 = zero
   * Note: don't forget the minus at A0 term! */
  dima = 2;
  nnzasum = 2;
  nblk = 1;
  idblk = 0;
  /* A0 */
  irowa[0] = 2;
  icola[0] = 2;
  a[0] = -1.0;
  nnza[0] = 1;
  /* A1 */
  irowa[1] = 1;
  icola[1] = 1;
  a[1] = 1.0;
  nnza[1] = 1;
  /* A2 */
  nnza[2] = 0;
  /* nag_opt_handle_set_linmatineq (e04rnc).
   * Add another linear matrix inequality which will be extended to BMI. */
  nag_opt_handle_set_linmatineq(h, nvar, dima, nnza, nnzasum, irowa, icola,
                                a, nblk, NULL, &idblk, NAGERR_DEFAULT);

  /* It is possible to add or extend existing matrix constraints. */
  printf("\nHandle after partial definition of the 2nd matrix constraint\n");
  fflush(stdout);
  nag_opt_handle_print(h, fileid, "Matrix Constraints", NAGERR_DEFAULT);

  /* Extending current matrix constraint (with IDBLK) by bilinear term
   *       ( 0  -1 )
   * Q12 = ( 0   0 ). */
  dima = 2;
  nnzasum = 1;
  nnza[0] = 1;
  irowa[0] = 1;
  icola[0] = 2;
  a[0] = -1.0;
  qi[0] = 1;
  qj[0] = 2;
  /* nag_opt_handle_set_quadmatineq (e04rpc).
   * Add a bilinear matrix term. */
  nag_opt_handle_set_quadmatineq(h, 1, qi, qj, dima, nnza, nnzasum, irowa,
                                 icola, a, &idblk, NAGERR_DEFAULT);

  /* The problem is completely defined. */
  printf("\nHandle with the complete problem formulation\n");
  fflush(stdout);
  nag_opt_handle_print(h, fileid,
                       "Overview,Matrix Constraints,Multipliers Sizes",
                       NAGERR_DEFAULT);
  nag_opt_handle_print(h, fileid, "Matrix Constraints Detailed",
                       NAGERR_DEFAULT);

  /* Set optional arguments of the solver by calling
   * nag_opt_handle_opt_set (e04zmc). */
  nag_opt_handle_opt_set(h, "Print Options = No", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(h, "Initial X = Automatic", NAGERR_DEFAULT);

  /* Options can be printed even outside the solver. */
  printf("\n");
  fflush(stdout);
  nag_opt_handle_print(h, fileid, "Options", NAGERR_DEFAULT);

  /* Pass the handle to the solver
   * nag_opt_handle_solve_pennon (e04svc). */
  INIT_FAIL(fail);
  nag_opt_handle_solve_pennon(h, nvar, x, 0, NULL, 0, NULL, 0, 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;
  }

  /* After the solver finished. */
  printf("\nProblem solved\n");
  fflush(stdout);
  nag_opt_handle_print(h, fileid, "Overview", NAGERR_DEFAULT);

  /* Print the result. */
  printf("\nFinal objective function = %f\n", rinfo[0]);
  printf("Final x = [%f, %f].\n", x[0], x[1]);

  /* nag_opt_handle_free (e04rzc).
   * Destroy the problem handle and deallocate all the memory. */
  nag_opt_handle_free(&h, NAGERR_DEFAULT);

END:
  return exit_status;
}