# SparseSparse for inverting structured sparse matrices

**URL:** <https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715>\
**Category:** Numerics\
**Created:** [November 2, 2023, 6:08pm UTC](https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715 "2023-11-02T18:08:10Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![loisel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/loisel/32/50626_2.png) [@loisel](https://discourse.julialang.org/u/loisel)\
**Post date:** [November 2, 2023, 6:08pm UTC](https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715/1 "2023-11-02T18:08:10Z")

</div>

Hello,

A few days ago, I discovered in [this discussion](https://discourse.julialang.org/t/block-diagonal-and-sparse-matrices/105551/14) 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](https://github.com/sloisel/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
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 2, 2023, 6:53pm UTC](https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715/2 "2023-11-02T18:53:17Z")

</div>

> [@loisel](#):
>
> This is what module [SparseSparse](https://github.com/sloisel/SparseSparse) does.

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 2, 2023, 6:57pm UTC](https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715/3 "2023-11-02T18:57:49Z")

</div>

> [@loisel](#):
>
> 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?

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

---

<div class="post-metadata">

**Author:** ![loisel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/loisel/32/50626_2.png) [@loisel](https://discourse.julialang.org/u/loisel)\
**Post date:** [November 2, 2023, 7:11pm UTC](https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715/4 "2023-11-02T19:11:29Z")

</div>

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? 🙂

Cheers,

S

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 27, 2023, 2:48pm UTC](https://discourse.julialang.org/t/sparsesparse-for-inverting-structured-sparse-matrices/105715/5 "2023-12-27T14:48:33Z")

</div>

> [@loisel](#):
>
> 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? 🙂

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:

> <https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md>