spmv: low-level API based on existing high-level API

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.
+*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.

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

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
+end submodule stdlib_sparse_spmv_coo