#:include "common.fypp"
#:set RANKS = range(1, 2+1)
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
#:set KINDS_TYPES = R_KINDS_TYPES+C_KINDS_TYPES
#! define ranks without parentheses
#:def rksfx2(rank)
#{if rank > 0}#${":," + ":," * (rank - 1)}$#{endif}#
#:enddef
submodule (stdlib_sparse_spmv) stdlib_sparse_spmv_coo
contains
!! spmv_coo
#:for k1, t1, s1 in (KINDS_TYPES)
#:for rank in RANKS
module subroutine spmv_coo_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
type(COO_${s1}$_type), intent(in) :: matrix
${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_, 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%data, matrix%index, matrix%storage, &
vec_x, beta_, vec_y)
end subroutine
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, intent(in) :: storage !! storage
${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
integer(ilp) :: col_index, k, row_index
integer(ilp) :: nnz !! number of non-zero values
nnz = size(index, 2)
vec_y = beta * vec_y
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)
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)
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)
end do
end if
case(sparse_op_transpose)
if(storage == sparse_full) then
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)
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)
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)
end do
end if
#:if t1.startswith('complex')
case(sparse_op_hermitian)
if(storage == sparse_full) then
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)
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)
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)
end do
end if
#:endif
end select
end subroutine
#:endfor
#:endfor
end submodule stdlib_sparse_spmv_coo
