#pragma once #include #include #include #include #include #include #include #include #include namespace AMDiS { /// \brief The basic container that stores a base matrix template class MatrixBackend { public: using Traits = TraitsType; /// The matrix type of the underlying base matrix using BaseMatrix = typename Traits::Mat; /// The type of the elements of the DOFMatrix using value_type = typename BaseMatrix::value_type; /// The index/size - type using size_type = typename BaseMatrix::size_type; /// The type of the inserter used when filling the matrix. NOTE: need not be public using Inserter = mtl::mat::inserter>; public: /// Constructor. Constructs new BaseMatrix. MatrixBackend(std::shared_ptr const&) {} /// Return a reference to the data-matrix \ref matrix BaseMatrix& matrix() { assert( !inserter_ ); return matrix_; } /// Return a reference to the data-matrix \ref matrix BaseMatrix const& matrix() const { assert( !inserter_ ); return matrix_; } /// Create inserter. Assumes that no inserter is currently active on this matrix. template void init(RowBasis const& rowBasis, ColBasis const& colBasis, SymmetryStructure symmetry) { test_exit(!inserter_, "Matrix already in insertion mode!"); calculateSlotSize(); matrix_.change_dim(rowBasis.dimension(), colBasis.dimension()); set_to_zero(matrix_); inserter_ = new Inserter(matrix_, slotSize_); symmetry_ = symmetry; } /// Delete inserter -> finish insertion. Must be called in order to fill the /// final construction of the matrix. void finish() { delete inserter_; inserter_ = nullptr; } /// \brief Returns an update-proxy of the inserter, to insert/update a value at /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can /// be created by \ref init and must be closed by \ref finish after insertion. void insert(size_type r, size_type c, value_type const& value) { test_exit_dbg(inserter_, "Inserter not initilized!"); test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_), "Indices out of range [0,{})x[0,{})", num_rows(matrix_), num_cols(matrix_)); if (value != value_type(0) || r == c) (*inserter_)[r][c] += value; } template void scatter(Ind const& idx, LocalMat const& mat) { scatter(idx, idx, mat); } template void scatter(RowInd const& rows, ColInd const& cols, LocalMat const& mat) { test_exit_dbg(inserter_, "Inserter not initilized!"); for (size_type i = 0; i < size_type(rows.size()); ++i) for (size_type j = 0; j < size_type(cols.size()); ++j) if (mat[i][j] != value_type(0) || i == j) (*inserter_)[rows[i]][cols[j]] += mat[i][j]; } /// Return the number of nonzeros in the matrix std::size_t nnz() const { return matrix_.nnz(); } SymmetryStructure symmetry() const { return symmetry_; } protected: // Estimates the slot-size used in the inserter to reserve enough space per row. void calculateSlotSize() { slotSize_ = 0; if (num_rows(matrix_) != 0) slotSize_ = int(double(matrix_.nnz()) / num_rows(matrix_) * 1.2); if (slotSize_ < MIN_NNZ_PER_ROW) slotSize_ = MIN_NNZ_PER_ROW; } private: /// The minimal number of nonzeros per row static constexpr int MIN_NNZ_PER_ROW = 50; /// The data-matrix (can hold a new BaseMatrix or a pointer to external data BaseMatrix matrix_; /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish Inserter* inserter_ = nullptr; /// The size of the slots used to insert values per row int slotSize_ = MIN_NNZ_PER_ROW; SymmetryStructure symmetry_ = SymmetryStructure::unknown; }; } // end namespace AMDiS