diff --git a/dune/solvers/iterationsteps/CMakeLists.txt b/dune/solvers/iterationsteps/CMakeLists.txt index f4e462c3c216ee2ea3316086f22fbe273091a4a5..b1d9ec1f6aa201314580eb532b3f6309ef5a785a 100644 --- a/dune/solvers/iterationsteps/CMakeLists.txt +++ b/dune/solvers/iterationsteps/CMakeLists.txt @@ -2,6 +2,7 @@ install(FILES amgstep.hh blockgsstep.cc blockgsstep.hh + blockgssteps.hh cgstep.cc cgstep.hh istlseqilu0step.hh diff --git a/dune/solvers/iterationsteps/blockgssteps.hh b/dune/solvers/iterationsteps/blockgssteps.hh new file mode 100644 index 0000000000000000000000000000000000000000..d26afeec22e216f9fdc8737d026b4cfa633f7092 --- /dev/null +++ b/dune/solvers/iterationsteps/blockgssteps.hh @@ -0,0 +1,372 @@ +#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() && 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 how to solve 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(it.ignoreNodes_ != nullptr); + assert(it.mat_ != nullptr); + assert(it.x_ != nullptr); + assert(it.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 diff --git a/dune/solvers/test/gssteptest.cc b/dune/solvers/test/gssteptest.cc index ab4b5b5327b7ae101bb378d8519968e118544bf6..89ce86973e6711fd9d552b54062a2cc7c58504b4 100644 --- a/dune/solvers/test/gssteptest.cc +++ b/dune/solvers/test/gssteptest.cc @@ -13,6 +13,7 @@ #include <dune/solvers/iterationsteps/blockgsstep.hh> #include <dune/solvers/iterationsteps/projectedblockgsstep.hh> #include <dune/solvers/iterationsteps/semidefiniteblockgsstep.hh> +#include <dune/solvers/iterationsteps/blockgssteps.hh> #include "common.hh" @@ -22,11 +23,9 @@ * is solved correctly for a random rhs by a LoopSolver employing * a GSStep. */ -struct GSTestSuite -{ +struct GSTestSuite { template <class GridType> - bool check(const GridType& grid) - { + bool check(const GridType& grid) { bool passed = true; using Problem = SymmetricSampleProblem<1, typename GridType::LevelGridView>; @@ -37,6 +36,17 @@ struct GSTestSuite using Matrix = typename Problem::Matrix; 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 Step = LinearIterationStep<Matrix, Vector>; @@ -46,8 +56,8 @@ struct GSTestSuite step->setProblem(p.A, u_copy, p.rhs); step->ignoreNodes_ = &p.ignore; - LoopSolver solver(step, maxIterations, stepTol, - &p.energyNorm, verbosity, relativeErrors); + LoopSolver solver(step, maxIterations, stepTol, &p.energyNorm, verbosity, + relativeErrors); solver.check(); solver.preprocess(); solver.solve(); @@ -55,11 +65,12 @@ struct GSTestSuite 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 std::string s = Dune::formatString("%*s", -60, testCase.c_str()); - std::cout << "error for solver \"" << s - << "\": " << normDiff << std::endl; + std::cout << "error for solver \"" << s << "\": " << normDiff + << std::endl; if (normDiff > tolerance) { std::cerr << "Error too large for solver \"" << testCase << "\"." << std::endl; @@ -73,25 +84,63 @@ struct GSTestSuite 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; 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 - for(auto reg : {false, true}) { - for(auto type : {"Direct", "LDLt", "CG"}) { + for (auto reg : {false, true}) { + for (auto type : {"Direct", "LDLt", "CG"}) { Dune::ParameterTree config; config["type"] = type; - config["regularize_diagonal"] = reg ? "1" : "0"; // for SemiDefiniteBlockGS - config["regularize_diag"] = reg ? "1e-10" : "0.0"; // for Dune::Solvers::*BlockGS - config["maxIter"] = "100"; // for Dune::Solvers::CGBlockGS - config["tol"] = "1e-10"; // for Dune::Solvers::CGBlockGS + config["regularize_diagonal"] = + reg ? "1" : "0"; // for SemiDefiniteBlockGS + config["regularize_diag"] = + 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); - test(&gsStep, Dune::formatString("SemiDefiniteBlockGS %s %s", - type, reg ? "regularized" : "")); + 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 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 test(&gsStep, "ProjectedBlockGS free"); } 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; gsStep.hasObstacle_ = &hasObstacle; @@ -127,8 +168,7 @@ struct GSTestSuite } }; -int main(int argc, char** argv) -{ +int main(int argc, char** argv) { Dune::MPIHelper::instance(argc, argv); bool passed(true);