ex_f01sb.cpp

// 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"
intmain(void){
std::cout<<"nagcpp::matop::real_nmf_rcomm Example"<<std::endl;
typedefstd::vector<double>::size_typevsize_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_tm=7;
size_tn=6;
size_tk=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::OptionalF01SBopt;
opt.seed(234);
opt.errtol(1.0e-3);
// we will terminate after 200 iterations
vsize_tmaxit=200;
// we are going to pretend that A is q x q, rather than m x n
size_tq=std::max(m,n);
std::vector<double>tmp_x(q,0);
std::vector<double>tmp_y(q);
vsize_ticount;
doublerel_resid_norm;
// find a non-negative matrix factorixation A ~= WH
try{
nagcpp::utility::CopyableCommcomm;
nagcpp::types::f77_integerirevcm=0;
for(icount=0;icount<maxit;){
nagcpp::matop::real_nmf_rcomm(irevcm,w,h,ht,comm,opt);
if(irevcm==1){
++icount;
}elseif(irevcm==2){
// compute a^t * w and store in ht
for(vsize_tj=0;j<k;++j){
for(vsize_ti=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_ti=0;i<n;++i){
ht(i,j)=tmp_y[i];
}
}
}elseif(irevcm==3){
// compute a * ht and store in w (f11mk)
for(vsize_tj=0;j<k;++j){
for(vsize_ti=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_ti=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_ti=0;i<a.size();++i){
da(irowix[i]-1,icolix[i]-1)=a[i];
}
// get ||A|| (f06ra)
doublenorma=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)
doublenorm=nagcpp::blas::dlange("F",da);
rel_resid_norm=norm/norma;
}catch(nagcpp::error_handler::Exception&e){
std::cout<<e.msg<<std::endl;
return1;
}
std::cout<<"W"<<std::endl;
std::cout.precision(4);
std::cout.setf(std::ios::scientific);
std::cout<<" ";
for(vsize_tj=0;j<k;++j){
std::cout<<std::setw(15)<<j+1;
}
std::cout<<std::endl;
for(vsize_ti=0;i<m;++i){
std::cout<<std::left<<std::setw(4)<<i+1;
for(vsize_tj=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_ti=0;i<k;++i){
std::cout<<std::setw(15)<<i+1;
}
std::cout<<std::endl;
for(vsize_tj=0;j<n;++j){
std::cout<<std::left<<std::setw(4)<<j+1;
for(vsize_ti=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;
return0;
}