// The easiest way to pull in all necessary definitions is to include oct.h
#include <octave/oct.h>
// Overwrite type names duplicated in Octave and NAG header files
#define Complex NagComplex
#define MatrixType NagMatrixType
// and include the appropriate NAG header files
#include <nag.h>
#include <nage04.h>

extern "C" {
  static void objfun(Integer n, double x[], double *objf,
                     double objgrd[], Nag_Comm *comm);
  static void confun(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_opt, 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 != 8)
    {
      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());
      NDArray x(args(7).array_value());
      dim_vector dv (1); dv(0) = n;
      NDArray objgrd(dv);
      if (! error_state)
        {
          // Declare local variables
          double objf;
          Nag_Comm comm;
          NagError fail;
          INIT_FAIL(fail);

          // 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(), objfun, confun, 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 objfun(Integer n, double x[], double *objf,
                   double objgrd[], Nag_Comm *comm)
{
  /* Routine to evaluate objective function and its 1st derivatives. */
  if (comm->flag == 0 || comm->flag == 2)
    *objf = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];

  if (comm->flag == 2)
    {
      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]);
    }
} /* objfun */

static void confun(Integer n, Integer ncnlin, Integer needc[],
                   double x[], double conf[], double conjac[],
                   Nag_Comm *comm)
{
  /* Routine to evaluate the nonlinear constraints and
   * their 1st derivatives.
   */
#define CONJAC(I,J) conjac[(I-1)*n + J - 1]

  Integer i, j;

  if (comm->first)
    {
      /* First call to confun.  Set all Jacobian elements to zero.
       * Note that this will only work when con_deriv = TRUE
       * (the default; see Section 4 (confun) and 8.2 (con_deriv)).
       */
      for (j = 1; j <= n; ++j)
        for (i = 1; i <= ncnlin; ++i)
          CONJAC(i,j) = 0.0;
    }
  if (needc[0] > 0)
    {
      if (comm->flag == 0 || comm->flag == 2)
        conf[0] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];

      if (comm->flag == 2)
        {
          CONJAC(1,1) = x[0] * 2.0;
          CONJAC(1,2) = x[1] * 2.0;
          CONJAC(1,3) = x[2] * 2.0;
          CONJAC(1,4) = x[3] * 2.0;
        }
    }
  if (needc[1] > 0)
    {
      if (comm->flag == 0 || comm->flag == 2)
        conf[1] = x[0] * x[1] * x[2] * x[3];

      if (comm->flag == 2)
        {
          CONJAC(2,1) = x[1] * x[2] * x[3];
          CONJAC(2,2) = x[0] * x[2] * x[3];
          CONJAC(2,3) = x[0] * x[1] * x[3];
          CONJAC(2,4) = x[0] * x[1] * x[2];
        }
    }
} /* confun */