/* nag_opt_nlp_revcomm_option_set_file (e04udc) 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)
{
  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;
  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("%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;
      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;

          /* Initialise 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 */
      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 (icomm) NAG_FREE(icomm);
      if (rcomm) NAG_FREE(rcomm);
    }
  return exit_status;
}