Pourquoi les codes parallélisés Rcpp openmp donnaient parfois: Erreur de segmentation (core dumped)

Question

J'ai écrit un code Rcpp qui a fonctionné avant de faire ce parallèle. J'ai parallélisé mon code avec openmp. voici mon cppcode:

#include <Rcpp.h> 
#include <omp.h> 

using namespace Rcpp; 
// [[Rcpp::export]] 
float fsvalue(NumericMatrix X1, NumericMatrix X2,int n_cpu=2) { 
    int p = X1.ncol(); 
    if (p !=X2.ncol()) 
    stop ("Number of Dimentions should be equal"); 
    int n1 = X1.nrow(); 
    int n2 = X2.nrow(); 
    float gamma = (float) n1/ (float) n2; 
    int step = 0; int max=p*n1*(n1-1)*n2*(n2-1); 
    int sums = 0; 
    float FS_Sum = 0; 
    omp_set_dynamic(0);   // Explicitly disable dynamic teams 
    omp_set_num_threads(n_cpu); // Use n_cpu threads for all 
    /* 
    * 
    */ 
    #pragma omp parallel 
    { 
    float sum_thread = 0; 
    #pragma omp for 
    for (int k=0; k<p;k++){ 
     for (int i=0; i<n1;i++){ 
     for (int j=0;j<n1;j++){ 
      for (int s=0;s<n2;s++){ 
      for (int t=0;t<n2;t++){ 
       if (i==j || s==t) 
       continue; 
       step++; 
       sum_thread += ((X1[i*p + k] - X2[s*p + k])*(X1[j*p + k]) - X2[t*p + k])/(gamma); 
      } 
      } 
     } 
     } 
    } 
    #pragma omp critical 
    { 
     FS_Sum += sum_thread; 
    } 
    } 
    return (FS_Sum/(n1*(n1-1)*n2*(n2-1))); 
    //return (X1[3*10 + 1]); 
} 

et dans mon code R:

library("Rcpp") 
#generate two matrixes for test 
source("MovingAverage_data.r") 
Sys.setenv("PKG_CXXFLAGS"="-fopenmp") 
sourceCpp("teststat.cpp") 
fsvalue(fsdata1,fsdata2) 

Après que quelques fois mon code donnent des résultats et parfois donnent cette erreur:

Segmentation fault (core dumped) 

Je suis tout à fait newbi à openmp et appris here. Quelle peut être la raison? Toute aide sera appréciée.

Modifier Quand je reçois cette erreur, la sortie de sessionInfo() est une erreur aussi:

> sessionInfo() 
Error in eval(formal.args[[deparse(substitute(arg))]]) : there is no .Internal function '�l�' 

mais quand je peux obtenir ouput, sortie de sessionInfo est:

> sessionInfo() 
R version 3.2.3 (2015-12-10) 
Platform: x86_64-pc-linux-gnu (64-bit) 
Running under: Ubuntu 16.04.3 LTS 

locale: 
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C   LC_TIME=fa_IR   
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=fa_IR  LC_MESSAGES=en_US.UTF-8 
[7] LC_PAPER=fa_IR   LC_NAME=C    LC_ADDRESS=C   
[10] LC_TELEPHONE=C   LC_MEASUREMENT=fa_IR LC_IDENTIFICATION=C  

attached base packages: 
[1] stats  graphics grDevices utils  datasets methods base  

other attached packages: 
[1] Rcpp_0.12.12 

loaded via a namespace (and not attached): 
[1] tools_3.2.3 

Et je génère mes données avec cette fonction:

# @getdata function 
# @params i, p, i \in \{1,2\} and k in \{i,\ldots,p\} 
# @return X_{ijk} vector - X{ijk} has n_i length 
get_data_of_dim_k <- function(i,n_i,k,p,lambda=10,scenario=1,randomize=TRUE,seed=342373){ 
    if (k>p){ 
     stop("your dimention (k) should not be greater than maximum number of dimentions (p)"); 
    } 
    l_i = c(3,4) 
    rho_generator <- function(){ 
     #rho11 rho12 rho13 - rho21 rho22 rho23 rho24 
     rho = c(); 
     set.seed(seed) 
     for (i in c(1:8)){ 
      rho[i]= runif(1,2,3) 
     } 
     rho = matrix(rho,nrow=2); 
     rho[1,4] = 0 
     if (randomize){ 
      set.seed(as.numeric(format(Sys.time(), "%OS3"))) 
     } 
     return (rho) 
    } 
    mu_generator = function(p,lambda=10,Seed){ 
     set.seed(Seed) 
     mu = runif(p,0,lambda) 
     set.seed(as.numeric(format(Sys.time(), "%OS3"))) 
     return(mu) 
    } 
    mu_ij = mu_generator(p,Seed = seed) 
    rho = rho_generator() 

    #Check this on windows and macintash 
    set.seed(as.numeric(format(Sys.time(), "%OS3"))) 
    #innov = runif(100, -0.5, 0.5) -------> http://stackoverflow.com/questions/39925470/simulate-an-ar1-process-with-uniform-innovations 
    rho = rho_generator() 
    #TODO add \mu_{ij} 
    switch(scenario, 
     "1" = { 
      X_ij = arima.sim(list(ma=rho[i,c(1,l_i[i])]),n=n_i)+mu_ij[k] 
     }, 
     "2" = { 
      if (k<=(p/2)){ 
       X_ij = arima.sim(list(ma=rho[i,c(1,l_i[i])]),n=n_i)+mu_ij[k] 
      }else{ 
       X_ij = arima.sim(list(ma=rho[i,c(1,l_i[i])]),n=n_i,innov = rgamma(n_i,1,rate=4))+mu_ij[k]   
      } 
     } 
    ) 
    return (as.vector(X_ij)) 
} 
#get_data_of_dim_k(1,35,1,25) 
get_data = function(i,n_i,p){ 
    data = matrix(get_data_of_dim_k(i,n_i,1,p),ncol=1) 
    for (j in c(2:p)){ 
     data = cbind(data,get_data_of_dim_k(i,n_i,j,p)) 
    } 
    return(data) 
} 
    # @var fsdata1 just for test 
fsdata1 = get_data(i = 1,n_i = 8,p = 10) 

# @var fsdata2 just for test 
fsdata2 = get_data(i= 2, n_i=8,p = 10) 

Answer 1

Ici une utilisation possible de RMatrix qui devrait aider les erreurs de segmentation:

// [[Rcpp::depends(RcppParallel)]] 
#include <RcppParallel.h> 
#include <Rcpp.h> 
#include <omp.h> 

using namespace Rcpp; 
// [[Rcpp::export]] 
float fsvalue(NumericMatrix x1, NumericMatrix x2,int n_cpu=2) { 
    RcppParallel::RMatrix<double> X1(x1); 
    RcppParallel::RMatrix<double> X2(x2); 
    ... [your code as above] ...