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

[Gregor GORJANC]

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>