首页 > 解决方案 > 如何将 OpenMp 添加到三重嵌套 for 循环

问题描述

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

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;
}

标签: c++multithreadingperformanceparallel-processingopenmp

解决方案


我认为你不能用这个算法进行更多的并行化,因为i外部循环的i - 1第 th 次k迭代取决于第 th 次迭代的结果,而内部循环的第 th 次迭代取决于第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;
}

推荐阅读