From 7e3b0a9cec7f8ef12292cecfe1829350daf428ee Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 2 May 2026 20:50:49 +0200 Subject: [PATCH 1/8] Low-level interfaces for COO --- src/sparse/stdlib_sparse_spmv.fypp | 12 ++++++++++++ src/sparse/stdlib_sparse_spmv_coo.fypp | 19 ++++++++++++++++--- 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 5524ef7c0..1445a7dc7 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -32,6 +32,18 @@ module stdlib_sparse_spmv character(1), intent(in), optional :: op end subroutine + module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, storage, nnz, vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: index(:,:) + integer, intent(in) :: storage + integer(ilp), intent(in) :: nnz + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + module subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op) type(CSR_${s1}$_type), intent(in) :: matrix ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 7aba6e9d8..0d009299b 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -20,6 +20,21 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op + + call spmv_coo_sub_${rank}$d_${s1}$(matrix%data, matrix%index, matrix%storage, matrix%nnz, vec_x,vec_y,alpha,beta,op) + + end subroutine + + module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, storage, nnz, vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: index(:,:) !! Matrix coordinates index(2,nnz) + integer, intent(in) :: storage !! storage + integer(ilp), intent(in) :: nnz !! number of non-zero values + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op ${t1}$ :: alpha_ character(1) :: op_ integer(ilp) :: col_index, k, row_index @@ -33,7 +48,6 @@ contains vec_y = zero_${s1}$ endif - associate( data => matrix%data, index => matrix%index, storage => matrix%storage, nnz => matrix%nnz ) select case(op_) case(sparse_op_none) if(storage == sparse_full) then @@ -92,10 +106,9 @@ contains end if #:endif end select - end associate end subroutine #:endfor #:endfor -end submodule stdlib_sparse_spmv_coo \ No newline at end of file +end submodule stdlib_sparse_spmv_coo From 51780555855940ade51bd5425d5f0f6d82d2c06a Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 2 May 2026 20:54:44 +0200 Subject: [PATCH 2/8] COO: low-level interfaces --- src/sparse/stdlib_sparse_spmv.fypp | 37 ++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 1445a7dc7..43b2522e2 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -16,6 +16,7 @@ module stdlib_sparse_spmv implicit none private + !! Version experimental !! !! Apply the sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ @@ -31,18 +32,6 @@ module stdlib_sparse_spmv ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op end subroutine - - module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, storage, nnz, vec_x,vec_y,alpha,beta,op) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: index(:,:) - integer, intent(in) :: storage - integer(ilp), intent(in) :: nnz - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ - ${t1}$, intent(in), optional :: alpha - ${t1}$, intent(in), optional :: beta - character(1), intent(in), optional :: op - end subroutine module subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op) type(CSR_${s1}$_type), intent(in) :: matrix @@ -82,6 +71,30 @@ module stdlib_sparse_spmv end subroutine #:endfor end interface + + !! Version experimental + !! + !! Apply the COO sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! [Specifications](../page/specs/stdlib_sparse.html#spmv) + interface spmv_coo + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, storage, nnz, vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: index(:,:) + integer, intent(in) :: storage + integer(ilp), intent(in) :: nnz + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + #:endfor + #:endfor + end interface + public :: spmv + public :: spmv_coo end module From 71519b9950bff58557c430406e4014751def764e Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 2 May 2026 21:09:14 +0200 Subject: [PATCH 3/8] COO: test for low-level spmv --- test/linalg/test_linalg_sparse.fypp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 0ebe235b5..35b16eac3 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -63,15 +63,28 @@ contains call check(error, all(vec_y1 == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return + ! High-level API call spmv( COO, vec_x, vec_y2 ) call check(error, all(vec_y1 == vec_y2) ) if (allocated(error)) return + ! Low-level API + call spmv_coo( COO%data ,COO%index, COO%storage, COO%nnz, vec_x, vec_y2 ) + call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' ) + if (allocated(error)) return + ! Test in-place transpose vec_y1 = 1._wp + ! High-level API call spmv( COO, vec_y1, vec_x, op=sparse_op_transpose ) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return + + ! Low-level API + call spmv_coo( COO%data ,COO%index, COO%storage, COO%nnz, vec_y1, vec_x, op=sparse_op_transpose ) + call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'COO - transpose - low level API: wrong output' ) + if (allocated(error)) return + end block #:endfor end subroutine @@ -715,4 +728,4 @@ program tester write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" error stop end if -end program \ No newline at end of file +end program From 0ccedc844c0bac4f13a4c169759783cccb2078d2 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 2 May 2026 21:23:28 +0200 Subject: [PATCH 4/8] cleaning --- src/sparse/stdlib_sparse_spmv.fypp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 43b2522e2..aeb5113a6 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -16,7 +16,6 @@ module stdlib_sparse_spmv implicit none private - !! Version experimental !! !! Apply the sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ @@ -32,7 +31,7 @@ module stdlib_sparse_spmv ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op end subroutine - + module subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op) type(CSR_${s1}$_type), intent(in) :: matrix ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ From 2be029fb84291761d1bbee54a2d023ec4f33d43b Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Mon, 4 May 2026 19:38:08 +0200 Subject: [PATCH 5/8] CRS: low-level API --- src/sparse/stdlib_sparse_spmv.fypp | 28 +++++++++++++++++++++++++- src/sparse/stdlib_sparse_spmv_coo.fypp | 9 +++++---- src/sparse/stdlib_sparse_spmv_csr.fypp | 27 +++++++++++++++++++------ test/linalg/test_linalg_sparse.fypp | 22 ++++++++++++++++---- 4 files changed, 71 insertions(+), 15 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index aeb5113a6..cea11dd0c 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -78,7 +78,7 @@ module stdlib_sparse_spmv interface spmv_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, storage, nnz, vec_x,vec_y,alpha,beta,op) + module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, nnz, storage, vec_x,vec_y,alpha,beta,op) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) integer, intent(in) :: storage @@ -93,7 +93,33 @@ module stdlib_sparse_spmv #:endfor end interface + !! Version experimental + !! + !! Apply the CRS sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! [Specifications](../page/specs/stdlib_sparse.html#spmv) + interface spmv_csr + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + module subroutine spmv_csr_sub_${rank}$d_${s1}$(data,col,rowptr,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: col(:) !! matrix column pointer + integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer + integer(ilp), intent(in) :: nnz !! number of non-zero values + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage !! storage + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + #:endfor + #:endfor + end interface + public :: spmv public :: spmv_coo + public :: spmv_csr end module diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 0d009299b..a7a362b08 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -21,15 +21,16 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_coo_sub_${rank}$d_${s1}$(matrix%data, matrix%index, matrix%storage, matrix%nnz, vec_x,vec_y,alpha,beta,op) + call spmv_coo_sub_${rank}$d_${s1}$(matrix%data, matrix%index, matrix%nnz, matrix%storage, & + vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, storage, nnz, vec_x,vec_y,alpha,beta,op) - ${t1}$, intent(in) :: data(:) + module subroutine spmv_coo_sub_${rank}$d_${s1}$(data,index,nnz,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) !! Matrix coordinates index(2,nnz) - integer, intent(in) :: storage !! storage integer(ilp), intent(in) :: nnz !! number of non-zero values + integer, intent(in) :: storage !! storage ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in), optional :: alpha diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index 929449969..f1e4bee22 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -20,6 +20,25 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op + + call spmv_csr_sub_${rank}$d_${s1}$(matrix%data, matrix%col, matrix%rowptr, matrix%nnz, matrix%nrows, matrix%ncols, matrix%storage, & + vec_x,vec_y,alpha,beta,op) + + end subroutine + + module subroutine spmv_csr_sub_${rank}$d_${s1}$(data,col,rowptr,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: col(:) !! matrix column pointer + integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer + integer(ilp), intent(in) :: nnz !! number of non-zero values + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage !! storage + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op ${t1}$ :: alpha_ character(1) :: op_ integer(ilp) :: i, j @@ -37,9 +56,6 @@ contains else vec_y = zero_${s1}$ endif - - associate( data => matrix%data, col => matrix%col, rowptr => matrix%rowptr, & - & nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage ) if( storage == sparse_full .and. op_==sparse_op_none ) then do i = 1, nrows @@ -114,10 +130,9 @@ contains end do #:endif end if - end associate end subroutine - + #:endfor #:endfor -end submodule stdlib_sparse_spmv_csr \ No newline at end of file +end submodule stdlib_sparse_spmv_csr diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 35b16eac3..fe191d0e6 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -69,7 +69,7 @@ contains if (allocated(error)) return ! Low-level API - call spmv_coo( COO%data ,COO%index, COO%storage, COO%nnz, vec_x, vec_y2 ) + call spmv_coo( COO%data ,COO%index, COO%nnz, COO%storage, vec_x, vec_y2 ) call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' ) if (allocated(error)) return @@ -81,7 +81,7 @@ contains if (allocated(error)) return ! Low-level API - call spmv_coo( COO%data ,COO%index, COO%storage, COO%nnz, vec_y1, vec_x, op=sparse_op_transpose ) + call spmv_coo( COO%data ,COO%index, COO%nnz, COO%storage, vec_y1, vec_x, op=sparse_op_transpose ) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'COO - transpose - low level API: wrong output' ) if (allocated(error)) return @@ -134,15 +134,29 @@ contains allocate( vec_x(5) , source = 1._wp ) allocate( vec_y(4) , source = 0._wp ) + ! High-level API call spmv( CSR, vec_x, vec_y ) - call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) + call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSR - high-level API - wrong output' ) + if (allocated(error)) return + + ! Low-level API + call spmv_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nnz, CSR%nrows, CSR%ncols, CSR%storage, & + vec_x, vec_y ) + + call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSR - low-level API - wrong output' ) if (allocated(error)) return ! Test in-place transpose vec_y = 1._wp + ! High-level API call spmv( CSR, vec_y, vec_x, op=sparse_op_transpose ) - call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) + call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'CSR transpose - high-level API - wrong output' ) + if (allocated(error)) return + + call spmv_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nnz, CSR%nrows, CSR%ncols, CSR%storage, & + vec_y, vec_x, op=sparse_op_transpose ) + call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'CSR transpose - low-level API - wrong output' ) if (allocated(error)) return end block #:endfor From f0b85fa5080cb1fae5b63ad4279890d50e19cbd0 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Mon, 4 May 2026 19:53:26 +0200 Subject: [PATCH 6/8] CSC: low-level API --- src/sparse/stdlib_sparse_spmv.fypp | 28 +++++++++++++++++++++++++- src/sparse/stdlib_sparse_spmv_csc.fypp | 25 ++++++++++++++++++----- test/linalg/test_linalg_sparse.fypp | 15 ++++++++++++++ 3 files changed, 62 insertions(+), 6 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index cea11dd0c..a68b47028 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -95,7 +95,32 @@ module stdlib_sparse_spmv !! Version experimental !! - !! Apply the CRS sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! Apply the CSC sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! [Specifications](../page/specs/stdlib_sparse.html#spmv) + interface spmv_csc + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + module subroutine spmv_csc_sub_${rank}$d_${s1}$(data,colptr,row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: colptr(:) !! matrix column pointer + integer(ilp), intent(in) :: row(:) !! matrix row pointer + integer(ilp), intent(in) :: nnz !! number of non-zero values + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage !! storage + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + #:endfor + #:endfor + end interface + + !! Version experimental + !! + !! Apply the CSR sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_csr #:for k1, t1, s1 in (KINDS_TYPES) @@ -120,6 +145,7 @@ module stdlib_sparse_spmv public :: spmv public :: spmv_coo + public :: spmv_csc public :: spmv_csr end module diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 83b4ed41d..95a173c10 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -20,6 +20,25 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op + + call spmv_csc_sub_${rank}$d_${s1}$(matrix%data,matrix%colptr,matrix%row,matrix%nnz,matrix%nrows,matrix%ncols,matrix%storage, & + vec_x,vec_y,alpha,beta,op) + + end subroutine + + module subroutine spmv_csc_sub_${rank}$d_${s1}$(data,colptr,row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: colptr(:) !! matrix column pointer + integer(ilp), intent(in) :: row(:) !! matrix row pointer + integer(ilp), intent(in) :: nnz !! number of non-zero values + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage !! storage + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op ${t1}$ :: alpha_ character(1) :: op_ integer(ilp) :: i, j @@ -38,8 +57,6 @@ contains vec_y = zero_${s1}$ endif - associate( data => matrix%data, colptr => matrix%colptr, row => matrix%row, & - & nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage ) if( storage == sparse_full .and. op_==sparse_op_none ) then do concurrent(j=1:ncols) aux = alpha_ * vec_x(${rksfx2(rank-1)}$j) @@ -110,10 +127,8 @@ contains end do #:endif end if - end associate end subroutine - #:endfor #:endfor -end submodule stdlib_sparse_spmv_csc \ No newline at end of file +end submodule stdlib_sparse_spmv_csc diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index fe191d0e6..b0e907c69 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -179,16 +179,31 @@ contains allocate( vec_x(5) , source = 1._wp ) allocate( vec_y(4) , source = 0._wp ) + !High-level API call spmv( CSC, vec_x, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return + !High-level API + call spmv_csc( CSC%data, CSC%colptr, CSC%row, CSC%nnz, CSC%nrows, CSC%ncols, CSC%storage, vec_x, vec_y ) + + call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) + if (allocated(error)) return + ! Test in-place transpose vec_y = 1._wp + !High-level API call spmv( CSC, vec_y, vec_x, op=sparse_op_transpose ) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return + + !Low-level API + call spmv_csc( CSC%data, CSC%colptr, CSC%row, CSC%nnz, CSC%nrows, CSC%ncols, CSC%storage, vec_y, vec_x, & + op=sparse_op_transpose ) + call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) + if (allocated(error)) return + end block #:endfor end subroutine From e4c796abde211d6833b86167575296fabae1d3ab Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Mon, 4 May 2026 20:38:45 +0200 Subject: [PATCH 7/8] ELL: low-level API --- src/sparse/stdlib_sparse_spmv.fypp | 26 ++++++++++++++++++++++++++ src/sparse/stdlib_sparse_spmv_ell.fypp | 26 +++++++++++++++++++++----- test/linalg/test_linalg_sparse.fypp | 15 +++++++++++++++ 3 files changed, 62 insertions(+), 5 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index a68b47028..f7a36e9cb 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -143,9 +143,35 @@ module stdlib_sparse_spmv #:endfor end interface + !! Version experimental + !! + !! Apply the ELL sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! [Specifications](../page/specs/stdlib_sparse.html#spmv) + interface spmv_ell + #:for k1, t1, s1 in (KINDS_TYPES) + #:for rank in RANKS + module subroutine spmv_ell_sub_${rank}$d_${s1}$(data,index,mnz_p_row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:,:) + integer(ilp), intent(in) :: index(:,:) + integer(ilp), intent(in) :: mnz_p_row + integer(ilp), intent(in) :: nnz + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + #:endfor + #:endfor + end interface + public :: spmv public :: spmv_coo public :: spmv_csc public :: spmv_csr + public :: spmv_ell end module diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index 243d1cbd4..68e8b32d6 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -20,6 +20,25 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op + + call spmv_ell_sub_${rank}$d_${s1}$(matrix%data,matrix%index,matrix%K, & + matrix%nnz,matrix%nrows,matrix%ncols,matrix%storage, & + vec_x,vec_y,alpha,beta,op) + end subroutine + + module subroutine spmv_ell_sub_${rank}$d_${s1}$(data,index,mnz_p_row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + ${t1}$, intent(in) :: data(:,:) + integer(ilp), intent(in) :: index(:,:) + integer, intent(in) :: mnz_p_row + integer(ilp), intent(in) :: nnz + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage + ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op ${t1}$ :: alpha_ character(1) :: op_ integer(ilp) :: i, j, k @@ -32,8 +51,6 @@ contains else vec_y = zero_${s1}$ endif - associate( data => matrix%data, index => matrix%index, MNZ_P_ROW => matrix%K, & - & nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage ) if( storage == sparse_full .and. op_==sparse_op_none ) then do i = 1, nrows do k = 1, MNZ_P_ROW @@ -78,10 +95,9 @@ contains end do #:endif end if - end associate end subroutine - + #:endfor #:endfor -end submodule stdlib_sparse_spmv_ell \ No newline at end of file +end submodule stdlib_sparse_spmv_ell diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index b0e907c69..77fadeb6d 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -231,16 +231,31 @@ contains allocate( vec_x(5) , source = 1._wp ) allocate( vec_y(4) , source = 0._wp ) + !High-level API call spmv( ELL, vec_x, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return + !Low-level API + call spmv_ell( ELL%data, ELL%index, ELL%K, & + ELL%nnz, ELL%nrows, ELL%ncols, ELL%storage, vec_x, vec_y ) + + call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) + if (allocated(error)) return + ! Test in-place transpose vec_y = 1._wp + !High-level API call spmv( ELL, vec_y, vec_x, op=sparse_op_transpose ) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return + + !Low-level API + call spmv_ell( ELL%data, ELL%index, ELL%K, & + ELL%nnz, ELL%nrows, ELL%ncols, ELL%storage, vec_y, vec_x, op=sparse_op_transpose ) + call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) + if (allocated(error)) return end block #:endfor From 7f4309aaca2fff74a3f9d4012084a25d2f699c38 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Mon, 4 May 2026 21:01:42 +0200 Subject: [PATCH 8/8] SELLC: low-level API --- src/sparse/stdlib_sparse_spmv.fypp | 26 +++++++++++++++++++++ src/sparse/stdlib_sparse_spmv_sellc.fypp | 29 +++++++++++++++++++----- test/linalg/test_linalg_sparse.fypp | 18 ++++++++++++++- 3 files changed, 66 insertions(+), 7 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index f7a36e9cb..6f408f2b0 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -168,10 +168,36 @@ module stdlib_sparse_spmv #:endfor end interface + !! Version experimental + !! + !! Apply the SELLC sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ + !! [Specifications](../page/specs/stdlib_sparse.html#spmv) + interface spmv_sellc + #:for k1, t1, s1 in (KINDS_TYPES) + module subroutine spmv_sellc_sub_${s1}$(data,ia,ja,cs,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves + ${t1}$, intent(in) :: data(:,:) + integer(ilp), intent(in) :: ia(:) + integer(ilp), intent(in) :: ja(:,:) + integer, intent(in) :: cs + integer(ilp), intent(in) :: nnz + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage + ${t1}$, intent(in) :: vec_x(:) + ${t1}$, intent(inout) :: vec_y(:) + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op + end subroutine + #:endfor + end interface + public :: spmv public :: spmv_coo public :: spmv_csc public :: spmv_csr public :: spmv_ell + public :: spmv_sellc end module diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 52205debd..c4982e3fd 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -16,6 +16,28 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op + + call spmv_sellc_sub_${s1}$(matrix%data, matrix%rowptr, matrix%col, matrix%chunk_size, & + matrix%nnz, matrix%nrows, matrix%ncols, matrix%storage, & + vec_x,vec_y,alpha,beta,op) + + end subroutine + + module subroutine spmv_sellc_sub_${s1}$(data,ia,ja,cs,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves + ${t1}$, intent(in) :: data(:,:) + integer(ilp), intent(in) :: ia(:) + integer(ilp), intent(in) :: ja(:,:) + integer, intent(in) :: cs + integer(ilp), intent(in) :: nnz + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols + integer, intent(in) :: storage + ${t1}$, intent(in) :: vec_x(:) + ${t1}$, intent(inout) :: vec_y(:) + ${t1}$, intent(in), optional :: alpha + ${t1}$, intent(in), optional :: beta + character(1), intent(in), optional :: op ${t1}$ :: alpha_ character(1) :: op_ integer(ilp) :: i, nz, rowidx, num_chunks, rm @@ -29,9 +51,6 @@ contains vec_y = zero_${s1}$ endif - associate( data => matrix%data, ia => matrix%rowptr , ja => matrix%col, cs => matrix%chunk_size, & - & nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage ) - if( .not.any( ${CHUNKS}$ == cs ) ) then print *, "error: sellc chunk size not supported." return @@ -107,7 +126,6 @@ contains print *, "error: sellc format for spmv operation not yet supported." return end if - end associate contains #:for chunk in CHUNKS @@ -187,7 +205,6 @@ contains #:endif end subroutine - #:endfor -end submodule stdlib_sparse_spmv_sellc \ No newline at end of file +end submodule stdlib_sparse_spmv_sellc diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 77fadeb6d..396bae95e 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -289,17 +289,33 @@ contains allocate( vec_x(6) , source = 1._wp ) allocate( vec_y(6) , source = 0._wp ) + !High-level API call spmv( SELLC, vec_x, vec_y ) call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)) ) if (allocated(error)) return + !Low-level API + call spmv_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & + SELLC%nnz, SELLC%nrows, SELLC%ncols, SELLC%storage, vec_x, vec_y ) + + call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)) ) + if (allocated(error)) return + ! Test in-place transpose vec_x = real( [1,2,3,4,5,6] , kind=wp ) call spmv( CSR, vec_x, vec_y , op = sparse_op_transpose ) allocate( vec_y2(6) , source = 0._wp ) + !High-level API call spmv( SELLC, vec_x, vec_y2 , op = sparse_op_transpose ) - + + call check(error, all(vec_y == vec_y2)) + if (allocated(error)) return + + !Low-level API + call spmv_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & + SELLC%nnz, SELLC%nrows, SELLC%ncols, SELLC%storage, vec_x, vec_y2 , op = sparse_op_transpose ) + call check(error, all(vec_y == vec_y2)) if (allocated(error)) return