I haven’t read your message in detail, but the second example from here might be helpful: https://cran.r-project.org/web/packages/dqrng/vignettes/parallel.html
There a matrix is filled in parallel by splitting the columns among the threads. Splitting by columns is helpful since matrices in R are stored that way. You could use the same method to obtain the transpose of your matrix. Greetings Ralf — Ralf Stubner Senior Software Engineer / Trainer daqana GmbH Dortustraße 48 14467 Potsdam T: +49 331 23 70 81 66 F: +49 331 23 70 81 67 M: +49 162 20 91 196 Mail: ralf.stub...@r-institute.com Sitz: Potsdam Register: AG Potsdam HRB 27966 P Ust.-IdNr: DE300072622 Geschäftsführer: Prof. Dr. Dr. Karl-Kuno Kunze > Am 26.09.2018 um 18:17 schrieb Joseph Wood <jwood...@gmail.com>: > > As the subject states, this question is regarding populating a matrix > in parallel. I am currently reading "C++ Concurrency in Action: > Practical Multithreading" as I'd like to take some algorithms I have > to the next level. I have looked at the RcppParallel package, but the > features offered there do not seem to apply to this situation. I will > explain my reasoning further down. First, we set up our scenario. > > 1. We have an empty matrix with the number of rows equal to numRows > 2. We are able to generate the ith row of the matrix at will > 3. Our underlying subroutine populates the matrix from any > particular starting point one by one. > > This scenario easily extends to a parallel setup. We have each entry > in our matrix being visited exactly one time by only one thread. The > idea is that if we have m threads, we can split up the work so that > each thread is responsible for populating roughly (numRows / m) number > of rows of our matrix. > > Here is a simplified example that represents my real situation (In my > project I don't have the cpp11 plugin as I take care of this in > Makevars file with CXX_STD = CXX11): > > #include <Rcpp.h> > #include <thread> > > // [[Rcpp::plugins(cpp11)]] > > int myFactorial(int n) { > int res = 1; > for (int i = 1; i <= n; ++i) > res *= i; > > return res; > } > > std::vector<int> nthPerm(int n, int index) { > int temp = myFactorial(n); > std::vector<int> indexVec(n), res(n); > std::iota(indexVec.begin(), indexVec.end(), 1); > > for (int k = 0, r = n; k < n; ++k, --r) { > temp /= r; > int j = (int) index / temp; > res[k] = indexVec[j]; > index -= (temp * j); > indexVec.erase(indexVec.begin() + j); > } > > return res; > } > > // Simplified version for demonstration only. The real subroutines > // that carry out this task are more complicated > void PopulateMatrix(Rcpp::IntegerMatrix permuteMatrix, > std::vector<int> z, int count, int nRows, int n) { > for (; count < nRows; ++count) { > for (int j = 0; j < n; ++j) > permuteMatrix(count, j) = z[j]; > > std::next_permutation(z.begin(), z.end()); > } > } > > // [[Rcpp::export]] > SEXP ParallelPerms(int n, int userThrds = 0) { > int nThreads = std::thread::hardware_concurrency() - 1; > std::vector<std::thread> myThreads; > nThreads = (userThrds > 0) ? std::min(userThrds, nThreads) : nThreads; > > int step = 0, numRows = myFactorial(n); > int stepSize = numRows / nThreads; > int nextStep = stepSize; > std::vector<int> z(n); > std::iota(z.begin(), z.end(), 1); > > Rcpp::IntegerMatrix myMat = Rcpp::no_init_matrix(numRows, n); > > for (std::size_t j = 0; j < (nThreads - 1); ++j) { > myThreads.emplace_back(PopulateMatrix, std::ref(myMat), z, > step, nextStep, n); > step += stepSize; > nextStep += stepSize; > z = nthPerm(n, step); > } > > // Guarantee that we get all the rows. N.B. We are passing numRows > // instead of nextStep... they aren't guaranteed to be the same > myThreads.emplace_back(PopulateMatrix, std::ref(myMat), z, step, > numRows, n); > > for (auto& thr: myThreads) > thr.join(); > > return myMat; > } > > I have read that Rcpp objects are not thread safe as they make > unpredictable calls to the garbage collector > (https://github.com/RcppCore/RcppParallel/issues/17), however Romain > Francois states: > > "As soon as you don't use references for Rcpp types, you are not safe. > If you use references, it all depends on what you do with them." > > I have a couple of questions regarding this. > > My initial thought was I thought Rcpp objects were passed by reference > by default. Secondly, if this isn't the case, is it as simple as > adding an ampersand to all of my Rcpp objects in the function > parameters? > > The project that I'm implementing this in (RcppAlgos) builds fine on > win-builder as well as all of the various rhub checks with no errors > (even check_with_sanitizers() and check_with_valgrind()). When I > submitted v 2.1.0 to CRAN, there were sporadic build errors on some of > the linux platforms. By sporadic, I mean sometimes it passes, and > other times it would fail with the error : segfault from C stack > overflow. The current version (v 2.2.0) still has the argument for > parallel computing, but it doesn't do anything. It is only there for > backwards compatibility. > > When I initially submitted, it should be noted that I did not have my > matrices wrapped with std::ref in the call to create new threads, so > I'm not sure what effect this will have on those builds if I were to > submit to CRAN now. I will say that after I saw the errors on the > CRAN checks, I immediately sought to replicate them. I was successful > in extreme situations. For example, if I called the parallel > implementation thousands of times I could generate the error. I would > put these extreme tests in my tests folder and check them in a unit > test environment so as not to crash my r session. > > I then sought out a more robust solution to my situation and found > that thread function arguments are by default pass by value, and if > you have a particular variable that is expected to be passed by > reference, then you must add std::ref (See > https://en.cppreference.com/w/cpp/thread/thread/thread). I have done > this and have noticed that I can't generate the issues with the > extreme tests above. HOWEVER, If I call it say 50000 times 10 times in > a row, I can sometimes generate an issue (not necessarily the segfault > above.. most of the time it is a stack imbalance warning... the > warning you get when the number of UNPROTECTS doesn't match the number > of PROTECTS in the R C API). > > I then revisited the RcppParallel package to see if there were any > solutions there. I know that RcppParallel implements fully thread safe > objects like RMatrix<T>, however I can't see a way to set up a Worker > to populate my matrix. I guess the issue I see is that if you look at > my set-up above, we first get the starter vector with a call to > nthPerm from the parent function, then I pass this to PopulateMatrix > which proceeds to populate rows of my matrix for a given number of > rows. With the examples I've seen, there is no dependency on an > external function to get a specific entry point. I thought about > bypassing the RcppParallel::Worker altogether and simply use the > RMatrix<T> object, however I don't think this is how that object is to > be utilized. For example, I have not seen any RcppParallel examples > that preallocate an RMatrix<T> object. > > My question is, is there a way to make my current set up thread-safe? > If not, is it possible at all to simply populate an object in parallel > safely? > > Alternatively, if my question seems a bit naive, a nudge in the right > direction would be greatly appreciated. I don't mind going back to > square one, if need be. > > Thanks, > Joseph Wood > _______________________________________________ > 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
_______________________________________________ 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