Sparse matrix support

[certik]

Prior art in other languages:

• SciPy: https://docs.scipy.org/doc/scipy/reference/sparse.html
• Matlab: https://www.mathworks.com/help/matlab/ref/sparse.html
• Julia: https://docs.julialang.org/en/v1/stdlib/SparseArrays/index.html

In Fortran (I keep this list updated with all implementations posted in this issue):

• @certik: https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/sparse.f90
• @jvdp1: https://github.com/jvdp1/libsparse
• @sfilippone: https://github.com/sfilippone/psblas3
• @victorsndvg: https://github.com/fempar/fempar/tree/experimental/Sources/Lib/LinearAlgebra/SparseMatrix
• @danshapero: https://github.com/danshapero/sigma
• https://github.com/cp2k/dbcsr
• https://github.com/TRIBO-Pprime/MSOLV
• Yale sparse matrix package: http://www.netlib.no/netlib/ode/yale.f

I really like the SciPy simple non-OO implementation, and I have ported it to modern Fortran in the link above. If people also want an OO implementation, then it can be build on top as an option.

One thing that I found out is that one must sort the indices and the overall speed very much depends on how quickly one can sort it. I ended up using quicksort, but it might be even faster to use some specialized sorting algorithm (such as Timsort) because in practice, the indices have subsections that are already sorted (typically coming from some local to global mapping as in finite elements), but overall it is not sorted.

[jvdp1]

A non-OO implementation could be an easier start indeed.
An issue we may need to discuss about sparse matrices is that NNZ may be larger than integer(4).
Also we will need to discuss the formats to support. COO and CRS3 could be a good start.

One thing that I found out is that one must sort the indices and the overall speed very much depends on how quickly one can sort it. I ended up using quicksort, but it might be even faster to use some specialized sorting algorithm (such as Timsort) because in practice, the indices have subsections that are already sorted (typically coming from some local to global mapping as in finite elements), but overall it is not sorted.

In my field, entries are added to a sparse matrix in a quite random way. Therefore, the indices are far to be sorted, even in a row. Therefore, I use quicksort as implemented in LAPACK.

[certik]

@jvdp1 eventually we can support all the formats like SciPy does. It's quite a bit of work, and so for my own work I just did COO and CSR, since I didn't have time to implement the others, but for stdlib, I think we have the manpower to implement all useful formats.

[rweed]

For sparse matrix support I suggest you look at the PSBLAS3 project of Filippone and Buttari. Its Object-Oriented Fortran 2003/2008. Go to

https://github.com/sfilippone/psblas3

I also recommend you read the papers listed there for ideas about implementing sparse matrix support in OO Fortran and the various "design patterns" associated with sparse matrices.

I like there approach to supporting different storage types. COO is used as input and the "bridge" type to all other storage formats. You only write converters to/from COO and your chosen storage format. Reduces the explosion of code trying to go directly from say CRS and other formats. For a "simple" sparse grid implementation, I use an approach similar to Java Sparse Arrays which is sorta like a CRS. I find it makes Matvecs simple

[zbeekman]

@rweed Thanks for the suggestion! I agree, and have been using PSBLAS in a number of projects for work. I added the suggestion to the list of prior art at the top.

[zbeekman]

Also, I'm not really sure where this class of sparse linear-algebra routines would belong, but I often see very large tridiagonal, and pentadiagonal systems. To my (admittedly limited) knowledge, there are not many good parallel algorithms for the solutions of these systems and there are not parallel algorithms in extant linear algebra libraries. The best/most promising algorithm that I know of is the SPIKE algorithm, but I've never seen a library that included it. I have a half-decent version specialized for tridiagonal matrices, and would love to see this and any other approaches implemented in the linear algebra portion of the stdlib.

Should I open a new issue for this, or does it fit within the general sparse matrix support category?

[sfilippone]

As the author of PSBLAS, I am very much aware of the need to sort indices, and this is taken care *internally* in my library, so that there is no burden on the user. And no, quicksort is not the best choice. I normally use a special version of merge sort, adapted from TAOCP. Hope this helps Salvatore

…

On Thu, Jan 2, 2020 at 5:04 PM Ondřej Čertík ***@***.***> wrote: Prior art: - SciPy: https://docs.scipy.org/doc/scipy/reference/sparse.html - Matlab: https://www.mathworks.com/help/matlab/ref/sparse.html - Julia: https://docs.julialang.org/en/v1/stdlib/SparseArrays/index.html - @certik <https://github.com/certik>: https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/sparse.f90 - @jvdp1 <https://github.com/jvdp1>: https://github.com/jvdp1/libsparse - @sfilippone <https://github.com/sfilippone>: https://github.com/sfilippone/psblas3 There are probably many more implementations. I really like the SciPy simple non-OO implementation, and I have ported it to modern Fortran in the link above. If people also want an OO implementation, then it can be build on top as an option. One thing that I found out is that one must sort the indices <https://github.com/certik/hfsolver/blob/b4c50c1979fb7e468b1852b144ba756f5a51788d/src/sparse.f90#L185> and the overall speed very much depends on how quickly one can sort it. I ended up using quicksort, but it might be even faster to use some specialized sorting algorithm (such as Timsort <https://en.wikipedia.org/wiki/Timsort>) because in practice, the indices have subsections that are already sorted (typically coming from some local to global mapping as in finite elements), but overall it is not sorted. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#38?email_source=notifications&email_token=AD274TY6ANZGMIDAYNGDX5DQ3YNANA5CNFSM4J6NGRLKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4ICGREGQ>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AD274TZUQE4LA53J5SW6ASTQ3YNANANCNFSM4J6NGRLA> .

[certik]

Thanks @sfilippone for the feedback! Indeed, I also suspect quicksort is not the best choice.

@zbeekman, I would start with COO, CSR that we all agree we need. Then as we are figuring out the API, let's think how tridiagonal (Lapack has those) and pentadiagonal systems fit in. Also, initially I was thinking of sticking to serial implementations.

We'll have to figure out how to best tackle parallel algorithms in stdlib, but I think that's an issue on its own, so I created #66 for it.

[zbeekman]

@certik good idea. Do you happen to know the name of the Lapack procedure for diagonal matrices? I remember looking for a good serial implementation too and not finding it, but perhaps I was only focused on parallel implementations.

[certik]

@zbeekman yes, I have used dstevd for eigenvalues of a tridiagonal matrix. If you just need a regular solve, then for example dptsv will do it.

[victorsndvg]

I don't know if this could be useful, but in FEMPAR project we have implemented an OO extendible sparse matrix based on Filipone PSBLAS sparse matrix. It's complicated ... but performs quite well.

https://github.com/fempar/fempar/tree/experimental/Sources/Lib/LinearAlgebra/SparseMatrix

[ivan-pi]

Then as we are figuring out the API, let's think how tridiagonal (Lapack has those) and pentadiagonal systems fit in. Also, initially I was thinking of sticking to serial implementations.

Pentadiagonal matrices can be treated with the banded solvers in LAPACK, but they require a special storage scheme. I have some simple examples available:

• Banded (four diagonals): https://gist.github.com/ivan-pi/a9c905065af75362c786c2032ac48c56
• Tridiagonal: https://gist.github.com/ivan-pi/9d7af8257392ec7340075d63c738fd79

The ancient Yale sparse matrix package has routines for ordering and solving both systems of symmetric and asymmetric sparse systems of equations. I don't know if the algorithms they used were any good though.

A few other Fortran sparse matrix libraries include:

• sigma by @danshapero
• dbcsr
• msolv is a common API to several direct sparse matrix solvers (UMFPACK, SuperLU, MUMPS)