Stochastic Simulation by OpenMP

[option72]

Hello, I am writing because I am experiencing troubles with parallel programming.
The following code is a parallel version of a stochastic model for an actuarial simulation.
Even with OpenMP there are unsolvable errors I am unable to fix (may be useful: win11, R 4.5.1, Rstudio 2025.05.1 Build 513)
Very briefly, I would like to repeat in parallel the same stochastic simulation model for evolving a popolution.
The sequential version works (EvoluzionePopolazione_v01.cpp), while the parallel version crashes
(EvoluzionePopolazione_par_v03_Gemini.cpp).

Here is my parallel code...

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


#include <omp.h>


// [[Rcpp::plugins(openmp)]]

using namespace Rcpp;

inline Rcpp::NumericMatrix random_mtx(int n, int m) {
  Rcpp::NumericVector res = dqrng::dqrunif(n * m, 0.0, 1.0);
  res.attr("dim") = Rcpp::Dimension(n, m);
  return Rcpp::as<Rcpp::NumericMatrix>(res);
}

inline Rcpp::IntegerMatrix permutazioni(int n,Rcpp::IntegerVector x)  {
  Rcpp::IntegerMatrix perm(n,x.size());
  for(int i = 0; i < n; ++i) {
    Rcpp::IntegerVector indices = dqrng::dqsample_int(x.size(), x.size());
    for(int j = 0; j <x.size();j++){
      perm(i,j)=x[indices[j]];
    }
  }
  return perm;
}

inline Rcpp::NumericVector DeterminazioniEvento(  
    int evento, int stato, int tp_sesso,
    double prob_q, double prob_tavM

){
  
  NumericVector mDetEvento(3);
  // 0. evento,       1.  accadimento,    2. stato
  
  mDetEvento[0] = evento;
  
  switch(evento) {
  
  case 1: 
  { 
    
    if (prob_q <= prob_tavM) {
    mDetEvento[1] = 1; //morte
    mDetEvento[2] = 5; 
    
    
  } else {
    mDetEvento[1] = 0; //sopravvivenza
    mDetEvento[2] = stato; //non cambia lo stato
  }
  }
    break;
    
  default: 
  {
    
  } 
  break;
  } 
  
  return mDetEvento;
}

inline Rcpp::NumericMatrix SintesiSimulazione(NumericMatrix mSim, int nR, int nAnni){
  
  Rcpp::NumericMatrix Sintesi(1, 4 * nAnni);
  
  for(int i = 0; i < nR; i++) {
    
    if(mSim(i,12) != 0) {Sintesi(0,            mSim(i,12)-1)++; } //conteggio decessi per epoca
    if(mSim(i,13) != 0) {Sintesi(0,    nAnni + mSim(i,13)-1)++; } //conteggio reversibilità per epoca
    if(mSim(i,14) != 0) {Sintesi(0,2 * nAnni + mSim(i,14)-1)++; } //conteggio indirette per epoca
    if(mSim(i,15) != 0) {Sintesi(0,3 * nAnni + mSim(i,15)-1)++; } //conteggio estinzioni per epoca
  } //next i
  
  return Sintesi;
  
}


// [[Rcpp::export]]
Rcpp::NumericMatrix popolazione_par_v02(NumericMatrix Iscritti, NumericMatrix Tmorte, NumericMatrix Teta, 
                                        int nRip, int nAnni, int nThreads) {
  
  int nR = Iscritti.nrow();
  int nC = Iscritti.ncol();
  
  Rcpp::NumericMatrix mSintesi(nRip, 1 + 4 * nAnni); //4 eventi
  if (nThreads > 0) {
    omp_set_num_threads(nThreads);
  }
  
  
  
#pragma omp parallel for
  for(int r = 0; r < nRip; r++) {
      
    Rcpp::NumericMatrix mSim(nR, 16);
    Rcpp::NumericMatrix mpcMorte = random_mtx(nR, nAnni);
    Rcpp::NumericMatrix SintesiAnno(1, 4 * nAnni);
    
    double prob_q=0.0; 
    double prob_tavM=0.0; 
    int evento = 0;
    
    for(int s = 0; s < nC; s++) {
      for(int z = 0; z < nR; z++) {
        mSim(z,s) = Iscritti(z,s);
      }
    }
    
    int eta_anag =0;
    int eta_q=0;
    int tp_sesso=0;
    
    int stato=0;
    int inC= -1; //incremento epoca
    
    for(int t = 1; t <= nAnni; t++) {
      inC = inC + 1;
      for(int i = 0; i < nR; i++) {
        NumericVector eventoRis(3);
        tp_sesso = mSim(i,5);
        stato = mSim(i,6);
        if (t == 1) {
          eta_anag = mSim(i,1);
          eta_q = mSim(i,4);     
        } else {
          eta_anag = mSim(i,1) + inC; //incremento età anagrafica       
          eta_q = mSim(i,4) + inC; //incremento età aggiustata
        } 
        
        evento = 1; //unico evento previsto
        
        if(mSim(i,6) != 5) {
          prob_q = mpcMorte(i, t - 1); 
          prob_tavM = Tmorte(eta_q, tp_sesso);   
          eventoRis = DeterminazioniEvento (evento, stato, tp_sesso, prob_q, prob_tavM);
          
          if (eventoRis[0] == 1 && eventoRis[1] == 1) {  
            mSim(i,6) = eventoRis[2]; //nuovo stato
            mSim(i,12) = t; //epoca decesso 
            mSim(i,15) = t; //epoca estinzione
          } 
        }
      } //next i
    }  //next t
    
    SintesiAnno = SintesiSimulazione(mSim, nR, nAnni);
    
    mSintesi(r,0) = r;
    for(int g = 0; g < SintesiAnno.ncol(); g++) {
      mSintesi(r,g + 1) = SintesiAnno(0,g);
    }
    
  } //next r
  
  return mSintesi;
}

it is possible to test the simulation and to reproduce the errors using the following R code:

library(dqrng)
library(Rcpp)

sourceCpp("~/R/MAGIS/CPA/codice/EvoluzionePopolazione_par_v03_Gemini.cpp")


Teta_i <- data.frame(eta= seq(1900, 1900+121,1), delta = sample(-2:2, 122, replace=T))
Tmorte_i <- data.frame(eta= seq(0, 121,1), m = runif(122, 0, 1), f = runif(122, 0, 1))

campiVuoti01 <- data.frame(CAMPO_01 = rep(NA,1000),
                           ETA_ULTIMO_FIGLIO = rep(0,1000), SESSO_ULTIMO_FIGLIO = rep(0,1000), EPOCA_INGRESSO = rep(0,1000), 
                           CAMPO_02 = rep(NA,1000), 
                           EPOCA_DECESSO = rep(0,1000), EPOCA_REVERSIBILITA = rep(0,1000), EPOCA_INDIRETTA = rep(0,1000), EPOCA_ESTINZIONE = rep(0,1000)
)

PopolazioneSub1 <- data.frame(NUM_ISCRIZIONE = seq(1000, 1999,1), 
                              ETA_DT_INTERESSE = sample(20:100, 1000, replace=T), 
                              TP_SESSO = sample(1:2, 1000, replace=T), 
                              STATO_CONTRIBUZIONE = 1)

PopolazioneSub1$ETA_SIM <- PopolazioneSub1$ETA_DT_INTERESSE
PopolazioneSub1$TP_SESSO_SIM <- PopolazioneSub1$TP_SESSO # Indice 1 o 2
PopolazioneSub1$STATO_SIM <- PopolazioneSub1$STATO_CONTRIBUZIONE

PopolazioneEvo <- data.frame(cbind(PopolazioneSub1, campiVuoti01 ) )



mTmorte <- as.matrix(Tmorte_i)
mTeta <- as.matrix(Teta_i)
mPopolazione <- as.matrix(PopolazioneEvo)

nRip =  10 #numero ripetizioni
nAnni = 30 #numero epoche massimo

dqrng::dqRNGkind("Xoshiro256++")
dqrng::register_methods()
set.seed(20250428)

ncores = 4

popolazione_par_res <- as.matrix(popolazione_par_v02(mPopolazione, mTmorte, mTeta, 
                                                     nRip, nAnni, ncores))


print(head(popolazione_par_res))
print(dim(popolazione_par_res))