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)
--
Khanh To Duc
PhD student, Department of Statistical Science, XXIX cycle.
Office: Room 130
Telephone: +39 049 827 4147
Mobile: +39 345 667 0711
http://www.stat.unipd.it/it/fare-ricerca/duc-khanh
_______________________________________________
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
