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

NAG CL Interface Introduction
Example description
/* nag_opt_handle_solve_ssqp (e04src) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.4, 2022.
 */

/* NLP example: Nonlinear objective, linear constraint and two nonlinear
 * constraints.
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
                            Integer *flag, Nag_Comm *comm);
static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
                            double fdx[], Integer *flag, Nag_Comm *comm);
static void NAG_CALL confun(Integer nvar, const double x[], Integer ncnln,
                            double gx[], Integer *flag, Nag_Comm *comm);
static void NAG_CALL congrd(Integer nvar, const double x[], Integer nnzgd,
                            double gdx[], Integer *flag, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void) {
#define BIGBND 1.0E20
  /* Scalars */
  Integer exit_status = 0;
  Integer i, idlc, nnzu, nvar, nclin, ncnln, nnzfd, nnzgd, flag;
  double lmtest[4];

  /* Arrays */
  double rinfo[100], stats[100], *x = 0, *u = 0;
  double *lambda = 0, *fdx = 0, *gdx = 0;
  Integer idxfd[4] = {1, 2, 3, 4};
  double bl[4] = {-BIGBND, -BIGBND, 0, 0};
  double bu[4] = {BIGBND, BIGBND, BIGBND, BIGBND};
  double linbl[1] = {0.0};
  double linbu[1] = {BIGBND};
  double nlnbl[2] = {2.0, 4.0};
  double nlnbu[2] = {2.0, 4.0};
  Integer irowb[2] = {1, 1};
  Integer icolb[2] = {1, 2};
  Integer irowgd[5] = {1, 1, 1, 2, 2};
  Integer icolgd[5] = {1, 2, 3, 2, 4};
  double b[2] = {2.0, 4.0};
  Integer xlin[1] = {4};
  void *handle = 0;

  /* Nag Types */
  NagError fail;
  Nag_FileID nout;
  Nag_Comm comm;

  /* nag_file_open (x04acc).
   *  Open unit number for reading, writing or appending, and
   *  associate unit with named file
   */
  nag_file_open("", 1, &nout, NAGERR_DEFAULT);

  /* nag_file_line_write (x04bac).
   * Write formatted record to external file
   */
  nag_file_line_write(nout, "nag_opt_handle_solve_ssqp (e04src) "
                            "Example Program Results");

  /* Problem size */
  nvar = 4;
  /* Counter for Lagrange multipliers */
  nnzu = 0;
  /* Objective gradient nonzero elements quantity */
  nnzfd = 4;
  /* Constraint jacobian nonzero elements quantity */
  nnzgd = 5;

  /* 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_simplebounds (e04rhc).
   * Define bounds on the variables
   */
  nag_opt_handle_set_simplebounds(handle, nvar, bl, bu, NAGERR_DEFAULT);
  /* Specify the amount of Lagrange mult. required */
  nnzu += 2 * nvar;

  /* nag_opt_handle_set_nlnobj (e04rgc).
   * Define nonlinear objective
   */
  nag_opt_handle_set_nlnobj(handle, nnzfd, idxfd, NAGERR_DEFAULT);

  /* Add one linear constraint */
  nclin = 1;
  idlc = 0;
  /* nag_opt_handle_set_linconstr (e04rjc).
   * Define a block of linear constraints
   */
  nag_opt_handle_set_linconstr(handle, nclin, linbl, linbu, 2, irowb,
                               icolb, b, &idlc, NAGERR_DEFAULT);
  /* Update the Lagrange mult. count */
  nnzu += 2 * nclin;

  /* Add two nonlinear constraints */
  ncnln = 2;

  /* nag_opt_handle_set_nlnconstr (e04rkc).
   * Define a block of nonlinear constraints
   */
  nag_opt_handle_set_nlnconstr(handle, ncnln, nlnbl, nlnbu, nnzgd, irowgd,
                               icolgd, NAGERR_DEFAULT);
  /* Update the Lagrange mult. count */
  nnzu += 2 * ncnln;

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

  /* Define initial guess point */
  for (i = 0; i < nvar; i++)
    x[i] = i+1;


  /* Optionally, define variable x(4) to be linear
   * Hinting which variables are linear in problems with many
   * variables can speed up performance
   */
  nag_opt_handle_set_property(handle,"linear",1,xlin,NAGERR_DEFAULT);

  /* Disable printing */
  nag_opt_handle_opt_set(handle, "Print Level = 0", NAGERR_DEFAULT);
  /* Do not print options */
  nag_opt_handle_opt_set(handle, "Print Options = No", NAGERR_DEFAULT);
  /* It is recommended on new problems to always verify the derivatives */
  nag_opt_handle_opt_set(handle, "Verify Derivatives = Yes", NAGERR_DEFAULT);
  /* Do not print solution, x and f(x) will be printed afterwards */
  nag_opt_handle_opt_set(handle, "Print Solution = No", NAGERR_DEFAULT);

  INIT_FAIL(fail);
  /* nag_opt_handle_solve_ssqp (e04src).
   * Call the solver
   */
  nag_opt_handle_solve_ssqp(handle, objfun, objgrd, confun, congrd, NULLFN, 
      NULLFN, nvar, x, nnzu, u, rinfo, stats, &comm, &fail);
  if (fail.code == NE_NOERROR) {
    char msg[80];
    nag_file_line_write(nout, "\nVariables");
    for (i = 0; i < nvar; i++) {
      sprintf(msg, "     x(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "", x[i]);
      nag_file_line_write(nout, msg);
    }
    nag_file_line_write(nout, "Variable bound Lagrange multipliers");
    for (i = 0; i < nvar; i++) {
      sprintf(msg, "    zL(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "", u[2 * i]);
      nag_file_line_write(nout, msg);
      sprintf(msg, "    zU(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
              u[2 * i + 1]);
      nag_file_line_write(nout, msg);
    }
    nag_file_line_write(nout, "Linear constraint Lagrange multipliers");
    for (i = 0; i < nclin; i++) {
      sprintf(msg, "    zL(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
              u[2 * nvar + 2 * i]);
      nag_file_line_write(nout, msg);
      sprintf(msg, "    zU(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
              u[2 * nvar + 2 * i + 1]);
      nag_file_line_write(nout, msg);
    }
    nag_file_line_write(nout, "Nonlinear constraints Lagrange multipliers");
    for (i = 0; i < ncnln; i++) {
      sprintf(msg, "    zL(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
              u[2 * (nvar + nclin) + 2 * i]);
      nag_file_line_write(nout, msg);
      sprintf(msg, "    zU(%10" NAG_IFMT ")%17s=%16.2E", i + 1, "",
              u[2 * (nvar + nclin) + 2 * i + 1]);
      nag_file_line_write(nout, msg);
    }

    if (!(lambda = NAG_ALLOC((Integer)nnzu/2, double)) ||
        !(fdx = NAG_ALLOC(nnzfd, double)) ||
        !(gdx = NAG_ALLOC(nnzgd, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

    /* Check of Lagrange multipliers (stationarity)
     * Gradient of the objective (also available in data%coeffs(1:nnzfd))
     */
    flag = 0;
    objgrd(nvar, x, nnzfd, fdx, &flag, &comm);
    if (flag != 0) {
      exit_status = -1;
      goto END;
    }

    /* Gradient of the nonlinear constraint */
    flag = 0;
    congrd(nvar, x, nnzgd, gdx, &flag, &comm);
    if (flag != 0) {
      exit_status = -1;
      goto END;
    }

    /* Obtain the multipliers with correct sign
     * 4 bound constraints, 1 linear constraint, and 2 nonlinear constraints
     */
    for (i = 0; i < 7; i++) {
      lambda[i] = u[2 * i + 1] - u[2 * i];
    }
    /* lambda (0..3) mult. related to bounds
     * lambda (4) mult. related to linear constraints
     * lambda (5..6) mult. related to the nonlinear constraint
     */
    nag_file_line_write(nout, "Stationarity test for Lagrange multipliers");
    lmtest[0] = fdx[0] + lambda[0] + lambda[4]*b[0] + lambda[5]*gdx[0];
    lmtest[1] = fdx[1] + lambda[1] + lambda[4]*b[1] + lambda[5]*gdx[1]
      + lambda[6]*gdx[3];
    lmtest[2] = fdx[2] + lambda[2] +                  lambda[5]*gdx[2];
    lmtest[3] = fdx[3] + lambda[3]
      + lambda[6]*gdx[4];
    for (i = 0; i < nvar; i++) {
      if (lmtest[i] < 2.5e-8)
        sprintf(msg, "    lx(%10" NAG_IFMT ")%17s=%16.2E%5s%6s", i + 1, "",
                lmtest[i], "", "Ok    ");
      else
        sprintf(msg, "    lx(%10" NAG_IFMT ")%17s=%16.2E%5s%6s", i + 1, "",
                lmtest[i], "", "Failed");
      nag_file_line_write(nout, msg);
    }

    NAG_FREE(lambda);
    NAG_FREE(fdx);
    NAG_FREE(gdx);

    sprintf(msg, "\nSolution found. Objective minimum  =%16.4E", rinfo[0]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "             Constraint violation  =%6s%10.2E", "", rinfo[1]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "             Dual infeasibility    =%6s%10.2E", "", rinfo[2]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "   Iterations                      =     %11" NAG_IFMT,
            (Integer) stats[0]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "   Objective evaluations           =     %11" NAG_IFMT,
            (Integer) stats[18]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "   Objective gradient evaluations  =     %11" NAG_IFMT,
            (Integer) stats[4]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "   Constraint evaluations          =     %11" NAG_IFMT,
            (Integer) stats[19]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "   Constraint gradient evaluations =     %11" NAG_IFMT,
            (Integer) stats[20]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "   Hessian evaluations             =     %11" NAG_IFMT,
            (Integer) stats[3]);
    nag_file_line_write(nout, msg);
  } else {
    printf("Error from nag_opt_handle_solve_ssqp (e04src).\n%s\n",
           fail.message);
    exit_status = 1;
  }
  if (handle)
    /* nag_opt_handle_free (e04rzc).
     * Destroy the problem handle and deallocate all the memory used
     */
    nag_opt_handle_free(&handle, NAGERR_DEFAULT);

  NAG_FREE(x);
  NAG_FREE(u);

END:
  return exit_status;
}

/* Subroutines */
static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
                            Integer *flag, Nag_Comm *comm) {
  *fx = (x[0]+x[1]+x[2])*(x[0]+x[1]+x[2]) + 3.0*x[2] + 5.0*x[3] + cos(0.01*x[0]) - 1.0;

  *flag = 0;
}
static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
                            double fdx[], Integer *flag, Nag_Comm *comm) {
  double sum;
  sum = 2.0*(x[0]+x[1]+x[2]);

  fdx[0] = sum - 0.01*sin(0.01*x[0]);
  fdx[1] = sum;
  fdx[2] = sum + 3.0;
  fdx[3] = 5.0;

  *flag = 0;
}
static void NAG_CALL confun(Integer nvar, const double x[], Integer ncnln,
                            double gx[], Integer *flag, Nag_Comm *comm) {
  gx[0] = x[0]*x[0] + x[1]*x[1] + x[2];
  gx[1] = x[1]*x[1]*x[1]*x[1] + x[3];

  *flag = 0;
}
static void NAG_CALL congrd(Integer nvar, const double x[], Integer nnzgd,
                            double gdx[], Integer *flag, Nag_Comm *comm) {
  gdx[0] = 2.0*x[0];
  gdx[1] = 2.0*x[1];
  gdx[2] = 1.0;
  gdx[3] = 4.0*x[1]*x[1]*x[1];
  gdx[4] = 1.0;
  
  *flag = 0;
}