Sparse matrix library

1

Can someone recommend a sparse matrix library in modern Fortran (i.e. from Fortran 90 onwards)? I need only basic stuff, like matrix multiplication or matrix vector multiplication.

To explain what I want to do:

To be more specific, I have to find the stationary distribution of a discrete state Markov chain. The stationary distribution is a row vector “mu” with N elements, the transition matrix “P” is a square matrix with N elements, all positive, with rows summing up to one. The stationary distribution mu satisfies the equation:

mu = mu*P

One way to find the stationary mu is to use fixed point iteration. So I will have to do a vector-matrix multiplication hundreds of times within a do while loop. I would like to take advantage of the fact that P is sparse.

Any suggestion is greatly appreciated!

2

Intel oneAPI MKL is not that difficult to use (https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2024-0/inspector-executor-sparse-blas-routines.html); the interface is kind of ugly though:


use mkl_spblas

! Sparse matrix data 
integer :: n, nnz
real(dp), allocatable, target :: values(:)
integer, allocatable, target :: ip(:), jp(:) 
type(SPARSE_MATRIX_T) :: p
type(MATRIX_DESCR) :: descr
integer :: stat

! Vectors
real(dp), allocatable :: mu1(:), mu2(:)

! Coefficients
real(dp) :: alpha, beta

n = ...
nnz = ... 

allocate(ip(n+1), jp(nnz))

! ... initialize CSR matrix graph ...

! Loop over rows to assemble matrix
do i = 1, n
    ipl = ip(i)
    ipu = ip(i+1) - 1
    values(ipl:ipu) = ...
end do

stat = mkl_sparse_d_create_csr (p, SPARSE_INDEX_BASE_ONE, &
                                n, n, ip(1:n), ip(2:n+1), jp, values)

descr%type = SPARSE_MATRIX_TYPE_GENERAL

alpha = 1.0_dp
beta = 0.0_dp

do while (.not. converged)

   ! mu_2 := P^T*mu_1
   stat = mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, alpha, P, descr, mu1, beta, mu2)
   ! mu_1 := P^T*mu_2
   stat = mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, alpha, P, descr, mu2, beta, mu1)

end do

! When done ...
stat = mkl_sparse_destroy(p)

Note I made use of the property, \mu^T = (\mu P)^T = P^T \mu^T. Two buffers are needed because the SpMV arguments should not alias.

3 Likes
3

Thank you so much for answer, @ivanpribec! I tried out your suggestion, but I’m not sure how to include “mkl_spblas” in my compile line. I’m using a windows pc with Intel oneAPI 2024.0. My fortran compiler is ifort and I have these flags:

/Qmkl /Qopenmp

. I can see the “mkl_spblas” files on my computer in “C:\Program Files (x86)\Intel\oneAPI\2024.0\include”. When I compile my code with the “use mkl_spblas” statement, I get the following error message:

error #7002: Error in opening the compiled module file. Check INCLUDE paths. [MKL_SPBLAS]
use mkl_spblas
--------^

4

Intel distributes the interface modules in source form, so they can be used also with other compilers.

The solution is to include it like this

include "mkl_spblas.f90"

module foo
use mkl_spblas
!...
end module

(This should only be done in one translation unit.)

2 Likes