# Sparse matrix library

**URL:** <https://fortran-lang.discourse.group/t/sparse-matrix-library/7653>\
**Category:** Uncategorized\
**Created:** [March 19, 2024, 9:27pm UTC](https://fortran-lang.discourse.group/t/sparse-matrix-library/7653 "2024-03-19T21:27:13Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [March 19, 2024, 9:27pm UTC](https://fortran-lang.discourse.group/t/sparse-matrix-library/7653/1 "2024-03-19T21:27:13Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [March 19, 2024, 10:01pm UTC](https://fortran-lang.discourse.group/t/sparse-matrix-library/7653/2 "2024-03-19T22:01:25Z")

</div>

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](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:

```fortran

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.

---

<div class="post-metadata">

**Author:** ![aledinola](https://avatars.discourse-cdn.com/v4/letter/a/3ec8ea/32.png) [@aledinola](https://fortran-lang.discourse.group/u/aledinola)\
**Post date:** [March 20, 2024, 12:10pm UTC](https://fortran-lang.discourse.group/t/sparse-matrix-library/7653/3 "2024-03-20T12:10:24Z")

</div>

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  
> --------^

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [March 20, 2024, 12:30pm UTC](https://fortran-lang.discourse.group/t/sparse-matrix-library/7653/4 "2024-03-20T12:30:38Z")

</div>

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

```auto
include "mkl_spblas.f90"

module foo
use mkl_spblas
!...
end module

```

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