[Rcpp-devel] bounds checking disabled for RcppEigen?

[Douglas Bates]

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,
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)