/* nag_opt_miqp_mps_read (e04mxc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.0, 2024.
*/
#include <nag.h>
#include <stdio.h>
#include <string.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL qphx(Integer ncolh, const double x[], double hx[],
Integer nstate, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
/* Make a typedef for convenience when allocating crname. */
typedef char e04mx_name[9];
int main(void) {
Integer exit_status = 0;
double obj, objadd, sinf;
Integer i, iobj, j, lenc, lintvar, m, maxlintvar, maxm, maxn, maxncolh,
maxnnz, maxnnzh, minmax, n, ncolh, ninf, nname, nnz, nnzh, ns;
Integer mpslst = 1, verbose_output;
double *a = 0, *bl = 0, *bu = 0, *h = 0, *pi = 0, *rc = 0, *ruser = 0, *x = 0;
Integer *helast = 0, *hs = 0, *iccola = 0, *iccolh = 0, *intvar = 0,
*irowa = 0, *irowh = 0, *iuser = 0;
char pnames[5][9] = {"", "", "", "", ""};
char(*crname)[9] = 0;
char **names = 0;
char fname[] = "e04mxce.opt";
/* Nag Types */
Nag_Boolean readints = Nag_FALSE;
Nag_E04State state;
NagError fail;
Nag_Comm comm;
Nag_FileID fileid;
INIT_FAIL(fail);
printf("nag_opt_miqp_mps_read (e04mxc) Example Program Results\n");
fflush(stdout);
/* nag_file_open (x04acc).
Open unit number for reading and associate unit with named file. */
nag_file_open(fname, 0, &fileid, NAGERR_DEFAULT);
/* nag_opt_miqp_mps_read (e04mxc).
Reads MPS data file defining LP, QP, MILP or MIQP problem.
Query call. */
nag_opt_miqp_mps_read(fileid, 0, 0, 0, 0, 0, 0, mpslst, &n, &m, &nnz, &ncolh,
&nnzh, &lintvar, NULL, NULL, NULL, NULL, NULL, NULL,
NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
NAGERR_DEFAULT);
/* nag_file_close (x04adc).
Close file associated with given unit number. */
nag_file_close(fileid, NAGERR_DEFAULT);
maxm = m;
maxn = n;
maxnnz = nnz;
maxnnzh = nnzh;
maxncolh = ncolh;
maxlintvar = (readints && lintvar > 0) ? lintvar : -1;
if (!(irowa = NAG_ALLOC(maxnnz, Integer)) ||
!(iccola = NAG_ALLOC(maxn + 1, Integer)) ||
!(a = NAG_ALLOC(maxnnz, double)) ||
!(bl = NAG_ALLOC(maxn + maxm, double)) ||
!(bu = NAG_ALLOC(maxn + maxm, double)) ||
!(irowh = NAG_ALLOC(maxnnzh, Integer)) ||
!(iccolh = NAG_ALLOC(maxncolh + 1, Integer)) ||
!(h = NAG_ALLOC(maxnnzh, double)) ||
(maxlintvar > 0 && !(intvar = NAG_ALLOC(maxlintvar, Integer))) ||
!(crname = NAG_ALLOC(maxn + maxm, e04mx_name))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
nag_file_open(fname, 0, &fileid, NAGERR_DEFAULT);
/* Full call to the reader. */
nag_opt_miqp_mps_read(fileid, maxn, maxm, maxnnz, maxncolh, maxnnzh,
maxlintvar, mpslst, &n, &m, &nnz, &ncolh, &nnzh,
&lintvar, &iobj, a, irowa, iccola, bl, bu, pnames,
&nname, crname, h, irowh, iccolh, &minmax,
(maxlintvar > 0) ? intvar : NULL, NAGERR_DEFAULT);
nag_file_close(fileid, NAGERR_DEFAULT);
/* Data has been read. Set up and run the solver.
We have no explicit objective vector so set lenc = 0; the
objective vector is stored in row iobj of a. */
lenc = 0;
objadd = 0.0;
if (!(helast = NAG_ALLOC(n + m, Integer)) ||
!(x = NAG_ALLOC(n + m, double)) || !(pi = NAG_ALLOC(m, double)) ||
!(rc = NAG_ALLOC(n + m, double)) || !(hs = NAG_ALLOC(n + m, Integer)) ||
!(iuser = NAG_ALLOC(ncolh + 1 + nnzh, Integer)) ||
!(ruser = NAG_ALLOC(nnzh, double)) ||
!(names = NAG_ALLOC(n + m, char *))) {
printf("Allocation failure\n");
exit_status = -3;
goto END;
}
/* Transform char (*crname)[9] to char *names[] to be compatible with the
solver. */
for (i = 0; i < n + m; i++)
names[i] = crname[i];
for (i = 0; i < n + m; i++)
helast[i] = 0;
for (i = 0; i < n + m; i++)
hs[i] = 0;
for (i = 0; i < n + m; i++)
x[i] = MIN(MAX(0.0, bl[i]), bu[i]);
if (ncolh > 0) {
/* Store the nonzeros of H in ruser for use by qphx. */
for (i = 0; i < nnzh; i++)
ruser[i] = h[i];
/* Store iccolh and irowh in iuser for use by qphx. */
for (i = 0; i < ncolh + 1; i++)
iuser[i] = iccolh[i];
for (i = ncolh + 1, j = 0; i < nnzh + ncolh + 1; i++, j++)
iuser[i] = irowh[j];
comm.iuser = iuser;
comm.user = ruser;
}
/* nag_opt_qpconvex2_sparse_init (e04npc).
Initialization function for nag_opt_qpconvex2_sparse_solve (e04nqc). */
nag_opt_qpconvex2_sparse_init(&state, NAGERR_DEFAULT);
/* Use nag_opt_qpconvex2_sparse_option_string (e04nsc) to change
the direction of optimization. Minimization is assumed by default. */
if (minmax == 1)
nag_opt_qpconvex2_sparse_option_string("Maximize", &state, NAGERR_DEFAULT);
else if (minmax == 0)
nag_opt_qpconvex2_sparse_option_string("Feasible Point", &state,
NAGERR_DEFAULT);
/* Set this to 1 to cause intermediate progress output to be produced. */
verbose_output = 0;
if (verbose_output == 1) {
/* By default nag_opt_qpconvex2_sparse_solve (e04nqc) does not print
* monitoring information. Call nag_file_open (x04acc) to set the print
* file fileid and then call
* nag_opt_qpconvex2_sparse_option_integer_set (e04ntc) to register that
* setting with the solver.
*/
nag_file_open("", 2, &fileid, NAGERR_DEFAULT);
nag_opt_qpconvex2_sparse_option_integer_set("Print file", fileid, &state,
NAGERR_DEFAULT);
} else {
printf("\nProblem contains %3" NAG_IFMT " variables and %3" NAG_IFMT
" linear constraints\n",
n, m);
}
/* nag_opt_qpconvex2_sparse_solve (e04nqc).
LP or QP problem (suitable for sparse problems). */
nag_opt_qpconvex2_sparse_solve(
Nag_Cold, qphx, m, n, nnz, nname, lenc, ncolh, iobj, objadd, pnames[0], a,
irowa, iccola, bl, bu, NULL, (const char **)names, helast, hs, x, pi, rc,
&ns, &ninf, &sinf, &obj, &state, &comm, &fail);
if (verbose_output == 0) {
printf("\n");
printf("Final objective value = %12.3e\n", obj);
printf("Optimal x = ");
for (i = 0; i < n; ++i) {
printf("%9.2f%s", x[i],
i % 6 == 5 || i == n - 1 ? "\n " : " ");
}
}
END:
NAG_FREE(a);
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(h);
NAG_FREE(pi);
NAG_FREE(rc);
NAG_FREE(ruser);
NAG_FREE(x);
NAG_FREE(helast);
NAG_FREE(hs);
NAG_FREE(iccola);
NAG_FREE(iccolh);
NAG_FREE(intvar);
NAG_FREE(irowa);
NAG_FREE(irowh);
NAG_FREE(iuser);
NAG_FREE(crname);
NAG_FREE(names);
return exit_status;
}
static void NAG_CALL qphx(Integer ncolh, const double x[], double hx[],
Integer nstate, Nag_Comm *comm) {
/* Function to compute H*x.
Note: comm->iuser and comm->user contain the following data:
comm->user[0:nnzh-1] = h[0:nnzh-1]
comm->iuser[0:ncolh] = iccolh[0:ncolh]
comm->iuser[ncolh+1:nnzh+ncolh] = irowh[0:nnzh-1] */
Integer i, end, icol, idx, irow, start;
for (i = 0; i < ncolh; i++)
hx[i] = 0.0;
for (icol = 0; icol < ncolh; icol++) {
start = comm->iuser[icol];
end = comm->iuser[icol + 1] - 1;
for (idx = start - 1; idx < end; idx++) {
irow = comm->iuser[ncolh + 1 + idx] - 1;
hx[irow] = hx[irow] + x[icol] * comm->user[idx];
if (irow != icol)
hx[icol] = hx[icol] + x[irow] * comm->user[idx];
}
}
}