Skip to content
Snippets Groups Projects
Commit 81a54245 authored by Elias Pipping's avatar Elias Pipping
Browse files

Merge branch 'generalized-blockgsstep-rebased' into 'master'

Add generalized BlockGS steps with modular local solvers.

This continues !1 .

See merge request !8
parents 65a885f0 33f01d08
No related branches found
No related tags found
1 merge request!8Add generalized BlockGS steps with modular local solvers.
Pipeline #
...@@ -2,6 +2,7 @@ install(FILES ...@@ -2,6 +2,7 @@ install(FILES
amgstep.hh amgstep.hh
blockgsstep.cc blockgsstep.cc
blockgsstep.hh blockgsstep.hh
blockgssteps.hh
cgstep.cc cgstep.cc
cgstep.hh cgstep.hh
istlseqilu0step.hh istlseqilu0step.hh
......
#ifndef DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
#define DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
#include <algorithm>
#include <functional>
#include <dune/common/parametertree.hh>
#include <dune/solvers/common/staticmatrixtools.hh>
#include <dune/solvers/common/defaultbitvector.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
namespace Dune {
namespace Solvers {
namespace BlockGS {
enum class Direction { FORWARD, BACKWARD, SYMMETRIC };
/**
* @brief Provides support to parse the @Direction enum.
*/
std::istream& operator>>(std::istream& lhs, Direction& t) {
std::string s;
lhs >> s;
std::transform(s.begin(), s.end(), s.begin(), ::tolower);
if (s == "forward")
t = Direction::FORWARD;
else if (s == "backward")
t = Direction::BACKWARD;
else if (s == "symmetric")
t = Direction::SYMMETRIC;
else
lhs.setstate(std::ios_base::failbit);
return lhs;
}
/**
* \brief Iterates over all rows i and updates the i-th component by solving
* the i-th correction problem by means of localSolver.
* \param m Matrix
* \param x Solution vector
* \param b Inhomogeneity
* \param ignore Ignored DOFs
* \param localSolver
* \param direction Iteration direction (e.g. forward, backward, ...)
*/
template <class M, class V, class BitVector, class LocalSolver>
void linearStep(const M& m, V& x, const V& b, const BitVector& ignore,
LocalSolver&& localSolver,
Direction direction = Direction::FORWARD) {
// Note: move-capture requires C++14
auto blockStep = [&, localSolver = std::move(localSolver) ](size_t i) {
const auto& row_i = m[i];
// Compute residual
auto ri = b[i];
for (auto cIt = row_i.begin(); cIt != row_i.end(); ++cIt)
cIt->mmv(x[cIt.index()], ri);
// Update iterate with correction
x[i] += localSolver(row_i[i], std::move(ri), ignore[i]);
};
if (direction != Direction::BACKWARD)
for (size_t i = 0; i < x.size(); ++i)
blockStep(i);
if (direction != Direction::FORWARD)
for (size_t i = x.size(); i > 0; --i)
blockStep(i - 1);
}
/// \brief Convenience method for LinearIterationSteps.
template <class M, class V, class BitVector, class LocalSolver>
void linearStep(LinearIterationStep<M, V, BitVector>& lis,
LocalSolver&& localSolver,
Direction direction = Direction::FORWARD) {
return linearStep(*lis.mat_, *lis.x_, *lis.rhs_, *lis.ignoreNodes_,
std::forward<LocalSolver>(localSolver), direction);
}
/**
* \brief All functions in this namespace are assumed to compute the solution
* to a correction problem. Hence they do not require an initial guess
* (and assume it to be zero) and return the (approximate) solution.
*/
namespace LocalSolvers {
namespace LinearSolvers {
/** \brief Solves a block problem using the MBlock-inherent solve method. See
* e.g. FieldMatrix::solve.
*/
template <class MBlock, class VBlock>
VBlock direct(const MBlock& m, const VBlock& b) {
VBlock x;
m.solve(x, b);
return x;
}
//! \brief Solves a block problem using LDL^t decomposition.
template <class MBlock, class VBlock>
VBlock ldlt(MBlock m, const VBlock& b) {
VBlock x;
StaticMatrix::ldlt(m, m, m);
StaticMatrix::ldltSolve(m, m, b, x);
return x;
}
static constexpr size_t defaultCgMaxIter = 100;
static constexpr double defaultCgTol = 1e-15;
//! \brief Solves a block problem using the conjugate gradient method.
template <class MBlock, class VBlock>
VBlock cg(const MBlock& m, VBlock r, size_t maxIter = defaultCgMaxIter,
double tol = defaultCgTol) {
VBlock x = 0;
VBlock p = r;
auto q = r;
auto rOldSquared = r * r;
auto tolSquared = tol * tol;
for (size_t j = 0; j < maxIter; ++j) {
m.mv(p, q);
if (rOldSquared < tolSquared)
break;
auto alpha = rOldSquared / (q * p);
x.axpy(alpha, p);
r.axpy(-alpha, q);
auto rSquared = r * r;
auto beta = rSquared / rOldSquared;
p *= beta;
p += r;
rOldSquared = rSquared;
}
return x;
}
} // end namespace LinearSolvers
namespace LocalSolverFromLinearSolver {
// Note: move-capture, auto-return and auto-lambda-arguments require C++14
/**
* \brief Applies truncation to the block problem according to the local ignore
* nodes. The system is modified in such a way that for nodes that are to be
* ignored, the corresponding rhs is set to 0 while the corresponding matrix row
* is set to be the appropriate euclidean basis vector, essentially enforcing an
* exact solution for this node to be zero. This leads to zero correction
* contributions for the ignored DOFs.
* \param ignore The global ignore nodes where the local ignore nodes can be
* obtained from.
* \param localSolver The block solver used to solve the (modified) local
* problem.
*/
template <class LinearSolver>
auto truncate(LinearSolver&& linearSolver) {
return [&, localSolver = std::move(linearSolver) ](
const auto& m, const auto& b, const auto& ignore) {
if (ignore.all())
return typename std::result_of<
LinearSolver(decltype(m), decltype(b))>::type{0};
if (ignore.none())
return linearSolver(m, b);
auto mTruncated = m;
auto bTruncated = b;
assert(b.size() == m.N());
assert(m.N() == m.M());
size_t blockSize = b.size();
for (size_t j = 0; j < blockSize; ++j) {
if (not ignore[j])
continue;
bTruncated[j] = 0;
for (size_t k = 0; k < blockSize; ++k)
mTruncated[j][k] = (k == j);
}
return linearSolver(std::move(mTruncated), std::move(bTruncated));
};
}
} // end namespace LocalSolverFromLinearSolver
namespace LocalSolverRegularizer {
static constexpr double defaultDiagRegularizeParameter = 1e-10;
// Note: move-capture, auto-return and auto-lambda-arguments require C++14
/**
* \brief Applies regularization to the diagonal of the block problem, to
* render possibly indefinite problems definite (and hence enable the
* application of more block problem solve techniques).
* \param p Diagonal entries that are non-zero are scaled by (1+p) while
* zero diagonal entries are set to be p.
*/
template <class LocalSolver>
auto diagRegularize(double p, LocalSolver&& localSolver) {
return [ p, localSolver = std::move(localSolver) ](
const auto& m, const auto& b, const auto& ignore) {
if (p == 0)
return localSolver(m, b, ignore);
auto m_regularized = m;
for (size_t j = 0; j < m_regularized.N(); ++j)
if (!ignore[j]) {
if (m_regularized[j][j] == 0)
m_regularized[j][j] = p;
else
m_regularized[j][j] *= 1.0 + p;
}
return localSolver(std::move(m_regularized), b, ignore);
};
}
} // end namespace LocalSolverRegularizer
//! \brief is @LinearSolvers::direct with ignore nodes.
//! \param r Regularization parameter. Set to 0.0 to switch off regularization.
auto direct(double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
return [r](const auto& m, const auto& b, const auto& ignore) {
auto directSolver = LinearSolvers::direct<std::decay_t<decltype(m)>,
std::decay_t<decltype(b)>>;
auto directTruncatedSolver =
LocalSolverFromLinearSolver::truncate(directSolver);
auto directTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r, directTruncatedSolver);
return directTruncatedRegularizedSolver(m, b, ignore);
};
}
//! \brief is @LinearSolvers::ldlt with ignore nodes.
auto ldlt(double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
return [r](const auto& m, const auto& b, const auto& ignore) {
auto ldltSolver = LinearSolvers::ldlt<std::decay_t<decltype(m)>,
std::decay_t<decltype(b)>>;
auto ldltTruncatedSolver =
LocalSolverFromLinearSolver::truncate(ldltSolver);
auto ldltTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r, ldltTruncatedSolver);
return ldltTruncatedRegularizedSolver(m, b, ignore);
};
}
//! \brief is @LinearSolvers::cg with ignore nodes.
auto cg(size_t maxIter = LinearSolvers::defaultCgMaxIter,
double tol = LinearSolvers::defaultCgTol,
double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
return [maxIter, tol, r](const auto& m, const auto& b, const auto& ignore) {
using namespace std::placeholders;
auto cgSolver = std::bind(
LinearSolvers::cg<std::decay_t<decltype(m)>, std::decay_t<decltype(b)>>,
_1, _2, maxIter, tol);
auto cgTruncatedSolver = LocalSolverFromLinearSolver::truncate(cgSolver);
auto cgTruncatedRegularizedSolver =
LocalSolverRegularizer::diagRegularize(r, cgTruncatedSolver);
return cgTruncatedRegularizedSolver(m, b, ignore);
};
}
} // end namespace LocalSolvers
} // end namespace BlockGS
/**
* \brief A Gauss--Seidel-type linear iteration step.
* \param localSolver The solver for the linear block correction problems.
*/
template <class Matrix, class Vector, class BitVector, class LocalSolver>
struct BlockGSStep : public LinearIterationStep<Matrix, Vector, BitVector> {
BlockGSStep(LocalSolver&& localSolver,
BlockGS::Direction direction = BlockGS::Direction::FORWARD)
: localSolver_(localSolver)
, direction_(direction) {}
void iterate() {
assert(this->ignoreNodes_ != nullptr);
assert(this->mat_ != nullptr);
assert(this->x_ != nullptr);
assert(this->rhs_ != nullptr);
BlockGS::linearStep(*this, localSolver_, direction_);
}
private:
LocalSolver localSolver_;
BlockGS::Direction direction_;
};
/**
* @brief The Type enum provides representatives for (new) block Gauss-Seidel
* steps that can be constructed via the @BlockGSStepFactory.
*/
enum class BlockGSType { Direct, LDLt, CG };
/**
* @brief Provides support to parse the @BlockGSType enum.
*/
std::istream& operator>>(std::istream& lhs, BlockGSType& t) {
std::string s;
lhs >> s;
std::transform(s.begin(), s.end(), s.begin(), ::tolower);
if (s == "direct")
t = BlockGSType::Direct;
else if (s == "ldlt")
t = BlockGSType::LDLt;
else if (s == "cg")
t = BlockGSType::CG;
else
lhs.setstate(std::ios_base::failbit);
return lhs;
}
/**
* @brief The BlockGSStepFactory struct allows to conveniently construct various
* (new) block Gauss-Seidel steps through a ParameterTree config.
*/
template <class Matrix, class Vector,
class BitVector = DefaultBitVector_t<Vector>>
struct BlockGSStepFactory {
using BlockGSStepPtr =
std::shared_ptr<LinearIterationStep<Matrix, Vector, BitVector>>;
// TODO add ProjectedBlockGSStep (not linear!)
// TODO add CachedLDLtBlockGSStep
//! \brief Create a BlockGSStep with custom local solver.
template <class LocalSolver>
static auto create(
LocalSolver&& localSolver,
BlockGS::Direction direction = BlockGS::Direction::FORWARD) {
return BlockGSStep<Matrix, Vector, BitVector, LocalSolver>(
std::forward<LocalSolver>(localSolver), direction);
}
template <class LocalSolver>
static auto createPtr(
LocalSolver&& localSolver,
BlockGS::Direction direction = BlockGS::Direction::FORWARD) {
return std::make_shared<
BlockGSStep<Matrix, Vector, BitVector, LocalSolver>>(
std::forward<LocalSolver>(localSolver), direction);
}
static BlockGSStepPtr createPtrFromConfig(const Dune::ParameterTree& config) {
auto type = config.get<BlockGSType>("type");
auto dir = config.get<BlockGS::Direction>("direction",
BlockGS::Direction::FORWARD);
auto reg = config.get("regularize_diag",
BlockGS::LocalSolvers::LocalSolverRegularizer::
defaultDiagRegularizeParameter);
switch (type) {
case BlockGSType::Direct:
return createPtr(BlockGS::LocalSolvers::direct(reg), dir);
case BlockGSType::LDLt:
return createPtr(BlockGS::LocalSolvers::ldlt(reg), dir);
case BlockGSType::CG: {
auto maxIter = config.get(
"maxIter", BlockGS::LocalSolvers::LinearSolvers::defaultCgMaxIter);
auto tol =
config.get("tol", BlockGS::LocalSolvers::LinearSolvers::defaultCgTol);
return createPtr(BlockGS::LocalSolvers::cg(maxIter, tol, reg), dir);
}
default:
DUNE_THROW(Dune::NotImplemented,
"Factory cannot create this BlockGSType.");
}
}
};
} // end namespace Solvers
} // end namespace Dune
#endif // DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
...@@ -13,6 +13,7 @@ ...@@ -13,6 +13,7 @@
#include <dune/solvers/iterationsteps/blockgsstep.hh> #include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/projectedblockgsstep.hh> #include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
#include <dune/solvers/iterationsteps/semidefiniteblockgsstep.hh> #include <dune/solvers/iterationsteps/semidefiniteblockgsstep.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include "common.hh" #include "common.hh"
...@@ -22,11 +23,9 @@ ...@@ -22,11 +23,9 @@
* is solved correctly for a random rhs by a LoopSolver employing * is solved correctly for a random rhs by a LoopSolver employing
* a GSStep. * a GSStep.
*/ */
struct GSTestSuite struct GSTestSuite {
{
template <class GridType> template <class GridType>
bool check(const GridType& grid) bool check(const GridType& grid) {
{
bool passed = true; bool passed = true;
using Problem = SymmetricSampleProblem<1, typename GridType::LevelGridView>; using Problem = SymmetricSampleProblem<1, typename GridType::LevelGridView>;
...@@ -37,6 +36,17 @@ struct GSTestSuite ...@@ -37,6 +36,17 @@ struct GSTestSuite
using Matrix = typename Problem::Matrix; using Matrix = typename Problem::Matrix;
using Vector = typename Problem::Vector; using Vector = typename Problem::Vector;
enum { BlockSize = Vector::block_type::dimension };
// create some obstacle
size_t size = p.u.size();
std::vector<BoxConstraint<double, BlockSize>> obstacles(size);
for (size_t i = 0; i < size; ++i) {
for (size_t j = 0; j < BlockSize; ++j) {
obstacles[i].lower(j) = j;
obstacles[i].upper(j) = j + 1;
}
}
using LoopSolver = ::LoopSolver<Vector>; using LoopSolver = ::LoopSolver<Vector>;
using Step = LinearIterationStep<Matrix, Vector>; using Step = LinearIterationStep<Matrix, Vector>;
...@@ -46,8 +56,8 @@ struct GSTestSuite ...@@ -46,8 +56,8 @@ struct GSTestSuite
step->setProblem(p.A, u_copy, p.rhs); step->setProblem(p.A, u_copy, p.rhs);
step->ignoreNodes_ = &p.ignore; step->ignoreNodes_ = &p.ignore;
LoopSolver solver(step, maxIterations, stepTol, LoopSolver solver(step, maxIterations, stepTol, &p.energyNorm, verbosity,
&p.energyNorm, verbosity, relativeErrors); relativeErrors);
solver.check(); solver.check();
solver.preprocess(); solver.preprocess();
solver.solve(); solver.solve();
...@@ -55,11 +65,12 @@ struct GSTestSuite ...@@ -55,11 +65,12 @@ struct GSTestSuite
return u_copy; return u_copy;
}; };
auto analyse = [&](const Vector& v, std::string testCase, double tolerance) { auto analyse =
[&](const Vector& v, std::string testCase, double tolerance) {
const auto normDiff = p.energyNorm.diff(p.u_ex, v); const auto normDiff = p.energyNorm.diff(p.u_ex, v);
const std::string s = Dune::formatString("%*s", -60, testCase.c_str()); const std::string s = Dune::formatString("%*s", -60, testCase.c_str());
std::cout << "error for solver \"" << s std::cout << "error for solver \"" << s << "\": " << normDiff
<< "\": " << normDiff << std::endl; << std::endl;
if (normDiff > tolerance) { if (normDiff > tolerance) {
std::cerr << "Error too large for solver \"" << testCase << "\"." std::cerr << "Error too large for solver \"" << testCase << "\"."
<< std::endl; << std::endl;
...@@ -73,25 +84,63 @@ struct GSTestSuite ...@@ -73,25 +84,63 @@ struct GSTestSuite
passed = passed and analyse(result, name, 1e-7); passed = passed and analyse(result, name, 1e-7);
}; };
auto diffTest = [&](Step* step1, std::string name1, Step* step2,
std::string name2, double maxDiff = 0.0) {
auto result1 = solve(step1, 0.0, 5);
auto result2 = solve(step2, 0.0, 5);
result1 -= result2;
if (result1.two_norm() > maxDiff) {
passed = false;
std::cerr << name1 << " and " << name2 << " differ with norm "
<< result1.two_norm() << "." << std::endl;
return;
}
std::cout << Dune::formatString("%*s", -60, name1.c_str()) << " and "
<< Dune::formatString("%*s", -60, name2.c_str()) << " coincide"
<< ((maxDiff > 0.0)
? " up to " + Dune::formatString("%e", maxDiff)
: ".") << std::endl;
passed = passed and true;
};
{ {
BlockGSStep<Matrix, Vector> gsStep; BlockGSStep<Matrix, Vector> gsStep;
test(&gsStep, "BlockGS"); test(&gsStep, "BlockGS");
} }
{
BlockGSStep<Matrix, Vector> gsStep1;
auto gsStep2 =
Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
diffTest(&gsStep1, "BlockGS", &gsStep2, "Dune::Solvers::DirectBlockGS");
}
// test regularizable block GS // test regularizable block GS
for(auto reg : {false, true}) { for (auto reg : {false, true}) {
for(auto type : {"Direct", "LDLt", "CG"}) { for (auto type : {"Direct", "LDLt", "CG"}) {
Dune::ParameterTree config; Dune::ParameterTree config;
config["type"] = type; config["type"] = type;
config["regularize_diagonal"] = reg ? "1" : "0"; // for SemiDefiniteBlockGS config["regularize_diagonal"] =
config["regularize_diag"] = reg ? "1e-10" : "0.0"; // for Dune::Solvers::*BlockGS reg ? "1" : "0"; // for SemiDefiniteBlockGS
config["maxIter"] = "100"; // for Dune::Solvers::CGBlockGS config["regularize_diag"] =
config["tol"] = "1e-10"; // for Dune::Solvers::CGBlockGS reg ? "1e-10" : "0.0"; // for Dune::Solvers::*BlockGS
config["cg_maxIter"] = "100"; // for SemiDefiniteBlockGS with CG
config["cg_tol"] = "1e-10"; // for Dune::Solvers::CGBlockGS with CG
config["maxIter"] = "100"; // for Dune::Solvers::CGBlockGS
config["tol"] = "1e-10"; // for Dune::Solvers::CGBlockGS
{ {
SemiDefiniteBlockGSStep<Matrix, Vector> gsStep(config); SemiDefiniteBlockGSStep<Matrix, Vector> gsStep1(config);
test(&gsStep, Dune::formatString("SemiDefiniteBlockGS %s %s", auto gsString1 = Dune::formatString("SemiDefiniteBlockGS %s %s", type,
type, reg ? "regularized" : "")); reg ? "regularized" : "");
auto gsStep2Ptr =
Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtrFromConfig(config);
auto gsString2 = Dune::formatString("Dune::Solvers::BlockGS(%s) %s",
type, reg ? "regularized" : "");
test(&gsStep1, gsString1);
test(gsStep2Ptr.get(), gsString2);
diffTest(&gsStep1, gsString1, gsStep2Ptr.get(), gsString2, 1e-7);
} }
} }
} }
...@@ -107,14 +156,6 @@ struct GSTestSuite ...@@ -107,14 +156,6 @@ struct GSTestSuite
test(&gsStep, "ProjectedBlockGS free"); test(&gsStep, "ProjectedBlockGS free");
} }
hasObstacle.setAll(); hasObstacle.setAll();
std::vector<typename GSStep::Obstacle> obstacles(size);
for(size_t i=0; i<size; ++i) {
hasObstacle[i] = true;
for(size_t j=0; j<Vector::block_type::dimension; ++j) {
obstacles[i].lower(j) = j;
obstacles[i].upper(j) = j+1;
}
}
{ {
ProjectedBlockGSStep<Matrix, Vector> gsStep; ProjectedBlockGSStep<Matrix, Vector> gsStep;
gsStep.hasObstacle_ = &hasObstacle; gsStep.hasObstacle_ = &hasObstacle;
...@@ -127,8 +168,7 @@ struct GSTestSuite ...@@ -127,8 +168,7 @@ struct GSTestSuite
} }
}; };
int main(int argc, char** argv) int main(int argc, char** argv) {
{
Dune::MPIHelper::instance(argc, argv); Dune::MPIHelper::instance(argc, argv);
bool passed(true); bool passed(true);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment