使用 OpenMP 并行化一维矩阵乘法

Question

我是 OpenMP 的新手，我正在尝试并行化一维矩阵乘法（仅乘以矩阵的上三角）。我现在的解决方案很快，但它目前只使用了大约一半的可用线程，所以我应该进一步并行化它但我不知道从哪里开始（我知道在结果数组中保存数据是最耗时的，这个数组是外部的，只能将数据存入其中）。

这是代码（我同时进行乘法和转置）：

#pragma omp parallel
{
#pragma omp for
for (int j=0; j < ny; ++j){
    for (int i=0; i < ny; ++i){
        if (i < j){
            continue;
        }
        double mat_mul_sum{0.};
        for (int k=0; k < nx; ++k){
             double mat_mul = new_data[k + j*nx] * new_data[k + i*nx];
             mat_mul_sum += mat_mul;
        }
        array_results[i + j*ny] = mat_mul_sum;
    }
}
}

非常感谢帮助。

Answer 1

首先，您可以从删除循环中的 continue 开始，然后 在值 j[=28 处开始基于 i 的循环=].然后，您可以专注于线程之间的负载平衡问题。

平衡线程之间的工作的一个简单解决方案是同时处理矩阵的开头和结尾（在每个线程中）。从视觉上看，就像是把三角形切掉一个角，然后移动，形成一个ny x (ny/2)大小的矩形。请注意 ny 为奇数的情况应仔细考虑（因为形状不会是纯矩形）。

最后，您可以使用 SIMD 减少 来优化 k 循环。

结果代码应如下所示（未经测试）：

#pragma omp parallel
{
    #pragma omp for
    for (int j=0; j < (ny+1)/2; ++j)
    {
        const int jUp = j;
        const int jDown = ny - 1 - j;

        for (int i=jUp; i < ny; ++i)
        {
            double mat_mul_sum{0.};
            #pragma omp simd reduction(+:mat_mul_sum)
            for (int k=0; k < nx; ++k){
                 double mat_mul = new_data[k + jUp*nx] * new_data[k + i*nx];
                 mat_mul_sum += mat_mul;
            }
            array_results[i + jUp*ny] = mat_mul_sum;
        }

        // Compute the center only once
        if(jUp != jDown)
        {
            for (int i=jDown; i < ny; ++i)
            {
                double mat_mul_sum{0.};
                #pragma omp simd reduction(+:mat_mul_sum)
                for (int k=0; k < nx; ++k){
                     double mat_mul = new_data[k + jDown*nx] * new_data[k + i*nx];
                     mat_mul_sum += mat_mul;
                }
                array_results[i + jDown*ny] = mat_mul_sum;
            }
        }
    }
}