/* nag_opt_handle_set_simplebounds (e04rhc) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.3, 2023.
*/
/* Matrix completion problem (rank minimization) solved approximatelly
* by SDP via nuclear norm minimization formulated as:
* min trace(X1) + trace(X2)
* s.t. [ X1, Y; Y', X2 ] >=0
* 0 <= Y_ij <= 1
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
int main(void) {
double const stol = 1e-5;
double const time_limit = 180.0;
Integer exit_status = 0;
Integer dima, i, idblk, idx, idxend, idxobj, idxx, inform, j, lwork, m, maxs,
n, nblk, nnz, nnzasum, nnzc, nnzu, nnzua, nnzuc, nvar, rank;
double rdummy[1], rinfo[32], stats[32];
double *a = 0, *bl = 0, *bu = 0, *c = 0, *s = 0, *work = 0, *x = 0, *y = 0;
Integer *icola = 0, *idxc = 0, *irowa = 0, *nnza = 0;
void *handle = 0;
/* Nag Types */
NagError fail;
#define Y(I, J) y[(J - 1) * m + I - 1]
printf("nag_opt_handle_set_simplebounds (e04rhc) Example Program Results"
"\n\n");
fflush(stdout);
/* Skip heading in data file. */
scanf("%*[^\n]");
/* Read in the problem size and ignore the rest of the line. */
scanf("%" NAG_IFMT " %" NAG_IFMT "%*[^\n]", &m, &n);
/* Allocate memory for matrix Y and read it in. */
if (!(y = NAG_ALLOC(m * n, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
for (i = 1; i <= m; i++)
for (j = 1; j <= n; j++)
scanf("%lf", &Y(i, j));
scanf("%*[^\n]");
/* Count the number of specified elements (i.e., nonnegative) */
nnz = 0;
for (i = 1; i <= m; i++)
for (j = 1; j <= n; j++)
if (Y(i, j) >= 0.0)
nnz++;
/* There are as many variables as missing entries in the Y matrix
* plus two full symmetric matrices m x m and n x n. */
nvar = m * (m + 1) / 2 + n * (n + 1) / 2 + m * n - nnz;
if (!(x = NAG_ALLOC(nvar, double)) || !(bl = NAG_ALLOC(nvar, double)) ||
!(bu = NAG_ALLOC(nvar, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* nag_opt_handle_init (e04rac).
* Initialize an empty problem handle with NVAR variables. */
nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);
/* Create bounds for the missing entries in Y matrix to be between 0 and 1 */
idxend = m * (m + 1) / 2 + n * (n + 1) / 2;
for (idx = 0; idx < idxend; idx++) {
bl[idx] = -1e+20;
bu[idx] = 1e+20;
}
for (; idx < nvar; idx++) {
bl[idx] = 0.0;
bu[idx] = 1.0;
}
/* nag_opt_handle_set_simplebounds (e04rhc).
* Define bounds on the variables. */
nag_opt_handle_set_simplebounds(handle, nvar, bl, bu, NAGERR_DEFAULT);
/* Allocate space for the objective - minimize trace of the matrix
* constraint. There is no quadratic part in the objective. */
nnzc = m + n;
if (!(idxc = NAG_ALLOC(nnzc, Integer)) || !(c = NAG_ALLOC(nnzc, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Construct linear matrix inequality to request that
* [ X1, Y; Y', X2] is positive semidefinite. */
/* How many nonzeros do we need? As many as number of variables
* and the number of specified elements together. */
nnzasum = m * (m + 1) / 2 + n * (n + 1) / 2 + m * n;
if (!(nnza = NAG_ALLOC(nvar + 1, Integer)) ||
!(irowa = NAG_ALLOC(nnzasum, Integer)) ||
!(icola = NAG_ALLOC(nnzasum, Integer)) ||
!(a = NAG_ALLOC(nnzasum, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
nnza[0] = nnz;
for (i = 1; i <= nvar; i++)
nnza[i] = 1;
/* Copy Y to the upper block of A_0 with different sign (because -A_0)
* (upper triangle) */
idx = 0;
for (i = 1; i <= m; i++)
for (j = 1; j <= n; j++)
if (Y(i, j) >= 0.0) {
irowa[idx] = i;
icola[idx] = m + j;
a[idx] = -Y(i, j);
idx++;
}
/* One matrix for each variable, A_i has just one nonzero - it binds
* x_i with its position in the matrix constraint. Set also the objective. */
/* 1,1 - block, X1 matrix (mxm) */
idxobj = 0;
idxx = 1;
for (i = 1; i <= m; i++) {
/* the next element is diagonal ==> part of the objective as a trace() */
idxc[idxobj] = idxx;
c[idxobj] = 1.0;
idxobj++;
for (j = i; j <= m; j++) {
irowa[idx] = i;
icola[idx] = j;
a[idx] = 1.0;
idx++;
idxx++;
}
}
/* 2,2 - block, X2 matrix (nxn) */
for (i = 1; i <= n; i++) {
/* the next element is diagonal ==> part of the objective as a trace() */
idxc[idxobj] = idxx;
c[idxobj] = 1.0;
idxobj++;
for (j = i; j <= n; j++) {
irowa[idx] = m + i;
icola[idx] = m + j;
a[idx] = 1.0;
idx++;
idxx++;
}
}
/* 1,2 - block, missing element in Y we are after */
for (i = 1; i <= m; i++)
for (j = 1; j <= n; j++)
if (Y(i, j) < 0.0) {
irowa[idx] = i;
icola[idx] = m + j;
a[idx] = 1.0;
idx++;
}
/* nag_opt_handle_set_quadobj (e04rfc).
* Add the sparse linear objective to the handle.*/
nag_opt_handle_set_quadobj(handle, nnzc, idxc, c, 0, NULL, NULL, NULL,
NAGERR_DEFAULT);
/* Just one matrix inequality of the dimension of the extended matrix. */
nblk = 1;
dima = m + n;
idblk = 0;
/* nag_opt_handle_set_linconstr (e04rnc).
* Add the linear matrix constraint to the problem formulation. */
nag_opt_handle_set_linmatineq(handle, nvar, dima, nnza, nnzasum, irowa, icola,
a, nblk, NULL, &idblk, NAGERR_DEFAULT);
/* nag_opt_handle_opt_set (e04zmc).
* Set optional arguments of the solver:
* Completely turn off printing, allow timing and
* turn on the monitor mode to stop every iteration. */
nag_opt_handle_opt_set(handle, "Print File = -1", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Stats Time = Yes", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Monitor Frequency = 1", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Initial X = Automatic", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Dimacs = Check", NAGERR_DEFAULT);
/* Pass the handle to the solver, we are not interested in
* Lagrangian multipliers. */
nnzu = 0;
nnzuc = 0;
nnzua = 0;
while (1) {
INIT_FAIL(fail);
/* nag_opt_handle_solve_pennon (e04svc). */
nag_opt_handle_solve_pennon(handle, nvar, x, nnzu, NULL, nnzuc, NULL, nnzua,
NULL, rinfo, stats, &inform, &fail);
if (inform == 0 || fail.code != NE_NOERROR) {
/* Final exit, solver finished. */
printf("Finished at iteration %2" NAG_IFMT
": objective %7.2f, avg.error %9.2e\n",
(Integer)stats[0], rinfo[0],
(rinfo[1] + rinfo[2] + rinfo[3]) / 3.0);
fflush(stdout);
break;
} else {
/* Monitor stop */
printf("Monitor at iteration %2" NAG_IFMT
": objective %7.2f, avg.error %9.2e\n",
(Integer)stats[0], rinfo[0],
(rinfo[1] + rinfo[2] + rinfo[3]) / 3.0);
fflush(stdout);
/* Check time limit and possibly stop the solver. */
if (stats[7] > time_limit)
inform = 0;
}
}
if (fail.code == NE_NOERROR || fail.code == NW_NOT_CONVERGED) {
/* Successful run, fill the missing elements in the matrix Y. */
idx = m * (m + 1) / 2 + n * (n + 1) / 2;
for (i = 1; i <= m; i++)
for (j = 1; j <= n; j++)
if (Y(i, j) < 0.0)
Y(i, j) = x[idx++];
/* nag_file_print_matrix_real_gen_comp (x04cbc).
* Print the matrix. */
nag_file_print_matrix_real_gen_comp(
Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, m, n, y, m, "%7.1f",
"Completed Matrix", Nag_IntegerLabels, NULL, Nag_IntegerLabels, NULL,
80, 0, NULL, NAGERR_DEFAULT);
/* Compute rank of the matrix via SVD, use the fact that the order
* of the singular values is descending. */
lwork = 20 * (m > n ? m : n);
maxs = m < n ? m : n;
if (!(s = NAG_ALLOC(maxs, double)) ||
!(work = NAG_ALLOC((lwork), double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* nag_lapackeig_dgesvd (f08kbc).
* Compute the singular values */
nag_lapackeig_dgesvd(Nag_ColMajor, Nag_NotU, Nag_NotVT, m, n, y, m, s,
rdummy, 1, rdummy, 1, work, NAGERR_DEFAULT);
for (rank = 0; rank < maxs; rank++)
if (s[rank] <= stol)
break;
printf("Rank is %" NAG_IFMT "\n", rank);
} else if (fail.code == NE_USER_STOP) {
printf("The given time limit was reached, run aborted.\n");
} else {
printf("Error from nag_opt_handle_solve_pennon (e04svc).\n%s\n",
fail.message);
exit_status = 1;
}
END:
/* nag_opt_handle_free (e04rzc).
* Destroy the problem handle and deallocate all the memory. */
if (handle)
nag_opt_handle_free(&handle, NAGERR_DEFAULT);
NAG_FREE(a);
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(c);
NAG_FREE(s);
NAG_FREE(work);
NAG_FREE(x);
NAG_FREE(y);
NAG_FREE(icola);
NAG_FREE(irowa);
NAG_FREE(nnza);
NAG_FREE(idxc);
return exit_status;
}