(Based on this Q: http://stackoverflow.com/q/36206524/164611 , and
available from this gist:
https://gist.github.com/shabbychef/9c2ed6b2ee92f104a8b1 )
I am trying to write some code using the parallelReduce framework from
RcppParallel. Somehow I must be violating thread safety. I reduced my real
code down to a simple MWE that does not work. I'll include it below. This
very silly example just computes the length of an input vector. I thought
it was working correctly, but noticed that if I ran it many times quickly
(via replicate in R), I would get different answers. Per Dirk's advice, I
tried to remove any hint of Rcpp functions to maintain thread safety, but I
am not sure what is left. Any help appreciated.
yunoparallel.cpp:
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;
// accumulate # of elements seen on the input iterator.
void accum_nel_int(RVector<double>::const_iterator vecit,
RVector<double>::const_iterator end,
RVector<int>::iterator xret) {
for (;vecit != end;++vecit) { *xret += 1; }
}
void join_em_int(RVector<int> lhs, const RVector<int> rhs) {
lhs[0] += rhs[0];
}
struct NEL_int : public Worker
{
// source vector
const RVector<double> input;
// accumulated value
RVector<int> output;
// constructors
NEL_int(const NumericVector input, IntegerVector output) : input(input),
output(output) {}
NEL_int(const NEL_int& othr, Split) : input(othr.input),
output(othr.output) {}
// accumulate just the element of the range I've been asked to
void operator()(std::size_t begin, std::size_t end) {
accum_nel_int(input.begin() + begin, input.begin() + end,
output.begin());
}
// join my value with that of another NEL_int
void join(const NEL_int& rhs) {
join_em_int(output,rhs.output);
}
};
// dumb example counting # of elements on input:
// [[Rcpp::export]]
IntegerVector dumbCount_int(NumericVector x) {
// allocate the output
IntegerVector output(1);
// declare the instance
NEL_int wrkr(x,output);
// call parallel_reduce to start the work
parallelReduce(0, x.length(), wrkr, 5000);
// return the computed number of non-na elements
return output;
}
//for vim modeline: (do not edit)
//
vim:ts=4:sw=4:tw=79:fdm=marker:fmr=FOLDUP,UNFOLD:cms=//%s:syn=cpp:ft=cpp:ai:si:et:cin:nu:fo=croql:cino=p0t0c5(0:
In R:
library(Rcpp)
source('yunoparallel.cpp')
x <- rnorm(1e5)
dumbCount_int(x)
dumbCount_int(x)
# fine, right? no.
replicate(10,dumbCount_int(x))
sd(replicate(100,dumbCount_int(x)))
sessionInfo:
> sessionInfo()
R version 3.2.4 Revised (2016-03-16 r70336)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=C
attached base packages:
[1] graphics grDevices datasets stats utils methods base
other attached packages:
[1] Rcpp_0.12.3 colorout_1.1-2 fortunes_1.5-2 Quandl_2.7.0
xts_0.9-7 zoo_1.7-12 devtools_1.10.0 drat_0.1.0
loaded via a namespace (and not attached):
[1] httr_1.1.0 R6_2.1.2 tools_3.2.4
memoise_1.0.0 grid_3.2.4 jsonlite_0.9.19 digest_0.6.9
RcppParallel_4.3.15
[9] lattice_0.20-33
--
--sep
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
