/* nag_opt_nlp_revcomm (e04ufc) Example Program.
 *
 * Copyright 2011 Numerical Algorithms Group.
 *
 * Mark 23, 2011
 *
 */
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>

int main(void)
{
  /* Scalars */
  double        objf, nctotal;
  Integer       exit_status=0, i, irevcm, iter, j, n, nclin, ncnln;
  Integer       tda, tdcj, tdr, liwork, lwork;

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

  /* Nag Types */
  NagError     fail;

#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 (e04ufc) Example Program Results\n");
  fflush(stdout);

  /* Skip heading in data file */
  scanf("%*[^\n]");

  scanf("%ld%ld%ld%*[^\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;
      liwork = 3*n + nclin + 2*ncnln;
      lwork = 21*n + 2;
      if (ncnln || nclin) lwork += 2*n*n + 11*nclin;
      if (ncnln) lwork += 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)) ||
          !(iwork = NAG_ALLOC(liwork, Integer)) ||
          !(work = NAG_ALLOC(lwork, 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;

          /* Initialise nag_opt_nlp_revcomm (e04ufc) */
          nag_opt_nlp_revcomm_init("e04ufc",cwsav,5,lwsav,120,iwsav,610,
                                   rwsav,475,&fail);

          /* Set print level */
          nag_opt_nlp_revcomm_option_set_string("Print Level = 10",lwsav,iwsav,
                                                rwsav,&fail);

          /* Solve the problem */
          irevcm = 0;

          do
            {
              /* nag_opt_nlp_revcomm (e04ufc).
               * Solves the nonlinear programming (NP) problem using
               * reverse communication for evaluation of functions.
               */
              nag_opt_nlp_revcomm(&irevcm,n,nclin,ncnln,tda,tdcj,tdr,a,bl,bu,
                                  &iter,istate,c,cjac,clamda,&objf,objgrd,r,x,
                                  needc,iwork,liwork,work,lwork,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("nag_opt_nlp_revcomm (e04ufc) failed.\n%s\n",fail.message);
              exit_status = 1;
            }
        }

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