OpenMP 点积和指针答案

【问题标题】：OpenMP Dot Product and PointersOpenMP 点积和指针
【发布时间】：2015-07-25 00:56:46
【问题描述】：

我正在尝试在 OpenMP 中使用 malloc 分配的大型数组来实现 dotproduct。但是，当我使用 reduction(+:result) 时，它会为每个程序运行产生不同的结果。为什么我会得到不同的结果？我该如何补救？以及如何优化这个例子？这是我的代码：

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <omp.h>

const int N = 1e1;

int main ()
{
  int    i, nthreads, tid;
  double x_seq, x_par, *y, *z, cpu_time_used;
  clock_t start, end;

  y = (double*)malloc(sizeof(double)*N);
  z = (double*)malloc(sizeof(double)*N);

  for (i=0; i<N; i++) {
      y[i] = i * 1.0;
      z[i] = i * 2.0;
  }

  x_seq = 0;
  x_par = 0;
  for (i=0; i<N; i++) x_seq += y[i] * z[i];

  #pragma omp parallel shared(y, z) private(i, tid)
  {
      #pragma omp single
      {
          nthreads = omp_get_num_threads();
      }
      tid = omp_get_thread_num();

      #pragma omp parallel for reduction(+:x_par)
      for (i=tid; i<N; i+=nthreads)
      {
          x_par += y[i] * z[i];
      }

  }
  return 0;
}

【问题讨论】：

请用您使用的编程语言标记您的问题。（我假设这是 C？）另外，请注意，如果这是 C（不是 C++），那么 don't cast the return value of malloc().
回答你的问题how can this example be optimized?答案是it's memory bandwidth bound so it won't scale with the number of threads though that does not mean multiple threads won't help。
@ace，这是有效的 C/C++ 代码，但不强制转换会使 C++ 代码无效。关于这个问题，你有什么有用的要说的吗？
@Zboson 检查修订历史。第一个版本没有标记语言，所以我要求 OP 相应地标记它。而且由于我注意到 OP 使用的是<stdlib.h> 而不是<cstdlib>，所以我猜它是 C 而不是 C++，所以我补充说，如果这是 C，它不应该被强制转换。如果你仍然认为我做了什么，请解释更多错了。此外，现在它被标记为 C，它使我之前的评论更加有用。
@ace 投不投是见仁见智的事，没有对错之分。

标签： c pointers for-loop openmp reduction

【解决方案1】：

这里有几个问题。

让我们看看当前的循环：

#pragma omp parallel shared(y, z) private(i, tid)
{
  #pragma omp single
  {
      nthreads = omp_get_num_threads();
  }
  tid = omp_get_thread_num();

  #pragma omp parallel for reduction(+:x_par)
  for (i=tid; i<N; i+=nthreads)
  {
      x_par += y[i] * z[i];
  }
}

所以 (1) 请注意，您（大概）希望在该区域之外可以访问 x_par。所以你需要reduction(+:x_par) 在外部，而不是内部。如果您还添加了非常有用的default(none) 子句，您还会发现没有描述nthreads 共享的子句；让我们明确地共享它。

让我们再看一遍：

#pragma omp parallel shared(y, z, nthreads) private(i, tid) reduction(+:x_par) default(none)
{
  #pragma omp single
  {
      nthreads = omp_get_num_threads();
  }
  tid = omp_get_thread_num();

  #pragma omp parallel for 
  for (i=tid; i<N; i+=nthreads)
  {
      x_par += y[i] * z[i];
  }
}

所以仔细观察，我们现在看到您有两个 omp parallel 部分。这意味着，如果启用了嵌套并行性，您将拥有 nthreads 任务，每个任务都会启动 nthreads 任务来执行该循环；因此，如果一切正常，循环最终会得到 nthreads 次正确答案。所以让我们摆脱并行，只使用for：

 #pragma omp parallel shared(y, z, nthreads) private(i, tid) reduction(+:x_par) default(none)
{
  #pragma omp single
  {
      nthreads = omp_get_num_threads();
  }
  tid = omp_get_thread_num();

  #pragma omp for 
  for (i=tid; i<N; i+=nthreads)
  {
      x_par += y[i] * z[i];
  }
}

所以共享正确并且不是嵌套并行性，但它仍然没有给出正确的答案；它给出的结果要小得多。怎么了？让我们看一下 for 循环。每个线程都想从 tid 开始并跳过 nthreads，很好；但是为什么我们要omp for 呢？

让我们看一个更简单的版本：

#pragma omp parallel shared(y, z) reduction(+:x_par) default(none)
{
  #pragma omp for
  for (i=0; i<N; i++)
  {
      x_par += y[i] * z[i];
  }
}

请注意，这里我们没有使用 tid 和 nthreads 显式分解循环 - 我们不必这样做，因为 omp for 为我们分解了循环；它将循环迭代分配给线程。

所以回顾我们所拥有的，我们对循环进行了手动分解——这很好，有时这就是你需要做的；和一个omp for，它试图获取该循环并将其拆分到多个线程中。但我们已经在这样做了； omp for 只是让我们在这里跳过迭代！

所以摆脱omp for

#pragma omp parallel shared(y, z, nthreads) private(i, tid) reduction(+:x_par) default(none)
{
  #pragma omp single
  {
      nthreads = omp_get_num_threads();
  }
  tid = omp_get_thread_num();

  for (i=tid; i<N; i+=nthreads)
  {
      x_par += y[i] * z[i];
  }
}

为我们提供正确的答案。

【讨论】：