[Rcpp-devel] Rcpp Parallel is not working

[Serguei Sokol]

Hi To,

Le 22/12/2016 à 12:55, To Duc Khanh a écrit :
> Dear all
>
> I am trying to apply RcppParallel for my function, but my code seems to be
> not working. Could you please help me to find any mistakes in my code or
> explain me why it is not working?
As your result is of type "aggregation" (res[0] += ...) wouldn't it be more
appropriate to use parallelReduce() instead of parallelFor()?

Another point to pay attention for is that index domain to cover is
square in your case (i and j) while parallelFor (and Reduce) are
intended to split plain linear index domains. So some additional effort
is needed to treat it properly.

Hoping it helps,
Serguei.

>
> I include my C/C++ code and R code at the end of this mail. Thank you for
> the help.
>
>
> // parallel version
> // [[Rcpp::depends(RcppParallel)]]
> #include<RcppParallel.h>
> #include <Rcpp.h>
>
> using namespace Rcpp;
> using namespace RcppParallel;
>
> inline double indvus(double a, double b, double c){
>   if((a < b) && (b < c)){
>     return 1.0;
>   } else if((a < b) && (b == c)){
>     return 0.5;
>   } else if((a == b) && (b < c)){
>     return 0.5;
>   } else if((a == b) && (b == c)){
>     return 1.0/6;
>   } else{
>     return 0.0;
>   }
> }
>
> struct VUS : public Worker {
>   // input matrix
>   const RVector<double> test;
>   const RMatrix<double> dise;
>   // output value
>   RVector<double> rvus;
>     // initialize from Rcpp input and output
>   VUS(const NumericVector test, const NumericMatrix dise, NumericVector rvus)
>     : test(test), dise(dise), rvus(rvus) {}
>   // function call operator that work for the specified range (begin/end)
>   void operator()(std::size_t begin, std::size_t end) {
>     for(std::size_t i = begin; i < end; i++) {
>       for(std::size_t j = begin; j < end; j++) {
>         if(j != i){
>           for(std::size_t k = begin; k < end; k++){
>             if((k != j) && (k != i)){
>               double temp = dise(i,0)*dise(j,1)*dise(k,2);
>               rvus[0] += temp*indvus(test[i], test[j], test[k]);
>               rvus[1] += temp;
>             }
>           }
>         }
>       }
>     }
>   }
> };
>
> // [[Rcpp::export]]
> NumericVector parallel_vus(NumericVector tt, NumericMatrix dd) {
>   // allocate the matrix we will return
>   NumericVector rvus(2);
>   // create the worker
>   VUS pvus(tt, dd, rvus);
>   // call it with parallelFor
>   parallelFor(0, tt.size(), pvus);
>   return rvus;
> }
>
>
> //serial version
> // [[Rcpp::export]]
> NumericVector serial_vus(NumericVector tt, NumericMatrix dd) {
>   NumericVector res(2);
>   int n = tt.size();
>   double temp = 0.0;
>   for(int i = 0; i < n; i++){
>     for(int j = 0; j < n; j++){
>       if(j != i){
>         for(int k = 0; k < n; k++){
>           if((k != j) && (k != i)){
>             temp = dd(i, 0)*dd(j, 1)*dd(k, 2);
>             res[0] += temp*indvus(tt[i], tt[j], tt[k]);
>             res[1] += temp;
>           }
>         }
>       }
>     }
>   }
>   return res;
> }
>
>
> ### R code
> dd = t(rmultinom(1000, 1, c(0.4, 0.35, 0.25)))
> tt <- numeric(1000)
> tt[dd[,1] == 1] <- rnorm(sum(dd[,1] == 1), 3, sqrt(1.2))
> tt[dd[,2] == 1] <- rnorm(sum(dd[,2] == 1), 6, sqrt(1.2))
> tt[dd[,3] == 1] <- rnorm(sum(dd[,3] == 1), 9, sqrt(1.2))
>
> serial_vus(tt,dd)
>
> parallel_vus(tt, dd)
>
>