[Rcpp-devel] bounds checking disabled for RcppEigen?

[Dirk Eddelbuettel]

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 un-C++-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 pre-inline 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)
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com