Sparse matrix support about stdlib HOT 57 OPEN

from stdlib.

from stdlib.

Before trying to reinvent the wheel check out PSBLAS 3 which is a parallel Fortran 2003 OOP
implementation of sparse matrixes (see http://people.uniroma2.it/salvatore.filippone/psblas and the
selected publications listed).

from stdlib.

Hi @certik and @jvdp1 . Thanks @ivan-pi for making me aware of the effort.
I'm a former PhD student of @sfilippone (and author of http://librsb.sourceforge.net/) and can only second him.

Said that, seeing that the discussion is taking time, please let me share a bit of my experience, too.

Concerning the API:

• I see https://stdlib.fortran-lang.org/module/stdlib_linalg.html ; you are striving for modern OO-oriented Fortran API for "Sparse Basic Linear Algebra" functionality; I'm not sure if this has been tackled as such after 1999; surely it has been tackled in @sfilippone 's work within PSBLAS though (so, with the "MPI-Parallel" additional aspect) (and countless other projects, but without the "Sparse BLAS" generality aim)
• wrt the 1999 Sparse BLAS Technical Forum specification -- http://www.netlib.org/blas/blast-forum/chapter3.pdf -- it contains a lot of wisdom (maybe foremost: format independence)
• I've worked on allowing usage of that particular API (http://librsb.sourceforge.net/doc/1.3.0.0/namespaceblas__sparse.html) -- as a wrapper around the LIBRSB's own API (http://librsb.sourceforge.net/doc/1.3.0.0/rsb_8F90.html)
• the Sparse BLAS API is node-level (so, more limited in scope than PSBLAS) and format-independent, though it's trivial to have a format specified (may use ussp "set property" function for extension -- see http://librsb.sourceforge.net/doc/1.3.0.0/namespaceblas__sparse.html#a469df92a4d25a9554fb1d79cdac1de84 -- this is the trick that even in LIBRSB, allows some control over the format)
• reason why I cite the 1999 Sparse BLAS API is that in principle, one may wrap around it if needed

Whatever you end up choosing for the most external API later, I recommend caution with respect to the internals, especially if you start an effort from scratch (see above @certik @sfilippone 's comments about sorting :-) ).
It's 2022, and if at some point the user discovers that at least multithreaded parallelism would be desirable to have, then things get complicated easily.
For instance, one may discover that "simple" COO/CSR kernels are not enough to be used within an OpenMP-parallel region, and more offsets or strides are needed to be specified (see e.g. any of the functions here
https://salsa.debian.org/science-team/librsb/-/blob/master/rsb_krnl_bcss_spmv_u.c?expanded=true&viewer=simple --- even if expressed in modern C++ -- see https://salsa.debian.org/science-team/librsb/-/blob/master/librsbpp/rsbpp.hpp#L2064 -- it remains complicated business).

Staying flexible with respect to the internals, perhaps engineering a good API+wrappers code generator that would interface with an existing library is maybe the best thing.

I'm not discouraging from offering a sane COO/CSR interface to the user -- but maybe joining forces with one or more existing projects, maybe even reusing the internals, may help getting something usable very quickly.

And maybe, kick more library implementors in making their libraries substitutable one-to-one for stdlib_linalg..

Michele

from stdlib.

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.

from stdlib.

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

from stdlib.

Ok @certik, the representations of the sparse matrices are certainly in the core of this work, I will start writing some text about the objectives of the API, its ambition and philosophy (e.g. to have a minimum number of functions/subroutines instead of being extensive as I envision now) and in a more technical section I will start writing about these two formats to start with.
Tomorrow I will have something to show.

from stdlib.

from stdlib.

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.

from stdlib.

@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.

from stdlib.

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

from stdlib.

@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.

from stdlib.

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?

from stdlib.

@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.

from stdlib.

@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.

from stdlib.

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)

from stdlib.

@certik What is the status of this discussion? From the previous posts I understand that there is a number of related projects, some discussion about technical aspects but it is not clear to me if there is a work started to create this library in stdlib (or even have people decided that this proposal for a sparse matrix library has passed step 1 of the Workflow in workflow.md).

Sorry for what I am missing.
Williams

from stdlib.

@ghwilliams the status of the discussion at this issue as I understand it is that there is a general agreement that we want this, and so now we need to start implementing it, say the COO / CSR subset to get started, then submit a PR with a draft of a specification and then people will comment if they like it or not.

It should be designed taking into account the discussion above and all the other libraries as linked at the top post.

If you want to go ahead to take a lead on this one, that would be great. I'll help out as much as I can.

from stdlib.

Ok @certik . I will try to direct the discussion toward more concrete conclusions. So I will follow the workglow.md guide line and move on to Step 2. Defining an API.
So, I will post here later some ideas for the API. Anyone else interested could contribute as well. Let's just keep the focus now, as I said, on step 2 of the guideline because I think the step one was already confirmed unless someone comes and say "Not too fast Williams".
Let's move.

regards
Williams

from stdlib.

@ghwilliams I agree, let's move! Good to go forward with step 2 (API) as far as I'm concerned.

I'm happy you're here.

from stdlib.

from stdlib.

Good point @rweed.

@certik What's missing in the 4 existing Fortran implementations that you listed? Or do you suggest that one of the existing implementations is provided under stdlib?

Perhaps naive questions, but I'm not too familiar with this domain.

from stdlib.

from stdlib.

I'm going to keep a document for the specification of the API. I will collect all the suggestions posted here for the API. In a first round I will include many suggestions for data structures/algorithms in the previous posts, also taking the other libraries as a reference. After some time I will propose to filter this list to get a final list of features to go in a version 1 of the API.
It will be good also to decide later on a deadline for the conclusion of step 2. In this way we can have something material after a finite amount of time.

Just as a side note, initially I'm planning to keep Object Oriented solutions out and focus on a clean, procedural design. The reason for that is to notenter an infinite loop of discussions about the pros and cons of one design or another. Let's choose one and move on.

I plan to write the API document using MS Word but if anyone has another suggestion for writing the text (latex, notepad, whatever), for me it is ok, just let me know.

Williams

from stdlib.

In general most people agree to do both, have a low level procedural (no side effects) design, and then have an optional higher level OO design.

I would recommend to start with a serial design for COO and CSR matrices and let's start the discussion on the API for those.

I updated all the Fortran implementations in the comment above that were suggested so far.

Most of them are either fully OO, or partially OO. The only procedural implementations that I found are:

• hfsolver (mine)
• Yale
• psblas3 (only some parts, e.g. sorting seems to be implemented in a procedural way)

So those three can serve as an inspiration for the low level API. The rest of the packages can serve as an inspiration for the OO API.

from stdlib.

@ghwilliams perfect, looking forward!

from stdlib.

Hello everybody.
I started writing the API document. I'm sending both a MS Word and pdf file. Although it doesn't have too much written it certainly will serve as a concrete place to develop the ideas for the API. anyone who wants to contribute can edit directly the Word file.
Discussing ideas for the API having this document as a convergence point certainly will make the discussions more oriented toward our goals.
Stdlib Sparse matrix API.docx
Stdlib Sparse matrix API.pdf

Of course any suggestions, critiques are welcome.

from stdlib.

Thanks @ghwilliams for starting it. Would it be ok to use Markdown and put the document on our wiki? That way others can easily edit and it's easier to view, and later it will be easier to create the specification document, which also uses Markdown.

Regarding indexing, in Fortran the standard is to start from 1, so that's what I would recommend.

from stdlib.

Of course, having the document in the wiki is much better just let me know how to do that.

from stdlib.

Here is the API that I propose for serial COO and CSR functionality: #189

from stdlib.

@ghwilliams click on the Wiki tab (https://github.com/fortran-lang/stdlib/wiki) and click on the green "New Page" button to create a new page, and create a page and then keep editing it.

from stdlib.

Ok. Sohuld I merge with the information in #189?

from stdlib.

@ghwilliams in #189 we have to discuss as a community if this is an API that we would like to use, or if not, how to improve it. I submitted the code so that we have something more concrete to discuss and so that people can compile and run it and play with it.

from stdlib.

Allright, how #189 then relates to this document that Im creating? Because in my document we are going to propose an API with the same purpose right?

from stdlib.

@ghwilliams yes, exactly. The PR #189 provides code implementation of an API that you are writing the specs for in the document. In order to efficiently discuss the API, it is helpful to have code that implements it, so that we can play with it.

from stdlib.

Ok I will continue then writing the API document in the wiki. It is good to keep track of all the posts here related to it.

from stdlib.

from stdlib.

I created the wiki page with all the information in the previous Word file. But there are some formatting issues to be corrected. I've never used Markdown before and I don't know about its capabilities to render math.
Anyway, I will continue working on the API from there.

from stdlib.

I would like to suggest that as policy COO will be used as the bridge format for implementing new sparse formats, for initializing an existing format and conversion to/from other formats. I believe
IMSL does this for their sparse matrix support and if I remember correctly so does PSBLAS. The
advantage of this is people implementing a new format only have to write code to go to/from COO
(or whatever format is selected) instead of having to support direct conversion to all the other formats that might eventually be impemented. I'm guessing this is standard practice in most other sparse matrix packages.

from stdlib.

from stdlib.

Yes SPARSKIT does the same. I will add some text about it.

from stdlib.

@certik I would like to include an image in the wiki page. Clicking the wiki edit tool asks me for a url but the file is in my local computer. The alternative, I think, is to have the image file uploaded to the repository but I don't have write access to it, I guess.

from stdlib.

from stdlib.

from stdlib.

I added (wiki page) some ideas for the relation between the storage format and the routines in the API.
I'm still working on that.

from stdlib.

I don't think it has been linked previously, but on netlib there is a folder of templates for sparse solvers:

• https://www.netlib.org/templates/

Section 4.3.1 of the documentation provides a survey of sparse storage formats.

from stdlib.

I know that book. If you can look at the wiki and leave some comments it will be greatly appreciated.

Thanks.

from stdlib.

from stdlib.

Thanks @sfilippone .

from stdlib.

Today I found this interesting work: https://www.math.uzh.ch/pages/spam/
Good to add to the list of references.

from stdlib.

@sfilippone thanks for the information. As I wrote in #189, I also suggest we do not default to quicksort.

from stdlib.

I've been editing the wiki (https://github.com/fortran-lang/stdlib/wiki/Stdlib-Sparse-matrix-API) and I would like to have some feedback to know if I am in the right direction. Any comments are welcome.
Among many things that need to be incorporated I would like to add state variables for the derived types representing the matrices, similar to what is done in PSBLAS (build, assembled, update).

from stdlib.

@ghwilliams thanks a lot for creating the page and suggesting an API.

My main feedback is that it looks like this would be the high level API -- or middle level, depending on how we want to view it, but not the low level API. Because it uses a derived type to store the CSR matrix, which I also ended up doing in my application, but I suggest for the low level API to not do that, unless absolutely necessary (which it isn't in this case), because you are then forcing your derived type on all applications. As an example, I would not be able to just use a routine from stdlib in my code, because I already use a different derived type for the CSR matrix. In the same way, unless PSBLAS moves to this derived type, it would not be able to use the subroutines either.

So for that reason, I propose for the low level API to be roughly in the sense of #189, and then the high level API roughly what you wrote in your wiki, or what PSBLAS and other libraries are using.

from stdlib.

Thanks @certik . I will work on your suggestions.

from stdlib.

I separated the design now in low and middle level parts. Low level subroutines take plain arrays of integers/real/complex values for the elements in the matrix representation.
But I would like to discuss here one particular aspect of the design of this API. It was suggested, and I'm following this suggestion, to have a preferred matrix format against with all the API core subroutines would be written and then have conversion routines for the other formats. The drawback of this approach, as far as I can see now, is that the data will be duplicated requiring, broadly speaking, two times RAM during the API calls.

from stdlib.

@sfilippone I will look at psblas (and other libs) to see how this problem is approached.
I don't like the idea of duplicating data/code so I will try to find a good solution that will not blamed by users regarding its performance nor by developers for having to write many versions of a routine for each data representation..

from stdlib.

Thank you @michelemartone for your comments.

Staying flexible with respect to the internals, perhaps engineering a good API+wrappers code generator that would interface >with an existing library is maybe the best thing.

I'm not discouraging from offering a sane COO/CSR interface to the user -- but maybe joining forces with one or more existing >projects, maybe even reusing the internals, may help getting something usable very quickly.

I mostly agree with all your comments. At this stage, I believe that a good (low-level) API is the first thing to define. And the sparse BLAS API might be a good start.
Re: the internals, joining forces would be really nice. I was not aware of your library. I will look at it.

from stdlib.