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

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

/* NLP example: linear objective + linear constraint + nonlinear constraint.
 * For illustrative purposes, the linear objective will be treated as a
 * nonlinear function to show the usage of objfun and objgrd user call-backs.
 */

#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);
static void NAG_CALL hess(Integer nvar, const double x[], Integer ncnln,
                          Integer idf, double sigma, const double lambda[],
                          Integer nnzh, double hx[], Integer *flag,
                          Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

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

  /* Arrays */
  Integer iuser[2];
  Integer idxfd[4] = {1, 2, 3, 4};
  double ruser[4] = {24.55, 26.75, 39.00, 40.50};
  double rinfo[100], stats[100], *x = 0, *u = 0;
  double *lambda = 0, *fdx = 0, *gdx = 0;
  double bl[4] = {0, 0, 0, 0};
  double bu[4] = {BIGBND, BIGBND, BIGBND, BIGBND};
  double linbl[2] = {5.0, 1.0};
  double linbu[2] = {BIGBND, 1.0};
  double nlnbl[1] = {21.0};
  double nlnbu[1] = {BIGBND};
  Integer irowb[8] = {1, 1, 1, 1, 2, 2, 2, 2};
  Integer icolb[8] = {1, 2, 3, 4, 1, 2, 3, 4};
  Integer irowgd[4] = {1, 1, 1, 1};
  Integer icolgd[4] = {1, 2, 3, 4};
  Integer irowh[10], icolh[10];
  double b[8] = {2.3, 5.6, 11.1, 1.3, 1.0, 1.0, 1.0, 1.0};
  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_ipopt (e04stc) "
                            "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 = 4;

  /* 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 two linear constraints */
  nclin = 2;
  idlc = 0;
  /* nag_opt_handle_set_linconstr (e04rjc).
   * Define a block of linear constraints
   */
  nag_opt_handle_set_linconstr(handle, nclin, linbl, linbu, nclin * nvar, irowb,
                               icolb, b, &idlc, NAGERR_DEFAULT);
  /* Update the Lagrange mult. count */
  nnzu += 2 * nclin;

  /* Add one nonlinear constraint */
  ncnln = 1;

  /* 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;

  /* Define structure of the Hessian, dense upper triangle */
  for (idx = 0, i = 1; i <= nvar; i++)
    for (j = i; j <= nvar; j++, idx++) {
      icolh[idx] = j;
      irowh[idx] = i;
    }
  /* nag_opt_handle_set_nlnhess (e04rlc).
   * Define structure of Hessian of objective, constraints or the Lagrangian
   */
  nag_opt_handle_set_nlnhess(handle, 1, idx, irowh, icolh, NAGERR_DEFAULT);
  nag_opt_handle_set_nlnhess(handle, 0, idx, irowh, icolh, NAGERR_DEFAULT);

  if (!(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] = 1.0;
  iuser[0] = 0;

  nag_opt_handle_opt_set(handle, "Print Level = 0", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Outer Iteration Limit = 26", NAGERR_DEFAULT);
  nag_opt_handle_opt_set(handle, "Stop Tolerance 1 = 2.5e-8", NAGERR_DEFAULT);

  comm.iuser = iuser;
  comm.user = ruser;

  INIT_FAIL(fail);
  /* nag_opt_handle_solve_ipopt (e04stc).
   * Call the solver
   */
  nag_opt_handle_solve_ipopt(handle, objfun, objgrd, confun, congrd, hess,
                             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 constraints 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(7, double)) || !(fdx = NAG_ALLOC(nnzfd, double)) ||
        !(gdx = NAG_ALLOC(nnzgd, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

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

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

    /* Obtain the multipliers with correct sign
     * 4 bound constraints, 2 linear constraints, and 1 nonlinear constraint
     */
    for (i = 0; i < 7; i++) {
      lambda[i] = u[2 * i + 1] - u[2 * i];
    }
    /* lambda (0..3) mult. related to bounds
     * lambda (4..5) mult. related to linear constraints
     * lambda (6) mult. related to the nonlinear constraint
     */

    nag_file_line_write(nout, "Stationarity test for Lagrange multipliers");
    for (i = 0; i < nvar; i++) {
      lmtest = fdx[i] + lambda[i] + lambda[4] * b[i] + lambda[5] * b[4 + i] +
               lambda[6] * gdx[i];
      if (lmtest < 2.5e-8)
        sprintf(msg, "    lx(%10" NAG_IFMT ")%17s=%16.2E%5s%6s", i + 1, "",
                lmtest, "", "Ok    ");

      else
        sprintf(msg, "    lx(%10" NAG_IFMT ")%17s=%16.2E%5s%6s", i + 1, "",
                lmtest, "", "Failed");
      nag_file_line_write(nout, msg);
    }

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

    sprintf(msg, "\nAt solution, 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, "             Complementarity       =%6s%10.2E", "", rinfo[3]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "             KKT error             =%6s%10.2E", "", rinfo[4]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "    after iterations                        :%11.0f",
            stats[0]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "    after objective evaluations             :%11.0f",
            stats[18]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "    after objective gradient evaluations    :%11.0f",
            stats[4]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "    after constraint evaluations            :%11.0f",
            stats[19]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "    after constraint gradient evaluations   :%11.0f",
            stats[20]);
    nag_file_line_write(nout, msg);
    sprintf(msg, "    after hessian evaluations               :%11.0f",
            stats[3]);
    nag_file_line_write(nout, msg);
  } else {
    printf("Error from nag_opt_handle_solve_ipopt (e04stc).\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;
}

/* Subroutine */
static void NAG_CALL objfun(Integer nvar, const double x[], double *fx,
                            Integer *flag, Nag_Comm *comm) {
  *flag = 0;
  /* nag_blast_ddot (f16eac).
   * Dot product of two vectors, allows scaling and accumulation
   */
  nag_blast_ddot(Nag_NoConj, 4, 1.0, x, 1, 0.0, comm->user, 1, fx,
                 NAGERR_DEFAULT);
}
static void NAG_CALL objgrd(Integer nvar, const double x[], Integer nnzfd,
                            double fdx[], Integer *flag, Nag_Comm *comm) {
  *flag = 0;
  /* nag_blast_dge_copy (f16qfc).
   * Matrix copy, real rectangular matrix
   */
  nag_blast_dge_copy(Nag_ColMajor, Nag_NoTrans, nnzfd, 1, comm->user, nnzfd,
                     fdx, nnzfd, NAGERR_DEFAULT);
}
static void NAG_CALL confun(Integer nvar, const double x[], Integer ncnln,
                            double gx[], Integer *flag, Nag_Comm *comm) {
  *flag = 0;
  gx[0] = 12.0 * x[0] + 11.9 * x[1] + 41.8 * x[2] + 52.1 * x[3] -
          1.645 * sqrt(.28 * x[0] * x[0] + .19 * x[1] * x[1] +
                       20.5 * x[2] * x[2] + .62 * x[3] * x[3]);
}
static void NAG_CALL congrd(Integer nvar, const double x[], Integer nnzgd,
                            double gdx[], Integer *flag, Nag_Comm *comm) {
  double tmp;
  *flag = 0;
  tmp = sqrt(0.62 * x[3] * x[3] + 20.5 * x[2] * x[2] + 0.19 * x[1] * x[1] +
             0.28 * x[0] * x[0]);
  gdx[0] = (12.0 * tmp - 0.4606 * x[0]) / tmp;
  gdx[1] = (11.9 * tmp - 0.31255 * x[1]) / tmp;
  gdx[2] = (41.8 * tmp - 33.7225 * x[2]) / tmp;
  gdx[3] = (52.1 * tmp - 1.0199 * x[3]) / tmp;
}
static void NAG_CALL hess(Integer nvar, const double x[], Integer ncnln,
                          Integer idf, double sigma, const double lambda[],
                          Integer nnzh, double hx[], Integer *flag,
                          Nag_Comm *comm) {
  double tmp;
  *flag = 0;
  /* nag_blast_dload (f16fbc).
   * Broadcast scalar into real vector
   */
  nag_blast_dload(nnzh, 0.0, hx, 1, NAGERR_DEFAULT);

  if (idf == 0)
    return; /* objective is linear */
  tmp = sqrt(0.62 * x[3] * x[3] + 20.5 * x[2] * x[2] + 0.19 * x[1] * x[1] +
             0.28 * x[0] * x[0]);
  tmp = tmp * (x[3] * x[3] + 33.064516129032258064 * x[2] * x[2] +
               0.30645161290322580645 * x[1] * x[1] +
               0.45161290322580645161 * x[0] * x[0]);
  /* Col 1 */
  hx[0] = (-0.4606 * x[3] * x[3] - 15.229516129032258064 * x[2] * x[2] -
           0.14115161290322580645 * x[1] * x[1]) /
          tmp;
  hx[1] = (0.14115161290322580645 * x[0] * x[1]) / tmp;
  hx[2] = (15.229516129032258064 * x[0] * x[2]) / tmp;
  hx[3] = (0.4606 * x[0] * x[3]) / tmp;
  /* Col 2 */
  hx[4] = (-0.31255 * x[3] * x[3] - 10.334314516129032258 * x[2] * x[2] -
           0.14115161290322580645 * x[0] * x[0]) /
          tmp;
  hx[5] = (10.334314516129032258 * x[1] * x[2]) / tmp;
  hx[6] = (0.31255 * x[1] * x[3]) / tmp;
  /* Col 3 */
  hx[7] = (-33.7225 * x[3] * x[3] - 10.334314516129032258 * x[1] * x[1] -
           15.229516129032258065 * x[0] * x[0]) /
          tmp;
  hx[8] = (33.7225 * x[2] * x[3]) / tmp;
  /* Col 4 */
  hx[9] =
      (-33.7225 * x[2] * x[2] - 0.31255 * x[1] * x[1] - 0.4606 * x[0] * x[0]) /
      tmp;
  /* nag_blast_daxpby (f16ecc).
   * Real weighted vector addition
   */
  if (idf == -1)
    nag_blast_daxpby(nnzh, 0.0, hx, 1, lambda[0], hx, 1, NAGERR_DEFAULT);
}