From 7e3b0a9cec7f8ef12292cecfe1829350daf428ee Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 2 May 2026 20:50:49 +0200 Subject: [PATCH 01/27] 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 02/27] 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 03/27] 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 04/27] 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 05/27] 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 06/27] 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 07/27] 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 08/27] 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 From 2aa51174d055a9e2011559b2372b84805687a137 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 24 Jul 2026 16:32:05 +0200 Subject: [PATCH 09/27] Rename specific spmv by spmv_kernel --- src/sparse/stdlib_sparse_spmv.fypp | 30 ++++++++++++------------ src/sparse/stdlib_sparse_spmv_coo.fypp | 4 ++-- src/sparse/stdlib_sparse_spmv_csc.fypp | 4 ++-- src/sparse/stdlib_sparse_spmv_csr.fypp | 4 ++-- src/sparse/stdlib_sparse_spmv_ell.fypp | 4 ++-- src/sparse/stdlib_sparse_spmv_sellc.fypp | 4 ++-- test/linalg/test_linalg_sparse.fypp | 20 ++++++++-------- 7 files changed, 35 insertions(+), 35 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 6f408f2b0..5c89327a1 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -75,10 +75,10 @@ module stdlib_sparse_spmv !! !! Apply the COO sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ !! [Specifications](../page/specs/stdlib_sparse.html#spmv) - interface spmv_coo + interface spmv_kernel_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_coo_sub_${rank}$d_${s1}$(data, index, nnz, storage, vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_coo_${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 @@ -97,10 +97,10 @@ module stdlib_sparse_spmv !! !! Apply the CSC sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ !! [Specifications](../page/specs/stdlib_sparse.html#spmv) - interface spmv_csc + interface spmv_kernel_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) + module subroutine spmv_kernel_csc_${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 @@ -122,10 +122,10 @@ module stdlib_sparse_spmv !! !! Apply the CSR sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ !! [Specifications](../page/specs/stdlib_sparse.html#spmv) - interface spmv_csr + interface spmv_kernel_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) + module subroutine spmv_kernel_csr_${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 @@ -147,10 +147,10 @@ module stdlib_sparse_spmv !! !! Apply the ELL sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ !! [Specifications](../page/specs/stdlib_sparse.html#spmv) - interface spmv_ell + interface spmv_kernel_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) + module subroutine spmv_kernel_ell_${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 @@ -172,9 +172,9 @@ module stdlib_sparse_spmv !! !! Apply the SELLC sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$ !! [Specifications](../page/specs/stdlib_sparse.html#spmv) - interface spmv_sellc + interface spmv_kernel_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) + module subroutine spmv_kernel_sellc_${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(:) @@ -194,10 +194,10 @@ module stdlib_sparse_spmv end interface public :: spmv - public :: spmv_coo - public :: spmv_csc - public :: spmv_csr - public :: spmv_ell - public :: spmv_sellc + public :: spmv_kernel_coo + public :: spmv_kernel_csc + public :: spmv_kernel_csr + public :: spmv_kernel_ell + public :: spmv_kernel_sellc end module diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index a7a362b08..31e3b10bd 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -21,12 +21,12 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_coo_sub_${rank}$d_${s1}$(matrix%data, matrix%index, matrix%nnz, matrix%storage, & + call spmv_kernel_coo_${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,nnz,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_coo_${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(ilp), intent(in) :: nnz !! number of non-zero values diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 95a173c10..8100fd571 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -21,12 +21,12 @@ contains ${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, & + call spmv_kernel_csc_${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) + module subroutine spmv_kernel_csc_${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 diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index f1e4bee22..8e1101533 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -21,12 +21,12 @@ contains ${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, & + call spmv_kernel_csr_${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) + module subroutine spmv_kernel_csr_${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 diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index 68e8b32d6..b492860da 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -21,12 +21,12 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_ell_sub_${rank}$d_${s1}$(matrix%data,matrix%index,matrix%K, & + call spmv_kernel_ell_${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) + module subroutine spmv_kernel_ell_${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 diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index c4982e3fd..714c63c4e 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -17,13 +17,13 @@ contains ${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, & + call spmv_kernel_sellc_${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) + module subroutine spmv_kernel_sellc_${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(:) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 396bae95e..e267b053e 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%nnz, COO%storage, vec_x, vec_y2 ) + call spmv_kernel_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%nnz, COO%storage, vec_y1, vec_x, op=sparse_op_transpose ) + call spmv_kernel_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 @@ -141,7 +141,7 @@ contains if (allocated(error)) return ! Low-level API - call spmv_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nnz, CSR%nrows, CSR%ncols, CSR%storage, & + call spmv_kernel_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' ) @@ -154,7 +154,7 @@ contains 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, & + call spmv_kernel_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 @@ -186,7 +186,7 @@ contains 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 spmv_kernel_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 @@ -199,7 +199,7 @@ contains 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, & + call spmv_kernel_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 @@ -238,7 +238,7 @@ contains if (allocated(error)) return !Low-level API - call spmv_ell( ELL%data, ELL%index, ELL%K, & + call spmv_kernel_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)) ) @@ -252,7 +252,7 @@ contains if (allocated(error)) return !Low-level API - call spmv_ell( ELL%data, ELL%index, ELL%K, & + call spmv_kernel_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 @@ -296,7 +296,7 @@ contains if (allocated(error)) return !Low-level API - call spmv_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & + call spmv_kernel_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)) ) @@ -313,7 +313,7 @@ contains if (allocated(error)) return !Low-level API - call spmv_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & + call spmv_kernel_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)) From e7c195d6b315d49776d855d621323f8971717aaf Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 24 Jul 2026 21:06:17 +0200 Subject: [PATCH 10/27] remove nnz API --- doc/specs/stdlib_sparse.md | 105 ++++++++++++++++++++++- src/sparse/stdlib_sparse_spmv.fypp | 12 +-- src/sparse/stdlib_sparse_spmv_csc.fypp | 5 +- src/sparse/stdlib_sparse_spmv_csr.fypp | 5 +- src/sparse/stdlib_sparse_spmv_ell.fypp | 5 +- src/sparse/stdlib_sparse_spmv_sellc.fypp | 5 +- test/linalg/test_linalg_sparse.fypp | 16 ++-- 7 files changed, 124 insertions(+), 29 deletions(-) diff --git a/doc/specs/stdlib_sparse.md b/doc/specs/stdlib_sparse.md index 7ce53008f..fb0f2fc5e 100644 --- a/doc/specs/stdlib_sparse.md +++ b/doc/specs/stdlib_sparse.md @@ -217,6 +217,109 @@ $$y=\alpha*op(M)*x+\beta*y$$ `op`, `optional`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian` + +## `spmv_kernel` - Non-object-oriented sparse matrix-vector product + +### Status + +Experimental + +### Description + +Provide non-object-oriented sparse matrix-vector product kernels for the current supported sparse matrix types. + +$$y=\alpha*op(M)*x+\beta*y$$ + +### Syntax + +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_coo(interface)]] `(data,index,nnz,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(data,colptr,row,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(data,col,rowptr,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(data,index,mnz_p_row,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(data,ia,ja,cs,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` + +### Arguments + +For the `COO` format: + +`data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument. + +`index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. + +`nnz`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. integer(ilp), intent(in) :: nnz + +For the `CSC` format: + +`data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument. + +`colptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix column pointer + +`row`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix row pointer + +`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. !! number of non-zero values + +`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +For the `CSR` format: + +`data`: Shall be a rank-1 array of `real` or `complex` type array. It is an `intent(in)` argument. + +`col`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix column pointer + +`rowptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix row pointer + +`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. !! number of non-zero values + +`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +For the `ELL` format: + +`data`: Shall be a rank-2 array of `real` or `complex` type. It is an `intent(in)` argument. + +`index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. + +`mnz_p_row`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +For the `SELLC` format: + +`data`: Shall be a rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. + +`ia`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. + +`ja`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. + +`cs`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. + +`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +For the all formats: + +`storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. + +`vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. + +`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument. + +`alpha`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument. + +`beta`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument. + +`op`, `optional`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian` + ## Sparse matrix to matrix conversions @@ -411,4 +514,4 @@ The definition of all standard arithmetic operators have been overloaded to be a `B = A / alpha` or `B = alpha / A` -*Note*: scalar addition and subtraction operators perform element-wise operations only on the stored (non-zero) values, not on the full mathematical matrix. Meaning, the sparsity pattern is preserved. \ No newline at end of file +*Note*: scalar addition and subtraction operators perform element-wise operations only on the stored (non-zero) values, not on the full mathematical matrix. Meaning, the sparsity pattern is preserved. diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 5c89327a1..62e146fa7 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -100,11 +100,10 @@ module stdlib_sparse_spmv interface spmv_kernel_csc #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(data,colptr,row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(data,colptr,row,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 @@ -125,11 +124,10 @@ module stdlib_sparse_spmv interface spmv_kernel_csr #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(data,col,rowptr,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(data,col,rowptr,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 @@ -150,11 +148,10 @@ module stdlib_sparse_spmv interface spmv_kernel_ell #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(data,index,mnz_p_row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(data,index,mnz_p_row,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 @@ -174,13 +171,12 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(data,ia,ja,cs,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(data,ia,ja,cs,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 diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 8100fd571..e050c3ac0 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -21,16 +21,15 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_csc_${rank}$d_${s1}$(matrix%data,matrix%colptr,matrix%row,matrix%nnz,matrix%nrows,matrix%ncols,matrix%storage, & + call spmv_kernel_csc_${rank}$d_${s1}$(matrix%data,matrix%colptr,matrix%row,matrix%nrows,matrix%ncols,matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(data,colptr,row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(data,colptr,row,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 diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index 8e1101533..4acfadf70 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -21,16 +21,15 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_csr_${rank}$d_${s1}$(matrix%data, matrix%col, matrix%rowptr, matrix%nnz, matrix%nrows, matrix%ncols, matrix%storage, & + call spmv_kernel_csr_${rank}$d_${s1}$(matrix%data, matrix%col, matrix%rowptr, matrix%nrows, matrix%ncols, matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(data,col,rowptr,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(data,col,rowptr,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 diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index b492860da..8860df6ee 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -22,15 +22,14 @@ contains character(1), intent(in), optional :: op call spmv_kernel_ell_${rank}$d_${s1}$(matrix%data,matrix%index,matrix%K, & - matrix%nnz,matrix%nrows,matrix%ncols,matrix%storage, & + matrix%nrows,matrix%ncols,matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(data,index,mnz_p_row,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(data,index,mnz_p_row,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 diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 714c63c4e..134e420d3 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -18,18 +18,17 @@ contains character(1), intent(in), optional :: op call spmv_kernel_sellc_${s1}$(matrix%data, matrix%rowptr, matrix%col, matrix%chunk_size, & - matrix%nnz, matrix%nrows, matrix%ncols, matrix%storage, & + matrix%nrows, matrix%ncols, matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(data,ia,ja,cs,nnz,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(data,ia,ja,cs,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 diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index e267b053e..ef422832a 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -141,7 +141,7 @@ contains if (allocated(error)) return ! Low-level API - call spmv_kernel_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nnz, CSR%nrows, CSR%ncols, CSR%storage, & + call spmv_kernel_csr( CSR%data, CSR%col, CSR%rowptr, 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' ) @@ -154,7 +154,7 @@ contains 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_kernel_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nnz, CSR%nrows, CSR%ncols, CSR%storage, & + call spmv_kernel_csr( CSR%data, CSR%col, CSR%rowptr, 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 @@ -186,7 +186,7 @@ contains if (allocated(error)) return !High-level API - call spmv_kernel_csc( CSC%data, CSC%colptr, CSC%row, CSC%nnz, CSC%nrows, CSC%ncols, CSC%storage, vec_x, vec_y ) + call spmv_kernel_csc( CSC%data, CSC%colptr, CSC%row, 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 @@ -199,7 +199,7 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_csc( CSC%data, CSC%colptr, CSC%row, CSC%nnz, CSC%nrows, CSC%ncols, CSC%storage, vec_y, vec_x, & + call spmv_kernel_csc( CSC%data, CSC%colptr, CSC%row, 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 @@ -239,7 +239,7 @@ contains !Low-level API call spmv_kernel_ell( ELL%data, ELL%index, ELL%K, & - ELL%nnz, ELL%nrows, ELL%ncols, ELL%storage, vec_x, vec_y ) + 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 @@ -253,7 +253,7 @@ contains !Low-level API call spmv_kernel_ell( ELL%data, ELL%index, ELL%K, & - ELL%nnz, ELL%nrows, ELL%ncols, ELL%storage, vec_y, vec_x, op=sparse_op_transpose ) + 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 @@ -297,7 +297,7 @@ contains !Low-level API call spmv_kernel_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & - SELLC%nnz, SELLC%nrows, SELLC%ncols, SELLC%storage, vec_x, vec_y ) + 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 @@ -314,7 +314,7 @@ contains !Low-level API call spmv_kernel_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 ) + 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 From bb4ba13199a9e8a8a2895924a4d6d0b08bb95d4e Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 24 Jul 2026 21:36:04 +0200 Subject: [PATCH 11/27] mv nrows and ncols in the api of spmv csc csr --- src/sparse/stdlib_sparse_spmv.fypp | 12 ++++++------ src/sparse/stdlib_sparse_spmv_csc.fypp | 8 ++++---- src/sparse/stdlib_sparse_spmv_csr.fypp | 8 ++++---- test/linalg/test_linalg_sparse.fypp | 8 ++++---- 4 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 62e146fa7..2414ca2e9 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -100,12 +100,12 @@ module stdlib_sparse_spmv interface spmv_kernel_csc #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(data,colptr,row,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(nrows,ncols,data,colptr,row,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: colptr(:) !! matrix column pointer integer(ilp), intent(in) :: row(:) !! matrix row pointer - 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)}$ @@ -124,12 +124,12 @@ module stdlib_sparse_spmv interface spmv_kernel_csr #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(data,col,rowptr,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: col(:) !! matrix column pointer integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer - 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)}$ diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index e050c3ac0..22db7d1a6 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -21,17 +21,17 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_csc_${rank}$d_${s1}$(matrix%data,matrix%colptr,matrix%row,matrix%nrows,matrix%ncols,matrix%storage, & + call spmv_kernel_csc_${rank}$d_${s1}$(matrix%nrows,matrix%ncols,matrix%data,matrix%colptr,matrix%row,matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(data,colptr,row,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(nrows,ncols,data,colptr,row,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: colptr(:) !! matrix column pointer integer(ilp), intent(in) :: row(:) !! matrix row pointer - 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)}$ diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index 4acfadf70..677ccb0e0 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -21,17 +21,17 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_csr_${rank}$d_${s1}$(matrix%data, matrix%col, matrix%rowptr, matrix%nrows, matrix%ncols, matrix%storage, & + call spmv_kernel_csr_${rank}$d_${s1}$(matrix%nrows, matrix%ncols, matrix%data, matrix%col, matrix%rowptr, matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(data,col,rowptr,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: col(:) !! matrix column pointer integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer - 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)}$ diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index ef422832a..2d74cbafe 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -141,7 +141,7 @@ contains if (allocated(error)) return ! Low-level API - call spmv_kernel_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nrows, CSR%ncols, CSR%storage, & + call spmv_kernel_csr(CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, 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' ) @@ -154,7 +154,7 @@ contains 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_kernel_csr( CSR%data, CSR%col, CSR%rowptr, CSR%nrows, CSR%ncols, CSR%storage, & + call spmv_kernel_csr(CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, 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 @@ -186,7 +186,7 @@ contains if (allocated(error)) return !High-level API - call spmv_kernel_csc( CSC%data, CSC%colptr, CSC%row, CSC%nrows, CSC%ncols, CSC%storage, vec_x, vec_y ) + call spmv_kernel_csc(CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, vec_x, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return @@ -199,7 +199,7 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_csc( CSC%data, CSC%colptr, CSC%row, CSC%nrows, CSC%ncols, CSC%storage, vec_y, vec_x, & + call spmv_kernel_csc(CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, 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 From e3e31155cf34b458f1bd6ff48d2009cc8367ec9c Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 24 Jul 2026 21:52:19 +0200 Subject: [PATCH 12/27] spmv_kernel: clean api --- doc/specs/stdlib_sparse.md | 64 +++++++++--------------- src/sparse/stdlib_sparse_spmv.fypp | 16 +++--- src/sparse/stdlib_sparse_spmv_coo.fypp | 6 ++- src/sparse/stdlib_sparse_spmv_ell.fypp | 10 ++-- src/sparse/stdlib_sparse_spmv_sellc.fypp | 10 ++-- test/linalg/test_linalg_sparse.fypp | 21 ++++---- 6 files changed, 59 insertions(+), 68 deletions(-) diff --git a/doc/specs/stdlib_sparse.md b/doc/specs/stdlib_sparse.md index fb0f2fc5e..6d42e0a6e 100644 --- a/doc/specs/stdlib_sparse.md +++ b/doc/specs/stdlib_sparse.md @@ -232,21 +232,40 @@ $$y=\alpha*op(M)*x+\beta*y$$ ### Syntax -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_coo(interface)]] `(data,index,nnz,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(data,colptr,row,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(data,col,rowptr,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(data,index,mnz_p_row,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(data,ia,ja,cs,nnz,nrows,ncols,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_coo(interface)]] `(nrows,ncols,data,index,nnz,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(nrows,ncols,data,colptr,row,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(nrows,ncols,data,ia,ja,cs,storage,vec_x,vec_y [,alpha,beta,op])` ### Arguments +For the all formats: + +`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. + +`storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. + +`vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. + +`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument. + +`alpha`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument. + +`beta`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument. + +`op`, `optional`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian` + For the `COO` format: `data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument. `index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`nnz`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. integer(ilp), intent(in) :: nnz +`nnz`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. !! number of non-zero values + For the `CSC` format: @@ -256,11 +275,6 @@ For the `CSC` format: `row`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix row pointer -`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. !! number of non-zero values - -`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. For the `CSR` format: @@ -270,11 +284,6 @@ For the `CSR` format: `rowptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix row pointer -`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. !! number of non-zero values - -`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. For the `ELL` format: @@ -284,11 +293,6 @@ For the `ELL` format: `mnz_p_row`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. -`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. For the `SELLC` format: @@ -300,25 +304,7 @@ For the `SELLC` format: `cs`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. -`nnz`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. - -For the all formats: -`storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. - -`vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. - -`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument. - -`alpha`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument. - -`beta`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument. - -`op`, `optional`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian` ## Sparse matrix to matrix conversions diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 2414ca2e9..27976bec5 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -78,7 +78,9 @@ module stdlib_sparse_spmv interface spmv_kernel_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(data, index, nnz, storage, vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(nrows, ncols, data, index, nnz, storage, vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) integer, intent(in) :: storage @@ -148,12 +150,12 @@ module stdlib_sparse_spmv interface spmv_kernel_ell #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(data,index,mnz_p_row,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: index(:,:) integer(ilp), intent(in) :: mnz_p_row - 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)}$ @@ -171,14 +173,14 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(data,ia,ja,cs,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,cs,storage,vec_x,vec_y,alpha,beta,op) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) integer(ilp), intent(in) :: ja(:,:) integer, intent(in) :: cs - integer(ilp), intent(in) :: nrows - integer(ilp), intent(in) :: ncols integer, intent(in) :: storage ${t1}$, intent(in) :: vec_x(:) ${t1}$, intent(inout) :: vec_y(:) diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 31e3b10bd..45682e3f9 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -21,12 +21,14 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_coo_${rank}$d_${s1}$(matrix%data, matrix%index, matrix%nnz, matrix%storage, & + call spmv_kernel_coo_${rank}$d_${s1}$(matrix%nrows, matrix%ncols, matrix%data, matrix%index, matrix%nnz, matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(data,index,nnz,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(nrows,ncols,data,index,nnz,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) !! Matrix coordinates index(2,nnz) integer(ilp), intent(in) :: nnz !! number of non-zero values diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index 8860df6ee..f1f988897 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -21,17 +21,17 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_ell_${rank}$d_${s1}$(matrix%data,matrix%index,matrix%K, & - matrix%nrows,matrix%ncols,matrix%storage, & + call spmv_kernel_ell_${rank}$d_${s1}$(matrix%nrows,matrix%ncols,matrix%data,matrix%index,matrix%K, & + matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(data,index,mnz_p_row,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y,alpha,beta,op) + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: index(:,:) integer, intent(in) :: mnz_p_row - 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)}$ diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 134e420d3..7b0d4b88b 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -17,20 +17,20 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_sellc_${s1}$(matrix%data, matrix%rowptr, matrix%col, matrix%chunk_size, & - matrix%nrows, matrix%ncols, matrix%storage, & + call spmv_kernel_sellc_${s1}$(matrix%nrows, matrix%ncols, matrix%data, matrix%rowptr, matrix%col, & + matrix%chunk_size, matrix%storage, & vec_x,vec_y,alpha,beta,op) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(data,ia,ja,cs,nrows,ncols,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,cs,storage,vec_x,vec_y,alpha,beta,op) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves + integer(ilp), intent(in) :: nrows + integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) integer(ilp), intent(in) :: ja(:,:) integer, intent(in) :: cs - integer(ilp), intent(in) :: nrows - integer(ilp), intent(in) :: ncols integer, intent(in) :: storage ${t1}$, intent(in) :: vec_x(:) ${t1}$, intent(inout) :: vec_y(:) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 2d74cbafe..1e174ce15 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_kernel_coo( COO%data ,COO%index, COO%nnz, COO%storage, vec_x, vec_y2 ) + call spmv_kernel_coo(COO%nrows, COO%ncols, 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,8 @@ contains if (allocated(error)) return ! Low-level API - call spmv_kernel_coo( COO%data ,COO%index, COO%nnz, COO%storage, vec_y1, vec_x, op=sparse_op_transpose ) + call spmv_kernel_coo(COO%nrows, COO%ncols, 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 @@ -238,8 +239,8 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_ell( ELL%data, ELL%index, ELL%K, & - ELL%nrows, ELL%ncols, ELL%storage, vec_x, vec_y ) + call spmv_kernel_ell(ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, & + ELL%storage, vec_x, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return @@ -252,8 +253,8 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_ell( ELL%data, ELL%index, ELL%K, & - ELL%nrows, ELL%ncols, ELL%storage, vec_y, vec_x, op=sparse_op_transpose ) + call spmv_kernel_ell(ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, & + 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 @@ -296,8 +297,8 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & - SELLC%nrows, SELLC%ncols, SELLC%storage, vec_x, vec_y ) + call spmv_kernel_sellc(SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & + 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 @@ -313,8 +314,8 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_sellc( SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & - SELLC%nrows, SELLC%ncols, SELLC%storage, vec_x, vec_y2 , op = sparse_op_transpose ) + call spmv_kernel_sellc(SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & + SELLC%storage, vec_x, vec_y2 , op = sparse_op_transpose ) call check(error, all(vec_y == vec_y2)) if (allocated(error)) return From efcfa9d44368d2b7744a56b898f599ea4cc6e79c Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Sat, 25 Jul 2026 08:28:32 +0200 Subject: [PATCH 13/27] Update specs --- doc/specs/stdlib_sparse.md | 25 ++++++++++----------- src/sparse/stdlib_sparse_spmv.fypp | 4 ++-- src/sparse/stdlib_sparse_spmv_sellc.fypp | 28 ++++++++++++------------ 3 files changed, 28 insertions(+), 29 deletions(-) diff --git a/doc/specs/stdlib_sparse.md b/doc/specs/stdlib_sparse.md index 6d42e0a6e..bb749d068 100644 --- a/doc/specs/stdlib_sparse.md +++ b/doc/specs/stdlib_sparse.md @@ -94,7 +94,7 @@ CSR%rowptr(:) = [1,3,5,8,11] Experimental #### Description -The Compressed Sparse Colum `CSC` is similar to the `CSR` format but values are accesed first by column, thus an index counter is given by `colptr` which enables to know the first and last non-zero row index of a given colum. +The Compressed Sparse Colum `CSC` is similar to the `CSR` format but values are accessed first by column, thus an index counter is given by `colptr` which enables to know the first and last non-zero row index of a given colum. ```Fortran type(CSC_sp_type) :: CSC @@ -236,17 +236,17 @@ $$y=\alpha*op(M)*x+\beta*y$$ `call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(nrows,ncols,data,colptr,row,storage,vec_x,vec_y [,alpha,beta,op])` `call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y [,alpha,beta,op])` `call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(nrows,ncols,data,ia,ja,cs,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,vec_y [,alpha,beta,op])` ### Arguments For the all formats: -`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. +`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the number of rows in the sparse matrix. -`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. +`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the number of columns in the sparse matrix. -`storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. +`storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. It defines the symmetry storage of the sparse matrix and must be one of `sparse_full`, `sparse_lower`, or `sparse_upper`. `vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. @@ -264,25 +264,25 @@ For the `COO` format: `index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`nnz`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. !! number of non-zero values +`nnz`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. It specifies the number of non-zero values in `index`. For the `CSC` format: `data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument. -`colptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix column pointer +`colptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. -`row`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix row pointer +`row`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. For the `CSR` format: `data`: Shall be a rank-1 array of `real` or `complex` type array. It is an `intent(in)` argument. -`col`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix column pointer +`col`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. -`rowptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. !! matrix row pointer +`rowptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. For the `ELL` format: @@ -291,7 +291,7 @@ For the `ELL` format: `index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`mnz_p_row`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. +`mnz_p_row`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the constant number of elements per row $K$. For the `SELLC` format: @@ -302,8 +302,7 @@ For the `SELLC` format: `ja`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`cs`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. - +`chunk_size`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 27976bec5..a315de3dc 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -173,14 +173,14 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,cs,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,vec_y,alpha,beta,op) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) integer(ilp), intent(in) :: ja(:,:) - integer, intent(in) :: cs + integer, intent(in) :: chunk_size integer, intent(in) :: storage ${t1}$, intent(in) :: vec_x(:) ${t1}$, intent(inout) :: vec_y(:) diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 7b0d4b88b..35fad6e39 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -23,14 +23,14 @@ contains end subroutine - module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,cs,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,vec_y,alpha,beta,op) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) integer(ilp), intent(in) :: ja(:,:) - integer, intent(in) :: cs + integer, intent(in) :: chunk_size integer, intent(in) :: storage ${t1}$, intent(in) :: vec_x(:) ${t1}$, intent(inout) :: vec_y(:) @@ -50,16 +50,16 @@ contains vec_y = zero_${s1}$ endif - if( .not.any( ${CHUNKS}$ == cs ) ) then + if( .not.any( ${CHUNKS}$ == chunk_size ) ) then print *, "error: sellc chunk size not supported." return end if - num_chunks = nrows / cs - rm = nrows - num_chunks * cs + num_chunks = nrows / chunk_size + rm = nrows - num_chunks * chunk_size if( storage == sparse_full .and. op_==sparse_op_none ) then - select case(cs) + select case(chunk_size) #:for chunk in CHUNKS case(${chunk}$) do i = 1, num_chunks @@ -74,13 +74,13 @@ contains if(rm>0)then i = num_chunks + 1 nz = ia(i+1) - ia(i) - rowidx = (i - 1)*cs + 1 - call chunk_kernel_rm(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:)) + rowidx = (i - 1)*chunk_size + 1 + call chunk_kernel_rm(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:)) end if else if( storage == sparse_full .and. op_==sparse_op_transpose ) then - select case(cs) + select case(chunk_size) #:for chunk in CHUNKS case(${chunk}$) do i = 1, num_chunks @@ -95,14 +95,14 @@ contains if(rm>0)then i = num_chunks + 1 nz = ia(i+1) - ia(i) - rowidx = (i - 1)*cs + 1 - call chunk_kernel_rm_trans(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y) + rowidx = (i - 1)*chunk_size + 1 + call chunk_kernel_rm_trans(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y) end if #:if t1.startswith('complex') else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then - select case(cs) + select case(chunk_size) #:for chunk in CHUNKS case(${chunk}$) do i = 1, num_chunks @@ -117,8 +117,8 @@ contains if(rm>0)then i = num_chunks + 1 nz = ia(i+1) - ia(i) - rowidx = (i - 1)*cs + 1 - call chunk_kernel_rm_herm(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y) + rowidx = (i - 1)*chunk_size + 1 + call chunk_kernel_rm_herm(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y) end if #:endif else From 9bef87843469d2b7cb03ca8a505931a2f20837fd Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 09:57:05 +0200 Subject: [PATCH 14/27] Set op, alpha and beta non-optional --- src/sparse/stdlib_sparse_spmv.fypp | 40 ++++++------- src/sparse/stdlib_sparse_spmv_coo.fypp | 55 +++++++++-------- src/sparse/stdlib_sparse_spmv_csc.fypp | 76 +++++++++++++----------- src/sparse/stdlib_sparse_spmv_csr.fypp | 73 ++++++++++++----------- src/sparse/stdlib_sparse_spmv_ell.fypp | 61 ++++++++++--------- src/sparse/stdlib_sparse_spmv_sellc.fypp | 55 +++++++++-------- test/linalg/test_linalg_sparse.fypp | 58 ++++++++++++------ 7 files changed, 231 insertions(+), 187 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index a315de3dc..8371f8058 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_kernel_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(nrows, ncols, data, index, nnz, storage, vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,nnz,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) @@ -87,9 +87,9 @@ module stdlib_sparse_spmv 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 + ${t1}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op end subroutine #:endfor #:endfor @@ -102,7 +102,7 @@ module stdlib_sparse_spmv interface spmv_kernel_csc #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(nrows,ncols,data,colptr,row,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,colptr,row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) @@ -111,9 +111,9 @@ module stdlib_sparse_spmv 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}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op end subroutine #:endfor #:endfor @@ -126,7 +126,7 @@ module stdlib_sparse_spmv interface spmv_kernel_csr #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) @@ -135,9 +135,9 @@ module stdlib_sparse_spmv 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}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op end subroutine #:endfor #:endfor @@ -150,7 +150,7 @@ module stdlib_sparse_spmv interface spmv_kernel_ell #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) @@ -159,9 +159,9 @@ module stdlib_sparse_spmv 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}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op end subroutine #:endfor #:endfor @@ -173,7 +173,7 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols @@ -184,9 +184,9 @@ module stdlib_sparse_spmv 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}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op end subroutine #:endfor end interface diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 45682e3f9..ac4f2613f 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -21,12 +21,24 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_coo_${rank}$d_${s1}$(matrix%nrows, matrix%ncols, matrix%data, matrix%index, matrix%nnz, matrix%storage, & - vec_x,vec_y,alpha,beta,op) + ${t1}$ :: alpha_, beta_ + character(1) :: op_ + + op_ = sparse_op_none; if(present(op)) op_ = op + + alpha_ = one_${k1}$ + if(present(alpha)) alpha_ = alpha + + beta_ = zero_${s1}$ + if(present(beta)) beta_ = beta + + call spmv_kernel_coo_${rank}$d_${s1}$(op_, alpha_, & + matrix%nrows, matrix%ncols, matrix%data, matrix%index, matrix%nnz, matrix%storage, & + vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(nrows,ncols,data,index,nnz,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,nnz,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) @@ -35,38 +47,29 @@ contains 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_ + ${t1}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op integer(ilp) :: col_index, k, row_index - op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ - if(present(alpha)) alpha_ = alpha - if(present(beta)) then - vec_y = beta * vec_y - else - vec_y = zero_${s1}$ - endif + vec_y = beta * vec_y - select case(op_) + select case(op) case(sparse_op_none) if(storage == sparse_full) then do k = 1, nnz row_index = index(1,k) col_index = index(2,k) - vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index) + vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha*data(k) * vec_x(${rksfx2(rank-1)}$col_index) end do else do k = 1, nnz row_index = index(1,k) col_index = index(2,k) - vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index) + vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha*data(k) * vec_x(${rksfx2(rank-1)}$col_index) if( row_index==col_index ) cycle - vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$row_index) + vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha*data(k) * vec_x(${rksfx2(rank-1)}$row_index) end do end if @@ -75,16 +78,16 @@ contains do k = 1, nnz col_index = index(1,k) row_index = index(2,k) - vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index) + vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha*data(k) * vec_x(${rksfx2(rank-1)}$col_index) end do else do k = 1, nnz col_index = index(1,k) row_index = index(2,k) - vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index) + vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha*data(k) * vec_x(${rksfx2(rank-1)}$col_index) if( row_index==col_index ) cycle - vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$row_index) + vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha*data(k) * vec_x(${rksfx2(rank-1)}$row_index) end do end if @@ -94,16 +97,16 @@ contains do k = 1, nnz col_index = index(1,k) row_index = index(2,k) - vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$col_index) + vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$col_index) end do else do k = 1, nnz col_index = index(1,k) row_index = index(2,k) - vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$col_index) + vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$col_index) if( row_index==col_index ) cycle - vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha_*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$row_index) + vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$row_index) end do end if diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 22db7d1a6..6a4f83ae2 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -21,12 +21,29 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_csc_${rank}$d_${s1}$(matrix%nrows,matrix%ncols,matrix%data,matrix%colptr,matrix%row,matrix%storage, & - vec_x,vec_y,alpha,beta,op) + ${t1}$ :: alpha_, beta_ + character(1) :: op_ + + op_ = sparse_op_none; if(present(op)) op_ = op + + alpha_ = one_${k1}$ + if(present(alpha)) alpha_ = alpha + + beta_ = zero_${k1}$ + if(present(beta)) beta_ = beta + + if(present(beta)) then + else + vec_y = zero_${s1}$ + endif + + call spmv_kernel_csc_${rank}$d_${s1}$(op_,alpha_, & + matrix%nrows,matrix%ncols,matrix%data,matrix%colptr,matrix%row,matrix%storage, & + vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(nrows,ncols,data,colptr,row,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,colptr,row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) @@ -35,11 +52,9 @@ contains 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_ + ${t1}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op integer(ilp) :: i, j #:if rank == 1 ${t1}$ :: aux @@ -47,82 +62,75 @@ contains ${t1}$ :: aux(size(vec_x,dim=1)) #:endif - op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ - if(present(alpha)) alpha_ = alpha - if(present(beta)) then - vec_y = beta * vec_y - else - vec_y = zero_${s1}$ - endif + vec_y = beta * vec_y - if( storage == sparse_full .and. op_==sparse_op_none ) then + if( storage == sparse_full .and. op==sparse_op_none ) then do concurrent(j=1:ncols) - aux = alpha_ * vec_x(${rksfx2(rank-1)}$j) + aux = alpha * vec_x(${rksfx2(rank-1)}$j) do i = colptr(j), colptr(j+1)-1 vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + data(i) * aux end do end do - else if( storage == sparse_full .and. op_==sparse_op_transpose ) then + else if( storage == sparse_full .and. op==sparse_op_transpose ) then do concurrent(j=1:ncols) aux = zero_${k1}$ do i = colptr(j), colptr(j+1)-1 aux = aux + data(i) * vec_x(${rksfx2(rank-1)}$row(i)) end do - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha * aux end do - else if( storage == sparse_lower .and. op_/=sparse_op_hermitian )then + else if( storage == sparse_lower .and. op/=sparse_op_hermitian )then do j = 1 , ncols aux = vec_x(${rksfx2(rank-1)}$j) * data(colptr(j)) do i = colptr(j)+1, colptr(j+1)-1 aux = aux + data(i) * vec_x(${rksfx2(rank-1)}$row(i)) - vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * data(i) * vec_x(${rksfx2(rank-1)}$j) + vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha * data(i) * vec_x(${rksfx2(rank-1)}$j) end do - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha * aux end do - else if( storage == sparse_upper .and. op_/=sparse_op_hermitian )then + else if( storage == sparse_upper .and. op/=sparse_op_hermitian )then do j = 1 , ncols aux = zero_${s1}$ do i = colptr(j), colptr(j+1)-2 aux = aux + data(i) * vec_x(${rksfx2(rank-1)}$row(i)) - vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * data(i) * vec_x(${rksfx2(rank-1)}$j) + vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha * data(i) * vec_x(${rksfx2(rank-1)}$j) end do aux = aux + data(colptr(j+1)-1) * vec_x(${rksfx2(rank-1)}$j) - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha * aux end do #:if t1.startswith('complex') - else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then + else if( storage == sparse_full .and. op==sparse_op_hermitian ) then do concurrent(j=1:ncols) aux = zero_${k1}$ do i = colptr(j), colptr(j+1)-1 aux = aux + conjg(data(i)) * vec_x(${rksfx2(rank-1)}$row(i)) end do - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha * aux end do - else if( storage == sparse_lower .and. op_==sparse_op_hermitian )then + else if( storage == sparse_lower .and. op==sparse_op_hermitian )then do j = 1 , ncols aux = vec_x(${rksfx2(rank-1)}$j) * conjg(data(colptr(j))) do i = colptr(j)+1, colptr(j+1)-1 aux = aux + conjg(data(i)) * vec_x(${rksfx2(rank-1)}$row(i)) - vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j) + vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha * conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j) end do - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha * aux end do - else if( storage == sparse_upper .and. op_==sparse_op_hermitian )then + else if( storage == sparse_upper .and. op==sparse_op_hermitian )then do j = 1 , ncols aux = zero_${s1}$ do i = colptr(j), colptr(j+1)-2 aux = aux + conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j) - vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j) + vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha * conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j) end do aux = aux + conjg(data(colptr(j+1)-1)) * vec_x(${rksfx2(rank-1)}$j) - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha * aux end do #:endif end if diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index 677ccb0e0..b6a93a0cb 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -21,12 +21,24 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_csr_${rank}$d_${s1}$(matrix%nrows, matrix%ncols, matrix%data, matrix%col, matrix%rowptr, matrix%storage, & - vec_x,vec_y,alpha,beta,op) + ${t1}$ :: alpha_, beta_ + character(1) :: op_ + + op_ = sparse_op_none; if(present(op)) op_ = op + + alpha_ = one_${k1}$ + if(present(alpha)) alpha_ = alpha + + beta_ = zero_${k1}$ + if(present(beta)) beta_ = beta + + call spmv_kernel_csr_${rank}$d_${s1}$(op_, alpha_, & + matrix%nrows, matrix%ncols, matrix%data, matrix%col, matrix%rowptr, matrix%storage, & + vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) @@ -35,11 +47,9 @@ contains 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_ + ${t1}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op integer(ilp) :: i, j #:if rank == 1 ${t1}$ :: aux, aux2 @@ -47,85 +57,78 @@ contains ${t1}$ :: aux(size(vec_x,dim=1)), aux2(size(vec_x,dim=1)) #:endif - op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ - if(present(alpha)) alpha_ = alpha - if(present(beta)) then - vec_y = beta * vec_y - else - vec_y = zero_${s1}$ - endif - - if( storage == sparse_full .and. op_==sparse_op_none ) then + vec_y = beta * vec_y + + if( storage == sparse_full .and. op==sparse_op_none ) then do i = 1, nrows aux = zero_${k1}$ do j = rowptr(i), rowptr(i+1)-1 aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j)) end do - vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha * aux end do - else if( storage == sparse_full .and. op_==sparse_op_transpose ) then + else if( storage == sparse_full .and. op==sparse_op_transpose ) then do i = 1, nrows - aux = alpha_ * vec_x(${rksfx2(rank-1)}$i) + aux = alpha * vec_x(${rksfx2(rank-1)}$i) do j = rowptr(i), rowptr(i+1)-1 vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux end do end do - else if( storage == sparse_lower .and. op_/=sparse_op_hermitian )then + else if( storage == sparse_lower .and. op/=sparse_op_hermitian )then do i = 1 , nrows aux = zero_${s1}$ - aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i) + aux2 = alpha * vec_x(${rksfx2(rank-1)}$i) do j = rowptr(i), rowptr(i+1)-2 aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j)) vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux2 end do - aux = alpha_ * aux + data(j) * aux2 + aux = alpha * aux + data(j) * aux2 vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + aux end do - else if( storage == sparse_upper .and. op_/=sparse_op_hermitian )then + else if( storage == sparse_upper .and. op/=sparse_op_hermitian )then do i = 1 , nrows aux = vec_x(${rksfx2(rank-1)}$i) * data(rowptr(i)) - aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i) + aux2 = alpha * vec_x(${rksfx2(rank-1)}$i) do j = rowptr(i)+1, rowptr(i+1)-1 aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j)) vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux2 end do - vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha * aux end do #:if t1.startswith('complex') - else if( storage == sparse_full .and. op_==sparse_op_hermitian) then + else if( storage == sparse_full .and. op==sparse_op_hermitian) then do i = 1, nrows - aux = alpha_ * vec_x(${rksfx2(rank-1)}$i) + aux = alpha * vec_x(${rksfx2(rank-1)}$i) do j = rowptr(i), rowptr(i+1)-1 vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux end do end do - else if( storage == sparse_lower .and. op_==sparse_op_hermitian )then + else if( storage == sparse_lower .and. op==sparse_op_hermitian )then do i = 1 , nrows aux = zero_${s1}$ - aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i) + aux2 = alpha * vec_x(${rksfx2(rank-1)}$i) do j = rowptr(i), rowptr(i+1)-2 aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j)) vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2 end do - aux = alpha_ * aux + conjg(data(j)) * aux2 + aux = alpha * aux + conjg(data(j)) * aux2 vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + aux end do - else if( storage == sparse_upper .and. op_==sparse_op_hermitian )then + else if( storage == sparse_upper .and. op==sparse_op_hermitian )then do i = 1 , nrows aux = vec_x(${rksfx2(rank-1)}$i) * conjg(data(rowptr(i))) - aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i) + aux2 = alpha * vec_x(${rksfx2(rank-1)}$i) do j = rowptr(i)+1, rowptr(i+1)-1 aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j)) vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2 end do - vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux + vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha * aux end do #:endif end if diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index f1f988897..a6e3d7fa1 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -21,12 +21,25 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_ell_${rank}$d_${s1}$(matrix%nrows,matrix%ncols,matrix%data,matrix%index,matrix%K, & + ${t1}$ :: alpha_, beta_ + character(1) :: op_ + + op_ = sparse_op_none; if(present(op)) op_ = op + + alpha_ = one_${k1}$ + if(present(alpha)) alpha_ = alpha + + beta_ = zero_${k1}$ + if(present(beta)) beta_ = alpha + + call spmv_kernel_ell_${rank}$d_${s1}$(op_,alpha_, & + matrix%nrows,matrix%ncols,matrix%data,matrix%index,matrix%K, & matrix%storage, & - vec_x,vec_y,alpha,beta,op) + vec_x,beta_,vec_y) + end subroutine - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) @@ -35,61 +48,53 @@ contains 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_ + ${t1}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op integer(ilp) :: i, j, k - op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ - if(present(alpha)) alpha_ = alpha - if(present(beta)) then - vec_y = beta * vec_y - else - vec_y = zero_${s1}$ - endif - if( storage == sparse_full .and. op_==sparse_op_none ) then + vec_y = beta * vec_y + + if( storage == sparse_full .and. op==sparse_op_none ) then do i = 1, nrows do k = 1, MNZ_P_ROW j = index(i,k) - if(j>0) vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$j) + if(j>0) vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$j) end do end do - else if( storage == sparse_full .and. op_==sparse_op_transpose ) then + else if( storage == sparse_full .and. op==sparse_op_transpose ) then do i = 1, nrows do k = 1, MNZ_P_ROW j = index(i,k) - if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$i) + if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$i) end do end do - else if( storage /= sparse_full .and. op_/=sparse_op_hermitian ) then + else if( storage /= sparse_full .and. op/=sparse_op_hermitian ) then do i = 1, nrows do k = 1, MNZ_P_ROW j = index(i,k) if(j<=0) cycle - vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$j) + vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$j) if(i==j) cycle - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$i) + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$i) end do end do #:if t1.startswith('complex') - else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then + else if( storage == sparse_full .and. op==sparse_op_hermitian ) then do i = 1, nrows do k = 1, MNZ_P_ROW j = index(i,k) - if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i) + if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i) end do end do - else if( storage /= sparse_full .and. op_==sparse_op_hermitian ) then + else if( storage /= sparse_full .and. op==sparse_op_hermitian ) then do i = 1, nrows do k = 1, MNZ_P_ROW j = index(i,k) if(j<=0) cycle - vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$j) + vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$j) if(i==j) cycle - vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i) + vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i) end do end do #:endif diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 35fad6e39..a1b830b09 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -17,13 +17,25 @@ contains ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - call spmv_kernel_sellc_${s1}$(matrix%nrows, matrix%ncols, matrix%data, matrix%rowptr, matrix%col, & + ${t1}$ :: alpha_, beta_ + character(1) :: op_ + + op_ = sparse_op_none; if(present(op)) op_ = op + + alpha_ = one_${s1}$ + if(present(alpha)) alpha_ = alpha + + beta_ = zero_${s1}$ + if(present(beta)) beta_ = beta + + call spmv_kernel_sellc_${s1}$(op_,alpha_, & + matrix%nrows, matrix%ncols, matrix%data, matrix%rowptr, matrix%col, & matrix%chunk_size, matrix%storage, & - vec_x,vec_y,alpha,beta,op) + vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,vec_y,alpha,beta,op) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols @@ -34,30 +46,21 @@ contains 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_ + ${t1}$, intent(in) :: alpha + ${t1}$, intent(in) :: beta + character(1), intent(in) :: op integer(ilp) :: i, nz, rowidx, num_chunks, rm - op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${s1}$ - if(present(alpha)) alpha_ = alpha - if(present(beta)) then - vec_y = beta * vec_y - else - vec_y = zero_${s1}$ - endif - if( .not.any( ${CHUNKS}$ == chunk_size ) ) then print *, "error: sellc chunk size not supported." return end if + vec_y = beta * vec_y + num_chunks = nrows / chunk_size rm = nrows - num_chunks * chunk_size - if( storage == sparse_full .and. op_==sparse_op_none ) then + if( storage == sparse_full .and. op==sparse_op_none ) then select case(chunk_size) #:for chunk in CHUNKS @@ -78,7 +81,7 @@ contains call chunk_kernel_rm(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:)) end if - else if( storage == sparse_full .and. op_==sparse_op_transpose ) then + else if( storage == sparse_full .and. op==sparse_op_transpose ) then select case(chunk_size) #:for chunk in CHUNKS @@ -100,7 +103,7 @@ contains end if #:if t1.startswith('complex') - else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then + else if( storage == sparse_full .and. op==sparse_op_hermitian ) then select case(chunk_size) #:for chunk in CHUNKS @@ -135,7 +138,7 @@ contains ${t1}$, intent(inout) :: y(${chunk}$) integer :: j do j = 1, n - y(:) = y(:) + alpha_ * a(:,j) * x(col(:,j)) + y(:) = y(:) + alpha * a(:,j) * x(col(:,j)) end do end subroutine pure subroutine chunk_kernel_trans_${chunk}$(n,a,col,x,y) @@ -146,7 +149,7 @@ contains integer :: j, k do j = 1, n do k = 1, ${chunk}$ - y(col(k,j)) = y(col(k,j)) + alpha_ * a(k,j) * x(k) + y(col(k,j)) = y(col(k,j)) + alpha * a(k,j) * x(k) end do end do end subroutine @@ -159,7 +162,7 @@ contains integer :: j, k do j = 1, n do k = 1, ${chunk}$ - y(col(k,j)) = y(col(k,j)) + alpha_ * conjg(a(k,j)) * x(k) + y(col(k,j)) = y(col(k,j)) + alpha * conjg(a(k,j)) * x(k) end do end do end subroutine @@ -173,7 +176,7 @@ contains ${t1}$, intent(inout) :: y(r) integer :: j do j = 1, n - y(1:r) = y(1:r) + alpha_ * a(1:r,j) * x(col(1:r,j)) + y(1:r) = y(1:r) + alpha * a(1:r,j) * x(col(1:r,j)) end do end subroutine pure subroutine chunk_kernel_rm_trans(n,cs,r,a,col,x,y) @@ -184,7 +187,7 @@ contains integer :: j, k do j = 1, n do k = 1, r - y(col(k,j)) = y(col(k,j)) + alpha_ * a(k,j) * x(k) + y(col(k,j)) = y(col(k,j)) + alpha * a(k,j) * x(k) end do end do end subroutine @@ -197,7 +200,7 @@ contains integer :: j, k do j = 1, n do k = 1, r - y(col(k,j)) = y(col(k,j)) + alpha_ * conjg(a(k,j)) * x(k) + y(col(k,j)) = y(col(k,j)) + alpha * conjg(a(k,j)) * x(k) end do end do end subroutine diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 1e174ce15..7698538cf 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -41,6 +41,8 @@ contains block integer, parameter :: wp = ${k1}$ type(COO_${s1}$_type) :: COO + ${t1}$, parameter :: alpha = 1._wp + ${t1}$, parameter :: beta = 0._wp ${t1}$, allocatable :: dense(:,:) ${t1}$, allocatable :: vec_x(:) ${t1}$, allocatable :: vec_y1(:), vec_y2(:) @@ -69,7 +71,9 @@ contains if (allocated(error)) return ! Low-level API - call spmv_kernel_coo(COO%nrows, COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, vec_x, vec_y2 ) + call spmv_kernel_coo(sparse_op_none, alpha, & + COO%nrows, COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & + vec_x, beta, vec_y2 ) call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' ) if (allocated(error)) return @@ -81,8 +85,9 @@ contains if (allocated(error)) return ! Low-level API - call spmv_kernel_coo(COO%nrows, COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, vec_y1, vec_x, & - op=sparse_op_transpose ) + call spmv_kernel_coo(sparse_op_transpose, alpha, & + COO%nrows, COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & + vec_y1, beta, vec_x) 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 @@ -125,6 +130,8 @@ contains block integer, parameter :: wp = ${k1}$ type(CSR_${s1}$_type) :: CSR + ${t1}$, parameter :: alpha = 1._wp + ${t1}$, parameter :: beta = 0._wp ${t1}$, allocatable :: vec_x(:) ${t1}$, allocatable :: vec_y(:) @@ -142,8 +149,9 @@ contains if (allocated(error)) return ! Low-level API - call spmv_kernel_csr(CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & - vec_x, vec_y ) + call spmv_kernel_csr(sparse_op_none, alpha, & + CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & + vec_x, beta, 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 @@ -155,8 +163,9 @@ contains 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_kernel_csr(CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & - vec_y, vec_x, op=sparse_op_transpose ) + call spmv_kernel_csr(sparse_op_transpose, alpha, & + CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & + vec_y, beta, vec_x) 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 @@ -170,6 +179,8 @@ contains block integer, parameter :: wp = ${k1}$ type(CSC_${s1}$_type) :: CSC + ${t1}$, parameter :: alpha = 1._wp + ${t1}$, parameter :: beta = 0._wp ${t1}$, allocatable :: vec_x(:) ${t1}$, allocatable :: vec_y(:) @@ -187,7 +198,9 @@ contains if (allocated(error)) return !High-level API - call spmv_kernel_csc(CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, vec_x, vec_y ) + call spmv_kernel_csc(sparse_op_none, alpha, & + CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & + vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return @@ -200,8 +213,9 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_csc(CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, vec_y, vec_x, & - op=sparse_op_transpose ) + call spmv_kernel_csc(sparse_op_transpose, alpha, & + CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & + vec_y, beta, vec_x) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return @@ -216,6 +230,8 @@ contains block integer, parameter :: wp = ${k1}$ type(ELL_${s1}$_type) :: ELL + ${t1}$, parameter :: alpha = 1._wp + ${t1}$, parameter :: beta = 0._wp ${t1}$, allocatable :: vec_x(:) ${t1}$, allocatable :: vec_y(:) @@ -239,8 +255,9 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_ell(ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, & - ELL%storage, vec_x, vec_y ) + call spmv_kernel_ell(sparse_op_none, alpha, & + ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & + vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) if (allocated(error)) return @@ -253,8 +270,9 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_ell(ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, & - ELL%storage, vec_y, vec_x, op=sparse_op_transpose ) + call spmv_kernel_ell(sparse_op_transpose, alpha, & + ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & + vec_y, beta, vec_x) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return end block @@ -270,6 +288,8 @@ contains integer, parameter :: wp = ${k1}$ type(SELLC_${s1}$_type) :: SELLC type(CSR_${s1}$_type) :: CSR + ${t1}$, parameter :: alpha = 1._wp + ${t1}$, parameter :: beta = 0._wp ${t1}$, allocatable :: vec_x(:) ${t1}$, allocatable :: vec_y(:), vec_y2(:) integer :: i @@ -297,8 +317,9 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_sellc(SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & - SELLC%storage, vec_x, vec_y ) + call spmv_kernel_sellc(sparse_op_none, alpha, & + SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & + vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)) ) if (allocated(error)) return @@ -314,8 +335,9 @@ contains if (allocated(error)) return !Low-level API - call spmv_kernel_sellc(SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, & - SELLC%storage, vec_x, vec_y2 , op = sparse_op_transpose ) + call spmv_kernel_sellc(sparse_op_transpose, alpha, & + SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & + vec_x, beta, vec_y2) call check(error, all(vec_y == vec_y2)) if (allocated(error)) return From 067f991f537c9fc080ad0874f67701bd0ba2dfb6 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 10:11:50 +0200 Subject: [PATCH 15/27] kernel: remove nrows --- src/sparse/stdlib_sparse_spmv.fypp | 15 +++++---------- src/sparse/stdlib_sparse_spmv_coo.fypp | 5 ++--- src/sparse/stdlib_sparse_spmv_csc.fypp | 5 ++--- src/sparse/stdlib_sparse_spmv_csr.fypp | 8 +++++--- src/sparse/stdlib_sparse_spmv_ell.fypp | 8 +++++--- src/sparse/stdlib_sparse_spmv_sellc.fypp | 8 +++++--- test/linalg/test_linalg_sparse.fypp | 20 ++++++++++---------- 7 files changed, 34 insertions(+), 35 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 8371f8058..d32ef7595 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -78,8 +78,7 @@ module stdlib_sparse_spmv interface spmv_kernel_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,nnz,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,ncols,data,index,nnz,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) @@ -102,8 +101,7 @@ module stdlib_sparse_spmv interface spmv_kernel_csc #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,colptr,row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,ncols,data,colptr,row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: colptr(:) !! matrix column pointer @@ -126,8 +124,7 @@ module stdlib_sparse_spmv interface spmv_kernel_csr #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: col(:) !! matrix column pointer @@ -150,8 +147,7 @@ module stdlib_sparse_spmv interface spmv_kernel_ell #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: index(:,:) @@ -173,9 +169,8 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves - integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index ac4f2613f..20d0f523e 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -33,13 +33,12 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_coo_${rank}$d_${s1}$(op_, alpha_, & - matrix%nrows, matrix%ncols, matrix%data, matrix%index, matrix%nnz, matrix%storage, & + matrix%ncols, matrix%data, matrix%index, matrix%nnz, matrix%storage, & vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,nnz,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,ncols,data,index,nnz,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) !! Matrix coordinates index(2,nnz) diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 6a4f83ae2..28f30568b 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -38,13 +38,12 @@ contains endif call spmv_kernel_csc_${rank}$d_${s1}$(op_,alpha_, & - matrix%nrows,matrix%ncols,matrix%data,matrix%colptr,matrix%row,matrix%storage, & + matrix%ncols,matrix%data,matrix%colptr,matrix%row,matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,colptr,row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,ncols,data,colptr,row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: colptr(:) !! matrix column pointer diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index b6a93a0cb..53d8a162c 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -33,13 +33,12 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_csr_${rank}$d_${s1}$(op_, alpha_, & - matrix%nrows, matrix%ncols, matrix%data, matrix%col, matrix%rowptr, matrix%storage, & + matrix%ncols, matrix%data, matrix%col, matrix%rowptr, matrix%storage, & vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: col(:) !! matrix column pointer @@ -51,12 +50,15 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, j + integer(ilp) :: nrows #:if rank == 1 ${t1}$ :: aux, aux2 #:else ${t1}$ :: aux(size(vec_x,dim=1)), aux2(size(vec_x,dim=1)) #:endif + nrows = size(rowptr)-1 + vec_y = beta * vec_y if( storage == sparse_full .and. op==sparse_op_none ) then diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index a6e3d7fa1..c957cbe94 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -33,14 +33,13 @@ contains if(present(beta)) beta_ = alpha call spmv_kernel_ell_${rank}$d_${s1}$(op_,alpha_, & - matrix%nrows,matrix%ncols,matrix%data,matrix%index,matrix%K, & + matrix%ncols,matrix%data,matrix%index,matrix%K, & matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,nrows,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: nrows + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: index(:,:) @@ -52,7 +51,10 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, j, k + integer(ilp) :: nrows + nrows = size(index, 1) + vec_y = beta * vec_y if( storage == sparse_full .and. op==sparse_op_none ) then diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index a1b830b09..97a12b52f 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -29,15 +29,14 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_sellc_${s1}$(op_,alpha_, & - matrix%nrows, matrix%ncols, matrix%data, matrix%rowptr, matrix%col, & + matrix%ncols, matrix%data, matrix%rowptr, matrix%col, & matrix%chunk_size, matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves - integer(ilp), intent(in) :: nrows integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) @@ -50,12 +49,15 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, nz, rowidx, num_chunks, rm + integer(ilp) :: nrows if( .not.any( ${CHUNKS}$ == chunk_size ) ) then print *, "error: sellc chunk size not supported." return end if + nrows = size(ia) - 1 + vec_y = beta * vec_y num_chunks = nrows / chunk_size diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 7698538cf..8ad7734fb 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -72,7 +72,7 @@ contains ! Low-level API call spmv_kernel_coo(sparse_op_none, alpha, & - COO%nrows, COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & + COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & vec_x, beta, vec_y2 ) call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' ) if (allocated(error)) return @@ -86,7 +86,7 @@ contains ! Low-level API call spmv_kernel_coo(sparse_op_transpose, alpha, & - COO%nrows, COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & + COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & vec_y1, beta, vec_x) 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 @@ -150,7 +150,7 @@ contains ! Low-level API call spmv_kernel_csr(sparse_op_none, alpha, & - CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & + CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSR - low-level API - wrong output' ) @@ -164,7 +164,7 @@ contains if (allocated(error)) return call spmv_kernel_csr(sparse_op_transpose, alpha, & - CSR%nrows, CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & + CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & vec_y, beta, vec_x) 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 @@ -199,7 +199,7 @@ contains !High-level API call spmv_kernel_csc(sparse_op_none, alpha, & - CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & + CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) @@ -214,7 +214,7 @@ contains !Low-level API call spmv_kernel_csc(sparse_op_transpose, alpha, & - CSC%nrows, CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & + CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & vec_y, beta, vec_x) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return @@ -256,7 +256,7 @@ contains !Low-level API call spmv_kernel_ell(sparse_op_none, alpha, & - ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & + ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)) ) @@ -271,7 +271,7 @@ contains !Low-level API call spmv_kernel_ell(sparse_op_transpose, alpha, & - ELL%nrows, ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & + ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & vec_y, beta, vec_x) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) ) if (allocated(error)) return @@ -318,7 +318,7 @@ contains !Low-level API call spmv_kernel_sellc(sparse_op_none, alpha, & - SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & + SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)) ) @@ -336,7 +336,7 @@ contains !Low-level API call spmv_kernel_sellc(sparse_op_transpose, alpha, & - SELLC%nrows, SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & + SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & vec_x, beta, vec_y2) call check(error, all(vec_y == vec_y2)) From 5a2390768fce78196a33c7b64603ee7cbfaab2ac Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:25:59 +0200 Subject: [PATCH 16/27] remove nrows and ncol --- src/sparse/stdlib_sparse_spmv.fypp | 16 ++++------- src/sparse/stdlib_sparse_spmv_coo.fypp | 5 ++-- src/sparse/stdlib_sparse_spmv_csc.fypp | 8 ++++-- src/sparse/stdlib_sparse_spmv_csr.fypp | 5 ++-- src/sparse/stdlib_sparse_spmv_ell.fypp | 6 ++-- src/sparse/stdlib_sparse_spmv_sellc.fypp | 15 +++++----- test/linalg/test_linalg_sparse.fypp | 36 ++++++++++++------------ 7 files changed, 42 insertions(+), 49 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index d32ef7595..e67a6027c 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -78,8 +78,7 @@ module stdlib_sparse_spmv interface spmv_kernel_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,ncols,data,index,nnz,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,nnz,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) integer, intent(in) :: storage @@ -101,8 +100,7 @@ module stdlib_sparse_spmv interface spmv_kernel_csc #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,ncols,data,colptr,row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: colptr(:) !! matrix column pointer integer(ilp), intent(in) :: row(:) !! matrix row pointer @@ -124,8 +122,7 @@ module stdlib_sparse_spmv interface spmv_kernel_csr #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: col(:) !! matrix column pointer integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer @@ -147,8 +144,7 @@ module stdlib_sparse_spmv interface spmv_kernel_ell #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,mnz_p_row,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: index(:,:) integer(ilp), intent(in) :: mnz_p_row @@ -169,12 +165,12 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,chunk_size,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves - integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) integer(ilp), intent(in) :: ja(:,:) + integer, intent(in) :: nrows integer, intent(in) :: chunk_size integer, intent(in) :: storage ${t1}$, intent(in) :: vec_x(:) diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 20d0f523e..6f7476f70 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -33,13 +33,12 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_coo_${rank}$d_${s1}$(op_, alpha_, & - matrix%ncols, matrix%data, matrix%index, matrix%nnz, matrix%storage, & + matrix%data, matrix%index, matrix%nnz, matrix%storage, & vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,ncols,data,index,nnz,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,nnz,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: index(:,:) !! Matrix coordinates index(2,nnz) integer(ilp), intent(in) :: nnz !! number of non-zero values diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 28f30568b..e10102771 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -38,13 +38,12 @@ contains endif call spmv_kernel_csc_${rank}$d_${s1}$(op_,alpha_, & - matrix%ncols,matrix%data,matrix%colptr,matrix%row,matrix%storage, & + matrix%data,matrix%colptr,matrix%row,matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,ncols,data,colptr,row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: colptr(:) !! matrix column pointer integer(ilp), intent(in) :: row(:) !! matrix row pointer @@ -55,12 +54,15 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, j + integer(ilp) :: ncols #:if rank == 1 ${t1}$ :: aux #:else ${t1}$ :: aux(size(vec_x,dim=1)) #:endif + ncols = size(colptr)-1 + vec_y = beta * vec_y if( storage == sparse_full .and. op==sparse_op_none ) then diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index 53d8a162c..da303ee3c 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -33,13 +33,12 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_csr_${rank}$d_${s1}$(op_, alpha_, & - matrix%ncols, matrix%data, matrix%col, matrix%rowptr, matrix%storage, & + matrix%data, matrix%col, matrix%rowptr, matrix%storage, & vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,ncols,data,col,rowptr,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:) integer(ilp), intent(in) :: col(:) !! matrix column pointer integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index c957cbe94..82f773efd 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -33,14 +33,12 @@ contains if(present(beta)) beta_ = alpha call spmv_kernel_ell_${rank}$d_${s1}$(op_,alpha_, & - matrix%ncols,matrix%data,matrix%index,matrix%K, & - matrix%storage, & + matrix%data,matrix%index,matrix%K,matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,ncols,data,index,mnz_p_row,storage,vec_x,beta,vec_y) - integer(ilp), intent(in) :: ncols + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,mnz_p_row,storage,vec_x,beta,vec_y) ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: index(:,:) integer, intent(in) :: mnz_p_row diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 97a12b52f..68e35ae26 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -29,18 +29,18 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_sellc_${s1}$(op_,alpha_, & - matrix%ncols, matrix%data, matrix%rowptr, matrix%col, & - matrix%chunk_size, matrix%storage, & + matrix%data, matrix%rowptr, matrix%col, & + matrix%nrows, matrix%chunk_size, matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,ncols,data,ia,ja,chunk_size,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,chunk_size,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves - integer(ilp), intent(in) :: ncols ${t1}$, intent(in) :: data(:,:) integer(ilp), intent(in) :: ia(:) integer(ilp), intent(in) :: ja(:,:) + integer(ilp), intent(in) :: nrows integer, intent(in) :: chunk_size integer, intent(in) :: storage ${t1}$, intent(in) :: vec_x(:) @@ -49,19 +49,18 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, nz, rowidx, num_chunks, rm - integer(ilp) :: nrows if( .not.any( ${CHUNKS}$ == chunk_size ) ) then print *, "error: sellc chunk size not supported." return end if - nrows = size(ia) - 1 - vec_y = beta * vec_y - num_chunks = nrows / chunk_size + num_chunks = size(ia)-1 !nrows / chunk_size + rm = nrows - num_chunks * chunk_size + if( storage == sparse_full .and. op==sparse_op_none ) then select case(chunk_size) diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 8ad7734fb..c3ccd552e 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -72,7 +72,7 @@ contains ! Low-level API call spmv_kernel_coo(sparse_op_none, alpha, & - COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & + COO%data ,COO%index, COO%nnz, COO%storage, & vec_x, beta, vec_y2 ) call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' ) if (allocated(error)) return @@ -86,7 +86,7 @@ contains ! Low-level API call spmv_kernel_coo(sparse_op_transpose, alpha, & - COO%ncols, COO%data ,COO%index, COO%nnz, COO%storage, & + COO%data ,COO%index, COO%nnz, COO%storage, & vec_y1, beta, vec_x) 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 @@ -150,7 +150,7 @@ contains ! Low-level API call spmv_kernel_csr(sparse_op_none, alpha, & - CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & + CSR%data, CSR%col, CSR%rowptr, CSR%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSR - low-level API - wrong output' ) @@ -164,7 +164,7 @@ contains if (allocated(error)) return call spmv_kernel_csr(sparse_op_transpose, alpha, & - CSR%ncols, CSR%data, CSR%col, CSR%rowptr, CSR%storage, & + CSR%data, CSR%col, CSR%rowptr, CSR%storage, & vec_y, beta, vec_x) 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 @@ -199,10 +199,10 @@ contains !High-level API call spmv_kernel_csc(sparse_op_none, alpha, & - CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & + CSC%data, CSC%colptr, CSC%row, CSC%storage, & vec_x, beta, 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)), 'CSC - low-level API - wrong output' ) if (allocated(error)) return ! Test in-place transpose @@ -214,9 +214,9 @@ contains !Low-level API call spmv_kernel_csc(sparse_op_transpose, alpha, & - CSC%ncols, CSC%data, CSC%colptr, CSC%row, CSC%storage, & + CSC%data, CSC%colptr, CSC%row, CSC%storage, & vec_y, beta, vec_x) - 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)), 'CSC tranpose - low-level API - wrong output' ) if (allocated(error)) return end block @@ -256,10 +256,10 @@ contains !Low-level API call spmv_kernel_ell(sparse_op_none, alpha, & - ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & + ELL%data, ELL%index, ELL%K, ELL%storage, & vec_x, beta, 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)), 'ELL - low-level API - wrong output' ) if (allocated(error)) return ! Test in-place transpose @@ -271,9 +271,9 @@ contains !Low-level API call spmv_kernel_ell(sparse_op_transpose, alpha, & - ELL%ncols, ELL%data, ELL%index, ELL%K, ELL%storage, & + ELL%data, ELL%index, ELL%K, ELL%storage, & vec_y, beta, vec_x) - 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)), 'ELL tranpose - low-level API - wrong output' ) if (allocated(error)) return end block #:endfor @@ -313,15 +313,15 @@ contains !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)) ) + call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)), 'SELLC - high-level API - wrong output' ) if (allocated(error)) return !Low-level API call spmv_kernel_sellc(sparse_op_none, alpha, & - SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & + SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%chunk_size, SELLC%storage, & vec_x, beta, vec_y ) - call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)) ) + call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)), 'SELLC - low-level API - wrong output' ) if (allocated(error)) return ! Test in-place transpose @@ -331,15 +331,15 @@ contains !High-level API call spmv( SELLC, vec_x, vec_y2 , op = sparse_op_transpose ) - call check(error, all(vec_y == vec_y2)) + call check(error, all(vec_y == vec_y2), 'SELLC tranpose - high-level API - wrong output' ) if (allocated(error)) return !Low-level API call spmv_kernel_sellc(sparse_op_transpose, alpha, & - SELLC%ncols, SELLC%data, SELLC%rowptr, SELLC%col, SELLC%chunk_size, SELLC%storage, & + SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%chunk_size, SELLC%storage, & vec_x, beta, vec_y2) - call check(error, all(vec_y == vec_y2)) + call check(error, all(vec_y == vec_y2), 'SELLC tranpose - low-level API - wrong output' ) if (allocated(error)) return end block From d1b7b4c189c02466305228a5425dae31cd3967b7 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:31:30 +0200 Subject: [PATCH 17/27] fix spmv kernel sellc --- src/sparse/stdlib_sparse_spmv_sellc.fypp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 68e35ae26..2558d7a25 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -57,8 +57,7 @@ contains vec_y = beta * vec_y - num_chunks = size(ia)-1 !nrows / chunk_size - + num_chunks = nrows / chunk_size rm = nrows - num_chunks * chunk_size if( storage == sparse_full .and. op==sparse_op_none ) then From ffb0f262a3c36cd492c9bdad31e5c48eee8091a0 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:40:36 +0200 Subject: [PATCH 18/27] addition of contiguous --- src/sparse/stdlib_sparse_spmv.fypp | 46 ++++++++++++------------ src/sparse/stdlib_sparse_spmv_coo.fypp | 8 ++--- src/sparse/stdlib_sparse_spmv_csc.fypp | 10 +++--- src/sparse/stdlib_sparse_spmv_csr.fypp | 10 +++--- src/sparse/stdlib_sparse_spmv_ell.fypp | 8 ++--- src/sparse/stdlib_sparse_spmv_sellc.fypp | 10 +++--- 6 files changed, 46 insertions(+), 46 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index e67a6027c..0a7f89af2 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -79,12 +79,12 @@ module stdlib_sparse_spmv #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,nnz,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: index(:,:) + ${t1}$, intent(in), contiguous :: data(:) + integer(ilp), intent(in), contiguous :: 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), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op @@ -101,12 +101,12 @@ module stdlib_sparse_spmv #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: colptr(:) !! matrix column pointer - integer(ilp), intent(in) :: row(:) !! matrix row pointer + ${t1}$, intent(in), contiguous :: data(:) + integer(ilp), intent(in), contiguous :: colptr(:) !! matrix column pointer + integer(ilp), intent(in), contiguous :: row(:) !! matrix row pointer integer, intent(in) :: storage !! storage - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op @@ -123,12 +123,12 @@ module stdlib_sparse_spmv #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: col(:) !! matrix column pointer - integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer + ${t1}$, intent(in), contiguous :: data(:) + integer(ilp), intent(in), contiguous :: col(:) !! matrix column pointer + integer(ilp), intent(in), contiguous :: rowptr(:) !! matrix row pointer integer, intent(in) :: storage !! storage - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op @@ -145,12 +145,12 @@ module stdlib_sparse_spmv #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,mnz_p_row,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:,:) - integer(ilp), intent(in) :: index(:,:) + ${t1}$, intent(in), contiguous :: data(:,:) + integer(ilp), intent(in), contiguous :: index(:,:) integer(ilp), intent(in) :: mnz_p_row integer, intent(in) :: storage - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op @@ -167,14 +167,14 @@ module stdlib_sparse_spmv #:for k1, t1, s1 in (KINDS_TYPES) module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,chunk_size,storage,vec_x,beta,vec_y) !! 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(:,:) + ${t1}$, intent(in), contiguous :: data(:,:) + integer(ilp), intent(in), contiguous :: ia(:) + integer(ilp), intent(in), contiguous :: ja(:,:) integer, intent(in) :: nrows integer, intent(in) :: chunk_size integer, intent(in) :: storage - ${t1}$, intent(in) :: vec_x(:) - ${t1}$, intent(inout) :: vec_y(:) + ${t1}$, intent(in), contiguous :: vec_x(:) + ${t1}$, intent(inout), contiguous :: vec_y(:) ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 6f7476f70..922a59614 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -39,12 +39,12 @@ contains end subroutine module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,nnz,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: index(:,:) !! Matrix coordinates index(2,nnz) + ${t1}$, intent(in), contiguous :: data(:) + integer(ilp), intent(in), contiguous :: index(:,:) !! Matrix coordinates index(2,nnz) 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), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index e10102771..23eb9e900 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -44,12 +44,12 @@ contains end subroutine module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: colptr(:) !! matrix column pointer - integer(ilp), intent(in) :: row(:) !! matrix row pointer + ${t1}$, intent(in), contiguous :: data(:) + integer(ilp), intent(in), contiguous :: colptr(:) !! matrix column pointer + integer(ilp), intent(in), contiguous :: row(:) !! matrix row pointer integer, intent(in) :: storage !! storage - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index da303ee3c..f9ae7fd88 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -39,12 +39,12 @@ contains end subroutine module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:) - integer(ilp), intent(in) :: col(:) !! matrix column pointer - integer(ilp), intent(in) :: rowptr(:) !! matrix row pointer + ${t1}$, intent(in), contiguous :: data(:) + integer(ilp), intent(in), contiguous :: col(:) !! matrix column pointer + integer(ilp), intent(in), contiguous :: rowptr(:) !! matrix row pointer integer, intent(in) :: storage !! storage - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index 82f773efd..4a11ae478 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -39,12 +39,12 @@ contains end subroutine module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,mnz_p_row,storage,vec_x,beta,vec_y) - ${t1}$, intent(in) :: data(:,:) - integer(ilp), intent(in) :: index(:,:) + ${t1}$, intent(in), contiguous :: data(:,:) + integer(ilp), intent(in), contiguous :: index(:,:) integer, intent(in) :: mnz_p_row integer, intent(in) :: storage - ${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$ - ${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$ + ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ + ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 2558d7a25..85daf89b8 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -37,14 +37,14 @@ contains module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,chunk_size,storage,vec_x,beta,vec_y) !! 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(:,:) + ${t1}$, intent(in), contiguous :: data(:,:) + integer(ilp), intent(in), contiguous :: ia(:) + integer(ilp), intent(in), contiguous :: ja(:,:) integer(ilp), intent(in) :: nrows integer, intent(in) :: chunk_size integer, intent(in) :: storage - ${t1}$, intent(in) :: vec_x(:) - ${t1}$, intent(inout) :: vec_y(:) + ${t1}$, intent(in), contiguous :: vec_x(:) + ${t1}$, intent(inout), contiguous :: vec_y(:) ${t1}$, intent(in) :: alpha ${t1}$, intent(in) :: beta character(1), intent(in) :: op From 6b8e8776e60ba3a6c9a1ac5de4fa9896c423d3b7 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:41:29 +0200 Subject: [PATCH 19/27] Update doc/specs/stdlib_sparse.md MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: José Alves <102541118+jalvesz@users.noreply.github.com> --- doc/specs/stdlib_sparse.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_sparse.md b/doc/specs/stdlib_sparse.md index bb749d068..3a1d5cec9 100644 --- a/doc/specs/stdlib_sparse.md +++ b/doc/specs/stdlib_sparse.md @@ -240,7 +240,7 @@ $$y=\alpha*op(M)*x+\beta*y$$ ### Arguments -For the all formats: +Common arguments for all formats `nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the number of rows in the sparse matrix. From 856a27d2ad0b31e0c99e5023fa87ebdc489e9d32 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:46:41 +0200 Subject: [PATCH 20/27] fix format --- src/sparse/stdlib_sparse_spmv_coo.fypp | 2 -- src/sparse/stdlib_sparse_spmv_csc.fypp | 7 +------ src/sparse/stdlib_sparse_spmv_csr.fypp | 2 -- src/sparse/stdlib_sparse_spmv_ell.fypp | 2 -- src/sparse/stdlib_sparse_spmv_sellc.fypp | 2 -- 5 files changed, 1 insertion(+), 14 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 922a59614..8aa229d5d 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -20,12 +20,10 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - ${t1}$ :: alpha_, beta_ character(1) :: op_ op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ if(present(alpha)) alpha_ = alpha diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 23eb9e900..2d1846779 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -20,22 +20,17 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - ${t1}$ :: alpha_, beta_ character(1) :: op_ op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ if(present(alpha)) alpha_ = alpha beta_ = zero_${k1}$ if(present(beta)) beta_ = beta - if(present(beta)) then - else - vec_y = zero_${s1}$ - endif + if(present(beta))vec_y = zero_${s1}$ call spmv_kernel_csc_${rank}$d_${s1}$(op_,alpha_, & matrix%data,matrix%colptr,matrix%row,matrix%storage, & diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp index f9ae7fd88..cfe3dee61 100644 --- a/src/sparse/stdlib_sparse_spmv_csr.fypp +++ b/src/sparse/stdlib_sparse_spmv_csr.fypp @@ -20,12 +20,10 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - ${t1}$ :: alpha_, beta_ character(1) :: op_ op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ if(present(alpha)) alpha_ = alpha diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index 4a11ae478..8cc636ec5 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -20,12 +20,10 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - ${t1}$ :: alpha_, beta_ character(1) :: op_ op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${k1}$ if(present(alpha)) alpha_ = alpha diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 85daf89b8..837f7075f 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -16,12 +16,10 @@ contains ${t1}$, intent(in), optional :: alpha ${t1}$, intent(in), optional :: beta character(1), intent(in), optional :: op - ${t1}$ :: alpha_, beta_ character(1) :: op_ op_ = sparse_op_none; if(present(op)) op_ = op - alpha_ = one_${s1}$ if(present(alpha)) alpha_ = alpha From dc2e00b2e13120bca844f3a66acb535e8c81ae3f Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:48:48 +0200 Subject: [PATCH 21/27] fix csc --- src/sparse/stdlib_sparse_spmv_csc.fypp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp index 2d1846779..52ee6d3e8 100644 --- a/src/sparse/stdlib_sparse_spmv_csc.fypp +++ b/src/sparse/stdlib_sparse_spmv_csc.fypp @@ -30,8 +30,6 @@ contains beta_ = zero_${k1}$ if(present(beta)) beta_ = beta - if(present(beta))vec_y = zero_${s1}$ - call spmv_kernel_csc_${rank}$d_${s1}$(op_,alpha_, & matrix%data,matrix%colptr,matrix%row,matrix%storage, & vec_x,beta_,vec_y) From 12d1a62f2c3456e7c5fecefeff62d5a82db3b075 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 11:52:36 +0200 Subject: [PATCH 22/27] fix spmv ell --- src/sparse/stdlib_sparse_spmv_ell.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index 8cc636ec5..a15ea1a0c 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -28,7 +28,7 @@ contains if(present(alpha)) alpha_ = alpha beta_ = zero_${k1}$ - if(present(beta)) beta_ = alpha + if(present(beta)) beta_ = beta call spmv_kernel_ell_${rank}$d_${s1}$(op_,alpha_, & matrix%data,matrix%index,matrix%K,matrix%storage, & From 4a4ed9bce823c6d16b3340483991279c9ffe2205 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 12:04:55 +0200 Subject: [PATCH 23/27] remove nnz from spmv coo --- src/sparse/stdlib_sparse_spmv.fypp | 3 +-- src/sparse/stdlib_sparse_spmv_coo.fypp | 8 +++++--- test/linalg/test_linalg_sparse.fypp | 4 ++-- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 0a7f89af2..679f65e06 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -78,11 +78,10 @@ module stdlib_sparse_spmv interface spmv_kernel_coo #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,nnz,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,storage,vec_x,beta,vec_y) ${t1}$, intent(in), contiguous :: data(:) integer(ilp), intent(in), contiguous :: index(:,:) integer, intent(in) :: storage - integer(ilp), intent(in) :: nnz ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ ${t1}$, intent(in) :: alpha diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 8aa229d5d..e755a8ffc 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -31,15 +31,14 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_coo_${rank}$d_${s1}$(op_, alpha_, & - matrix%data, matrix%index, matrix%nnz, matrix%storage, & + matrix%data, matrix%index, matrix%storage, & vec_x, beta_, vec_y) end subroutine - module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,nnz,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_coo_${rank}$d_${s1}$(op,alpha,data,index,storage,vec_x,beta,vec_y) ${t1}$, intent(in), contiguous :: data(:) integer(ilp), intent(in), contiguous :: index(:,:) !! Matrix coordinates index(2,nnz) - integer(ilp), intent(in) :: nnz !! number of non-zero values integer, intent(in) :: storage !! storage ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ @@ -47,6 +46,9 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: col_index, k, row_index + integer(ilp) :: nnz !! number of non-zero values + + nnz = size(index, 2) vec_y = beta * vec_y diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index c3ccd552e..3ea27d10a 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -72,7 +72,7 @@ contains ! Low-level API call spmv_kernel_coo(sparse_op_none, alpha, & - COO%data ,COO%index, COO%nnz, COO%storage, & + COO%data ,COO%index, COO%storage, & vec_x, beta, vec_y2 ) call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' ) if (allocated(error)) return @@ -86,7 +86,7 @@ contains ! Low-level API call spmv_kernel_coo(sparse_op_transpose, alpha, & - COO%data ,COO%index, COO%nnz, COO%storage, & + COO%data ,COO%index, COO%storage, & vec_y1, beta, vec_x) 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 From 0c33ea8d4950955be712bd836a195060ef7420da Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 15:11:37 +0200 Subject: [PATCH 24/27] remove chunck_size --- src/sparse/stdlib_sparse_spmv.fypp | 3 +-- src/sparse/stdlib_sparse_spmv_sellc.fypp | 8 +++++--- test/linalg/test_linalg_sparse.fypp | 4 ++-- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 679f65e06..e1f92139b 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -164,13 +164,12 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,chunk_size,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves ${t1}$, intent(in), contiguous :: data(:,:) integer(ilp), intent(in), contiguous :: ia(:) integer(ilp), intent(in), contiguous :: ja(:,:) integer, intent(in) :: nrows - integer, intent(in) :: chunk_size integer, intent(in) :: storage ${t1}$, intent(in), contiguous :: vec_x(:) ${t1}$, intent(inout), contiguous :: vec_y(:) diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index 837f7075f..c7ec9b1a9 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -28,18 +28,17 @@ contains call spmv_kernel_sellc_${s1}$(op_,alpha_, & matrix%data, matrix%rowptr, matrix%col, & - matrix%nrows, matrix%chunk_size, matrix%storage, & + matrix%nrows, matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,chunk_size,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves ${t1}$, intent(in), contiguous :: data(:,:) integer(ilp), intent(in), contiguous :: ia(:) integer(ilp), intent(in), contiguous :: ja(:,:) integer(ilp), intent(in) :: nrows - integer, intent(in) :: chunk_size integer, intent(in) :: storage ${t1}$, intent(in), contiguous :: vec_x(:) ${t1}$, intent(inout), contiguous :: vec_y(:) @@ -47,6 +46,9 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, nz, rowidx, num_chunks, rm + integer :: chunk_size + + chunk_size = size(data, 1) if( .not.any( ${CHUNKS}$ == chunk_size ) ) then print *, "error: sellc chunk size not supported." diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 3ea27d10a..5bdf27d47 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -318,7 +318,7 @@ contains !Low-level API call spmv_kernel_sellc(sparse_op_none, alpha, & - SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%chunk_size, SELLC%storage, & + SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)), 'SELLC - low-level API - wrong output' ) @@ -336,7 +336,7 @@ contains !Low-level API call spmv_kernel_sellc(sparse_op_transpose, alpha, & - SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%chunk_size, SELLC%storage, & + SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%storage, & vec_x, beta, vec_y2) call check(error, all(vec_y == vec_y2), 'SELLC tranpose - low-level API - wrong output' ) From 98d89814586d7ecdc04ff8c7d4111489bc41ae8b Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 15:21:15 +0200 Subject: [PATCH 25/27] sellc spmv: remove nrows from API --- src/sparse/stdlib_sparse_spmv.fypp | 3 +-- src/sparse/stdlib_sparse_spmv_sellc.fypp | 8 +++++--- test/linalg/test_linalg_sparse.fypp | 4 ++-- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index e1f92139b..5f56e7602 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -164,12 +164,11 @@ module stdlib_sparse_spmv !! [Specifications](../page/specs/stdlib_sparse.html#spmv) interface spmv_kernel_sellc #:for k1, t1, s1 in (KINDS_TYPES) - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves ${t1}$, intent(in), contiguous :: data(:,:) integer(ilp), intent(in), contiguous :: ia(:) integer(ilp), intent(in), contiguous :: ja(:,:) - integer, intent(in) :: nrows integer, intent(in) :: storage ${t1}$, intent(in), contiguous :: vec_x(:) ${t1}$, intent(inout), contiguous :: vec_y(:) diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp index c7ec9b1a9..203336cd2 100644 --- a/src/sparse/stdlib_sparse_spmv_sellc.fypp +++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp @@ -28,17 +28,16 @@ contains call spmv_kernel_sellc_${s1}$(op_,alpha_, & matrix%data, matrix%rowptr, matrix%col, & - matrix%nrows, matrix%storage, & + matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,nrows,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,storage,vec_x,beta,vec_y) !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves ${t1}$, intent(in), contiguous :: data(:,:) integer(ilp), intent(in), contiguous :: ia(:) integer(ilp), intent(in), contiguous :: ja(:,:) - integer(ilp), intent(in) :: nrows integer, intent(in) :: storage ${t1}$, intent(in), contiguous :: vec_x(:) ${t1}$, intent(inout), contiguous :: vec_y(:) @@ -46,6 +45,7 @@ contains ${t1}$, intent(in) :: beta character(1), intent(in) :: op integer(ilp) :: i, nz, rowidx, num_chunks, rm + integer(ilp) :: nrows integer :: chunk_size chunk_size = size(data, 1) @@ -57,6 +57,8 @@ contains vec_y = beta * vec_y + nrows = merge(size(vec_y), size(vec_x), op==sparse_op_none) + num_chunks = nrows / chunk_size rm = nrows - num_chunks * chunk_size diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 5bdf27d47..541550b52 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -318,7 +318,7 @@ contains !Low-level API call spmv_kernel_sellc(sparse_op_none, alpha, & - SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%storage, & + SELLC%data, SELLC%rowptr, SELLC%col, SELLC%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)), 'SELLC - low-level API - wrong output' ) @@ -336,7 +336,7 @@ contains !Low-level API call spmv_kernel_sellc(sparse_op_transpose, alpha, & - SELLC%data, SELLC%rowptr, SELLC%col, SELLC%nrows, SELLC%storage, & + SELLC%data, SELLC%rowptr, SELLC%col, SELLC%storage, & vec_x, beta, vec_y2) call check(error, all(vec_y == vec_y2), 'SELLC tranpose - low-level API - wrong output' ) From 041e3eff2ff5515285271e4fda1fbd5b8433ea92 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Fri, 31 Jul 2026 15:30:55 +0200 Subject: [PATCH 26/27] spmv_kernel ELL: remove mnz_p_row --- src/sparse/stdlib_sparse_spmv.fypp | 3 +-- src/sparse/stdlib_sparse_spmv_ell.fypp | 7 ++++--- test/linalg/test_linalg_sparse.fypp | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 5f56e7602..c21054024 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -143,10 +143,9 @@ module stdlib_sparse_spmv interface spmv_kernel_ell #:for k1, t1, s1 in (KINDS_TYPES) #:for rank in RANKS - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,mnz_p_row,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,storage,vec_x,beta,vec_y) ${t1}$, intent(in), contiguous :: data(:,:) integer(ilp), intent(in), contiguous :: index(:,:) - integer(ilp), intent(in) :: mnz_p_row integer, intent(in) :: storage ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp index a15ea1a0c..9f268c155 100644 --- a/src/sparse/stdlib_sparse_spmv_ell.fypp +++ b/src/sparse/stdlib_sparse_spmv_ell.fypp @@ -31,15 +31,14 @@ contains if(present(beta)) beta_ = beta call spmv_kernel_ell_${rank}$d_${s1}$(op_,alpha_, & - matrix%data,matrix%index,matrix%K,matrix%storage, & + matrix%data,matrix%index,matrix%storage, & vec_x,beta_,vec_y) end subroutine - module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,mnz_p_row,storage,vec_x,beta,vec_y) + module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,storage,vec_x,beta,vec_y) ${t1}$, intent(in), contiguous :: data(:,:) integer(ilp), intent(in), contiguous :: index(:,:) - integer, intent(in) :: mnz_p_row integer, intent(in) :: storage ${t1}$, intent(in), contiguous :: vec_x${ranksuffix(rank)}$ ${t1}$, intent(inout), contiguous :: vec_y${ranksuffix(rank)}$ @@ -48,8 +47,10 @@ contains character(1), intent(in) :: op integer(ilp) :: i, j, k integer(ilp) :: nrows + integer :: mnz_p_row nrows = size(index, 1) + mnz_p_row = size(index, 2) vec_y = beta * vec_y diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp index 541550b52..d03c1cc5e 100644 --- a/test/linalg/test_linalg_sparse.fypp +++ b/test/linalg/test_linalg_sparse.fypp @@ -256,7 +256,7 @@ contains !Low-level API call spmv_kernel_ell(sparse_op_none, alpha, & - ELL%data, ELL%index, ELL%K, ELL%storage, & + ELL%data, ELL%index, ELL%storage, & vec_x, beta, vec_y ) call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'ELL - low-level API - wrong output' ) @@ -271,7 +271,7 @@ contains !Low-level API call spmv_kernel_ell(sparse_op_transpose, alpha, & - ELL%data, ELL%index, ELL%K, ELL%storage, & + ELL%data, ELL%index, ELL%storage, & vec_y, beta, vec_x) call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'ELL tranpose - low-level API - wrong output' ) if (allocated(error)) return From 259b368f7a7dceed501d8d50ef595654c539d491 Mon Sep 17 00:00:00 2001 From: Jeremie Vandenplas Date: Tue, 1 Sep 2026 14:16:41 +0200 Subject: [PATCH 27/27] Update spmv_kernel specs in stdlib_sparse.md --- doc/specs/stdlib_sparse.md | 31 +++++++++---------------------- 1 file changed, 9 insertions(+), 22 deletions(-) diff --git a/doc/specs/stdlib_sparse.md b/doc/specs/stdlib_sparse.md index 3a1d5cec9..369c2b792 100644 --- a/doc/specs/stdlib_sparse.md +++ b/doc/specs/stdlib_sparse.md @@ -232,31 +232,27 @@ $$y=\alpha*op(M)*x+\beta*y$$ ### Syntax -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_coo(interface)]] `(nrows,ncols,data,index,nnz,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(nrows,ncols,data,colptr,row,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(nrows,ncols,data,col,rowptr,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(nrows,ncols,data,index,mnz_p_row,storage,vec_x,vec_y [,alpha,beta,op])` -`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(nrows,ncols,data,ia,ja,chunk_size,storage,vec_x,vec_y [,alpha,beta,op])` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_coo(interface)]] `(op,alpha,data,index,storage,vec_x,beta,vec_y)` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y)` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y)` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(op,alpha,data,index,storage,vec_x,beta,vec_y)` +`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(op,alpha,data,ia,ja,storage,vec_x,beta,vec_y)` ### Arguments Common arguments for all formats -`nrows`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the number of rows in the sparse matrix. +`op`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian` -`ncols`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the number of columns in the sparse matrix. +`alpha`: Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument. `storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. It defines the symmetry storage of the sparse matrix and must be one of `sparse_full`, `sparse_lower`, or `sparse_upper`. `vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. -`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument. - -`alpha`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument. +`beta`: Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument. -`beta`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument. - -`op`, `optional`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian` +`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument. For the `COO` format: @@ -264,9 +260,6 @@ For the `COO` format: `index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`nnz`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. It specifies the number of non-zero values in `index`. - - For the `CSC` format: `data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument. @@ -284,16 +277,12 @@ For the `CSR` format: `rowptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument. - For the `ELL` format: `data`: Shall be a rank-2 array of `real` or `complex` type. It is an `intent(in)` argument. `index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`mnz_p_row`: Shall be a scalar of `integer(ilp)` type. It is an `intent(in)` argument. It specifies the constant number of elements per row $K$. - - For the `SELLC` format: `data`: Shall be a rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument. @@ -302,8 +291,6 @@ For the `SELLC` format: `ja`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument. -`chunk_size`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. - ## Sparse matrix to matrix conversions