Morpheus_DiaMatrix.hpp

 
#ifndef MORPHEUS_DIAMATRIX_HPP
#define MORPHEUS_DIAMATRIX_HPP
 
#include <Morpheus_Exceptions.hpp>
#include <Morpheus_FormatTags.hpp>
#include <Morpheus_DenseVector.hpp>
#include <Morpheus_DenseMatrix.hpp>
#include <Morpheus_DynamicMatrix.hpp>
#include <Morpheus_MatrixBase.hpp>
 
#include <impl/Dia/Morpheus_Utils_Impl.hpp>
 
namespace Morpheus {
 
template <class ValueType, class... Properties>
class DiaMatrix : public MatrixBase<DiaMatrix, ValueType, Properties...> {
 public:
  using traits = ContainerTraits<DiaMatrix, ValueType, Properties...>;
  using type = typename traits::type;
  using base = MatrixBase<DiaMatrix, ValueType, Properties...>;
  using tag = typename MatrixFormatTag<DiaFormatTag>::tag;
 
  using value_type = typename traits::value_type;
  using non_const_value_type = typename traits::non_const_value_type;
  using size_type            = typename traits::size_type;
  using index_type = typename traits::index_type;
  using non_const_index_type = typename traits::non_const_index_type;
 
  using array_layout    = typename traits::array_layout;
  using backend         = typename traits::backend;
  using memory_space    = typename traits::memory_space;
  using execution_space = typename traits::execution_space;
  using device_type     = typename traits::device_type;
  using memory_traits   = typename traits::memory_traits;
  using HostMirror      = typename traits::HostMirror;
 
  using pointer         = typename traits::pointer;
  using const_pointer   = typename traits::const_pointer;
  using reference       = typename traits::reference;
  using const_reference = typename traits::const_reference;
 
  using index_array_type =
      Morpheus::DenseVector<index_type, size_type, array_layout, backend,
                            memory_traits>;
  using const_index_array_type = const index_array_type;
  using index_array_pointer    = typename index_array_type::value_array_pointer;
  using index_array_reference =
      typename index_array_type::value_array_reference;
  using const_index_array_reference = const index_array_reference;
 
  using value_array_type =
      Morpheus::DenseMatrix<value_type, size_type, array_layout, backend,
                            memory_traits>;
  using const_value_array_type = const value_array_type;
  using value_array_pointer    = typename value_array_type::value_array_pointer;
  using value_array_reference =
      typename value_array_type::value_array_reference;
  using const_value_array_reference = const value_array_reference;
 
  ~DiaMatrix() = default;
  DiaMatrix(const DiaMatrix &) = default;
  DiaMatrix(DiaMatrix &&) = default;
  DiaMatrix &operator=(const DiaMatrix &) = default;
  DiaMatrix &operator=(DiaMatrix &&) = default;
 
  inline DiaMatrix()
      : base(), _ndiags(0), _alignment(0), _diagonal_offsets(), _values() {}
 
  inline DiaMatrix(const size_type num_rows, const size_type num_cols,
                   const size_type num_entries, const size_type num_diagonals,
                   const size_type alignment = 32)
      : base(num_rows, num_cols, num_entries),
        _ndiags(num_diagonals),
        _alignment(alignment),
        _diagonal_offsets(num_diagonals) {
    _values.resize(Impl::get_pad_size<size_type>(num_rows, alignment),
                   num_diagonals);
  }
 
  template <typename ValueArray, typename IndexArray>
  explicit inline DiaMatrix(
      const size_type num_rows, const size_type num_cols,
      const size_type num_entries, const IndexArray &diag_offsets,
      const ValueArray &vals,
      typename std::enable_if<
          is_dense_matrix_format_container<ValueArray>::value &&
          is_dense_vector_format_container<IndexArray>::value &&
          is_compatible<typename DiaMatrix::value_array_type,
                        ValueArray>::value &&
          is_compatible<typename DiaMatrix::index_array_type,
                        IndexArray>::value &&
          !ValueArray::memory_traits::is_unmanaged &&
          !IndexArray::memory_traits::is_unmanaged>::type * = nullptr)
      : base(num_rows, num_cols, num_entries),
        _diagonal_offsets(diag_offsets),
        _values(vals) {
    _ndiags    = _diagonal_offsets.size();
    _alignment = _values.nrows();
  }
 
  // Construct from raw pointers
  template <typename ValuePtr, typename IndexPtr>
  explicit inline DiaMatrix(
      const size_type num_rows, const size_type num_cols,
      const size_type num_entries, IndexPtr diag_offsets_ptr, ValuePtr vals_ptr,
      const size_type num_diagonals, const size_type alignment = 32,
      typename std::enable_if<
          (std::is_pointer<ValuePtr>::value &&
           is_same_value_type<value_type, ValuePtr>::value &&
           memory_traits::is_unmanaged) &&
          (std::is_pointer<IndexPtr>::value &&
           is_same_index_type<index_type, IndexPtr>::value &&
           memory_traits::is_unmanaged)>::type * = nullptr)
      : base(num_rows, num_cols, num_entries),
        _diagonal_offsets(num_diagonals, diag_offsets_ptr),
        _values(Impl::get_pad_size<size_type>(num_rows, alignment),
                num_diagonals, vals_ptr) {}
 
  template <class VR, class... PR>
  DiaMatrix(const DiaMatrix<VR, PR...> &src,
            typename std::enable_if<is_format_compatible<
                DiaMatrix, DiaMatrix<VR, PR...>>::value>::type * = nullptr)
      : base(src.nrows(), src.ncols(), src.nnnz()),
        _ndiags(src.ndiags()),
        _alignment(src.alignment()),
        _diagonal_offsets(src.cdiagonal_offsets()),
        _values(src.cvalues()) {}
 
  template <class VR, class... PR>
  typename std::enable_if<
      is_format_compatible<DiaMatrix, DiaMatrix<VR, PR...>>::value,
      DiaMatrix &>::type
  operator=(const DiaMatrix<VR, PR...> &src) {
    this->set_nrows(src.nrows());
    this->set_ncols(src.ncols());
    this->set_nnnz(src.nnnz());
 
    _ndiags           = src.ndiags();
    _alignment        = src.alignment();
    _diagonal_offsets = src.cdiagonal_offsets();
    _values           = src.cvalues();
 
    return *this;
  }
 
  template <class VR, class... PR>
  DiaMatrix(const DynamicMatrix<VR, PR...> &src,
            typename std::enable_if<is_dynamically_compatible<
                DiaMatrix, DynamicMatrix<VR, PR...>>::value>::type * = nullptr)
      : base(src.nrows(), src.ncols(), src.nnnz()) {
    auto f = std::bind(Impl::any_type_assign(), std::placeholders::_1,
                       std::ref(*this));
 
    std::visit(f, src.const_formats());
  }
 
  template <class VR, class... PR>
  typename std::enable_if<
      is_dynamically_compatible<DiaMatrix, DynamicMatrix<VR, PR...>>::value,
      DiaMatrix &>::type
  operator=(const DynamicMatrix<VR, PR...> &src) {
    auto f = std::bind(Impl::any_type_assign(), std::placeholders::_1,
                       std::ref(*this));
 
    std::visit(f, src.const_formats());
 
    return *this;
  }
 
  template <typename MatrixType>
  DiaMatrix(const MatrixType &src) = delete;
 
  template <typename MatrixType>
  reference operator=(const MatrixType &src) = delete;
 
  inline void resize(const size_type num_rows, const size_type num_cols,
                     const size_type num_entries, const size_type num_diagonals,
                     const size_type alignment = 32) {
    base::resize(num_rows, num_cols, num_entries);
    _ndiags    = num_diagonals;
    _alignment = alignment;
 
    if (Impl::exceeds_tolerance(num_rows, num_entries, _ndiags)) {
      throw Morpheus::FormatConversionException(
          "DiaMatrix fill-in would exceed maximum tolerance");
    }
 
    _diagonal_offsets.resize(num_diagonals);
    _values.resize(Impl::get_pad_size<size_type>(num_rows, alignment),
                   num_diagonals);
  }
 
  template <class VR, class... PR>
  inline void resize(const DiaMatrix<VR, PR...> &src) {
    resize(src.nrows(), src.ncols(), src.nnnz(), src.ndiags(), src.alignment());
  }
 
  template <class VR, class... PR>
  inline DiaMatrix &allocate(const DiaMatrix<VR, PR...> &src) {
    resize(src.nrows(), src.ncols(), src.nnnz(), src.ndiags(), src.alignment());
    return *this;
  }
 
  formats_e format_enum() const { return _id; }
 
  int format_index() const { return static_cast<int>(_id); }
 
  MORPHEUS_FORCEINLINE_FUNCTION index_array_reference
  diagonal_offsets(size_type n) {
    return _diagonal_offsets(n);
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION value_array_reference values(size_type i,
                                                             size_type j) {
    return _values(i, j);
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION const_index_array_reference
  cdiagonal_offsets(size_type n) const {
    return _diagonal_offsets(n);
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION const_value_array_reference
  cvalues(size_type i, size_type j) const {
    return _values(i, j);
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION index_array_type &diagonal_offsets() {
    return _diagonal_offsets;
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION value_array_type &values() { return _values; }
 
  MORPHEUS_FORCEINLINE_FUNCTION const_index_array_type &cdiagonal_offsets()
      const {
    return _diagonal_offsets;
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION const_value_array_type &cvalues() const {
    return _values;
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION size_type ndiags() const { return _ndiags; }
 
  MORPHEUS_FORCEINLINE_FUNCTION size_type alignment() const {
    return _alignment;
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION void set_ndiags(const size_type num_diagonals) {
    _ndiags = num_diagonals;
  }
 
  MORPHEUS_FORCEINLINE_FUNCTION void set_alignment(const size_type alignment) {
    _alignment = alignment;
  }
 
 private:
  size_type _ndiags, _alignment;
  index_array_type _diagonal_offsets;
  value_array_type _values;
  static constexpr formats_e _id = Morpheus::DIA_FORMAT;
};
 
}  // namespace Morpheus
 
#endif  // MORPHEUS_DIAMATRIX_HPP