C++ Armadillo 和 OpenMp：外积求和的并行化 - 定义 Armadillo 矩阵的约简

Posted 2023-02-17

【中文标题】C++ Armadillo 和 OpenMp：外积求和的并行化 - 定义 Armadillo 矩阵的约简

【英文标题】：C++ Armadillo and OpenMp: Parallelization of summation of outer products - define reduction for Armadillo matrix

【发布时间】：2017-07-09 05:32:31

【问题描述】：

我正在尝试使用对犰狳矩阵求和的 OpenMP 并行化 for 循环。我有以下代码：

#include <armadillo>
#include <omp.h>

int main()


        arma::mat A = arma::randu<arma::mat>(1000,700);
        arma::mat X = arma::zeros(700,700);
        arma::rowvec point = A.row(0);

        # pragma omp parallel for shared(A) reduction(+:X)
        for(unsigned int i = 0; i < A.n_rows; i++)
                arma::rowvec diff = point - A.row(i);
                X += diff.t() * diff; // Adding the matrices to X here
        

我收到此错误：

[Legendre@localhost ~]$ g++ test2.cpp -o test2 -O2 -larmadillo -fopenmp
test2.cpp: In function ‘int main()’:
test2.cpp:11:52: error: user defined reduction not found for ‘X’

我阅读了有关定义约简的内容，但我还没有找到使用犰狳矩阵的示例。在我的案例中，定义犰狳矩阵减少的最佳方法是什么？

【参考方案1】：

这些缩减仅适用于内置类型（double、int 等）。因此，您必须自己进行还原，这很简单。只需将每个线程的结果累积在一个线程局部变量中，然后将其添加到临界区中的全局结果中。

#include <armadillo>
#include <omp.h>

int main()


  arma::mat A = arma::randu<arma::mat>(1000,700);
  arma::mat X = arma::zeros(700,700);
  arma::rowvec point = A.row(0);

  #pragma omp parallel shared(A)
  
    arma::mat X_local = arma::zeros(700,700);

    #pragma omp for
    for(unsigned int i = 0; i < A.n_rows; i++)
    
      arma::rowvec diff = point - A.row(i);
      X_local += diff.t() * diff; // Adding the matrices to X here
    

    #pragma omp critical
    X += X_local;
  

使用更新的 OpenMP（我认为是 4.5？），您还可以为您的类型声明用户定义的归约。

#include <armadillo>
#include <omp.h>

#pragma omp declare reduction( + : arma::mat : omp_out += omp_in ) \
  initializer( omp_priv = omp_orig )

int main()


  arma::mat A = arma::randu<arma::mat>(1000,700);
  arma::mat X = arma::zeros(700,700);
  arma::rowvec point = A.row(0);

  #pragma omp parallel shared(A) reduction(+:X)
  for(unsigned int i = 0; i < A.n_rows; i++)
  
    arma::rowvec diff = point - A.row(i);
    X += diff.t() * diff; // Adding the matrices to X here
  

【讨论】：

谢谢。为什么在添加到 X_local 之前不需要关键子句？我也尝试了第二个答案，它给了我这个错误（我知道我的矩阵的尺寸没有问题）：error: addition: incompatible matrix dimensions: 0x0 and 700x700 terminate called after throwing an instance of 'std::logic_error' what(): addition: incompatible matrix dimensions: 0x0 and 700x700 terminate called recursively Aborted (core dumped)

@Cauchy (1) X_local 对于每个线程都是本地的。指令 #pragma omp for 将 for 循环拆分为每个线程的分区，因此没有两个线程写入同一个 X_local。 (2) 是的，这是一个错误。我忘记了初始化程序，矩阵默认初始化为 0x0。查看更新后的答案。

@Cauchy 默认情况下，外部作用域中的所有变量都作为共享进入并行区域（循环变量除外），因此您不必显式标记任何共享。