simplify sellc spmv kernel

Author: jalvesz

src/stdlib_sparse_spmv.fypp

@@ -418,7 +418,7 @@ contains
         character(1), intent(in), optional :: op
         ${t1}$ :: alpha_
         character(1) :: op_
-        integer(ilp) :: i, nz, rowidx, num_chunks, rm
+        integer(ilp) :: i, j, k, nz, rowidx, num_chunks, rm
 
         op_ = sparse_op_none; if(present(op)) op_ = op
         alpha_ = one_${s1}$
@@ -447,7 +447,12 @@ contains
                 do i = 1, num_chunks
                     nz = ia(i+1) - ia(i)
                     rowidx = (i - 1)*${chunk}$ + 1 
-                    call chunk_kernel_${chunk}$(nz,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:))
+                    associate(col => ja(1:${chunk}$,ia(i):ia(i)+nz-1), mat => data(1:${chunk}$,ia(i):ia(i)+nz-1), &
+                    &           x => vec_x, y => vec_y(rowidx:rowidx+${chunk}$-1) )
+                    do j = 1, nz
+                        where(col(:,j) > 0) y = y + alpha_ * mat(:,j) * x(col(:,j))
+                    end do
+                    end associate
                 end do
             #:endfor
             end select
@@ -457,7 +462,12 @@ contains
                 i = num_chunks + 1 
                 nz = ia(i+1) - ia(i)
                 rowidx = (i - 1)*cs + 1
-                call chunk_kernel_remainder(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x,vec_y(rowidx:))
+                associate(col => ja(1:${chunk}$,ia(i):ia(i)+nz-1), mat => data(1:${chunk}$,ia(i):ia(i)+nz-1), &
+                &           x => vec_x, y => vec_y(rowidx:rowidx+rm-1) )
+                do j = 1, nz
+                    where(col(1:rm,j) > 0) y = y + alpha_ * mat(1:rm,j) * x(col(1:rm,j))
+                end do
+                end associate
             end if
 
         else if( storage == sparse_full .and. op_==sparse_op_transpose ) then
@@ -468,7 +478,14 @@ contains
                 do i = 1, num_chunks
                     nz = ia(i+1) - ia(i)
                     rowidx = (i - 1)*${chunk}$ + 1 
-                    call chunk_kernel_trans_${chunk}$(nz,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y)
+                    associate(col => ja(1:${chunk}$,ia(i):ia(i)+nz-1), mat => data(1:${chunk}$,ia(i):ia(i)+nz-1), &
+                    &           x => vec_x(rowidx:rowidx+${chunk}$-1), y => vec_y )
+                    do j = 1, nz
+                        do k = 1, ${chunk}$
+                            if(col(k,j) > 0) y(col(k,j)) = y(col(k,j)) + alpha_ * mat(k,j) * x(k)
+                        end do
+                    end do
+                    end associate
                 end do
             #:endfor
             end select
@@ -478,63 +495,21 @@ contains
                 i = num_chunks + 1 
                 nz = ia(i+1) - ia(i)
                 rowidx = (i - 1)*cs + 1
-                call chunk_kernel_remainder_trans(nz,cs,rm,data(:,ia(i)),ja(:,ia(i)),vec_x(rowidx:),vec_y)
+                associate(col => ja(1:${chunk}$,ia(i):ia(i)+nz-1), mat => data(1:${chunk}$,ia(i):ia(i)+nz-1), &
+                &           x => vec_x(rowidx:rowidx+rm-1), y => vec_y )
+                do j = 1, nz
+                    do k = 1, rm
+                        if(col(k,j) > 0) y(col(k,j)) = y(col(k,j)) + alpha_ * mat(k,j) * x(k)
+                    end do
+                end do
+                end associate
             end if
         else
             print *, "error: sellc format for spmv operation not yet supported."
             return
         end if
         end associate
 
-    contains
-    #:for chunk in CHUNKS
-    pure subroutine chunk_kernel_${chunk}$(n,a,col,x,y)
-        integer, value      :: n
-        ${t1}$, intent(in)  :: a(${chunk}$,n), x(*)
-        integer(ilp), intent(in) :: col(${chunk}$,n)
-        ${t1}$, intent(inout) :: y(${chunk}$)
-        integer :: j
-        do j = 1, n
-            where(col(:,j) > 0) y = y + alpha_ * a(:,j) * x(col(:,j))
-        end do
-    end subroutine
-    pure subroutine chunk_kernel_trans_${chunk}$(n,a,col,x,y)
-        integer, value      :: n
-        ${t1}$, intent(in)  :: a(${chunk}$,n), x(${chunk}$)
-        integer(ilp), intent(in) :: col(${chunk}$,n)
-        ${t1}$, intent(inout) :: y(*)
-        integer :: j, k
-        do j = 1, n
-            do k = 1, ${chunk}$
-                if(col(k,j) > 0) y(col(k,j)) = y(col(k,j)) + alpha_ * a(k,j) * x(k)
-            end do
-        end do
-    end subroutine
-    #:endfor
-
-    pure subroutine chunk_kernel_remainder(n,cs,rm,a,col,x,y)
-        integer, value      :: n, cs, rm
-        ${t1}$, intent(in)  :: a(cs,n), x(*)
-        integer(ilp), intent(in) :: col(cs,n)
-        ${t1}$, intent(inout) :: y(rm)
-        integer :: j
-        do j = 1, n
-            where(col(1:rm,j) > 0) y = y + alpha_ * a(1:rm,j) * x(col(1:rm,j))
-        end do
-    end subroutine
-    pure subroutine chunk_kernel_remainder_trans(n,cs,rm,a,col,x,y)
-        integer, value      :: n, cs, rm
-        ${t1}$, intent(in)  :: a(cs,n), x(rm)
-        integer(ilp), intent(in) :: col(cs,n)
-        ${t1}$, intent(inout) :: y(*)
-        integer :: j, k
-        do j = 1, n
-            do k = 1, rm
-                if(col(k,j) > 0) y(col(k,j)) = y(col(k,j)) + alpha_ * a(k,j) * x(k)
-            end do
-        end do
-    end subroutine  
-
     end subroutine
 
     #:endfor