Mistake in sparse linear system solution
MarshalLeeeeee opened this issue · 1 comments
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.
My bad, the ground truth I provided is not correct. I will check my code.