spgemv_csr Subroutine

subroutine spgemv_csr(m, a, ia, ja, x, y, transa)

   !> Number of rows of the matrix A
   integer, intent(in) :: m

   !> Non-zero elements of the matrix A
   real(wp), intent(in) :: a(:)

   !> Row offsets of A into a and ja, size m + 1
   integer(i8), intent(in) :: ia(:)

   !> Column indices of the non-zero elements of A
   integer, intent(in) :: ja(:)

   !> Input vector
   real(wp), intent(in) :: x(:)

   !> Output vector, its first m elements are overwritten
   real(wp), intent(inout) :: y(:)

   !> Operation selector, 'N' for A*x (default), 'T' or 'C' for A^T*x
   character(len=1), intent(in), optional :: transa

   integer :: i, j
   integer(i8) :: k
   logical :: trans
   real(wp) :: y_tmp_i
   real(wp), allocatable :: y_priv(:)

   trans = .false.
   if (present(transa)) trans = scan(transa, "TtCc") > 0

   if (.not. trans) then
      !$omp parallel do default(none) schedule(guided) &
      !$omp& private(i, k, y_tmp_i) shared(m, a, ia, ja, x, y)
      do i = 1, m
         y_tmp_i = 0.0_wp
         do k = ia(i), ia(i+1) - 1
            y_tmp_i = y_tmp_i + a(k) * x(ja(k))
         end do
         y(i) = y_tmp_i
      end do
      !$omp end parallel do
   else
      y(:m) = 0.0_wp

      !$omp parallel default(none) &
      !$omp& private(y_priv, i, k, j) shared(m, a, ia, ja, x, y)
      allocate(y_priv(m), source=0.0_wp)

      !$omp do schedule(guided)
      do i = 1, m
         do k = ia(i), ia(i+1) - 1
            j = ja(k)
            y_priv(j) = y_priv(j) + a(k) * x(i)
         end do
      end do
      !$omp end do

      !$omp critical
      y(:m) = y(:m) + y_priv
      !$omp end critical

      deallocate(y_priv)
      !$omp end parallel
   end if

end subroutine spgemv_csr