<html><body><div style="color:#000; background-color:#fff; font-family:verdana, helvetica, sans-serif;font-size:13px"><div class="" style=""><span class="" style="">Hi,</span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><span class="" style=""><br class="" style=""></span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><span class="" style="">I want to post the following post to the list. </span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><span class="" style=""><br class="" style=""></span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style:
normal; background-color: transparent;" class=""><span class="" style="">I want to use conjugate gradient algorithm implemented in RcppEigen package for solving large sparse matrices. </span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><span class="" style=""><br class="" style=""></span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><span class="" style="">Since I am new to Rcpp and C++, I just started with the dense matrices.</span></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><span class="" style=""><br class="" style=""></span></div><div style="background-color: transparent;" class="">//
[[Rcpp::depends(RcppEigen)]]</div><div style="background-color: transparent;" class="">#include <Rcpp.h></div><div style="background-color: transparent;" class="">#include <RcppEigen.h></div><div style="background-color: transparent;" class="">#include <Eigen/IterativeLinearSolvers></div><div style="background-color: transparent;" class="">using Eigen::SparseMatrix;</div><div style="background-color: transparent;" class="">using Eigen::MappedSparseMatrix;</div><div style="background-color: transparent;" class="">using Eigen::Map;</div><div style="background-color: transparent;" class="">using Eigen::MatrixXd;</div><div style="background-color: transparent;" class="">using Eigen::VectorXd;</div><div style="background-color: transparent;" class="">using Rcpp::as;</div><div style="background-color: transparent;" class="">using Eigen::ConjugateGradient;</div><div style="background-color: transparent;" class="">typedef
Eigen::MappedSparseMatrix<double> MSpMat;</div><div style="background-color: transparent;" class=""><br class="" style=""></div><div style="background-color: transparent;" class="">// [[Rcpp::export]]</div><div style="background-color: transparent;" class="">VectorXd getEigenValues(SEXP As, SEXP bs) {</div><div style="background-color: transparent;" class=""> <span style="background-color: transparent;" class="">const Map<MatrixXd> A(as<Map<MatrixXd> > (As));</span></div><div style="background-color: transparent;" class=""> const Map<VectorXd> b(as<Map<VectorXd> > (bs));</div><div style="background-color: transparent;" class=""> ConjugateGradient<MatrixXd> cg;</div><div style="background-color: transparent;" class=""> cg.compute(A);</div><div style="background-color: transparent;" class=""> VectorXd x=cg.solve(b);</div><div style="background-color: transparent;"
class=""> return x;</div><div style="background-color: transparent;" class=""><span class="" style=""></span></div><div style="background-color: transparent;" class="">}</div><div style="background-color: transparent;" class=""><br class="" style=""></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class="">This works as expected. Therefore, I thought to extend this to suit to sparse matrices.</div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><br class="" style=""></div><div style="background-color: transparent;" class="">// [[Rcpp::depends(RcppEigen)]]</div><div style="background-color: transparent;" class="">#include <Rcpp.h></div><div style="background-color: transparent;" class="">#include <RcppEigen.h></div><div
style="background-color: transparent;" class="">#include <Eigen/IterativeLinearSolvers></div><div style="background-color: transparent;" class="">using Eigen::SparseMatrix;</div><div style="background-color: transparent;" class="">using Eigen::MappedSparseMatrix;</div><div style="background-color: transparent;" class="">using Eigen::Map;</div><div style="background-color: transparent;" class="">using Eigen::MatrixXd;</div><div style="background-color: transparent;" class="">using Eigen::VectorXd;</div><div style="background-color: transparent;" class="">using Rcpp::as;</div><div style="background-color: transparent;" class="">using Eigen::ConjugateGradient;</div><div style="background-color: transparent;" class="">typedef Eigen::MappedSparseMatrix<double> MSpMat;</div><div style="background-color: transparent;" class=""><br class="" style=""></div><div style="background-color: transparent;" class="">// [[Rcpp::export]]</div><div
style="background-color: transparent;" class="">VectorXd getEigenValues(SEXP As, SEXP bs) {</div><div style="background-color: transparent;" class=""> const MSpMat A = as<MSpMat>(As);</div><div style="background-color: transparent;" class=""> const Map<VectorXd> b(as<Map<VectorXd> > (bs));</div><div style="background-color: transparent;" class=""> ConjugateGradient<SparseMatrix<double> > cg;</div><div style="background-color: transparent;" class=""><span style="background-color: transparent;"> cg.compute(A);</span><br></div><div style="background-color: transparent;" class=""> VectorXd x=cg.solve(b);</div><div style="background-color: transparent;" class=""> return x;</div><div style="background-color: transparent;" class="">}</div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;"
class=""><br class="" style=""></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class="">However, this tends to give really strange results. The code in itself is also not giving any errors. Hope someone could help me to correct this error.</div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><br></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class="">Thanks.</div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""><br></div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal;
background-color: transparent;" class="">Regards,</div><div style="color: rgb(0, 0, 0); font-size: 13px; font-family: verdana, helvetica, sans-serif; font-style: normal; background-color: transparent;" class=""> <br></div><div class="" style="">Shanika Wickramasuriya</div></div></body></html>