/* nag_opt_qp_dense_solve (e04nfc) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.0, 2023.
*
*/
#include <nag.h>
#include <stdio.h>
#include <string.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[2] = {-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_dense_solve (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));
/* Read the bounds. */
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;
/* nag_opt_read (e04xyc).
* Read options from a text file
*/
nag_opt_read("e04nfc", optionsfile, &options, print, "stdout", &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_opt_read (e04xyc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Solve the problem from a cold start.
* The Hessian is defined implicitly by function qphess1.
*/
/* nag_opt_qp_dense_solve (e04nfc), see above. */
nag_opt_qp_dense_solve(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_dense_solve (e04nfc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* The following is for illustrative purposes only. We do a warm
* start with the final working set of the previous run.
* 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_dense_solve (e04nfc), see above. */
nag_opt_qp_dense_solve(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_dense_solve (e04nfc).\n%s\n", fail.message);
exit_status = 1;
}
/* Free memory allocated by nag_opt_qp_dense_solve (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[0] == -1.0) {
printf("(User-supplied callback qphess1, first invocation.)\n");
fflush(stdout);
comm->user[0] = 0.0;
}
hx[0] = 2.0 * x[0];
hx[1] = 2.0 * x[1];
hx[2] = 2.0 * (x[2] + x[3]);
hx[3] = hx[2];
hx[4] = 2.0 * x[4];
hx[5] = -2.0 * (x[5] + x[6]);
hx[6] = hx[5];
} /* 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] == -1.0) {
printf("(User-supplied callback qphess2, first invocation.)\n");
fflush(stdout);
comm->user[1] = 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 */