/* nag_opt_nlp_solve (e04wdc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 8, 2004.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL confun(Integer *mode, Integer ncnln, Integer n,
Integer ldcj, const Integer needc[],
const double x[], double ccon[], double cjac[],
Integer nstate, Nag_Comm *comm);
static void NAG_CALL objfun(Integer *mode, Integer n, const double x[],
double *objf, double grad[], Integer nstate,
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{

/* Scalars */
double       objf;
Integer      exit_status, i, j, majits, n, nclin, ncnln, nctotal, pda, pdcj,
pdh;

/* Arrays */
static double ruser[2] = {-1.0, -1.0};
double       *a = 0, *bl = 0, *bu = 0, *ccon = 0, *cjac = 0, *clamda = 0;
double       *grad = 0, *h = 0, *x = 0;
Integer      *istate = 0;

/* Nag Types */
Nag_E04State state;
NagError     fail;
Nag_Comm     comm;
Nag_FileID   fileid;

#define A(I, J) a[(I-1)*pda + J - 1]

exit_status = 0;
INIT_FAIL(fail);

printf("nag_opt_nlp_solve (e04wdc) Example Program Results\n");

/* For communication with user-supplied functions: */
comm.user = ruser;

/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%ld%ld%ld%*[^\n] ", &n, &nclin, &ncnln);
if (n > 0 && nclin >= 0 && ncnln >= 0)
{
/* Allocate memory */
nctotal = n + nclin + ncnln;
if (!(a = NAG_ALLOC(nclin*n, double)) ||
!(bl = NAG_ALLOC(nctotal, double)) ||
!(bu = NAG_ALLOC(nctotal, double)) ||
!(ccon = NAG_ALLOC(ncnln, double)) ||
!(cjac = NAG_ALLOC(ncnln*n, double)) ||
!(clamda = NAG_ALLOC(nctotal, double)) ||
!(h = NAG_ALLOC(n*n, double)) ||
!(x = NAG_ALLOC(n, double)) ||
!(istate = NAG_ALLOC(nctotal, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
pda = n;
pdcj = n;
pdh = n;

/* 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 = 1; i <= n+nclin+ncnln; ++i)
{
scanf("%lf", &bl[i - 1]);
}
scanf("%*[^\n] ");

for (i = 1; i <= n+nclin+ncnln; ++i)
{
scanf("%lf", &bu[i - 1]);
}
scanf("%*[^\n] ");

for (i = 1; i <= n; ++i)
{
scanf("%lf", &x[i - 1]);
}
scanf("%*[^\n] ");

/* nag_opt_nlp_init (e04wcc).
* Initialization function for nag_opt_nlp_solve (e04wdc)
*/
nag_opt_nlp_init(&state, &fail);
if (fail.code != NE_NOERROR)
{
printf("Initialisation of nag_opt_nlp_init (e04wcc) failed.\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* By default nag_opt_nlp_solve (e04wdc) does not print monitoring
* information. Call nag_open_file (x04acc) to set the print file fileid.
*/
/* nag_open_file (x04acc).
* Open unit number for reading, writing or appending, and
* associate unit with named file
*/
nag_open_file("", 2, &fileid, &fail);
if (fail.code != NE_NOERROR)
{
exit_status = 2;
goto END;
}
/* nag_opt_nlp_option_set_integer (e04wgc).
* Set a single option for nag_opt_nlp_solve (e04wdc) from
* an integer argument
*/
fflush(stdout);
nag_opt_nlp_option_set_integer("Print file", fileid, &state, &fail);

/* Solve the problem. */
/* nag_opt_nlp_solve (e04wdc).
* Solves the nonlinear programming (NP) problem
*/
nag_opt_nlp_solve(n, nclin, ncnln, pda, pdcj, pdh, a, bl, bu,
confun, objfun, &majits, istate, ccon, cjac, clamda,
&objf, grad, h, x, &state, &comm, &fail);
fflush(stdout);

if (fail.code == NE_NOERROR)
{
printf("\n\nFinal objective value = %11.3f\n", objf);

printf("Optimal X = ");

for (i = 1; i <= n; ++i)
printf("%9.2f%s", x[i - 1], i%7 == 0 || i == n?"\n":" ");
}
else
{
printf(
"Error message from nag_opt_nlp_solve (e04wdc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

if (fail.code != NE_NOERROR)
exit_status = 2;

}
END:
NAG_FREE(a);
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(ccon);
NAG_FREE(cjac);
NAG_FREE(clamda);
NAG_FREE(h);
NAG_FREE(x);
NAG_FREE(istate);

return exit_status;
}

#undef A

static void NAG_CALL objfun(Integer *mode, Integer n, const double x[],
double *objf, double grad[], Integer nstate,
Nag_Comm *comm)
{
/* Routine to evaluate objective function and its 1st derivatives. */

/* Function Body */
if (comm->user[0] == -1.0)
{
fflush(stdout);
printf("(User-supplied callback objfun, first invocation.)\n");
comm->user[0] = 0.0;
fflush(stdout);
}
if (*mode == 0 || *mode == 2)
{
*objf = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];
}

if (*mode == 1 || *mode == 2)
{
grad[0] = x[3] * (x[0] * 2. + x[1] + x[2]);
grad[2] = x[0] * x[3] + 1.;
grad[3] = x[0] * (x[0] + x[1] + x[2]);
}

return;
} /* objfun */

static void NAG_CALL confun(Integer *mode, Integer ncnln, Integer n,
Integer pdcj, const Integer needc[],
const double x[],
double ccon[], double cjac[], Integer nstate,
Nag_Comm *comm)
{
/* Scalars */
Integer i, j;

#define CJAC(I, J) cjac[(I-1)*pdcj + J-1]

/* Routine to evaluate the nonlinear constraints and their 1st */
/* derivatives. */

/* Function Body */
if (comm->user[1] == -1.0)
{
fflush(stdout);
printf("(User-supplied callback confun, first invocation.)\n");
comm->user[1] = 0.0;
fflush(stdout);
}
if (nstate == 1)
{
/* First call to confun.  Set all 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.;
}
}
}

if (needc[0] > 0)
{
if (*mode == 0 || *mode == 2)
{
ccon[0] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
}
if (*mode == 1 || *mode == 2)
{
CJAC(1, 1) = x[0] * 2.;
CJAC(1, 2) = x[1] * 2.;
CJAC(1, 3) = x[2] * 2.;
CJAC(1, 4) = x[3] * 2.;
}
}

if (needc[1] > 0)
{
if (*mode == 0 || *mode == 2)
{
ccon[1] = x[0] * x[1] * x[2] * x[3];
}
if (*mode == 1 || *mode == 2)
{
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];
}
}

return;
} /* confun */

#undef CJAC