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)
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
