使用 pthreads 和 OpenMp 在 C 中并行化 Simpson 方法

Question

我对整个多处理和并行化世界还很陌生。我目前正在研究一种实现辛普森方法和辛普森 3/8 的算法。我已经以串行形式在 C 中实现了该算法。

我需要一些帮助来使用 OpenMp 和 Pthreads 并行化这两种算法。欢迎和赞赏任何建议。

//Simpson.c
#include <stdio.h>
#include <stdlib.h>

double function(double x) //función a la cual se le aplicará el método
{
  return 0.5*x;    
}

int main(int argc, char*argv[])
{
  double resultado = 0.0; //resultado de la integral tras aplicar el método de simpson
  double a = 0.0;  //límite inferior de la integral
  double b = 0.0;  //límite superior de la integral
  int n = 0; //cantidad de particiones entre los intervalos, tiene que ser un número par
  int i = 0; //variable para poder iterar sobre las diferentes particiones
  double h = 0.0; //variable para guardar el incremento que habrá entre cada valor de "x"
  double calc = 0.0; //variable intermedia para ir acumulando el resultado de evaluar las diferentes "x" en la función
  
  //variables para almacenar los diferentes valores de "x", es decir los límites inferiores y superiores y el valor intermedio de cada partición  
  double x0 = 0.0;
  double x1 = 0.0;
  double x2 = 0.0;  
    
  printf("Introduce el límite inferior: ");
  scanf("%lf", &a);
  printf("Introduce el límite superior: ");
  scanf("%lf", &b);
  printf("Introduce la cantidad de particiones (número par): ");
  scanf("%d", &n);  
    
  h = (b-a)/n; //se calcula el incremento para poder iterar sobre todas las particiones
    
      //Para el cálculo de la integral se utilizan xj, x(j+1) y x(j+2) en el subintervalo [xj, x(j+2)] para cada i = 1, 3, 5, ..., n
    
  for (i=0; i<=n; i+=2)
  {   
    x0 = a + h*i; //límite inferior
    calc += function(x0); //se evalua la x en la función
    
    x1 = a + h*(i+1); //valor intermedio
    calc += 4*function(x1);
      
    x2 = a + h*(i+2); //límite superior
    calc += function(x2);
      
    calc *= ((x2-x0)/6) ; //se vuelve a calcular el incremento entre nuestros límites, y se divide entre 6
      
    resultado += calc; //variable para ir acumulando los cálculos intermedios de cada "x"
  }
    
  printf("La aproximación a la integral es: %f \n", resultado);

}

#include<stdio.h>
#include<math.h>

// función a integrar
#define f(x) (0.5*x)

int main()
{
  float b; //Limite superior
  float a; //Limite inferior
  float resultado=0.0; 
  float dx;
  float k;
  int i; 
  int N;  //numero de intervalos o iteraciones

  
 // Entradas
 printf("Dame el limite inferior: ");
 scanf("%f", &a);
 printf("Dame el limite superior: ");
 scanf("%f", &b);
 printf("Número de iteraciones/intervalos (número par): ");
 scanf("%d", &N);

  //Calculando delta de x
  dx = (b - a)/N;
  
  //Sumas de la integracion
  resultado = f(a) + f(b); //primera suma
  
  for(i=1; i<= N-1; i++)
  {
    k = a + i*dx; //Suma del intervalo
    
    if(i%3 == 0)
    {
      resultado = resultado + 2 * f(k);
    }
    else
    {
      resultado = resultado + 3 * f(k);
    }
  }
  
  resultado = resultado * dx*3/8;
  
  printf("\nEl valor de la integral es: %f \n", resultado);
  printf("\n");
  
  return 0;
}

Answer 1

我的第一个建议是首先阅读有关共享内存范式及其问题的大量信息。然后阅读 Pthreads 和 OpenMP，然后再看看具体的 OpenMP 和 Pthread 并行化，以及它们如何处理共享内存范式的一些问题。如果你必须同时学习 Pthreads 和 OpenMP，我认为，如果我是你，我会先研究 Pthread，然后再研究 OpenMP。后者有助于解决使用前者编码的一些繁琐细节和陷阱。

我将尝试为您提供一个非常笼统（不是超级详细）公式关于如何处理代码并行化的方法，我将用作示例您的第二个代码和 OpenMP。

如果您选择并行化算法，首先您需要查看值得并行化的代码部分（即，并行化的开销被它带来的加速所掩盖） .因为这种分析是基本的，但是，对于像你这样的较小的代码，它很容易弄清楚，它通常是循环，例如：

for(i=1; i<= N-1; i++)
  {
    k = a + i*dx; //Suma del intervalo
    
    if(i%3 == 0)
    {
      resultado = resultado + 2 * f(k);
    }
    else
    {
      resultado = resultado + 3 * f(k);
    }
  }

要并行化这段代码，我应该使用哪个构造函数 tasks ？ parallel for ?您可以在 SO thread 上查看该主题的精彩答案。在您的代码中，很明显我们应该使用parallel for（即#pragma omp parallel for）。因此，我们告诉 OpenMP 在线程之间划分该循环的迭代。

现在您应该考虑在并行区域内使用的变量应该是线程间的private 还是shared。为此，openMP 提供了 private、firstprivate、lastprivate 和 shared 等构造函数。关于这个主题的另一个很棒的 SO 线程可以在 here 找到。要对此进行评估，您确实需要了解 shared memory 范式以及 OpenMP 如何处理它。

然后看看潜在的竞争条件、变量之间的相互依赖关系等等。 OpenMP 提供了 critical、atomic operations、reductions 等构造函数来解决这些问题（其他出色的 SO 线程 atomic vs critical 和 reduction vs atomic）。 TL;DR：您应该选择在不影响正确性的情况下提供最佳性能的解决方案。

在您的情况下，变量“resultado”的reduction 显然是最佳选择：

 #pragma omp parallel reduction ( +: resultado) 
 for(i=1; i<= N-1; i++)
  {
    k = a + i*dx; //Suma del intervalo
    
    ... 
  }

保证代码的正确性后，您可以查看负载不平衡问题（即每个线程执行的工作之间的差异）。为了解决这个问题，openMP 提供了并行 for 调度策略，例如 dynamic 和 guided（另一个 SO thread）以及调整这些分布的块大小的可能性。

在您的代码中，线程将执行相同数量的并行工作，因此无需使用dynamic 或guided 调度，您可以选择static 替代。 static 相对于其他的好处是它不会引入必须在运行时协调线程之间的任务分配的开销。

要从运行代码的架构中提取最大性能，还有很多东西需要研究。但就目前而言，我认为这是一个足够好的开始。