diff --git a/doc/specs/stdlib_sparse.md b/doc/specs/stdlib_sparse.md
index 7ce53008f..369c2b792 100644
--- a/doc/specs/stdlib_sparse.md
+++ b/doc/specs/stdlib_sparse.md
@@ -94,7 +94,7 @@ CSR%rowptr(:) = [1,3,5,8,11]
 Experimental
 
 #### Description
-The Compressed Sparse Colum `CSC` is similar to the `CSR` format but values are accesed first by column, thus an index counter is given by `colptr` which enables to know the first and last non-zero row index of a given colum. 
+The Compressed Sparse Colum `CSC` is similar to the `CSR` format but values are accessed first by column, thus an index counter is given by `colptr` which enables to know the first and last non-zero row index of a given colum.
 
 ```Fortran
 type(CSC_sp_type) :: CSC
@@ -217,6 +217,81 @@ $$y=\alpha*op(M)*x+\beta*y$$
 
 `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`
 
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+## `spmv_kernel` - Non-object-oriented sparse matrix-vector product
+
+### Status
+
+Experimental
+
+### Description
+
+Provide non-object-oriented sparse matrix-vector product kernels for the current supported sparse matrix types.
+
+$$y=\alpha*op(M)*x+\beta*y$$
+
+### Syntax
+
+`call ` [[stdlib_sparse_spmv(module):spmv_kernel_coo(interface)]] `(op,alpha,data,index,storage,vec_x,beta,vec_y)`
+`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csc(interface)]] `(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y)`
+`call ` [[stdlib_sparse_spmv(module):spmv_kernel_csr(interface)]] `(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y)`
+`call ` [[stdlib_sparse_spmv(module):spmv_kernel_ell(interface)]] `(op,alpha,data,index,storage,vec_x,beta,vec_y)`
+`call ` [[stdlib_sparse_spmv(module):spmv_kernel_sellc(interface)]] `(op,alpha,data,ia,ja,storage,vec_x,beta,vec_y)`
+
+### Arguments
+
+Common arguments for all formats
+
+`op`: 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`
+
+`alpha`: Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument.
+
+`storage`: Shall be a scalar of `integer` type. It is an `intent(in)` argument. It defines the symmetry storage of the sparse matrix and must be one of `sparse_full`, `sparse_lower`, or `sparse_upper`.
+
+`vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument.
+
+`beta`: Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument.
+
+`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument.
+
+For the `COO` format:
+
+`data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument.
+
+`index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+For the `CSC` format:
+
+`data`: Shall be a rank-1 array of `real` or `complex` type. It is an `intent(in)` argument.
+
+`colptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+`row`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+
+For the `CSR` format:
+
+`data`: Shall be a rank-1 array of `real` or `complex` type array. It is an `intent(in)` argument.
+
+`col`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+`rowptr`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+For the `ELL` format:
+
+`data`: Shall be a rank-2 array of `real` or `complex` type. It is an `intent(in)` argument.
+
+`index`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+For the `SELLC` format:
+
+`data`: Shall be a rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument.
+
+`ia`: Shall be a rank-1 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+`ja`: Shall be a rank-2 array of `integer(ilp)` type. It is an `intent(in)` argument.
+
+
 <!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
 ## Sparse matrix to matrix conversions
 
@@ -411,4 +486,4 @@ The definition of all standard arithmetic operators have been overloaded to be a
 
 `B = A / alpha` or `B = alpha / A`
 
-*Note*: scalar addition and subtraction operators perform element-wise operations only on the stored (non-zero) values, not on the full mathematical matrix. Meaning, the sparsity pattern is preserved.
\ No newline at end of file
+*Note*: scalar addition and subtraction operators perform element-wise operations only on the stored (non-zero) values, not on the full mathematical matrix. Meaning, the sparsity pattern is preserved.
diff --git a/src/sparse/stdlib_sparse_spmv.fypp b/src/sparse/stdlib_sparse_spmv.fypp
index 5524ef7c0..c21054024 100644
--- a/src/sparse/stdlib_sparse_spmv.fypp
+++ b/src/sparse/stdlib_sparse_spmv.fypp
@@ -70,6 +70,119 @@ module stdlib_sparse_spmv
         end subroutine
         #:endfor
     end interface
+
+    !! Version experimental
+    !!
+    !! Apply the COO sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
+    !! [Specifications](../page/specs/stdlib_sparse.html#spmv)
+    interface spmv_kernel_coo
+        #:for k1, t1, s1 in (KINDS_TYPES)
+        #:for rank in RANKS
+        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(:,:)
+            integer, intent(in) :: 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
+        end subroutine
+        #:endfor
+        #:endfor
+    end interface
+
+    !! Version experimental
+    !!
+    !! Apply the CSC sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
+    !! [Specifications](../page/specs/stdlib_sparse.html#spmv)
+    interface spmv_kernel_csc
+        #:for k1, t1, s1 in (KINDS_TYPES)
+        #:for rank in RANKS
+        module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y)
+            ${t1}$, intent(in), contiguous :: data(:)
+            integer(ilp), intent(in), contiguous :: colptr(:) !! matrix column pointer
+            integer(ilp), intent(in), contiguous :: row(:)  !! matrix row pointer
+            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
+        end subroutine
+        #:endfor
+        #:endfor
+    end interface
+
+    !! Version experimental
+    !!
+    !! Apply the CSR sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
+    !! [Specifications](../page/specs/stdlib_sparse.html#spmv)
+    interface spmv_kernel_csr
+        #:for k1, t1, s1 in (KINDS_TYPES)
+        #:for rank in RANKS
+        module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y)
+            ${t1}$, intent(in), contiguous :: data(:)
+            integer(ilp), intent(in), contiguous :: col(:) !! matrix column pointer
+            integer(ilp), intent(in), contiguous :: rowptr(:)  !! matrix row pointer
+            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
+        end subroutine
+        #:endfor
+        #:endfor
+    end interface
+
+    !! Version experimental
+    !!
+    !! Apply the ELL sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
+    !! [Specifications](../page/specs/stdlib_sparse.html#spmv)
+    interface spmv_kernel_ell
+        #:for k1, t1, s1 in (KINDS_TYPES)
+        #:for rank in RANKS
+        module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,storage,vec_x,beta,vec_y)
+            ${t1}$, intent(in), contiguous :: data(:,:)
+            integer(ilp), intent(in), contiguous :: index(:,:)
+            integer, intent(in) :: 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
+        end subroutine
+        #:endfor
+        #:endfor
+    end interface
+
+    !! Version experimental
+    !!
+    !! Apply the SELLC sparse matrix-vector product $$y = \alpha * op(M) * x + \beta * y $$
+    !! [Specifications](../page/specs/stdlib_sparse.html#spmv)
+    interface spmv_kernel_sellc
+        #:for k1, t1, s1 in (KINDS_TYPES)
+        module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,storage,vec_x,beta,vec_y)
+            !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves
+            ${t1}$, intent(in), contiguous :: data(:,:)
+            integer(ilp), intent(in), contiguous :: ia(:)
+            integer(ilp), intent(in), contiguous :: ja(:,:)
+            integer, intent(in) :: storage
+            ${t1}$, intent(in), contiguous    :: vec_x(:)
+            ${t1}$, intent(inout), contiguous :: vec_y(:)
+            ${t1}$, intent(in) :: alpha
+            ${t1}$, intent(in) :: beta
+            character(1), intent(in) :: op
+        end subroutine
+        #:endfor
+    end interface
+
     public :: spmv
+    public :: spmv_kernel_coo
+    public :: spmv_kernel_csc
+    public :: spmv_kernel_csr
+    public :: spmv_kernel_ell
+    public :: spmv_kernel_sellc
 
 end module
diff --git a/src/sparse/stdlib_sparse_spmv_coo.fypp b/src/sparse/stdlib_sparse_spmv_coo.fypp
index 7aba6e9d8..e755a8ffc 100644
--- a/src/sparse/stdlib_sparse_spmv_coo.fypp
+++ b/src/sparse/stdlib_sparse_spmv_coo.fypp
@@ -20,36 +20,54 @@ contains
         ${t1}$, intent(in), optional :: alpha
         ${t1}$, intent(in), optional :: beta
         character(1), intent(in), optional :: op
-        ${t1}$ :: alpha_
+        ${t1}$ :: alpha_, beta_
         character(1) :: op_
-        integer(ilp) :: col_index, k, row_index
 
         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, index => matrix%index, storage => matrix%storage, nnz => matrix%nnz )
-        select case(op_)
+
+        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)
+                    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)
+                    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)
+                    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
@@ -58,16 +76,16 @@ contains
                 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)
+                    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)
+                    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)
+                    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
@@ -77,25 +95,24 @@ contains
                 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)
+                    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)
+                    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)
+                    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 associate
     end subroutine
 
     #:endfor
     #:endfor
 
-end submodule stdlib_sparse_spmv_coo
\ No newline at end of file
+end submodule stdlib_sparse_spmv_coo
diff --git a/src/sparse/stdlib_sparse_spmv_csc.fypp b/src/sparse/stdlib_sparse_spmv_csc.fypp
index 83b4ed41d..52ee6d3e8 100644
--- a/src/sparse/stdlib_sparse_spmv_csc.fypp
+++ b/src/sparse/stdlib_sparse_spmv_csc.fypp
@@ -20,100 +20,116 @@ contains
         ${t1}$, intent(in), optional :: alpha
         ${t1}$, intent(in), optional :: beta
         character(1), intent(in), optional :: op
-        ${t1}$ :: alpha_
+        ${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_${k1}$
+        if(present(beta)) beta_ = beta
+
+        call spmv_kernel_csc_${rank}$d_${s1}$(op_,alpha_, &
+            matrix%data,matrix%colptr,matrix%row,matrix%storage, &
+            vec_x,beta_,vec_y)
+
+    end subroutine
+    
+    module subroutine spmv_kernel_csc_${rank}$d_${s1}$(op,alpha,data,colptr,row,storage,vec_x,beta,vec_y)
+        ${t1}$, intent(in), contiguous :: data(:)
+        integer(ilp), intent(in), contiguous :: colptr(:) !! matrix column pointer
+        integer(ilp), intent(in), contiguous :: row(:)  !! matrix row pointer
+        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) :: i, j
+        integer(ilp) :: ncols
         #: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
+        ncols = size(colptr)-1
+
+        vec_y = beta * vec_y
 
-        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
+            if( storage == sparse_full .and. op==sparse_op_none ) then
                 do concurrent(j=1:ncols)
-                    aux = alpha_ * vec_x(${rksfx2(rank-1)}$j)
+                    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
+            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
+                    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
+            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)
+                        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
+                    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
+            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)
+                        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
+                    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
+            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
+                    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
+            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)
+                        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
+                    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
+            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)
+                        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
+                    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
\ No newline at end of file
+end submodule stdlib_sparse_spmv_csc
diff --git a/src/sparse/stdlib_sparse_spmv_csr.fypp b/src/sparse/stdlib_sparse_spmv_csr.fypp
index 929449969..cfe3dee61 100644
--- a/src/sparse/stdlib_sparse_spmv_csr.fypp
+++ b/src/sparse/stdlib_sparse_spmv_csr.fypp
@@ -20,104 +20,120 @@ contains
         ${t1}$, intent(in), optional :: alpha
         ${t1}$, intent(in), optional :: beta
         character(1), intent(in), optional :: op
-        ${t1}$ :: alpha_
+        ${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_${k1}$
+        if(present(beta)) beta_ = beta
+
+        call spmv_kernel_csr_${rank}$d_${s1}$(op_, alpha_, &
+            matrix%data, matrix%col, matrix%rowptr, matrix%storage, &
+            vec_x, beta_, vec_y)
+
+    end subroutine
+
+    module subroutine spmv_kernel_csr_${rank}$d_${s1}$(op,alpha,data,col,rowptr,storage,vec_x,beta,vec_y)
+        ${t1}$, intent(in), contiguous :: data(:)
+        integer(ilp), intent(in), contiguous :: col(:) !! matrix column pointer
+        integer(ilp), intent(in), contiguous :: rowptr(:)  !! matrix row pointer
+        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) :: i, j
+        integer(ilp) :: nrows
         #:if rank == 1
         ${t1}$ :: aux, aux2
         #:else
         ${t1}$ :: aux(size(vec_x,dim=1)), aux2(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
+        nrows = size(rowptr)-1
 
-        associate( data => matrix%data, col => matrix%col, rowptr => matrix%rowptr, &
-            & nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage )
-    
-            if( storage == sparse_full .and. op_==sparse_op_none ) then
+        vec_y = beta * vec_y
+   
+            if( storage == sparse_full .and. op==sparse_op_none ) then
                 do i = 1, nrows
                     aux = zero_${k1}$
                     do j = rowptr(i), rowptr(i+1)-1
                         aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j))
                     end do
-                    vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
+                    vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha * aux
                 end do
 
-            else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
+            else if( storage == sparse_full .and. op==sparse_op_transpose ) then
                 do i = 1, nrows
-                    aux = alpha_ * vec_x(${rksfx2(rank-1)}$i)
+                    aux = alpha * vec_x(${rksfx2(rank-1)}$i)
                     do j = rowptr(i), rowptr(i+1)-1
                         vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux
                     end do
                 end do
                 
-            else if( storage == sparse_lower .and. op_/=sparse_op_hermitian )then
+            else if( storage == sparse_lower .and. op/=sparse_op_hermitian )then
                 do i = 1 , nrows
                     aux  = zero_${s1}$
-                    aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
+                    aux2 = alpha * vec_x(${rksfx2(rank-1)}$i)
                     do j = rowptr(i), rowptr(i+1)-2
                         aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j))
                         vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux2
                     end do
-                    aux = alpha_ * aux + data(j) * aux2
+                    aux = alpha * aux + data(j) * aux2
                     vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + aux
                 end do
 
-            else if( storage == sparse_upper .and. op_/=sparse_op_hermitian )then
+            else if( storage == sparse_upper .and. op/=sparse_op_hermitian )then
                 do i = 1 , nrows
                     aux  = vec_x(${rksfx2(rank-1)}$i) * data(rowptr(i))
-                    aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
+                    aux2 = alpha * vec_x(${rksfx2(rank-1)}$i)
                     do j = rowptr(i)+1, rowptr(i+1)-1
                         aux = aux + data(j) * vec_x(${rksfx2(rank-1)}$col(j))
                         vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + data(j) * aux2
                     end do
-                    vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
+                    vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha * aux
                 end do
                 
             #:if t1.startswith('complex')
-            else if( storage == sparse_full .and. op_==sparse_op_hermitian) then
+            else if( storage == sparse_full .and. op==sparse_op_hermitian) then
                 do i = 1, nrows
-                    aux = alpha_ * vec_x(${rksfx2(rank-1)}$i)
+                    aux = alpha * vec_x(${rksfx2(rank-1)}$i)
                     do j = rowptr(i), rowptr(i+1)-1
                         vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux
                     end do
                 end do
 
-            else if( storage == sparse_lower .and. op_==sparse_op_hermitian )then
+            else if( storage == sparse_lower .and. op==sparse_op_hermitian )then
                 do i = 1 , nrows
                     aux  = zero_${s1}$
-                    aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
+                    aux2 = alpha * vec_x(${rksfx2(rank-1)}$i)
                     do j = rowptr(i), rowptr(i+1)-2
                         aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j))
                         vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2
                     end do
-                    aux = alpha_ * aux + conjg(data(j)) * aux2
+                    aux = alpha * aux + conjg(data(j)) * aux2
                     vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + aux
                 end do
 
-            else if( storage == sparse_upper .and. op_==sparse_op_hermitian )then
+            else if( storage == sparse_upper .and. op==sparse_op_hermitian )then
                 do i = 1 , nrows
                     aux  = vec_x(${rksfx2(rank-1)}$i) * conjg(data(rowptr(i)))
-                    aux2 = alpha_ * vec_x(${rksfx2(rank-1)}$i)
+                    aux2 = alpha * vec_x(${rksfx2(rank-1)}$i)
                     do j = rowptr(i)+1, rowptr(i+1)-1
                         aux = aux + conjg(data(j)) * vec_x(${rksfx2(rank-1)}$col(j))
                         vec_y(${rksfx2(rank-1)}$col(j)) = vec_y(${rksfx2(rank-1)}$col(j)) + conjg(data(j)) * aux2
                     end do
-                    vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_ * aux
+                    vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha * aux
                 end do
             #:endif
             end if
-        end associate
     end subroutine
-    
+
     #:endfor
     #:endfor
 
-end submodule stdlib_sparse_spmv_csr
\ No newline at end of file
+end submodule stdlib_sparse_spmv_csr
diff --git a/src/sparse/stdlib_sparse_spmv_ell.fypp b/src/sparse/stdlib_sparse_spmv_ell.fypp
index 243d1cbd4..9f268c155 100644
--- a/src/sparse/stdlib_sparse_spmv_ell.fypp
+++ b/src/sparse/stdlib_sparse_spmv_ell.fypp
@@ -20,68 +20,87 @@ contains
         ${t1}$, intent(in), optional :: alpha
         ${t1}$, intent(in), optional :: beta
         character(1), intent(in), optional :: op
-        ${t1}$ :: alpha_
+        ${t1}$ :: alpha_, beta_
         character(1) :: op_
-        integer(ilp) :: i, j, k
-        
+
         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, index => matrix%index, MNZ_P_ROW => matrix%K, &
-            & nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage )
-            if( storage == sparse_full .and. op_==sparse_op_none ) then
+
+        beta_ = zero_${k1}$
+        if(present(beta)) beta_ = beta
+
+        call spmv_kernel_ell_${rank}$d_${s1}$(op_,alpha_, &
+            matrix%data,matrix%index,matrix%storage, &
+            vec_x,beta_,vec_y)
+
+    end subroutine
+    
+    module subroutine spmv_kernel_ell_${rank}$d_${s1}$(op,alpha,data,index,storage,vec_x,beta,vec_y)
+        ${t1}$, intent(in), contiguous :: data(:,:)
+        integer(ilp), intent(in), contiguous :: index(:,:)
+        integer, intent(in) :: 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) :: i, j, k
+        integer(ilp) :: nrows
+        integer :: mnz_p_row
+        
+        nrows = size(index, 1)
+        mnz_p_row = size(index, 2)
+
+        vec_y = beta * vec_y
+
+            if( storage == sparse_full .and. op==sparse_op_none ) then
                 do i = 1, nrows
                     do k = 1, MNZ_P_ROW
                         j = index(i,k)
-                        if(j>0) vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$j)
+                        if(j>0) vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$j)
                     end do
                 end do
-            else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
+            else if( storage == sparse_full .and. op==sparse_op_transpose ) then
                 do i = 1, nrows
                     do k = 1, MNZ_P_ROW
                         j = index(i,k)
-                        if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$i)
+                        if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$i)
                     end do
                 end do
-            else if( storage /= sparse_full .and. op_/=sparse_op_hermitian ) then
+            else if( storage /= sparse_full .and. op/=sparse_op_hermitian ) then
                 do i = 1, nrows
                     do k = 1, MNZ_P_ROW
                         j = index(i,k)
                         if(j<=0) cycle
-                        vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$j)
+                        vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$j)
                         if(i==j) cycle 
-                        vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*data(i,k) * vec_x(${rksfx2(rank-1)}$i)
+                        vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*data(i,k) * vec_x(${rksfx2(rank-1)}$i)
                     end do
                 end do
             #:if t1.startswith('complex')
-            else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then
+            else if( storage == sparse_full .and. op==sparse_op_hermitian ) then
                 do i = 1, nrows
                     do k = 1, MNZ_P_ROW
                         j = index(i,k)
-                        if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i)
+                        if(j>0) vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i)
                     end do
                 end do
-            else if( storage /= sparse_full .and. op_==sparse_op_hermitian ) then
+            else if( storage /= sparse_full .and. op==sparse_op_hermitian ) then
                 do i = 1, nrows
                     do k = 1, MNZ_P_ROW
                         j = index(i,k)
                         if(j<=0) cycle
-                        vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha_*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$j)
+                        vec_y(${rksfx2(rank-1)}$i) = vec_y(${rksfx2(rank-1)}$i) + alpha*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$j)
                         if(i==j) cycle 
-                        vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha_*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i)
+                        vec_y(${rksfx2(rank-1)}$j) = vec_y(${rksfx2(rank-1)}$j) + alpha*conjg(data(i,k)) * vec_x(${rksfx2(rank-1)}$i)
                     end do
                 end do
             #:endif
             end if
-        end associate
     end subroutine
-    
+
     #:endfor
     #:endfor
 
-end submodule stdlib_sparse_spmv_ell
\ No newline at end of file
+end submodule stdlib_sparse_spmv_ell
diff --git a/src/sparse/stdlib_sparse_spmv_sellc.fypp b/src/sparse/stdlib_sparse_spmv_sellc.fypp
index 52205debd..203336cd2 100644
--- a/src/sparse/stdlib_sparse_spmv_sellc.fypp
+++ b/src/sparse/stdlib_sparse_spmv_sellc.fypp
@@ -16,32 +16,55 @@ contains
         ${t1}$, intent(in), optional :: alpha
         ${t1}$, intent(in), optional :: beta
         character(1), intent(in), optional :: op
-        ${t1}$ :: alpha_
+        ${t1}$ :: alpha_, beta_
         character(1) :: op_
-        integer(ilp) :: i, nz, rowidx, num_chunks, rm
 
         op_ = sparse_op_none; if(present(op)) op_ = op
         alpha_ = one_${s1}$
         if(present(alpha)) alpha_ = alpha
-        if(present(beta)) then
-            vec_y = beta * vec_y
-        else 
-            vec_y = zero_${s1}$
-        endif
 
-        associate( data => matrix%data, ia => matrix%rowptr , ja => matrix%col, cs => matrix%chunk_size, &
-        &   nnz => matrix%nnz, nrows => matrix%nrows, ncols => matrix%ncols, storage => matrix%storage  )
+        beta_ = zero_${s1}$
+        if(present(beta)) beta_ = beta
+
+        call spmv_kernel_sellc_${s1}$(op_,alpha_, &
+            matrix%data, matrix%rowptr, matrix%col, &
+            matrix%storage, &
+            vec_x,beta_,vec_y)
+
+    end subroutine
+
+    module subroutine spmv_kernel_sellc_${s1}$(op,alpha,data,ia,ja,storage,vec_x,beta,vec_y)
+        !! This algorithm was gracefully provided by Ivan Privec and adapted by Jose Alves
+        ${t1}$, intent(in), contiguous :: data(:,:)
+        integer(ilp), intent(in), contiguous :: ia(:)
+        integer(ilp), intent(in), contiguous :: ja(:,:)
+        integer, intent(in) :: storage
+        ${t1}$, intent(in), contiguous    :: vec_x(:)
+        ${t1}$, intent(inout), contiguous :: vec_y(:)
+        ${t1}$, intent(in) :: alpha
+        ${t1}$, intent(in) :: beta
+        character(1), intent(in) :: op
+        integer(ilp) :: i, nz, rowidx, num_chunks, rm
+        integer(ilp) :: nrows
+        integer :: chunk_size
+
+        chunk_size = size(data, 1)
 
-        if( .not.any( ${CHUNKS}$ == cs ) ) then
+        if( .not.any( ${CHUNKS}$ == chunk_size ) ) then
             print *, "error: sellc chunk size not supported."
             return
         end if
 
-        num_chunks = nrows / cs
-        rm = nrows - num_chunks * cs
-        if( storage == sparse_full .and. op_==sparse_op_none ) then
+        vec_y = beta * vec_y
+
+        nrows = merge(size(vec_y), size(vec_x), op==sparse_op_none)
 
-            select case(cs)
+        num_chunks = nrows / chunk_size
+        rm = nrows - num_chunks * chunk_size
+
+        if( storage == sparse_full .and. op==sparse_op_none ) then
+
+            select case(chunk_size)
             #:for chunk in CHUNKS
             case(${chunk}$)
                 do i = 1, num_chunks
@@ -56,13 +79,13 @@ contains
             if(rm>0)then 
                 i = num_chunks + 1 
                 nz = ia(i+1) - ia(i)
-                rowidx = (i - 1)*cs + 1
-                call chunk_kernel_rm(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:))
+                rowidx = (i - 1)*chunk_size + 1
+                call chunk_kernel_rm(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:))
             end if
             
-        else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
+        else if( storage == sparse_full .and. op==sparse_op_transpose ) then
 
-            select case(cs)
+            select case(chunk_size)
             #:for chunk in CHUNKS
             case(${chunk}$)
                 do i = 1, num_chunks
@@ -77,14 +100,14 @@ contains
             if(rm>0)then 
                 i = num_chunks + 1 
                 nz = ia(i+1) - ia(i)
-                rowidx = (i - 1)*cs + 1
-                call chunk_kernel_rm_trans(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y)
+                rowidx = (i - 1)*chunk_size + 1
+                call chunk_kernel_rm_trans(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y)
             end if
         
         #:if t1.startswith('complex')
-        else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then
+        else if( storage == sparse_full .and. op==sparse_op_hermitian ) then
 
-            select case(cs)
+            select case(chunk_size)
             #:for chunk in CHUNKS
             case(${chunk}$)
                 do i = 1, num_chunks
@@ -99,15 +122,14 @@ contains
             if(rm>0)then 
                 i = num_chunks + 1 
                 nz = ia(i+1) - ia(i)
-                rowidx = (i - 1)*cs + 1
-                call chunk_kernel_rm_herm(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y)
+                rowidx = (i - 1)*chunk_size + 1
+                call chunk_kernel_rm_herm(nz,chunk_size,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y)
             end if
         #:endif
         else
             print *, "error: sellc format for spmv operation not yet supported."
             return
         end if
-        end associate
 
     contains
         #:for chunk in CHUNKS
@@ -118,7 +140,7 @@ contains
             ${t1}$, intent(inout) :: y(${chunk}$)
             integer :: j
             do j = 1, n
-                y(:) = y(:) + alpha_ * a(:,j) * x(col(:,j))
+                y(:) = y(:) + alpha * a(:,j) * x(col(:,j))
             end do
         end subroutine
         pure subroutine chunk_kernel_trans_${chunk}$(n,a,col,x,y)
@@ -129,7 +151,7 @@ contains
             integer :: j, k
             do j = 1, n
                 do k = 1, ${chunk}$
-                    y(col(k,j)) = y(col(k,j)) + alpha_ * a(k,j) * x(k)
+                    y(col(k,j)) = y(col(k,j)) + alpha * a(k,j) * x(k)
                 end do
             end do
         end subroutine
@@ -142,7 +164,7 @@ contains
             integer :: j, k
             do j = 1, n
                 do k = 1, ${chunk}$
-                    y(col(k,j)) = y(col(k,j)) + alpha_ * conjg(a(k,j)) * x(k)
+                    y(col(k,j)) = y(col(k,j)) + alpha * conjg(a(k,j)) * x(k)
                 end do
             end do
         end subroutine
@@ -156,7 +178,7 @@ contains
             ${t1}$, intent(inout) :: y(r)
             integer :: j
             do j = 1, n
-                y(1:r) = y(1:r) + alpha_ * a(1:r,j) * x(col(1:r,j))
+                y(1:r) = y(1:r) + alpha * a(1:r,j) * x(col(1:r,j))
             end do
         end subroutine
         pure subroutine chunk_kernel_rm_trans(n,cs,r,a,col,x,y)
@@ -167,7 +189,7 @@ contains
             integer :: j, k
             do j = 1, n
                 do k = 1, r
-                    y(col(k,j)) = y(col(k,j)) + alpha_ * a(k,j) * x(k)
+                    y(col(k,j)) = y(col(k,j)) + alpha * a(k,j) * x(k)
                 end do
             end do
         end subroutine
@@ -180,14 +202,13 @@ contains
             integer :: j, k
             do j = 1, n
                 do k = 1, r
-                    y(col(k,j)) = y(col(k,j)) + alpha_ * conjg(a(k,j)) * x(k)
+                    y(col(k,j)) = y(col(k,j)) + alpha * conjg(a(k,j)) * x(k)
                 end do
             end do
         end subroutine
         #:endif
 
     end subroutine
-    
     #:endfor
 
-end submodule stdlib_sparse_spmv_sellc
\ No newline at end of file
+end submodule stdlib_sparse_spmv_sellc
diff --git a/test/linalg/test_linalg_sparse.fypp b/test/linalg/test_linalg_sparse.fypp
index 0ebe235b5..d03c1cc5e 100644
--- a/test/linalg/test_linalg_sparse.fypp
+++ b/test/linalg/test_linalg_sparse.fypp
@@ -41,6 +41,8 @@ contains
         block
             integer, parameter :: wp = ${k1}$
             type(COO_${s1}$_type) :: COO
+            ${t1}$, parameter :: alpha = 1._wp
+            ${t1}$, parameter :: beta = 0._wp
             ${t1}$, allocatable :: dense(:,:)
             ${t1}$, allocatable :: vec_x(:)
             ${t1}$, allocatable :: vec_y1(:), vec_y2(:)
@@ -63,15 +65,32 @@ contains
             call check(error, all(vec_y1 == real([6,11,15,15],kind=wp)) )
             if (allocated(error)) return
 
+            ! High-level API
             call spmv( COO, vec_x, vec_y2 )
             call check(error, all(vec_y1 == vec_y2) )
             if (allocated(error)) return
 
+            ! Low-level API
+            call spmv_kernel_coo(sparse_op_none, alpha, &
+                COO%data ,COO%index, COO%storage, &
+                vec_x, beta, vec_y2 )
+            call check(error, all(vec_y1 == vec_y2), 'COO - low level API: wrong output' )
+            if (allocated(error)) return
+
             ! Test in-place transpose
             vec_y1 = 1._wp
+            ! High-level API
             call spmv( COO, vec_y1, vec_x, op=sparse_op_transpose )
             call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) )
             if (allocated(error)) return
+
+            ! Low-level API
+            call spmv_kernel_coo(sparse_op_transpose, alpha, &
+                COO%data ,COO%index, COO%storage, &
+                vec_y1, beta, vec_x)
+            call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'COO - transpose - low level API: wrong output' )
+            if (allocated(error)) return
+ 
         end block
         #:endfor
     end subroutine
@@ -111,6 +130,8 @@ contains
         block
             integer, parameter :: wp = ${k1}$
             type(CSR_${s1}$_type) :: CSR
+            ${t1}$, parameter :: alpha = 1._wp
+            ${t1}$, parameter :: beta = 0._wp
             ${t1}$, allocatable :: vec_x(:)
             ${t1}$, allocatable :: vec_y(:)
 
@@ -121,15 +142,31 @@ contains
             
             allocate( vec_x(5) , source = 1._wp )
             allocate( vec_y(4) , source = 0._wp )
+            ! High-level API
             call spmv( CSR, vec_x, vec_y )
             
-            call check(error, all(vec_y == real([6,11,15,15],kind=wp)) )
+            call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSR - high-level API - wrong output' )
+            if (allocated(error)) return
+
+            ! Low-level API
+            call spmv_kernel_csr(sparse_op_none, alpha, &
+                CSR%data, CSR%col, CSR%rowptr, CSR%storage, &
+                vec_x, beta, vec_y )
+            
+            call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSR - low-level API - wrong output' )
             if (allocated(error)) return
 
             ! Test in-place transpose
             vec_y = 1._wp
+            ! High-level API
             call spmv( CSR, vec_y, vec_x, op=sparse_op_transpose )
-            call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) )
+            call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'CSR transpose - high-level API - wrong output' )
+            if (allocated(error)) return
+
+            call spmv_kernel_csr(sparse_op_transpose, alpha, &
+                CSR%data, CSR%col, CSR%rowptr, CSR%storage, &
+                vec_y, beta, vec_x)
+            call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'CSR transpose - low-level API - wrong output' )
             if (allocated(error)) return
         end block
         #:endfor
@@ -142,6 +179,8 @@ contains
         block
             integer, parameter :: wp = ${k1}$
             type(CSC_${s1}$_type) :: CSC
+            ${t1}$, parameter :: alpha = 1._wp
+            ${t1}$, parameter :: beta = 0._wp
             ${t1}$, allocatable :: vec_x(:)
             ${t1}$, allocatable :: vec_y(:)
 
@@ -152,16 +191,34 @@ contains
             
             allocate( vec_x(5) , source = 1._wp )
             allocate( vec_y(4) , source = 0._wp )
+            !High-level API
             call spmv( CSC, vec_x, vec_y )
             
             call check(error, all(vec_y == real([6,11,15,15],kind=wp)) )
             if (allocated(error)) return
 
+            !High-level API
+            call spmv_kernel_csc(sparse_op_none, alpha, &
+                CSC%data, CSC%colptr, CSC%row, CSC%storage, &
+                vec_x, beta, vec_y )
+            
+            call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'CSC - low-level API - wrong output' )
+            if (allocated(error)) return
+
             ! Test in-place transpose
             vec_y = 1._wp
+            !High-level API
             call spmv( CSC, vec_y, vec_x, op=sparse_op_transpose )
             call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) )
             if (allocated(error)) return
+
+            !Low-level API
+            call spmv_kernel_csc(sparse_op_transpose, alpha, &
+                CSC%data, CSC%colptr, CSC%row, CSC%storage, &
+                vec_y, beta, vec_x)
+            call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'CSC tranpose - low-level API - wrong output' )
+            if (allocated(error)) return
+
         end block
         #:endfor
     end subroutine
@@ -173,6 +230,8 @@ contains
         block
             integer, parameter :: wp = ${k1}$
             type(ELL_${s1}$_type) :: ELL
+            ${t1}$, parameter :: alpha = 1._wp
+            ${t1}$, parameter :: beta = 0._wp
             ${t1}$, allocatable :: vec_x(:)
             ${t1}$, allocatable :: vec_y(:)
 
@@ -189,16 +248,33 @@ contains
             
             allocate( vec_x(5) , source = 1._wp )
             allocate( vec_y(4) , source = 0._wp )
+            !High-level API
             call spmv( ELL, vec_x, vec_y )
             
             call check(error, all(vec_y == real([6,11,15,15],kind=wp)) )
             if (allocated(error)) return
 
+            !Low-level API
+            call spmv_kernel_ell(sparse_op_none, alpha, &
+                ELL%data, ELL%index, ELL%storage, &
+                vec_x, beta, vec_y )
+
+            call check(error, all(vec_y == real([6,11,15,15],kind=wp)), 'ELL - low-level API - wrong output' )
+            if (allocated(error)) return
+
             ! Test in-place transpose
             vec_y = 1._wp
+            !High-level API
             call spmv( ELL, vec_y, vec_x, op=sparse_op_transpose )
             call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)) )
             if (allocated(error)) return
+
+            !Low-level API
+            call spmv_kernel_ell(sparse_op_transpose, alpha, &
+                ELL%data, ELL%index, ELL%storage, &
+                vec_y, beta, vec_x)
+            call check(error, all(vec_x == real([17,15,4,14,-3],kind=wp)), 'ELL tranpose - low-level API - wrong output' )
+            if (allocated(error)) return
         end block
         #:endfor
         
@@ -212,6 +288,8 @@ contains
             integer, parameter :: wp = ${k1}$
             type(SELLC_${s1}$_type) :: SELLC
             type(CSR_${s1}$_type)   :: CSR
+            ${t1}$, parameter :: alpha = 1._wp
+            ${t1}$, parameter :: beta = 0._wp
             ${t1}$, allocatable :: vec_x(:)
             ${t1}$, allocatable :: vec_y(:), vec_y2(:)
             integer :: i
@@ -232,18 +310,36 @@ contains
             allocate( vec_x(6) , source = 1._wp )
             allocate( vec_y(6) , source = 0._wp )
             
+            !High-level API
             call spmv( SELLC, vec_x, vec_y )
             
-            call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)) )
+            call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)), 'SELLC - high-level API - wrong output' )
+            if (allocated(error)) return
+
+            !Low-level API
+            call spmv_kernel_sellc(sparse_op_none, alpha, &
+                SELLC%data, SELLC%rowptr, SELLC%col, SELLC%storage, &
+                vec_x, beta, vec_y )
+
+            call check(error, all(vec_y == real([6,22,27,23,27,48],kind=wp)), 'SELLC - low-level API - wrong output' )
             if (allocated(error)) return
 
             ! Test in-place transpose
             vec_x = real( [1,2,3,4,5,6] , kind=wp )
             call spmv( CSR, vec_x, vec_y , op = sparse_op_transpose )
             allocate( vec_y2(6) , source = 0._wp )
+            !High-level API
             call spmv( SELLC, vec_x, vec_y2 , op = sparse_op_transpose )
-            
-            call check(error, all(vec_y == vec_y2))
+
+            call check(error, all(vec_y == vec_y2), 'SELLC tranpose - high-level API - wrong output' )
+            if (allocated(error)) return
+
+            !Low-level API
+            call spmv_kernel_sellc(sparse_op_transpose, alpha, &
+                SELLC%data, SELLC%rowptr, SELLC%col, SELLC%storage, &
+                vec_x, beta, vec_y2)
+
+            call check(error, all(vec_y == vec_y2), 'SELLC tranpose - low-level API - wrong output' )
             if (allocated(error)) return
 
         end block
@@ -715,4 +811,4 @@ program tester
         write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
         error stop
     end if
-end program
\ No newline at end of file
+end program