multigrid.hpp

// SPDX-FileCopyrightText: 2017 - 2026 The Ginkgo authors
//
// SPDX-License-Identifier: BSD-3-Clause
 
#ifndef GKO_PUBLIC_CORE_SOLVER_MULTIGRID_HPP_
#define GKO_PUBLIC_CORE_SOLVER_MULTIGRID_HPP_
 
 
#include <complex>
#include <functional>
#include <memory>
#include <utility>
 
#include <ginkgo/core/base/abstract_factory.hpp>
#include <ginkgo/core/base/array.hpp>
#include <ginkgo/core/base/exception_helpers.hpp>
#include <ginkgo/core/base/lin_op.hpp>
#include <ginkgo/core/base/math.hpp>
#include <ginkgo/core/base/precision_dispatch.hpp>
#include <ginkgo/core/base/types.hpp>
#include <ginkgo/core/config/config.hpp>
#include <ginkgo/core/config/registry.hpp>
#include <ginkgo/core/config/type_descriptor.hpp>
#include <ginkgo/core/log/logger.hpp>
#include <ginkgo/core/matrix/dense.hpp>
#include <ginkgo/core/matrix/identity.hpp>
#include <ginkgo/core/multigrid/multigrid_level.hpp>
#include <ginkgo/core/solver/solver_base.hpp>
#include <ginkgo/core/stop/combined.hpp>
#include <ginkgo/core/stop/criterion.hpp>
 
 
namespace gko {
namespace solver {
 
 
namespace multigrid {
 
 
enum class cycle { v, f, w };
 
 
enum class mid_smooth_type { both, post_smoother, pre_smoother, standalone };
 
 
namespace detail {
 
 
// It should only be used internally
class MultigridState;
 
 
}  // namespace detail
}  // namespace multigrid
 
 
class Multigrid : public LinOp,
                  public EnableSolverBase<Multigrid>,
                  public EnableIterativeBase<Multigrid>,
                  public EnableApplyWithInitialGuess<Multigrid> {
    friend class EnableApplyWithInitialGuess<Multigrid>;
 
public:
    bool apply_uses_initial_guess() const override
    {
        return this->get_default_initial_guess() ==
               initial_guess_mode::provided;
    }
 
    std::vector<std::shared_ptr<const gko::multigrid::MultigridLevel>>
    get_mg_level_list() const
    {
        return mg_level_list_;
    }
 
    std::vector<std::shared_ptr<const LinOp>> get_pre_smoother_list() const
    {
        return pre_smoother_list_;
    }
 
    std::vector<std::shared_ptr<const LinOp>> get_mid_smoother_list() const
    {
        return mid_smoother_list_;
    }
 
    std::vector<std::shared_ptr<const LinOp>> get_post_smoother_list() const
    {
        return post_smoother_list_;
    }
 
    std::shared_ptr<const LinOp> get_coarsest_solver() const
    {
        return coarsest_solver_;
    }
 
    multigrid::cycle get_cycle() const { return parameters_.cycle; }
 
    void set_cycle(multigrid::cycle cycle) { parameters_.cycle = cycle; }
 
 
    class Factory;
 
    struct parameters_type
        : public enable_iterative_solver_factory_parameters<parameters_type,
                                                            Factory> {
        std::vector<std::shared_ptr<const LinOpFactory>>
            GKO_DEFERRED_FACTORY_VECTOR_PARAMETER(mg_level);
 
        std::function<size_type(const size_type, const LinOp*)>
            GKO_FACTORY_PARAMETER_SCALAR(level_selector, nullptr);
 
        std::vector<std::shared_ptr<const LinOpFactory>>
            GKO_DEFERRED_FACTORY_VECTOR_PARAMETER(pre_smoother);
 
        std::vector<std::shared_ptr<const LinOpFactory>>
            GKO_DEFERRED_FACTORY_VECTOR_PARAMETER(post_smoother);
 
        std::vector<std::shared_ptr<const LinOpFactory>>
            GKO_DEFERRED_FACTORY_VECTOR_PARAMETER(mid_smoother);
 
        bool GKO_FACTORY_PARAMETER_SCALAR(post_uses_pre, true);
 
        multigrid::mid_smooth_type GKO_FACTORY_PARAMETER_SCALAR(
            mid_case, multigrid::mid_smooth_type::standalone);
 
        size_type GKO_FACTORY_PARAMETER_SCALAR(max_levels, 10);
 
        size_type GKO_FACTORY_PARAMETER_SCALAR(min_coarse_rows, 64);
 
        std::vector<std::shared_ptr<const LinOpFactory>>
            GKO_DEFERRED_FACTORY_VECTOR_PARAMETER(coarsest_solver);
 
        std::function<size_type(const size_type, const LinOp*)>
            GKO_FACTORY_PARAMETER_SCALAR(solver_selector, nullptr);
 
        multigrid::cycle GKO_FACTORY_PARAMETER_SCALAR(cycle,
                                                      multigrid::cycle::v);
 
        size_type GKO_FACTORY_PARAMETER_SCALAR(kcycle_base, 1);
 
        double GKO_FACTORY_PARAMETER_SCALAR(kcycle_rel_tol, 0.25);
 
        std::complex<double> GKO_FACTORY_PARAMETER_SCALAR(smoother_relax, 0.9);
 
        size_type GKO_FACTORY_PARAMETER_SCALAR(smoother_iters, 1);
 
        initial_guess_mode GKO_FACTORY_PARAMETER_SCALAR(
            default_initial_guess, initial_guess_mode::zero);
    };
    class reuse_data_type;
    GKO_ENABLE_LIN_OP_FACTORY_WITH_REUSE(Multigrid, parameters, Factory,
                                         reuse_data_type);
    GKO_ENABLE_BUILD_METHOD(Factory);
 
    class reuse_data_type : public LinOpFactory::ReuseData {
        friend class Multigrid;
 
        struct level_data {
            size_type index{};
            dim<2> coarse_size{};
            std::unique_ptr<LinOpFactory::ReuseData> level;
            std::unique_ptr<LinOpFactory::ReuseData> pre_smoother;
            std::unique_ptr<LinOpFactory::ReuseData> mid_smoother;
            std::unique_ptr<LinOpFactory::ReuseData> post_smoother;
        };
 
        bool initialized_ = false;
        dim<2> size_{};
        std::vector<level_data> levels_;
        size_type coarsest_index_{};
        std::unique_ptr<LinOpFactory::ReuseData> coarsest_solver_;
    };
 
    static parameters_type parse(const config::pnode& config,
                                 const config::registry& context,
                                 const config::type_descriptor& td_for_child =
                                     config::make_type_descriptor<>());
 
protected:
    void apply_impl(const LinOp* b, LinOp* x) const override;
 
    void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta,
                    LinOp* x) const override;
 
    void apply_with_initial_guess_impl(const LinOp* b, LinOp* x,
                                       initial_guess_mode guess) const override;
 
    void apply_with_initial_guess_impl(const LinOp* alpha, const LinOp* b,
                                       const LinOp* beta, LinOp* x,
                                       initial_guess_mode guess) const override;
 
    template <typename VectorType>
    void apply_dense_impl(const VectorType* b, VectorType* x,
                          initial_guess_mode guess) const;
 
    void generate(reuse_data_type* reuse_data);
 
    void generate_smoothers(
        std::shared_ptr<const gko::multigrid::MultigridLevel> mg_level,
        size_type index,
        std::vector<std::shared_ptr<const LinOp>>& pre_smoother_list,
        std::vector<std::shared_ptr<const LinOp>>& mid_smoother_list,
        std::vector<std::shared_ptr<const LinOp>>& post_smoother_list,
        reuse_data_type::level_data* level_reuse);
 
    std::unique_ptr<LinOp> generate_coarsest_solver(
        std::shared_ptr<const gko::multigrid::MultigridLevel> last_mg_level,
        size_type level, std::shared_ptr<const LinOp> matrix,
        std::unique_ptr<LinOpFactory::ReuseData>* reuse_data,
        const size_type* recorded_index, size_type* selected_index);
 
    explicit Multigrid(std::shared_ptr<const Executor> exec);
 
    explicit Multigrid(const Factory* factory,
                       std::shared_ptr<const LinOp> system_matrix);
 
    Multigrid(const Factory* factory,
              std::shared_ptr<const LinOp> system_matrix,
              reuse_data_type& reuse_data);
 
    Multigrid(const Factory* factory,
              std::shared_ptr<const LinOp> system_matrix,
              reuse_data_type* reuse_data);
 
    static void check_reuse_consistent(const Factory* factory,
                                       const LinOp* input,
                                       const reuse_data_type& reuse_data);
 
    void validate();
 
    void verify_legal_length(bool checked, size_type len, size_type ref_len);
 
    void create_state() const;
 
private:
    std::vector<std::shared_ptr<const gko::multigrid::MultigridLevel>>
        mg_level_list_{};
    std::vector<std::shared_ptr<const LinOp>> pre_smoother_list_{};
    std::vector<std::shared_ptr<const LinOp>> mid_smoother_list_{};
    std::vector<std::shared_ptr<const LinOp>> post_smoother_list_{};
    std::shared_ptr<const LinOp> coarsest_solver_{};
    std::function<size_type(const size_type, const LinOp*)> level_selector_;
    std::function<size_type(const size_type, const LinOp*)> solver_selector_;
 
    mutable struct cache_struct {
        cache_struct() = default;
 
        ~cache_struct() = default;
 
        cache_struct(const cache_struct&) {}
 
        cache_struct(cache_struct&&) noexcept {}
 
        cache_struct& operator=(const cache_struct&) { return *this; }
 
        cache_struct& operator=(cache_struct&&) noexcept { return *this; }
 
        // unique_ptr with default destructor does not work with the incomplete
        // type.
        std::shared_ptr<multigrid::detail::MultigridState> state{};
    } cache_;
};
 
template <>
struct workspace_traits<Multigrid> {
    using Solver = Multigrid;
    // number of vectors used by this workspace
    static int num_vectors(const Solver&);
    // number of arrays used by this workspace
    static int num_arrays(const Solver&);
    // array containing the num_vectors names for the workspace vectors
    static std::vector<std::string> op_names(const Solver&);
    // array containing the num_arrays names for the workspace vectors
    static std::vector<std::string> array_names(const Solver&);
    // array containing all varying scalar vectors (independent of problem size)
    static std::vector<int> scalars(const Solver&);
    // array containing all varying vectors (dependent on problem size)
    static std::vector<int> vectors(const Solver&);
 
    // stopping status array
    constexpr static int stop = 0;
};
 
 
}  // namespace solver
}  // namespace gko
 
 
#endif  // GKO_PUBLIC_CORE_SOLVER_MULTIGRID_HPP_