#: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