c++ - 使用 OpenMP 和减少优化外循环

Question

我对一个函数有点纠结。如果我尝试用

并行化外循环，计算是错误的

#pragma omp parallel reduction(+:det). 

有人可以告诉我如何解决它以及它失败的原因吗？

// template<class T> using vector2D = std::vector<std::vector<T>>;

float Det(vector2DF &a, int n)
{
  vector2DF m(n - 1, vector1DF(n - 1, 0));

  if (n == 1) return a[0][0];
  if (n == 2) return a[0][0] * a[1][1] - a[1][0] * a[0][1];

  float det = 0;
  for (int i = 0; i < n; i++)
  {
    int l = 0;
#pragma omp parallel for private(l)
    for (int j = 1; j < n; j++)
    {
      l = 0;
      for (int k = 0; k < n; k++)
      {
        if (k == i) continue;
        m[j - 1][l] = a[j][k];
        l++;
      }
    }
    det += std::pow(-1.0, 1.0 + i + 1.0) * a[0][i] * Det(m, n - 1);
  }

  return det;
}

Answer 1

如果将外循环并行化，则此行存在竞争条件:

m[j - 1][l] = a[j][k];

您可能还需要一个parallel for reduction 而不仅仅是一个parallel reduction。

问题是，m 是共享的，尽管考虑到它在内部循环中被完全覆盖，这不是必需的。始终尽可能在本地声明变量，这避免了错误共享变量的问题，例如:

float Det(vector2DF &a, int n)
{
  if (n == 1) return a[0][0];
  if (n == 2) return a[0][0] * a[1][1] - a[1][0] * a[0][1];

  float det = 0;
  #pragma omp parallel reduction(+:det)
  for (int i = 0; i < n; i++)
  {
    vector2DF m(n - 1, vector1DF(n - 1, 0));
    for (int j = 1; j < n; j++)
    {
      int l = 0;
      for (int k = 0; k < n; k++)
      {
        if (k == i) continue;
        m[j - 1][l] = a[j][k];
        l++;
      }
    }
    det += std::pow(-1.0, 1.0 + i + 1.0) * a[0][i] * Det(m, n - 1);
  }
  return det;
}

现在这是正确的，但是由于 m 的分配成本很高，如果不在每次迭代中都这样做，性能可能会有所提高。这可以通过拆分 parallel 和 for 指令来完成:

float Det(vector2DF &a, int n)
{
  if (n == 1) return a[0][0];
  if (n == 2) return a[0][0] * a[1][1] - a[1][0] * a[0][1];

  float det = 0;
  #pragma omp parallel reduction(+:det)
  {
    vector2DF m(n - 1, vector1DF(n - 1, 0));
    #pragma omp parallel for
    for (int i = 0; i < n; i++)
    {
      for (int j = 1; j < n; j++)
      {
        int l = 0;
        for (int k = 0; k < n; k++)
        {
          if (k == i) continue;
          m[j - 1][l] = a[j][k];
          l++;
        }
      }
      det += std::pow(-1.0, 1.0 + i + 1.0) * a[0][i] * Det(m, n - 1);
    }
  }
  return det;
}

现在您也可以将 m 声明为 firstprivate，但这将假设复制构造函数进行完全独立的深度复制，从而使代码更难原因。

请注意，您应该始终包括预期输出、实际输出和 minimal complete and verifiable example .