Suitesparse equivalent in Fortran

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.

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

UMFPACK is part of SuiteSparse and has a Fortran interface.

1 Like

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.

1 Like

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 :wink:
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