jeudi 10 août 2017

Parallelizing distance computation: segfault

I have a matrix, for which I want to compute the distance (let's say Euclidean) between the ith row and every other row(i.e. I want the ith row of the pairwise distance matrix).

#include <Rcpp.h>
#include <cmath>
#include <algorithm>
#include <RcppParallel.h>
//#include <RcppArmadillo.h>
#include <queue>

using namespace std;
using namespace Rcpp;
using namespace RcppParallel;

// [[Rcpp::export]]
double dist_fun(NumericVector row1, NumericVector row2){
   double rval = 0;
   for (int i = 0; i < row1.length(); i++){
      rval += (row1[i] - row2[i]) * (row1[i] - row2[i]);
   }
   return rval;
}

// [[Rcpp::export]]
NumericVector dist_row(NumericMatrix mat, int i){
   NumericVector row(mat.nrow());

   NumericMatrix::Row row1 = mat.row(i - 1);
   for (int j = 0; j < mat.nrow(); j++){
      NumericMatrix::Row row2 = mat.row(j);

      row(j) = dist_fun(row1, row2);
   }
   return row;
}

// [[Rcpp::depends(RcppParallel)]]
struct JsDistance: public Worker {

   // input matrix to read from
   const NumericMatrix mat;
   int i;

   // output vector to write to
   NumericVector output;

   // initialize from Rcpp input and output matrixes (the RMatrix class
   // can be automatically converted to from the Rcpp matrix type)
   JsDistance(const NumericMatrix mat, int i, NumericVector output)
      : mat(mat), i(i), output(output) {}


   // function call operator that work for the specified range (begin/end)
    void operator()(std::size_t begin, std::size_t end) {
      NumericVector row1 = mat.row(i);
      for (std::size_t j = begin; j < end; j++) {

         NumericVector row2 = mat.row(j);
         output[j] = dist_fun(row1, row2);
      }
    }
};

// [[Rcpp::export]]
NumericVector parallel_dist_row(NumericMatrix mat, int i) {

   // allocate the matrix we will return
   NumericVector output(mat.nrow());

   // create the worker
   JsDistance JsDistance(mat, i, output);

   // call it with parallelFor
   parallelFor(0, mat.nrow(), JsDistance);
   return output;
}

The sequential way using Rcpp is the function 'row_dist' as written above. Yet the matrix I want to work with is very large so I want to parallelize it. But then I will run into a segfault error which I don't quite understand why. To trigger the error you can run the following code:

library(Rcpp)
library(RcppParallel)

setThreadOptions(numThreads = 20)


set.seed(42)
X = matrix(rnorm(10000 * 400), 10000, 400)
sourceCpp("question.cpp")


start1 = proc.time()
print(dist_row(X, 2)[1:30])
print(proc.time() - start1)

start2 = proc.time()
print(parallel_dist_row(X, 2)[1:30])
print(proc.time() - start2)

enter image description here

Can someone give me some hint about what I did wrong? Thanks in advance for your time!

Aucun commentaire:

Enregistrer un commentaire