Skip to content

Commit a4d9306

Browse files
jalveszjalvesz
committed
Add in place operator for coo and csr spmv
1 parent 7279461 commit a4d9306

File tree

4 files changed

+140
-15
lines changed

4 files changed

+140
-15
lines changed

‎doc/specs/stdlib_sparse.md‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -194,11 +194,11 @@ Experimental
194194

195195
Provide sparse matrix-vector product kernels for the current supported sparse matrix types.
196196

197-
$$y=\alpha*M*x+\beta*y$$
197+
$$y=\alpha*op(M)*x+\beta*y$$
198198

199199
### Syntax
200200

201-
`call ` [[stdlib_sparse_spmv(module):spmv(interface)]] `(matrix,vec_x,vec_y [,alpha,beta])`
201+
`call ` [[stdlib_sparse_spmv(module):spmv(interface)]] `(matrix,vec_x,vec_y [,alpha,beta,op])`
202202

203203
### Arguments
204204

@@ -212,6 +212,8 @@ $$y=\alpha*M*x+\beta*y$$
212212

213213
`beta`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument.
214214

215+
`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`
216+
215217
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
216218
## Sparse matrix to matrix conversions
217219

‎src/stdlib_sparse_constants.fypp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,10 @@ module stdlib_sparse_constants
1212
enumerator :: sparse_lower !! Symmetric Sparse matrix with triangular inferior storage
1313
enumerator :: sparse_upper !! Symmetric Sparse matrix with triangular supperior storage
1414
end enum
15+
16+
character(1), parameter :: sparse_op_none = 'N' !! no transpose
17+
character(1), parameter :: sparse_op_transpose = 'T' !! transpose
18+
character(1), parameter :: sparse_op_hermitian = 'H' !! conjugate or hermitian transpose
1519

1620
! Integer size support for ILP64 builds should be done here
1721
integer, parameter :: ilp = int32

‎src/stdlib_sparse_kinds.fypp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ module stdlib_sparse_kinds
1212
implicit none
1313
private
1414
public :: sparse_full, sparse_lower, sparse_upper
15+
public :: sparse_op_none, sparse_op_transpose, sparse_op_hermitian
1516
!! version: experimental
1617
!!
1718
!! Base sparse type holding the meta data related to the storage capacity of a matrix.

‎src/stdlib_sparse_spmv.fypp‎

Lines changed: 131 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ module stdlib_sparse_spmv
1818

1919
!! Version experimental
2020
!!
21-
!! Applay the sparse matrix-vector product $$y = \alpha * M * x + \beta * y $$
21+
!! Applay the sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
2222
!! [Specifications](../page/specs/stdlib_sparse.html#spmv)
2323
interface spmv
2424
#:for k1, t1, s1 in (KINDS_TYPES)
@@ -38,15 +38,18 @@ contains
3838
!! spmv_coo
3939
#:for k1, t1, s1 in (KINDS_TYPES)
4040
#:for rank in RANKS
41-
subroutine spmv_coo_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta)
41+
subroutine spmv_coo_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
4242
type(COO_${s1}$_type), intent(in) :: matrix
4343
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
4444
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
4545
${t1}$, intent(in), optional :: alpha
4646
${t1}$, intent(in), optional :: beta
47+
character(1), intent(in), optional :: op
4748
${t1}$ :: alpha_, beta_
49+
character(1) :: op_
4850
integer(ilp) :: k, ik, jk
4951

52+
op_ = sparse_op_none; if(present(op)) op_ = op
5053
alpha_ = one_${k1}$
5154
if(present(alpha)) alpha_ = alpha
5255
if(present(beta)) then
@@ -55,7 +58,9 @@ contains
5558
vec_y = zero_${s1}$
5659
endif
5760
associate( data => matrix%data, index => matrix%index, storage => matrix%storage, nnz => matrix%nnz )
58-
if( storage == sparse_full) then
61+
select case(op_)
62+
case(sparse_op_none)
63+
if(storage == sparse_full) then
5964
do concurrent (k = 1:nnz)
6065
ik = index(1,k)
6166
jk = index(2,k)
@@ -72,6 +77,45 @@ contains
7277
end do
7378

7479
end if
80+
case(sparse_op_transpose)
81+
if(storage == sparse_full) then
82+
do concurrent (k = 1:nnz)
83+
jk = index(1,k)
84+
ik = index(2,k)
85+
vec_y(${rksfx2(rank-1)}$ik) = vec_y(${rksfx2(rank-1)}$ik) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$jk)
86+
end do
87+
88+
else
89+
do concurrent (k = 1:nnz)
90+
jk = index(1,k)
91+
ik = index(2,k)
92+
vec_y(${rksfx2(rank-1)}$ik) = vec_y(${rksfx2(rank-1)}$ik) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$jk)
93+
if( ik==jk ) cycle
94+
vec_y(${rksfx2(rank-1)}$jk) = vec_y(${rksfx2(rank-1)}$jk) + alpha_*data(k) * vec_x(${rksfx2(rank-1)}$ik)
95+
end do
96+
97+
end if
98+
#:if t1.startswith('complex')
99+
case(sparse_op_hermitian)
100+
if(storage == sparse_full) then
101+
do concurrent (k = 1:nnz)
102+
jk = index(1,k)
103+
ik = index(2,k)
104+
vec_y(${rksfx2(rank-1)}$ik) = vec_y(${rksfx2(rank-1)}$ik) + alpha_*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$jk)
105+
end do
106+
107+
else
108+
do concurrent (k = 1:nnz)
109+
jk = index(1,k)
110+
ik = index(2,k)
111+
vec_y(${rksfx2(rank-1)}$ik) = vec_y(${rksfx2(rank-1)}$ik) + alpha_*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$jk)
112+
if( ik==jk ) cycle
113+
vec_y(${rksfx2(rank-1)}$jk) = vec_y(${rksfx2(rank-1)}$jk) + alpha_*conjg(data(k)) * vec_x(${rksfx2(rank-1)}$ik)
114+
end do
115+
116+
end if
117+
#:endif
118+
end select
75119
end associate
76120
end subroutine
77121

@@ -81,28 +125,32 @@ contains
81125
!! spmv_csr
82126
#:for k1, t1, s1 in (KINDS_TYPES)
83127
#:for rank in RANKS
84-
subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta)
128+
subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
85129
type(CSR_${s1}$_type), intent(in) :: matrix
86130
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
87131
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
88132
${t1}$, intent(in), optional :: alpha
89133
${t1}$, intent(in), optional :: beta
134+
character(1), intent(in), optional :: op
90135
${t1}$ :: alpha_, beta_
136+
character(1) :: op_
91137
integer(ilp) :: i, j
92138
#:if rank == 1
93139
${t1}$ :: aux, aux2
94140
#:else
95141
${t1}$ :: aux(size(vec_x,dim=1)), aux2(size(vec_x,dim=1))
96142
#:endif
97143

144+
op_ = sparse_op_none; if(present(op)) op_ = op
98145
alpha_ = one_${k1}$
99146
if(present(alpha)) alpha_ = alpha
100147
beta_ = zero_${k1}$
101148
if(present(beta)) beta_ = beta
102149

103150
associate( data => matrix%data, col => matrix%col, rowptr => matrix%rowptr, &
104151
& nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage )
105-
if( storage == sparse_full) then
152+
153+
if( storage == sparse_full .and. op_==sparse_op_none ) then
106154
do i = 1, nrows
107155
aux = zero_${k1}$
108156
do j = rowptr(i), rowptr(i+1)-1
@@ -114,8 +162,21 @@ contains
114162
vec_y(${rksfx2(rank-1)}$i) = alpha_ * aux
115163
end if
116164
end do
165+
166+
else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
167+
if(present(beta)) then
168+
vec_y = beta * vec_y
169+
else
170+
vec_y = zero_${s1}$
171+
endif
172+
do i = 1, nrows
173+
aux = alpha_ * vec_x(${rksfx2(rank-1)}$i)
174+
do j = rowptr(i), rowptr(i+1)-1
175+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux
176+
end do
177+
end do
117178

118-
else if( storage == sparse_lower )then
179+
else if( storage == sparse_lower .and. op_/=sparse_op_hermitian )then
119180
do i = 1 , nrows
120181
aux = zero_${s1}$
121182
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
@@ -132,7 +193,7 @@ contains
132193
end if
133194
end do
134195

135-
else if( storage == sparse_upper )then
196+
else if( storage == sparse_upper .and. op_/=sparse_op_hermitian )then
136197
do i = 1 , nrows
137198
aux = vec_x(${rksfx2(rank-1)}$i) * data(rowptr(i))
138199
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
@@ -150,7 +211,57 @@ contains
150211
end if
151212
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
152213
end do
214+
215+
#:if t1.startswith('complex')
216+
else if( storage == sparse_full .and. op_==sparse_op_hermitian) then
217+
if(present(beta)) then
218+
vec_y = beta * vec_y
219+
else
220+
vec_y = zero_${s1}$
221+
endif
222+
do i = 1, nrows
223+
aux = alpha_ * vec_x(${rksfx2(rank-1)}$i)
224+
do j = rowptr(i), rowptr(i+1)-1
225+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux
226+
end do
227+
end do
228+
229+
else if( storage == sparse_lower .and. op_==sparse_op_hermitian )then
230+
do i = 1 , nrows
231+
aux = zero_${s1}$
232+
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
233+
do j = rowptr(i), rowptr(i+1)-2
234+
aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j))
235+
vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2
236+
end do
237+
aux = alpha_ * aux + conjg(data(j)) * aux2
238+
239+
if(present(beta)) then
240+
vec_y(${rksfx2(rank-1)}$i) = beta_ * vec_y(${rksfx2(rank-1)}$i) + aux
241+
else
242+
vec_y(${rksfx2(rank-1)}$i) = aux
243+
end if
244+
end do
153245

246+
else if( storage == sparse_upper .and. op_==sparse_op_hermitian )then
247+
do i = 1 , nrows
248+
aux = vec_x(${rksfx2(rank-1)}$i) * conjg(data(rowptr(i)))
249+
aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
250+
do j = rowptr(i)+1, rowptr(i+1)-1
251+
aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j))
252+
end do
253+
if(present(beta)) then
254+
do j = rowptr(i)+1, rowptr(i+1)-1
255+
vec_y(${rksfx2(rank-1)}$col(j)) = beta_ * vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2
256+
end do
257+
else
258+
do j = rowptr(i)+1, rowptr(i+1)-1
259+
vec_y(${rksfx2(rank-1)}$col(j)) = conjg(data(j)) * aux2
260+
end do
261+
end if
262+
vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
263+
end do
264+
#:endif
154265
end if
155266
end associate
156267
end subroutine
@@ -161,20 +272,23 @@ contains
161272
!! spmv_csc
162273
#:for k1, t1, s1 in (KINDS_TYPES)
163274
#:for rank in RANKS
164-
subroutine spmv_csc_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta)
275+
subroutine spmv_csc_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
165276
type(CSC_${s1}$_type), intent(in) :: matrix
166277
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
167278
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
168279
${t1}$, intent(in), optional :: alpha
169280
${t1}$, intent(in), optional :: beta
281+
character(1), intent(in), optional :: op
170282
${t1}$ :: alpha_, beta_
283+
character(1) :: op_
171284
integer(ilp) :: i, j
172285
#:if rank == 1
173286
${t1}$ :: aux
174287
#:else
175288
${t1}$ :: aux(size(vec_x,dim=1))
176289
#:endif
177290

291+
op_ = sparse_op_none; if(present(op)) op_ = op
178292
alpha_ = one_${k1}$
179293
if(present(alpha)) alpha_ = alpha
180294
if(present(beta)) then
@@ -193,7 +307,6 @@ contains
193307
end do
194308

195309
else if( storage == sparse_lower )then
196-
! NOT TESTED
197310
do j = 1 , ncols
198311
aux = vec_x(${rksfx2(rank-1)}$j) * data(colptr(j))
199312
do i = colptr(j)+1, colptr(j+1)-1
@@ -204,7 +317,6 @@ contains
204317
end do
205318

206319
else if( storage == sparse_upper )then
207-
! NOT TESTED
208320
do j = 1 , ncols
209321
aux = zero_${s1}$
210322
do i = colptr(j), colptr(i+1)-2
@@ -225,15 +337,18 @@ contains
225337
!! spmv_ell
226338
#:for k1, t1, s1 in (KINDS_TYPES)
227339
#:for rank in RANKS
228-
subroutine spmv_ell_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta)
340+
subroutine spmv_ell_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
229341
type(ELL_${s1}$_type), intent(in) :: matrix
230342
${t1}$, intent(in) :: vec_x${ranksuffix(rank)}$
231343
${t1}$, intent(inout) :: vec_y${ranksuffix(rank)}$
232344
${t1}$, intent(in), optional :: alpha
233345
${t1}$, intent(in), optional :: beta
346+
character(1), intent(in), optional :: op
234347
${t1}$ :: alpha_, beta_
348+
character(1) :: op_
235349
integer(ilp) :: i, j, k
236-
350+
351+
op_ = sparse_op_none; if(present(op)) op_ = op
237352
alpha_ = one_${k1}$
238353
if(present(alpha)) alpha_ = alpha
239354
if(present(beta)) then
@@ -259,16 +374,19 @@ contains
259374
!! spmv_sellc
260375
#:set CHUNKS = [4,8,16]
261376
#:for k1, t1, s1 in (KINDS_TYPES)
262-
subroutine spmv_sellc_${s1}$(matrix,vec_x,vec_y,alpha,beta)
377+
subroutine spmv_sellc_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
263378
!! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves
264379
type(SELLC_${s1}$_type), intent(in) :: matrix
265380
${t1}$, intent(in) :: vec_x(:)
266381
${t1}$, intent(inout) :: vec_y(:)
267382
${t1}$, intent(in), optional :: alpha
268383
${t1}$, intent(in), optional :: beta
384+
character(1), intent(in), optional :: op
269385
${t1}$ :: alpha_, beta_
386+
character(1) :: op_
270387
integer(ilp) :: i, nz, rowidx, num_chunks, rm
271388

389+
op_ = sparse_op_none; if(present(op)) op_ = op
272390
alpha_ = one_${s1}$
273391
if(present(alpha)) alpha_ = alpha
274392
if(present(beta)) then

0 commit comments

Comments
 (0)