如何用OpenMP reduction优化矩阵极值（含位置）并行计算？

我想借助OpenMP的reduction子句计算矩阵的总和、最小值与最大值（及其位置）。当前遇到的问题是无法对自定义的Extreme结构体应用min或max归约操作，尝试用declare reduction子句自定义操作也未成功，推测是归约不支持结构体所致。为临时解决问题，我采用critical section来安全更新minimum和maximum结构体的值，但这会导致大矩阵的并行执行速度极慢。请问：该如何优化此方案？是否可以对这类结构体使用归约？或者是否需要换用其他OpenMP子句来处理该问题？

/* 基于OpenMP的矩阵求和
gcc编译要求（版本4.2及以上）：
gcc -O -fopenmp -o matrixSum-openmp matrixSum-openmp.c
运行：./matrixSum-openmp size numWorkers
*/

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

#define MAXSIZE 10000           // 矩阵最大尺寸
#define MATRIX_LIMIT 100        // 矩阵元素范围 [0, MATRIX_LIMIT - 1]
#define MAXWORKERS 4            // 最大工作线程数
#define PROGRAM_EXECUTIONS 5    // 程序执行次数

int numWorkers, size;
int matrix[MAXSIZE][MAXSIZE];
double start_time, end_time;

struct Extreme {
    int value;
    int pos_i, pos_j;
};

struct Extreme min (struct Extreme e1, struct Extreme e2) {
    return e1.value < e2.value ? e1 : e2;
}

struct Extreme max(struct Extreme e1, struct Extreme e2) {
    return e1.value > e2.value ? e1 : e2;
}

void readCommandLine(int argc, char *argv[]) {
    size = (argc > 1) ? atoi(argv[1]) : MAXSIZE;
    numWorkers = (argc > 2) ? atoi(argv[2]) : MAXWORKERS;
    if (size > MAXSIZE)
        size = MAXSIZE;
    if (numWorkers > MAXWORKERS)
        numWorkers = MAXWORKERS;
}

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

double calculateAvg(const double v[PROGRAM_EXECUTIONS]) {
    double avg = 0;
    for(int i = 0; i < PROGRAM_EXECUTIONS; i++)
        avg += v[i];
    return avg / PROGRAM_EXECUTIONS;
}

// 读取命令行参数、初始化矩阵并创建线程
int main(int argc, char *argv[]) {

    // 读取命令行参数（若有）
    readCommandLine(argc, argv);

    // 初始化矩阵
    initializeMatrix();

    // 存储不同处理器数量下的每次执行时间
    // execution_times[i][j] = 使用i个处理器时第j次执行的耗时（秒）
    // 注：execution_times[0][i] 是第i次执行的串行时间
    double execution_times[numWorkers][PROGRAM_EXECUTIONS];

    for(int num_proc = 0; num_proc < numWorkers; num_proc++) {
        // 设置OpenMP线程数
        omp_set_num_threads(num_proc + 1);

        printf("\n============================| 处理器数量: %d |============================\n", num_proc + 1);
        for (int num_prog = 0; num_prog < PROGRAM_EXECUTIONS; num_prog++) {

            // 重置变量
            int total_sum = 0, i, j;
            struct Extreme minimum, maximum;
            minimum.value = MATRIX_LIMIT;
            maximum.value = -MATRIX_LIMIT;

            printf("\n程序执行次数 %d: \n", num_prog + 1);

            // 开始计时
            start_time = omp_get_wtime();

            //#pragma omp declare reduction(myMin : struct Extreme : combinerMin) initializer(omp_orig = initMin)
            //#pragma omp declare reduction(myMax : struct Extreme : combinerMax) initializer(omp_priv = initMax)
            #pragma omp parallel for reduction (+:total_sum) private(j) // (myMin: minimum) (myMax: maximum)
            for (i = 0; i < size; i++) {
                for (j = 0; j < size; j++) {
                    total_sum += matrix[i][j];
                    #pragma omp critical
                    {
                        struct Extreme candidate = {.value = matrix[i][j], .pos_i = i, .pos_j = j};
                        minimum = min(minimum, candidate);
                        maximum = max(maximum, candidate);
                    }
                }
            }

            // 隐式屏障
            end_time = omp_get_wtime();

            execution_times[num_proc][num_prog] = end_time - start_time;

            printf("总和为 %d\n", total_sum);
            printf("最小值为 %d，位置 (%d, %d)\n", minimum.value, minimum.pos_i, minimum.pos_j);
            printf("最大值为 %d，位置 (%d, %d)\n", maximum.value, maximum.pos_i, maximum.pos_j);
            printf("串行执行时间: %f s\n", execution_times[0][num_prog]);
            printf("并行执行时间: %f s\n", execution_times[num_proc][num_prog]);
        }

        double seq_avg = calculateAvg(execution_times[0]);
        double exec_avg = calculateAvg(execution_times[num_proc]);
        double speedup = (seq_avg / exec_avg) * 100;

        printf("\n平均串行时间（%d次执行）: %f", PROGRAM_EXECUTIONS, seq_avg);
        printf("\n使用 %d 个处理器的平均执行时间（%d次执行）: %f s\n", num_proc + 1, PROGRAM_EXECUTIONS, exec_avg);
        printf("加速比: %.2f%%", speedup);
    }

    printf("\n\n");
}

方案一：正确使用自定义结构体的OpenMP归约（推荐） #

OpenMP完全支持对自定义结构体使用reduction，之前失败是因为declare reduction的语法和初始化方式不正确。需要明确定义归约的合并逻辑和私有变量的初始化方式：

1. 添加正确的declare reduction语句
   在并行区域前添加以下代码，分别定义最小值和最大值的归约规则：

// 自定义最小值归约：合并时取value更小的Extreme
#pragma omp declare reduction(myMin : struct Extreme : \
    omp_out = min(omp_out, omp_in)) \
    initializer(omp_priv = (struct Extreme){MATRIX_LIMIT, -1, -1})

// 自定义最大值归约：合并时取value更大的Extreme
#pragma omp declare reduction(myMax : struct Extreme : \
    omp_out = max(omp_out, omp_in)) \
    initializer(omp_priv = (struct Extreme){-MATRIX_LIMIT, -1, -1})

• omp_out是归约后的结果变量，omp_in是当前线程的私有变量
• initializer指定每个线程私有变量的初始值，和全局的初始化逻辑一致

2. 修改并行区域的reduction子句
   去掉critical块，在parallel for中添加自定义的归约：

#pragma omp parallel for reduction (+:total_sum) reduction(myMin: minimum) reduction(myMax: maximum) private(j)
for (i = 0; i < size; i++) {
    for (j = 0; j < size; j++) {
        total_sum += matrix[i][j];
        struct Extreme candidate = {.value = matrix[i][j], .pos_i = i, .pos_j = j};
        minimum = min(minimum, candidate);
        maximum = max(maximum, candidate);
    }
}

这样每个线程会维护自己的minimum和maximum副本，并行结束后自动合并，完全避免critical的性能开销。

方案二：手动维护线程局部最值（兼容旧版本OpenMP） #

如果编译器对自定义结构体归约支持不好，可以手动让每个线程计算自己负责区域的局部最值，最后在主线程合并：

1. 并行区域内定义线程局部变量

#pragma omp parallel private(j)
{
    // 每个线程初始化自己的局部最值
    struct Extreme local_min = {MATRIX_LIMIT, -1, -1};
    struct Extreme local_max = {-MATRIX_LIMIT, -1, -1};
    int local_sum = 0;

    #pragma omp for
    for (i = 0; i < size; i++) {
        for (j = 0; j < size; j++) {
            local_sum += matrix[i][j];
            struct Extreme candidate = {.value = matrix[i][j], .pos_i = i, .pos_j = j};
            local_min = min(local_min, candidate);
            local_max = max(local_max, candidate);
        }
    }

    // 仅在合并全局变量时使用一次critical（大幅减少锁竞争）
    #pragma omp critical
    {
        total_sum += local_sum;
        minimum = min(minimum, local_min);
        maximum = max(maximum, local_max);
    }
}

这种方式只在最后合并时进入一次critical，相比原代码每个元素都加锁，性能提升非常明显。

内容的提问来源于stack exchange，提问作者Pablo