1

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.

  1. What is the flaw here ?
  2. How could I access safely to a row in a matrix shared by all the thread ?
user20650
  • 24,654
  • 5
  • 56
  • 91
Maxime2506
  • 11
  • 2

0 Answers0