[Rcpp-devel] Populate a Matrix in Parallel (Text Version)

[Ralf Stubner]

Hi Joseph,

On 9/28/18 3:26 AM, Joseph Wood wrote:
> I understand that parallelism is not a magic wand. Have you read my
> original post? I have managed to parallelize generating permutations
> by taking advantage of the fact that I can generate the ith
> permutation via nthPerm. My question is about making this thread safe
> not if it is possible.

unfortunately, you messages are a bit long and contradictory. I had the
same reaction as Jeff when I read

>> My situation is fundamentally different. The algorithm that fills the
>> matrix does so in a way that relies on the previous row and more
[...]

If you have code that generates the n-th row based on some global input,
then the example I refered to would be the right starting point. It uses
global state (seed) and column indices to fill a matrix by column. It
does it for multiple columns together, since it is more efficient to
have one thread process multiple columns. However, the default grain
size for parallelFor is one, so it is easy to create one thread per
column and perform some action based on some global input and on the
column index:


#include <Rcpp.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
// [[Rcpp::plugins(cpp11)]]

struct ParallelFill : public RcppParallel::Worker {
  RcppParallel::RMatrix<double> output;
  RcppParallel::RVector<double> global_input;

  ParallelFill(Rcpp::NumericMatrix output,
               Rcpp::NumericVector input) :
               output(output), global_input(input) {};

  std::vector<double> create_column(std::size_t index) {
    std::vector<double> result(global_input.size());
    std::transform(global_input.begin(),
                   global_input.end(),
                   result.begin(),
                   [&index] (double a) {return a + index;});
    return result;
  }

  // default grain size is 1, i.e. end == begin + 1
  void operator()(std::size_t begin, std::size_t end) {
    std::vector<double> column = create_column(begin);
    std::copy(column.begin(), column.end(),
              output.begin() + begin * output.nrow());
  }
};

// [[Rcpp::export]]
Rcpp::NumericMatrix parallel_matrix(const int n,
                                    Rcpp::NumericVector input) {
  Rcpp::NumericMatrix res(input.length(), n);
  ParallelFill parallelFill(res, input);
  RcppParallel::parallelFor(0,  n, parallelFill);
  return res;
}

/*** R
set.seed(42)
res <- parallel_matrix(8, runif(1e7))
head(res)
*/

All this is done *by column*, since matrices in R are stored that way.

cheerio
ralf


-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 833 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20180928/fb0a9e04/attachment.sig>