RcppParallel Распараллеливание вычисления расстояния: segfault

У меня есть матрица, для которой я хочу вычислить расстояние (скажем, евклидово) между яй ряд и каждый второй ряд (т. е. я хочу яСтрока матрицы парных расстояний).

#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;
}

Последовательным способом использования Rcpp является функция ‘row_dist’, как написано выше. Тем не менее, матрица, с которой я хочу работать, очень большая, поэтому я хочу распараллелить ее. Но тогда я столкнусь с ошибкой segfault, которую я не совсем понимаю, почему. Чтобы вызвать ошибку, вы можете запустить следующий код:

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)

Может кто-нибудь подсказать мне, что я сделал не так? Спасибо заранее за ваше время!

================================================== =====================

Редактировать:

inline double d(double a, double b){
return fabs(a - b);
}

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

// Input matrix to read from must be of the RMatrix<T> form
// if using Rcpp objects
const RMatrix<double> mat;
int i;

// Output vector to write to must be of the RVector<T> form
// if using Rcpp objects
RVector<double> output;

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

// Note the -1 ^^^^ to match results from prior function

// Function call operator to iterate over a specified range (begin/end)
void operator()(std::size_t begin, std::size_t end) {

RMatrix<double>::Row row1 = mat.row(i);

for (std::size_t j = begin; j < end; ++j) {

RMatrix<double>::Row row2 = mat.row(j);

size_t n = row1.length();
size_t m = row2.length();
NumericMatrix cost(n + 1, m + 1);

for (int ii = 1; ii <= n; ii++){
cost(i, 0) = numeric_limits<double>::infinity();
}

for (int jj = 1; jj <= m; jj++){
cost(0, j) = numeric_limits<double>::infinity();
}

for (int ii = 1; ii <= n; ii++){
for (int jj = 1; jj <= m; jj++){
double dist = d(row1[ii - 1], row2[jj - 1]);
cost(ii, jj) = dist + min(min(cost(ii - 1, jj), cost(ii, jj - 1)), cost(ii - 1, jj - 1));
//cout << ii << ", " << jj << ", " << cost(ii, jj) << "\n";
}
}
output[j] = cost(n, m);

}
}
};
// [[Rcpp::export]]
NumericVector parallel_dist_row_dtw(NumericMatrix mat, int i) {

// allocate the matrix we will return
//RMatrix<double> input(mat);
NumericVector y(mat.nrow());
//RVector<double> output(y);

// create the worker
dtwDistance dtwDistance(mat, i, y);

// call it with parallelFor
parallelFor(0, mat.nrow(), dtwDistance);
return y;
}

Расстояние, которое мне нужно было рассчитать, — это динамическое время деформации. Я реализовал это, как указано выше. Тем не менее, при запуске он выдаст предупреждение «дисбаланс стека». И после нескольких запусков будет сегфо. Мне интересно, в чем проблема сейчас.

Чтобы вызвать проблему, я сделал:

library(Rcpp)
library(RcppParallel)
setThreadOptions(numThreads = 4)
sourceCpp("scripts/chisq_dtw.cpp")
set.seed(42)
X = matrix(rnorm(1000), 100, 10)

parallel_dist_row_dtw(X, 1)
parallel_dist_row_dtw(X, 2)
parallel_dist_row_dtw(X, 3)
parallel_dist_row_dtw(X, 4)
parallel_dist_row_dtw(X, 5)