C++11使用OpenMP并行构造距离矩阵的7条假设正确性核验

我希望在C++11中使用OpenMP并行构造距离矩阵，已查阅各类文档、教程和示例，但仍存在疑问，现将疑问整理为编号1至7的假设，请指出各假设是否正确。

稠密距离矩阵实现 #

串行实现 #

// [[Rcpp::export]]
arma::mat compute_dist_mat(arma::mat &coordinates, unsigned int n_points) {
  arma::mat dist_mat(n_points, n_points, arma::fill::zeros);
  double dist {};
  for(unsigned int i {0}; i < n_points; i++) {
    for(unsigned int j = i + 1; j < n_points; j++) {
      dist = compute_dist(coordinates(i, 1), coordinates(j, 1), coordinates(i, 0), coordinates(j, 0));
      dist_mat.at(i, j) = dist;
      dist_mat.at(j, i) = dist;
    }
  }
  return dist_mat;
}

头文件与编译配置 #

该函数通过Rcpp接口供R调用，文件头部包含以下配置：

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

using namespace Rcpp;
using namespace arma;

即使没有R接口，该函数也可正常运行。

并行改造代码 #

为实现并行化，新增n_threads入参，将循环替换为如下形式：

unsigned int i {};
unsigned int j {};
# pragma omp parallel for private(dist, i, j) num_threads(n_threads) if(n_threads > 1)
for(i = 0; i < n_points; i++) {
  for(j = i + 1; j < n_points; j++) {
    dist = compute_dist(coordinates(i, 1), coordinates(j, 1), coordinates(i, 0), coordinates(j, 0));
    dist_mat.at(i, j) = dist;
    dist_mat.at(j, i) = dist;
  }
}

稀疏距离矩阵实现 #

因为需要将低于指定阈值的距离设为0，因此采用稀疏矩阵存储，先将非零值存入vector再批量构造Armadillo稀疏矩阵。

串行实现 #

// [[Rcpp::export]]
arma::sp_mat compute_dist_spmat(arma::mat &coordinates, unsigned int n_points, double dist_threshold) {
  std::vector<double> dists;
  std::vector<unsigned int> dist_i;
  std::vector<unsigned int> dist_j;
  double dist {};
  for(unsigned long int i {0}; i < n_points; i++) {
    for(unsigned long int j = i + 1; j < n_points; j++) {
      dist = compute_dist(coordinates(i, 1), coordinates(j, 1), coordinates(i, 0), coordinates(j, 0));
      if(dist >= dist_threshold) {
        dists.push_back(dist);
        dist_i.push_back(i);
        dist_j.push_back(j);
      }
    }
  }
  unsigned int mat_size = dist_i.size();
  arma::umat index_mat(2, mat_size * 2);
  arma::vec dists_vec(mat_size * 2);
  unsigned int j {};
  for(unsigned int i {0}; i < mat_size; i++) {
    j = i * 2;
    index_mat.at(0, j) = dist_i[i];
    index_mat.at(1, j) = dist_j[i];
    index_mat.at(0, j + 1) = dist_j[i];
    index_mat.at(1, j + 1) = dist_i[i];
    dists_vec.at(j) = dists[i];
    dists_vec.at(j + 1) = dists[i];
  }
  arma::sp_mat dist_mat(index_mat, dists_vec, n_points, n_points);
  return dist_mat;
}

并行改造代码 #

// [[Rcpp::export]]
arma::sp_mat compute_dist_spmat(arma::mat &coordinates, unsigned int n_points, double dist_threshold, unsigned short int n_threads) {
  std::vector<std::vector<double>> dists(n_points);
  std::vector<std::vector<unsigned int>> dist_j(n_points);
  double dist {};
  unsigned int i {};
  unsigned int j {};
  # pragma omp parallel for private(dist, i, j) num_threads(n_threads) if(n_threads > 1)
  for(i = 0; i < n_points; i++) {
    for(j = i + 1; j < n_points; j++) {
      dist = compute_dist(coordinates(i, 1), coordinates(j, 1), coordinates(i, 0), coordinates(j, 0));
      if(dist >= dist_threshold) {
        dists[i].push_back(dist);
        dist_j[i].push_back(j);
      }
    }
  }
  unsigned int vec_intervals[n_points + 1];
  vec_intervals[0] = 0;
  for (i = 0; i < n_points; i++) {
    vec_intervals[i + 1] = vec_intervals[i] + dist_j[i].size();
  }
  unsigned int mat_size {vec_intervals[n_points]};
  arma::umat index_mat(2, mat_size * 2);
  arma::vec dists_vec(mat_size * 2);
  unsigned int vec_begins_i {};
  unsigned int vec_length_i {};
  unsigned int k {};
  # pragma omp parallel for private(i, j, k, vec_begins_i, vec_length_i) num_threads(n_threads) if(n_threads > 1)
  for(i = 0; i < n_points; i++) {
    vec_begins_i = vec_intervals[i];
    vec_length_i = vec_intervals[i + 1] - vec_begins_i;
    for(j = 0; j < vec_length_i; j++) {
      k = (vec_begins_i + j) * 2;
      index_mat.at(0, k) = i;
      index_mat.at(1, k) = dist_j[i][j];
      index_mat.at(0, k + 1) = dist_j[i][j];
      index_mat.at(1, k + 1) = i;
      dists_vec.at(k) = dists[i][j];
      dists_vec.at(k + 1) = dists[i][j];
    }
  }
  arma::sp_mat dist_mat(index_mat, dists_vec, n_points, n_points);
  return dist_mat;
}

7个假设的验证结果 #

• 假设1：该写法将外层循环的迭代分配给不同线程，内层循环的迭代由负责对应外层迭代的线程执行。
  结论：正确。OpenMP的parallel for默认作用于紧跟的第一层for循环，会将外层i的迭代拆分给不同线程，每个线程独立执行分配到的i对应的全部内层j循环。
• 假设2：两层循环无法合并，因为内层循环依赖外层循环的取值。
  结论：错误。可以将三角遍历的所有(i,j)对映射为全局线性索引，总共有n_points*(n_points-1)/2个计算任务，单个全局id可以反推出对应的i和j，因此两层循环可以合并为单层循环做并行。
• 假设3：dist、i、j都需要在#pragma行上方初始化，然后声明为私有变量，而非在循环内部初始化。
  结论：错误。C++支持在for循环头内定义循环变量，比如for(unsigned int i=0; i<n_points; i++)，这种场景下循环变量会自动被OpenMP识别为私有变量，无需提前定义再手动声明private。
• 假设4：当n_threads = 1时#pragma行不生效，代码串行执行。
  结论：正确。pragma语句中指定了if(n_threads > 1)的判断条件，当条件不成立时，该并行区域会自动退化为串行执行。
• 假设5：在循环中使用动态vector是线程安全的。
  结论：正确。此处采用了按行拆分的私有存储设计，每个线程仅操作dists[i]和dist_j[i]中与自己分配到的i对应的子vector，不存在多线程同时读写同一个vector的场景，因此线程安全。
• 假设6：dist、i、j、k、vec_begins_i、vec_length_i都需要在#pragma行上方初始化，然后声明为私有变量，而非在循环内部初始化。
  结论：错误。和假设3同理，这些变量都可以定义在循环内部或并行区域内部，自动成为线程私有变量，无需提前定义再手动声明私有。
• 假设7：不需要将任何内容标记为section。
  结论：正确。sections指令用于手动划分不同的独立代码块给不同线程执行，当前场景全部为循环迭代并行，仅用parallel for即可满足需求，无需使用sections。

内容的提问来源于stack exchange，提问作者user