Commit 5d5f2a1a authored by Müller, Felix's avatar Müller, Felix Committed by Praetorius, Simon
Browse files

Only create SparsityPattern when required

Remove BiLinearForm::pattern_
Change Pattern to be created in BiLinearForm::init()
Change BiLinearForm Observer to set a flag instead of directly rebuilding the pattern
Change Pattern ctor to replace call to init
Make Pattern::init protected
Add MatrixFacade::init() and implementations MatrixBackend::init() with no argument
Change MTL matrix backend slot size calculation
parent 1acbfeb0
......@@ -55,8 +55,8 @@ namespace AMDiS
, symmetry_(SymmetryStructure::unknown)
, rowBasis_(rowBasis)
, colBasis_(colBasis)
, updatePattern_(true)
{
pattern_.init(*rowBasis_, *colBasis_, symmetry_);
operators_.init(*rowBasis_, *colBasis_);
auto const rowSize = rowBasis_->localView().maxSize();
......@@ -129,13 +129,14 @@ namespace AMDiS
}
/// @}
/// Set the symmetry property of the bilinear form
/// Provide some additional information about the structure of the values of the pattern
void setSymmetryStructure(SymmetryStructure symm)
{
updatePattern_ = (symmetry_ != symm); // Force rebuild if symmetry structure changes
symmetry_ = symm;
}
/// Set the symmetry property using a string
/// Set the symmetry structure using a string
void setSymmetryStructure(std::string str) { setSymmetryStructure(symmetryStructure(str)); }
/// Assemble the matrix operators on the bound element.
......@@ -148,32 +149,25 @@ namespace AMDiS
/// Assemble all matrix operators, TODO: incorporate boundary conditions
void assemble();
void init()
/// Prepare the underlying matrix for insertion and rebuild the sparsity pattern if required
void init(bool forcePatternRebuild = false)
{
Super::init(pattern_);
}
void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
void updateImpl(event::adapt e, index_t<1> i) override { updateImpl2(e,i); }
private:
// Track for each basis wether updated and if both are, update the sparsity pattern
template <std::size_t I>
void updateImpl2(event::adapt, index_t<I>)
{
assert(!updateCounter_.test(I));
updateCounter_.set(I);
if (updateCounter_.all()) {
pattern_.init(*rowBasis_, *colBasis_, symmetry_);
updateCounter_.reset();
if (updatePattern_ || forcePatternRebuild)
{
Super::init(SparsityPattern{*rowBasis_, *colBasis_, symmetry_});
updatePattern_ = false;
}
else
{
Super::init();
}
}
protected:
/// The structure of the non-zeros in the matrix
SparsityPattern pattern_;
/// Set the updatePattern_ flag since the pattern must be rebuilt after a basis change
void updateImpl(event::adapt e, index_t<0> i) override { updatePattern_ = true; }
void updateImpl(event::adapt e, index_t<1> i) override { updatePattern_ = true; }
protected:
/// The symmetry property if the bilinear form
SymmetryStructure symmetry_;
......@@ -187,7 +181,8 @@ namespace AMDiS
std::shared_ptr<ColBasis const> colBasis_;
private:
std::bitset<2> updateCounter_ = 0;
/// The pattern is rebuilt on calling init if this flag is set
bool updatePattern_;
};
......
......@@ -29,6 +29,7 @@ addOperator(ContextTag contextTag, Expr const& expr,
auto localAssembler = makeUniquePtr(makeAssembler<Tr>(std::move(op), i, j));
operators_[i][j].push(contextTag, std::move(localAssembler));
updatePattern_ = true;
}
......
......@@ -39,18 +39,19 @@ namespace AMDiS
Impl const& impl() const { return impl_; }
Impl& impl() { return impl_; }
/// \brief Initialize the matrix for insertion, i.e. allocate the non-zero pattern
/**
* With the optional parameter \p symmetry some additional information about the
* structure of the values or the sparsity pattern can be provided.
* See \ref SymmetryStructure.
**/
/// Initialize the matrix for insertion and allocate the non-zero pattern
template <class SparsityPattern>
void init(SparsityPattern const& pattern)
{
impl_.init(pattern);
}
/// Initialize the matrix for insertion while keeping the pattern unchanged
void init()
{
impl_.init();
}
/// Finish the matrix insertion, e.g. cleanup or final insertion
void finish()
{
......
......@@ -14,6 +14,28 @@ namespace AMDiS
class SparsityPattern
{
public:
template <class Basis>
SparsityPattern(Basis const& basis, SymmetryStructure symmetry = SymmetryStructure::unknown)
: rows_(basis.dimension())
, cols_(rows_)
, pattern_(rows_, cols_)
{
init(basis);
}
template <class RowBasis, class ColBasis>
SparsityPattern(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
: rows_(rowBasis.dimension())
, cols_(colBasis.dimension())
, pattern_(rows_, cols_)
{
if (uintptr_t(&rowBasis) == uintptr_t(&colBasis))
init(rowBasis);
else
init(rowBasis, colBasis);
}
/// Number of rows in the matrix
std::size_t rows() const
{
......@@ -55,15 +77,12 @@ namespace AMDiS
pattern_.exportIdx(matrix);
}
// Update pattern when basis is updated. This method is called if rowBasis == colBasis.
protected:
// Create pattern from basis. This method is called if rowBasis == colBasis.
template <class Basis>
void init(Basis const& basis, SymmetryStructure symmetry = SymmetryStructure::unknown)
void init(Basis const& basis)
{
rows_ = basis.dimension();
cols_ = rows_;
pattern_.resize(0, 0); // clear the old pattern
pattern_.resize(rows_, cols_);
auto localView = basis.localView();
for (const auto& element : elements(basis.gridView())) {
localView.bind(element);
......@@ -80,19 +99,10 @@ namespace AMDiS
}
}
// Update pattern when basis is updated. This method is called if rowBasis != colBasis.
// Create pattern from basis. This method is called if rowBasis != colBasis.
template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
void init(RowBasis const& rowBasis, ColBasis const& colBasis)
{
if (uintptr_t(&rowBasis) == uintptr_t(&colBasis))
return init(rowBasis);
rows_ = rowBasis.dimension();
cols_ = colBasis.dimension();
pattern_.resize(0, 0); // clear the old pattern
pattern_.resize(rows_, cols_);
auto rowLocalView = rowBasis.localView();
auto colLocalView = colBasis.localView();
for (const auto& element : elements(rowBasis.gridView())) {
......
......@@ -72,7 +72,7 @@ namespace AMDiS
}
/// Create inserter. Assumes that no inserter is currently active on this matrix.
/// Resize the matrix according to the pattern provided and set all entries to zero.
template <class Pattern>
void init(Pattern const& pattern)
{
......@@ -80,8 +80,13 @@ namespace AMDiS
matrix_.setZero();
}
/// Delete inserter -> finish insertion. Must be called in order to fill the
/// final construction of the matrix.
/// Set all matrix entries to zero while keeping the size unchanged.
void init()
{
matrix_.setZero();
}
/// Set the matrix entries from the triplet list and compress the matrix afterwards.
void finish()
{
matrix_.setFromTriplets(tripletList_.begin(), tripletList_.end());
......
......@@ -11,6 +11,13 @@ namespace AMDiS
class MatrixSize
{
public:
template <class RowBasis, class ColBasis>
MatrixSize(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
: rows_(rowBasis.dimension())
, cols_(colBasis.dimension())
{}
/// Number of rows in the matrix
std::size_t rows() const
{
......@@ -23,17 +30,9 @@ namespace AMDiS
return cols_;
}
template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
{
rows_ = rowBasis.dimension();
cols_ = colBasis.dimension();
}
private:
std::size_t rows_ = 0;
std::size_t cols_ = 0;
std::size_t rows_;
std::size_t cols_;
};
} // end namespace AMDiS
......@@ -61,7 +61,7 @@ namespace AMDiS
Comm const& comm() const { return *comm_; }
/// create occupation pattern and apply it to the matrix
/// Create occupation pattern and apply it to the matrix
template <class Pattern>
void init(Pattern const& pattern)
{
......@@ -69,6 +69,13 @@ namespace AMDiS
initialized_ = true;
}
/// Set all entries to zero while keeping the occupation pattern intact
void init()
{
matrix_ = value_type(0);
initialized_ = true;
}
void finish()
{
initialized_ = false;
......
......@@ -56,19 +56,34 @@ namespace AMDiS
return matrix_;
}
/// Create inserter. Assumes that no inserter is currently active on this matrix.
/// Create inserter. Assumes that no inserter is currently active on this matrix. Resizes the
/// matrix according to the provided pattern.
void init(Pattern const& pattern)
{
test_exit(!inserter_, "Matrix already in insertion mode!");
std::size_t slotSize = pattern.avgRowSize(matrix_);
std::size_t slotSize = nnz() > 0 ? 6*nnz() / (5*num_rows(matrix_))
: pattern.rowSizeEstimate();
matrix_.change_dim(pattern.rows(), pattern.cols());
set_to_zero(matrix_);
inserter_ = new Inserter(matrix_, slotSize);
symmetry_ = pattern.symmetry();
inserter_ = new Inserter(matrix_, slotSize);
}
/// Create inserter. Assumes that no inserter is currently active on this matrix. Does not
/// change matrix dimensions.
void init()
{
test_exit(!inserter_, "Matrix already in insertion mode!");
std::size_t slotSize = 6*nnz() / (5*num_rows(matrix_));
set_to_zero(matrix_);
inserter_ = new Inserter(matrix_, slotSize);
}
/// Delete inserter -> finish insertion. Must be called in order to fill the
/// final construction of the matrix.
void finish()
......
#pragma once
#include <algorithm>
#include <cassert>
#include <functional>
#include <memory>
#include <dune/common/math.hh>
#include <amdis/Observer.hpp>
......@@ -15,6 +15,16 @@ namespace AMDiS
class SlotSize
{
public:
template <class RowBasis, class ColBasis>
SlotSize(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
: rows_(rowBasis.dimension())
, cols_(colBasis.dimension())
, symmetry_(symmetry)
{
init(rowBasis, colBasis);
}
/// Number of rows in the matrix
std::size_t rows() const
{
......@@ -27,14 +37,10 @@ namespace AMDiS
return cols_;
}
/// Average number of non-zeros per row
/// In the first time this function is called, use the estRowSize,
/// otherwise take `1.2 * nnz/rows`
template <class Matrix>
std::size_t avgRowSize(Matrix const& matrix) const
/// Estimate of the non-zeros per row
std::size_t rowSizeEstimate() const
{
assert(rows_ > 0 && cols_ > 0);
return matrix.nnz() > 0 ? 6*matrix.nnz() / (5*rows_) : estRowSize_;
return estRowSize_;
}
/// Symmetry of the matrix entries
......@@ -43,29 +49,26 @@ namespace AMDiS
return symmetry_;
}
protected:
// estimate the number of columns by multiplying the maximal node size with the
// number of element surrounding a vertex. This number is approximated by the
// number of simplices surrounding a vertex in a kuhn triangulation
template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
void init(RowBasis const& rowBasis, ColBasis const& colBasis)
{
using GridView = typename RowBasis::GridView;
// TODO(FM): Simplify using Dune::factorial with 2.7
static std::size_t surrounding
= Math::pow<GridView::dimension>(2) * Dune::Factorial<int(GridView::dimension)>::factorial;
rows_ = rowBasis.dimension();
cols_ = colBasis.dimension();
estRowSize_ = std::min(cols_, colBasis.localView().maxSize() * surrounding);
symmetry_ = symmetry;
}
private:
std::size_t rows_ = 0;
std::size_t cols_ = 0;
std::size_t estRowSize_ = 50;
SymmetryStructure symmetry_ = SymmetryStructure::unknown;
std::size_t rows_;
std::size_t cols_;
std::size_t estRowSize_;
SymmetryStructure symmetry_;
};
} // end namespace AMDiS
......@@ -147,6 +147,13 @@ namespace AMDiS
initialized_ = true;
}
/// Reuse the matrix pattern and set all entries to zero
void init()
{
MatZeroEntries(matrix_);
initialized_ = true;
}
/// Finish assembly. Must be called before matrix can be used in a KSP
void finish()
{
......
......@@ -16,6 +16,14 @@ namespace AMDiS
class MatrixNnzStructure
{
public:
template <class RowBasis, class ColBasis>
MatrixNnzStructure(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
: symmetry_(symmetry)
{
init(rowBasis, colBasis);
}
// Return Number of nonzeros in the diagonal part (owner part)
std::vector<PetscInt> const& d_nnz() const
{
......@@ -34,19 +42,15 @@ namespace AMDiS
return symmetry_;
}
// Update pattern when basis is updated
protected:
// Create pattern from basis
template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown);
void init(RowBasis const& rowBasis, ColBasis const& colBasis);
private:
std::vector<PetscInt> dnnz_; //< number of nonzeros in the diagonal part (owner part)
std::vector<PetscInt> onnz_; //< number of nonzeros in the off-diagonal part (overlap part)
SymmetryStructure symmetry_ = SymmetryStructure::unknown;
// use a grid change index to decide whether to update the matrix pattern or not.
unsigned long changeIndex_ = 0;
bool initialized_ = false;
SymmetryStructure symmetry_;
#if HAVE_MPI
const Mpi::Tag tag_{916821}; //< Communication tag used internally
......
......@@ -17,12 +17,10 @@
namespace AMDiS {
template <class RowBasis, class ColBasis>
void MatrixNnzStructure::init(RowBasis const& rowBasis, ColBasis const& /*colBasis*/,
SymmetryStructure symmetry)
void MatrixNnzStructure::init(RowBasis const& rowBasis, ColBasis const& /*colBasis*/)
{
Dune::Timer t;
symmetry_ = symmetry;
auto const& basis = rowBasis; // TODO: generalize to row != col basis
auto const& dofMap = basis.communicator().dofMap();
std::size_t localSize = dofMap.localSize();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment