/* nag_opt_handle_solve_lp_ipm (e04mtc) Example Program.
 *
 * Copyright 2019 Numerical Algorithms Group.
 *
 * Mark 27.0, 2019.
 */

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

#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL monit(void *handle, const double rinfo[],
                           const double stats[], Nag_Comm *comm,
                           Integer *inform);
#ifdef __cplusplus
}
#endif

int main(void){

  Integer nclin, nvar, nnza, nnzc, nnzu, exit_status, i;
  Integer idlc;
  Integer *irowa = 0, *icola = 0;
  Integer iuser[2];
  double *cvec = 0, *a = 0, *bla = 0, *bua = 0, *xl = 0, *xu = 0,
    *x = 0, *u = 0;
  double rinfo[100], stats[100];
  void *handle = 0;
  /* Nag Types */
  Nag_Comm comm;

  exit_status = 0;

  printf("nag_opt_handle_solve_lp_ipm (e04mtc) Example Program Results\n\n");
  fflush(stdout);

  /* Read the data file and allocate memory */
  scanf(" %*[^\n]"); /* Skip heading in data file */
  scanf("%" NAG_IFMT" %" NAG_IFMT" %" NAG_IFMT" %" NAG_IFMT" %*[^\n]",&nclin,&nvar,&nnza,&nnzc);
  /* Allocate memory */
  nnzu = 2*nvar + 2*nclin;
  if (!(irowa = NAG_ALLOC(nnza, Integer))  ||
      !(icola = NAG_ALLOC(nnza, Integer))  ||
      !(cvec = NAG_ALLOC(nnzc,double))     ||
      !(a = NAG_ALLOC(nnza,double))        ||
      !(bla = NAG_ALLOC(nclin,double))     ||
      !(bua = NAG_ALLOC(nclin,double))     ||
      !(xl = NAG_ALLOC(nvar,double))       ||
      !(xu = NAG_ALLOC(nvar,double))       ||
      !(x = NAG_ALLOC(nvar,double))        ||
      !(u = NAG_ALLOC(nnzu,double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  for (i=0; i< nvar; i++){
    x[i] = 0.0;
  }

  /* Read objective */
  for (i=0; i<nnzc; i++){
    scanf("%lf",&cvec[i]);
  }
  scanf("%*[^\n]");
  /* Read constraint matrix row indices */
  for (i=0; i<nnza; i++){
    scanf("%" NAG_IFMT,&irowa[i]);
  }
  scanf("%*[^\n]");
  /* Read constraint matrix col indices */
  for (i=0; i<nnza; i++){
    scanf("%" NAG_IFMT,&icola[i]);
  }
  scanf("%*[^\n]");
  /* Read constraint matrix values */
  for (i=0; i<nnza; i++){
    scanf("%lf",&a[i]);
  }
  scanf("%*[^\n]");
  /* Read linear constraints lower bounds */
  for (i=0; i<nclin; i++){
    scanf("%lf ",&bla[i]);
  }
  scanf("%*[^\n]");
  /* Read linear constraints upper bounds */
  for (i=0; i<nclin; i++){
    scanf("%lf ",&bua[i]);
  }
  scanf("%*[^\n]");
  /* Read variables lower bounds */
  for (i=0; i<nvar; i++){
    scanf("%lf ",&xl[i]);
  }
  scanf("%*[^\n]");
  /* Read variables upper bounds */
  for (i=0; i<nvar; i++){
    scanf("%lf ",&xu[i]);
  }
  scanf("%*[^\n]");

  /* Create the problem handle */
  /* 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)
   * Define a linear objective */
  nag_opt_handle_set_linobj(handle,nvar,cvec,NAGERR_DEFAULT);

  /* nag_opt_handle_set_simplebounds (e04rhc)
   * Define bounds on the variables */
  nag_opt_handle_set_simplebounds(handle,nvar,xl,xu,NAGERR_DEFAULT);

  /* nag_opt_handle_set_linconstr (e04rjc)
   * Define linear constraints */
  idlc = 0;
  nag_opt_handle_set_linconstr(handle,nclin,bla,bua,nnza,irowa,
                               icola,a,&idlc,NAGERR_DEFAULT);

  /* nag_opt_handle_opt_set (e04zmc)
   * Require a high accuracy solution */
  nag_opt_handle_opt_set(handle, "LPIPM Stop Tolerance = 1.0e-10",
                         NAGERR_DEFAULT);
  /* Require printing of the solution at the end of the solve */
  nag_opt_handle_opt_set(handle, "Print Solution = Yes",
                         NAGERR_DEFAULT);
  /* Use a constant number of centrality correctors steps */
  nag_opt_handle_opt_set(handle, "LPIPM Centrality Correctors = -6",
                         NAGERR_DEFAULT);
  /* Turn on monitoring */
  nag_opt_handle_opt_set(handle, "LPIPM Monitor Frequency = 1",
                         NAGERR_DEFAULT);

  comm.iuser = iuser;
  iuser[0] = 1;

  /* nag_opt_handle_solve_lp_ipm (e04mtc)
   * Call LP interior point solver with the default (primal-dual) algorithm */
  printf("\n++++++++++ Use the Primal-Dual algorithm ++++++++++\n");
  fflush(stdout);
  nag_opt_handle_solve_lp_ipm(handle,nvar,x,nnzu,u,rinfo,stats,monit,
                              &comm,NAGERR_DEFAULT);

  iuser[0] = 2;
  /* Solve the same problem with the self-dual algorithm */
  printf("\n++++++++++ Use the Self-Dual algorithm ++++++++++\n");
  fflush(stdout);
  nag_opt_handle_opt_set(handle, "LPIPM Algorithm = Self-Dual",
                         NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "LPIPM Stop Tolerance 2 = 1.0e-11",
                         NAGERR_DEFAULT);
  nag_opt_handle_solve_lp_ipm(handle,nvar,x,nnzu,u,rinfo,stats,monit,
                              &comm,NAGERR_DEFAULT);

 END:
  NAG_FREE(cvec);
  NAG_FREE(irowa);
  NAG_FREE(icola);
  NAG_FREE(a);
  NAG_FREE(bla);
  NAG_FREE(bua);
  NAG_FREE(xl);
  NAG_FREE(xu);
  NAG_FREE(x);
  NAG_FREE(u);
  /* nag_opt_handle_free (e04rzc).
   * Destroy the problem handle and deallocate all the memory. */
  if (handle)
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);

  return exit_status;
}


static void NAG_CALL monit(void *handle, const double rinfo[],
                      const double stats[], Nag_Comm *comm,
                      Integer *inform){
  /* Monitoring function */
  double tol = 1.2e-08;

  if (!comm || !comm->iuser){
    /* The communication structure is not correctly allocated, abort solve */
    *inform = -1;
    return;
  }
  /* If x is close to the solution, print a message */
  if (comm->iuser[0]==1){
    if (rinfo[4]<tol && rinfo[5]<tol &&rinfo[6]<tol){
      printf("     Iteration %" NAG_IFMT"\n", (Integer)stats[0]);
      printf("     monit() reports good approximate solution "
             "(tol =, %8.2e):\n",tol);
    }
  }
  else {
    if (rinfo[14]<tol && rinfo[15]<tol &&rinfo[16]<tol){
      printf("     Iteration %" NAG_IFMT"\n", (Integer)stats[0]);
      printf("     monit() reports good approximate solution "
             "(tol =, %8.2e):\n",tol);
    }
  }
  fflush(stdout);
}