<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">I'm not sure I understand the problem. You've found RcppParallel, have you looked at the examples? There's one with handling matrices here:<div><br></div><div><a href="http://gallery.rcpp.org/articles/parallel-matrix-transform/">http://gallery.rcpp.org/articles/parallel-matrix-transform/</a><br></div><div><br></div><div>In the call to parallelFor() you can change the values you iterate over and change your logic based on it.</div><div>In my package I do a lot of cross-distance matrix calculations with RcppParallel, see e.g.:</div><div><br></div><div><a href="https://github.com/asardaes/dtwclust/blob/master/src/distmat/fillers.cpp">https://github.com/asardaes/dtwclust/blob/master/src/distmat/fillers.cpp</a><br></div><div><a href="https://github.com/asardaes/dtwclust/blob/master/CONTRIBUTING.md">https://github.com/asardaes/dtwclust/blob/master/CONTRIBUTING.md</a><br></div><div><br></div><div>Regards,</div><div>Alexis.</div><div><br></div></div></div></div></div><br><div class="gmail_quote"><div dir="ltr">On Wed, Sep 26, 2018 at 6:37 PM Ralf Stubner <<a href="mailto:ralf.stubner@daqana.com">ralf.stubner@daqana.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div 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" target="_blank">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="m_-5270309003638191281AppleMailSignature"><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt">— </p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">Ralf Stubner<u></u><u></u></span></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">Senior Software Engineer / Trainer</span></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><u></u> <u></u></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">daqana GmbH<u></u><u></u></span></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><font color="#000000"><span style="background-color:rgba(255,255,255,0)"><a dir="ltr">Dortustraße 48</a><u></u><u></u></span></font></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><font color="#000000"><span style="background-color:rgba(255,255,255,0)"><a dir="ltr">14467 Potsdam</a><u></u><u></u></span></font></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><u></u> <u></u></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">T: +49 331 23 70 81 66<u></u><u></u></span></p><p class="m_-5270309003638191281MsoPlainText" 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="m_-5270309003638191281MsoPlainText" 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" target="_blank">+49 162 20 91 196</a><u></u><u></u></span></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><u></u> <u></u><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" target="_blank">ralf.stubner@r-institute.com</a></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><u></u> <u></u></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">Sitz: Potsdam<u></u><u></u></span></p><p class="m_-5270309003638191281MsoPlainText" 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="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">Ust.-IdNr: DE300072622<u></u><u></u></span></p><p class="m_-5270309003638191281MsoPlainText" style="margin:0cm 0cm 0.0001pt"><span style="background-color:rgba(255,255,255,0)">Geschäftsführer: Prof. Dr. Dr. Karl-Kuno Kunze<u></u><u></u></span></p><p class="MsoNormal" style="margin:0cm 0cm 0.0001pt"><u></u> <u></u></p></div><div><br>Am 26.09.2018 um 18:17 schrieb Joseph Wood <<a href="mailto:jwood000@gmail.com" target="_blank">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" target="_blank">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" target="_blank">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" target="_blank">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" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a></span><br></div></blockquote></div></div>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" rel="noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a></blockquote></div>