大型稀疏Armadillo矩阵多线程构建慢于单线程的问题及优化

我正在编写代码，需要以位置独立的方式构建超大型矩阵——矩阵任意位置的值不依赖其他位置的值，所以我觉得这个场景非常适合并行化处理。具体逻辑如下：

1. 根据线程数划分矩阵位置（确保无竞争），为每个线程分配内存存储局部中间矩阵（尺寸3x3至8x8）来计算所需值；
2. 执行计算确定要输入的值，再把这些值赋值到稀疏Armadillo矩阵中。

但实际运行时多线程速度反而比单线程慢，我怀疑和稀疏矩阵格式有关。请问这是什么原因导致的？有哪些方法可以提升矩阵构建的速度？

以下是我用的多线程代码：

void PopulateGlobalStiffnessMatrixMultiThreaded(arma::sp_mat* global_stiffness_matrix, Mesh* mesh, Material* material,
                                                int number_of_cpus){
    int step = mesh->elements.n_rows / number_of_cpus;
    int start_element = 0;

    std::vector<std::thread> threads;

    for (int i = 0; i < number_of_cpus; i++){
        threads.push_back(std::thread(PopulateGlobalStiffnessMatrix, global_stiffness_matrix, mesh, material, start_element, step));
        start_element = start_element + step;
    }

    threads.push_back(std::thread(PopulateGlobalStiffnessMatrix, global_stiffness_matrix, mesh, material, start_element,
                                  mesh->elements.n_rows - start_element - 1));

    for(std::thread &t : threads){
        if (t.joinable()){
            t.join();
        }
    }

}

为什么多线程反而更慢？ #

• 稀疏矩阵的全局锁竞争：Armadillo的sp_mat并非线程安全的写入结构，底层存在隐式的全局锁或内存同步逻辑。多个线程同时往同一个全局稀疏矩阵写入数据时，会频繁触发锁等待，并行计算的收益被完全抵消，甚至因为锁开销拖慢整体速度。
• 线程调度与创建开销：如果待处理的元素总数不算极大，创建多线程的初始化、上下文切换开销可能超过并行计算节省的时间。另外你代码里额外多创建了一个线程，导致线程总数为number_of_cpus+1，容易引发过度调度，进一步增加开销。
• 任务划分错误：最后一个线程的任务量计算有误（应该是mesh->elements.n_rows - start_element，而非mesh->elements.n_rows - start_element -1），会导致少处理一个元素，同时任务划分不均可能让某个线程负载过高，其他线程处于等待状态，降低整体效率。

提升矩阵构建速度的方法 #

1. 局部构建+全局合并（最核心优化） #

让每个线程独立构建自己的局部稀疏矩阵，完全避免全局写入竞争，最后再合并所有局部矩阵到全局：

// 修改线程函数为构建并返回局部稀疏矩阵
arma::sp_mat PopulateLocalStiffnessMatrix(Mesh* mesh, Material* material, int start_element, int count) {
    arma::sp_mat local_mat;
    // 原有计算逻辑，仅填充当前线程负责的局部矩阵
    return local_mat;
}

// 多线程主控函数修改
void PopulateGlobalStiffnessMatrixMultiThreaded(arma::sp_mat* global_stiffness_matrix, Mesh* mesh, Material* material, int number_of_cpus) {
    int total_elements = mesh->elements.n_rows;
    int step = total_elements / number_of_cpus;
    int start_element = 0;
    std::vector<std::future<arma::sp_mat>> futures;

    for (int i = 0; i < number_of_cpus; i++) {
        // 最后一个线程处理剩余所有元素，避免划分不均
        int current_count = (i == number_of_cpus - 1) ? (total_elements - start_element) : step;
        futures.push_back(std::async(std::launch::async, PopulateLocalStiffnessMatrix, mesh, material, start_element, current_count));
        start_element += step;
    }

    // 合并所有局部矩阵到全局
    for (auto& fut : futures) {
        *global_stiffness_matrix += fut.get();
    }
}

2. 优化线程数量与调度 #

• 线程数直接使用std::thread::hardware_concurrency()获取的系统物理核心数，避免过度调度；
• 不要额外多创建线程，让最后一个线程处理剩余所有元素，保证任务划分均匀。

3. 预分配稀疏矩阵内存 #

提前预估非零元素的数量，调用reserve()预分配内存，避免稀疏矩阵插入元素时频繁的内存重分配与拷贝：

// 全局矩阵预分配示例
int estimated_nonzeros = total_elements * 16; // 根据你的元素非零值数量调整
global_stiffness_matrix->reserve(estimated_nonzeros);

// 局部矩阵预分配示例
int local_estimated = current_count * 16;
local_mat.reserve(local_estimated);

4. 适配稀疏矩阵存储格式 #

Armadillo的sp_mat默认用*压缩稀疏列（CSC）*格式，如果你是按行写入数据，效率会较低。可以考虑：

• 局部矩阵用*压缩稀疏行（CSR）*格式构建，最后再转换为CSC格式；
• 直接指定全局矩阵的存储格式为更适配写入模式的类型。

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