<font face="courier new,monospace">Dear Rcpp developers and users,<br clear="all"></font><div><div><font face="'courier new', monospace"><br></font></div><div><font face="'courier new', monospace">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.</font></div>
<div><font face="'courier new', monospace"><br></font></div><div><font face="'courier new', monospace">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?</font></div>
<div><font class="Apple-style-span" face="'courier new', monospace"><br></font></div><div><font class="Apple-style-span" face="'courier new', monospace">Thx!</font></div><div><font class="Apple-style-span" face="'courier new', monospace"><br>
</font></div><div><font class="Apple-style-span" face="'courier new', monospace">Gregor Gorjanc</font></div>
</div><div><font face="'courier new', monospace"><br></font></div><div><font face="'courier new', monospace"><div>SEXP fx(SEXP c1_, SEXP c2_, SEXP nI_, SEXP nP_, SEXP nT_, SEXP y_, SEXP P_, SEXP Px_) {</div>
<div> using namespace Rcpp ;</div><div> int i, j, k, t, p;</div><div> double c1 = Rcpp::as<double>(c1_); </div><div> double c2 = Rcpp::as<double>(c2_); </div><div> int nI = Rcpp::as<int>(nI_); </div>
<div> int nP = Rcpp::as<int>(nP_);</div><div> int nT = Rcpp::as<int>(nT_);</div><div> Rcpp::NumericMatrix ped(y_);</div><div> Rcpp::IntegerVector P(P_); </div><div> Rcpp::IntegerVector Px(Px_);</div><div>
Rcpp::NumericMatrix pa(nI+1, nT);</div><div> Rcpp::NumericMatrix w(nI+1, nT);</div><div> Rcpp::NumericMatrix xa(nI+1, nP*nT);</div><div> </div><div> // --- Compute ---</div><div> </div><div> for(i = 1; i < nI+1; i++) {</div>
<div> for(t = 0; t < nT; t++) {</div><div> pa(i, t) = c1 * ped(ped(i, 1), t+3) +</div><div> c2 * ped(ped(i, 2), t+3);</div><div> </div><div> w(i, t) = ped(i, t+3) - pa(i, t);</div><div>
</div>
<div> for(p = 0; p < nP; p++) {</div><div> j = Px[t] + p;</div><div> xa(i, j) = c1 * xa(ped(i, 1), j) +</div><div> c2 * xa(ped(i, 2), j); </div><div> }</div><div> k = Px[t] + P[i];</div>
<div> xa(i, k) += w(i, t);</div><div> }</div><div> }</div><div> </div><div> // --- Return ---</div><div> </div><div> Rcpp::List ret;</div><div> ret["pa"] = pa;</div><div> ret["w"] = w;</div>
<div> ret["xa"] = xa;</div><div> return(ret);</div><div>}</div><div><br></div></font></div>