1

I need to compute singular value decomposition (SVD) of an array of matrices stored in a RcppArmadillo cube. I would like to do the SVDs in parallel for each matrix and want to enable multi-threading in calculating the SVDs via OpenBLAS. For that I use the following Cpp script:

field<vec> sing_vals(S);
#pragma omp parallel for
    for(unsigned s=0;s<S;++s){
      svd( sing_vals(s), X.slice(s) );
    }

Now before compiling this Cpp file from R via sourceCpp, I set the number of OpenMP and OpenBLAS threads with the following script:

library(OpenMPController )
omp_set_num_threads(nomp)

Sys.setenv("PKG_CXXFLAGS"="-fopenmp")
Sys.setenv("PKG_LIBS"="-fopenmp")


#########################set OPENBLAS nthreads
require(inline)
openblas.set.num.threads <- cfunction( signature(ipt="integer"),
                                       body = 'openblas_set_num_threads(*ipt);',
                                       otherdefs = c ('extern void openblas_set_num_threads(int);'),
                                       libargs = c ('-L/opt/openblas/lib -lopenblas'),
                                       language = "C",
                                       convention = ".C"
)
openblas.set.num.threads(nblas)

I have to set either of nblas and nomp to 1. Else I get the message 

OpenBLAS Warning : Detect OpenMP Loop and this application may hang. Please rebuild the library with USE_OPENMP=1 option

How can I do this nested parallelization in this exact problem? I read some materials on nested parallelization but could not figure out how to do it in my case.

noirritchandra
  • 103
  • 9
  • 1
    Short answer is "yes you likely can (with care)". Medium length answer is "please do not use `inline::cfunction` but Rcpp Attributes. Second medium term answer is try to search CRAN, other packages combine both too. – Dirk Eddelbuettel Jun 19 '20 at 16:01
  • 1
    And, as always, before you to OpenMP from R ensure that you _copy the data into non-R data structures_. But default RcppArmadillo uses zero-copy meaning you are using R memory blobs which is a no-no with OpenMP as R may do a `gc()` while you're at work. – Dirk Eddelbuettel Jun 19 '20 at 16:02
  • @DirkEddelbuettel What else should I use instead of ``inline::cfunction``? My main concern is the nested parallelization. How can I do that? – noirritchandra Jun 19 '20 at 16:30
  • 2
    Please read https://cloud.r-project.org/web/packages/Rcpp/vignettes/Rcpp-attributes.pdf -- Rcpp Attributes and its combo of `evalCpp()`, `cppFunction()` (which is what you should use) and `sourceCpp()` are preferred since 2013 (!!). See the relevant stories at the Rcpp Gallery on how to turn the OpenMP plugin in via Attributes. – Dirk Eddelbuettel Jun 19 '20 at 16:42
  • Where did your OpenBLAS library come from? Do you have alternates? On Debian and Ubuntu, I can pick between OpenBLAS built with pthread and OpenBLAS built with OpenMP (as well as a serial variant). There was a recent bug in the Ubuntu 20.04 builds which is why you may see the message. (The bug is fixed at Ubuntu too.) All that is of course a level down from RcppArmadillo and really more of a sys.admin question. – Dirk Eddelbuettel Jun 21 '20 at 15:23
  • @DirkEddelbuettel I am not sure how to check that. If you can help me find it, I can check. – noirritchandra Jun 22 '20 at 16:27
  • Back to my first question: How did you install OpenBLAS? I am mostly familiar with the binary package(s) for Debian and Ubuntu which _give me this very choice_ of switching between pthread and OpenMP backends. – Dirk Eddelbuettel Jun 22 '20 at 16:30
  • @DirkEddelbuettel I did not install OpenBLAS or OpenMP. I am running codes in my university servers and the system admins installed those libraries and the OS is Fedora. – noirritchandra Jun 23 '20 at 17:43

0 Answers0