Rcpp omp_set_num_threads внутри пакета

Question

Я написал следующий простой пример с Rcpp и OpenMP, который прекрасно работает, когда я получаю файл cpp из RStudio:

#include <Rcpp.h>
#include <omp.h>

// [[Rcpp::plugins(openmp)]]

using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix my_matrix(int I, int J, int nthreads) {
  NumericMatrix A(I,J);
  int i,j,tid;
  omp_set_num_threads(nthreads);
#pragma omp parallel for private(i, j, tid)
  for(int i = 0; i < I; i++) {
    for(int j = 0; j < J; j++) {
      tid = omp_get_thread_num();
      A(i,j) = tid ;
    }
  }

  return A;
}


/*** R
set.seed(42)
  my_matrix(10,10,5)
*/
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    1    1    1    1    1    1    1    1    1     1
 [4,]    1    1    1    1    1    1    1    1    1     1
 [5,]    2    2    2    2    2    2    2    2    2     2
 [6,]    2    2    2    2    2    2    2    2    2     2
 [7,]    3    3    3    3    3    3    3    3    3     3
 [8,]    3    3    3    3    3    3    3    3    3     3
 [9,]    4    4    4    4    4    4    4    4    4     4
[10,]    4    4    4    4    4    4    4    4    4     4

Однако тот же код не работает должным образом, если я создаю пакет:

> rcpphello::my_matrix(10,10,5)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    0    0    0    0    0    0     0
 [4,]    0    0    0    0    0    0    0    0    0     0
 [5,]    0    0    0    0    0    0    0    0    0     0
 [6,]    0    0    0    0    0    0    0    0    0     0
 [7,]    0    0    0    0    0    0    0    0    0     0
 [8,]    0    0    0    0    0    0    0    0    0     0
 [9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

Почему один и тот же код использует только один поток, если я вызываю его из пакета?На случай, если это поможет, я отправил весь код в этот репозиторий github

Answer 1

@ coatless уже ответил на вопрос.Я хотел бы добавить одну оговорку: не используйте структуры данных из R или Rcpp в параллельном коде.Вы можете использовать RcppParallel , однако:

#include <Rcpp.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix my_matrix(int I, int J, int nthreads) {
  NumericMatrix A(I,J);
  // create a thread safe accessor for A
  RcppParallel::RMatrix<double> a(A);
  int tid;
  omp_set_num_threads(nthreads);
#pragma omp parallel for private(tid)
  for(int j = 0; j < J; j++) {
    for(int i = 0; i < I; i++) {
      tid = omp_get_thread_num();
      a(i, j) = tid ;
    }
  }

  return A;
}


/*** R
set.seed(42)
my_matrix(12,10,5)
*/

Обратите внимание, что я также изменил доступ к главному столбцу и удалил дополнительное объявление i и j.Обратите внимание, что переменные, объявленные внутри секции parallel, автоматически становятся закрытыми.

И в случае, если вы хотите использовать RNG R (поскольку вы устанавливаете начальное число), есть еще одно «не делайте этого».Посмотрите пакеты типа sitmo или dqrng для RNG, которые можно использовать с параллельным кодом.

Answer 2

Добавить к src/Makevars и src/Makevars.win:

PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS)
PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS)

Включает флаг -fopenmp.В противном случае вы не будете в конечном итоге включать OpenMP в вашем пакете.

Примечание: При использовании:

// [[Rcpp::plugins(openmp)]]

Устанавливает параметр -fopenmp только при работе с sourceCpp().Эта опция не переносится в пакет.Следовательно, мы должны установить настройку в Makevars и Makevars.win.

Краткий пример можно найти здесь:

https://github.com/r-pkg-examples/rcpp-and-openmp

Хотя я будунужно немного почистить.