SparseSparse for inverting structured sparse matrices

Specific Domains Numerics
1

Hello,

A few days ago, I discovered in this discussion that there is seemingly no way of inverting structured sparse matrices in Julia. We say that a matrix A is morally block diagonal if there are permutation matrices P,Q such that PAQ is block diagonal. Such matrices have sparse inverses and commonly occur in PDEs, particularly in the area of H-matrix preconditioners, etc…

Manually tracking the block structure (which may be “encrypted” by permutations P and Q) is tedious and error-prone, and unnecessary, since the usual sparse solvers can be adapted to invert such matrices efficiently. This is what module SparseSparse does. The name is because it’s a Sparse system with a Sparse right-hand-side.

I was surprised to see that my implementation is fairly efficient. On the problem A\b with a dense vector b, converting b to a sparse vector and then using completely sparse algorithms is about half the speed of the builtin solver for a dense vector b.

I’m not sure my code has the same quality as whatever is builtin to Julia but maybe it can be an inspiration and encouragement to have builtin inverses and solvers for sparse right-hand-sides?

Here is an example:

julia> using SparseSparse, SparseArrays
       A = sparse([2 3 0 0
                   4 5 0 0
                   0 0 6 7
                   0 0 8 9.0])
       inv(A)
4×4 SparseMatrixCSC{Float64, Int64} with 8 stored entries:
 -2.5   1.5    ⋅     ⋅ 
  2.0  -1.0    ⋅     ⋅ 
   ⋅     ⋅   -4.5   3.5
   ⋅     ⋅    4.0  -3.0

Hopefully this is helpful to someone.

6 Likes
2

Great! Note that the git repository should conventionally end with “.jl” for Julia modules, i.e. SparseSparse.jl

3

Since your implementation is fairly short, it might be worth submitting a PR to the SparseArrays package and get some feedback.

4

Hi Steven,

I think that’s a lovely idea. I’m not quite as dumb as I portray myself to be, but still I’m not 100% sure on how to do this. Do I have to spend several weeks setting up a development environment for Julia, like I had to do when I joined SGI in 1995, or how does it work? :slight_smile:

Cheers,

S

5

Shouldn’t take long to set up; the julia build might take an hour the first time, but generally doesn’t require manual intervention (assuming you have compilers etc. installed). See: