Skip to content

Commit 4c8521d

Browse files
jalveszCopilotCopilotjvdp1
authored
housekeeping (sparse): split spmv kernels into submodules (#1180)
* split sparse_spmv into module submodules * cleanup unrelated changes * Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> * remove residual change * revert to plain do loops * fix indexing issue * Remove extra blank line in stdlib_sparse_kinds.fypp * add fixes to CSC format * make reference arrays in test static Co-authored-by: Copilot <copilot@github.com> * Update test/linalg/test_linalg_sparse.fypp Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com> * Update test/linalg/test_linalg_sparse.fypp Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com> --------- Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> Co-authored-by: Copilot <copilot@github.com> Co-authored-by: Jeremie Vandenplas <jeremie.vandenplas@gmail.com>
1 parent b1670f0 commit 4c8521d

9 files changed

Lines changed: 781 additions & 590 deletions

‎doc/specs/stdlib_sparse.md‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -146,7 +146,10 @@ Type-bound procedures to enable adding data in a sparse matrix.
146146

147147
### Syntax
148148

149-
`call matrix%add(i,j,v)` or
149+
* Add single value
150+
`call matrix%add(i,j,v)`
151+
152+
* Add a block of values
150153
`call matrix%add(i(:),j(:),v(:,:))`
151154

152155
### Arguments

‎src/sparse/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,11 @@ set(sparse_fppFiles
44
stdlib_sparse_kinds.fypp
55
stdlib_sparse_operators.fypp
66
stdlib_sparse_spmv.fypp
7+
stdlib_sparse_spmv_coo.fypp
8+
stdlib_sparse_spmv_csr.fypp
9+
stdlib_sparse_spmv_csc.fypp
10+
stdlib_sparse_spmv_ell.fypp
11+
stdlib_sparse_spmv_sellc.fypp
712
)
813

914
set(sparse_cppFiles

‎src/sparse/stdlib_sparse_spmv.fypp‎

Lines changed: 42 additions & 560 deletions
Large diffs are not rendered by default.

‎src/sparse/stdlib_sparse_spmv_coo.fypp‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#:include "common.fypp"
2+
#:set RANKS = range(1, 2+1)
3+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
4+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
5+
#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
6+
#! define ranks without parentheses
7+
#:def rksfx2(rank)
8+
#{if rank > 0}#${":," + ":," * (rank - 1)}$#{endif}#
9+
#:enddef
10+
submodule (stdlib_sparse_spmv) stdlib_sparse_spmv_coo
11+
contains
12+
13+
!! spmv_coo
14+
#:for k1, t1, s1 in (KINDS_TYPES)
15+
#:for rank in RANKS
16+
module subroutine spmv_coo_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
17+
type(COO_${s1}$_type), intent(in) :: matrix
18+
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
19+
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
20+
${t1}$, intent(in), optional :: alpha
21+
${t1}$, intent(in), optional :: beta
22+
character(1), intent(in), optional :: op
23+
${t1}$ :: alpha_
24+
character(1) :: op_
25+
integer(ilp) :: col_index, k, row_index
26+
27+
op_ = sparse_op_none; if(present(op)) op_ = op
28+
alpha_ = one_${k1}$
29+
if(present(alpha)) alpha_ = alpha
30+
if(present(beta)) then
31+
vec_y = beta * vec_y
32+
else
33+
vec_y = zero_${s1}$
34+
endif
35+
36+
associate( data => matrix%data, index => matrix%index, storage => matrix%storage, nnz => matrix%nnz )
37+
select case(op_)
38+
case(sparse_op_none)
39+
if(storage == sparse_full) then
40+
do k = 1, nnz
41+
row_index = index(1,k)
42+
col_index = index(2,k)
43+
vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index)
44+
end do
45+
46+
else
47+
do k = 1, nnz
48+
row_index = index(1,k)
49+
col_index = index(2,k)
50+
vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index)
51+
if( row_index==col_index ) cycle
52+
vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$row_index)
53+
end do
54+
55+
end if
56+
case(sparse_op_transpose)
57+
if(storage == sparse_full) then
58+
do k = 1, nnz
59+
col_index = index(1,k)
60+
row_index = index(2,k)
61+
vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index)
62+
end do
63+
64+
else
65+
do k = 1, nnz
66+
col_index = index(1,k)
67+
row_index = index(2,k)
68+
vec_y(${rksfx2(rank-1)}$row_index) = vec_y(${rksfx2(rank-1)}$row_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$col_index)
69+
if( row_index==col_index ) cycle
70+
vec_y(${rksfx2(rank-1)}$col_index) = vec_y(${rksfx2(rank-1)}$col_index) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$row_index)
71+
end do
72+
73+
end if
74+
#:if t1.startswith('complex')
75+
case(sparse_op_hermitian)
76+
if(storage == sparse_full) then
77+
do k = 1, nnz
78+
col_index = index(1,k)
79+
row_index = index(2,k)
80+
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)
81+
end do
82+
83+
else
84+
do k = 1, nnz
85+
col_index = index(1,k)
86+
row_index = index(2,k)
87+
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)
88+
if( row_index==col_index ) cycle
89+
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)
90+
end do
91+
92+
end if
93+
#:endif
94+
end select
95+
end associate
96+
end subroutine
97+
98+
#:endfor
99+
#:endfor
100+
101+
end submodule stdlib_sparse_spmv_coo

‎src/sparse/stdlib_sparse_spmv_csc.fypp‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
#:include "common.fypp"
2+
#:set RANKS = range(1, 2+1)
3+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
4+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
5+
#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
6+
#! define ranks without parentheses
7+
#:def rksfx2(rank)
8+
#{if rank > 0}#${":," + ":," * (rank - 1)}$#{endif}#
9+
#:enddef
10+
submodule (stdlib_sparse_spmv) stdlib_sparse_spmv_csc
11+
contains
12+
13+
!! spmv_csc
14+
#:for k1, t1, s1 in (KINDS_TYPES)
15+
#:for rank in RANKS
16+
module subroutine spmv_csc_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
17+
type(CSC_${s1}$_type), intent(in) :: matrix
18+
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
19+
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
20+
${t1}$, intent(in), optional :: alpha
21+
${t1}$, intent(in), optional :: beta
22+
character(1), intent(in), optional :: op
23+
${t1}$ :: alpha_
24+
character(1) :: op_
25+
integer(ilp) :: i, j
26+
#:if rank == 1
27+
${t1}$ :: aux
28+
#:else
29+
${t1}$ :: aux(size(vec_x,dim=1))
30+
#:endif
31+
32+
op_ = sparse_op_none; if(present(op)) op_ = op
33+
alpha_ = one_${k1}$
34+
if(present(alpha)) alpha_ = alpha
35+
if(present(beta)) then
36+
vec_y = beta * vec_y
37+
else
38+
vec_y = zero_${s1}$
39+
endif
40+
41+
associate( data => matrix%data, colptr => matrix%colptr, row => matrix%row, &
42+
& nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage )
43+
if( storage == sparse_full .and. op_==sparse_op_none ) then
44+
do concurrent(j=1:ncols)
45+
aux = alpha_ * vec_x(${rksfx2(rank-1)}$j)
46+
do i = colptr(j), colptr(j+1)-1
47+
vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + data(i) * aux
48+
end do
49+
end do
50+
51+
else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
52+
do concurrent(j=1:ncols)
53+
aux = zero_${k1}$
54+
do i = colptr(j), colptr(j+1)-1
55+
aux = aux + data(i) * vec_x(${rksfx2(rank-1)}$row(i))
56+
end do
57+
vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux
58+
end do
59+
60+
else if( storage == sparse_lower .and. op_/=sparse_op_hermitian )then
61+
do j = 1 , ncols
62+
aux = vec_x(${rksfx2(rank-1)}$j) * data(colptr(j))
63+
do i = colptr(j)+1, colptr(j+1)-1
64+
aux = aux + data(i) * vec_x(${rksfx2(rank-1)}$row(i))
65+
vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * data(i) * vec_x(${rksfx2(rank-1)}$j)
66+
end do
67+
vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux
68+
end do
69+
70+
else if( storage == sparse_upper .and. op_/=sparse_op_hermitian )then
71+
do j = 1 , ncols
72+
aux = zero_${s1}$
73+
do i = colptr(j), colptr(j+1)-2
74+
aux = aux + data(i) * vec_x(${rksfx2(rank-1)}$row(i))
75+
vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * data(i) * vec_x(${rksfx2(rank-1)}$j)
76+
end do
77+
aux = aux + data(colptr(j+1)-1) * vec_x(${rksfx2(rank-1)}$j)
78+
vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux
79+
end do
80+
81+
#:if t1.startswith('complex')
82+
else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then
83+
do concurrent(j=1:ncols)
84+
aux = zero_${k1}$
85+
do i = colptr(j), colptr(j+1)-1
86+
aux = aux + conjg(data(i)) * vec_x(${rksfx2(rank-1)}$row(i))
87+
end do
88+
vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux
89+
end do
90+
91+
else if( storage == sparse_lower .and. op_==sparse_op_hermitian )then
92+
do j = 1 , ncols
93+
aux = vec_x(${rksfx2(rank-1)}$j) * conjg(data(colptr(j)))
94+
do i = colptr(j)+1, colptr(j+1)-1
95+
aux = aux + conjg(data(i)) * vec_x(${rksfx2(rank-1)}$row(i))
96+
vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j)
97+
end do
98+
vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux
99+
end do
100+
101+
else if( storage == sparse_upper .and. op_==sparse_op_hermitian )then
102+
do j = 1 , ncols
103+
aux = zero_${s1}$
104+
do i = colptr(j), colptr(j+1)-2
105+
aux = aux + conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j)
106+
vec_y(${rksfx2(rank-1)}$row(i)) = vec_y(${rksfx2(rank-1)}$row(i)) + alpha_ * conjg(data(i)) * vec_x(${rksfx2(rank-1)}$j)
107+
end do
108+
aux = aux + conjg(data(colptr(j+1)-1)) * vec_x(${rksfx2(rank-1)}$j)
109+
vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_ * aux
110+
end do
111+
#:endif
112+
end if
113+
end associate
114+
end subroutine
115+
116+
#:endfor
117+
#:endfor
118+
119+
end submodule stdlib_sparse_spmv_csc

‎src/sparse/stdlib_sparse_spmv_csr.fypp‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,123 @@
1+
#:include "common.fypp"
2+
#:set RANKS = range(1, 2+1)
3+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
4+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
5+
#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
6+
#! define ranks without parentheses
7+
#:def rksfx2(rank)
8+
#{if rank > 0}#${":," + ":," * (rank - 1)}$#{endif}#
9+
#:enddef
10+
submodule (stdlib_sparse_spmv) stdlib_sparse_spmv_csr
11+
contains
12+
13+
!! spmv_csr
14+
#:for k1, t1, s1 in (KINDS_TYPES)
15+
#:for rank in RANKS
16+
module subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
17+
type(CSR_${s1}$_type), intent(in) :: matrix
18+
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
19+
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
20+
${t1}$, intent(in), optional :: alpha
21+
${t1}$, intent(in), optional :: beta
22+
character(1), intent(in), optional :: op
23+
${t1}$ :: alpha_
24+
character(1) :: op_
25+
integer(ilp) :: i, j
26+
#:if rank == 1
27+
${t1}$ :: aux, aux2
28+
#:else
29+
${t1}$ :: aux(size(vec_x,dim=1)), aux2(size(vec_x,dim=1))
30+
#:endif
31+
32+
op_ = sparse_op_none; if(present(op)) op_ = op
33+
alpha_ = one_${k1}$
34+
if(present(alpha)) alpha_ = alpha
35+
if(present(beta)) then
36+
vec_y = beta * vec_y
37+
else
38+
vec_y = zero_${s1}$
39+
endif
40+
41+
associate( data => matrix%data, col => matrix%col, rowptr => matrix%rowptr, &
42+
& nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage )
43+
44+
if( storage == sparse_full .and. op_==sparse_op_none ) then
45+
do i = 1, nrows
46+
aux = zero_${k1}$
47+
do j = rowptr(i), rowptr(i+1)-1
48+
aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j))
49+
end do
50+
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
51+
end do
52+
53+
else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
54+
do i = 1, nrows
55+
aux = alpha_ * vec_x(${rksfx2(rank-1)}$i)
56+
do j = rowptr(i), rowptr(i+1)-1
57+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux
58+
end do
59+
end do
60+
61+
else if( storage == sparse_lower .and. op_/=sparse_op_hermitian )then
62+
do i = 1 , nrows
63+
aux = zero_${s1}$
64+
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
65+
do j = rowptr(i), rowptr(i+1)-2
66+
aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j))
67+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux2
68+
end do
69+
aux = alpha_ * aux + data(j) * aux2
70+
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + aux
71+
end do
72+
73+
else if( storage == sparse_upper .and. op_/=sparse_op_hermitian )then
74+
do i = 1 , nrows
75+
aux = vec_x(${rksfx2(rank-1)}$i) * data(rowptr(i))
76+
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
77+
do j = rowptr(i)+1, rowptr(i+1)-1
78+
aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j))
79+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux2
80+
end do
81+
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
82+
end do
83+
84+
#:if t1.startswith('complex')
85+
else if( storage == sparse_full .and. op_==sparse_op_hermitian) then
86+
do i = 1, nrows
87+
aux = alpha_ * vec_x(${rksfx2(rank-1)}$i)
88+
do j = rowptr(i), rowptr(i+1)-1
89+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux
90+
end do
91+
end do
92+
93+
else if( storage == sparse_lower .and. op_==sparse_op_hermitian )then
94+
do i = 1 , nrows
95+
aux = zero_${s1}$
96+
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
97+
do j = rowptr(i), rowptr(i+1)-2
98+
aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j))
99+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2
100+
end do
101+
aux = alpha_ * aux + conjg(data(j)) * aux2
102+
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + aux
103+
end do
104+
105+
else if( storage == sparse_upper .and. op_==sparse_op_hermitian )then
106+
do i = 1 , nrows
107+
aux = vec_x(${rksfx2(rank-1)}$i) * conjg(data(rowptr(i)))
108+
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
109+
do j = rowptr(i)+1, rowptr(i+1)-1
110+
aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j))
111+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2
112+
end do
113+
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
114+
end do
115+
#:endif
116+
end if
117+
end associate
118+
end subroutine
119+
120+
#:endfor
121+
#:endfor
122+
123+
end submodule stdlib_sparse_spmv_csr

0 commit comments

Comments
 (0)