NAG Library Manual, Mark 29.2
```// example of using:
//   nagcpp::matop::real_nmf_rcomm (f01sb)
//   nagcpp::sparse::real_gen_matvec (f11xa)
//   nagcpp::blas::dlange (f06ra)
//   nagcpp::blas::dgemm (f06ya)
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <vector>

#include "f01/nagcpp_f01sb.hpp"
#include "f06/nagcpp_f06ra.hpp"
#include "f06/nagcpp_f06ya.hpp"
#include "f11/nagcpp_f11xa.hpp"

// a simple sample matrix class
#include "include/nag_my_matrix.hpp"

int main(void) {
std::cout << "nagcpp::matop::real_nmf_rcomm Example" << std::endl;
typedef std::vector<double>::size_type vsize_t;
std::vector<double> a = {6.0, 3.0, 6.0, 3.0, 1.0, 7.0, 2.0, 2.0, 4.0,
2.0, 2.0, 6.0, 3.0, 2.0, 1.0, 2.0, 1.0, 6.0,
7.0, 3.0, 1.0, 2.0, 1.0, 1.0, 3.0};
std::vector<nagcpp::types::f77_integer> icolix = {1, 1, 1, 1, 2, 2, 2, 2, 3,
3, 3, 3, 3, 4, 4, 4, 4, 5,
5, 5, 5, 6, 6, 6, 6};
std::vector<nagcpp::types::f77_integer> irowix = {2, 4, 5, 6, 1, 3, 5, 7, 1,
3, 4, 6, 7, 2, 4, 5, 6, 3,
4, 5, 7, 1, 3, 4, 6};
size_t m = 7;
size_t n = 6;
size_t k = 4;
MyMatrix<double> w(m, k);
MyMatrix<double> h(k, n);
MyMatrix<double> ht(n, k);

// choose values for the arguments errtol and seed
// rather than using the defaults
nagcpp::matop::OptionalF01SB opt;
opt.seed(234);
opt.errtol(1.0e-3);

// we will terminate after 200 iterations
vsize_t maxit = 200;

// we are going to pretend that A is q x q, rather than m x n
size_t q = std::max(m, n);
std::vector<double> tmp_x(q, 0);
std::vector<double> tmp_y(q);

vsize_t icount;
double rel_resid_norm;

// find a non-negative matrix factorixation A ~= WH
try {
nagcpp::utility::CopyableComm comm;
nagcpp::types::f77_integer irevcm = 0;
for (icount = 0; icount < maxit;) {
nagcpp::matop::real_nmf_rcomm(irevcm, w, h, ht, comm, opt);

if (irevcm == 1) {
++icount;

} else if (irevcm == 2) {
// compute a^t * w and store in ht
for (vsize_t j = 0; j < k; ++j) {
for (vsize_t i = 0; i < m; ++i) {
tmp_x[i] = w(i, j);
}
// sparse matrix vector product (f11xa)
nagcpp::sparse::real_gen_matvec("Transpose", a, irowix, icolix, tmp_x,
tmp_y);
for (vsize_t i = 0; i < n; ++i) {
ht(i, j) = tmp_y[i];
}
}

} else if (irevcm == 3) {
// compute a * ht and store in w (f11mk)
for (vsize_t j = 0; j < k; ++j) {
for (vsize_t i = 0; i < n; ++i) {
tmp_x[i] = ht(i, j);
}
// sparse matrix vector product (f11xa)
nagcpp::sparse::real_gen_matvec("NoTranspose", a, irowix, icolix,
tmp_x, tmp_y);
for (vsize_t i = 0; i < m; ++i) {
w(i, j) = tmp_y[i];
}
}

} else {
// nagcpp::matop::real_nmf_rcomm has finished
break;
}
}

std::cout << "Exited after ";
if (icount >= maxit) {
std::cout << "the maximum number of iterations (" << maxit << ")"
<< std::endl;
} else {
std::cout << icount << " iterations" << std::endl;
}

// get a dense version of the matrix A
MyMatrix<double> da(m, n, std::vector<double>(m * n, 0));
for (vsize_t i = 0; i < a.size(); ++i) {
da(irowix[i] - 1, icolix[i] - 1) = a[i];
}

// get ||A|| (f06ra)
double norma = nagcpp::blas::dlange("F", da);

// calculate A = A-WH (f06ya)
nagcpp::blas::dgemm("NoTranspose", "NoTranspose", -1.0, w, h, 1.0, da);

// get ||A-WH|| (f06ra)
double norm = nagcpp::blas::dlange("F", da);
rel_resid_norm = norm / norma;

} catch (nagcpp::error_handler::Exception &e) {
std::cout << e.msg << std::endl;
return 1;
}

std::cout << "W" << std::endl;
std::cout.precision(4);
std::cout.setf(std::ios::scientific);
std::cout << "         ";
for (vsize_t j = 0; j < k; ++j) {
std::cout << std::setw(15) << j + 1;
}
std::cout << std::endl;
for (vsize_t i = 0; i < m; ++i) {
std::cout << std::left << std::setw(4) << i + 1;
for (vsize_t j = 0; j < k; ++j) {
std::cout << std::setw(15) << w(i, j);
}
std::cout << std::endl;
}

std::cout << std::endl << "H^T" << std::endl;
std::cout.precision(4);
std::cout.setf(std::ios::scientific);
std::cout << "         ";
for (vsize_t i = 0; i < k; ++i) {
std::cout << std::setw(15) << i + 1;
}
std::cout << std::endl;
for (vsize_t j = 0; j < n; ++j) {
std::cout << std::left << std::setw(4) << j + 1;
for (vsize_t i = 0; i < k; ++i) {
std::cout << std::setw(15) << h(i, j);
}
std::cout << std::endl;
}

std::cout << std::endl;
std::cout << "The relative residual norm, ||A-WH|| / ||A||, is: "
<< rel_resid_norm << std::endl;

return 0;
}
```