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)