complete spmv for the ellpack format including symmetric representations

Author: jalvesz

doc/specs/stdlib_sparse.md

@@ -338,6 +338,18 @@ If the `diagonal` array has not been previously allocated, the `diag` subroutine
 
 ### Syntax
 
+`call ` [[stdlib_sparse_conversion(module):csr2sellc(interface)]] `(csr,ell[,num_nz_rows])`
+
+### Arguments
+
+`csr` : Shall be a `CSR` type of `real` or `complex` type. It is an `intent(in)` argument.
+
+`ell` : Shall be a `ELL` type of `real` or `complex` type. It is an `intent(out)` argument.
+
+`num_nz_rows`, `optional`: number of non zeros per row. If not give, it will correspond to the size of the longest row in the `CSR` matrix. It is an `intent(in)` argument.
+
+### Syntax
+
 `call ` [[stdlib_sparse_conversion(module):csc2coo(interface)]] `(csc,coo)`
 
 ### Arguments

src/stdlib_sparse_conversion.fypp

@@ -378,6 +378,7 @@ contains
             end do
         end if
         call ELL%malloc(CSR%nrows,CSR%ncols,num_nz_rows_)
+        ELL%storage = CSR%storage
         !-------------------------------------------
         do i = 1, CSR%nrows
             adr1 = CSR%rowptr(i)

src/stdlib_sparse_spmv.fypp

@@ -409,12 +409,30 @@ contains
                     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)
                 end do
+            else if( storage /= sparse_full .and. op_/=sparse_op_hermitian ) then
+                do concurrent (i = 1:nrows, k = 1:MNZ_P_ROW)
+                    j = index(i,k)
+                    if(j>0) then
+                        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)
+                    end if
+                end do
             #:if t1.startswith('complex')
             else if( storage == sparse_full .and. op_==sparse_op_hermitian ) then
                 do concurrent (i = 1:nrows, 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)
                 end do
+            else if( storage /= sparse_full .and. op_==sparse_op_hermitian ) then
+                do concurrent (i = 1:nrows, k = 1:MNZ_P_ROW)
+                    j = index(i,k)
+                    if(j>0) then
+                        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)
+                    end if
+                end do
             #:endif
             end if
         end associate

test/linalg/test_sparse_spmv.fypp

@@ -243,14 +243,16 @@ contains
             integer, parameter :: wp = ${k1}$
             type(COO_${s1}$_type) :: COO
             type(CSR_${s1}$_type) :: CSR
+            type(ELL_${s1}$_type) :: ELL
             ${t1}$, allocatable :: dense(:,:)
             ${t1}$, allocatable :: vec_x(:)
-            ${t1}$, allocatable :: vec_y1(:), vec_y2(:), vec_y3(:)
+            ${t1}$, allocatable :: vec_y1(:), vec_y2(:), vec_y3(:), vec_y4(:)
 
             allocate( vec_x(4)  , source = 1._wp )
             allocate( vec_y1(4) , source = 0._wp )
             allocate( vec_y2(4) , source = 0._wp )
             allocate( vec_y3(4) , source = 0._wp )
+            allocate( vec_y4(4) , source = 0._wp )
 
             allocate( dense(4,4) , source = &
                     reshape(real([1,0,0,0, &
@@ -261,6 +263,7 @@ contains
             call dense2coo( dense , COO )
             COO%storage = sparse_upper
             call coo2csr(COO, CSR)
+            call csr2ell(CSR, ELL)
 
             dense(2,1) = 2._wp; dense(3,2) = 2._wp; dense(4,3) = 2._wp
             vec_y1 = matmul( dense, vec_x )
@@ -274,6 +277,10 @@ contains
             call spmv( CSR , vec_x, vec_y3 )
             call check(error, all(vec_y1 == vec_y3) )
             if (allocated(error)) return
+
+            call spmv( ELL , vec_x, vec_y4 )
+            call check(error, all(vec_y1 == vec_y4) )
+            if (allocated(error)) return
         end block
         #:endfor
     end subroutine