OpenMP并行编程

OpenMP并行编程入门

OpenMP是一个编译器指令和库函数的集合，主要是为共享式存储计算机上的并行程序设计使用的。目前支持OpenMP的语言主要有Fortran，C/C++。

简单入门

#include <stdio.h>
#include <omp.h>
​
int main(void)
{
    omp_set_num_threads(10);
#pragma omp parallel
    {
        printf("Hello World\n");
    }
    return 0;
}

编译 gcc 1.c -o 1 -fopenmp

使用线程求解积分问题

我们将编写一个程序，该程序通过将曲线下的面积近似为矩形面积的和来估计一个定积分的值。选择积分和积分范围，使得这个积分的结果等于圆周率。

简单的原始代码：

/**
 * 求解圆周率
 */
#include <stdio.h>
#include <omp.h>
​
static long num_steps = 1000000000;
double step;
int main(void)
{
    step = 1.0 / (double)num_steps;
    // 使用 OpenMP 的计时器
    double start_time = omp_get_wtime();
    double sum = 0.0;
    for (int i = 0; i < num_steps; i++)
    {
        double x = (i + 0.5) * step;
        sum += 4.0 / (1.0 + x * x);
    }
    double pi = step * sum;
    double end_time = omp_get_wtime();
    double run_time = end_time - start_time;
    printf("pi = %lf, run_time = %lf seconds\n", pi, run_time);
    return 0;
}

多线程版本：

/**
 * 多线程求解圆周率
 */
#include <stdio.h>
#include <omp.h>
​
#define NUM_THREADS 4
​
static long num_steps = 1000000000;
double step;
int main(void)
{
    step = 1.0 / (double)num_steps;
    double results[NUM_THREADS] = {0.0};
    // 设置线程数量
    omp_set_num_threads(NUM_THREADS);
    double start_time = omp_get_wtime(); // 获取开始时间
    
    
    
#pragma omp parallel
    {
        // 当前线程的id
        int thread_id = omp_get_thread_num();
        for (int i = thread_id; i < num_steps; i += NUM_THREADS)
        {
            double x = (i + 0.5) * step;
            results[thread_id] += 4.0 / (1.0 + x * x);
        }
    }
    // 收集计算结果
    double sum = 0.0;
    for (int i = 0; i < NUM_THREADS; i++)
    {
        sum += results[i];
    }
    double pi = sum * step;
    double end_time = omp_get_wtime();       // 获取结束时间
    double run_time = end_time - start_time; // 计算运行时间
    printf("pi = %lf, run_time = %lf seconds\n", pi, run_time);
    return 0;
}
​

伪内存共享

上面的代码，可能会出现多线程比单线程还要慢的情况，这是因为什么呢？

看下图：

多个物理核心是共享一级缓存的，而对缓存操作的单位是缓存行，根据缓存一致性协议，同一时刻只能有一个线程去操纵缓存行，将会导致性能的下降。为什么同时只能有一个线程去操控缓存行呢？

如果线程A去修改Core1缓存中的x变量，由于缓存一致性协议的存在，Core2中对应的缓存行缓存的x，y变量将会失效，它会强制从主内存中去重新加载变量，这样会频繁的访问主内存，缓存基本都失效了，将会导致性能的下降。

如何解决

伪共享的原因就是变量在同一个缓存行上面，解决方案也很简答，就不要让变量在同一个缓存行上面即可。

字节填充：

代码如下：

/**
 * 字节填充解决伪共享
 */
#include <stdio.h>
#include <omp.h>
​
#define NUM_THREADS 4
​
static long num_steps = 1000000000;
double step;
​
int main(void)
{
    step = 1.0 / (double)num_steps;
    // 现代计算机缓存行一般是32-256字节
    double results[NUM_THREADS][16] = {{0.0}};
    // 设置线程数量
    omp_set_num_threads(NUM_THREADS);
    double start_time = omp_get_wtime(); // 获取开始时间
#pragma omp parallel
    {
        int thread_id = omp_get_thread_num();
        for (int i = thread_id; i < num_steps; i += NUM_THREADS)
        {
            double x = (i + 0.5) * step;
            results[thread_id][0] += 4.0 / (1.0 + x * x);
        }
    }
    // 收集计算结果
    double sum = 0.0;
    for (int i = 0; i < NUM_THREADS; i++)
    {
        sum += results[i][0] * step;
    }
    double pi = sum * step;
    double end_time = omp_get_wtime();
    double run_time = end_time - start_time;
    printf("pi = %lf, run_time = %lf seconds\n", pi, run_time);
    return 0;
}

临界区

最基本的同步构造为多线程运行的代码定义了一种相互排斥关系。相互排斥规定，如果一个线程正在执行一个代码块，而第二个线程试图执行同样的代码块，那么第二个线程将暂停并等待，直到第一个线程执行完成。在OpenMP中，用critical构造来定义相互排斥执行的代码块。

带有临界区的数值积分。每个线程分配的私有变量用来存储部分和，这些私有变量极不可能存储在相同的L1缓存行上，因此，不会出现伪共享。部分和在临界区内部合并，因此不会出现数据竞争:

/**
 * 临界区
 */
#include <stdio.h>
#include <omp.h>
​
#define NUM_THREADS 4
​
static long num_steps = 1000000000;
double step;
int main(void)
{
    step = 1.0 / (double)num_steps;
    double results[NUM_THREADS] = {0.0};
    // 设置线程数量
    omp_set_num_threads(NUM_THREADS);
    double sum = 0.0;
    double start_time = omp_get_wtime(); // 获取开始时间
#pragma omp parallel
    {
        // 当前线程的id
        int thread_id = omp_get_thread_num();
        double local_sum = 0.0;
        for (int i = thread_id; i < num_steps; i += NUM_THREADS)
        {
            double x = (i + 0.5) * step;
            local_sum += 4.0 / (1.0 + x * x);
        }
#pragma omp critical // 临界区
        sum += local_sum;
    }
    double pi = sum * step;
    double end_time = omp_get_wtime();       // 获取结束时间
    double run_time = end_time - start_time; // 计算运行时间
    printf("pi = %lf, run_time = %lf seconds\n", pi, run_time);
    return 0;
}

栅栏

OpenMP中最常用的同步结构是栅栏（barrier）。栅栏定义了程序中的一个点，在这个点上，所有线程必须在任意线程通过栅栏之前到达。我们在讨论并行构造结束处的行为时已经遇到了栅栏。当一个线程组到达并行构造内的结构化代码块的末端时，每个线程都会在结构化代码块的末端等待，直到所有线程都到达。然后线程组中的主线程，也就是最初遇到并行构造的线程继续执行，而其他线程则关闭。使用同步构造的术语，我们说一个并行区域的结束隐含着一个栅栏。

barrier指令是一个独立的、可执行的指令。

一个简单的例子

/**
 * 栅栏
 */
#include <stdio.h>
#include <omp.h>
#define NUM_THREADS 4
int main(void)
{
    omp_set_num_threads(NUM_THREADS);
#pragma omp parallel num_threads(4)
    {
        printf("Hello world! I am thread %d\n", omp_get_thread_num());
#pragma omp barrier
        printf("Thread %d is out of the barrier\n", omp_get_thread_num());
    }
    return 0;
}

并行化循环

考虑一个将两个向量a和b相加的简单程序。

for (int i = 0; i < N; i++) {
    a[i] = a[i] + b[i];
}

我们可以使用一个简单的指令omp for实现循环并行化：

#pragma omp parallel {
    #pragma omp for
    for (int i = 0; i < N; i++) {
        a[i] = a[i] + b[i];
    }
}

上面的指令也可以写成如下格式：

#pragma omp parallel for
for (int i = 0; i < N; i++) {
	a[i] = a[i] + b[i];
}

对于所有的循环并行，在循环代码块的末尾都有一个隐式栅栏。所有的线程都会在共享工作循环构造的结尾处等待，直到所有在该构造上运行的线程组结束。

规约

OpenMP中的规约（Reduction）是一种用于在并行区域中合并数据的操作，它允许多个线程协同工作，将多个值合并为一个。规约操作通常用于累加和（Summation）、乘积（Product）等操作。

比如：

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

int main() {
    int arr[100];
    int sum = 0;

    // 初始化数组
    for (int i = 0; i < 100; i++) {
        arr[i] = i;
    }

    // 使用OpenMP的并行for循环和reduction子句进行求和
    #pragma omp parallel for reduction(+:sum)
    for (int i = 0; i < 100; i++) {
        sum += arr[i];
    }

    printf("Sum is: %d\n", sum);
    return 0;
}

在这个例子中，每个线程都会计算自己的局部和，然后OpenMP运行时会在所有线程完成工作后，将这些局部和合并为最终的总和。这样可以避免在并行区域中使用锁或临界区，从而提高程序的并行效率。

需要注意的是，reduction子句中的变量在并行区域中被视为私有的，每个线程都有自己的副本。在并行区域开始时，这些私有变量会被初始化为指定的运算符所要求的初始值（例如，对于加法和减法是0，对于乘法是1）。在并行区域结束时，所有私有变量的值会根据reduction子句中指定的操作符进行合并，并将最终结果存储在主线程的原始变量中。

使用循环并行求解积分

通过循环并行化，我们可以使用一种简单而优雅的方法来并行化循环：

/**
 * 循环并行求解积分
 */
#include <stdio.h>
#include <omp.h>

#define NUM_THREADS 4

static long num_steps = 1000000000;
double step;
int main(void)
{
    step = 1.0 / (double)num_steps;
    double sum = 0.0;
    omp_set_num_threads(NUM_THREADS);
    double start_time = omp_get_wtime();
#pragma omp parallel for reduction(+ : sum)
    for (int i = 0; i < num_steps; i++)
    {
        double x = (i + 0.5) * step;
        sum += 4.0 / (1.0 + x * x);
    }
    double pi = step * sum;
    double end_time = omp_get_wtime();
    double run_time = end_time - start_time;
    printf("pi = %lf, run_time = %lf seconds\n", pi, run_time);
    return 0;
}

矩阵乘法

串行算法：

void matrixMultiply(int n, int m, int p, int A[n][m], int B[m][p], int C[n][p]) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < p; j++) {
            C[i][j] = 0;
            for (int k = 0; k < m; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

使用openmp：

注意事项：

• 在使用这种并行化方法时，需要确保循环的迭代之间没有依赖关系，这样才能安全地并行执行而不会产生竞争条件或数据不一致的问题。

• 由于 A, B, 和 C 被声明为共享变量，需要确保线程安全，特别是当多个线程尝试同时写入同一个内存位置时。在这个特定的例子中，每个线程写入的是 C 矩阵的不同部分，因此不会产生冲突。

/**
 * 矩阵乘法
 */
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#define N 4 // 定义矩阵大小

// 初始化矩阵
void initializeMatrix(int n, int matrix[n][n]) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            matrix[i][j] = rand() % 10; // 随机生成0-9之间的数字
        }
    }
}


// 并行分块矩阵乘法
void MatrixMultiply(int n, int A[n][n], int B[n][n], int C[n][n]) {
    #pragma omp parallel for collapse(2) shared(A, B, C)
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            C[i][j] = 0;
            for (int k = 0; k < n; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

void printMatrix(int n, int matrix[n][n]) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%d ", matrix[i][j]);
        }
        printf("\n");
    }
}

int main() {
    int A[N][N], B[N][N], C[N][N];

    // 初始化矩阵A和B
    initializeMatrix(N, A);
    initializeMatrix(N, B);

    // 执行并行分块矩阵乘法
    MatrixMultiply(N, A, B, C);

    // 打印结果
    printf("Matrix A:\n");
    printMatrix(N, A);
    printf("Matrix B:\n");
    printMatrix(N, B);
    printf("Matrix C (Product):\n");
    printMatrix(N, C);

    return 0;
}

#pragma omp parallel for

• #pragma omp 是告诉编译器下一行代码是一个OpenMP指令。

• parallel for 结合使用表示这是一个并行区域，其中的for循环将被多个线程并行执行。

collapse(2)

• collapse(2) 是一个扩展指令，用来并行化多层嵌套的循环。

• 数字 2 表示应该将接下来的两层循环合并为一个大的迭代空间，并在不同线程之间分配这些迭代。

• 这意味着在给定的例子中，i 循环和 j 循环会被合并，然后这些合并后的迭代会被分配给不同的线程执行。

shared(A, B, C)

• shared 说明在并行区域内，括号内的变量（在这里是 A, B, 和 C）是共享的。

• 共享变量意味着所有线程都使用相同的内存位置访问这些变量，而不是每个线程有自己的副本。

分块加速并行

上述过程没有使用“分块”（block-wise）矩阵乘法。它只是使用了OpenMP来并行化标准的三重嵌套循环矩阵乘法。在这个例子里，每个元素的计算都是独立的，但没有对矩阵进行物理上的分块处理。分块（包括行列分块、行行分块、列列分块）是并行化的基本思路。

分块矩阵乘法是一种算法，它将大矩阵划分成较小的子块（blocks），然后对这些子块执行乘法。这种方法可以提高缓存效率，因为小块更有可能完全适配在CPU缓存中，从而减少内存访问次数。

使用openmp实现矩阵分块乘法：

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

#define N 2<<8          // 矩阵大小
#define BLOCK_SIZE 32 // 块大小

// 初始化矩阵
void initializeMatrix(int n, int matrix[N][N]) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            matrix[i][j] = rand() % 10;
        }
    }
}

// 并行分块矩阵乘法
void blockMatrixMultiply(int n, int A[N][N], int B[N][N], int C[N][N]) {
    #pragma omp parallel for shared(A, B, C) collapse(2)
    for (int i = 0; i < n; i += BLOCK_SIZE) {
        for (int j = 0; j < n; j += BLOCK_SIZE) {
            int temp[BLOCK_SIZE][BLOCK_SIZE] = {0};
            for (int k = 0; k < n; k += BLOCK_SIZE) {
                for (int ii = i; ii < i + BLOCK_SIZE && ii < n; ii++) {
                    for (int jj = j; jj < j + BLOCK_SIZE && jj < n; jj++) {
                        for (int kk = k; kk < k + BLOCK_SIZE && kk < n; kk++) {
                            temp[ii - i][jj - j] += A[ii][kk] * B[kk][jj];
                        }
                    }
                }
            }
            for (int ii = i; ii < i + BLOCK_SIZE && ii < n; ii++) {
                for (int jj = j; jj < j + BLOCK_SIZE && jj < n; jj++) {
                    C[ii][jj] = temp[ii - i][jj - j];
                }
            }
        }
    }
}

int main() {
    int A[N][N], B[N][N], C[N][N] = {0};

    // 初始化矩阵A和B
    initializeMatrix(N, A);
    initializeMatrix(N, B);

    // 执行并行分块矩阵乘法
    blockMatrixMultiply(N, A, B, C);

    // 打印结果
    printf("Matrix A:\n");
    printMatrix(N, A);
    printf("Matrix B:\n");
    printMatrix(N, B);
    printf("Matrix C (Product):\n");
    printMatrix(N, C);

    return 0;
}