Skip to content

[Rcpp-devel] Rcpp Parallel is not working

2 messages · To Duc Khanh, Serguei Sokol

To Duc Khanh
# 
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?

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)
Serguei Sokol
# 
Hi To,

Le 22/12/2016 ? 12:55, To Duc Khanh a ?crit :
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.