Dia/OpenMP/Morpheus_Multiply_Impl.hpp

 
#ifndef MORPHEUS_DIA_OPENMP_MULTIPLY_IMPL_HPP
#define MORPHEUS_DIA_OPENMP_MULTIPLY_IMPL_HPP
 
#include <Morpheus_Macros.hpp>
#if defined(MORPHEUS_ENABLE_OPENMP)
 
#include <Morpheus_SpaceTraits.hpp>
#include <Morpheus_FormatTraits.hpp>
#include <Morpheus_FormatTags.hpp>
#include <Morpheus_Spaces.hpp>
 
namespace Morpheus {
namespace Impl {
 
template <typename ExecSpace, typename Matrix, typename Vector>
inline void multiply(
    const Matrix& A, const Vector& x, Vector& y, const bool init,
    typename std::enable_if_t<
        Morpheus::is_dia_matrix_format_container_v<Matrix> &&
        Morpheus::is_dense_vector_format_container_v<Vector> &&
        Morpheus::has_custom_backend_v<ExecSpace> &&
        Morpheus::has_openmp_execution_space_v<ExecSpace> &&
        Morpheus::has_access_v<ExecSpace, Matrix, Vector>>* = nullptr) {
  using size_type       = typename Matrix::size_type;
  using index_type      = typename Matrix::index_type;
  using value_type      = typename Matrix::value_type;
  const size_type ndiag = A.cdiagonal_offsets().size();
 
#pragma omp parallel for
  for (size_type row = 0; row < A.nrows(); row++) {
    value_type sum = init ? value_type(0) : y[row];
 
    for (size_type n = 0; n < ndiag; n++) {
      const index_type col = row + A.cdiagonal_offsets(n);
 
      if (col >= 0 && col < (index_type)A.ncols()) {
        sum += A.cvalues(row, n) * x[col];
      }
    }
    y[row] = sum;
  }
}
 
}  // namespace Impl
}  // namespace Morpheus
 
#endif  // MORPHEUS_ENABLE_OPENMP
#endif  // MORPHEUS_DIA_OPENMP_MULTIPLY_IMPL_HPP