sparsemat / sprs

sparse linear algebra library for rust
Apache License 2.0
381 stars 45 forks source link

Example on solving a linear system? #335

Closed Makogan closed 9 months ago

Makogan commented 9 months ago

I have been trying to figure this one out for a while and it is not obvious from either the code nor the examples. I think either sprs or sprs_ldl would greatly benefit from having a single example where a simple linear system is setup and solved, to walk new users as to how to do this.

This seems important given that sparse matrices are primarily use for solving large systems of equations in the first place.

jlogan03 commented 9 months ago

I'd recommend taking a look at the suitesparse umfpack bindings. I haven't added an ergonomics layer for them yet, but the tests give an example of usage to solve a general square Ax=b system. I plan to add a BiCGSTAB implementation in the next month or so, and it will probably make sense to add more ergonomics features and examples at the same time.

https://github.com/sparsemat/sprs/blob/f03627f2517823293e189167d3b5b581748aca6a/suitesparse_bindings/sprs_suitesparse_umfpack/src/lib.rs#L328

You will need suitesparse installed in order to use that. The install process is pretty easy on ubuntu and mac.

Makogan commented 9 months ago

I will piggy back off of this comment, I think there might be a bug in how matrix multiplication is performed by the crate?

I have this snippet:

 let laplacian: CsMat<_> = triplets.to_csr();

    let mut dbg_check = HashMap::new();
    for (v, (r, c)) in laplacian.iter()
    {
        if dbg_check.contains_key(&(c, r))
        {
            let o = *dbg_check.get(&(c, r)).unwrap();
            assert!(o == v, "{} {} {} {}", v, o, r, c);
        }
        dbg_check.insert((r, c), v);
    }

    let mass: CsMat<_> = mass.to_csr();

    for (_, (r, c)) in mass.iter()
    {
        assert!(r == c);
    }

    let high_laplacian = &mass * &laplacian;

    let mut dbg_check = HashMap::new();
    for (v, (r, c)) in high_laplacian.iter()
    {
        if dbg_check.contains_key(&(c, r))
        {
            let o = *dbg_check.get(&(c, r)).unwrap();
            assert!(o == v, "{} {} {} {}", v, o, r, c);
        }
        dbg_check.insert((r, c), v);
    }

According to the first check, the matrix laplacian is symmetric. Accordfing to the mid check, the matrix mass is diagonal. But the last check fails with:

 thread 'convex_polytope::tests::test_topologize_skeleton' panicked at '0.2576181530716957 0.12738811865592223 512 0', src/differential_geometry.rs:270:13
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
error: test failed, to rerun pass `--lib`

So the product &mass * &laplacian is producing a non-symmetric matrix despite being the product of a diagonal and a syemmtric matrix?

How is that possible?

jlogan03 commented 9 months ago

Nothing immediately comes to mind. I created #336 from your comment to track this separately.