/* nag_opt_handle_set_qconstr_fac (e04rtc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.6, 2022.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL monit(void *handle, const double rinfo[],
const double stats[], Nag_Comm *comm,
Integer *inform);
#ifdef __cplusplus
}
#endif
int main(void) {
Integer n, m, nnza;
Integer nnzu, nnzuc, exit_status, i, x_idx;
Integer idqc, idlc;
Integer verbose_output;
Integer tol_reached;
Integer *icola = 0, *irowa = 0, *idxr0 = 0;
Integer irowlc[3] = {1, 1, 1};
Integer icollc[3] = {1, 2, 3};
double *x = 0, *u = 0, *uc = 0;
double *a = 0, *b = 0, *r0 = 0, *xl = 0, *xu = 0;
double rinfo[100], stats[100];
double lc[3] = {1.0, 1.0, 1.0};
double lc_rhs[1] = {1.0};
double tol_monit;
void *handle = 0;
/* Nag Types */
Nag_Comm comm;
NagError fail;
exit_status = 0;
printf("nag_opt_handle_set_qconstr_fac (e04rtc) Example Program Results\n\n");
fflush(stdout);
/* Read the data file and allocate memory */
scanf(" %*[^\n]"); /* Skip heading in data file */
scanf("%" NAG_IFMT " %" NAG_IFMT " %" NAG_IFMT " %*[^\n]", &n, &m, &nnza);
/* Allocate memory to read data */
if (!(icola = NAG_ALLOC(nnza, Integer)) ||
!(irowa = NAG_ALLOC(nnza, Integer)) || !(idxr0 = NAG_ALLOC(n, Integer)) ||
!(a = NAG_ALLOC(nnza, double)) || !(b = NAG_ALLOC(m, double)) ||
!(r0 = NAG_ALLOC(n, double)) || !(xl = NAG_ALLOC(n, double)) ||
!(xu = NAG_ALLOC(n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read A matrix row indices */
for (i = 0; i < nnza; i++) {
scanf("%" NAG_IFMT, &irowa[i]);
}
scanf("%*[^\n]");
/* Read A matrix column indices */
for (i = 0; i < nnza; i++) {
scanf("%" NAG_IFMT, &icola[i]);
}
scanf("%*[^\n]");
/* Read A values*/
for (i = 0; i < nnza; i++) {
scanf("%lf", &a[i]);
}
scanf("%*[^\n]");
/* Read b values*/
for (i = 0; i < m; i++) {
scanf("%lf", &b[i]);
}
scanf("%*[^\n]");
/* Read xl values*/
for (i = 0; i < n; i++) {
scanf("%lf", &xl[i]);
}
scanf("%*[^\n]");
/* Read xu values*/
for (i = 0; i < n; i++) {
scanf("%lf", &xu[i]);
}
scanf("%*[^\n]");
/* Compute -2*b'A as linear term in quadratic function */
for (i = 0; i < n; i++) {
r0[i] = 0.0;
idxr0[i] = i + 1;
}
for (i = 0; i < nnza; i++) {
r0[icola[i] - 1] = r0[icola[i] - 1] + a[i] * b[irowa[i] - 1];
}
/* Compute size of multipliers */
/* One linear constraint in the model will have
* 2 multipliers for both bounds */
nnzu = 2 * n + 2;
/* No cone constraint in the model, so set nnzuc to 0 */
nnzuc = 0;
/* Allocate memory */
if (!(x = NAG_ALLOC(n, double)) || !(u = NAG_ALLOC(nnzu, double)) ||
!(uc = NAG_ALLOC(nnzuc, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 0; i < n; i++) {
x[i] = 0.0;
}
/* Create the problem handle */
/* nag_opt_handle_init (e04rac).
* Initialize an empty problem handle with n variables. */
nag_opt_handle_init(&handle, n, NAGERR_DEFAULT);
/* nag_opt_handle_set_qconstr_fac (e04rtc)
* Define a quadratic objective */
idqc = -1;
nag_opt_handle_set_qconstr_fac(handle, 0.0, n, idxr0, r0, m, nnza, irowa,
icola, a, &idqc, NAGERR_DEFAULT);
/* nag_opt_handle_set_simplebounds (e04rhc)
* Define bounds on the variables */
nag_opt_handle_set_simplebounds(handle, n, xl, xu, NAGERR_DEFAULT);
/* nag_opt_handle_set_linconstr (e04rjc)
* Define linear constraints: x1 + x2 + x3 = 1 */
idlc = 0;
nag_opt_handle_set_linconstr(handle, 1, lc_rhs, lc_rhs, 3, irowlc, icollc, lc,
&idlc, NAGERR_DEFAULT);
/* nag_opt_handle_opt_set (e04zmc) */
/* Turn on monitoring */
nag_opt_handle_opt_set(handle, "SOCP Monitor Frequency = 1", NAGERR_DEFAULT);
/* Set this to 1 to cause nag_opt_handle_solve_socp_ipm (e04ptc)
* to produce intermediate progress output */
verbose_output = 0;
if (verbose_output) {
/* Require printing of primal and dual solutions at the end of the solve */
nag_opt_handle_opt_set(handle, "Print Solution = Yes", NAGERR_DEFAULT);
} else {
/* Turn off printing of intermediate progress output */
nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
}
tol_reached = 0;
tol_monit = 1.0e-7;
comm.iuser = &tol_reached;
comm.user = &tol_monit;
/* nag_opt_handle_solve_socp_ipm (e04ptc) */
INIT_FAIL(fail);
nag_opt_handle_solve_socp_ipm(handle, n, x, nnzu, u, nnzuc, uc, rinfo, stats,
monit, &comm, &fail);
if (fail.code != NE_NOERROR && fail.code != NW_NOT_CONVERGED) {
printf("nag_opt_handle_solve_socp_ipm (e04ptc) failed.\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Print solution if optimal or suboptimal solution found */
printf(" Optimal X:\n");
printf(" x_idx Value\n");
for (x_idx = 0; x_idx < n; x_idx++) {
printf(" %5" NAG_IFMT " %11.3E\n", x_idx + 1, x[x_idx]);
}
END:
NAG_FREE(x);
NAG_FREE(u);
NAG_FREE(uc);
NAG_FREE(irowa);
NAG_FREE(icola);
NAG_FREE(a);
NAG_FREE(b);
NAG_FREE(r0);
NAG_FREE(idxr0);
NAG_FREE(xl);
NAG_FREE(xu);
/* nag_opt_handle_free (e04rzc).
* Destroy the problem handle and deallocate all the memory. */
if (handle)
nag_opt_handle_free(&handle, NAGERR_DEFAULT);
return exit_status;
}
static void NAG_CALL monit(void *handle, const double rinfo[],
const double stats[], Nag_Comm *comm,
Integer *inform) {
/* Monitoring function can be used to monitor the progress
* or, for example, to implement bespoke stopping criteria */
double tol = comm->user[0];
Integer *tol_reached = comm->iuser;
/* If x is close to the solution, print a message */
if (rinfo[14] < tol && rinfo[15] < tol && rinfo[16] < tol &&
rinfo[17] < tol) {
if (!*tol_reached) {
printf("\n monit() reports good approximate solution "
"(tol = %8.2E)\n",
tol);
*tol_reached = 1;
}
}
fflush(stdout);
}