/* The easiest way to pull in all necessary 
   definitions is to include OCT/src/oct.h */
#include <octave/oct.h>
#include <octave/parse.h>

/* Overwrite type names duplicated in Octave and NAG header files */
#define Complex NagComplex
#define MatrixType NagMatrixType
// and to include the appropriate NAG header files
#include <nag.h>
#include <nage04.h>

extern "C" {
  static void objfunC(Integer n, double x[], double *objf,
                      double objgrd[], Nag_Comm *comm);
  static void confunC(Integer n, Integer ncnlin, Integer needc[],
                      double x[], double conf[], double conjac[],
                      Nag_Comm *comm);
}

// The third parameter (nargout) is not used, so it is
// omitted from the list of arguments to DEFUN_DLD in order to avoid
// the warning from gcc about an unused function parameter. 
DEFUN_DLD (nag_opt2, args, ,
  "Calls nag_opt_nlp (e04ucc), which minimizes an arbitrary smooth function \n\
subject to constraints using a sequential quadratic programming (SQP) method.\n\
As many first derivatives as possible should be supplied by you; \n\
any unspecified derivatives are approximated by finite differences. \n\
It is not intended for large sparse problems.\n\
nag_opt_nlp (e04ucc) may also be used for unconstrained, bound-constrained \n\
and linearly constrained optimization.\n")
{
  // Variable to store function output values
  octave_value_list retval;
  int nargin = args.length ();
  if (nargin != 10)
    {
      error ("Insufficient input arguments.");
    }
  else
    {
      // Retrieve input arguments from args
      Integer n = args(0).int_value();
      Integer nclin = args(1).int_value();
      Integer ncnlin = args(2).int_value();
      Matrix a(args(3).matrix_value());
      Integer tda = args(4).int_value();
      NDArray bl(args(5).array_value());
      NDArray bu(args(6).array_value());
      struct fcnptrs {octave_function *obj, *con;} fcns;
      fcns.obj = args(7).function_value();
      fcns.con = args(8).function_value();
      NDArray x(args(9).array_value());
      dim_vector dv (1); dv(0) = n;
      NDArray objgrd(dv);
      if (! error_state)
        {
          // Declare local variables
          double objf=0.0;
          Nag_Comm comm;
          NagError fail;
          INIT_FAIL(fail);
          // Store Octave functions handles in communication structure comm
          comm.p = (Pointer) &fcns;
          
          // Call NAG routine
          /* nag_opt_nlp (e04ucc).
           * Minimization with nonlinear constraints using a
           * sequential QP method
           */
          nag_opt_nlp(n, nclin, ncnlin, a.fortran_vec(), tda, bl.fortran_vec(), 
                      bu.fortran_vec(),objfunC, confunC, x.fortran_vec(),&objf,
                      objgrd.fortran_vec(), E04_DEFAULT, &comm, &fail);
          /* T* fortran_vec (void) method returns a pointer to the underlying data of 
             the matrix or a array so that it can be manipulated directly
             by an external library. */
          
          if (fail.code != NE_NOERROR)
            {
              error("Error from nag_opt_nlp (e04ucc).\n%s\n", fail.message);
              goto END;
            }
          // Assign output arguments to retval
          retval(0) = objf;
          retval(1) = objgrd;
          retval(2) = fail.code;
        } 
    }

 END:
  return retval;
}

static void objfunC(Integer n, double x[], double *objf,
                    double objgrd[], Nag_Comm *comm)
{
  // Retrieve Octave objfun handle
  struct fcnptrs {octave_function *obj, *con;} *fcns;
  fcns = (struct fcnptrs *) comm->p;
  octave_function *fcn = (octave_function *) fcns->obj;
  // Copy C arrays into the Octave type ones
  dim_vector dv (1); dv(0) = n;
  NDArray x_oct(dv), objgrd_oct(dv);
  octave_idx_type i;
  for (i=0;i<n;i++) {
    x_oct.elem(i) = x[i];
    objgrd_oct.elem(i) = objgrd[i];
  }
  // Assign objfun input arguments
  octave_value_list newargs, retval;
  newargs(0)=octave_value(comm->flag);
  newargs(1)=octave_value(x_oct);
  newargs(2)=octave_value(*objf);
  newargs(3)=octave_value(objgrd_oct);
  int nargout=2;
  // Call Octave objfun function
  retval=feval(fcn,newargs,nargout);
  // Retrieve output arguments from retval
  *objf=retval(0).double_value();
  objgrd_oct = retval(1).array_value();
  // Copy output arrays back into C type arrays
  for (i=0;i<n;i++) {
    objgrd[i]=objgrd_oct.elem(i);
    x[i]=x_oct.elem(i);
  }
} /* objfunC */

static void confunC(Integer n, Integer ncnlin, Integer needc[],
                    double x[], double conf[], double conjac[],
                    Nag_Comm *comm)
{
#define CONJAC(I,J) conjac[(I-1)*n + J - 1]

  // Retrieve Octave confun handle
  struct fcnptrs {octave_function *obj, *con;} *fcns;
  fcns = (struct fcnptrs *) comm->p;
  octave_function *fcn = (octave_function *) fcns->con;
  // Copy C arrays into the Octave type arrays and matrices
  dim_vector dv (1); dv(0) = ncnlin;
  dim_vector dv2 (1); dv2(0) = n;
  NDArray needc_oct(dv), conf_oct(dv), x_oct(dv2);
  Matrix conjac_oct(ncnlin,n);
  octave_idx_type i,j;
  for (i=0;i<ncnlin;i++) {
    needc_oct.elem(i) = needc[i];
    conf_oct.elem(i) = conf[i];
    for (j=0;j<n;j++) {
      x_oct.elem(j) = x[j];
      conjac_oct.elem(i,j) = CONJAC(i+1,j+1);
    }
  }
  // Assign confun input arguments
  octave_value_list newargs, retval;
  newargs(0)=octave_value(comm->flag);
  newargs(1)=octave_value(ncnlin);
  newargs(2)=octave_value(n);
  newargs(3)=octave_value(needc_oct);
  newargs(4)=octave_value(x_oct);
  newargs(5)=octave_value(conjac_oct);
  newargs(6)=octave_value(comm->first);
  int nargout=2;
  // Call Octave confun function
  retval=feval(fcn,newargs,nargout);
  // Retrieve output arguments from retval
  conf_oct=retval(0).array_value();
  conjac_oct = retval(1).matrix_value();
  // Copy output arrays and matrices back into C type arrays
  for (i=0;i<ncnlin;i++) {
    needc[i] = needc_oct.elem(i);
    conf[i] = conf_oct.elem(i);
    for (j=0;j<n;j++) {
      x[j]=x_oct.elem(j);
      CONJAC(i+1,j+1) = conjac_oct.elem(i,j);
    }
  }
} /* confunC */