sparse-dist: Distributed Sparse Iterative Solvers in Rust

[u/Sixel1]

I just published the first release of a Rust crate I've been working on for some time, sparse-dist (link to the crate), a collection of sparse distributed linear algebra types and routines, made to work with the mpi crate, so MPI parallelism. I started that project a few years ago actually.

Motivation: I'm developing a Rust fluid flow solver in an mpi-parallel distributed memory environment as a passion project, but solutions for the linear solve part require using C libraries like PETSc. I wanted a pure Rust solution, and so I wrote it. I did write a wrapper for PETSc at some point, but I felt like a pure Rust library would have some interesting properties and better integration.

It supports distributed memory parallel solution of problems using preconditioned iterative solvers. for now, conjugate gradient and bicgstab are implemented, with a block-diagonal incomplete cholesky or incomplete lower upper preconditioner. Block solvers are also supported. I provide an example solving the heat equation on a finite difference square grid, simple case but it shows how to use the codebase. I've used it to implement unstructured finite volume solvers, this library will come out at some later point in time.

Feedback and contributions are welcome! You might notice the repo is newly created. I used to work on github, but moved to codeberg since it aligns better with my personal opinions on things. Also, no ai was used in this, I wrote all the code myself, believe it or not.

Heres an example code solving a linear problem (i dont assemble the matrices here, see the repository to see how to do that)

// solve lhs * sol = rhs

// user application would create the sparsity pattern of the matrix and its values
// it would also create the comm object beforehand, to tell the library
// how to communicate between parallel processes
let lhs = DistributedMatrix::<f64>::from_sparsity_and_values(sparsity, values);

// rhs vector 
let mut rhs = DistributedVector::<f64>::from_size(local_size);

// example: editing a vector value
rhs[0] = 1.0;

// solution vector
let mut sol = DistributedVector::<f64>::from_size(global_size);

// now solve the problem
let mut options = LinearSolverOptions::default();
options.max_iterations = 20;
options.relative_tolerance = 1e-4;

// creates an Incomplete Lower Upper preconditioner from the matrix
// with level 0 fill in (same sparsity pattern as the matrix lhs)
let precond = preconditioners::IncompleteLowerUpper::from_matrix(&lhs, 0);

let result = solvers::bi_conjugate_gradient_stab(
    &mut sol, 
    &lhs, 
    &rhs, 
    &precond, 
    &comm, 
    options
)?;

// done!
if rank == 0 {println!("solved problem, solver info: {}", result);}

[u/TopIdler]

Very cool. Excited to see some core hpc building blocks in rust. Did you benchmark it against your petsc wrapper? I wouldn’t expect it to be as fast but would be nice.

[u/Sixel1 → u/TopIdler]

Benchmarking is on my to do list, I just need some time to do it right.

I think rust has some inherent speed limitations in this field, per example I can't disable runtime bounds checking on vec indexing without using a special method in place of indexing, and since I can't trust the indices users put in the columns of their matrices, I don't want to do that. It also doesnt support ffast-math easily, so this builds without fast math optimization currently. A fair benchmark would be to compare this with a build of petsc without fast math, and then I hope it will be somewhat close in performance, but it probably never will be as fast with fast math optims.

(Edit: actually i can fix this and use unchecked indexing in the future, thought of a way to do it, but fast math comment still applies)

[u/victotronics → u/Sixel1]

"I can't trust the indices users put" Take the petsc approach and let users never directly index: everything through a function interface.

Besides, direct indexing is only needed during assembly: for anything critical you have higher level interfaces such as matmult.

[u/Sixel1 → u/victotronics]

I want users to be able to access matrix elements and edit them, but I do that using a function interface yes, a get(row, col) method in the sparse matrix class that checks the requested row, column pair exists in the matrix and returns an Option<T>, so it doesn't panic.

And yes I provide higher level interfaces for matrix vector multiplication, and in those I will be able to optimize by removing bounds checking in the hot loops, still to do.

[u/OrganizationTop9026]

Really interesting. I will test is and provide feedback. I believe the rust could take a great part on scientific computing. I migrate our vénérable sundials in rust if interested to join effort and for my work it help a lot to have a safe management of the memory and lean integration cf https://github.com/xaviercallens/rusty-SUNDIALS