[Rcpp-devel] MappedSparseMatrix in RcppEigen package

[Shanika Wickramasuriya]

Hi,

I want to post the following post to the list. 

I want to use conjugate gradient algorithm implemented in RcppEigen package for solving large sparse matrices. 

Since I am new to Rcpp and C++, I just started with the dense matrices.

// [[Rcpp::depends(RcppEigen)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <Eigen/IterativeLinearSolvers>
using Eigen::SparseMatrix;
using Eigen::MappedSparseMatrix;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;
using Rcpp::as;
using Eigen::ConjugateGradient;
typedef Eigen::MappedSparseMatrix<double> MSpMat;

// [[Rcpp::export]]
VectorXd getEigenValues(SEXP As, SEXP bs) {
  const Map<MatrixXd> A(as<Map<MatrixXd> > (As));
  const Map<VectorXd> b(as<Map<VectorXd> > (bs));
  ConjugateGradient<MatrixXd> cg;
  cg.compute(A);
  VectorXd x=cg.solve(b);
  return x;
}

This works as expected. Therefore, I thought to extend this to suit to sparse matrices.

// [[Rcpp::depends(RcppEigen)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <Eigen/IterativeLinearSolvers>
using Eigen::SparseMatrix;
using Eigen::MappedSparseMatrix;
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::VectorXd;
using Rcpp::as;
using Eigen::ConjugateGradient;
typedef Eigen::MappedSparseMatrix<double> MSpMat;

// [[Rcpp::export]]
VectorXd getEigenValues(SEXP As, SEXP bs) {
  const MSpMat A = as<MSpMat>(As);
  const Map<VectorXd> b(as<Map<VectorXd> > (bs));
  ConjugateGradient<SparseMatrix<double> > cg;
  cg.compute(A);

  VectorXd x=cg.solve(b);
  return x;
}

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.

Thanks.

Regards,
 

Shanika Wickramasuriya
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141002/398f64e4/attachment.html>