diff --git a/benchmark/utils/preconditioners.hpp b/benchmark/utils/preconditioners.hpp index a2166cf0270..03bbc808dca 100644 --- a/benchmark/utils/preconditioners.hpp +++ b/benchmark/utils/preconditioners.hpp @@ -78,6 +78,18 @@ DEFINE_double(mg_tolerance, false, "The tolerance for the coarse solver"); DEFINE_uint32(mg_max_iters, false, "The max number of iterations for the coarse solver"); +DEFINE_string(mg_scale_correction, "none", + "OpenFOAM-style Rayleigh scale correction mode for the Multigrid " + "preconditioner: 'none' (off), 'post' (coarse-correction scaling " + "only), or 'both' (pre-smooth + coarse-correction scaling)"); + +DEFINE_uint32(mg_smoother_iters, 2, + "Number of block-Jacobi Richardson sweeps used for the Multigrid " + "pre/post smoother (the coarsest solver uses 4x this)"); + +DEFINE_double(mg_smoother_relax, 0.8, + "Relaxation factor of the Multigrid block-Jacobi smoother"); + // parses the Jacobi storage optimization command line argument gko::precision_reduction parse_storage_optimization(const std::string& flag) @@ -315,6 +327,7 @@ const std::map( {"mg", [](std::shared_ptr exec) { using ir = gko::solver::Ir; + using jacobi = gko::preconditioner::Jacobi; auto iter_stop = gko::share(gko::stop::Iteration::build() .with_max_iters(FLAGS_mg_max_iters) .on(exec)); @@ -323,12 +336,45 @@ const std::map( .with_baseline(gko::stop::mode::absolute) .with_reduction_factor(FLAGS_mg_tolerance) .on(exec)); + // damped block-Jacobi Richardson smoother (and a heavier variant + // as the coarsest solver) -- required for the multigrid V-cycle to + // be an effective preconditioner, and hence for scale correction + // to have a measurable effect. + auto smoother = gko::share( + ir::build() + .with_solver(jacobi::build().with_max_block_size(1u)) + .with_relaxation_factor( + static_cast(FLAGS_mg_smoother_relax)) + .with_criteria(gko::stop::Iteration::build().with_max_iters( + FLAGS_mg_smoother_iters)) + .on(exec)); + auto coarsest = gko::share( + ir::build() + .with_solver(jacobi::build().with_max_block_size(1u)) + .with_relaxation_factor( + static_cast(FLAGS_mg_smoother_relax)) + .with_criteria(gko::stop::Iteration::build().with_max_iters( + 4u * FLAGS_mg_smoother_iters)) + .on(exec)); + const auto& sc_mode = FLAGS_mg_scale_correction; + if (sc_mode != "none" && sc_mode != "post" && sc_mode != "both") { + throw std::runtime_error( + "Unknown -mg_scale_correction mode '" + sc_mode + + "', expected 'none', 'post' or 'both'"); + } + const bool scale_correction = sc_mode != "none"; + const bool scale_correction_pre = sc_mode == "both"; return gko::solver::Multigrid::build() .with_mg_level( gko::multigrid::Pgm::build() .with_deterministic(FLAGS_pgm_deterministic)) + .with_pre_smoother(smoother) + .with_post_smoother(smoother) + .with_coarsest_solver(coarsest) .with_criteria(iter_stop, tol_stop) .with_max_levels(FLAGS_mg_max_num_levels) + .with_scale_correction(scale_correction) + .with_scale_correction_pre_smooth(scale_correction_pre) .on(exec); }} #if GINKGO_BUILD_MPI diff --git a/core/solver/multigrid.cpp b/core/solver/multigrid.cpp index 20649831f1f..435453bb7be 100644 --- a/core/solver/multigrid.cpp +++ b/core/solver/multigrid.cpp @@ -5,6 +5,7 @@ #include "ginkgo/core/solver/multigrid.hpp" #include +#include #include #include @@ -291,6 +292,13 @@ class MultigridState { std::vector> one_list; std::vector> next_one_list; std::vector> neg_one_list; + // scale correction workspace (only allocated when scale_correction=true) + std::vector> acf_list; // A*δ scratch, nrows×nrhs + std::vector> + delta_pre_list; // A*δ_c / smoother scratch + std::vector> alpha_list; // Rayleigh scalar, 1×nrhs + std::vector> + denom_list; // denominator scalar, 1×nrhs const LinOp* system_matrix; const Multigrid* multigrid; size_type nrhs; @@ -314,6 +322,12 @@ void MultigridState::generate(const LinOp* system_matrix_in, clear_and_reserve(one_list, list_size); clear_and_reserve(next_one_list, list_size); clear_and_reserve(neg_one_list, list_size); + if (multigrid_in->get_parameters().scale_correction) { + clear_and_reserve(acf_list, list_size); + clear_and_reserve(delta_pre_list, list_size); + clear_and_reserve(alpha_list, list_size); + clear_and_reserve(denom_list, list_size); + } // Allocate memory first such that reusing allocation in each iter. for (int i = 0; i < mg_level_list.size(); i++) { auto next_nrows = mg_level_list.at(i)->get_coarse_op()->get_size()[0]; @@ -401,6 +415,13 @@ void MultigridState::allocate_memory(int level, multigrid::cycle cycle, } one_list.emplace_back(initialize({one()}, exec)); neg_one_list.emplace_back(initialize({-one()}, exec)); + if (multigrid->get_parameters().scale_correction) { + acf_list.emplace_back(vec::create(exec, dim<2>{current_nrows, nrhs})); + delta_pre_list.emplace_back( + vec::create(exec, dim<2>{current_nrows, nrhs})); + alpha_list.emplace_back(vec::create(exec, dim<2>{1, nrhs})); + denom_list.emplace_back(vec::create(exec, dim<2>{1, nrhs})); + } } @@ -449,6 +470,16 @@ void MultigridState::allocate_memory( one_list.emplace_back(initialize({one()}, exec)); neg_one_list.emplace_back( initialize({-one()}, exec)); + if (multigrid->get_parameters().scale_correction) { + acf_list.emplace_back(vec::create(exec, current_comm, + dim<2>{current_nrows, nrhs}, + dim<2>{current_local_nrows, nrhs})); + delta_pre_list.emplace_back( + vec::create(exec, current_comm, dim<2>{current_nrows, nrhs}, + dim<2>{current_local_nrows, nrhs})); + alpha_list.emplace_back(dense_vec::create(exec, dim<2>{1, nrhs})); + denom_list.emplace_back(dense_vec::create(exec, dim<2>{1, nrhs})); + } } @@ -503,22 +534,47 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, auto r = r_list.at(level); auto g = g_list.at(level); auto e = e_list.at(level); - // get mg_level auto mg_level = multigrid->get_mg_level_list().at(level); - // get the pre_smoother auto pre_smoother = multigrid->get_pre_smoother_list().at(level); - // get the mid_smoother std::shared_ptr mid_smoother{nullptr}; auto mid_case = multigrid->get_parameters().mid_case; if (mid_case == multigrid::mid_smooth_type::standalone) { mid_smoother = multigrid->get_mid_smoother_list().at(level); } - // get the post_smoother auto post_smoother = multigrid->get_post_smoother_list().at(level); auto one = one_list.at(level).get(); auto next_one = next_one_list.at(level).get(); auto neg_one = neg_one_list.at(level).get(); - // origin or next or first + + // scale correction applies at all levels except immediately above coarsest + // (at that level the coarse solver provides a near-exact result, sf ≈ 1) + bool do_scale = + multigrid->get_parameters().scale_correction && level < total_level - 1; + // the pre-smooth (downward) scaling can be disabled independently, leaving + // only the post-smooth (coarse-correction) scaling ("post-only" mode) + bool do_pre_scale = + do_scale && multigrid->get_parameters().scale_correction_pre_smooth; + + // [NeoN patch] device-side guarded Rayleigh reciprocal: sf = num / (denom + eps), + // computed entirely on-device to avoid the per-correction-point copy_val_to_host(denom) + // D2H sync (that guard fires a synchronizing device->host copy at every level every + // V-cycle -- ~90/solve on the occDrivAer pressure solve, a leading cost when the solve is + // synchronization-bound). Safe because num = delta.b and denom = delta.A.delta are BOTH + // exactly zero iff delta = 0, so at delta=0 sf = 0/eps = 0 (the correction reduces to + // smoother(b) -- benign, and does not occur after a pre-smoother sweep on nonzero b); + // eps = smallest positive normal is negligible vs any nonzero denom. Replaces the old + // `if (copy_val_to_host(denom) != 0) { inv_scale(...); }`. + auto safe_inv_scale = [&](matrix::Dense* alpha_d, + matrix::Dense* denom_d) { + using real_type = gko::remove_complex; + auto exec_l = multigrid->get_executor(); + auto eps = matrix::Dense::create( + exec_l, dim<2>{1, alpha_d->get_size()[1]}); + eps->fill(value_type{std::numeric_limits::min()}); + denom_d->add_scaled(one, eps.get()); // denom += eps (device, no sync) + alpha_d->inv_scale(denom_d); // sf = num / (denom + eps) + }; + bool use_pre = has_property(mode, cycle_mode::first_of_cycle) || mid_case == multigrid::mid_smooth_type::both || mid_case == multigrid::mid_smooth_type::pre_smoother; @@ -540,15 +596,49 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, pre_smoother->apply(b, x); } } + + // Pre-smooth scale correction (OpenFOAM GAMGSolverSolve.C downward pass): + // Rayleigh-scale δ_pre = x, deflating r before restriction. + // Aδ = A * δ_pre + // sf = (δ_pre · b) / (δ_pre · Aδ) + // δ_pre = sf * δ_pre + smoother(b − sf * Aδ) [reuses acf, r scratch] + if (do_pre_scale && use_pre && pre_smoother) { + auto acf = acf_list.at(level); + auto dp = delta_pre_list.at(level); + auto alpha_dense = as>(alpha_list.at(level)); + auto denom_dense = as>(denom_list.at(level)); + + matrix->apply(x, acf); // acf = A * δ_pre + as(x)->compute_dot(b, alpha_list.at(level)); + as(x)->compute_dot(acf, denom_list.at(level)); + // [NeoN patch] device-side guarded reciprocal (no copy_val_to_host D2H sync). + safe_inv_scale(alpha_dense.get(), denom_dense.get()); // sf = (δ·b)/(δ·Aδ) + { + // r temporarily holds r_scaled = b − sf * Aδ + as(acf)->scale(alpha_dense); // acf = sf * Aδ + as(r)->copy_from(as(b)); + as(r)->add_scaled(neg_one, acf); // r = b − sf*Aδ + + // dp = smoother(r_scaled) starting from zero + as(dp)->fill(zero()); + pre_smoother->apply(r, dp); + + // x = sf * δ_pre + smoother(b − sf*Aδ) + as(x)->scale(alpha_dense); + as(x)->add_scaled(one, dp); + } + // r is overwritten with the actual (deflated) residual below + } + // The common smoother is wrapped by IR and IR already split the iter and // residual check. Thus, when the IR only contains iter limit, there's no - // additional residual computation + // additional residual computation. // TODO: if already computes the residual outside, the first level may not // need this residual computation when no presmoother in the first level. - as(r)->copy_from(as(b)); // n * b - matrix->apply(neg_one, x, one, r); + as(r)->copy_from(as(b)); + matrix->apply(neg_one, x, one, r); // r = b − A*x (deflated if scaled) - // first cycle + // restrict mg_level->get_restrict_op()->apply(r, g); // next level if (level + 1 == total_level) { @@ -568,9 +658,8 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, next_mode); if (level < multigrid->get_mg_level_list().size() - 1) { // additional work for non-v_cycle - // next level if (cycle == multigrid::cycle::f) { - // f_cycle call v_cycle in the second cycle + // f_cycle calls v_cycle in the second cycle this->run_mg_cycle(multigrid::cycle::v, level + 1, next_level_matrix, g.get(), e.get(), cycle_mode::end_of_cycle); @@ -579,20 +668,61 @@ void MultigridState::run_cycle(multigrid::cycle cycle, size_type level, e.get(), cycle_mode::end_of_cycle); } } - // prolong - mg_level->get_prolong_op()->apply(next_one, e, next_one, x); - // end or origin previous + // Post-smooth scale correction (OpenFOAM GAMGSolverSolve.C upward pass): + // Prolong coarse correction δ_c into acf, Rayleigh-scale it w.r.t. the + // deflated residual r, then merge with δ_pre (= current x). + // δ_c = prolong(e) + // Aδ = A * δ_c [stored in delta_pre scratch] + // sf = (δ_c · r) / (δ_c · Aδ) + // δ_c = sf * δ_c + smoother(r − sf * Aδ) + // x += δ_c [x = δ_pre + scale-corrected δ_c] + if (do_scale) { + auto acf = acf_list.at(level); + auto dp = delta_pre_list.at(level); + auto alpha_dense = as>(alpha_list.at(level)); + auto denom_dense = as>(denom_list.at(level)); + + // prolong e into acf (δ_c = prolong(e)) + as(acf)->fill(zero()); + mg_level->get_prolong_op()->apply(next_one, e, next_one, acf); + + matrix->apply(acf, dp); // dp = A * δ_c + as(acf)->compute_dot(r, alpha_list.at(level)); + as(acf)->compute_dot(dp, denom_list.at(level)); + // [NeoN patch] device-side guarded reciprocal (no copy_val_to_host D2H sync). + safe_inv_scale(alpha_dense.get(), + denom_dense.get()); // sf = (δ_c·r)/(δ_c·Aδ_c) + { + // r temporarily holds r_scaled = r − sf * Aδ_c + as(dp)->scale(alpha_dense); // dp = sf * Aδ_c + as(r)->add_scaled(neg_one, dp); // r = r − sf*Aδ_c + + // dp = smoother(r_scaled) starting from zero + as(dp)->fill(zero()); + if (pre_smoother) { + pre_smoother->apply(r, dp); + } + + // acf = sf * δ_c + smoother(r − sf*Aδ_c) + as(acf)->scale(alpha_dense); + as(acf)->add_scaled(one, dp); + } + // x = δ_pre + scale-corrected δ_c + as(x)->add_scaled(one, acf); + } else { + // standard prolongation: x += prolong(e) + mg_level->get_prolong_op()->apply(next_one, e, next_one, x); + } + bool use_post = has_property(mode, cycle_mode::end_of_cycle) || mid_case == multigrid::mid_smooth_type::both || mid_case == multigrid::mid_smooth_type::post_smoother; - // post-smooth if (use_post && post_smoother) { post_smoother->apply(b, x); } - // put the mid smoother into the end of previous cycle - // only W/F cycle + // put the mid smoother into the end of previous cycle (W/F cycle only) bool use_mid = (cycle == multigrid::cycle::w || cycle == multigrid::cycle::f) && !has_property(mode, cycle_mode::end_of_cycle) && @@ -695,6 +825,12 @@ typename Multigrid::parameters_type Multigrid::parse( params.with_default_initial_guess( config::get_value(obj)); } + if (auto& obj = config_check.get("scale_correction")) { + params.with_scale_correction(config::get_value(obj)); + } + if (auto& obj = config_check.get("scale_correction_pre_smooth")) { + params.with_scale_correction_pre_smooth(config::get_value(obj)); + } return params; } diff --git a/core/test/solver/multigrid.cpp b/core/test/solver/multigrid.cpp index 6eb279584bd..88a339bdd16 100644 --- a/core/test/solver/multigrid.cpp +++ b/core/test/solver/multigrid.cpp @@ -837,4 +837,70 @@ TYPED_TEST(Multigrid, DeferredFactoryParameter) } +TYPED_TEST(Multigrid, ScaleCorrectionIsDisabledByDefault) +{ + using Solver = typename TestFixture::Solver; + + auto factory = Solver::build().on(this->exec); + + ASSERT_FALSE(factory->get_parameters().scale_correction); +} + + +TYPED_TEST(Multigrid, ScaleCorrectionCanBeEnabled) +{ + using Solver = typename TestFixture::Solver; + + auto factory = Solver::build().with_scale_correction(true).on(this->exec); + + ASSERT_TRUE(factory->get_parameters().scale_correction); +} + + +TYPED_TEST(Multigrid, ScaleCorrectionPreSmoothDefaultsTrue) +{ + using Solver = typename TestFixture::Solver; + + auto factory = Solver::build().on(this->exec); + + ASSERT_TRUE(factory->get_parameters().scale_correction_pre_smooth); +} + + +TYPED_TEST(Multigrid, ScaleCorrectionPreSmoothCanBeDisabled) +{ + using Solver = typename TestFixture::Solver; + + auto factory = + Solver::build().with_scale_correction_pre_smooth(false).on(this->exec); + + ASSERT_FALSE(factory->get_parameters().scale_correction_pre_smooth); +} + + +TYPED_TEST(Multigrid, ScaleCorrectionIsPropagatedToGeneratedSolver) +{ + using Solver = typename TestFixture::Solver; + using DummyRPFactory = typename TestFixture::DummyRPFactory; + using DummyFactory = typename TestFixture::DummyFactory; + + auto solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(1u)) + .with_max_levels(2u) + .with_min_coarse_rows(2u) + .with_mg_level(this->rp_factory) + .with_pre_smoother(this->lo_factory) + .with_post_smoother(this->lo_factory) + .with_coarsest_solver(this->lo_factory) + .with_scale_correction(true) + .on(this->exec) + ->generate(this->mtx); + + ASSERT_TRUE(static_cast(solver.get()) + ->get_parameters() + .scale_correction); +} + + } // namespace diff --git a/include/ginkgo/core/solver/multigrid.hpp b/include/ginkgo/core/solver/multigrid.hpp index ede8990d000..2d9cb2a05e5 100644 --- a/include/ginkgo/core/solver/multigrid.hpp +++ b/include/ginkgo/core/solver/multigrid.hpp @@ -375,6 +375,42 @@ class Multigrid : public LinOp, */ initial_guess_mode GKO_FACTORY_PARAMETER_SCALAR( default_initial_guess, initial_guess_mode::zero); + + /** + * Per-level Rayleigh-quotient scale correction (mirrors OpenFOAM's + * scaleCorrection_ flag). When enabled, at each level l up to but not + * including the level immediately above the coarsest: + * + * Pre-smooth: the pre-smoother correction δ_pre is Rayleigh-scaled + * before restriction, deflating the component it already captures. + * Aδ = A_l * δ_pre + * sf = (δ_pre · r_l) / (δ_pre · Aδ) + * δ_pre = sf * δ_pre + smoother_l(r_l − sf * Aδ) + * r_l -= A_l * δ_pre [restrict deflated r down] + * + * Post-smooth: the prolonged coarse correction δ_c is Rayleigh-scaled + * before being merged with δ_pre and post-smoothed. + * Aδ = A_l * δ_c + * sf = (δ_c · r_l) / (δ_c · Aδ) [r_l = deflated residual] + * δ_c = sf * δ_c + smoother_l(r_l − sf * Aδ) + * x = δ_pre + δ_c + * + * Mirrors OpenFOAM GAMGSolverSolve.C::Vcycle() + + * GAMGSolverScale.C::scale(), replacing the hardcoded Jacobi (÷D) + * step with the configured pre_smoother. + */ + bool GKO_FACTORY_PARAMETER_SCALAR(scale_correction, false); + + /** + * When scale_correction is enabled, also apply the *pre-smooth* + * (downward pass) Rayleigh scaling described above. When false, only + * the post-smooth (coarse-correction) scaling is applied -- this is the + * "post-only" mode, which in practice captures nearly all of the + * benefit at half the extra work. Has no effect when scale_correction + * is false. Enabled by default so that scale_correction alone + * reproduces the full OpenFOAM-style pre+post behavior. + */ + bool GKO_FACTORY_PARAMETER_SCALAR(scale_correction_pre_smooth, true); }; GKO_ENABLE_LIN_OP_FACTORY(Multigrid, parameters, Factory); GKO_ENABLE_BUILD_METHOD(Factory); diff --git a/reference/test/solver/CMakeLists.txt b/reference/test/solver/CMakeLists.txt index b70a3299a44..df948a765a4 100644 --- a/reference/test/solver/CMakeLists.txt +++ b/reference/test/solver/CMakeLists.txt @@ -17,5 +17,6 @@ ginkgo_create_test(lower_trs) ginkgo_create_test(lower_trs_kernels) ginkgo_create_test(minres_kernels) ginkgo_create_test(multigrid_kernels) +ginkgo_create_test(multigrid_scale_correction) ginkgo_create_test(upper_trs) ginkgo_create_test(upper_trs_kernels) diff --git a/reference/test/solver/multigrid_scale_correction.cpp b/reference/test/solver/multigrid_scale_correction.cpp new file mode 100644 index 00000000000..7bf0ceaaefc --- /dev/null +++ b/reference/test/solver/multigrid_scale_correction.cpp @@ -0,0 +1,244 @@ +// SPDX-FileCopyrightText: 2017 - 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" + + +namespace { + + +// Exercises the multigrid `scale_correction` feature (OpenFOAM-style Rayleigh +// scaling of the pre-smooth and coarse corrections, implemented in +// MultigridState::run_cycle). These are behavioral tests on the reference +// executor: a real Pgm V-cycle multigrid is run as a stand-alone solver on an +// SPD system, so the whole scale-correction path (compute_dot -> device-side +// guarded reciprocal `safe_inv_scale` -> scaled correction) is executed. +class MultigridScaleCorrection : public ::testing::Test { +protected: + using value_type = double; + using index_type = int; + using Csr = gko::matrix::Csr; + using Vec = gko::matrix::Dense; + using Pgm = gko::multigrid::Pgm; + using Ir = gko::solver::Ir; + using Cg = gko::solver::Cg; + using Jacobi = gko::preconditioner::Jacobi; + using Mg = gko::solver::Multigrid; + + MultigridScaleCorrection() + : exec(gko::ReferenceExecutor::create()), + n(300), + mtx(gko::share(Csr::create(exec))), + x_exact(Vec::create(exec, gko::dim<2>(n, 1))), + b(Vec::create(exec, gko::dim<2>(n, 1))) + { + // 1D Laplacian tridiag(-1, 2, -1): symmetric positive definite, and + // > min_coarse_rows so Pgm produces several levels (scale correction + // is only active on levels above the coarsest). + gko::matrix_data data(gko::dim<2>(n, n)); + for (index_type i = 0; i < n; i++) { + if (i > 0) { + data.nonzeros.emplace_back(i, i - 1, value_type{-1.0}); + } + data.nonzeros.emplace_back(i, i, value_type{2.0}); + if (i < n - 1) { + data.nonzeros.emplace_back(i, i + 1, value_type{-1.0}); + } + } + mtx->read(data); + + x_exact->fill(value_type{1.0}); + mtx->apply(x_exact, b); // b = A * x_exact + } + + // Smoother / coarsest-solver factories shared by every multigrid we build. + std::shared_ptr smoother_factory() + { + return Ir::build() + .with_solver(Jacobi::build().with_max_block_size(1u)) + .with_relaxation_factor(value_type{2.0 / 3.0}) + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .on(exec); + } + + std::shared_ptr coarsest_factory() + { + return Cg::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(50u), + gko::stop::ResidualNorm::build() + .with_reduction_factor(value_type{1e-4})) + .on(exec); + } + + // The three scale-correction configurations exposed by the feature. + enum class sc_mode { none, post_only, pre_and_post }; + + // A single-V-cycle multigrid factory used as a preconditioner (the way the + // OpenFOAM-style scale-corrected multigrid is deployed in production: it is + // a *nonlinear* correction and is driven by an outer Krylov solver rather + // than iterated stand-alone). + std::shared_ptr mg_precond_factory(sc_mode mode) + { + return Mg::build() + .with_mg_level(Pgm::build().with_deterministic(true)) + .with_pre_smoother(smoother_factory()) + .with_post_smoother(smoother_factory()) + .with_coarsest_solver(coarsest_factory()) + .with_min_coarse_rows(16u) + .with_scale_correction(mode != sc_mode::none) + .with_scale_correction_pre_smooth(mode == sc_mode::pre_and_post) + .with_criteria(gko::stop::Iteration::build().with_max_iters(1u)) + .on(exec); + } + + // CG preconditioned by a (possibly scale-corrected) multigrid V-cycle. + std::unique_ptr build_cg_mg(sc_mode mode) + { + return Cg::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(300u), + gko::stop::ResidualNorm::build() + .with_baseline(gko::stop::mode::rhs_norm) + .with_reduction_factor(value_type{1e-11})) + .with_preconditioner(mg_precond_factory(mode)) + .on(exec) + ->generate(mtx); + } + + // A stand-alone V-cycle multigrid solver running a fixed number of cycles + // (Iteration criterion only) -- used by the zero-rhs guard test, which + // needs the scale-correction path to execute even though the residual is + // exactly zero. + std::unique_ptr build_standalone_mg(bool scale_correction, + unsigned num_cycles) + { + return Mg::build() + .with_mg_level(Pgm::build().with_deterministic(true)) + .with_pre_smoother(smoother_factory()) + .with_post_smoother(smoother_factory()) + .with_coarsest_solver(coarsest_factory()) + .with_min_coarse_rows(16u) + .with_scale_correction(scale_correction) + .with_default_initial_guess( + gko::solver::initial_guess_mode::provided) + .with_criteria( + gko::stop::Iteration::build().with_max_iters(num_cycles)) + .on(exec) + ->generate(mtx); + } + + // relative residual ||b - A x|| / ||b|| + value_type relative_residual(const Vec* x) + { + auto res = gko::clone(b); + auto one = gko::initialize({value_type{1.0}}, exec); + auto neg_one = gko::initialize({value_type{-1.0}}, exec); + mtx->apply(neg_one, x, one, res); // res = b - A x + auto rnorm = Vec::create(exec, gko::dim<2>(1, 1)); + auto bnorm = Vec::create(exec, gko::dim<2>(1, 1)); + res->compute_norm2(rnorm); + b->compute_norm2(bnorm); + return rnorm->at(0, 0) / bnorm->at(0, 0); + } + + std::shared_ptr exec; + index_type n; + std::shared_ptr mtx; + std::unique_ptr x_exact; + std::unique_ptr b; +}; + + +TEST_F(MultigridScaleCorrection, ScaleCorrectedVCyclePreconditionsCgToSolution) +{ + // A CG solve preconditioned by the scale-corrected multigrid V-cycle must + // converge to the true solution of the SPD system. This drives every part + // of the scale-correction code path (pre-smooth + coarse Rayleigh scaling + // via the device-side safe_inv_scale) once per outer iteration. + auto solver = build_cg_mg(sc_mode::pre_and_post); + auto x = Vec::create(exec, gko::dim<2>(n, 1)); + x->fill(value_type{0.0}); + + solver->apply(b, x); + + ASSERT_LT(relative_residual(x.get()), value_type{1e-9}); + GKO_ASSERT_MTX_NEAR(x, x_exact, value_type{1e-5}); +} + + +TEST_F(MultigridScaleCorrection, PostOnlyScaleCorrectionPreconditionsCgToSolution) +{ + // The "post-only" mode (coarse-correction scaling, pre-smooth scaling + // disabled) must also yield a correct CG solve. Exercises the do_pre_scale + // == false path while do_scale == true. + auto solver = build_cg_mg(sc_mode::post_only); + auto x = Vec::create(exec, gko::dim<2>(n, 1)); + x->fill(value_type{0.0}); + + solver->apply(b, x); + + ASSERT_LT(relative_residual(x.get()), value_type{1e-9}); + GKO_ASSERT_MTX_NEAR(x, x_exact, value_type{1e-5}); +} + + +TEST_F(MultigridScaleCorrection, AllModesReachTheSameSolution) +{ + // Scale correction (in either mode) changes the multigrid convergence path + // but not the fixed point: all three configurations drive CG to the same + // solution of the SPD system. + auto x_none = Vec::create(exec, gko::dim<2>(n, 1)); + auto x_post = Vec::create(exec, gko::dim<2>(n, 1)); + auto x_both = Vec::create(exec, gko::dim<2>(n, 1)); + x_none->fill(value_type{0.0}); + x_post->fill(value_type{0.0}); + x_both->fill(value_type{0.0}); + + build_cg_mg(sc_mode::none)->apply(b, x_none); + build_cg_mg(sc_mode::post_only)->apply(b, x_post); + build_cg_mg(sc_mode::pre_and_post)->apply(b, x_both); + + GKO_ASSERT_MTX_NEAR(x_post, x_none, value_type{1e-5}); + GKO_ASSERT_MTX_NEAR(x_both, x_none, value_type{1e-5}); +} + + +TEST_F(MultigridScaleCorrection, ZeroRhsStaysZeroWithoutNan) +{ + // Regression guard for the device-side reciprocal `safe_inv_scale`: with a + // zero rhs and zero initial guess every correction delta is exactly zero, + // so the Rayleigh denominator delta.A.delta is exactly zero. The guard must + // yield a scale factor of 0 (num/(0 + eps) = 0), NOT a 0/0 NaN. Iteration- + // only criteria force the V-cycles to actually run. + auto solver = build_standalone_mg(true, 3u); + auto x = Vec::create(exec, gko::dim<2>(n, 1)); + x->fill(value_type{0.0}); + auto zero = Vec::create(exec, gko::dim<2>(n, 1)); + zero->fill(value_type{0.0}); + + solver->apply(zero, x); + + for (index_type i = 0; i < n; i++) { + ASSERT_TRUE(gko::is_finite(x->at(i, 0))); + } + GKO_ASSERT_MTX_NEAR(x, zero, value_type{0.0}); +} + + +} // namespace