包内的 Rcpp omp_set_num_threads

Posted 2023-02-21

【中文标题】包内的 Rcpp omp_set_num_threads

【英文标题】：Rcpp omp_set_num_threads within a package

【发布时间】：2019-02-15 19:09:09

【问题描述】：

我使用 Rcpp 和 OpenMP 编写了以下简单示例，当我从 RStudio 获取 cpp 文件时，它可以正常工作：

#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

如果我从包中调用相同的代码，为什么只使用一个线程？如果有帮助，我将所有代码推送到this github repo

【问题讨论】：

在我看来，标准的 OpenMP 环境变量设置控制着它。当然，还有帮助程序包RhpcBLASctl。哦，我看你有电话。那么@coatless 是完全正确的。

两个警告：不要使用来自Rcpp 的数据结构进行并行处理。不过可以使用来自RcppParallel 的那些。 R 对矩阵使用列主要布局。一般来说，最好在您的代码中遵循这一点。

@ralf-stubner 你能告诉我你对RcppParallel的答案吗？我刚开始学习这些东西，很高兴看到您的答案得到实施。

【参考方案1】：

添加到src/Makevars 和src/Makevars.win：

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

这会启用-fopenmp 标志。否则，您最终不会在您的包中启用 OpenMP。

注意：使用时：

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

当使用sourceCpp() 运行时，这仅设置-fopenmp 参数。此选项不转移到包中。因此，我们必须在Makevars和Makevars.win中建立设置。

一个简短的例子可以在这里找到：

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

不过，我需要稍微清理一下。

【参考方案2】：

@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 部分中声明的变量自动是私有的。

如果你想使用 R 的 RNG（因为你正在设置种子），还有另一个“不要那样做”。查看类似 sitmo 或 dqrng 的包，这些包可以用于并行代码的 RNG。

【讨论】：

非常感谢@ralf-stubner。如果函数声明了NumericMatrix A(I,J);，您能否解释一下RcppParallel::RMatrix<double> a(A); 行在做什么以及为什么要执行a(j, i)？如果我使用您的代码 R 调用 my_matrix(12,10,5) 会崩溃

@Ignacio 该行为A 创建了一个线程安全的访问器。有关更多信息，请参阅链接的 RcppParallel 文档。感谢您检查非方阵。那确实是头脑简单。切换到列主要访问权限的问题现已修复。