<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">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.</div><div dir="ltr"><br></div><div dir="ltr">  1. We have an empty matrix with the number of rows equal to numRows</div><div dir="ltr">  2. We are able to generate the ith row of the matrix at will</div><div dir="ltr">  3. Our underlying subroutine populates the matrix from any particular starting point one by one.</div><div dir="ltr"><br></div><div dir="ltr">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.</div><div dir="ltr"><br></div><div dir="ltr">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):</div><div dir="ltr"><br></div><div>************************** Code Section Start ***************************************</div><div dir="ltr"><br></div><div dir="ltr">#include <Rcpp.h></div><div dir="ltr">#include <thread></div><div dir="ltr"><br></div><div dir="ltr">// [[Rcpp::plugins(cpp11)]]</div><div dir="ltr"><br></div><div dir="ltr">int myFactorial(int n) {</div><div dir="ltr">    int res = 1;</div><div dir="ltr">    for (int i = 1; i <= n; ++i)</div><div dir="ltr">        res *= i;</div><div dir="ltr">    </div><div dir="ltr">    return res;</div><div dir="ltr">}</div><div dir="ltr"><br></div><div dir="ltr">std::vector<int> nthPerm(int n, int index) {</div><div dir="ltr">    int temp = myFactorial(n);</div><div dir="ltr">    std::vector<int> indexVec(n), res(n);</div><div dir="ltr">    std::iota(indexVec.begin(), indexVec.end(), 1);</div><div dir="ltr">    </div><div dir="ltr">    for (int k = 0, r = n; k < n; ++k, --r) {</div><div dir="ltr">        temp /= r;</div><div dir="ltr">        int j = (int) index / temp;</div><div dir="ltr">        res[k] = indexVec[j];</div><div dir="ltr">        index -= (temp * j);</div><div dir="ltr">        indexVec.erase(indexVec.begin() + j);</div><div dir="ltr">    }</div><div dir="ltr">    </div><div dir="ltr">    return res;</div><div dir="ltr">}</div><div dir="ltr"><br></div><div dir="ltr">// Simplified version for demonstration only. The real subroutines</div><div dir="ltr">// that carry out this task are more complicated</div><div dir="ltr">void PopulateMatrix(Rcpp::IntegerMatrix permuteMatrix,</div><div dir="ltr">                    std::vector<int> z, int count, int nRows, int n) {</div><div dir="ltr">    for (; count < nRows; ++count) {</div><div dir="ltr">        for (int j = 0; j < n; ++j)</div><div dir="ltr">            permuteMatrix(count, j) = z[j];</div><div dir="ltr">        </div><div dir="ltr">        std::next_permutation(z.begin(), z.end());</div><div dir="ltr">    }</div><div dir="ltr">}</div><div dir="ltr"><br></div><div dir="ltr">// [[Rcpp::export]]</div><div dir="ltr">SEXP ParallelPerms(int n, int userThrds = 0) {</div><div dir="ltr">    int nThreads = std::thread::hardware_concurrency() - 1;</div><div dir="ltr">    std::vector<std::thread> myThreads;</div><div dir="ltr">    nThreads = (userThrds > 0) ? std::min(userThrds, nThreads) : nThreads;</div><div dir="ltr">    </div><div dir="ltr">    int step = 0, numRows = myFactorial(n);</div><div dir="ltr">    int stepSize = numRows / nThreads;</div><div dir="ltr">    int nextStep = stepSize;</div><div dir="ltr">    std::vector<int> z(n);</div><div dir="ltr">    std::iota(z.begin(), z.end(), 1);</div><div dir="ltr">    </div><div dir="ltr">    Rcpp::IntegerMatrix myMat = Rcpp::no_init_matrix(numRows, n);</div><div dir="ltr">    </div><div dir="ltr">    for (std::size_t j = 0; j < (nThreads - 1); ++j) {</div><div dir="ltr">        myThreads.emplace_back(PopulateMatrix, std::ref(myMat), z, step, nextStep, n);</div><div dir="ltr">        step += stepSize;</div><div dir="ltr">        nextStep += stepSize;</div><div dir="ltr">        z = nthPerm(n, step);</div><div dir="ltr">    }</div><div dir="ltr">    </div><div dir="ltr">    // Guarantee that we get all the rows. N.B. We are passing numRows</div><div dir="ltr">    // instead of nextStep... they aren't guaranteed to be the same</div><div dir="ltr">    myThreads.emplace_back(PopulateMatrix, std::ref(myMat), z, step, numRows, n);</div><div dir="ltr">    </div><div dir="ltr">    for (auto& thr: myThreads)</div><div dir="ltr">        thr.join();</div><div dir="ltr">    </div><div dir="ltr">    return myMat;</div><div dir="ltr">}</div><div dir="ltr"><br></div><div dir="ltr">************************** Code Section End ***************************************<br></div><div dir="ltr"><br></div><div dir="ltr">I have read that Rcpp objects are not thread safe as they make unpredictable calls to the garbage collector (<a href="https://github.com/RcppCore/RcppParallel/issues/17">https://github.com/RcppCore/RcppParallel/issues/17</a>), however Romain Francois states:</div><div dir="ltr"><br></div><div dir="ltr">"As soon as you don't use references for Rcpp types, you are not safe.</div><div dir="ltr">If you use references, it all depends on what you do with them."</div><div dir="ltr"><br></div><div dir="ltr">I have a couple of questions regarding this.</div><div dir="ltr"><br></div><div dir="ltr">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?</div><div dir="ltr"><br></div><div dir="ltr">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.</div><div dir="ltr"><br></div><div dir="ltr">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.</div><div dir="ltr"><br></div><div dir="ltr">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 <a href="https://en.cppreference.com/w/cpp/thread/thread/thread">https://en.cppreference.com/w/cpp/thread/thread/thread</a>). 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).</div><div dir="ltr"><br></div><div dir="ltr">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.</div><div dir="ltr"><br></div><div dir="ltr">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?</div><div dir="ltr"><br></div><div dir="ltr">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.</div><div dir="ltr"><br></div><div dir="ltr">Thanks,</div><div dir="ltr">Joseph Wood</div></div></div></div></div>