2016-03-02 78 views
-1

我有以下RccpArmadillo函数运行良好,如果我在一个CPU核心上执行它。但是如果我使用多个内核,那么R将会崩溃。到目前为止我创建的所有其他Rcpp函数在几个核心上运行良好(使用foreach),但只有RccpArmadillo似乎有问题。任何想法如何解决这个问题?RcppArmadillo在几个CPU核心

cppFunction('double augmentedDickeyFullerCpp(NumericVector a, NumericVector b, double gamma, double mu, int lags) { 

    if (gamma < 0) { 
     return 0; 
    } 

    int n = a.size()-1; 
    int lags2 = lags + 1; 
    // first rows, then columns 
    NumericMatrix x(n-lags2,lags2); 
    NumericMatrix zdifflag(n-lags2+1,lags2); 
    NumericVector diff(n); 
    NumericVector zdiff(n-lags2+1); 
    NumericVector residuals(n+1); 
    residuals[0] = a[0] - gamma * b[0] - mu; 


    // residuals a is y and b is x 
    for(int i = 1; i < n+1; i++) { 
     residuals[i] = a[i] - gamma * b[i] - mu; 
     diff[i-1] = residuals[i] - residuals[i-1]; 
    } 
    for(int i = 0; i < n-lags2+1; i++) { 
     zdifflag[0,i] = residuals[i+lags2-1]; 
    } 

    for(int j = 0; j < n-lags2+1; j++) { 
     for(int i = 0; i < lags2; i++) { 
      x(j,i) = diff[j+lags2-1-i]; 
      if (i > 0) { 
       zdifflag(j,i) = x(j,i); 
      } 
     } 
     zdiff[j] = x(j,0); 
    } 

    int length = zdifflag.nrow(), k = zdifflag.ncol(); 

    arma::mat X(zdifflag.begin(), length, k, false); // reuses memory and avoids extra copy 
    arma::colvec y(zdiff.begin(), zdiff.size(), false); 
    arma::colvec coef = arma::solve(X, y); // fit model y ~ X 
    arma::colvec res = y - X*coef;   // residuals 

    // std.errors of coefficients 
    //arma::colvec res = y - X*coef[0]; 
    // sqrt(sum(residuals^2)/(length - k)) 
    double s2 = std::inner_product(res.begin(), res.end(), res.begin(), 0.0)/(length - k);               
    arma::colvec std_err = arma::sqrt(s2 * arma::diagvec(arma::pinv(arma::trans(X)*X))); 

    return coef[0]/std_err[0];  

}',depends = "RcppArmadillo", includes="#include <RcppArmadillo.h>") 

回答

2

我通常建议将代码放入一个小包中,让每个并行工作者加载包。这已知可以串行和并行方式工作,而依靠cppFunction()作为专门功能可能对于并行执行来说太脆弱了。