LU分解多线程实现异常排查：C++标准线程版本结果不稳定

已有可正常运行的单线程LU分解实现，使用C++标准线程改造为多线程版本后，运行结果不稳定，有时正确有时错误。多线程版本代码如下：

#define MAX_THREADS std::thread::hardware_concurrency()

void LUDecompositionMT(std::vector<std::vector<double>> matrix){
    int n = matrix.size();

    std::vector<std::vector<double>> lower(n, std::vector<double>(n));
    std::vector<std::vector<double>> upper(n, std::vector<double>(n));
    std::vector<std::thread> threads;

    upper = matrix;

    for (int c = 0; c < MAX_THREADS; c++) {
        threads.emplace_back(std::thread([=, &lower, &upper] {
            int start = (n * c / MAX_THREADS);
            int end = (n * (c + 1) / MAX_THREADS);

            for (int k = start == 0 ? 1 : start; k < end; k++)
            {
                for (int i = k - 1; i < n; i++)
                    for (int j = i; j < n; j++)
                    {
                        lower[j][i] = upper[j][i] / upper[i][i];
                    }


                for (int i = k; i < n; i++)
                    for (int j = k - 1; j < n; j++)
                    {
                        upper[i][j] = upper[i][j] - lower[i][k - 1] * upper[k - 1][j];

                    }
            }

            }));
    }

    for (auto& t : threads) {
        t.join();
    }
}

错误原因分析 #

• 数据竞争问题：多个线程同时对lower和upper的同一元素进行读写操作，无任何同步机制。不同线程处理的k区间可能操作相同的行/列，导致并发读写时数据被覆盖，出现不可预测的结果。
• LU分解的顺序依赖被破坏：LU分解的每一步消元（k步）依赖于前一步（k-1步）的计算结果。当前代码将k的区间直接分给不同线程并行执行，线程可能在k-1步未完成时就开始执行k步，使用未正确更新的lower和upper数据，导致计算错误。
• 任务划分逻辑错误：将k循环的区间拆分给线程是不合理的，因为k步本身必须串行执行，并行化的正确位置应该是在每个k步内部，处理相互独立的行或列。
• 遗漏lower矩阵初始化：LU分解中lower矩阵的对角线元素应为1，原代码未做初始化，这也是结果异常的潜在原因。

修复后的多线程实现 #

我们需要保证k步的串行执行，在每个k步内部，将独立的行处理任务分配给不同线程，避免数据竞争：

#include <vector>
#include <thread>

#define MAX_THREADS std::thread::hardware_concurrency()

void LUDecompositionMT(std::vector<std::vector<double>> matrix) {
    int n = matrix.size();
    std::vector<std::vector<double>> lower(n, std::vector<double>(n, 0.0));
    std::vector<std::vector<double>> upper = matrix;
    std::vector<std::thread> threads;

    // 初始化lower矩阵的对角线为1
    for (int i = 0; i < n; ++i) {
        lower[i][i] = 1.0;
    }

    // 串行处理每一个消元步k
    for (int k = 0; k < n - 1; ++k) {
        // 并行计算当前k列的lower元素
        auto computeLower = [&](int start, int end) {
            for (int i = start; i < end; ++i) {
                lower[i][k] = upper[i][k] / upper[k][k];
            }
        };

        // 拆分lower计算任务到多个线程
        threads.clear();
        int chunk = (n - k - 1) / MAX_THREADS;
        for (int c = 0; c < MAX_THREADS; ++c) {
            int start = k + 1 + c * chunk;
            int end = (c == MAX_THREADS - 1) ? n : start + chunk;
            if (start < n) {
                threads.emplace_back(computeLower, start, end);
            }
        }
        for (auto& t : threads) {
            t.join();
        }

        // 并行计算upper矩阵的后续行
        auto computeUpper = [&](int start, int end) {
            for (int i = start; i < end; ++i) {
                for (int j = k + 1; j < n; ++j) {
                    upper[i][j] -= lower[i][k] * upper[k][j];
                }
            }
        };

        // 拆分upper计算任务到多个线程
        threads.clear();
        chunk = (n - k - 1) / MAX_THREADS;
        for (int c = 0; c < MAX_THREADS; ++c) {
            int start = k + 1 + c * chunk;
            int end = (c == MAX_THREADS - 1) ? n : start + chunk;
            if (start < n) {
                threads.emplace_back(computeUpper, start, end);
            }
        }
        for (auto& t : threads) {
            t.join();
        }
    }

    // 可添加lower/upper矩阵的输出或返回逻辑
}

修复说明 #

1. 保证k步串行执行：外层k循环串行，确保每一步消元都基于前一步的正确结果，不会出现依赖错误。
2. 并行处理独立任务：在每个k步内部，将lower列的计算和upper行的计算拆分为多个线程任务，这些任务之间无数据依赖，避免竞争。
3. 补全lower矩阵初始化：设置对角线元素为1，符合LU分解的定义。
4. 合理拆分任务：按行拆分任务，每个线程处理连续的行区间，提升缓存效率的同时避免线程间资源冲突。

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