Segmentation Fault/Memory Not Mapped

[linliu-stats]

Hi all, I'm trying to use RcppParallel to parallelize an iteration process iterating over all possible k (k = 2, 3, 4...) permutations over an n-dimensional integer vector. Since the number of iterations is huge, I'm trying to use RcppParallel to speed up my process. I test my code with small n/k, n = 9 and k = 3, so there's only 504 iterations in total. When stopping the iterations at 10/20, my code runs fine and the result is the same as the non-parallel version. However, when i tried to increase the iteration, I often ran into segmentation fault/memory not mapped/bad value errors. I'm not an expert on parallel computing and only started to learn such technique. Do you have any idea on how I should change my code so that the segmentation error can be avoided?
Thanks in advance!

[eddelbuettel]

No reproducible code here so how do think we should be able to help you?

[linliu-stats]

Thanks for your timely response. Sorry about this oversight. Here's the code that gives me trouble
`// [[Rcpp::depends(RcppArmadillo, RcppParallel)]]
#include <RcppArmadillo.h>
#include <RcppParallel.h>
using namespace std;
using arma::mat;
using arma::vec;
using arma::rowvec;

struct Parallel : public RcppParallel::Worker
{
const int& m;
const mat& Omega;
const mat& Precision;
const mat& BS;
const mat& all_indices;
double value;

Parallel(const int& m, const mat& Omega, const mat& Precision, const mat& BS, const mat& all_indices):
m(m), Omega(Omega), Precision(Precision), BS(BS), all_indices(all_indices), value(0) {}

Parallel(const Parallel& parallel, RcppParallel::Split):
m(parallel.m), Omega(parallel.Omega), Precision(parallel.Precision), BS(parallel.BS), all_indices(parallel.all_indices), value(0) {}

void operator()(size_t begin, size_t end)
{
for (int i = (int) begin; i != (int) end; ++i)
{
    Rcpp::NumericVector ind = Rcpp::wrap(all_indices.rows(i, i).t());
    mat kernel(1, 1);
    kernel = - pow(-1, (m - 1)) * BS.rows(ind[0], ind[0]) * Precision;
    for (unsigned s = 2; s != m; ++s)
    {
	kernel *= (BS.rows(ind[s], ind[s]).t() * BS.rows(ind[s], ind[s]) - Omega) * Precision;
    }
    kernel *= BS.rows(ind[1], ind[1]).t();
    value += kernel(0, 0);
}
}

void join(const Parallel& rhs)
{
value += rhs.value;
}

};

// [[Rcpp::export]]

double C_Parallel(int m, mat Omega, mat Precision, mat BS, mat all_indices)
{
size_t K = all_indices.n_rows;
Parallel parallel(m, Omega, Precision, BS, all_indices);
parallelReduce(0, K, parallel);
return (parallel.value / K);
}
And the R code for simulating data to run the C_Parallel function:library(Rcpp)
library(RcppParallel)
library(RcppArmadillo)
library(splines)
library(gtools)

sourceCpp('test.cpp')

set.seed(42)
n <- 9

X <- runif(n)

knots <- seq(0, 1, length.out = 4)
knots <- knots[c(-1, -length(knots))]
CovB <- bs(X, knots = knots, degree = 2, intercept = TRUE, Boundary.knots = c(0, 1))
CovB <- as.matrix(CovB)

Omega <- matrix(0, nrow = ncol(CovB), ncol = ncol(CovB))
for (i in 1:nrow(CovB)) {
Omega <- Omega + CovB[i, ] %*% t(CovB[i, ]) / nrow(CovB)
}

Precision <- solve(Omega)

m <- 3

all.indices <- permutations(n, m) - 1
C_Parallel(m, Omega, Precision, CovB, all.indices)

`

[linliu-stats]

And running the code above in R-3.3.0 gave me the error message below:
Error: bad value

Error: bad value

invalid second argument

*** caught bus error ***
address 0x10c380938, cause 'non-existent physical address'

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace

R(89195,0x700001810000) malloc: *** error for object 0x7fb3b06a05e0: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug
*** caught bus error ***
address 0x7fff920b01fa, cause 'non-existent physical address'

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace

*** caught segfault ***
Error: bad value
Error: bad value
R(89195,0x700000804000) malloc: *** error for object 0x7fb3b06a05e0: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug
Abort trap: 6

[eddelbuettel]

The documentation stresses you should not use Rcpp::NumericVector and friends, but rather RVector and RMatrix. R itself is not multi-threaded so you need to keep data at arm's length. That is likely your cause.

Also note the Styling with Markdown is supported link below to learn about formatting here to make it easier to read your code.

I will close this now. If you can substantiate an error while using the package as documented / recommended please feel free to reopen.

[jjallaire]

See the docs thread safety: http://rcppcore.github.io/RcppParallel/#thread_safety You can't call the R API inside your parallel code.

…

On Fri, Jul 14, 2017 at 1:04 PM, Dirk Eddelbuettel ***@***.*** > wrote: Closed #47 <#47>. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#47 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAGXx7KcuD-YceC199CZJqQ1qWOMJfJxks5sN595gaJpZM4OYbfp> .

[linliu-stats]

Thank you both. I just have a follow-up question - does that mean that I couldn't use all the matrix operation available in the RcppArmadillo if I'd like to use RcppParallel?

[dselivanov]

Nope. You can. Main issue in threaded code is R's memory allocation. 15 июл. 2017 г. 4:08 ПП пользователь "mamba-stats" <notifications@github.com> написал: Thank you both. I just have a follow-up question - does that mean that I couldn't use all the matrix operation available in the RcppArmadillo if I'd like to use RcppParallel? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#47 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AE4u3XyWGUbwWeI0cmJFv7Rxnec26uKcks5sOLnmgaJpZM4OYbfp> .

[linliu-stats]

Thanks! Do you mind giving me more specifics? I tried to connect RMatrix with arma::mat, but couldn't succeed. Or do you mind pointing to me a more detailed document on RcppParallel? I'm a rookie in this business and need to learn RcppParallel from bottom up. Unfortunately I couldn't find any easy-to-understand tutorial on how to use RcppParallel together with RcppArmadillo, which is quite essential to what I'd like to do. Thanks again!

[eddelbuettel]

Use the explicit constructors which create new instances rather than reusing R's memory.

[linliu-stats]

Thanks for the suggestion. I've tried something along this line. However, the memory issue is still there.
For instance

    #include <RcppArmadillo.h>
    #include <RcppParallel.h>
    using namespace std;
    using namespace Rcpp;
    using namespace arma;
    using namespace RcppParallel;
    using namespace RcppArmadillo;
    
    struct Parallel : public Worker
    {
        const RMatrix<double> Precision;
        const RMatrix<double> BS;
        const RMatrix<double> all_indices;
        double value;
        
        Parallel(const NumericMatrix Precision, const NumericMatrix BS, const NumericMatrix all_indices):
            Precision(Precision), BS(BS), all_indices(all_indices), value(0) {}

        Parallel(const Parallel& parallel, RcppParallel::Split):
            Precision(parallel.Precision), BS(parallel.BS), all_indices(parallel.all_indices), value(0) {}
        
        void operator()(size_t begin, size_t end)
        {
            mat all_indices_mat(all_indices.begin(), all_indices.nrow(), all_indices.ncol());
            mat Precision_mat(Precision.begin(), Precision.nrow(), Precision.ncol());
            mat BS_mat(BS.begin(), BS.nrow(), BS.ncol());
            for (int i = (int) begin; i != (int) end; ++i)
            {
                NumericVector ind = wrap(all_indices_mat.rows(i, i).t());
                mat kernel(1, 1);
                kernel = BS_mat.rows(ind[0], ind[0]) * Precision_mat;
                kernel *= BS_mat.rows(ind[1], ind[1]).t();
                value += kernel(0, 0);
            }
        }
        
        void join(const Parallel& rhs)
        {
            value += rhs.value;
        }
    };```


My guess is that the armadillo objects are still killing the memory with such format. Do I misunderstand your point on how to do this with Armadillo? Is there a better way to do this?

[dselivanov]

Don't use R's API inside threaded code. In particular NumericVector.