c - OpenMP 点积和指针

Question

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

#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 并行 部分。这意味着，如果启用了嵌套并行性，您将拥有 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 开始并跳过 n 个线程，没问题；但为什么我们需要 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];
  }
}

为我们提供正确答案。