OpenMP向量欧氏距离并行代码正确性及调度子句咨询

任务要求 #

编写OpenMP程序计算两个等长n维向量的欧氏距离，每个线程需生成部分向量并计算元素差平方的局部和，最终所有线程将局部和发送至线程0，由其汇总求和、取平方根并输出结果。

已编写代码 #

1. 串行代码 #

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

int main()
{
    double res=0;
    double d;
    int n=350000;
    int a[n];
    int b[n];
    clock_t start,end;
    
    for (int i=0;i<n;i++){
        a[i]=i+1;
        b[i]=i+3;
    }
    start = clock();
    
    for (int i=0;i<n;i++){
        d = pow((a[i] - b[i]), 2);
        res = res + d;
    }   
    res = sqrt(res);
    end = clock();
    
    printf("The result is: %f\nwork took: %f seconds\n",res,((double)(end-start))/CLOCKS_PER_SEC);

    return 0;
}

2. 初始OpenMP代码 #

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

int main(int argc, char *argv[]) {
    int n = 10000000;
    double start, end;
    double *a = (double *) malloc(n * sizeof(double));
    double *b = (double *) malloc(n * sizeof(double));
    double local_sum = 0, threads_sum = 0;

    // 初始化向量a和b
    for (int i = 0; i < n; i++) {
        a[i] = i;
        b[i] = i * 2;
    }

    start = omp_get_wtime(); 
    omp_set_dynamic(0);
    omp_set_num_threads(16);

    // 计算局部和
    #pragma omp parallel for reduction(+:local_sum)
    for (int i = 0; i < n; i++) {
        double diff = a[i] - b[i];
        local_sum += diff * diff;
    }

    // 计算线程总和
    #pragma omp parallel for reduction(+:threads_sum)
    for (int i = 0; i < omp_get_num_threads(); i++) {
        threads_sum += local_sum;
    }

    // 计算欧氏距离
    double euclidean_distance = sqrt(threads_sum);
    end = omp_get_wtime();

    printf("Euclidean distance between vectors a and b: %f\nTime taken: %f\n", euclidean_distance, end - start);

    free(a);
    free(b);

    return 0;
}

反馈意见 #

你的串行代码没问题，但OpenMP代码不符合要求。只需在串行代码的循环上方添加OpenMP语句即可实现并行执行。

简化版OpenMP代码 #

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

int main() {
    int n = 350000;
    double start, end;
    double res = 0;
    double d;
    int *a = (int *)malloc(n * sizeof(int));
    int *b = (int *)malloc(n * sizeof(int));

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

    start = omp_get_wtime();

    #pragma omp parallel for reduction(+:res)
    for (int i = 0; i < n; i++) {
        d = pow((a[i] - b[i]), 2);
        res += d;
    }

    res = sqrt(res);
    end = omp_get_wtime();

    printf("The result is: %f\nWork took: %f seconds\n", res, end - start);

    free(a);
    free(b);

    return 0;
}

问题 #

请问该简化版代码是否正确？或者是否需要为循环添加schedule(static)子句（如#pragma omp parallel for schedule(static) reduction(+:res)）？

代码正确性判断 #

你的简化版OpenMP代码是完全正确的，完美契合任务要求和反馈意见：

1. 它在串行代码的核心求和循环上方添加了#pragma omp parallel for reduction(+:res)：
   • parallel for会自动将循环迭代均匀分配给多个线程执行
   • reduction(+:res)机制会让每个线程维护自己的局部res副本，计算完成后自动将所有局部和汇总到主线程的res中，正好实现了“每个线程计算局部和，最终汇总到线程0”的需求
2. 代码逻辑和串行版本完全一致，计算结果准确，同时用omp_get_wtime()替代clock()，更适合测量并行程序的运行时间。

是否需要添加schedule(static)？ #

不需要额外添加schedule(static)。因为OpenMP中parallel for的默认调度方式就是static，它会将循环迭代平均分配给各个线程，对于这种每个迭代计算量完全相同的均匀循环，默认的static调度已经是最优选择，手动添加只是显式写出了默认行为，不会带来性能提升。

另外可以做个小优化：把pow((a[i]-b[i]),2)改成(a[i]-b[i])*(a[i]-b[i])，pow函数用于整数平方的开销比直接乘法大，能进一步提升程序性能。

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