Unexpected results using `faer_core::solve::solve_lower_triangular_in_place_with_conj`

First of all, thanks a lot for your work in faer, it is impressive. It has also been a lot of fun to play around with the library.

Describe the bug The results of solving a lower triangular system are unexpected.

To Reproduce The following code reproduces the error. It also shows how the equivalent code with nalgebra gets the correct solution:

# Cargo.toml

[package]
name = "test"
version = "0.1.0"
edition = "2021"

[dependencies]
faer-core = "0.9.1"
nalgebra = "0.32.2"

// main.rs

const RESPONSE: [f64; 3] = [-18.0, -60.0, -80.0];

// Toeplitz and lower triangular
fn r_element(i: usize, j: usize) -> f64 {
    if i >= j {
        let diff = i - j;
        if diff < RESPONSE.len() {
            RESPONSE[diff]
        } else {
            0.0
        }
    } else {
        0.0
    }
}

fn r_matrix_faer(n: usize) -> faer_core::Mat<f64> {
    faer_core::Mat::with_dims(n, n, |i, j| r_element(i, j))
}

fn r_matrix_nalgebra(n: usize) -> nalgebra::DMatrix<f64> {
    nalgebra::DMatrix::from_fn(n, n, |i, j| r_element(i, j))
}

// Solve for X in Y = R * X
fn solve_x_faer(r: &faer_core::Mat<f64>, y: &mut faer_core::Mat<f64>) {
    faer_core::solve::solve_lower_triangular_in_place_with_conj(
        r.as_ref(),
        faer_core::Conj::No,
        y.as_mut(),
        faer_core::Parallelism::None,
    );
}

// Solve for X in Y = R * X
fn solve_x_nalgebra(r: &nalgebra::DMatrix<f64>, y: &mut nalgebra::DMatrix<f64>) {
    r.solve_lower_triangular_mut(y);
}

fn main() {
    let (i, j) = (400, 1);

    let r = r_matrix_faer(i);
    let mut x = faer_core::Mat::zeros(i, j);
    x.write(0, 0, 150.0);

    let mut y = r.clone() * x;

    solve_x_faer(&r, &mut y); // The values in `y` blow up

    // Now, the same thing with nalgebra
    let r = r_matrix_nalgebra(i);
    let mut x = nalgebra::DMatrix::zeros(i, j);
    x[(0, 0)] = 150.0;

    let mut y = r.clone() * x;

    solve_x_nalgebra(&r, &mut y); // The values in `y` are correct
}

Expected behavior I would expect the same results as nalgebra

sarah-ek / faer-rs

Unexpected results using `faer_core::solve::solve_lower_triangular_in_place_with_conj` #45