Why Rcpp openmp parallized codes sometimes gave: Segmentation fault (core dumped)

Question

I wrote an Rcpp code which worked before I make that parallel. I parallelized my code with openmp. this is my 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]);
}

and in my R code:

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

After that some times my code give results and sometime give this error:

Segmentation fault (core dumped)

I am absolutely newbi to openmp and learned here. What Can be the reason? Any help will be appreciated.

Edit When I get this error, output of sessionInfo() is error too:

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

but when I can get ouput, Output of sessionInfo is:

> 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

And I generate my data with this function:

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

Here a possible use of RMatrix that should help with the segmentation faults:

// [[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] ...