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

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

#include <math.h>
#include <nag.h>
#include <stdio.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 n, nnzq0, nnzq1;
  Integer nnzu, nnzuc, exit_status, i, x_idx;
  Integer idqc;
  Integer verbose_output;
  Integer tol_reached;
  Integer *irowq0 = 0, *icolq0 = 0, *irowq1 = 0, *icolq1 = 0, *idxr0 = 0,
          *idxr1 = 0;
  double *x = 0, *u = 0, *uc = 0;
  double *q0 = 0, *q1 = 0, *r0 = 0, *r1 = 0;
  double rinfo[100], stats[100];
  double s, tol_monit;
  void *handle = 0;
  /* Nag Types */
  Nag_Comm comm;
  NagError fail;

  exit_status = 0;

  printf("nag_opt_handle_set_qconstr (e04rsc) 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 " %*[^\n]", &n, &nnzq0,
        &nnzq1);

  /* Allocate memory to read data */
  if (!(irowq0 = NAG_ALLOC(nnzq0, Integer)) ||
      !(icolq0 = NAG_ALLOC(nnzq0, Integer)) ||
      !(irowq1 = NAG_ALLOC(nnzq1, Integer)) ||
      !(icolq1 = NAG_ALLOC(nnzq1, Integer)) ||
      !(idxr0 = NAG_ALLOC(n, Integer)) || !(idxr1 = NAG_ALLOC(n, Integer)) ||
      !(q0 = NAG_ALLOC(nnzq0, double)) || !(q1 = NAG_ALLOC(nnzq1, double)) ||
      !(r0 = NAG_ALLOC(n, double)) || !(r1 = NAG_ALLOC(n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Read Q0 matrix row indices */
  for (i = 0; i < nnzq0; i++) {
    scanf("%" NAG_IFMT, &irowq0[i]);
  }
  scanf("%*[^\n]");
  /* Read Q0 matrix column indices */
  for (i = 0; i < nnzq0; i++) {
    scanf("%" NAG_IFMT, &icolq0[i]);
  }
  scanf("%*[^\n]");
  /* Read Q0 values*/
  for (i = 0; i < nnzq0; i++) {
    scanf("%lf", &q0[i]);
  }
  scanf("%*[^\n]");
  /* Read Q1 matrix row indices */
  for (i = 0; i < nnzq1; i++) {
    scanf("%" NAG_IFMT, &irowq1[i]);
  }
  scanf("%*[^\n]");
  /* Read Q1 matrix column indices */
  for (i = 0; i < nnzq1; i++) {
    scanf("%" NAG_IFMT, &icolq1[i]);
  }
  scanf("%*[^\n]");
  /* Read Q1 values*/
  for (i = 0; i < nnzq1; i++) {
    scanf("%lf", &q1[i]);
  }
  scanf("%*[^\n]");
  /* Read idxr0 values*/
  for (i = 0; i < n; i++) {
    scanf("%" NAG_IFMT, &idxr0[i]);
  }
  scanf("%*[^\n]");
  /* Read r0 values*/
  for (i = 0; i < n; i++) {
    scanf("%lf", &r0[i]);
  }
  scanf("%*[^\n]");
  /* Read idxr1 values*/
  for (i = 0; i < n; i++) {
    scanf("%" NAG_IFMT, &idxr1[i]);
  }
  scanf("%*[^\n]");
  /* Read r1 values*/
  for (i = 0; i < n; i++) {
    scanf("%lf", &r1[i]);
  }
  scanf("%*[^\n]");
  /* Read s */
  scanf("%lf", &s);

  /* Compute size of multipliers */
  /* One quadratic constraint in the model will have
   * 2 multipliers for both bounds */
  nnzu = 2;
  /* No cone constraint in the model, so set nnzuc to 0 */
  nnzuc = 0;

  /* Allocate memory */
  if (!(x = NAG_ALLOC(n, double)) || !(u = NAG_ALLOC(nnzu, double)) ||
      !(uc = NAG_ALLOC(nnzuc, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  for (i = 0; i < n; i++) {
    x[i] = 0.0;
  }

  /* Create the problem handle */
  /* nag_opt_handle_init (e04rac).
   * Initialize an empty problem handle with n variables. */
  nag_opt_handle_init(&handle, n, NAGERR_DEFAULT);

  /* nag_opt_handle_set_qconstr (e04rsc)
   * Define a quadratic objective */
  idqc = -1;
  nag_opt_handle_set_qconstr(handle, 0.0, n, idxr0, r0, nnzq0, irowq0, icolq0,
                             q0, &idqc, NAGERR_DEFAULT);

  /* nag_opt_handle_set_qconstr (e04rsc)
   * Define a quadratic constraint */
  idqc = 0;
  nag_opt_handle_set_qconstr(handle, s, n, idxr1, r1, nnzq1, irowq1, icolq1, q1,
                             &idqc, NAGERR_DEFAULT);

  /* nag_opt_handle_opt_set (e04zmc) */
  /* Turn on monitoring */
  nag_opt_handle_opt_set(handle, "SOCP Monitor Frequency = 1", NAGERR_DEFAULT);

  /* Set this to 1 to cause nag_opt_handle_solve_socp_ipm (e04ptc)
   * to produce intermediate progress output */
  verbose_output = 0;

  if (verbose_output) {
    /* Require printing of primal and dual solutions at the end of the solve */
    nag_opt_handle_opt_set(handle, "Print Solution = Yes", NAGERR_DEFAULT);
  } else {
    /* Turn off printing of intermediate progress output */
    nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
  }

  tol_reached = 0;
  tol_monit = 1.0e-7;
  comm.iuser = &tol_reached;
  comm.user = &tol_monit;
  /* nag_opt_handle_solve_socp_ipm (e04ptc) */
  INIT_FAIL(fail);
  nag_opt_handle_solve_socp_ipm(handle, n, x, nnzu, u, nnzuc, uc, rinfo, stats,
                                monit, &comm, &fail);
  if (fail.code != NE_NOERROR && fail.code != NW_NOT_CONVERGED) {
    printf("nag_opt_handle_solve_socp_ipm (e04ptc) failed.\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Print solution if optimal or suboptimal solution found */
  printf(" Optimal X:\n");
  printf("  x_idx       Value\n");
  for (x_idx = 0; x_idx < n; x_idx++) {
    printf("  %5" NAG_IFMT "   %11.3E\n", x_idx + 1, x[x_idx]);
  }

END:
  NAG_FREE(x);
  NAG_FREE(u);
  NAG_FREE(uc);
  NAG_FREE(irowq0);
  NAG_FREE(icolq0);
  NAG_FREE(irowq1);
  NAG_FREE(icolq1);
  NAG_FREE(q0);
  NAG_FREE(q1);
  NAG_FREE(r0);
  NAG_FREE(r1);
  NAG_FREE(idxr0);
  NAG_FREE(idxr1);

  /* 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 can be used to monitor the progress
   *  or, for example,  to implement bespoke stopping criteria */
  double tol = comm->user[0];
  Integer *tol_reached = comm->iuser;

  /* If x is close to the solution, print a message */
  if (rinfo[14] < tol && rinfo[15] < tol && rinfo[16] < tol &&
      rinfo[17] < tol) {
    if (!*tol_reached) {
      printf("\n     monit() reports good approximate solution "
             "(tol = %8.2E)\n",
             tol);
      *tol_reached = 1;
    }
  }
  fflush(stdout);
}