[Rcpp-devel] bounds checking disabled for RcppEigen?

[Douglas Bates]

On Wed, Aug 1, 2012 at 12:46 PM, Andrew Slaughter
<andrew.slaughter at gmail.com> wrote:
> 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.)

Actually I think the debug checking is, by default, turned off in R
package compilation with the -DNDEBUG flag unless you work out a way
to override it.  I'm not sure if that applies in code compiled through
the inline package.

One of the things I would definitely do in your RcppEigen is to use a
Map<MatrixXd> not a MatrixXd for the input value, as this avoids a
copy operation.  Another thing to check is whether the colwise or
rowwise methods can be used to advantage.  I'll take a look and see if
I can come up with a faster formulation in RcppEigen.

> 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
>
>
>
>
>
> _______________________________________________
> 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