#: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_csc
contains
!! spmv_csc
#:for k1, t1, s1 in (KINDS_TYPES)
#:for rank in RANKS
module subroutine spmv_csc_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
type(CSC_${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_
character(1) :: op_
integer(ilp) :: i, j
#:if rank == 1
${t1}$ :: aux
#:else
${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
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)
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
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
end do
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)
end do
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
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)
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
end do
#:if t1.startswith('complex')
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
end do
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)
end do
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
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)
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
end do
#:endif
end if
end associate
end subroutine
#:endfor
#:endfor
end submodule stdlib_sparse_spmv_csc
