spsymv_csr Subroutine

subroutine spsymv_csr(m, a, ia, ja, x, y, uplo)

   !> 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(:)

   !> Triangle selector, 'U' for the upper (default), 'L' for the lower one
   character(len=1), intent(in), optional :: uplo

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

   upper = .true.
   if (present(uplo)) upper = scan(uplo, "Uu") > 0

   y(:m) = 0.0_wp

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

   !$omp do schedule(guided)
   do i = 1, m
      y_tmp_i = 0.0_wp
      do k = ia(i), ia(i+1) - 1
         j = ja(k)
         if (j == i) then
            y_tmp_i = y_tmp_i + a(k) * x(i)
         else if ((j > i) .eqv. upper) then
            ! Contributions to row i and to the mirrored row j
            y_tmp_i = y_tmp_i + a(k) * x(j)
            y_priv(j) = y_priv(j) + a(k) * x(i)
         end if
      end do
      y_priv(i) = y_priv(i) + y_tmp_i
   end do
   !$omp end do

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

   deallocate(y_priv)
   !$omp end parallel

end subroutine spsymv_csr