spmspv_csr Subroutine

subroutine spmspv_csr(list, mlist, xptr, xval, yptr, yval, alpha)

   !> CSR neighbour-list structure
   type(csr_list), intent(in) :: list

   !> Matrix elements in CSR order
   real(wp), intent(in) :: mlist(:)

   !> Row indices of the non-zero elements of the input vector
   integer, intent(in) :: xptr(:)

   !> Values of the non-zero elements of the input vector
   real(wp), intent(in) :: xval(:)

   !> Row indices of the non-zero elements of the product vector, allocated here
   integer, allocatable, intent(out) :: yptr(:)

   !> Values of the non-zero elements of the product vector, allocated here
   real(wp), allocatable, intent(out) :: yval(:)

   !> Matrix scaling factor, defaults to one
   real(wp), intent(in), optional :: alpha

   integer :: i, j, n, m, nnz
   integer(i8) :: k
   real(wp) :: a, y_tmp_i
   real(wp), allocatable :: xdense(:), ydense(:)
   logical, allocatable :: yflag(:)

   a = 1.0_wp
   if (present(alpha)) a = alpha

   n = size(list%inl) - 1

   if (size(mlist) /= size(list%nlat) .or. n < 0) then
      allocate(yptr(0), yval(0))
      return
   end if

   allocate(xdense(n), source=0.0_wp)
   m = min(size(xptr), size(xval))
   do i = 1, m
      j = xptr(i)
      if (j >= 1 .and. j <= n) xdense(j) = xdense(j) + xval(i)
   end do

   allocate(ydense(n))
   allocate(yflag(n), source=.false.)

   !$omp parallel do default(none) schedule(guided) &
   !$omp& private(i, k, j, y_tmp_i) &
   !$omp& shared(list, mlist, xdense, ydense, yflag, a, n)
   do i = 1, n
      y_tmp_i = 0.0_wp
      do k = list%inl(i), list%inl(i+1) - 1
         j = list%nlat(k)
         y_tmp_i = y_tmp_i + mlist(k) * xdense(j)
      end do
      ydense(i) = a * y_tmp_i
      yflag(i) = ydense(i) /= 0.0_wp
   end do
   !$omp end parallel do

   nnz = count(yflag)
   allocate(yptr(nnz))
   allocate(yval(nnz))

   nnz = 0
   do i = 1, n
      if (yflag(i)) then
         nnz = nnz + 1
         yptr(nnz) = i
         yval(nnz) = ydense(i)
      end if
   end do

end subroutine spmspv_csr