C 使用 pragma omp 时间优化问题

Question

我正在尝试并行运行此代码，以优化“挂钟时间”

在不做任何更改的情况下运行它大约 0.00775 秒

原代码段：

for (it=1;it<=itmax;it++)
  {
    dphimax=0.;
    for (k=1;k<kmax;k++)
    {
      for (i=1;i<imax;i++)
      {
        dphi=(phi[i+1][k]+phi[i-1][k]-2.*phi[i][k])*dy2i
            +(phi[i][k+1]+phi[i][k-1]-2.*phi[i][k])*dx2i;
        dphi=dphi*dt;
        dphimax=max(dphimax,dphi);
        phin[i][k]=phi[i][k]+dphi;
      }
    }
/* save values */
    for (k=1;k<kmax;k++)
    {
      for (i=1;i<imax;i++)
      {
    phi[i][k]=phin[i][k];
      }
    }
    if(dphimax<eps) break;
  }

  t2=clock();
# ifdef _OPENMP
    wt2=omp_get_wtime();
# endif
  gettimeofday(&tv2, &tz);

  printf("\nphi after %d iterations\n",it);
  heatpr(phi);
  printf( "CPU time (clock)                = %12.4g sec\n", (t2-t1)/1000000.0 );
# ifdef _OPENMP
    printf( "wall clock time (omp_get_wtime) = %12.4g sec\n", wt2-wt1 );
# endif
  printf( "wall clock time (gettimeofday)  = %12.4g sec\n", (tv2.tv_sec-tv1.tv_sec) + (tv2.tv_usec-tv1.tv_usec)*1e-6 );

} 

当我尝试使用 openmp 进行优化时，它只会变得更糟。

挂钟时间上升到大约 0.01644 秒，我不知道我做错了什么？

这样做了：

  #pragma omp parallel private(it, k, i), shared(phi, phin, dy2i, dx2i, dphi) //tried also with shared(dphi, dphimax, phi, phin)
  {
  for (it=1;it<=itmax;it++)
  { 
    dphimax=0.;
    //pragma omp parallel for
    for (k=1;k<kmax;k++)
    {
      for (i=1;i<imax;i++)
      {
        dphi=(phi[i+1][k]+phi[i-1][k]-2.*phi[i][k])*dy2i
            +(phi[i][k+1]+phi[i][k-1]-2.*phi[i][k])*dx2i;
        dphi=dphi*dt;
        dphimax=max(dphimax,dphi);
        phin[i][k]=phi[i][k]+dphi;
      }
    }

如何优化？

Answer 1

在不做任何更改的情况下运行它大约 0.00775 秒

要真正看到并行性的好处，您需要增加输入大小。否则，并行执行的计算将不足以克服并行引入的额外开销（例如，线程创建和同步）

让我们首先声明更接近它们将使用的范围的变量。这将避免必须显式添加 OpenMP pragma 以使这些变量成为私有变量（即 private(it, k, i)）

double dphimax;
for (it=1;it<=itmax;it++){    
    dphimax=0.0;
    for (k=1;k<kmax;k++){
       for (i=1;i<imax;i++){
          dphi=(phi[i+1][k]+phi[i-1][k]-2.*phi[i][k])*dy2i
              +(phi[i][k+1]+phi[i][k-1]-2.*phi[i][k])*dx2i;
          dphi=dphi*dt;
          dphimax=max(dphimax,dphi);
          phin[i][k]=phi[i][k]+dphi;
      }
   }
    for (k=1;k<kmax;k++)
      for (i=1;i<imax;i++)
          phi[i][k]=phin[i][k];
    if(dphimax<eps) break;
}

对于并行区域，我们可以尝试添加整个代码块：

double dphimax;
#pragma omp parallel
for (it=1;it<=itmax;it++){    
    ....
}

现在让我们看看哪些循环迭代可以在线程之间划分。最外层的循环不起作用，因为它的迭代之间存在依赖关系。下一个候选者是循环：

for (k=1;k<kmax;k++){
   for (i=1;i<imax;i++){
      dphi=(phi[i+1][k]+phi[i-1][k]-2.*phi[i][k])*dy2i
          +(phi[i][k+1]+phi[i][k-1]-2.*phi[i][k])*dx2i;
      dphi=dphi*dt;
      dphimax=max(dphimax,dphi);
      phin[i][k]=phi[i][k]+dphi;
  }

}

为了并行化这个，我们需要修复变量dphimax更新的竞争条件，为此我们可以使用OpenMP的reduction子句：

  #pragma omp for reduction(max: dphimax)
  for (int k=1; k<kmax; k++){
     for (int i=1; i<imax ;i++){
       double dphi=(phi[i+1][k]+phi[i-1][k]-2.*phi[i][k])*dy2i +
                   (phi[i][k+1]+phi[i][k-1]-2.*phi[i][k])*dx2i;
       dphi = dphi * dt;
       dphimax = max(dphimax, dphi);
       phin[i][k] = phi[i][k] + dphi;
     }
 }

我们还可以并行化最后一个循环：

#pragma omp for
for (k=1;k<kmax;k++)
  for (i=1;i<imax;i++)
      phi[i][k]=phin[i][k];

您也可以尝试对其进行矢量化：

#pragma omp simd
for (k=1;k<kmax;k++)
  for (i=1;i<imax;i++)
      phi[i][k]=phin[i][k];

所有东西放在一起：

double dphimax;
#pragma omp parallel
for (it=1;it<=itmax;it++){    
    dphimax=0.0;
    #pragma omp for reduction(max: dphimax)
    for (k=1;k<kmax;k++){
       for (i=1;i<imax;i++){
          dphi=(phi[i+1][k]+phi[i-1][k]-2.*phi[i][k])*dy2i
              +(phi[i][k+1]+phi[i][k-1]-2.*phi[i][k])*dx2i;
          dphi=dphi*dt;
          dphimax=max(dphimax,dphi);
          phin[i][k]=phi[i][k]+dphi;
      }
   }
    #pragma omp for
    for (k=1;k<kmax;k++)
      for (i=1;i<imax;i++)
          phi[i][k]=phin[i][k];
    if(dphimax<eps) break;
}