/* nag_matop_real_nmf_rcomm (f01sbc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.5, 2022.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#define A_DENSE(I, J) a_dense[(J - 1) * pda + I - 1]
int main(void) {
Integer exit_status = 0, i, j, irevcm = 0, m, n, k, icount = 0, maxit,
exit_loop = 0, pdw, pdh, pdht, pda, q, seed, nnz;
double *comm = 0, *a_dense = 0, *a = 0, *w = 0, *h = 0, *ht = 0;
Integer icomm[9], *icolzp = 0, *irowix = 0, element;
double errtol, norm, norma;
/* Nag Types */
Nag_OrderType order = Nag_ColMajor;
NagError fail;
INIT_FAIL(fail);
printf("nag_matop_real_nmf_rcomm (f01sbc) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
/* Read values of m, n and k */
scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &m, &n, &k);
/* Pretend A is a square q x q matrix to allow use of
* nag_sparse_direct_real_gen_matmul for matrix multiplications
*/
q = MAX(m, n);
/* Allocate memory */
pdw = q;
pdh = k;
pdht = q;
pda = m;
if (!(w = NAG_ALLOC(pdw * k, double)) || !(h = NAG_ALLOC(pdh * n, double)) ||
!(ht = NAG_ALLOC(pdht * k, double)) ||
!(a_dense = NAG_ALLOC(pda * n, double)) ||
!(icolzp = NAG_ALLOC(q + 1, Integer)) ||
!(comm = NAG_ALLOC((2 * m + n) * k + 3, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read the matrix A */
for (i = 0; i < n + 1; ++i)
scanf("%" NAG_IFMT "%*[^\n] ", &icolzp[i]);
nnz = icolzp[n] - 1;
if (!(a = NAG_ALLOC(nnz, double)) || !(irowix = NAG_ALLOC(nnz, Integer))) {
printf("Allocation failure\n");
exit_status = -2;
goto END;
}
for (i = 0; i < nnz; ++i)
scanf("%lf%" NAG_IFMT "%*[^\n] ", &a[i], &irowix[i]);
/* If m > n (q==m), nag_sparse_direct_real_gen_matmul will need
* a total of q columns, otherwise no adjustment is necessary
*/
for (i = n + 1; i < q + 1; ++i)
icolzp[i] = icolzp[n];
/* Initialize w and ht to zero */
for (i = 0; i < pdw * k; ++i)
w[i] = 0.0;
for (i = 0; i < pdht * k; ++i)
ht[i] = 0.0;
/* Choose the values of seed and errtol */
seed = 234;
errtol = 1.0e-3;
/* We will terminate after 200 iterations */
maxit = 200;
/* nag_matop_real_nmf_rcomm (f01sbc).
* Non-negative matrix factorisation
* (reverse communication)
*/
do {
nag_matop_real_nmf_rcomm(&irevcm, m, n, k, w, pdw, h, pdh, ht, pdht, seed,
errtol, comm, icomm, &fail);
switch (irevcm) {
case 1:
if (icount == maxit) {
printf("Exiting after the maximum number of iterations\n\n");
exit_loop = 1;
} else {
icount++;
}
/* w and h are available for printing */
break;
case 2:
/* compute a^t * w and store in ht */
nag_sparse_direct_real_gen_matmul(order, Nag_Trans, q, k, 1.0, icolzp,
irowix, a, w, pdw, 0.0, ht, pdht,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparse_direct_real_gen_matmul (f11mkc) \
\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
break;
case 3:
/* compute a * ht and store in w */
nag_sparse_direct_real_gen_matmul(order, Nag_NoTrans, q, k, 1.0, icolzp,
irowix, a, ht, pdht, 0.0, w, pdw,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparse_direct_real_gen_matmul (f11mkc) \
\n%s\n",
fail.message);
exit_status = 2;
goto END;
}
break;
}
} while (irevcm != 0 && exit_loop == 0);
if (fail.code != NE_NOERROR) {
printf("Error from nag_matop_real_nmf_rcomm (f01sbc).\n%s\n", fail.message);
exit_status = 3;
goto END;
}
/* Print matrix W using nag_file_print_matrix_real_gen (x04cac)
* Print real general matrix (easy-to-use)
*/
nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m,
k, w, pdw, "W", NULL, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen (x04cac)\n%s\n",
fail.message);
exit_status = 4;
goto END;
}
printf("\n");
/* Print matrix H using nag_file_print_matrix_real_gen (x04cac)
* Print real general matrix (easy-to-use)
*/
nag_file_print_matrix_real_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag, k,
n, h, pdh, "H", NULL, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_gen (x04cac)\n%s\n",
fail.message);
exit_status = 5;
goto END;
}
/* In order to compute the residual we convert a to dense format */
for (i = 1; i <= n; i++) {
for (j = 1; j <= m; j++) {
A_DENSE(j, i) = 0.0;
}
}
element = 0;
for (i = 1; i <= n; i++) {
for (j = 1; j <= icolzp[i] - icolzp[i - 1]; j++) {
A_DENSE(irowix[element], i) = a[element];
element++;
}
}
/* Compute ||A||_F using
* nag_blast_dge_norm (f16rac)
*/
nag_blast_dge_norm(order, Nag_FrobeniusNorm, m, n, a_dense, pda, &norma,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_blast_dge_norm (f16yac)\n%s\n", fail.message);
exit_status = 6;
goto END;
}
/* Compute residual A-WH using nag_blast_dgemm (f16yac)
* Matrix-matrix product, two rectangular matrices
*/
nag_blast_dgemm(order, Nag_NoTrans, Nag_NoTrans, m, n, k, -1.0, w, pdw, h,
pdh, 1.0, a_dense, pda, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_blast_dgemm (f16yac)\n%s\n", fail.message);
exit_status = 7;
goto END;
}
/* Compute ||A-WH||_F using
* nag_blast_dge_norm (f16rac)
*/
nag_blast_dge_norm(order, Nag_FrobeniusNorm, m, n, a_dense, pda, &norm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_blast_dge_norm (f16yac)\n%s\n", fail.message);
exit_status = 8;
goto END;
}
/* Print relative residual norm */
printf("\nThe relative residual norm, ||A-WH||/||A||, is: %9.2e\n",
norm / norma);
END:
NAG_FREE(a);
NAG_FREE(a_dense);
NAG_FREE(icolzp);
NAG_FREE(irowix);
NAG_FREE(comm);
NAG_FREE(w);
NAG_FREE(h);
NAG_FREE(ht);
return exit_status;
}