Hi,
Can somebody please tell me if there is any Fortran equivalent to Suitesparse by Dr.Tim Davis Suitesparse. This is basically a library for handling sparse matrices. It is written in ‘c’. I am confused seeing all the pointer magics they do in ‘c’. So I am looking for a fortran alternative. Matlab uses this library under the hood to handle large sparse matrices.
If you want to solve AX = B in matlab, you can simply use A/B and you will get the X matrix. The ‘/’ operator is heavily overloaded with ‘blas’, ‘lapack’ and ‘suitesparse’ libraries - so that user need not worry when to call what.
Thank you.
5 Likes
I’m not familiar with using suitesparse, but here are some possible Fortran options:
Not all of them provide linear solvers I don’t think and it probably goes without saying that the Fortran options are much more complicated interfaces than MATLAB’s convenient operators.
2 Likes
Thank you.
Yeah, I don’t expect them to be as convenient.
Hello @Ashok,
based on your initial post, it seems that you are looking for direct sparse solvers. The ones I am aware of include:
Most of these libraries offer a very low-level interface and expect a matrix in either CSC or CSR format. Both of these formats are typically represented by a collection of three arrays (typically named something like ia, ja, and values). For easier management it makes sense to encapsulate these in a derived type along the lines of:
type :: csr_matrix_t
integer :: nrows, ncols
!! Dimensions of the matrix
integer :: nnz, maxnnz
!! Number of nonzero values, and size of storage arrays
integer, allocatable :: ia(:)
!! Row pointer of size nrows + 1
integer, allocatable :: ja(:)
!! Column pointer; only nnz values in use
real(dp), allocatable :: values(:)
!! Sparse matrix values; only nnz values in use
end type
Personally, I have been using the Intel MKL library for sparse linear algebra since Intel MKL offers a complete implementation of the Sparse BLAS standard. You can build a small sparse matrix library around the more recent inspector-executor Sparse BLAS routines. The operations supported are matrix-vector products, matrix addition, solution of triangular systems, and matrix-matrix multiplication.
Jack Dongarra maintains a more comprehensive list of freely available software for linear algebra at: http://www.netlib.org/utk/people/JackDongarra/la-sw.html
2 Likes
The derived type shown above is very similar to the record structure used to store sparse matrices in MATLAB or the CSparse library (part of SuiteSparse).
Specifically, for MATLAB,
… A sparse matrix is a C record structure with the following constituents. The nonzero elements are stored in a one dimensional array of double precision reals, in column major order. (If the matrix is complex, the imaginary parts are stored in another array.) A second array of integers stores the row indices. A third array of n+1 integers stores the index into the first two array of the leading entry in each of the n columns, and a terminating index whose value is nnz. Thus a real m \times n matrix with nnz nonzeros uses nnz reals and nnz + n + 1 integers.
The sparse data structure is allowed to have unused elements after the end of the last column of the matrix. Thus an algorithm that builds up a matrix one column at a time can be implemented efficiently by allocating enough space for all the expected nonzeros at the outset.
The quotes are taken from the paper Sparse Matrices in MATLAB: Design and implementation.
Thanks a lot @ivanpribec .
The collection that Jack Dongarra has, many are written in Fortran. Unfortunately they are not modernized. That’s why c++ took up.
qr_mumps demonstrates superior performance to SuiteSparseQR and is written in modern Fortran.
2 Likes
Yousef Saad’s SPARSKIT2 software is still available. It’s only Fortran 77 but I think
John Burkhardt has an F90 translation on his site. It has the basic sparse BLAS routines along with translation to/from various sparse matrix formats. Also, some of his
Newton-Krylov (GMRES etc) routines if I remember correctly. See his book
Interative Methods for Sparse Linear Systems, PWS publishing
for more details.
PDFs of both the first and second editions are also available from his web site (for free)
https://www-users.cs.umn.edu/~saad
You can get Saad’s version from
https://www-users.cs.umn.edu/~saad/software/SPARSKIT/index.html
or
wget http://www-users.cs.umn.edu/~saad/software/SPARSKIT/SPARSKIT2.tar.gz
Modern does not automatically means better. I benchmarked many of libraries mentioned and the winner was y12m.
2 Likes
@Jamie , @rwmsu , @themos, @Beliavsky
Thank you. I will go through them and see what suits best to my use case.
Do you have a public repository with the benchmarks you tested?
I think the use of y12m could be simplified if the workspace arguments were allocated internally or hidden in a derived type. I’ve opened a few issues with ideas for improvements in my fpm-compatible mirror of the code.
Y12m is kind of spaghetti code but compilers tend to optimize it very well, like Linpack. Would be interesting to see if modernization kills performance as its was the case around 2000, when compilers poorly supported modern fortran.
1 Like
How about Y12m compared with MA48? Thanks!
For my own research, I developed this OO library for managing sparse matrices. Note that I never benchmarked, but it supports COO and CRS, as well as different solvers, like CG and PARDISO, and it is flexible enough for my work 
Here is an example:
...
type(coosparse)::coo
type(crssparse)::crs
!COO
open(newunit=iunit,file='nput_sparsematrix.ascii',status='old',action='read')
read(iunit,*) nrow
coo = coosparse(5, lupper=.true.)
call coo%setsymmetric()
do
read(iunit,*,iostat=istat) row,col,val
if(istat.ne.0)exit
call coo%add(row,col,val)
end do
close(iunit)
call coo%add(1,1,10._wp)
!CSR UPPER
crs = coo
call crs%printsquare()
! Get permutation vector using Metis
perm=crs%getordering(compress=1,ctype=METIS_CTYPE_SHEM)
...
call crs%setpermutation(perm)
...
call crs%printstats()
...
call coo%cg(x, y) !using own CG algorithm
...
call crs%solve(x,y) !using PARDISO
...
4 Likes