Commit 93c69d6e authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

several tests copied from wir-common repository, cleanup of CMakeLists files...

several tests copied from wir-common repository, cleanup of CMakeLists files and GenericIterator modified.
parent 7749968f
......@@ -4,6 +4,7 @@ dune_library_add_sources(dunedec SOURCES
)
install(FILES
CoordsGridFactory.hpp
Dec.hpp
DecGrid.hpp
FlatSharp.hpp
......
......@@ -168,14 +168,14 @@ namespace Dec
template <class T, std::size_t N>
auto max(std::array<T, N> const& x)
{
return aux::accumulate(x, operation::max{});
return aux::accumulate(x, operation::maximum{});
}
/// Minimum over all vector entries
template <class T, std::size_t N>
auto min(std::array<T, N> const& x)
{
return aux::accumulate(x, operation::min{});
return aux::accumulate(x, operation::minimum{});
}
/// Maximum of the absolute values of vector entries
......
......@@ -11,7 +11,6 @@ namespace Dec
namespace concepts
{
#ifndef DOXYGEN
namespace definition
{
template <class E, class T, std::size_t N>
......@@ -31,25 +30,12 @@ namespace Dec
std::false_type is_matrix(...);
} // end namespace definition
#endif
/** \addtogroup concept
* @{
**/
#ifdef HAVE_CONCEPTS
/// \brief Argument is a vector
template <class V>
concept bool Vector = requires(V v) { { definition::is_vector(v) } -> std::true_type };
/// \brief Argument is a matrix
template <class E>
concept bool Matrix = requires(E e) { { definition::is_matrix(e) } -> std::true_type };
#else // HAVE_CONCEPTS
/// \brief Argument is a vector
template <class V>
constexpr bool Vector = decltype(definition::is_vector(std::declval<V>()))::value;
......@@ -58,8 +44,7 @@ namespace Dec
template <class E>
constexpr bool Matrix = decltype(definition::is_matrix(std::declval<E>()))::value;
#endif // HAVE_CONCEPTS
/** @} **/
}
namespace aux
......@@ -72,22 +57,25 @@ namespace Dec
namespace concepts
{
/** \addtogroup concept
* @{
**/
/// \brief Argument is a vector of size N
template <class V, std::size_t N>
CONCEPT bool Vector_size = Vector<V> && Size<V> == N;
constexpr bool Vector_size = Vector<V> && Size<V> == N;
/// \brief Argument is a matrix of size NxM
template <class E, std::size_t N, std::size_t M>
CONCEPT bool Matrix_size = Matrix<E> && E::rows() == N && E::cols() = M;
constexpr bool Matrix_size = Matrix<E> && E::num_rows() == N && E::num_cols() = M;
/// \brief Arguments are vectors of the same size
template <class V0, class V1>
CONCEPT bool Vectors = Vector<V0> && Vector_size<V1, Size<V0>>;
constexpr bool Vectors = Vector<V0> && Vector_size<V1, Size<V0>>;
/// \brief Arguments are matrices of the same size
template <class E0, class E1>
CONCEPT bool Matrices = Matrix<E0> && Matrix_size<E1, E0::rows(), E0::cols()>;
constexpr bool Matrices = Matrix<E0> && Matrix_size<E1, E0::num_rows(), E0::num_cols()>;
/** @} **/
......
......@@ -162,21 +162,21 @@ namespace Dec
} // end namespace aux
template <class E>
struct Collection;
template <class Model, class T, std::size_t N, std::size_t M>
struct Collection<matrix_expression_c<Model, T, N, M>>
{
using Expr = Model;
using value_type = T;
static constexpr std::size_t rows = N;
static constexpr std::size_t cols = M;
static constexpr std::size_t size = N*M;
};
template <class Expr>
struct Collection : Collection<aux::Base_t<Expr>> {};
// template <class E>
// struct Collection;
//
// template <class Model, class T, std::size_t N, std::size_t M>
// struct Collection<matrix_expression_c<Model, T, N, M>>
// {
// using Expr = Model;
// using value_type = T;
//
// static constexpr std::size_t rows = N;
// static constexpr std::size_t cols = M;
// static constexpr std::size_t size = N*M;
// };
//
// template <class Expr>
// struct Collection : Collection<aux::Base_t<Expr>> {};
} // end namespace Dec
......@@ -6,8 +6,11 @@ install(FILES
block.hpp
default.hpp
eigen.hpp
iteration.hpp
multigrid.hpp
petsc.hpp
preconditioner.hpp
smoother.hpp
DESTINATION include/dune/dec/linear_algebra/
)
......@@ -20,7 +23,9 @@ else ()
endif ()
add_subdirectory("block")
add_subdirectory("interface")
add_subdirectory("iteration")
add_subdirectory("krylov")
add_subdirectory("multigrid")
add_subdirectory("preconditioner")
add_subdirectory("smoother")
......@@ -4,6 +4,7 @@
#include <dune/dec/common/Compatible.hpp>
#include <dune/dec/linear_algebra/block.hpp>
#include <dune/dec/linear_algebra/interface/Solver.hpp>
namespace Dec
{
......@@ -22,6 +23,7 @@ namespace Dec
**/
template <class... Solver>
struct BlockJacobi
: public SolverBase<BlockJacobi<Solver...>>
{
static constexpr int N = sizeof...(Solver);
......
install(FILES
Preconditioner.hpp
Smoother.hpp
Solver.hpp
DESTINATION include/dune/dec/linear_algebra/interface/
)
......@@ -8,10 +8,10 @@ namespace Dec
**/
/// Interface of a preconditioner
struct Preconditioner
struct PreconBase
{
/// Prepare the preconditioner, based on the matrix \f$ A \f$
virtual Preconditioner& compute(DOFMatrix<float_type> const& A) = 0;
virtual PreconBase& compute(DOFMatrix<float_type> const& A) = 0;
/// \brief Apply the preconditioner \f$ P \simeq A^{-1} \f$ to the vector \f$ b \f$, i.e.
/// calculate \f$ u = P\cdot b \f$.
......
......@@ -8,7 +8,7 @@ namespace Dec
**/
/// Interface of a smoother
struct Smoother
struct SmootherBase
{
/// Apply one smoothing iteration to the system \f$ A\cdot u = b \f$
virtual void operator()(DOFMatrix<float_type> const& A, DOFVector<float_type>& u, DOFVector<float_type> const& b) const = 0;
......
......@@ -8,24 +8,78 @@ namespace Dec
**/
/// Interface of a linear solver
struct Solver
template <class Impl>
struct SolverBase
{
/// Prepare the solver, based on the matrix \f$ A \f$
virtual Solver& compute(DOFMatrix<float_type> const& A) = 0;
/// \brief Prepare the solver, based on the matrix \f$ A \f$
/**
* \note Should be implemented by derived class.
**/
template <class Matrix>
Impl& compute(Matrix const& A)
{
error_exit("Should be implemented by derived class.");
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial condition is zero.
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial
/// condition is zero.
/**
* Requirement:
* - \ref compute() must be called in advance
**/
template <class VectorB, class VectorU>
void solve(VectorB const& b, VectorU& u) const
{
u.setZero();
solveWithGuessImpl(b, u);
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ with initial
/// condition is zero using iteration.
/**
* Requirement:
* - \ref compute() must be called in advance
**/
virtual void solve(DOFVector<float_type> const& b, DOFVector<float_type>& u) const = 0;
template <class VectorB, class VectorU, class Iteration>
void solve(VectorB const& b, VectorU& u, Iteration& iter) const
{
u.setZero();
solveWithGuess(b, u, iter);
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ .
/// \brief Solve the linear system \f$ A\cdot u = b \f$
/**
* Requirement:
* - \ref compute() must be called in advance
*
* \note Should be implemented by derived class.
**/
template <class VectorB, class VectorU>
void solveWithGuess(VectorB const& b, VectorU& u) const
{
error_exit("Should be implemented by derived class.");
}
/// \brief Solve the linear system \f$ A\cdot u = b \f$ using iteration.
/**
* Requirement:
* - \ref compute() must be called in advance
**/
virtual void solveWithGuess(DOFVector<float_type> const& b, DOFVector<float_type>& u) const = 0;
template <class VectorB, class VectorU, class Iteration>
void solveWithGuess(VectorB const& b, VectorU& u, Iteration& iter) const
{
for (iter.init(b); !iter.finished(); ++iter)
solveWithGuessImpl(b, u);
}
private:
// redirect to implementation
template <class VectorB, class VectorU>
void solveWithGuessImpl(VectorB const& b, VectorU& u) const
{
static_cast<Impl const&>(*this).solveWithGuess(b, u);
}
};
/** @} */
......
#pragma once
/**
* \defgroup iteration Iteration
* \brief Repeated and controlled iteration of (linear) solvers
**/
#include "iteration/BasicIteration.hpp"
#include "iteration/ResidualIteration.hpp"
#include "iteration/GenericIteration.hpp"
#include "iteration/NoisyIteration.hpp"
#include "iteration/IteratedSolver.hpp"
#pragma once
#include <dune/dec/common/ConceptsBase.hpp>
namespace Dec
{
// forward declarations
template <class,int> class DOFMatrix;
template <class> class DOFVector;
/// A basic iteration scheme, returns finihed() == false while i < numIter
/**
* \addtogroup iteration
* @{
**/
/// \brief A basic iteration scheme, returns `finished() == false` while `i < numIter`
template <class Real>
class BasicIteration
{
using Self = BasicIteration;
public:
using real_type = Real;
/// Constructor, stores maximal number of iterations `numIter`.
BasicIteration(int numIter)
......@@ -57,6 +62,48 @@ namespace Dec
int i_ = 0;
};
/** @} */
namespace concepts
{
#ifndef DOXYGEN
namespace definition
{
struct Iteration
{
template <class I, class Real = typename I::real_type>
auto requires_(I&& iter) -> decltype(
concepts::valid_expr(
++iter,
iter++,
int( *iter ),
static_cast<BasicIteration<Real> const&>(iter)
));
};
} // end namespace definition
#endif
/** \addtogroup concept
* @{
**/
/// \brief Argument `I` is an iteration object.
/**
* Requirement:
* - pre-/post-increment operator (go to next iteration)
* - dereferencing returns integer (number of iteration)
* - method `init()` (resets iteration to initial state)
* - method `finished() -> bool` (tests whether to stop the iteration)
**/
template <class I>
constexpr bool Iteration = models<definition::Iteration(I)>;
/** @} */
} // end namespace concepts
} // end namespace Dec
install(FILES
BasicIteration.hpp
GenericIteration.hpp
IteratedSolver.hpp
NoisyIteration.hpp
ResidualIteration.hpp
......
#pragma once
#include <functional>
#include <dune/dec/LinearAlgebra.hpp>
#include "BasicIteration.hpp"
namespace Dec
{
/**
* \addtogroup iteration
* @{
**/
/// \brief Class representing an iteration that watches the residuum of a
/// linear system and breaks if residuum is below some tolerance.
/// The residuum is calculated by a functor passed to the class.
/**
* Requirement:
* - `IterBase` models `concepts::Iteration`.
*
* \see \ref BasicIteration.
**/
template <class IterBase,
REQUIRES( concepts::Iteration<IterBase> )>
class GenericIteration
: public IterBase
{
using Super = IterBase;
public:
using real_type = typename Super::real_type;
/// Constructor, stores function to calc the residuum and initializes the base-class iteration
template <class... Args>
GenericIteration(std::function<real_type(void)> residual, Args&&... args)
: Super(std::forward<Args>(args)...)
, residual_(residual)
{}
/// Calls the \ref Super::finished() with the evaluation of \ref residual_.
template <class... Args>
bool finished(Args&&...) const
{
return Super::finished(residual_());
}
protected:
std::function<real_type(void)> residual_;
};
/** @} */
} // end namespace Dec
......@@ -4,14 +4,18 @@
namespace Dec
{
/// \see Iteration: \ref BasicIteration, \ref ResidualIteration
/**
* \addtogroup linear_solver
* @{
**/
/// \brief Apply a solver repeatedly. \see Iteration: \ref BasicIteration, \ref ResidualIteration
template <class Solver, class Iteration>
struct IteratedSolver
{
template <class Solver_>
IteratedSolver(Solver_&& solver, DOFMatrix<float_type> const& mat, Iteration const& iter)
IteratedSolver(Solver_&& solver, Iteration& iter)
: solver_{std::forward<Solver_>(solver)}
, mat_(mat)
, iter_(iter)
{}
......@@ -32,25 +36,23 @@ namespace Dec
template <class VectorIn, class VectorOut>
void solveWithGuess(VectorIn const& b, VectorOut& u) const
{
Iteration iter{iter_};
for (iter.init(b); !iter.finished(mat_, u, b); ++iter)
for (iter_.init(b); !iter_.finished(); ++iter_)
solver_->solveWithGuess(b, u);
}
private:
SmartRef<Solver> solver_;
DOFMatrix<float_type> const& mat_;
Iteration iter_;
Iteration& iter_;
};
/// Generator for iterated solver
/// Generator for iterated solver, \relates IteratedSolver
template <class Solver, class Iteration>
IteratedSolver<Decay_t<Solver>, Decay_t<Iteration>>
makeIteratedSolver(Solver&& solver, DOFMatrix<float_type> const& mat, Iteration const& iter)
iterated(Solver&& solver, Iteration& iter)
{
return {std::forward<Solver>(solver), mat, iter};
return {std::forward<Solver>(solver), iter};
}
/** @} */
} // end namespace Dec
......@@ -4,6 +4,11 @@
namespace Dec
{
/**
* \addtogroup iteration
* @{
**/
/// An iteration that prints the actual residuum every few iterations
template <class Real>
class CyclicIteration
......@@ -51,4 +56,6 @@ namespace Dec
{}
};
/** @} */
} // end namespace Dec
......@@ -6,7 +6,12 @@
namespace Dec
{
/// Class representing an iteration that watches the residuum of a linear system
/**
* \addtogroup iteration
* @{
**/
/// \brief Class representing an iteration that watches the residuum of a linear system
/// and breaks if residuum is below some tolerance.
template <class Real>
class ResidualIteration
......@@ -33,7 +38,8 @@ namespace Dec
/// Calculates the norm of the vecotr `b`, used as relative value for the break
/// test.
void init(DOFVector<double> const& b)
template <class Vector>
void init(Vector const& b)
{
Super::init(b);
nrmRhs_ = two_norm(b);
......@@ -48,7 +54,8 @@ namespace Dec
}
/// Calculates the residuum of the linear system A*u=b
bool finished(DOFMatrix<Real> const& A, DOFVector<Real> const& u, DOFVector<Real> const& b) const
template <class Matrix, class VectorU, class VectorB>
bool finished(Matrix const& A, VectorU const& u, VectorB const& b) const
{
return finished( residuum(A, u, b) );
}
......@@ -78,5 +85,6 @@ namespace Dec
mutable Real resid_ = 0.0;
};
/** @} */
} // end namespace Dec
install(FILES
bcgs.hpp
cg.hpp
diagonal_pc.hpp
identity_pc.hpp
richardson.hpp
DESTINATION include/dune/dec/linear_algebra/krylov/
)
\ No newline at end of file
)
......@@ -12,6 +12,7 @@
#include "Hierarchy.hpp"
#include <dune/dec/LinearAlgebra.hpp>
#include <dune/dec/linear_algebra/interface/Solver.hpp>
namespace Dec
{
......@@ -43,7 +44,10 @@ namespace Dec
class Transfer,
class CoarseSpaceSolver = Eigen::UmfPackLU<DOFMatrix<float_type>::Matrix>>
class MGSolver
: public SolverBase<MGSolver<LinearOperator, Transfer, CoarseSpaceSolver>>
{
using Super = SolverBase<MGSolver<LinearOperator, Transfer, CoarseSpaceSolver>>;
public:
/// \brief Constructor, initializes incredients of multi-grid method.
......@@ -113,6 +117,9 @@ namespace Dec
int l = maxLevel();
cycle(l, matrix(l), sol, rhs);
}
using Super::solve;
using Super::solveWithGuess;
/// Return the minimal multi-grid level
int minLevel() const { return minLevel_; }
......
install(FILES
Diagonal.hpp
Identity.hpp
DESTINATION include/dune/dec/linear_algebra/preconditioner/