OpenMP 点积和指针

Question

我正在尝试在 OpenMP 中使用分配有 malloc 的大型数组来实现点积。但是，当我使用 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;
}

Answer 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];
  }
}

给我们正确的答案。

OpenMP 点积和指针

OpenMP Dot Product and Pointers

c

pointers

for-loop

openmp

reduction