/* nag_opt_nlp_revcomm_option_set_file (e04udc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 *
 */
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>

int main(void)
{
  const char *optionsfile = "e04udce.opt";

  /* Scalars */
  double objf, nctotal;
  Integer exit_status = 0, i, irevcm, iter, j, n, nclin, ncnln;
  Integer tda, tdcj, tdr, licomm, lrcomm;

  /* Arrays */
  double *a = 0, *bl = 0, *bu = 0, *c = 0, *cjac = 0, *clamda = 0, *objgrd =
         0;
  double *r = 0, *rcomm = 0, rwsav[475], *x = 0;
  Integer *icomm = 0, *istate = 0, iwsav[610], *needc = 0;
  char cwsav[5 * 80];

  /* Nag Types */
  Nag_Boolean lwsav[120];
  NagError fail;
  Nag_FileID optfileid;

#define A(I,J) a[(I-1)*tda + J - 1]
#define CJAC(I,J) cjac[(I-1)*tdcj + J - 1]

  INIT_FAIL(fail);

  printf("nag_opt_nlp_revcomm_option_set_file (e04udc) "
         "Example Program Results\n");
  fflush(stdout);

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &n, &nclin,
        &ncnln);

  if (n <= 0 || nclin < 0 || ncnln < 0) {
    printf("At least one of n, nclin, or ncnln is invalid\n");
    exit_status = 1;
  }
  else {
    tda = MAX(nclin, n);
    tdcj = MAX(ncnln, n);
    tdr = n;
    nctotal = n + nclin + ncnln;
    licomm = 3 * n + nclin + 2 * ncnln;
    lrcomm = 21 * n + 2;
    if (ncnln || nclin)
      lrcomm += 2 * n * n + 11 * nclin;
    if (ncnln)
      lrcomm += n * nclin + 2 * n * ncnln + 22 * ncnln - 1;

    /* Allocate memory */
    if (!(a = NAG_ALLOC(tda * MAX(1, nclin), double)) ||
        !(bl = NAG_ALLOC(nctotal, double)) ||
        !(bu = NAG_ALLOC(nctotal, double)) ||
        !(istate = NAG_ALLOC(nctotal, Integer)) ||
        !(c = NAG_ALLOC(ncnln, double)) ||
        !(cjac = NAG_ALLOC(tdcj * MAX(1, ncnln), double)) ||
        !(clamda = NAG_ALLOC(nctotal, double)) ||
        !(objgrd = NAG_ALLOC(n, double)) ||
        !(r = NAG_ALLOC(tdr * n, double)) ||
        !(x = NAG_ALLOC(n, double)) ||
        !(needc = NAG_ALLOC(ncnln, Integer)) ||
        !(icomm = NAG_ALLOC(licomm, Integer)) ||
        !(rcomm = NAG_ALLOC(lrcomm, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
    }
    else {
      /* Read A, BL, BU and X from data file */
      if (nclin > 0) {
        for (i = 1; i <= nclin; ++i)
          for (j = 1; j <= n; ++j)
            scanf("%lf", &A(i, j));
        scanf("%*[^\n] ");
      }

      for (i = 0; i < nctotal; ++i)
        scanf("%lf", &bl[i]);
      scanf("%*[^\n] ");
      for (i = 0; i < nctotal; ++i)
        scanf("%lf", &bu[i]);
      scanf("%*[^\n] ");
      for (i = 0; i < n; ++i)
        scanf("%lf", &x[i]);

      /* Set all constraint Jacobian elements to zero.
         Note that this will only work when 'Derivative Level = 3'
         (the default; see Section 11.2) */

      for (j = 1; j <= n; ++j)
        for (i = 1; i <= ncnln; ++i)
          CJAC(i, j) = 0.0;

      /* Initialize nag_opt_nlp_revcomm (e04ufc) and check for error
         exits */
      nag_opt_nlp_revcomm_init("e04ufc", cwsav, 5, lwsav, 120, iwsav, 610,
                               rwsav, 475, &fail);

      /* Set three options using
         nag_opt_nlp_revcomm_option_set_string (e04uec) */
      nag_opt_nlp_revcomm_option_set_string("Infinite Bound Size = 1.0d+25",
                                            lwsav, iwsav, rwsav, &fail);
      nag_opt_nlp_revcomm_option_set_string("Print Level = 1", lwsav, iwsav,
                                            rwsav, &fail);
      nag_opt_nlp_revcomm_option_set_string("Verify Level = 11", lwsav, iwsav,
                                            rwsav, &fail);

      /* Use nag_opt_sparse_nlp_option_set_file (e04vkc) to read some
       * options from the options file. Call nag_open_file (x04acc) to
       * set the options file optfileid */
      /* nag_open_file (x04acc), see above. */
      nag_open_file(optionsfile, 0, &optfileid, &fail);

      /* nag_opt_nlp_revcomm_option_set_file (e04udc).
       * Supply optional parameter values for
       * nag_opt_nlp_revcomm (e04ufc) from external file
       */
      nag_opt_nlp_revcomm_option_set_file(optfileid, lwsav, iwsav, rwsav,
                                          &fail);

      /* Solve the problem */
      irevcm = 0;

      do {
        nag_opt_nlp_revcomm(&irevcm, n, nclin, ncnln, tda, tdcj, tdr, a,
                            bl, bu, &iter, istate, c, cjac, clamda, &objf,
                            objgrd, r, x, needc, icomm, licomm, rcomm, lrcomm,
                            cwsav, lwsav, iwsav, rwsav, &fail);

        if (irevcm == 1 || irevcm == 3)
          /* Evaluate the objective function */
          objf = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];

        if (irevcm == 2 || irevcm == 3) {
          /* Evaluate the objective gradient */
          objgrd[0] = x[3] * (2.0 * x[0] + x[1] + x[2]);
          objgrd[1] = x[0] * x[3];
          objgrd[2] = x[0] * x[3] + 1.0;
          objgrd[3] = x[0] * (x[0] + x[1] + x[2]);
        }

        if (irevcm == 4 || irevcm == 6) {
          /* Evaluate the nonlinear constraint functions */
          if (needc[0] > 0)
            c[0] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
          if (needc[1] > 0)
            c[1] = x[0] * x[1] * x[2] * x[3];
        }

        if (irevcm == 5 || irevcm == 6) {
          /* Evaluate the constraint Jacobian */
          if (needc[0] > 0) {
            CJAC(1, 1) = 2.0 * x[0];
            CJAC(1, 2) = 2.0 * x[1];
            CJAC(1, 3) = 2.0 * x[2];
            CJAC(1, 4) = 2.0 * x[3];
          }

          if (needc[1] > 0) {
            CJAC(2, 1) = x[1] * x[2] * x[3];
            CJAC(2, 2) = x[0] * x[2] * x[3];
            CJAC(2, 3) = x[0] * x[1] * x[3];
            CJAC(2, 4) = x[0] * x[1] * x[2];
          }
        }
      } while (irevcm > 0);

      if (fail.code != NE_NOERROR) {
        printf("e04ufc failed.\n%s\n", fail.message);
        exit_status = 1;
      }
    }

    /* Deallocate any allocated arrays */
    NAG_FREE(a);
    NAG_FREE(bl);
    NAG_FREE(bu);
    NAG_FREE(istate);
    NAG_FREE(c);
    NAG_FREE(cjac);
    NAG_FREE(clamda);
    NAG_FREE(objgrd);
    NAG_FREE(r);
    NAG_FREE(x);
    NAG_FREE(needc);
    NAG_FREE(icomm);
    NAG_FREE(rcomm);
  }
  return exit_status;
}