Rust并行化牛顿-拉夫逊法无性能提升问题排查

问题1：为何并行版和非并行版速度一样慢？ #

核心原因是全局互斥锁完全阻塞了并行执行。你的代码中，每个线程启动后立即锁定整个矩阵atom_a和向量atom_b，导致所有线程必须串行执行——第一个线程拿到锁后，其他线程全部处于等待状态，直到它完成任务释放锁，下一个线程才能继续。这种情况下，并行代码的执行逻辑和串行完全一致，还额外增加了线程创建、切换和锁的开销，甚至可能比串行更慢。

问题2：实现中哪里出错导致性能无提升？ #

1. 错误的锁粒度：用Mutex包裹整个矩阵和向量，完全消除了并行性。高斯消元的并行化需要细粒度的并发控制，而非全局锁。
2. 并行划分逻辑不符合算法依赖：高斯消元的消去阶段，第k步的操作依赖第k行的结果，你直接将行范围平均分给线程，但不同线程处理的行可能存在未解决的依赖关系（比如线程A处理i=0的行，线程B处理i=1的行，但i=1的消去需要i=0的行已处理完成），导致并行逻辑无效。
3. 数据副本未同步回原结构：线程修改的是atom_a和atom_b的副本，但最后回代阶段使用的是原始的a和b，这直接导致结果异常（极端值就是因为用了未消元的原始矩阵进行回代计算）。
4. 不必要的内存开销：derivada_parcial每次调用都复制两次整个输入向量，对于大维度（如TAM=1000）的向量，会产生大量内存分配和复制开销。

问题3：如何优化才能获得性能收益？ #

1. 修复并行高斯消元的核心逻辑 #

高斯消元的可并行阶段是：对于第k步，所有行号j > k的行可以并行处理（它们仅依赖第k行的结果，彼此之间无依赖）。正确的做法是按步骤k循环，每一步先完成主元选择（必须串行），再并行处理所有j > k的行。

2. 使用Rayon库简化并行实现 #

手动管理线程容易出错，推荐使用Rust生态的rayon库，它的并行迭代器能自动处理负载均衡和线程调度，大幅简化代码。

首先在Cargo.toml添加依赖：

[dependencies]
rayon = "1.8"
rand = "0.8"

修正后的并行高斯消元函数：

use rayon::prelude::*;

pub fn par_elim_gauss(tam: usize, a: &mut [Vec<f64>], b: &mut [f64]) -> Vec<f64> {
    // 主元选择+消去阶段
    for k in 0..tam {
        // 串行选择主元行
        let max_linha = (k..tam)
            .max_by(|&i, &j| a[i][k].abs().partial_cmp(&a[j][k].abs()).unwrap())
            .unwrap();
        if max_linha != k {
            a.swap(k, max_linha);
            b.swap(k, max_linha);
        }

        // 提取第k行的快照，避免并行修改时的冲突
        let row_k = a[k].clone();
        let b_k = b[k];
        let pivot = a[k][k];

        // 并行处理所有j > k的行
        (k+1..tam).into_par_iter().for_each(|j| {
            let m = a[j][k] / pivot;
            for col in k..tam {
                a[j][col] -= m * row_k[col];
            }
            b[j] -= m * b_k;
        });
    }

    // 回代阶段（串行，因为x[i]依赖x[i+1..tam]的结果）
    let mut x = vec![0.0; tam];
    for i in (0..tam).rev() {
        let soma: f64 = (i+1..tam).map(|j| a[i][j] * x[j]).sum();
        x[i] = (b[i] - soma) / a[i][i];
    }
    x
}

3. 修复结果异常问题 #

之前的代码中线程修改的是数据副本，回代用原始数据导致错误。上面的修正代码直接修改传入的a和b，回代阶段使用消元后的矩阵，结果会正确。同时将f32替换为f64，提升数值精度，减少溢出和极端值的出现。

4. 优化导数计算的内存开销 #

修改derivada_parcial，只复制一次向量并复用，减少内存分配：

pub fn derivada_parcial(funcao: impl Fn(&[f64]) -> f64, indice: usize) -> impl Fn(&[f64]) -> f64 {
    let h: f64 = 0.001;
    move |x: &[f64]| {
        let mut x_temp = x.to_vec();
        // 计算x+h的函数值
        x_temp[indice] += h;
        let f_mais = funcao(&x_temp);
        // 调整为x-h的状态，复用同一个向量
        x_temp[indice] -= 2.0 * h;
        let f_menos = funcao(&x_temp);
        (f_mais - f_menos) / (2.0 * h)
    }
}

5. 调整主函数的类型匹配 #

将主函数中的f32全部替换为f64，与修正后的函数保持一致：

use rand::Rng;
use rayon::prelude::*;

use crate::teste::{derivada_parcial, par_elim_gauss};

mod teste;

fn main() {
    const TAM: usize = 1000;
    let a: &mut[Vec<f64>] = &mut vec![vec![0.0; TAM]; TAM];
    let x_0: &[f64] = &[1.0; TAM];
    let b: &mut [f64] = &mut [0.0; TAM];

    for i in 0..TAM {
        for j in 0..TAM {
            let potencia: f64 = rand::thread_rng().gen_range(1.0..=6.0);
            let constante: f64 = rand::thread_rng().gen_range(1.0..=10.0);
            let aleatorio_i: i32 = rand::thread_rng().gen_range(0..=1);

            let funcao = |x: &[f64]| (-1.0_f64).powi(aleatorio_i) * constante * x[j].powf(potencia);
            let derivada = derivada_parcial(funcao, j);

            let valor = derivada(x_0);
            a[i][j] += valor;
            b[i] += funcao(x_0) * (-1.0);
        }
    }

    let resultados = par_elim_gauss(TAM, a, b);
    println!("Resultado: {:#?}", resultados);
}

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