housekeeping (sparse): split spmv kernels into submodules

doc/specs/stdlib_sparse.md

@@ -146,7 +146,10 @@ Type-bound procedures to enable adding data in a sparse matrix.
 
 ### Syntax
 
-`call matrix%add(i,j,v)` or
+* Add single value
+`call matrix%add(i,j,v)` 
+
+* Add a block of values
 `call matrix%add(i(:),j(:),v(:,:))`
 
 ### Arguments

src/sparse/CMakeLists.txt

@@ -4,6 +4,11 @@ set(sparse_fppFiles
     stdlib_sparse_kinds.fypp
     stdlib_sparse_operators.fypp
     stdlib_sparse_spmv.fypp
+    stdlib_sparse_spmv_coo.fypp
+    stdlib_sparse_spmv_csr.fypp
+    stdlib_sparse_spmv_csc.fypp
+    stdlib_sparse_spmv_ell.fypp
+    stdlib_sparse_spmv_sellc.fypp
     )
 
 set(sparse_cppFiles

src/sparse/stdlib_sparse_spmv_coo.fypp

@@ -0,0 +1,101 @@
+#: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_
+        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_)
+        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 associate
+    end subroutine
+
+    #:endfor
+    #:endfor
+
+end submodule stdlib_sparse_spmv_coo

src/sparse/stdlib_sparse_spmv_csc.fypp

@@ -0,0 +1,119 @@
+#: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

src/sparse/stdlib_sparse_spmv_csr.fypp

@@ -0,0 +1,123 @@
+#: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_csr
+contains
+
+    !! spmv_csr
+    #:for k1, t1, s1 in (KINDS_TYPES)
+    #:for rank in RANKS
+    module subroutine spmv_csr_${rank}$d_${s1}$(matrix,vec_x,vec_y,alpha,beta,op)
+        type(CSR_${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, 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
+
+        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
+                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
+                end do
+
+            else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
+                do i = 1, nrows
+                    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
+                do i = 1 , nrows
+                    aux  = zero_${s1}$
+                    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
+                    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
+                do i = 1 , nrows
+                    aux  = vec_x(${rksfx2(rank-1)}$i) * data(rowptr(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
+                end do
+
+            #:if t1.startswith('complex')
+            else if( storage == sparse_full .and. op_==sparse_op_hermitian) then
+                do i = 1, nrows
+                    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
+                do i = 1 , nrows
+                    aux  = zero_${s1}$
+                    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
+                    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
+                do i = 1 , nrows
+                    aux  = vec_x(${rksfx2(rank-1)}$i) * conjg(data(rowptr(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
+                end do
+            #:endif
+            end if
+        end associate
+    end subroutine
+
+    #:endfor
+    #:endfor
+
+end submodule stdlib_sparse_spmv_csr