Mistake in sparse linear system solution

Hi, I attempt to use sprs_ldl to solve a sparse Laplacian system. However, in the above toy case, the solution does not seem to be correct.

The toy case is to solve x for Ax=b, where A is the Laplacian matrix for a 33 grid G, i.e., A is 99. Since A is a Laplacian matrix, we have

A[i,i] = - count of the neighbors of G_idx(i) = G[j,k], i = j*3+k
A[i, i'] = 1, if G_idx(i') is a neighbor of G_idx(i)
          = 0, otherwise

b is set as [0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0]

The ground truth should be [x, x, x+1, x, x, x+1, x, x, x+1], however sprs_ldl's solution is [-8388606, -8388606, -8388607, -8388606, -8388606, -8388607, -8388605.5, -8388606.5, -8388607]

The code is

fn pressure_projection(d: &Vec<f32>) {
    let w = 3;
    let h = 3;
    let sz = 9;

    let mut pa = TriMat::<f32>::new((sz, sz));
    let mut pbi = Vec::with_capacity(sz);
    let mut pbv: Vec<f32> = Vec::with_capacity(sz);
    for i in 0..w {
        for j in 0..h {
            let mut cnt = 0;
            if i > 0 {
                pa.add_triplet(i*h+j, (i-1)*h+j, 1_f32);
                cnt += 1;
            }
            if i < w-1 {
                pa.add_triplet(i*h+j, (i+1)*h+j, 1_f32);
                cnt += 1;
            }
            if j > 0 {
                pa.add_triplet(i*h+j, i*h+j-1, 1_f32);
                cnt += 1;
            }
            if j < h-1 {
                pa.add_triplet(i*h+j, i*h+j+1, 1_f32);
                cnt += 1;
            }
            pa.add_triplet(i*h+j, i*h+j, -cnt as f32);

            pbi.push(i*h+j);
            pbv.push(*d.get(i*h+j).unwrap());
        }
    }
    let pa: CsMat<_> = pa.to_csr();
    let pb = CsVec::new(sz, pbi, pbv);
    let ldl = Ldl::default();
    let system = ldl.numeric(pa.view()).unwrap();
    let p = system.solve(pb.to_dense());
    println!("solution: {}", p);
}

fn main() {
    let d: Vec<f32> = vec![0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0];
    pressure_projection(&d);
}

The question is if this error is caused by bug, or there are mistakes in my code, or there is a way to construct a system with higher precision. Thanks in advance.

sparsemat / sprs

Mistake in sparse linear system solution #346