Function solve_lu

src/polynomial.rs:112–150

(n: usize, lhs: &[f64], rhs: &[f64])

}

pub fn solve_lu(n: usize, lhs: &[f64], rhs: &[f64]) -> Vec<f64> {
    let mut lu = vec![0f64; n * n];
    let mut sum;
    for i in 0..n {
        for j in i..n {
            sum = 0.;
            for k in 0..i {
                sum += lu[i * n + k] * lu[k * n + j];
            }
            lu[i * n + j] = lhs[i * n + j] - sum;
        }
        for j in (i + 1)..n {
            sum = 0.;
            for k in 0..i {
                sum += lu[j * n + k] * lu[k * n + i];
            }
            lu[j * n + i] = (1. / lu[i * n + i]) * (lhs[j * n + i] - sum)
        }
    }

    let mut y = vec![0.; n];
    for i in 0..n {
        sum = 0.;
        for k in 0..i {
            sum += lu[i * n + k] * y[k];
        }
        y[i] = rhs[i] - sum;
    }

    let mut x = vec![0.; n];
    for i in (0..n).rev() {
        sum = 0.;
        for k in (i + 1)..n {
            sum += lu[i * n + k] * x[k];
        }
        x[i] = (1. / lu[i * n + i]) * (y[i] - sum);
    }
    x
}