Dear All, <br>I have modified the code as following. Since I found the old code had the memory allocation issue, I calculated each term of the matrix with dim 4. The speed is pretty much good now. But since I am not familiar with Rcpp, the code sometimes will get error, like<br>
<br> *** caught segfault ***<br>address 0x16, cause 'memory not mapped'<br>Error in betahat(ker, inv[l], cbind(inv[1:ll], inv[(ll + 1):(2 * ll)]),  :<br>  badly formed function expression<br>Calls: system.time ... apply -> FUN -> betahat -> .External -> cpp_exception<br>
Execution halted<br><br>I am not sure whether my following code has some big issue mentioned in the above error hint. I repeatedly call the following function in my entire code.<br><br>// [[Rcpp::depends(RcppArmadillo)]]<br>
<br>#include <RcppArmadillo.h><br><br>using namespace Rcpp;<br><br>// [[Rcpp::export]]<br>List betahat(Function ker, double t0, NumericMatrix Xr, NumericMatrix yr, NumericVector tr, double h, int m) {<br>  int n = Xr.nrow(), p = Xr.ncol();<br>
  arma::mat X(Xr.begin(), n, p, false);<br>  arma::mat y(yr.begin(), n, 1, false);<br>  arma::colvec t(tr.begin(), tr.size(), false);<br>  arma::mat T = X;<br>  T.each_col() %= (t-t0)/h;<br>  arma::vec K =as<arma::vec>(ker(tr-t0,h))/m;<br>
  double L1 = arma::accu(K%X.col(0)%X.col(0));<br>  double L2 = arma::accu(K%X.col(0)%X.col(1));<br>  double L3 = arma::accu(K%X.col(1)%X.col(1));<br>  double L4 = arma::accu(K%X.col(0)%T.col(0));<br>  double L5 = arma::accu(K%X.col(1)%T.col(0));<br>
  double L6 = arma::accu(K%X.col(1)%T.col(1));<br>  double L7 = arma::accu(K%T.col(0)%T.col(0));<br>  double L8 = arma::accu(K%T.col(0)%T.col(1));<br>  double L9 = arma::accu(K%T.col(1)%T.col(1));<br>  double R1 = arma::accu(K%X.col(0)%y);<br>
  double R2 = arma::accu(K%X.col(1)%y);<br>  double R3 = arma::accu(K%T.col(0)%y);<br>  double R4 = arma::accu(K%T.col(1)%y);<br>  arma::mat L(2*p,2*p);<br>  L(0,0)=L1;L(0,1)=L2;L(0,2)=L4;L(0,3)=L5;<br>  L(1,0)=L2;L(1,1)=L3;L(1,2)=L5;L(1,3)=L6;<br>
  L(2,0)=L4;L(2,1)=L5;L(2,2)=L7;L(2,3)=L8;<br>  L(3,0)=L5;L(3,1)=L6;L(3,2)=L8;L(3,3)=L9;<br>  arma::mat R(2*p,1);<br>  R(0,0)=R1;R(1,0)=R2;R(2,0)=R3;R(3,0)=R4;<br>  arma::vec betahat = arma::solve(L,R);<br>  arma::colvec betahat0(betahat.begin(),betahat.size()/2,false);<br>
  return List::create(Named("betahat") = betahat0);<br>}<br><br clear="all"><div><div>Best wishes!</div><div> </div><div>Honglang Wang</div><div> </div><div>Office C402 Wells Hall</div><div>Department of Statistics and Probability</div>
<div>Michigan State University</div><div>1579 I Spartan Village, East Lansing, MI 48823</div><div><a href="mailto:wangho16@msu.edu" target="_blank">wangho16@msu.edu</a></div></div><br>
<br><br><div class="gmail_quote">On Thu, Dec 6, 2012 at 2:10 AM, Darren Cook <span dir="ltr"><<a href="mailto:darren@dcook.org" target="_blank">darren@dcook.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="im">> this part will always make your code crawl along:<br>
<br>
</div><div class="im">>>   arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;<br>
<br>
</div>I'd be very interested to see the before/after version of your code<br>
code, Honglang. With timings. It'll make a good blog post or academic<br>
paper (depending on just how many learnings you get from the process) :-)<br>
<span class="HOEnZb"><font color="#888888"><br>
Darren<br>
<br>
<br>
--<br>
Darren Cook, Software Researcher/Developer<br>
<br>
<a href="http://dcook.org/work/" target="_blank">http://dcook.org/work/</a> (About me and my work)<br>
<a href="http://dcook.org/blogs.html" target="_blank">http://dcook.org/blogs.html</a> (My blogs and articles)<br>
</font></span><div class="HOEnZb"><div class="h5">_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</div></div></blockquote></div><br>