Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Is it allowed/possible to call an R function or fortran code within a pragma openmp parallel for loop in Rcpp?

In an Rcpp project I would like to be able to either call an R function (the cobs function from the cobs package to do a concave spline fit) or call the fortran code that it relies on (the cobs function uses quantreg's rq.fit.sfnc function to fit the constrained spline model, which in turn relies on the fortran coded srqfnc function in the quantreg package) within a pragma openmp parallel for loop (the rest of my code mainly requires some simple linear algebra, so that would be no problem, but sadly each inner loop iteration also requires me to do a concave spline fit). I was wondering if this would be allowed or possible though, as I presume such calls would not be thread safe? Would there be an easy fix for this, like to surround those calls by a #pragma omp critical ? Would anyone have any examples of this? Or would the only way in this case involve doing a full Rcpp port of the cobs & rq.fit.sfnc functions first using thread-safe Armadillo classes?

like image 919
Tom Wenseleers Avatar asked Mar 05 '23 07:03

Tom Wenseleers


1 Answers

Citing the manual:

Calling any of the R API from threaded code is ‘for experts only’ and strongly discouraged. Many functions in the R API modify internal R data structures and might corrupt these data structures if called simultaneously from multiple threads. Most R API functions can signal errors, which must only happen on the R main thread. Also, external libraries (e.g. LAPACK) may not be thread-safe.

I have always interpreted this as "one must not call an R API function from threaded code". Calling an R function, irregardless of what's used inside, from inside a omp parallel region would be just that. Using #pragma omp critical might work, but if it breaks you got to keep the pieces ...

It would be safer to either re-implement the code in question or look for an existing implementation in C++/C/Fortran and call that directly.

like image 129
Ralf Stubner Avatar answered Mar 06 '23 23:03

Ralf Stubner