# [Rcpp-devel] First foray into Rcpp and comparison of speed with R

Gregor GORJANC gregor.gorjanc at bf.uni-lj.si
Thu Sep 1 10:06:03 CEST 2011

```Dear Rcpp developers and users,

I wrote a simple application in R. Application is not very complex, but
since I have to deal with large data I have devoted a lot of attention to
vectorize R computations as much as possible. This lead me to a very decent
performance, but my computation involves a large outer loop, due to the
recursive nature of a problem. I looked into Rcpp and I was amazed that I
was able to come up with its use in my application even though I do not had
*any* practical experience with C++ before!!! I managed to half the
execution speed of my application, which is great. I kind of hoped for more,
but given that I vectorized much of R code and that computations are simple
additions and multiplications I wonder if I could get much more gains with
Rcpp foray.

Does anyone have time/will to quickly look at my code bellow to see if I am
missing something very obvious. I am basically doing one walk from top to
bottom of a rather large matrix (most of it is dense). In R I was able to
vectorize all within line computations, but in C++ I have to use explicit
looping, right?

Thx!

Gregor Gorjanc

SEXP fx(SEXP c1_, SEXP c2_, SEXP nI_, SEXP nP_, SEXP nT_, SEXP y_, SEXP P_,
SEXP Px_) {
using namespace Rcpp ;
int i, j, k, t, p;
double c1 = Rcpp::as<double>(c1_);
double c2 = Rcpp::as<double>(c2_);
int nI = Rcpp::as<int>(nI_);
int nP = Rcpp::as<int>(nP_);
int nT = Rcpp::as<int>(nT_);
Rcpp::NumericMatrix ped(y_);
Rcpp::IntegerVector P(P_);
Rcpp::IntegerVector Px(Px_);
Rcpp::NumericMatrix pa(nI+1, nT);
Rcpp::NumericMatrix  w(nI+1, nT);
Rcpp::NumericMatrix xa(nI+1, nP*nT);

// --- Compute ---

for(i = 1; i < nI+1; i++) {
for(t = 0; t < nT; t++) {
pa(i, t) = c1 * ped(ped(i, 1), t+3) +
c2 * ped(ped(i, 2), t+3);

w(i, t) = ped(i, t+3) - pa(i, t);

for(p = 0; p < nP; p++) {
j = Px[t] + p;
xa(i, j) = c1 * xa(ped(i, 1), j) +
c2 * xa(ped(i, 2), j);
}
k = Px[t] + P[i];
xa(i, k) += w(i, t);
}
}

// --- Return ---

Rcpp::List ret;
ret["pa"] = pa;
ret["w"]  = w;
ret["xa"] = xa;
return(ret);
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20110901/c0bf5c16/attachment.htm>
```