diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp index 5524ef7c0..6f408f2b0 100644 --- a/src/sparse/stdlib_sparse_spmv.fypp +++ b/src/sparse/stdlib_sparse_spmv.fypp @@ -70,6 +70,134 @@ 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, nnz, storage, 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 + + !! Version experimental + !! + !! 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) + #: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 + + !! 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 + + !! 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_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp index 7aba6e9d8..a7a362b08 100644 --- a/src/sparse/stdlib_sparse_spmv_coo.fypp +++ b/src/sparse/stdlib_sparse_spmv_coo.fypp @@ -20,6 +20,22 @@ 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%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) + ${t1}$, intent(in) :: data(:) + integer(ilp), intent(in) :: 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), 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 +49,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 +107,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 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/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/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/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 0ebe235b5..396bae95e 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%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 + ! 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%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 + end block #:endfor end subroutine @@ -121,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 @@ -152,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 @@ -189,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 @@ -232,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 @@ -715,4 +788,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