RcppArmadillo on several cpu cores

Question

I have the following RccpArmadillo function that runs fine if I execute it on one cpu core. But if I use several cores, then R will crash. All the other Rcpp functions I created so far run fine on several cores (with foreach), only RccpArmadillo seems to be problematic. Any ideas how to fix that?

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

Answer 1

I generally recommend putting the code into a small package, and having each parallel worker load the package. That is known to work, both in serial and parallel, whereas relying on cppFunction() for an ad-hoc function may be too fragile for parallel execution.