[Rcpp-devel] bounds checking disabled for RcppEigen?

[Andrew Slaughter]

First of all, great job with Rcpp - it's an awesome tool, and finally
gave me a reason to go back to my C++ books and get beyond "Hello World"!

I did have a question about RcppEigen, though: specifically, when I use
RcppEigen through inline, is bounds checking disabled? If not, how could I
go about disabling it? I did try adding the eigen_no_debug flag to my
includes, but it wouldn't compile. (I know the preferred way might be to
build a package, but inline is so much faster for quickly testing changes
to my code, I hate to give it up right now.)

As a separate but related question: is there agreement about good dynamic
(not fixed-size) C++ analogues to R's array? Especially anaologues that are
Rcpp-friendly. :) I really want to be able to quickly loop over and access
elements using something like "array(i,j,k,l,m)".

BTW, sorry if these are simple questions, but I'm learning C/C++ for this
project as I go along.

Thanks,
Andy


Background for the curious:

The vast majority of my code is looping over elements of multidimensional
numeric arrays (not fixed at compile time) and performing simple
multiplication, addition, and subtraction on those elements. Basically, the
time required to access elements of a numeric array is pretty important. A
much much smaller set of my code requires some matrix algebra and linear
algebra routines (matrix multiplication, cholesky decomp, etc.) on small
vectors and matrices. Right now I use BLAS/LAPACK, but it's harder for me
to read and experiment with, so I thought I'd give Eigen a try.

Since RcppArmadillo was much slower than using STL vectors last time I
tried it, I decided to benchmark Eigen before rewriting any code (see code
below comparing STL vector, NumericMatrix, and MatrixXd). It seems a good
bit slower, which I suspect might be due to bounds checking when I access
elements of the matrix.

Ultimately, I guess I could mix and match if I have to (use STL vectors for
the simple operations, and other classes for the matrix algebra stuff), but
for some reason that "feels" wrong..... but I'm no C++ person, so I'm
probably wrong.


# Here's a copy of the R syntax for the benchmark I mention above:

library(RcppEigen)
library(inline)
library(rbenchmark)

mat1<-matrix(rbinom(200*200, 1, 0.5), ncol=200)

header<-'
using namespace std;
using namespace Eigen;
'
RcppMatrixtest<-'
const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);
int p = X.nrow();
std::vector<double> V(p*p);
Eigen::MatrixXd EM(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){
        EM(i,j) = V[i + j*p] = X(i,j);
    }
}

double tran = 0.0;

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

return Rcpp::wrap(tran);
'
rcpp<-cxxfunction(
    signature(
    mat="numeric"
    ),
    body=RcppMatrixtest,
    includes=header,
    plugin="RcppEigen",
    verbose=T
    )

STLvectortest<-'
const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);
int p = X.nrow();
std::vector<double> V(p*p);
Eigen::MatrixXd EM(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){
        EM(i,j) = V[i + j*p] = X(i,j);
    }
}

double tran = 0.0;

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

return Rcpp::wrap(tran);
'
stl<-cxxfunction(
    signature(
    mat="numeric"
    ),
    body=STLvectortest,
    includes=header,
    plugin="RcppEigen",
    verbose=T
    )

Reigentest<-'
const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);
int p = X.nrow();
std::vector<double> V(p*p);
Eigen::MatrixXd EM(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){
        EM(i,j) = V[i + j*p] = X(i,j);
    }
}

double tran = 0.0;

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

return Rcpp::wrap(tran);
'
eig<-cxxfunction(
    signature(
    mat="numeric"
    ),
    body=Reigentest,
    includes=header,
    plugin="RcppEigen",
    verbose=T
    )

Rcppbase<-'
const NumericMatrix X = Rcpp::as<NumericMatrix>(mat);
int p = X.nrow();
std::vector<double> V(p*p);
Eigen::MatrixXd EM(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){
        EM(i,j) = V[i + j*p] = X(i,j);
    }
}

double tran = 0.0;

return Rcpp::wrap(tran);
'
rcppbase<-cxxfunction(
    signature(
    mat="numeric"
    ),
    body=Rcppbase,
    includes=header,
    plugin="RcppEigen",
    settings=settings,
    verbose=T
    )

# rcppbase is there to provide an estimate of how much time it takes to
copy, create
# necessary data structures, assuming it's not optimized away by the
compiler....

m1<-benchmark(
    rcpp(mat1),
    stl(mat1),
    eig(mat1),
    rcppbase(mat1),
    replications=500)

m1
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20120801/f0f15d9a/attachment.html>