I have rewritten the example given in the package RcppParallel
for the function ParallelReduce
, to understand how it works. Instead of summing a RVector
to a double, I've tried to sum each row of a RMatrix
to obtain a RVector
. It seems to work, but only on very small matrix ; actually, with a 10x15 matrix, it just gets wrong and I don't know why.
// [[Rcpp::depends(RcppParallel)]]
#include <Rcpp.h>
#include <RcppParallel.h>
using namespace Rcpp;
using namespace RcppParallel;
struct Sum : public Worker
{
const RMatrix<double> input;
// accumulated values
RVector<double> value;
// constructors.
//Init to zeros in {} is probably useless, anyway it doesn't solve my problem.
Sum(const NumericMatrix input, NumericVector value) : input(input), value(value) {
for (std::size_t i = 0; i < input.nrow(); ++i) {
value[i] = 0.0;
}
}
Sum(const Sum& sum, Split) : input(sum.input), value(sum.value) {
for (std::size_t i = 0; i < input.nrow(); ++i) {
value[i] = 0.0;
}
}
// accumulate just the elements of the range I've been asked to
void operator()(std::size_t begin, std::size_t end) {
for (std::size_t j = begin; j < end; ++j) {
for(std::size_t i = 0; i < input.nrow(); ++i) {
RMatrix<double>::Row ligne = input.row(i);
value[i] += std::accumulate(ligne.begin() + begin, ligne.begin() + end, 0.0);
}
}
}
// join my values with that of another Sum
void join(const Sum& rhs) {
for(std::size_t i = 0; i < input.nrow(); ++i) {
value[i] += rhs.value[i];
}
}
};
// [[Rcpp::export]]
NumericVector parallelVectorSum(NumericMatrix x) {
NumericVector resultat(x.nrow());
// declare the SumBody instance
Sum sum(x, resultat);
// call parallel_reduce to start the work
parallelReduce(0, x.ncol(), sum);
// return the computed sum
return resultat;
}
It seems ok
parallelVectorSum(matrix(1, nrow=5, ncol = 10))
[1] 10 10 10 10 10
But it isn't. I was expecting a simple vector of 50 times 100.
parallelVectorSum(matrix(1, nrow=50, ncol = 100))
[1] 4096 4096 4096 4096 2048 2048 6144 6144 6144 6144 6144 6144 6144 6144 6144 6144 10240 8192
[19] 6144 8192 8192 8192 8192 4096 11264 11264 13312 13312 13312 13312 13312 11264 11264 11264 11264 10240
[37] 10240 10240 10240 5120 7168 7168 7168 7168 5120 6144 6144 5632 9472 9472
My understanding of Rmatrix::Row is probably erroneous.
- What is the flaw here ?
- How could I access safely to a row in a matrix shared by all the thread ?