[Rcpp-devel] Sparse matrix and RcppEigen

Dear all,

I need to solve a linear system using only the lower triangular matrix from
the Cholesky decomposition of a SPD matrix. My goal is to compute C^{-1/2}
V^{-1/2} where V^{-1/2} and C are known matrices.
For general matrix class MatrixXd I can do using the code below. But, now I
would like to do the same using the SparseMatrix class. I tried to adapt my
code, but it did not work. Any help will be welcome.


src <- '
using namespace Rcpp;
using namespace Eigen;
const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas));
const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs));
MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV);
return wrap(LT);'


srcsparse <- '
using Eigen::Map;
using Eigen::SparseMatrix;
using Eigen::LLT;

const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));
const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));
Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver;
SparseMatrix<double>  LT = solver.compute(Omega).matrixL()
*.solve(invSqrtV);*
return wrap(LT);'

The problem is that there is no .solve() method associated with the
matrixL().
It is just a intermediary step when computing

solver.compute(Omega).solve(invSqrtV)

This code computing C^{-1}V^{-1/2}, but I want only C^{-1/2}V^{-1/2}. It
will be used on the Generalized Kronecker product in the context of
multivariate spatial models.

Thank you.

Wagner Hugo Bonat
----------------------------------------------------------------------------------------------
Department of Mathematics and Computer Science (IMADA)
University of Southern Denmark (SDU) and
Laborat?rio de Estat?stica e Geoinforma??o (LEG)
Universidade Federal do Paran? (UFPR)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150608/62182ec6/attachment.html>