diff --git a/dune/solvers/iterationsteps/blockgssteps.hh b/dune/solvers/iterationsteps/blockgssteps.hh index b00584c52923bb29fc87d9966658f348f417bc74..7b351987915e5ecb8282c3d17a6eab0821ef8adf 100644 --- a/dune/solvers/iterationsteps/blockgssteps.hh +++ b/dune/solvers/iterationsteps/blockgssteps.hh @@ -121,6 +121,29 @@ VBlock cg(const MBlock& m, VBlock r, size_t maxIter = defaultCgMaxIter, return x; } +static constexpr double defaultGsTol = 9e-11; + +//! \brief Solves a block problem using the conjugate gradient method. +template <class MBlock, class VBlock> +VBlock gs(const MBlock& m, const VBlock& b, double tol = defaultGsTol) { + VBlock x(0); + for (size_t i = 0; i < m.N(); ++i) { + const auto& mi = m[i]; + const auto& mii = mi[i]; + if (std::abs(mii) > tol) { + x[i] = b[i]; + const auto& end = mi.end(); + for (auto it = mi.begin(); it != end; ++it) { + auto j = it.index(); + if (j != i) + x[i] -= (*it) * x[j]; + } + x[i] /= mii; + } + } + return x; +} + } // end namespace LinearSolvers namespace LocalSolverFromLinearSolver { @@ -244,6 +267,24 @@ auto cg(size_t maxIter = LinearSolvers::defaultCgMaxIter, }; } +//! \brief is @LinearSolvers::direct with ignore nodes. +//! \param r Regularization parameter. Set to 0.0 to switch off regularization. +template <typename dummy = void> +auto gs(double tol = LinearSolvers::defaultGsTol, + double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) { + return [tol, r](const auto& m, const auto& b, const auto& ignore) { + using namespace std::placeholders; + auto gsSolver = std::bind( + LinearSolvers::gs<std::decay_t<decltype(m)>, std::decay_t<decltype(b)>>, + _1, _2, tol); + auto gsTruncatedSolver = + LocalSolverFromLinearSolver::truncateSymmetrically(gsSolver); + auto gsTruncatedRegularizedSolver = + LocalSolverRegularizer::diagRegularize(r, gsTruncatedSolver); + return gsTruncatedRegularizedSolver(m, b, ignore); + }; +} + } // end namespace LocalSolvers } // end namespace BlockGS diff --git a/dune/solvers/test/gssteptest.cc b/dune/solvers/test/gssteptest.cc index d6e56b98990ae9ce89a23a8a42790575febc4341..c6131e6d8e2f5134bd41ceca9f2da70adfaa668d 100644 --- a/dune/solvers/test/gssteptest.cc +++ b/dune/solvers/test/gssteptest.cc @@ -14,6 +14,7 @@ #include <dune/solvers/iterationsteps/projectedblockgsstep.hh> #include <dune/solvers/iterationsteps/semidefiniteblockgsstep.hh> #include <dune/solvers/iterationsteps/blockgssteps.hh> +#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh> #include "common.hh" @@ -108,10 +109,18 @@ struct GSTestSuite { } { 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"); + auto gsStep2 = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create( + Dune::Solvers::BlockGS::LocalSolvers::direct(0.0)); + diffTest(&gsStep1, "BlockGS", &gsStep2, "Dune::Solvers::BlockGS(Direct)"); + } + { + TruncatedBlockGSStep<Matrix, Vector> gsStep1(1); + auto gsStep2 = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create( + Dune::Solvers::BlockGS::LocalSolvers::gs(0.0, 0.0)); + test(&gsStep1, "TruncatedBlockGS recursive"); + test(&gsStep2, "Dune::Solvers::BlockGS(GS)"); + diffTest(&gsStep1, "TruncatedBlockGS recursive", &gsStep2, + "Dune::Solvers::BlockGS(GS)", 1e-9); } // test regularizable block GS @@ -132,8 +141,8 @@ struct GSTestSuite { SemiDefiniteBlockGSStep<Matrix, Vector> gsStep1(config); auto gsString1 = Dune::formatString("SemiDefiniteBlockGS %s %s", type, reg ? "regularized" : ""); - auto gsStep2Ptr = - Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtrFromConfig(config); + auto gsStep2Ptr = Dune::Solvers::BlockGSStepFactory< + Matrix, Vector>::createPtrFromConfig(config); auto gsString2 = Dune::formatString("Dune::Solvers::BlockGS(%s) %s", type, reg ? "regularized" : ""); test(&gsStep1, gsString1);