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