[Rcppdevel] bounds checking disabled for RcppEigen?
Dirk Eddelbuettel
edd at debian.org
Thu Aug 2 02:44:04 CEST 2012
On 1 August 2012 at 17:41, Douglas Bates wrote:
 I don't really think that bounds checking is slowing down this
 calculation. You can make the whole thing run faster by the simple
 expedient of reordering the loops. For me, this version

 STLvectortest<'
 const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);
 int p = X.nrow();
 std::vector<double> V(p*p);

 // Copy data into STL vector and Eigen dynamic matrix
 for (int i=0; i<p; ++i){
 for (int j=0; j<p; ++j){
 V[i + j*p] = X(i,j);
 }
 }

 double tran = 0.0;

 for (int k=0; k<p; ++k){
 for (int j=0; j<p; ++j){
 for (int i=0; i<p; ++i){
 if ((j != i) && (k != i) && (j != k)) {
 tran += V[i + j*p]*V[j + k*p]*V[i + k*p];
 }
 }
 }
 }

 return Rcpp::wrap(tran);
 '

 is about 40% faster than your version because of localization of access.

 As is pointed out in one of the Rcpp vignettes, if you want to make
 accesses to elements of an array fast, you should use pointers. Also,
Well ... it just looks somewhat unC++like to me.
Romain and I beat this to death in one of the early example still in the
package as examples/ConvolveBenchmark/  started via 'buildAndRun.sh' (as we
may have written this in the preinline days). It compares the vector
convolution (from 'Writing R Extensions') in just about every way.
Fastest is ... a solution with loop unrolling. Pointers are on par with
simple class wrapper which inlines operator[].
Dirk
 remove repeated calculations of loop invariants. A good compiler can
 do such loop lifting for you but you might as well do it yourself when
 you this happening. The fastest version I can come up with is

 Ptrtest<'
 const NumericMatrix X(mat);
 int p = X.nrow();
 double *xpt = &X[0];
 double tran = 0.0;

 for (int k=0; k<p; ++k){
 double *xk = xpt + k*p;
 for (int j=0; j<p; ++j){
 if (j != k) {
 double *xj = xpt + j*p;
 double xjk = xpt[j + k*p];
 for (int i=0; i<p; ++i){
 if ((i != j) && (i != k)) {
 tran += xj[i] * xjk * xk[i];
 }
 }
 }
 }
 }

 return Rcpp::wrap(tran);
 '
 ptr <cxxfunction(signature(mat="numeric"), Ptrtest, "RcppEigen",
 header, verbose=TRUE)
 _______________________________________________
 Rcppdevel mailing list
 Rcppdevel at lists.rforge.rproject.org
 https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel

Dirk Eddelbuettel  edd at debian.org  http://dirk.eddelbuettel.com
More information about the Rcppdevel
mailing list