Merge branch 'feature/keep_pattern' into 'master'

Keep unchanged sparsity pattern

See merge request !166

amdis/BiLinearForm.hpp

@@ -52,10 +52,11 @@ namespace AMDiS
     BiLinearForm(std::shared_ptr<RB> const& rowBasis, std::shared_ptr<CB> const& colBasis)
       : Super(*rowBasis, *colBasis)
       , ObserverSequence<event::adapt,2>(*rowBasis, *colBasis)
+      , symmetry_(SymmetryStructure::unknown)
       , rowBasis_(rowBasis)
       , colBasis_(colBasis)
+      , updatePattern_(true)
     {
-      pattern_.init(*rowBasis_, *colBasis_);
       operators_.init(*rowBasis_, *colBasis_);
 
       auto const rowSize = rowBasis_->localView().maxSize();
@@ -128,6 +129,23 @@ namespace AMDiS
     }
     /// @}
 
+    /// Provide some additional information about the structure of the matrix pattern
+    /**
+     * If the symmetry structure is changed, it forces the matrix to
+     * rebuild its patterns.
+     **/
+    void setSymmetryStructure(SymmetryStructure symm)
+    {
+      updatePattern_ = (symmetry_ != symm);
+      symmetry_ = symm;
+    }
+
+    /// Set the symmetry structure using a string
+    void setSymmetryStructure(std::string str)
+    {
+      setSymmetryStructure(symmetryStructure(str));
+    }
+
     /// Assemble the matrix operators on the bound element.
     template <class RowLocalView, class ColLocalView,
       REQUIRES(Concepts::LocalView<RowLocalView>),
@@ -136,35 +154,39 @@ namespace AMDiS
                   ColLocalView const& colLocalView);
 
     /// Assemble all matrix operators, TODO: incorporate boundary conditions
-    void assemble(SymmetryStructure symmetry = SymmetryStructure::unknown);
-
-    void init(SymmetryStructure symmetry = SymmetryStructure::unknown)
-    {
-      Super::init(pattern_, symmetry);
-    }
-
-    using Super::init;
-
-    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); }
+    void assemble();
 
-  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>)
+    /// Prepare the underlying matrix for insertion and rebuild the sparsity pattern
+    /// if required.
+    /**
+     * The pattern is rebuild if either the parameter `forceRebuild` is set to true
+     * or if for some other reasons the internal flag `updatePattern_` is set.
+     * This might be the case if the symmetry structure is changed or the basis
+     * is updated.
+     *
+     * If the pattern is not updated, just the values are set to zero.
+     **/
+    void init(bool forcePatternRebuild = false)
     {
-      assert(!updateCounter_.test(I));
-      updateCounter_.set(I);
-
-      if (updateCounter_.all()) {
-        pattern_.init(*rowBasis_, *colBasis_);
-        updateCounter_.reset();
+      if (updatePattern_ || forcePatternRebuild)
+      {
+        Super::init(SparsityPattern{*rowBasis_, *colBasis_, symmetry_});
+        updatePattern_ = false;
+      }
+      else
+      {
+        Super::init();
       }
     }
 
+    /// Set the flag that forces an update of the pattern since the underlying
+    /// basis that defines the indexset has been changed
+    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 structure of the non-zeros in the matrix
-    SparsityPattern pattern_;
+    /// The symmetry property if the bilinear form
+    SymmetryStructure symmetry_;
 
     /// Dense matrix to store coefficients during \ref assemble()
     ElementMatrix elementMatrix_;
@@ -176,7 +198,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_;
   };
 
 

amdis/BiLinearForm.inc.hpp

@@ -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;
 }
 
 
@@ -63,12 +64,12 @@ assemble(RowLocalView const& rowLocalView, ColLocalView const& colLocalView)
 
 template <class RB, class CB, class T, class Traits>
 void BiLinearForm<RB,CB,T,Traits>::
-assemble(SymmetryStructure symmetry)
+assemble()
 {
   auto rowLocalView = this->rowBasis()->localView();
   auto colLocalView = this->colBasis()->localView();
 
-  this->init(symmetry);
+  this->init();
   for (auto const& element : elements(this->rowBasis()->gridView(), typename Traits::PartitionSet{})) {
     rowLocalView.bind(element);
     if (this->rowBasis() == this->colBasis())

amdis/ProblemStat.inc.hpp

@@ -219,6 +219,10 @@ template <class Traits>
 void ProblemStat<Traits>::createMatricesAndVectors()
 {
   systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
+  std::string symmetryStr = "unknown";
+  Parameters::get(name_ + "->symmetry", symmetryStr);
+  systemMatrix_->setSymmetryStructure(symmetryStr);
+
   solution_ = std::make_shared<SolutionVector>(globalBasis_);
   rhs_ = std::make_shared<SystemVector>(globalBasis_);
 
@@ -476,10 +480,7 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
   t2.reset();
 
   // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
-  std::string symmetryStr = "unknown";
-  Parameters::get(name_ + "->symmetry", symmetryStr);
-
-  systemMatrix_->init(symmetryStructure(symmetryStr));
+  systemMatrix_->init();
   rhs_->init(sizeInfo(*globalBasis_), asmVector);
 
   // statistic about system size

amdis/linearalgebra/MatrixFacade.hpp

@@ -8,7 +8,6 @@
 #include <amdis/common/ConceptsBase.hpp>
 #include <amdis/common/TypeTraits.hpp>
 #include <amdis/functions/NodeIndices.hpp>
-#include <amdis/linearalgebra/SymmetryStructure.hpp>
 #include <amdis/typetree/MultiIndex.hpp>
 
 namespace AMDiS
@@ -40,16 +39,17 @@ 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, SymmetryStructure symmetry = SymmetryStructure::unknown)
+    void init(SparsityPattern const& pattern)
     {
-      impl_.init(pattern, symmetry);
+      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

amdis/linearalgebra/SparsityPattern.hpp

@@ -5,6 +5,7 @@
 #include <dune/istl/matrixindexset.hh>
 
 #include <amdis/common/Index.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 namespace AMDiS
 {
@@ -13,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
     {
@@ -54,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)
     {
-      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);
@@ -79,18 +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)
     {
-      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())) {

amdis/linearalgebra/eigen/MatrixBackend.hpp

@@ -10,7 +10,6 @@
 #include <dune/common/timer.hh>
 
 #include <amdis/Output.hpp>
-#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 namespace AMDiS
 {
@@ -73,16 +72,21 @@ 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, SymmetryStructure symmetry)
+    void init(Pattern const& pattern)
     {
       matrix_.resize(pattern.rows(), pattern.cols());
       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());

amdis/linearalgebra/eigen/MatrixSize.hpp

@@ -4,12 +4,20 @@
 #include <memory>
 
 #include <amdis/common/Index.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 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
     {
@@ -22,16 +30,9 @@ namespace AMDiS
       return cols_;
     }
 
-    template <class RowBasis, class ColBasis>
-    void init(RowBasis const& rowBasis, ColBasis const& colBasis)
-    {
-      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

amdis/linearalgebra/istl/MatrixBackend.hpp

@@ -61,12 +61,18 @@ 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, SymmetryStructure symmetry)
+    void init(Pattern const& pattern)
     {
       pattern.applyTo(matrix_);
-      symmetry_ = symmetry;
+      initialized_ = true;
+    }
+
+    /// Set all entries to zero while keeping the occupation pattern intact
+    void init()
+    {
+      matrix_ = value_type(0);
       initialized_ = true;
     }
 
@@ -100,12 +106,6 @@ namespace AMDiS
           matrix_[rows[i]][cols[j]] += mat[i][j];
     }
 
-
-    SymmetryStructure symmetry() const
-    {
-      return symmetry_;
-    }
-
     std::size_t nnz() const
     {
       return matrix_.nonzeroes();
@@ -116,7 +116,6 @@ namespace AMDiS
     Comm const* comm_;
 
     bool initialized_ = false;
-    SymmetryStructure symmetry_;
   };
 
 } // end namespace AMDiS

amdis/linearalgebra/mtl/MatrixBackend.hpp

@@ -56,19 +56,34 @@ namespace AMDiS
       return matrix_;
     }
 
-    /// Create inserter. Assumes that no inserter is currently active on this matrix.
-    void init(Pattern const& pattern, SymmetryStructure symmetry)
+    /// 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_);
 
+      symmetry_ = pattern.symmetry();
       inserter_ = new Inserter(matrix_, slotSize);
-      symmetry_ = symmetry;
     }
 
+    /// 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()

amdis/linearalgebra/mtl/SlotSize.hpp

 #pragma once
 
+#include <algorithm>
 #include <cassert>
-#include <functional>
-#include <memory>
+
 #include <dune/common/math.hh>
 
 #include <amdis/Observer.hpp>
 #include <amdis/common/Index.hpp>
 #include <amdis/common/Math.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 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
     {
@@ -26,16 +37,19 @@ 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
+    {
+      return estRowSize_;
+    }
+
+    /// Symmetry of the matrix entries
+    SymmetryStructure symmetry() const
     {
-      assert(rows_ > 0 && cols_ > 0);
-      return matrix.nnz() > 0 ? 6*matrix.nnz() / (5*rows_) : estRowSize_;
+      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
@@ -47,15 +61,14 @@ namespace AMDiS
       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);
     }
 
   private:
-    std::size_t rows_ = 0;
-    std::size_t cols_ = 0;
-    std::size_t estRowSize_ = 50;
+    std::size_t rows_;
+    std::size_t cols_;
+    std::size_t estRowSize_;
+    SymmetryStructure symmetry_;
   };
 
 } // end namespace AMDiS

amdis/linearalgebra/petsc/MatrixBackend.hpp

@@ -104,7 +104,7 @@ namespace AMDiS
 
     /// Create and initialize the matrix
     template <class Pattern>
-    void init(Pattern const& pattern, SymmetryStructure symmetry)
+    void init(Pattern const& pattern)
     {
       Dune::Timer t;
 
@@ -123,7 +123,7 @@ namespace AMDiS
       MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
 
       // set symmetry properties of the matrix
-      switch (symmetry) {
+      switch (pattern.symmetry()) {
         case SymmetryStructure::spd:
           MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
           break;
@@ -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()
     {

amdis/linearalgebra/petsc/MatrixNnzStructure.hpp

@@ -4,6 +4,7 @@
 #include <vector>
 
 #include <amdis/common/Index.hpp>
+#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
 #if HAVE_MPI
 #include <amdis/common/parallel/Communicator.hpp>
@@ -15,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
     {
@@ -27,17 +36,21 @@ namespace AMDiS
       return onnz_;
     }
 
-    // Update pattern when basis is updated
+    /// Symmetry of the matrix entries
+    SymmetryStructure symmetry() const
+    {
+      return symmetry_;
+    }
+
+  protected:
+    // Create pattern from basis
     template <class RowBasis, class ColBasis>
     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)
-
-    // 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