/* nag_opt_qp (e04nfc) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*
*/

#include <nag.h>
#include <stdio.h>
#include <string.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL qphess1(Integer n, Integer jthcol, const double h[],
Integer tdh, const double x[], double hx[],
Nag_Comm *comm);
static void NAG_CALL qphess2(Integer n, Integer jthcol, const double h[],
Integer tdh, const double x[], double hx[],
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

#define A(I, J) a[(I) *tda + J]
#define H(I, J) h[(I) *tdh + J]

int main(void)
{
const char *optionsfile = "e04nfce.opt";
static double ruser = { -1.0, -1.0 };
Nag_Boolean print;
Integer exit_status = 0, i, j, n, nbnd, nclin, tda, tdh;
Nag_E04_Opt options;
double *a = 0, *bl = 0, *bu = 0, *cvec = 0, *h = 0, objf, *x = 0;
Nag_Comm comm;
NagError fail;

INIT_FAIL(fail);

printf("nag_opt_qp (e04nfc) Example Program Results\n");

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

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

/* Set the actual problem dimensions.
* n = the number of variables.
* nclin = the number of general linear constraints (may be 0).
*/
n = 7;
nclin = 7;
if (n > 0 && nclin >= 0) {
nbnd = n + nclin;
if (!(x = NAG_ALLOC(n, double)) ||
!(cvec = NAG_ALLOC(n, double)) ||
!(a = NAG_ALLOC(nclin * n, double)) ||
!(h = NAG_ALLOC(n * n, double)) ||
!(bl = NAG_ALLOC(nbnd, double)) || !(bu = NAG_ALLOC(nbnd, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
tda = n;
tdh = n;
}
else {
printf("Invalid n or nclin.\n");
exit_status = 1;
return exit_status;
}
/* cvec   = the coefficients of the explicit linear term of f(x).
* a      = the linear constraint matrix.
* bl     = the lower bounds on x and A*x.
* bu     = the upper bounds on x and A*x.
* x      = the initial estimate of the solution.
*/

/* Read the coefficients of the explicit linear term of f(x). */
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < n; ++i)
scanf("%lf", &cvec[i]);

/* Read the linear constraint matrix A. */
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < nclin; ++i)
for (j = 0; j < n; ++j)
scanf("%lf", &A(i, j));

nbnd = n + nclin;
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < nbnd; ++i)
scanf("%lf", &bl[i]);
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < nbnd; ++i)
scanf("%lf", &bu[i]);

/* Read the initial estimate of x. */
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < n; ++i)
scanf("%lf", &x[i]);

/* nag_opt_init (e04xxc).
* Initialization function for option setting
*/
nag_opt_init(&options); /* Initialize options structure */
/* Set one option directly
* Bounds  >=   inf_bound will be treated as plus  infinity.
* Bounds  <=  -inf_bound will be treated as minus infinity.
*/
options.inf_bound = 1.0e21;

/* Read remaining option values from file */
print = Nag_TRUE;
* Read options from a text file
*/
nag_opt_read("e04nfc", optionsfile, &options, print, "stdout", &fail);
if (fail.code != NE_NOERROR) {
exit_status = 1;
goto END;
}

/* Solve the problem from a cold start.
* The Hessian is defined implicitly by function qphess1.
*/
/* nag_opt_qp (e04nfc), see above. */
nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *) 0, tdh,
qphess1, x, &objf, &options, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* The following is for illustrative purposes only. We do a warm
* This time we store the Hessian explicitly in h[][], and use
* the corresponding function qphess2().
* Only the final solution from the results is printed.
*/
printf("\nA run of the same example with a warm start:\n");
fflush(stdout);

options.start = Nag_Warm;
options.print_level = Nag_Soln;

for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j)
H(i, j) = 0.0;
if (i <= 4)
H(i, i) = 2.0;
else
H(i, i) = -2.0;
}
H(2, 3) = 2.0;
H(3, 2) = 2.0;
H(5, 6) = -2.0;
H(6, 5) = -2.0;

/* Solve the problem again. */
/* nag_opt_qp (e04nfc), see above. */
nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh,
qphess2, x, &objf, &options, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message);
exit_status = 1;
}
/* Free memory allocated by nag_opt_qp (e04nfc) to pointers in options */
/* nag_opt_free (e04xzc).
* Memory freeing function for use with option setting
*/
nag_opt_free(&options, "all", &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

END:
NAG_FREE(x);
NAG_FREE(cvec);
NAG_FREE(a);
NAG_FREE(h);
NAG_FREE(bl);
NAG_FREE(bu);

return exit_status;
}

static void NAG_CALL qphess1(Integer n, Integer jthcol, const double h[],
Integer tdh, const double x[], double hx[],
Nag_Comm *comm)
{
/* In this version of qphess the Hessian matrix is implicit.
* The array h[] is not accessed. There is no special coding
* for the case jthcol > 0.
*/

if (comm->user == -1.0) {
printf("(User-supplied callback qphess1, first invocation.)\n");
fflush(stdout);
comm->user = 0.0;
}

hx = 2.0 * x;
hx = 2.0 * x;
hx = 2.0 * (x + x);
hx = hx;
hx = 2.0 * x;
hx = -2.0 * (x + x);
hx = hx;
} /* qphess1 */

#undef H

static void NAG_CALL qphess2(Integer n, Integer jthcol, const double h[],
Integer tdh, const double x[], double hx[],
Nag_Comm *comm)
{
/* In this version of qphess, the matrix H is stored in h[]
* as a full two-dimensional array.
*/

#define H(I, J) h[(I) *tdh + (J)]

Integer i, j;

if (comm->user == -1.0) {
printf("(User-supplied callback qphess2, first invocation.)\n");
fflush(stdout);
comm->user = 0.0;
}

if (jthcol != 0) {
/* Special case -- extract one column of  H. */
j = jthcol - 1;
for (i = 0; i < n; ++i)
hx[i] = H(i, j);
}
else {
/* Normal Case. */
for (i = 0; i < n; ++i)
hx[i] = 0.0;

for (i = 0; i < n; ++i)
for (j = 0; j < n; ++j)
hx[i] += H(i, j) * x[j];
}
} /* qphess2 */