[Rcppdevel] Is there a way to pass references to R
Søren Højsgaard
sorenh at math.aau.dk
Sat Dec 22 18:47:57 CET 2012
Dear all,
To supplement Dirk's remarks about using RcppEigen in connection with sparse matrices: Here is an example where I use RcppEigen to moralized a directed acyclic graph represented by an adjacency matrix which is sparse (a dgCMatrix from R's Matrix package). The code can actually be made more efficient by replacing the forloops by a different type of iterator which ensures that only nonzero entries are visited. The "triplettrick" is the "recommended" way of building up a sparse matrix  as far as I can tell. There may be many other possible improvents which I am unaware of (as I don't know c++).
Best regards
Søren
RcppExport SEXP C_moralizeM ( SEXP XX_){
using Eigen::Map;
using namespace Rcpp;
typedef Eigen::MappedSparseMatrix<double> MSpMat;
typedef Eigen::SparseMatrix<double> SpMat;
SpMat X(as<MSpMat>(XX_));
typedef Eigen::Triplet<double> T;
std::vector<T> triplets;
triplets.reserve(X.nonZeros() * 2);
int nrX(X.rows());
int kk, ll, vv;
for (vv=0; vv<nrX; vv++){ /* consider vertex vv */
for (kk=0; kk<nrX; kk++){
if (X.coeff(kk, vv) != 0){ /* yes, kk>vv */
for (ll=kk+1; ll<nrX; ll++){
if (X.coeff(ll, vv) != 0){ /* yes, ll>vv */
if ((X.coeff(kk, ll)==0) && (X.coeff(ll, kk)==0)){ /* kk not~ ll */
triplets.push_back(T(kk, ll, 1));
triplets.push_back(T(ll, kk, 1));
}
}
}
}
}
}
SpMat ans(X.rows(), X.cols());
ans.setFromTriplets(triplets.begin(), triplets.end());
SpMat Xt(X.transpose());
ans = ans + Xt + X;
for (kk=0; kk<nrX; kk++){
for (ll=kk+1; ll<nrX; ll++){
if (ans.coeff(kk,ll)!=0){
ans.coeffRef(kk,ll)=1;
ans.coeffRef(ll,kk)=1;
}
}
}
ans.makeCompressed();
return(wrap(ans));
}
Original Message
From: rcppdevelbounces at lists.rforge.rproject.org [mailto:rcppdevelbounces at lists.rforge.rproject.org] On Behalf Of Dirk Eddelbuettel
Sent: 22. december 2012 01:41
To: matzke at berkeley.edu
Cc: rcppdevel at lists.rforge.rproject.org
Subject: Re: [Rcppdevel] Is there a way to pass references to R
Nick,
On 21 December 2012 at 16:16, Nick Matzke wrote:
 I'm working on a problem that involves calculating a very large,
 sparse transition probability matrix, P. It takes significant time to
 run through the whole matrix, check all the transitions, and put in
 the appropriate probabilities, based on multiplying parameter values
 j, d, e, etc.

 Each time I update j, or d, or e, I would rather not go reprocess
 which cells of P get j, j*d, d*e, etc. It would be better to just do
 this once, with j, d, and e each just referring to a memory address.
 I could then update the values at that memory address and reuse the
 matrix P without rebuilding it again and again.

 Currently, my script builds P via an Rcpp function, then passes it to
 R, where it gets used (with the values
 calculated) a few hundred times (via another Rcpp function); then j,
 d, and e are updated and I have to recalculate P.

 So, my question is, is there a way to make P, currently a vector of
 float vectors, into a vector of float pointers, which can be passed to
 R and then back to Rcpp? Or is this just ridiculous?
There will be lots of different ways to skin this cat. Sparse matrices may be one, I have no experience there. They now exist in Armadillo, though without [at current] glue code for as<> and wrap in RcppArmadillo; this does exist in Eigen and RcppEigen. Others may help you on sparse stuff.
You can also define you own structure to hold those vectors, and then simply pass one single pointer around using the Rcpp::XPtr wrapper class for R's external pointer. If data volume is an issue, external pointers are your friends. (Some of) the database class package uses this, and as does the bigmemory familiy of packages.
But you should do some measuring and profiling before you go off and rearchitect your application based on a hunch. Great slide deck seen today and retweeted:
The only good intuition: 'I should time time this'
 Andrei Alexandrescu
Via http://isocpp.org/blog/2012/12/threeoptimizationtipsalexandrescu
and I think he has that exactly right. And he has about a gazillion times as much street cred on C++ as I do ....
Dirk
PS You need to work on your signature. Not sure it's long and detailed enough.

Dirk Eddelbuettel  edd at debian.org  http://dirk.eddelbuettel.com _______________________________________________
Rcppdevel mailing list
Rcppdevel at lists.rforge.rproject.org
https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel
More information about the Rcppdevel
mailing list