Sparse algebra support with OOP API

doc/specs/stdlib_sparse.md

@@ -0,0 +1,243 @@
+---
+title: sparse
+---
+
+# The `stdlib_sparse` module
+
+[TOC]
+
+## Introduction
+
+The `stdlib_sparse` module provides several derived types defining known sparse matrix data structures. It also provides basic sparse kernels such as sparse matrix vector and conversion between matrix types.
+
+## Derived types provided
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### The `sparse_type` abstract derived type
+#### Status
+
+Experimental
+
+#### Description
+The `sparse_type` is defined as an abstract derived type holding the basic common meta data needed to define a sparse matrix. All other sparse types falvors are derived from the `sparse_type`.
+
+```Fortran
+type, public, abstract :: sparse_type
+    integer :: nrows   !> number of rows
+    integer :: ncols   !> number of columns
+    integer :: nnz     !> number of non-zero values
+    integer :: storage !> assumed storage symmetry
+    integer :: base    !> index base = 0 for (C) or 1 (Fortran)
+end type
+```
+
+The storage integer label should be assigned from the module's internal enumerator containing the following three enums:
+
+```Fortran
+enum, bind(C)
+    enumerator :: sparse_full  !> Full Sparse matrix (no symmetry considerations)
+    enumerator :: sparse_lower !> Symmetric Sparse matrix with triangular inferior storage
+    enumerator :: sparse_upper !> Symmetric Sparse matrix with triangular supperior storage
+end enum
+```
+In the following, all sparse kinds will be presented in two main flavors: a data-less type `<matrix>_type` useful for topological graph operations. And real/complex valued types `<matrix>_<kind>` containing the `data` buffer for the matrix values. The following rectangular matrix will be used to showcase how each sparse matrix holds the data internally:
+
+$$ M = \begin{bmatrix} 
+    9 & 0 & 0  & 0 & -3 \\
+    4 & 7 & 0  & 0 & 0 \\
+    0 & 8 & -1 & 8 & 0 \\
+    4 & 0 & 5  & 6 & 0 \\
+  \end{bmatrix} $$
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `COO`: The COOrdinates compressed sparse format
+#### Status
+
+Experimental
+
+#### Description
+The `COO`, triplet or `ijv` format defines all non-zero elements of the matrix by explicitly allocating the `i,j` index and the value of the matrix. While some implementations use separate `row` and `col` arrays for the index, here we use a 2D array in order to promote fast memory acces to `ij`.
+
+```Fortran
+type(COO_sp) :: COO
+call COO%malloc(4,5,10)
+COO%data(:)   = real([9,-3,4,7,8,-1,8,4,5,6])
+COO%index(1:2,1)  = [1,1]
+COO%index(1:2,2)  = [1,5]
+COO%index(1:2,3)  = [2,1]
+COO%index(1:2,4)  = [2,2]
+COO%index(1:2,5)  = [3,2]
+COO%index(1:2,6)  = [3,3]
+COO%index(1:2,7)  = [3,4]
+COO%index(1:2,8)  = [4,1]
+COO%index(1:2,9)  = [4,3]
+COO%index(1:2,10) = [4,4]
+```
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `CSR`: The Compressed Sparse Row or Yale format
+#### Status
+
+Experimental
+
+#### Description
+The Compressed Sparse Row or Yale format `CSR` stores the matrix index by compressing the row indeces with a counter pointer `rowptr` enabling to know the first and last non-zero colum index `col` of the given row. 
+
+```Fortran
+type(CSR_sp) :: CSR
+call CSR%malloc(4,5,10)
+CSR%data(:)   = real([9,-3,4,7,8,-1,8,4,5,6])
+CSR%col(:)    = [1,5,1,2,2,3,4,1,3,4]
+CSR%rowptr(:) = [1,3,5,8,11]
+```
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `CSC`: The Compressed Sparse Column format
+#### Status
+
+Experimental
+
+#### Description
+The Compressed Sparse Colum `CSC` is similar to the `CSR` format but values are accesed first by colum, thus an index counter is given by `colptr` which enables accessing the start and ending rows of a given colum in the `row` index table. 
+
+```Fortran
+type(CSC_sp) :: CSC
+call CSC%malloc(4,5,10)
+CSC%data(:)   = real([9,4,4,7,8,-1,5,8,6,-3])
+CSC%row(:)    = [1,2,4,2,3,3,4,3,4,1]
+CSC%colptr(:) = [1,4,6,8,10,11]
+```
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `ELLPACK`: ELL-pack storage format
+#### Status
+
+Experimental
+
+#### Description
+The `ELL` format stores the data in a dense matrix of $nrows \times K$ in column major order. By imposing a constant number of zeros per row $K$, this format will incure in additional zeros being stored, but it enables efficient vectorization as memory acces are carried out by constant sized strides. 
+
+```Fortran
+type(ELL_sp) :: ELL
+call ELL%malloc(num_rows=4,num_cols=5,num_nz_row=3)
+ELL%data(1,1:3)   = real([9,-3,0])
+ELL%data(2,1:3)   = real([4,7,0])
+ELL%data(3,1:3)   = real([8,-1,8])
+ELL%data(4,1:3)   = real([4,5,6])
+
+ELL%index(1,1:3) = [1,5,0]
+ELL%index(2,1:3) = [1,2,0]
+ELL%index(3,1:3) = [2,3,4]
+ELL%index(4,1:3) = [1,3,4]
+```
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `SELL-C`: The Sliced ELLPACK with Constant blocks format
+#### Status
+
+Experimental
+
+#### Description
+The Sliced ELLPACK format `SELLC` is a variation of the `ELLPACK` format. This modification reduces the storage size compared to the `ELLPACK` format but maintaining its efficient data access scheme. It can be seen as an intermediate format between `CSR` and `ELLPACK`. For more details read [here](https://arxiv.org/pdf/1307.6209v1)
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+## `spmv` - Sparse Matrix-Vector product
+
+### Status
+
+Experimental
+
+### Description
+
+Provide sparse matrix-vector product kernels for the current supported sparse matrix types.
+
+$$y=\alpha*M*x+\beta*y$$
+
+### Syntax
+
+`call ` [[stdlib_sparse_spmv(module):spmv(interface)]] `(matrix,vec_x,vec_y [,alpha,beta])`
+
+### Arguments
+
+`matrix`, `intent(in)`: Shall be a `real` or `complex` sparse type matrix.
+
+`vec_x`, `intent(in)`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array.
+
+`vec_y`, `intent(inout)`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array.
+
+`alpha`, `intent(in)`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`.
+
+`beta`, `intent(in)`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`.
+
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+## `sparse_conversion` - Sparse matrix to matrix conversions
+
+### Status
+
+Experimental
+
+### Description
+
+This module provides facility functions for converting between storage formats.
+
+### Syntax
+
+`call ` [[stdlib_sparse_conversion(module):coo2ordered(interface)]] `(coo)`
+
+### Arguments
+
+`COO`, `intent(inout)`: Shall be any `COO` type. The same object will be returned with the arrays reallocated to the correct size after removing duplicates.
+
+`sort_data`, `logical(in)`, `optional`:: Boolean optional to determine whether to sort the data in the COO graph while sorting the index array, defult `.false.`.
+
+### Syntax
+
+`call ` [[stdlib_sparse_conversion(module):dense2coo(interface)]] `(dense,coo)`
+
+### Arguments
+
+`dense`, `intent(in)`: Shall be a rank-2 array of `real` or `complex` type.
+
+`coo`, `intent(inout)`: Shall be a `COO` type of `real` or `complex` type.
+
+### Syntax
+
+`call ` [[stdlib_sparse_conversion(module):coo2dense(interface)]] `(coo,dense)`
+
+### Arguments
+
+`coo`, `intent(in)`: Shall be a `COO` type of `real` or `complex` type.
+
+`dense`, `intent(inout)`: Shall be a rank-2 array of `real` or `complex` type.
+
+### Syntax
+
+`call ` [[stdlib_sparse_conversion(module):coo2csr(interface)]] `(coo,csr)`
+
+### Arguments
+
+`coo`, `intent(in)`: Shall be a `COO` type of `real` or `complex` type.
+
+`csr`, `intent(inout)`: Shall be a `CSR` type of `real` or `complex` type.
+
+### Syntax
+
+`call ` [[stdlib_sparse_conversion(module):csr2coo(interface)]] `(csr,coo)`
+
+### Arguments
+
+`csr`, `intent(in)`: Shall be a `CSR` type of `real` or `complex` type.
+
+`coo`, `intent(inout)`: Shall be a `COO` type of `real` or `complex` type.
+
+### Syntax
+
+`call ` [[stdlib_sparse_conversion(module):csr2sellc(interface)]] `(csr,sellc[,chunk])`
+
+### Arguments
+
+`csr`, `intent(in)`: Shall be a `CSR` type of `real` or `complex` type.
+
+`sellc`, `intent(inout)`: Shall be a `SELLC` type of `real` or `complex` type.
+
+`chunk`, `intent(in)`, `optional`: chunk size for the Sliced ELLPACK format.
+
+## Example
+```fortran
+{!example/linalg/example_sparse_spmv.f90!}
+```

example/linalg/CMakeLists.txt

@@ -23,6 +23,7 @@ ADD_EXAMPLE(lstsq2)
 ADD_EXAMPLE(solve1)
 ADD_EXAMPLE(solve2)
 ADD_EXAMPLE(solve3)
+ADD_EXAMPLE(sparse_spmv)
 ADD_EXAMPLE(svd)
 ADD_EXAMPLE(svdvals)
 ADD_EXAMPLE(determinant)

example/linalg/example_sparse_spmv.f90

@@ -0,0 +1,36 @@
+program example_sparse_spmv
+    use stdlib_linalg_constants, only: dp
+    use stdlib_sparse
+    implicit none
+
+    integer, parameter :: m = 4, n = 2
+    real(dp) :: A(m,n), x(n)
+    real(dp) :: y_dense(m), y_coo(m), y_csr(m)
+    real(dp) :: alpha, beta
+    type(COO_dp) :: COO
+    type(CSR_dp) :: CSR
+
+    call random_number(A)
+    ! Convert from dense to COO and CSR matrices
+    call dense2coo( A , COO )
+    call coo2csr( COO , CSR )
+
+    ! Initialize vectors
+    x       = 1._dp
+    y_dense = 2._dp
+    y_coo   = y_dense
+    y_csr   = y_dense
+
+    ! Perform matrix-vector product
+    alpha = 3._dp; beta = 2._dp
+    y_dense = alpha * matmul(A,x) + beta * y_dense
+    call spmv( COO , x , y_coo , alpha = alpha, beta = beta )
+    call spmv( CSR , x , y_csr , alpha = alpha, beta = beta )
+
+    print *, 'dense :', y_dense
+    print *, 'coo   :', y_coo
+    print *, 'csr   :', y_csr
+
+  end program example_sparse_spmv
+
+

include/common.fypp

@@ -38,6 +38,7 @@
 
 #! Real types to be considered during templating
 #:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS]
+#:set REAL_SUFFIX = REAL_KINDS
 
 #! Collected (kind, type) tuples for real types
 #:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_INIT))
@@ -62,6 +63,7 @@
 
 #! Complex types to be considered during templating
 #:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS]
+#:set CMPLX_SUFFIX = ["c{}".format(k) for k in CMPLX_KINDS]
 
 #! Collected (kind, type, initial) tuples for complex types
 #:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_INIT))
@@ -102,6 +104,9 @@
 #! Bitset types to be considered during templating
 #:set BITSET_TYPES = ["type({})".format(k) for k in BITSET_KINDS]
 
+#! Sparse types to be considered during templating
+#:set SPARSE_KINDS = ["COO", "CSR", "CSC", "ELL"]
+
 #! Whether Fortran 90 compatible code should be generated
 #:set VERSION90 = defined('VERSION90')
 

src/CMakeLists.txt

@@ -37,6 +37,9 @@ set(fppFiles
     stdlib_sorting_ord_sort.fypp
     stdlib_sorting_sort.fypp
     stdlib_sorting_sort_index.fypp
+    stdlib_sparse_conversion.fypp
+    stdlib_sparse_kinds.fypp
+    stdlib_sparse_spmv.fypp
     stdlib_specialfunctions_gamma.fypp
     stdlib_stats.fypp
     stdlib_stats_corr.fypp
@@ -109,6 +112,7 @@ set(SRC
     stdlib_logger.f90
     stdlib_sorting_radix_sort.f90
     stdlib_system.F90
+    stdlib_sparse.f90
     stdlib_specialfunctions.f90
     stdlib_specialfunctions_legendre.f90
     stdlib_quadrature_gauss.f90

src/stdlib_sparse.f90

@@ -0,0 +1,6 @@
+!! public API
+module stdlib_sparse
+    use stdlib_sparse_kinds
+    use stdlib_sparse_spmv
+    use stdlib_sparse_conversion
+end module stdlib_sparse