如何将 OpenMp 添加到三重嵌套 for 循环答案

【问题标题】：How to add OpenMp to triple nested for-loop如何将 OpenMp 添加到三重嵌套 for 循环
【发布时间】：2021-06-21 20:09:08
【问题描述】：

目标是将尽可能多的 OpenMP 添加到以下 Cholesky 因子函数中以提高并行度。到目前为止，我只有一个 #pragma omp parallel for 正确实施。 vector<vector<double>> 代表一个二维矩阵。我已经尝试为
for (int i = 0; i < n; ++i)、for (int k = 0; k < i; ++k) 和for (int j = 0; j < k; ++j) 添加#pragma omp parallel for，但并行化出错了。 makeMatrix(n, n) 初始化大小为 nxn 的全零的 vector<vector<double>>。

vector<vector<double>> cholesky_factor(vector<vector<double>> input)
{
    int n = input.size();
    vector<vector<double>> result = makeMatrix(n, n);
        
    for (int i = 0; i < n; ++i) 
    {
        for (int k = 0; k < i; ++k)
        {
            double value = input[i][k];
            for (int j = 0; j < k; ++j)
            {
                value -= result[i][j] * result[k][j];
            }
            result[i][k] = value / result[k][k];
        }
        double value = input[i][i];
        #pragma omp parallel for
        for (int j = 0; j < i; ++j)
        {
            value -= result[i][j] * result[i][j];
        }
        result[i][i] = std::sqrt(value);
    }

    return result;
}

【问题讨论】：

“但并行化出错了” - 你需要更具体。
“向嵌套循环添加更多 OpenMP”只会在外部循环本身不能提供足够的并行度时提供更好的性能。在这种情况下，可以使用collapse 子句，但在这里我不认为有一种简单的方法可以做到这一点，因为(i,j,k) 的不同组合在读写时矩阵的区域重叠。要使第二个 j-loop 的当前并行化正确，您必须添加 reduction(-: value)。
您的内循环边界也取决于外循环索引。因此你不能使用collapse

标签： c++ multithreading performance parallel-processing openmp

【解决方案1】：

我不认为你可以用这个算法进行更多的并行化，因为外部循环的ith 迭代取决于i - 1th 迭代的结果和kth 迭代的结果内循环取决于k - 1th 迭代的结果。

vector<vector<double>> cholesky_factor(vector<vector<double>> input)
{
    int n = input.size();
    vector<vector<double>> result = makeMatrix(n, n);
        
    for (int i = 0; i < n; ++i) 
    {
        for (int k = 0; k < i; ++k)
        {
            double value = input[i][k];
            // reduction(-: value) does the same 
            // (private instances of value are initialized to zero and
            // added to the initial instance of value when the threads are joining
            #pragma omp parallel for reduction(+: value)
            for (int j = 0; j < k; ++j)
            {
                value -= result[i][j] * result[k][j];
            }
            result[i][k] = value / result[k][k];
        }
        double value = input[i][i];
        #pragma omp parallel for reduction(+: value)
        for (int j = 0; j < i; ++j)
        {
            value -= result[i][j] * result[i][j];
        }
        result[i][i] = std::sqrt(value);
    }

    return result;
}

【讨论】：