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