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

Move SymmetryStructure from MatrixFacade::init param to Pattern::init param

parent e53f5c6f
...@@ -56,7 +56,7 @@ namespace AMDiS ...@@ -56,7 +56,7 @@ namespace AMDiS
, rowBasis_(rowBasis) , rowBasis_(rowBasis)
, colBasis_(colBasis) , colBasis_(colBasis)
{ {
pattern_.init(*rowBasis_, *colBasis_); pattern_.init(*rowBasis_, *colBasis_, symmetry_);
operators_.init(*rowBasis_, *colBasis_); operators_.init(*rowBasis_, *colBasis_);
auto const rowSize = rowBasis_->localView().maxSize(); auto const rowSize = rowBasis_->localView().maxSize();
...@@ -150,7 +150,7 @@ namespace AMDiS ...@@ -150,7 +150,7 @@ namespace AMDiS
void init() void init()
{ {
Super::init(pattern_, symmetry_); Super::init(pattern_);
} }
void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); } void updateImpl(event::adapt e, index_t<0> i) override { updateImpl2(e,i); }
...@@ -165,7 +165,7 @@ namespace AMDiS ...@@ -165,7 +165,7 @@ namespace AMDiS
updateCounter_.set(I); updateCounter_.set(I);
if (updateCounter_.all()) { if (updateCounter_.all()) {
pattern_.init(*rowBasis_, *colBasis_); pattern_.init(*rowBasis_, *colBasis_, symmetry_);
updateCounter_.reset(); updateCounter_.reset();
} }
} }
......
...@@ -8,7 +8,6 @@ ...@@ -8,7 +8,6 @@
#include <amdis/common/ConceptsBase.hpp> #include <amdis/common/ConceptsBase.hpp>
#include <amdis/common/TypeTraits.hpp> #include <amdis/common/TypeTraits.hpp>
#include <amdis/functions/NodeIndices.hpp> #include <amdis/functions/NodeIndices.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/typetree/MultiIndex.hpp> #include <amdis/typetree/MultiIndex.hpp>
namespace AMDiS namespace AMDiS
...@@ -47,9 +46,9 @@ namespace AMDiS ...@@ -47,9 +46,9 @@ namespace AMDiS
* See \ref SymmetryStructure. * See \ref SymmetryStructure.
**/ **/
template <class SparsityPattern> template <class SparsityPattern>
void init(SparsityPattern const& pattern, SymmetryStructure symmetry = SymmetryStructure::unknown) void init(SparsityPattern const& pattern)
{ {
impl_.init(pattern, symmetry); impl_.init(pattern);
} }
/// Finish the matrix insertion, e.g. cleanup or final insertion /// Finish the matrix insertion, e.g. cleanup or final insertion
......
...@@ -5,6 +5,7 @@ ...@@ -5,6 +5,7 @@
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
#include <amdis/common/Index.hpp> #include <amdis/common/Index.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
namespace AMDiS namespace AMDiS
{ {
...@@ -56,7 +57,7 @@ namespace AMDiS ...@@ -56,7 +57,7 @@ namespace AMDiS
// Update pattern when basis is updated. This method is called if rowBasis == colBasis. // Update pattern when basis is updated. This method is called if rowBasis == colBasis.
template <class Basis> template <class Basis>
void init(Basis const& basis) void init(Basis const& basis, SymmetryStructure symmetry = SymmetryStructure::unknown)
{ {
rows_ = basis.dimension(); rows_ = basis.dimension();
cols_ = rows_; cols_ = rows_;
...@@ -81,7 +82,8 @@ namespace AMDiS ...@@ -81,7 +82,8 @@ namespace AMDiS
// Update pattern when basis is updated. This method is called if rowBasis != colBasis. // Update pattern when basis is updated. This method is called if rowBasis != colBasis.
template <class RowBasis, class ColBasis> template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis) void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
{ {
if (uintptr_t(&rowBasis) == uintptr_t(&colBasis)) if (uintptr_t(&rowBasis) == uintptr_t(&colBasis))
return init(rowBasis); return init(rowBasis);
......
...@@ -10,7 +10,6 @@ ...@@ -10,7 +10,6 @@
#include <dune/common/timer.hh> #include <dune/common/timer.hh>
#include <amdis/Output.hpp> #include <amdis/Output.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
namespace AMDiS namespace AMDiS
{ {
...@@ -75,7 +74,7 @@ namespace AMDiS ...@@ -75,7 +74,7 @@ namespace AMDiS
/// Create inserter. Assumes that no inserter is currently active on this matrix. /// Create inserter. Assumes that no inserter is currently active on this matrix.
template <class Pattern> template <class Pattern>
void init(Pattern const& pattern, SymmetryStructure symmetry) void init(Pattern const& pattern)
{ {
matrix_.resize(pattern.rows(), pattern.cols()); matrix_.resize(pattern.rows(), pattern.cols());
matrix_.setZero(); matrix_.setZero();
......
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#include <memory> #include <memory>
#include <amdis/common/Index.hpp> #include <amdis/common/Index.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
namespace AMDiS namespace AMDiS
{ {
...@@ -23,7 +24,8 @@ namespace AMDiS ...@@ -23,7 +24,8 @@ namespace AMDiS
} }
template <class RowBasis, class ColBasis> template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis) void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
{ {
rows_ = rowBasis.dimension(); rows_ = rowBasis.dimension();
cols_ = colBasis.dimension(); cols_ = colBasis.dimension();
......
...@@ -63,10 +63,9 @@ namespace AMDiS ...@@ -63,10 +63,9 @@ namespace AMDiS
/// create occupation pattern and apply it to the matrix /// create occupation pattern and apply it to the matrix
template <class Pattern> template <class Pattern>
void init(Pattern const& pattern, SymmetryStructure symmetry) void init(Pattern const& pattern)
{ {
pattern.applyTo(matrix_); pattern.applyTo(matrix_);
symmetry_ = symmetry;
initialized_ = true; initialized_ = true;
} }
...@@ -100,12 +99,6 @@ namespace AMDiS ...@@ -100,12 +99,6 @@ namespace AMDiS
matrix_[rows[i]][cols[j]] += mat[i][j]; matrix_[rows[i]][cols[j]] += mat[i][j];
} }
SymmetryStructure symmetry() const
{
return symmetry_;
}
std::size_t nnz() const std::size_t nnz() const
{ {
return matrix_.nonzeroes(); return matrix_.nonzeroes();
...@@ -116,7 +109,6 @@ namespace AMDiS ...@@ -116,7 +109,6 @@ namespace AMDiS
Comm const* comm_; Comm const* comm_;
bool initialized_ = false; bool initialized_ = false;
SymmetryStructure symmetry_;
}; };
} // end namespace AMDiS } // end namespace AMDiS
...@@ -57,7 +57,7 @@ namespace AMDiS ...@@ -57,7 +57,7 @@ namespace AMDiS
} }
/// Create inserter. Assumes that no inserter is currently active on this matrix. /// Create inserter. Assumes that no inserter is currently active on this matrix.
void init(Pattern const& pattern, SymmetryStructure symmetry) void init(Pattern const& pattern)
{ {
test_exit(!inserter_, "Matrix already in insertion mode!"); test_exit(!inserter_, "Matrix already in insertion mode!");
...@@ -66,7 +66,7 @@ namespace AMDiS ...@@ -66,7 +66,7 @@ namespace AMDiS
set_to_zero(matrix_); set_to_zero(matrix_);
inserter_ = new Inserter(matrix_, slotSize); inserter_ = new Inserter(matrix_, slotSize);
symmetry_ = symmetry; symmetry_ = pattern.symmetry();
} }
/// Delete inserter -> finish insertion. Must be called in order to fill the /// Delete inserter -> finish insertion. Must be called in order to fill the
......
...@@ -8,6 +8,7 @@ ...@@ -8,6 +8,7 @@
#include <amdis/Observer.hpp> #include <amdis/Observer.hpp>
#include <amdis/common/Index.hpp> #include <amdis/common/Index.hpp>
#include <amdis/common/Math.hpp> #include <amdis/common/Math.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
namespace AMDiS namespace AMDiS
{ {
...@@ -36,11 +37,18 @@ namespace AMDiS ...@@ -36,11 +37,18 @@ namespace AMDiS
return matrix.nnz() > 0 ? 6*matrix.nnz() / (5*rows_) : estRowSize_; return matrix.nnz() > 0 ? 6*matrix.nnz() / (5*rows_) : estRowSize_;
} }
/// Symmetry of the matrix entries
SymmetryStructure symmetry() const
{
return symmetry_;
}
// estimate the number of columns by multiplying the maximal node size with the // 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 element surrounding a vertex. This number is approximated by the
// number of simplices surrounding a vertex in a kuhn triangulation // number of simplices surrounding a vertex in a kuhn triangulation
template <class RowBasis, class ColBasis> template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis) void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown)
{ {
using GridView = typename RowBasis::GridView; using GridView = typename RowBasis::GridView;
// TODO(FM): Simplify using Dune::factorial with 2.7 // TODO(FM): Simplify using Dune::factorial with 2.7
...@@ -50,12 +58,14 @@ namespace AMDiS ...@@ -50,12 +58,14 @@ namespace AMDiS
rows_ = rowBasis.dimension(); rows_ = rowBasis.dimension();
cols_ = colBasis.dimension(); cols_ = colBasis.dimension();
estRowSize_ = std::min(cols_, colBasis.localView().maxSize() * surrounding); estRowSize_ = std::min(cols_, colBasis.localView().maxSize() * surrounding);
symmetry_ = symmetry;
} }
private: private:
std::size_t rows_ = 0; std::size_t rows_ = 0;
std::size_t cols_ = 0; std::size_t cols_ = 0;
std::size_t estRowSize_ = 50; std::size_t estRowSize_ = 50;
SymmetryStructure symmetry_ = SymmetryStructure::unknown;
}; };
} // end namespace AMDiS } // end namespace AMDiS
...@@ -104,7 +104,7 @@ namespace AMDiS ...@@ -104,7 +104,7 @@ namespace AMDiS
/// Create and initialize the matrix /// Create and initialize the matrix
template <class Pattern> template <class Pattern>
void init(Pattern const& pattern, SymmetryStructure symmetry) void init(Pattern const& pattern)
{ {
Dune::Timer t; Dune::Timer t;
...@@ -123,7 +123,7 @@ namespace AMDiS ...@@ -123,7 +123,7 @@ namespace AMDiS
MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); MatSetOption(matrix_, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
// set symmetry properties of the matrix // set symmetry properties of the matrix
switch (symmetry) { switch (pattern.symmetry()) {
case SymmetryStructure::spd: case SymmetryStructure::spd:
MatSetOption(matrix_, MAT_SPD, PETSC_TRUE); MatSetOption(matrix_, MAT_SPD, PETSC_TRUE);
break; break;
......
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#include <vector> #include <vector>
#include <amdis/common/Index.hpp> #include <amdis/common/Index.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#if HAVE_MPI #if HAVE_MPI
#include <amdis/common/parallel/Communicator.hpp> #include <amdis/common/parallel/Communicator.hpp>
...@@ -27,13 +28,21 @@ namespace AMDiS ...@@ -27,13 +28,21 @@ namespace AMDiS
return onnz_; return onnz_;
} }
/// Symmetry of the matrix entries
SymmetryStructure symmetry() const
{
return symmetry_;
}
// Update pattern when basis is updated // Update pattern when basis is updated
template <class RowBasis, class ColBasis> template <class RowBasis, class ColBasis>
void init(RowBasis const& rowBasis, ColBasis const& colBasis); void init(RowBasis const& rowBasis, ColBasis const& colBasis,
SymmetryStructure symmetry = SymmetryStructure::unknown);
private: private:
std::vector<PetscInt> dnnz_; //< number of nonzeros in the diagonal part (owner part) 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) 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. // use a grid change index to decide whether to update the matrix pattern or not.
unsigned long changeIndex_ = 0; unsigned long changeIndex_ = 0;
......
...@@ -17,10 +17,12 @@ ...@@ -17,10 +17,12 @@
namespace AMDiS { namespace AMDiS {
template <class RowBasis, class ColBasis> template <class RowBasis, class ColBasis>
void MatrixNnzStructure::init(RowBasis const& rowBasis, ColBasis const& /*colBasis*/) void MatrixNnzStructure::init(RowBasis const& rowBasis, ColBasis const& /*colBasis*/,
SymmetryStructure symmetry)
{ {
Dune::Timer t; Dune::Timer t;
symmetry_ = symmetry;
auto const& basis = rowBasis; // TODO: generalize to row != col basis auto const& basis = rowBasis; // TODO: generalize to row != col basis
auto const& dofMap = basis.communicator().dofMap(); auto const& dofMap = basis.communicator().dofMap();
std::size_t localSize = dofMap.localSize(); std::size_t localSize = dofMap.localSize();
......
Supports Markdown
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