diff --git a/common/unified/CMakeLists.txt b/common/unified/CMakeLists.txt index 847be9b333a..ccd9f7c69d8 100644 --- a/common/unified/CMakeLists.txt +++ b/common/unified/CMakeLists.txt @@ -25,6 +25,7 @@ set(UNIFIED_SOURCES solver/bicgstab_kernels.cpp solver/cg_kernels.cpp solver/cgs_kernels.cpp + solver/chebyshev_kernels.cpp solver/common_gmres_kernels.cpp solver/fcg_kernels.cpp solver/gcr_kernels.cpp diff --git a/common/unified/solver/chebyshev_kernels.cpp b/common/unified/solver/chebyshev_kernels.cpp new file mode 100644 index 00000000000..a2461a1e930 --- /dev/null +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -0,0 +1,114 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include + +#include + +#include "common/unified/base/kernel_launch.hpp" + + +namespace gko { +namespace kernels { +namespace GKO_DEVICE_NAMESPACE { +namespace chebyshev { + + +template +using coeff_type = gko::kernels::chebyshev::coeff_type; + + +#if GINKGO_DPCPP_SINGLE_MODE + + +// we only change type in device code to keep the interface is the same as the +// other backend. +template +using if_single_only_type = + std::conditional_t, float, + std::complex>; + + +#else + + +template +struct type_identity { + using type = T; +}; + + +template +using if_single_only_type = typename type_identity::type; + + +#endif + + +template +void init_update(std::shared_ptr exec, + const coeff_type alpha, + const matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + using actual_coeff_type = if_single_only_type>; + using type = device_type; + + auto alpha_val = static_cast(alpha); + + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto col, auto alpha, auto inner_sol, + auto update_sol, auto output) { + const auto inner_val = static_cast(inner_sol(row, col)); + update_sol(row, col) = + static_cast>(inner_val); + output(row, col) = static_cast>( + static_cast(output(row, col)) + + static_cast(alpha) * inner_val); + }, + output->get_size(), alpha_val, inner_sol, update_sol, output); +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); + + +template +void update(std::shared_ptr exec, const ScalarType alpha, + const ScalarType beta, matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + using actual_coeff_type = if_single_only_type>; + using type = device_type; + + auto alpha_val = static_cast(alpha); + auto beta_val = static_cast(beta); + + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto col, auto alpha, auto beta, auto inner_sol, + auto update_sol, auto output) { + const auto val = static_cast(inner_sol(row, col)) + + static_cast(beta) * + static_cast(update_sol(row, col)); + inner_sol(row, col) = static_cast>(val); + update_sol(row, col) = static_cast>(val); + output(row, col) = static_cast>( + static_cast(output(row, col)) + + static_cast(alpha) * val); + }, + output->get_size(), alpha_val, beta_val, inner_sol, update_sol, output); +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); + + +} // namespace chebyshev +} // namespace GKO_DEVICE_NAMESPACE +} // namespace kernels +} // namespace gko diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index 54ce8fb59e6..e4f6fce6dc7 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -92,8 +92,8 @@ target_sources(${ginkgo_core} matrix/scaled_permutation.cpp matrix/sellp.cpp matrix/sparsity_csr.cpp - multigrid/pgm.cpp multigrid/fixed_coarsening.cpp + multigrid/pgm.cpp preconditioner/batch_jacobi.cpp preconditioner/gauss_seidel.cpp preconditioner/sor.cpp @@ -112,6 +112,7 @@ target_sources(${ginkgo_core} solver/cb_gmres.cpp solver/cg.cpp solver/cgs.cpp + solver/chebyshev.cpp solver/direct.cpp solver/fcg.cpp solver/gcr.cpp diff --git a/core/config/config_helper.hpp b/core/config/config_helper.hpp index b1fa7bd69b5..7da9391e61b 100644 --- a/core/config/config_helper.hpp +++ b/core/config/config_helper.hpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -24,10 +24,9 @@ namespace gko { namespace config { -#define GKO_INVALID_CONFIG_VALUE(_entry, _value) \ - GKO_INVALID_STATE(std::string("The value >" + _value + \ - "< is invalid for the entry >" + _entry + \ - "<")) +#define GKO_INVALID_CONFIG_VALUE(_entry, _value) \ + GKO_INVALID_STATE(std::string("The value >") + _value + \ + "< is invalid for the entry >" + _entry + "<") #define GKO_MISSING_CONFIG_ENTRY(_entry) \ @@ -52,6 +51,7 @@ enum class LinOpFactoryType : int { Direct, LowerTrs, UpperTrs, + Chebyshev, Factorization_Ic, Factorization_Ilu, Cholesky, diff --git a/core/config/registry.cpp b/core/config/registry.cpp index 19da3ed2559..d85f9bf241b 100644 --- a/core/config/registry.cpp +++ b/core/config/registry.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -32,6 +32,7 @@ configuration_map generate_config_map() {"solver::Direct", parse}, {"solver::LowerTrs", parse}, {"solver::UpperTrs", parse}, + {"solver::Chebyshev", parse}, {"factorization::Ic", parse}, {"factorization::Ilu", parse}, {"factorization::Cholesky", parse}, diff --git a/core/config/solver_config.cpp b/core/config/solver_config.cpp index d13c20f901e..d6c88c868b5 100644 --- a/core/config/solver_config.cpp +++ b/core/config/solver_config.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -43,6 +44,7 @@ GKO_PARSE_VALUE_TYPE_BASE(CbGmres, gko::solver::CbGmres); GKO_PARSE_VALUE_AND_INDEX_TYPE(Direct, gko::experimental::solver::Direct); GKO_PARSE_VALUE_AND_INDEX_TYPE(LowerTrs, gko::solver::LowerTrs); GKO_PARSE_VALUE_AND_INDEX_TYPE(UpperTrs, gko::solver::UpperTrs); +GKO_PARSE_VALUE_TYPE(Chebyshev, gko::solver::Chebyshev); template <> diff --git a/core/device_hooks/common_kernels.inc.cpp b/core/device_hooks/common_kernels.inc.cpp index 5ed2daa540c..8faf4fd12ba 100644 --- a/core/device_hooks/common_kernels.inc.cpp +++ b/core/device_hooks/common_kernels.inc.cpp @@ -63,6 +63,7 @@ #include "core/solver/cb_gmres_kernels.hpp" #include "core/solver/cg_kernels.hpp" #include "core/solver/cgs_kernels.hpp" +#include "core/solver/chebyshev_kernels.hpp" #include "core/solver/common_gmres_kernels.hpp" #include "core/solver/fcg_kernels.hpp" #include "core/solver/gcr_kernels.hpp" @@ -674,6 +675,16 @@ GKO_STUB_CB_GMRES_CONST(GKO_DECLARE_CB_GMRES_SOLVE_KRYLOV_KERNEL); } // namespace cb_gmres +namespace chebyshev { + + +GKO_STUB_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); +GKO_STUB_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); + + +} // namespace chebyshev + + namespace ir { diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp new file mode 100644 index 00000000000..409204d60b5 --- /dev/null +++ b/core/solver/chebyshev.cpp @@ -0,0 +1,341 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "ginkgo/core/solver/chebyshev.hpp" + +#include +#include +#include + +#include "core/config/solver_config.hpp" +#include "core/distributed/helpers.hpp" +#include "core/solver/chebyshev_kernels.hpp" +#include "core/solver/ir_kernels.hpp" +#include "core/solver/solver_base.hpp" +#include "core/solver/solver_boilerplate.hpp" +#include "core/solver/update_residual.hpp" + + +namespace gko { +namespace solver { +namespace chebyshev { +namespace { + + +GKO_REGISTER_OPERATION(initialize, ir::initialize); +GKO_REGISTER_OPERATION(init_update, chebyshev::init_update); +GKO_REGISTER_OPERATION(update, chebyshev::update); + + +} // anonymous namespace +} // namespace chebyshev + +template +typename Chebyshev::parameters_type Chebyshev::parse( + const config::pnode& config, const config::registry& context, + const config::type_descriptor& td_for_child) +{ + auto params = solver::Chebyshev::build(); + common_solver_parse(params, config, context, td_for_child); + if (auto& obj = config.get("foci")) { + auto arr = obj.get_array(); + if (arr.size() != 2) { + GKO_INVALID_CONFIG_VALUE("foci", "must contain two elements"); + } + params.with_foci(gko::config::get_value(arr.at(0)), + gko::config::get_value(arr.at(1))); + } + if (auto& obj = config.get("default_initial_guess")) { + params.with_default_initial_guess( + gko::config::get_value(obj)); + } + return params; +} + + +template +Chebyshev::Chebyshev(const Factory* factory, + std::shared_ptr system_matrix) + : EnableLinOp(factory->get_executor(), + gko::transpose(system_matrix->get_size())), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{ + this->set_default_initial_guess(parameters_.default_initial_guess); + center_ = (std::get<0>(parameters_.foci) + std::get<1>(parameters_.foci)) / + ValueType{2}; + foci_direction_ = + (std::get<1>(parameters_.foci) - std::get<0>(parameters_.foci)) / + ValueType{2}; +} + + +template +Chebyshev& Chebyshev::operator=(const Chebyshev& other) +{ + if (&other != this) { + EnableLinOp::operator=(other); + EnablePreconditionedIterativeSolver< + ValueType, Chebyshev>::operator=(other); + this->parameters_ = other.parameters_; + } + return *this; +} + + +template +Chebyshev& Chebyshev::operator=(Chebyshev&& other) +{ + if (&other != this) { + EnableLinOp::operator=(std::move(other)); + EnablePreconditionedIterativeSolver< + ValueType, Chebyshev>::operator=(std::move(other)); + } + return *this; +} + + +template +Chebyshev::Chebyshev(const Chebyshev& other) + : Chebyshev(other.get_executor()) +{ + *this = other; +} + + +template +Chebyshev::Chebyshev(Chebyshev&& other) + : Chebyshev(other.get_executor()) +{ + *this = std::move(other); +} + + +template +std::unique_ptr Chebyshev::transpose() const +{ + return build() + .with_generated_preconditioner( + share(as(this->get_preconditioner())->transpose())) + .with_criteria(this->get_stop_criterion_factory()) + .with_foci(parameters_.foci) + .on(this->get_executor()) + ->generate( + share(as(this->get_system_matrix())->transpose())); +} + + +template +std::unique_ptr Chebyshev::conj_transpose() const +{ + return build() + .with_generated_preconditioner(share( + as(this->get_preconditioner())->conj_transpose())) + .with_criteria(this->get_stop_criterion_factory()) + .with_foci(conj(std::get<0>(parameters_.foci)), + conj(std::get<1>(parameters_.foci))) + .on(this->get_executor()) + ->generate(share( + as(this->get_system_matrix())->conj_transpose())); +} + + +template +void Chebyshev::apply_impl(const LinOp* b, LinOp* x) const +{ + this->apply_with_initial_guess(b, x, this->get_default_initial_guess()); +} + + +template +void Chebyshev::apply_with_initial_guess_impl( + const LinOp* b, LinOp* x, initial_guess_mode guess) const +{ + if (!this->get_system_matrix()) { + return; + } + experimental::precision_dispatch_real_complex_distributed( + [this, guess](auto dense_b, auto dense_x) { + prepare_initial_guess(dense_b, dense_x, guess); + this->apply_dense_impl(dense_b, dense_x, guess); + }, + b, x); +} + + +template +template +void Chebyshev::apply_dense_impl(const VectorType* dense_b, + VectorType* dense_x, + initial_guess_mode guess) const +{ + using Vector = matrix::Dense; + using ws = workspace_traits; + using coeff_type = detail::coeff_type; + + auto exec = this->get_executor(); + this->setup_workspace(); + + GKO_SOLVER_VECTOR(residual, dense_b); + GKO_SOLVER_VECTOR(inner_solution, dense_b); + GKO_SOLVER_VECTOR(update_solution, dense_b); + + GKO_SOLVER_ONE_MINUS_ONE(); + + auto alpha_ref = coeff_type{1} / center_; + auto beta_ref = coeff_type{0.5} * (foci_direction_ * alpha_ref) * + (foci_direction_ * alpha_ref); + + auto& stop_status = this->template create_workspace_array( + ws::stop, dense_b->get_size()[1]); + exec->run(chebyshev::make_initialize(&stop_status)); + if (guess != initial_guess_mode::zero) { + residual->copy_from(dense_b); + this->get_system_matrix()->apply(neg_one_op, dense_x, one_op, residual); + } + // zero input the residual is dense_b + const VectorType* residual_ptr = + guess == initial_guess_mode::zero ? dense_b : residual; + + auto stop_criterion = this->get_stop_criterion_factory()->generate( + this->get_system_matrix(), + std::shared_ptr(dense_b, [](const LinOp*) {}), dense_x, + residual_ptr); + + int iter = -1; + while (true) { + ++iter; + auto log_func = [this](auto solver, auto dense_b, auto dense_x, + auto iter, auto residual_ptr, + array& stop_status, + bool all_stopped) { + this->template log( + solver, dense_b, dense_x, iter, residual_ptr, nullptr, nullptr, + &stop_status, all_stopped); + }; + bool all_stopped = update_residual( + this, iter, dense_b, dense_x, residual, residual_ptr, + stop_criterion, stop_status, log_func); + if (all_stopped) { + break; + } + + if (this->get_preconditioner()->apply_uses_initial_guess()) { + // Use the inner solver to solve + // A * inner_solution = residual + // with residual as initial guess. + inner_solution->copy_from(residual_ptr); + } + this->get_preconditioner()->apply(residual_ptr, inner_solution); + if (iter == 0) { + // x = x + alpha * inner_solution + // update_solultion = inner_solution + exec->run(chebyshev::make_init_update( + alpha_ref, gko::detail::get_local(inner_solution), + gko::detail::get_local(update_solution), + gko::detail::get_local(dense_x))); + continue; + } + // beta_ref for iter == 1 is initialized in the beginning + if (iter > 1) { + beta_ref = (foci_direction_ * alpha_ref / coeff_type{2.0}) * + (foci_direction_ * alpha_ref / coeff_type{2.0}); + } + alpha_ref = coeff_type{1.0} / (center_ - beta_ref / alpha_ref); + // z = z + beta * p + // p = z + // x += alpha * p + exec->run(chebyshev::make_update( + alpha_ref, beta_ref, gko::detail::get_local(inner_solution), + gko::detail::get_local(update_solution), + gko::detail::get_local(dense_x))); + } +} + + +template +void Chebyshev::apply_impl(const LinOp* alpha, const LinOp* b, + const LinOp* beta, LinOp* x) const +{ + this->apply_with_initial_guess(alpha, b, beta, x, + this->get_default_initial_guess()); +} + +template +void Chebyshev::apply_with_initial_guess_impl( + const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x, + initial_guess_mode guess) const +{ + if (!this->get_system_matrix()) { + return; + } + experimental::precision_dispatch_real_complex_distributed( + [this, guess](auto dense_alpha, auto dense_b, auto dense_beta, + auto dense_x) { + prepare_initial_guess(dense_b, dense_x, guess); + auto x_clone = dense_x->clone(); + this->apply_dense_impl(dense_b, x_clone.get(), guess); + dense_x->scale(dense_beta); + dense_x->add_scaled(dense_alpha, x_clone.get()); + }, + alpha, b, beta, x); +} + + +template +int workspace_traits>::num_arrays(const Solver&) +{ + return 1; +} + + +template +int workspace_traits>::num_vectors(const Solver&) +{ + return 5; +} + + +template +std::vector workspace_traits>::op_names( + const Solver&) +{ + return { + "residual", "inner_solution", "update_solution", "one", "minus_one", + }; +} + + +template +std::vector workspace_traits>::array_names( + const Solver&) +{ + return {"stop"}; +} + + +template +std::vector workspace_traits>::scalars(const Solver&) +{ + return {}; +} + + +template +std::vector workspace_traits>::vectors(const Solver&) +{ + return {residual, inner_solution, update_solution}; +} + + +#define GKO_DECLARE_CHEBYSHEV(_type) class Chebyshev<_type> +#define GKO_DECLARE_CHEBYSHEV_TRAITS(_type) \ + struct workspace_traits> +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV); +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_TRAITS); + + +} // namespace solver +} // namespace gko diff --git a/core/solver/chebyshev_kernels.hpp b/core/solver/chebyshev_kernels.hpp new file mode 100644 index 00000000000..1f8f67c32d5 --- /dev/null +++ b/core/solver/chebyshev_kernels.hpp @@ -0,0 +1,64 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_CORE_SOLVER_CHEBYSHEV_KERNELS_HPP_ +#define GKO_CORE_SOLVER_CHEBYSHEV_KERNELS_HPP_ + + +#include + +#include +#include +#include + +#include "core/base/kernel_declaration.hpp" + + +namespace gko { +namespace kernels { +namespace chebyshev { + + +template +using coeff_type = + std::conditional_t, std::complex, double>; + + +#define GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType) \ + void init_update(std::shared_ptr exec, \ + const coeff_type alpha, \ + const matrix::Dense* inner_sol, \ + matrix::Dense* update_sol, \ + matrix::Dense* output) + +#define GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType) \ + void update(std::shared_ptr exec, \ + const coeff_type alpha, \ + const coeff_type beta, \ + matrix::Dense* inner_sol, \ + matrix::Dense* update_sol, \ + matrix::Dense* output) + +#define GKO_DECLARE_ALL_AS_TEMPLATES \ + template \ + GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType); \ + template \ + GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType) + + +} // namespace chebyshev + + +GKO_DECLARE_FOR_ALL_EXECUTOR_NAMESPACES(chebyshev, + GKO_DECLARE_ALL_AS_TEMPLATES); + + +#undef GKO_DECLARE_ALL_AS_TEMPLATES + + +} // namespace kernels +} // namespace gko + + +#endif // GKO_CORE_SOLVER_CHEBYSHEV_KERNELS_HPP_ diff --git a/core/solver/ir.cpp b/core/solver/ir.cpp index 75efac351f9..69a76c59aa8 100644 --- a/core/solver/ir.cpp +++ b/core/solver/ir.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -13,6 +13,7 @@ #include "core/solver/ir_kernels.hpp" #include "core/solver/solver_base.hpp" #include "core/solver/solver_boilerplate.hpp" +#include "core/solver/update_residual.hpp" namespace gko { @@ -194,7 +195,6 @@ void Ir::apply_dense_impl(const VectorType* dense_b, { using Vector = matrix::Dense; using ws = workspace_traits; - constexpr uint8 relative_stopping_id{1}; auto exec = this->get_executor(); this->setup_workspace(); @@ -204,7 +204,6 @@ void Ir::apply_dense_impl(const VectorType* dense_b, GKO_SOLVER_ONE_MINUS_ONE(); - bool one_changed{}; auto& stop_status = this->template create_workspace_array( ws::stop, dense_b->get_size()[1]); exec->run(ir::make_initialize(&stop_status)); @@ -225,52 +224,19 @@ void Ir::apply_dense_impl(const VectorType* dense_b, while (true) { ++iter; - if (iter == 0) { - // In iter 0, the iteration and residual are updated. - bool all_stopped = stop_criterion->update() - .num_iterations(iter) - .residual(residual_ptr) - .solution(dense_x) - .check(relative_stopping_id, true, - &stop_status, &one_changed); + auto log_func = [this](auto solver, auto dense_b, auto dense_x, + auto iter, auto residual_ptr, + array& stop_status, + bool all_stopped) { this->template log( - this, dense_b, dense_x, iter, residual_ptr, nullptr, nullptr, + solver, dense_b, dense_x, iter, residual_ptr, nullptr, nullptr, &stop_status, all_stopped); - if (all_stopped) { - break; - } - } else { - // In the other iterations, the residual can be updated separately. - bool all_stopped = stop_criterion->update() - .num_iterations(iter) - .solution(dense_x) - // we have the residual check later - .ignore_residual_check(true) - .check(relative_stopping_id, false, - &stop_status, &one_changed); - if (all_stopped) { - this->template log( - this, dense_b, dense_x, iter, nullptr, nullptr, nullptr, - &stop_status, all_stopped); - break; - } - residual_ptr = residual; - // residual = b - A * x - residual->copy_from(dense_b); - this->get_system_matrix()->apply(neg_one_op, dense_x, one_op, - residual); - all_stopped = stop_criterion->update() - .num_iterations(iter) - .residual(residual_ptr) - .solution(dense_x) - .check(relative_stopping_id, true, &stop_status, - &one_changed); - this->template log( - this, dense_b, dense_x, iter, residual_ptr, nullptr, nullptr, - &stop_status, all_stopped); - if (all_stopped) { - break; - } + }; + bool all_stopped = update_residual( + this, iter, dense_b, dense_x, residual, residual_ptr, + stop_criterion, stop_status, log_func); + if (all_stopped) { + break; } if (solver_->apply_uses_initial_guess()) { diff --git a/core/solver/update_residual.hpp b/core/solver/update_residual.hpp new file mode 100644 index 00000000000..b4f94083dcd --- /dev/null +++ b/core/solver/update_residual.hpp @@ -0,0 +1,79 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_CORE_SOLVER_UPDATE_RESIDUAL_HPP_ +#define GKO_CORE_SOLVER_UPDATE_RESIDUAL_HPP_ + + +#include +#include +#include + + +namespace gko { +namespace solver { + + +template +bool update_residual(SolverType* solver, int iter, const VectorType* dense_b, + VectorType* dense_x, VectorType* residual, + const VectorType*& residual_ptr, + std::unique_ptr& stop_criterion, + array& stop_status, LogFunc log) +{ + using ws = workspace_traits>; + constexpr uint8 relative_stopping_id{1}; + + // It's required to be initialized outside. + auto one_op = solver->get_workspace_op(ws::one); + auto neg_one_op = solver->get_workspace_op(ws::minus_one); + + bool one_changed{}; + if (iter == 0) { + // In iter 0, the iteration and residual are updated. + bool all_stopped = + stop_criterion->update() + .num_iterations(iter) + .residual(residual_ptr) + .solution(dense_x) + .check(relative_stopping_id, true, &stop_status, &one_changed); + log(solver, dense_b, dense_x, iter, residual_ptr, stop_status, + all_stopped); + return all_stopped; + } else { + // In the other iterations, the residual can be updated separately. + bool all_stopped = + stop_criterion->update() + .num_iterations(iter) + .solution(dense_x) + // we have the residual check later + .ignore_residual_check(true) + .check(relative_stopping_id, false, &stop_status, &one_changed); + if (all_stopped) { + log(solver, dense_b, dense_x, iter, nullptr, stop_status, + all_stopped); + return all_stopped; + } + residual_ptr = residual; + // residual = b - A * x + residual->copy_from(dense_b); + solver->get_system_matrix()->apply(neg_one_op, dense_x, one_op, + residual); + all_stopped = + stop_criterion->update() + .num_iterations(iter) + .residual(residual_ptr) + .solution(dense_x) + .check(relative_stopping_id, true, &stop_status, &one_changed); + log(solver, dense_b, dense_x, iter, residual_ptr, stop_status, + all_stopped); + return all_stopped; + } +} + + +} // namespace solver +} // namespace gko + +#endif // GKO_CORE_SOLVER_UPDATE_RESIDUAL_HPP_ diff --git a/core/test/config/solver.cpp b/core/test/config/solver.cpp index a170ebb1e04..02560b6c587 100644 --- a/core/test/config/solver.cpp +++ b/core/test/config/solver.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -439,6 +440,41 @@ struct UpperTrs : TrsHelper { }; +struct Chebyshev : SolverConfigTest, + gko::solver::Chebyshev> { + static pnode::map_type setup_base() + { + return {{"type", pnode{"solver::Chebyshev"}}}; + } + + template + static void set(pnode::map_type& config_map, ParamType& param, registry reg, + std::shared_ptr exec) + { + solver_config_test::template set(config_map, param, reg, + exec); + using fvt = typename decltype(param.foci)::first_type; + config_map["foci"] = + pnode{pnode::array_type{pnode{fvt{0.5}}, pnode{fvt{1.5}}}}; + param.with_foci(fvt{0.5}, fvt{1.5}); + config_map["default_initial_guess"] = pnode{"zero"}; + param.with_default_initial_guess(gko::solver::initial_guess_mode::zero); + } + + template + static void validate(gko::LinOpFactory* result, AnswerType* answer) + { + auto res_param = gko::as(result)->get_parameters(); + auto ans_param = answer->get_parameters(); + + solver_config_test::template validate(result, answer); + ASSERT_EQ(res_param.foci, ans_param.foci); + ASSERT_EQ(res_param.default_initial_guess, + ans_param.default_initial_guess); + } +}; + + template class Solver : public ::testing::Test { protected: diff --git a/core/test/solver/CMakeLists.txt b/core/test/solver/CMakeLists.txt index 712d323fcd7..58942b0f956 100644 --- a/core/test/solver/CMakeLists.txt +++ b/core/test/solver/CMakeLists.txt @@ -2,13 +2,14 @@ ginkgo_create_test(batch_bicgstab) ginkgo_create_test(batch_cg) ginkgo_create_test(bicg) ginkgo_create_test(bicgstab) +ginkgo_create_test(cb_gmres) ginkgo_create_test(cg) ginkgo_create_test(cgs) +ginkgo_create_test(chebyshev) ginkgo_create_test(direct) ginkgo_create_test(fcg) ginkgo_create_test(gcr) ginkgo_create_test(gmres) -ginkgo_create_test(cb_gmres) ginkgo_create_test(idr) ginkgo_create_test(ir) ginkgo_create_test(lower_trs) diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp new file mode 100644 index 00000000000..8e10d1f13ff --- /dev/null +++ b/core/test/solver/chebyshev.cpp @@ -0,0 +1,329 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" + + +namespace { + + +template +class Chebyshev : public ::testing::Test { +protected: + using value_type = T; + using Mtx = gko::matrix::Dense; + using Solver = gko::solver::Chebyshev; + + Chebyshev() + : exec(gko::ReferenceExecutor::create()), + mtx(gko::initialize( + {{2, -1.0, 0.0}, {-1.0, 2, -1.0}, {0.0, -1.0, 2}}, exec)), + chebyshev_factory( + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(3u), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value)) + .on(exec)), + solver(chebyshev_factory->generate(mtx)) + {} + + std::shared_ptr exec; + std::shared_ptr mtx; + std::shared_ptr chebyshev_factory; + std::unique_ptr solver; + + static void assert_same_matrices(const Mtx* m1, const Mtx* m2) + { + ASSERT_EQ(m1->get_size()[0], m2->get_size()[0]); + ASSERT_EQ(m1->get_size()[1], m2->get_size()[1]); + for (gko::size_type i = 0; i < m1->get_size()[0]; ++i) { + for (gko::size_type j = 0; j < m2->get_size()[1]; ++j) { + EXPECT_EQ(m1->at(i, j), m2->at(i, j)); + } + } + } +}; + +TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); + + +TYPED_TEST(Chebyshev, ChebyshevFactoryKnowsItsExecutor) +{ + ASSERT_EQ(this->chebyshev_factory->get_executor(), this->exec); +} + + +TYPED_TEST(Chebyshev, ChebyshevFactoryCreatesCorrectSolver) +{ + using Solver = typename TestFixture::Solver; + auto solver = static_cast(this->solver.get()); + + ASSERT_EQ(this->solver->get_size(), gko::dim<2>(3, 3)); + ASSERT_NE(solver->get_system_matrix(), nullptr); + ASSERT_EQ(solver->get_system_matrix(), this->mtx); +} + + +TYPED_TEST(Chebyshev, CanBeCopied) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + auto copy = this->chebyshev_factory->generate(Mtx::create(this->exec)); + + copy->copy_from(this->solver.get()); + + ASSERT_EQ(copy->get_size(), gko::dim<2>(3, 3)); + auto copy_mtx = static_cast(copy.get())->get_system_matrix(); + this->assert_same_matrices(static_cast(copy_mtx.get()), + this->mtx.get()); +} + + +TYPED_TEST(Chebyshev, CanBeMoved) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + auto copy = this->chebyshev_factory->generate(Mtx::create(this->exec)); + + copy->move_from(this->solver); + + ASSERT_EQ(copy->get_size(), gko::dim<2>(3, 3)); + auto copy_mtx = static_cast(copy.get())->get_system_matrix(); + this->assert_same_matrices(static_cast(copy_mtx.get()), + this->mtx.get()); +} + + +TYPED_TEST(Chebyshev, CanBeCloned) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + + auto clone = this->solver->clone(); + + ASSERT_EQ(clone->get_size(), gko::dim<2>(3, 3)); + auto clone_mtx = static_cast(clone.get())->get_system_matrix(); + this->assert_same_matrices(static_cast(clone_mtx.get()), + this->mtx.get()); +} + + +TYPED_TEST(Chebyshev, CanBeCleared) +{ + using Solver = typename TestFixture::Solver; + + this->solver->clear(); + + ASSERT_EQ(this->solver->get_size(), gko::dim<2>(0, 0)); + auto solver_mtx = + static_cast(this->solver.get())->get_system_matrix(); + ASSERT_EQ(solver_mtx, nullptr); +} + + +TYPED_TEST(Chebyshev, DefaultApplyUsesInitialGuess) +{ + ASSERT_TRUE(this->solver->apply_uses_initial_guess()); +} + + +TYPED_TEST(Chebyshev, CanSetEigenRegion) +{ + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + + std::shared_ptr chebyshev_solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .with_foci(value_type{0.2}, value_type{1.2}) + .on(this->exec) + ->generate(this->mtx); + + ASSERT_EQ(chebyshev_solver->get_parameters().foci, + std::make_pair(value_type{0.2}, value_type{1.2})); +} + + +TYPED_TEST(Chebyshev, CanSetInnerSolverInFactory) +{ + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + + auto chebyshev_factory = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .with_preconditioner(Solver::build().with_criteria( + gko::stop::Iteration::build().with_max_iters(3u))) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + auto preconditioner = dynamic_cast( + static_cast(solver.get())->get_preconditioner().get()); + + ASSERT_NE(preconditioner, nullptr); + ASSERT_EQ(preconditioner->get_size(), gko::dim<2>(3, 3)); + ASSERT_EQ(preconditioner->get_system_matrix(), this->mtx); +} + + +TYPED_TEST(Chebyshev, CanSetGeneratedInnerSolverInFactory) +{ + using Solver = typename TestFixture::Solver; + std::shared_ptr chebyshev_solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .on(this->exec) + ->generate(this->mtx); + + auto chebyshev_factory = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .with_generated_preconditioner(chebyshev_solver) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + auto preconditioner = solver->get_preconditioner(); + + ASSERT_NE(preconditioner.get(), nullptr); + ASSERT_EQ(preconditioner.get(), chebyshev_solver.get()); +} + + +TYPED_TEST(Chebyshev, CanSetCriteriaAgain) +{ + using Solver = typename TestFixture::Solver; + std::shared_ptr init_crit = + gko::stop::Iteration::build().with_max_iters(3u).on(this->exec); + auto chebyshev_factory = + Solver::build().with_criteria(init_crit).on(this->exec); + + ASSERT_EQ((chebyshev_factory->get_parameters().criteria).back(), init_crit); + + auto solver = chebyshev_factory->generate(this->mtx); + std::shared_ptr new_crit = + gko::stop::Iteration::build().with_max_iters(5u).on(this->exec); + + solver->set_stop_criterion_factory(new_crit); + auto new_crit_fac = solver->get_stop_criterion_factory(); + auto niter = + static_cast(new_crit_fac.get()) + ->get_parameters() + .max_iters; + + ASSERT_EQ(niter, 5); +} + + +TYPED_TEST(Chebyshev, ThrowsOnWrongInnerSolverInFactory) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + std::shared_ptr wrong_sized_mtx = + Mtx::create(this->exec, gko::dim<2>{2, 2}); + std::shared_ptr chebyshev_solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .on(this->exec) + ->generate(wrong_sized_mtx); + + auto chebyshev_factory = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .with_generated_preconditioner(chebyshev_solver) + .on(this->exec); + + ASSERT_THROW(chebyshev_factory->generate(this->mtx), + gko::DimensionMismatch); +} + + +TYPED_TEST(Chebyshev, CanSetInnerSolver) +{ + using Solver = typename TestFixture::Solver; + std::shared_ptr chebyshev_solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .on(this->exec) + ->generate(this->mtx); + + auto chebyshev_factory = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + solver->set_preconditioner(chebyshev_solver); + auto preconditioner = solver->get_preconditioner(); + + ASSERT_NE(preconditioner.get(), nullptr); + ASSERT_EQ(preconditioner.get(), chebyshev_solver.get()); +} + + +TYPED_TEST(Chebyshev, CanSetApplyWithInitialGuessMode) +{ + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + using initial_guess_mode = gko::solver::initial_guess_mode; + + for (auto guess : {initial_guess_mode::provided, initial_guess_mode::rhs, + initial_guess_mode::zero}) { + auto chebyshev_factory = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .with_default_initial_guess(guess) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + + ASSERT_EQ(solver->apply_uses_initial_guess(), + guess == gko::solver::initial_guess_mode::provided); + } +} + + +TYPED_TEST(Chebyshev, ThrowOnWrongInnerSolverSet) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + std::shared_ptr wrong_sized_mtx = + Mtx::create(this->exec, gko::dim<2>{2, 2}); + std::shared_ptr chebyshev_solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .on(this->exec) + ->generate(wrong_sized_mtx); + + auto chebyshev_factory = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + + ASSERT_THROW(solver->set_preconditioner(chebyshev_solver), + gko::DimensionMismatch); +} + + +TYPED_TEST(Chebyshev, ThrowsOnRectangularMatrixInFactory) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + std::shared_ptr rectangular_mtx = + Mtx::create(this->exec, gko::dim<2>{1, 2}); + + ASSERT_THROW(this->chebyshev_factory->generate(rectangular_mtx), + gko::DimensionMismatch); +} + + +} // namespace diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp new file mode 100644 index 00000000000..890becedde6 --- /dev/null +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -0,0 +1,234 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_PUBLIC_CORE_SOLVER_CHEBYSHEV_HPP_ +#define GKO_PUBLIC_CORE_SOLVER_CHEBYSHEV_HPP_ + + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace gko { +namespace solver { +namespace detail { + + +template +using coeff_type = + std::conditional_t, std::complex, double>; + + +} + +/** + * Chebyshev iteration is an iterative method for solving nonsymmetric problems + * based on some knowledge of the spectrum of the (preconditioned) system + * matrix. It avoids the computation of inner products, which may be a + * performance bottleneck for distributed system. Chebyshev iteration is + * developed based on Chebyshev polynomials of the first kind. + * This implementation follows the algorithm in "Templates for the + * Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd + * Edition". + * + * ``` + * solution = initial_guess + * while not converged: + * residual = b - A solution + * error = preconditioner(A) * residual + * solution = solution + alpha_i * error + beta_i * (solution_i - + * solution_{i-1}) + * ``` + * + * + * @tparam ValueType precision of matrix elements + * + * @ingroup solvers + * @ingroup LinOp + */ +template +class Chebyshev final + : public EnableLinOp>, + public EnablePreconditionedIterativeSolver>, + public EnableApplyWithInitialGuess>, + public Transposable { + friend class EnableLinOp; + friend class EnablePolymorphicObject; + friend class EnableApplyWithInitialGuess; + +public: + using value_type = ValueType; + using transposed_type = Chebyshev; + + std::unique_ptr transpose() const override; + + std::unique_ptr conj_transpose() const override; + + /** + * Return true as iterative solvers use the data in x as an initial guess. + * + * @return true as iterative solvers use the data in x as an initial guess. + */ + bool apply_uses_initial_guess() const override + { + return this->get_default_initial_guess() == + initial_guess_mode::provided; + } + + /** + * Copy-assigns a Chebyshev solver. Preserves the executor, shallow-copies + * inner solver, stopping criterion and system matrix. If the executors + * mismatch, clones inner solver, stopping criterion and system matrix onto + * this executor. + */ + Chebyshev& operator=(const Chebyshev&); + + /** + * Move-assigns a Chebyshev solver. Preserves the executor, moves inner + * solver, stopping criterion and system matrix. If the executors mismatch, + * clones inner solver, stopping criterion and system matrix onto this + * executor. The moved-from object is empty (0x0 and nullptr inner solver, + * stopping criterion and system matrix) + */ + Chebyshev& operator=(Chebyshev&&); + + /** + * Copy-constructs an Chebyshev solver. Inherits the executor, + * shallow-copies inner solver, stopping criterion and system matrix. + */ + Chebyshev(const Chebyshev&); + + /** + * Move-constructs an Chebyshev solver. Preserves the executor, moves inner + * solver, stopping criterion and system matrix. The moved-from object is + * empty (0x0 and nullptr inner solver, stopping criterion and system + * matrix) + */ + Chebyshev(Chebyshev&&); + + class Factory; + + struct parameters_type + : enable_preconditioned_iterative_solver_factory_parameters< + parameters_type, Factory> { + /** + * The pair of foci of ellipse, which covers the eigenvalues of + * preconditioned system. It is usually be {lower bound of eigval, upper + * bound of eigval} of preconditioned real matrices. + */ + std::pair, + detail::coeff_type> + GKO_FACTORY_PARAMETER_VECTOR(foci, + detail::coeff_type{0}, + detail::coeff_type{1}); + + /** + * Default initial guess mode. The available options are under + * initial_guess_mode. + */ + initial_guess_mode GKO_FACTORY_PARAMETER_SCALAR( + default_initial_guess, initial_guess_mode::provided); + }; + + GKO_ENABLE_LIN_OP_FACTORY(Chebyshev, parameters, Factory); + GKO_ENABLE_BUILD_METHOD(Factory); + + /** + * Create the parameters from the property_tree. + * Because this is directly tied to the specific type, the value/index type + * settings within config are ignored and type_descriptor is only used + * for children configs. + * + * @param config the property tree for setting + * @param context the registry + * @param td_for_child the type descriptor for children configs. The + * default uses the value type of this class. + * + * @return parameters + */ + 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; + + template + void apply_dense_impl(const VectorType* b, VectorType* x, + initial_guess_mode guess) const; + + 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; + + explicit Chebyshev(std::shared_ptr exec) + : EnableLinOp(std::move(exec)) + {} + + explicit Chebyshev(const Factory* factory, + std::shared_ptr system_matrix); + +private: + std::shared_ptr solver_{}; + detail::coeff_type center_; + detail::coeff_type foci_direction_; +}; + + +template +struct workspace_traits> { + using Solver = Chebyshev; + // 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 op_names(const Solver&); + // array containing the num_arrays names for the workspace vectors + static std::vector array_names(const Solver&); + // array containing all varying scalar vectors (independent of problem size) + static std::vector scalars(const Solver&); + // array containing all varying vectors (dependent on problem size) + static std::vector vectors(const Solver&); + + // residual vector + constexpr static int residual = 0; + // inner solution vector + constexpr static int inner_solution = 1; + // update solution + constexpr static int update_solution = 2; + // constant 1.0 scalar + constexpr static int one = 3; + // constant -1.0 scalar + constexpr static int minus_one = 4; + + // stopping status array + constexpr static int stop = 0; +}; + + +} // namespace solver +} // namespace gko + + +#endif // GKO_PUBLIC_CORE_SOLVER_CHEBYSHEV_HPP_ diff --git a/include/ginkgo/ginkgo.hpp b/include/ginkgo/ginkgo.hpp index 56c2c0ba6fb..3148df7c7e2 100644 --- a/include/ginkgo/ginkgo.hpp +++ b/include/ginkgo/ginkgo.hpp @@ -139,6 +139,7 @@ #include #include #include +#include #include #include #include diff --git a/reference/CMakeLists.txt b/reference/CMakeLists.txt index 9bea1190ef1..a9eba3f1f83 100644 --- a/reference/CMakeLists.txt +++ b/reference/CMakeLists.txt @@ -54,13 +54,14 @@ target_sources(ginkgo_reference solver/batch_cg_kernels.cpp solver/bicg_kernels.cpp solver/bicgstab_kernels.cpp + solver/cb_gmres_kernels.cpp solver/cg_kernels.cpp solver/cgs_kernels.cpp + solver/chebyshev_kernels.cpp + solver/common_gmres_kernels.cpp solver/fcg_kernels.cpp solver/gcr_kernels.cpp solver/gmres_kernels.cpp - solver/cb_gmres_kernels.cpp - solver/common_gmres_kernels.cpp solver/idr_kernels.cpp solver/ir_kernels.cpp solver/lower_trs_kernels.cpp diff --git a/reference/solver/chebyshev_kernels.cpp b/reference/solver/chebyshev_kernels.cpp new file mode 100644 index 00000000000..a6cfddec2f1 --- /dev/null +++ b/reference/solver/chebyshev_kernels.cpp @@ -0,0 +1,69 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include + +namespace gko { +namespace kernels { +namespace reference { +namespace chebyshev { + + +template +using coeff_type = gko::kernels::chebyshev::coeff_type; + + +template +void init_update(std::shared_ptr exec, + const coeff_type alpha, + const matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + using type = coeff_type; + for (size_t row = 0; row < output->get_size()[0]; row++) { + for (size_t col = 0; col < output->get_size()[1]; col++) { + const auto inner_val = static_cast(inner_sol->at(row, col)); + update_sol->at(row, col) = static_cast(inner_val); + output->at(row, col) = + static_cast(static_cast(output->at(row, col)) + + static_cast(alpha) * inner_val); + } + } +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); + + +template +void update(std::shared_ptr exec, + const coeff_type alpha, const coeff_type beta, + matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + using type = coeff_type; + for (size_t row = 0; row < output->get_size()[0]; row++) { + for (size_t col = 0; col < output->get_size()[1]; col++) { + const auto val = static_cast(inner_sol->at(row, col)) + + static_cast(beta) * + static_cast(update_sol->at(row, col)); + inner_sol->at(row, col) = static_cast(val); + update_sol->at(row, col) = static_cast(val); + output->at(row, col) = + static_cast(static_cast(output->at(row, col)) + + static_cast(alpha) * val); + } + } +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); + + +} // namespace chebyshev +} // namespace reference +} // namespace kernels +} // namespace gko diff --git a/reference/test/solver/CMakeLists.txt b/reference/test/solver/CMakeLists.txt index fa14330681d..269d406086f 100644 --- a/reference/test/solver/CMakeLists.txt +++ b/reference/test/solver/CMakeLists.txt @@ -2,13 +2,14 @@ ginkgo_create_test(batch_bicgstab_kernels) ginkgo_create_test(batch_cg_kernels) ginkgo_create_test(bicg_kernels) ginkgo_create_test(bicgstab_kernels) +ginkgo_create_test(cb_gmres_kernels) ginkgo_create_test(cg_kernels) ginkgo_create_test(cgs_kernels) +ginkgo_create_test(chebyshev_kernels) ginkgo_create_test(direct) ginkgo_create_test(fcg_kernels) ginkgo_create_test(gcr_kernels) ginkgo_create_test(gmres_kernels) -ginkgo_create_test(cb_gmres_kernels) ginkgo_create_test(idr_kernels) ginkgo_create_test(ir_kernels) ginkgo_create_test(lower_trs) diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp new file mode 100644 index 00000000000..d8bd602b86c --- /dev/null +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -0,0 +1,472 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" + +template +class Chebyshev : public ::testing::Test { +protected: + using value_type = T; + using Mtx = gko::matrix::Dense; + using Solver = gko::solver::Chebyshev; + Chebyshev() + : exec(gko::ReferenceExecutor::create()), + mtx(gko::initialize( + {{0.9, -1.0, 3.0}, {0.0, 1.0, 3.0}, {0.0, 0.0, 1.1}}, exec)), + // Eigenvalues of mtx are 0.9, 1.0 and 1.1 + chebyshev_factory( + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(30u), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value)) + .with_foci(value_type{0.9}, value_type{1.1}) + .on(exec)) + {} + + std::shared_ptr exec; + std::shared_ptr mtx; + std::unique_ptr chebyshev_factory; +}; + +TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); + + +TYPED_TEST(Chebyshev, KernelInitUpdate) +{ + using value_type = typename TestFixture::value_type; + using Mtx = typename TestFixture::Mtx; + value_type alpha(0.5); + auto inner_sol = gko::initialize( + {{0.5, 0.125, -0.125}, {0.25, 0.5, -1.0}, {1.5, -0.25, 1.5}}, + this->exec); + auto update_sol = gko::initialize( + {{0.125, 0.0, -0.5}, {0.5, 0.125, -1.0}, {-1.5, -1.25, 1.0}}, + this->exec); + auto output = gko::initialize( + {{-1.0, 0.5, -0.0}, {0.75, 0.25, -1.25}, {1.0, -1.25, 3.0}}, + this->exec); + + gko::kernels::reference::chebyshev::init_update( + this->exec, alpha, inner_sol.get(), update_sol.get(), output.get()); + + GKO_ASSERT_MTX_NEAR(update_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR(output, + gko::initialize({{-0.75, 0.5625, -0.0625}, + {0.875, 0.5, -1.75}, + {1.75, -1.375, 3.75}}, + this->exec), + r::value); +} + + +TYPED_TEST(Chebyshev, KernelUpdate) +{ + using value_type = typename TestFixture::value_type; + using Mtx = typename TestFixture::Mtx; + value_type alpha(0.5); + value_type beta(0.25); + auto inner_sol = gko::initialize( + {{0.5, 0.125, -0.125}, {0.25, 0.5, -1.0}, {1.5, -0.25, 1.5}}, + this->exec); + auto update_sol = gko::initialize( + {{1.0, 0.0, -0.5}, {0.5, 1.0, -1.0}, {-1.5, -0.5, 1.0}}, this->exec); + auto output = gko::initialize( + {{-1.0, 0.5, -0.0}, {0.75, 0.25, -1.25}, {1.0, -1.25, 3.0}}, + this->exec); + + gko::kernels::reference::chebyshev::update(this->exec, alpha, beta, + inner_sol.get(), + update_sol.get(), output.get()); + + GKO_ASSERT_MTX_NEAR(update_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR( + inner_sol, + gko::initialize( + {{0.75, 0.125, -0.25}, {0.375, 0.75, -1.25}, {1.125, -0.375, 1.75}}, + this->exec), + r::value); + GKO_ASSERT_MTX_NEAR(output, + gko::initialize({{-0.625, 0.5625, -0.125}, + {0.9375, 0.625, -1.875}, + {1.5625, -1.4375, 3.875}}, + this->exec), + r::value); +} + + +#ifdef GINKGO_MIXED_PRECISION + + +TYPED_TEST(Chebyshev, MixedKernelInitUpdate) +{ + using value_type = typename TestFixture::value_type; + using scalar_type = gko::next_precision; + using Mtx = typename TestFixture::Mtx; + scalar_type alpha(0.5); + auto inner_sol = gko::initialize( + {{0.5, 0.125, -0.125}, {0.25, 0.5, -1.0}, {1.5, -0.25, 1.5}}, + this->exec); + auto update_sol = gko::initialize( + {{0.125, 0.0, -0.5}, {0.5, 0.125, -1.0}, {-1.5, -1.25, 1.0}}, + this->exec); + auto output = gko::initialize( + {{-1.0, 0.5, -0.0}, {0.75, 0.25, -1.25}, {1.0, -1.25, 3.0}}, + this->exec); + + gko::kernels::reference::chebyshev::init_update( + this->exec, alpha, inner_sol.get(), update_sol.get(), output.get()); + + GKO_ASSERT_MTX_NEAR(update_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR(output, + gko::initialize({{-0.75, 0.5625, -0.0625}, + {0.875, 0.5, -1.75}, + {1.75, -1.375, 3.75}}, + this->exec), + (r_mixed())); +} + + +TYPED_TEST(Chebyshev, MixedKernelUpdate) +{ + using value_type = typename TestFixture::value_type; + using scalar_type = gko::next_precision; + using Mtx = typename TestFixture::Mtx; + value_type alpha(0.5); + value_type beta(0.25); + auto inner_sol = gko::initialize( + {{0.5, 0.125, -0.125}, {0.25, 0.5, -1.0}, {1.5, -0.25, 1.5}}, + this->exec); + auto update_sol = gko::initialize( + {{1.0, 0.0, -0.5}, {0.5, 1.0, -1.0}, {-1.5, -0.5, 1.0}}, this->exec); + auto output = gko::initialize( + {{-1.0, 0.5, -0.0}, {0.75, 0.25, -1.25}, {1.0, -1.25, 3.0}}, + this->exec); + + gko::kernels::reference::chebyshev::update(this->exec, alpha, beta, + inner_sol.get(), + update_sol.get(), output.get()); + + GKO_ASSERT_MTX_NEAR(update_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR( + inner_sol, + gko::initialize( + {{0.75, 0.125, -0.25}, {0.375, 0.75, -1.25}, {1.125, -0.375, 1.75}}, + this->exec), + (r_mixed())); + GKO_ASSERT_MTX_NEAR(output, + gko::initialize({{-0.625, 0.5625, -0.125}, + {0.9375, 0.625, -1.875}, + {1.5625, -1.4375, 3.875}}, + this->exec), + (r_mixed())); +} + + +#endif + + +TYPED_TEST(Chebyshev, SolvesTriangularSystem) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.0, 0.0, 0.0}, this->exec); + + solver->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.0, 3.0, 2.0}), r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemMixed) +{ + using mixed_type = gko::next_precision; + using MixedMtx = gko::matrix::Dense; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.0, 0.0, 0.0}, this->exec); + + solver->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.0, 3.0, 2.0}), + (r_mixed()) * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemComplex) +{ + using Mtx = gko::to_complex; + using value_type = typename Mtx::value_type; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto b = gko::initialize( + {value_type{3.9, -7.8}, value_type{9.0, -18.0}, value_type{2.2, -4.4}}, + this->exec); + auto x = gko::initialize( + {value_type{0.0, 0.0}, value_type{0.0, 0.0}, value_type{0.0, 0.0}}, + this->exec); + + solver->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, + l({value_type{1.0, -2.0}, value_type{3.0, -6.0}, + value_type{2.0, -4.0}}), + r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemMixedComplex) +{ + using mixed_complex_type = + gko::to_complex>; + using MixedMtx = gko::matrix::Dense; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto b = gko::initialize( + {mixed_complex_type{3.9, -7.8}, mixed_complex_type{9.0, -18.0}, + mixed_complex_type{2.2, -4.4}}, + this->exec); + auto x = gko::initialize( + {mixed_complex_type{0.0, 0.0}, mixed_complex_type{0.0, 0.0}, + mixed_complex_type{0.0, 0.0}}, + this->exec); + + solver->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR( + x, + l({mixed_complex_type{1.0, -2.0}, mixed_complex_type{3.0, -6.0}, + mixed_complex_type{2.0, -4.0}}), + (r_mixed()) * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemWithIterativeInnerSolver) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + const gko::remove_complex inner_reduction_factor = 1e-2; + auto precond_factory = gko::share( + gko::solver::Gmres::build() + .with_criteria(gko::stop::ResidualNorm::build() + .with_reduction_factor(inner_reduction_factor)) + .on(this->exec)); + auto solver_factory = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(30u), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value)) + .with_preconditioner(precond_factory) + .with_foci(value_type{0.9}, value_type{1.1}) + .on(this->exec); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.0, 0.0, 0.0}, this->exec); + + solver_factory->generate(this->mtx)->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.0, 3.0, 2.0}), r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesMultipleTriangularSystems) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + using T = value_type; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto b = gko::initialize( + {I{3.9, 2.9}, I{9.0, 4.0}, I{2.2, 1.1}}, this->exec); + auto x = gko::initialize( + {I{0.0, 0.0}, I{0.0, 0.0}, I{0.0, 0.0}}, this->exec); + + solver->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({{1.0, 1.0}, {3.0, 1.0}, {2.0, 1.0}}), + r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemUsingAdvancedApply) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto alpha = gko::initialize({2.0}, this->exec); + auto beta = gko::initialize({-1.0}, this->exec); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.5, 1.0, 2.0}, this->exec); + + solver->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.5, 5.0, 2.0}), r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemUsingAdvancedApplyMixed) +{ + using mixed_type = gko::next_precision; + using MixedMtx = gko::matrix::Dense; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto alpha = gko::initialize({2.0}, this->exec); + auto beta = gko::initialize({-1.0}, this->exec); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.5, 1.0, 2.0}, this->exec); + + solver->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.5, 5.0, 2.0}), + (r_mixed()) * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemUsingAdvancedApplyComplex) +{ + using Scalar = typename TestFixture::Mtx; + using Mtx = gko::to_complex; + using value_type = typename Mtx::value_type; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto alpha = gko::initialize({2.0}, this->exec); + auto beta = gko::initialize({-1.0}, this->exec); + auto b = gko::initialize( + {value_type{3.9, -7.8}, value_type{9.0, -18.0}, value_type{2.2, -4.4}}, + this->exec); + auto x = gko::initialize( + {value_type{0.5, -1.0}, value_type{1.0, -2.0}, value_type{2.0, -4.0}}, + this->exec); + + solver->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, + l({value_type{1.5, -3.0}, value_type{5.0, -10.0}, + value_type{2.0, -4.0}}), + (r_mixed()) * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTriangularSystemUsingAdvancedApplyMixedComplex) +{ + using mixed_type = gko::next_precision; + using mixed_complex_type = gko::to_complex; + using Scalar = gko::matrix::Dense; + using MixedMtx = gko::matrix::Dense; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto alpha = gko::initialize({2.0}, this->exec); + auto beta = gko::initialize({-1.0}, this->exec); + auto b = gko::initialize( + {mixed_complex_type{3.9, -7.8}, mixed_complex_type{9.0, -18.0}, + mixed_complex_type{2.2, -4.4}}, + this->exec); + auto x = gko::initialize( + {mixed_complex_type{0.5, -1.0}, mixed_complex_type{1.0, -2.0}, + mixed_complex_type{2.0, -4.0}}, + this->exec); + + solver->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR( + x, + l({mixed_complex_type{1.5, -3.0}, mixed_complex_type{5.0, -10.0}, + mixed_complex_type{2.0, -4.0}}), + (r_mixed()) * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesMultipleStencilSystemsUsingAdvancedApply) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + using T = value_type; + auto solver = this->chebyshev_factory->generate(this->mtx); + auto alpha = gko::initialize({2.0}, this->exec); + auto beta = gko::initialize({-1.0}, this->exec); + auto b = gko::initialize( + {I{3.9, 2.9}, I{9.0, 4.0}, I{2.2, 1.1}}, this->exec); + auto x = gko::initialize( + {I{0.5, 1.0}, I{1.0, 2.0}, I{2.0, 3.0}}, this->exec); + + solver->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({{1.5, 1.0}, {5.0, 0.0}, {2.0, -1.0}}), + r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesTransposedTriangularSystem) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + auto solver = this->chebyshev_factory->generate(this->mtx->transpose()); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.0, 0.0, 0.0}, this->exec); + + solver->transpose()->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.0, 3.0, 2.0}), r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, SolvesConjTransposedTriangularSystem) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + auto solver = + this->chebyshev_factory->generate(this->mtx->conj_transpose()); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.0, 0.0, 0.0}, this->exec); + + solver->conj_transpose()->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({1.0, 3.0, 2.0}), r::value * 1e1); +} + + +TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) +{ + using Mtx = typename TestFixture::Mtx; + using value_type = typename TestFixture::value_type; + using initial_guess_mode = gko::solver::initial_guess_mode; + auto ref_solver = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(1u)) + .with_foci(value_type{0.9}, value_type{1.1}) + .on(this->exec) + ->generate(this->mtx); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + for (auto guess : {initial_guess_mode::provided, initial_guess_mode::rhs, + initial_guess_mode::zero}) { + auto solver = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(1u)) + .with_foci(value_type{0.9}, value_type{1.1}) + .with_default_initial_guess(guess) + .on(this->exec) + ->generate(this->mtx); + auto x = gko::initialize({1.0, -1.0, 1.0}, this->exec); + std::shared_ptr ref_x = nullptr; + if (guess == initial_guess_mode::provided) { + ref_x = x->clone(); + } else if (guess == initial_guess_mode::rhs) { + ref_x = b->clone(); + } else { + ref_x = gko::initialize({0.0, 0.0, 0.0}, this->exec); + } + solver->apply(b, x); + ref_solver->apply(b, ref_x); + + GKO_ASSERT_MTX_NEAR(x, ref_x, 0.0); + } +} diff --git a/reference/test/solver/ir_kernels.cpp b/reference/test/solver/ir_kernels.cpp index b0c1029f693..cb0ff5b5751 100644 --- a/reference/test/solver/ir_kernels.cpp +++ b/reference/test/solver/ir_kernels.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -203,17 +203,18 @@ TYPED_TEST(Ir, SolvesTriangularSystemUsingAdvancedApply) TYPED_TEST(Ir, SolvesTriangularSystemUsingAdvancedApplyMixed) { - using Mtx = typename TestFixture::Mtx; - using value_type = typename TestFixture::value_type; + using mixed_type = gko::next_precision; + using MixedMtx = gko::matrix::Dense; auto solver = this->ir_factory->generate(this->mtx); - auto alpha = gko::initialize({2.0}, this->exec); - auto beta = gko::initialize({-1.0}, this->exec); - auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); - auto x = gko::initialize({0.5, 1.0, 2.0}, this->exec); + auto alpha = gko::initialize({2.0}, this->exec); + auto beta = gko::initialize({-1.0}, this->exec); + auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); + auto x = gko::initialize({0.5, 1.0, 2.0}, this->exec); solver->apply(alpha, b, beta, x); - GKO_ASSERT_MTX_NEAR(x, l({1.5, 5.0, 2.0}), r::value * 1e1); + GKO_ASSERT_MTX_NEAR(x, l({1.5, 5.0, 2.0}), + (r_mixed()) * 1e1); } @@ -243,26 +244,29 @@ TYPED_TEST(Ir, SolvesTriangularSystemUsingAdvancedApplyComplex) TYPED_TEST(Ir, SolvesTriangularSystemUsingAdvancedApplyMixedComplex) { - using Scalar = gko::matrix::Dense< - gko::next_precision>; - using Mtx = gko::to_complex; - using value_type = typename Mtx::value_type; + using mixed_type = gko::next_precision; + using mixed_complex_type = gko::to_complex; + using Scalar = gko::matrix::Dense; + using MixedMtx = gko::matrix::Dense; auto solver = this->ir_factory->generate(this->mtx); auto alpha = gko::initialize({2.0}, this->exec); auto beta = gko::initialize({-1.0}, this->exec); - auto b = gko::initialize( - {value_type{3.9, -7.8}, value_type{9.0, -18.0}, value_type{2.2, -4.4}}, + auto b = gko::initialize( + {mixed_complex_type{3.9, -7.8}, mixed_complex_type{9.0, -18.0}, + mixed_complex_type{2.2, -4.4}}, this->exec); - auto x = gko::initialize( - {value_type{0.5, -1.0}, value_type{1.0, -2.0}, value_type{2.0, -4.0}}, + auto x = gko::initialize( + {mixed_complex_type{0.5, -1.0}, mixed_complex_type{1.0, -2.0}, + mixed_complex_type{2.0, -4.0}}, this->exec); - solver->apply(alpha, b, beta, x); + solver->apply(alpha.get(), b.get(), beta.get(), x.get()); - GKO_ASSERT_MTX_NEAR(x, - l({value_type{1.5, -3.0}, value_type{5.0, -10.0}, - value_type{2.0, -4.0}}), - r::value * 1e1); + GKO_ASSERT_MTX_NEAR( + x, + l({mixed_complex_type{1.5, -3.0}, mixed_complex_type{5.0, -10.0}, + mixed_complex_type{2.0, -4.0}}), + (r_mixed()) * 1e1); } diff --git a/test/solver/CMakeLists.txt b/test/solver/CMakeLists.txt index b24e063bc6d..db58d0b56b4 100644 --- a/test/solver/CMakeLists.txt +++ b/test/solver/CMakeLists.txt @@ -5,6 +5,7 @@ ginkgo_create_common_test(bicgstab_kernels) ginkgo_create_common_test(cb_gmres_kernels) ginkgo_create_common_test(cg_kernels) ginkgo_create_common_test(cgs_kernels) +ginkgo_create_common_test(chebyshev_kernels) ginkgo_create_common_test(direct DISABLE_EXECUTORS dpcpp) ginkgo_create_common_test(fcg_kernels) ginkgo_create_common_test(gcr_kernels) diff --git a/test/solver/chebyshev_kernels.cpp b/test/solver/chebyshev_kernels.cpp new file mode 100644 index 00000000000..55c71bdc68b --- /dev/null +++ b/test/solver/chebyshev_kernels.cpp @@ -0,0 +1,237 @@ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" +#include "test/utils/common_fixture.hpp" +#include "test/utils/executor.hpp" + + +class Chebyshev : public CommonTestFixture { +protected: + using Mtx = gko::matrix::Dense; + + Chebyshev() : rand_engine(30) {} + + std::unique_ptr gen_mtx(gko::size_type num_rows, + gko::size_type num_cols, gko::size_type stride) + { + auto tmp_mtx = gko::test::generate_random_matrix( + num_rows, num_cols, + std::uniform_int_distribution<>(num_cols, num_cols), + std::normal_distribution(-1.0, 1.0), rand_engine, ref); + auto result = Mtx::create(ref, gko::dim<2>{num_rows, num_cols}, stride); + result->copy_from(tmp_mtx); + return result; + } + + std::default_random_engine rand_engine; +}; + + +TEST_F(Chebyshev, KernelInitUpdate) +{ + value_type alpha(0.5); + auto inner_sol = gen_mtx(30, 3, 3); + auto update_sol = gen_mtx(30, 3, 3); + auto output = gen_mtx(30, 3, 3); + auto d_inner_sol = gko::clone(exec, inner_sol); + auto d_update_sol = gko::clone(exec, update_sol); + auto d_output = gko::clone(exec, output); + + gko::kernels::reference::chebyshev::init_update( + ref, alpha, inner_sol.get(), update_sol.get(), output.get()); + gko::kernels::GKO_DEVICE_NAMESPACE::chebyshev::init_update( + exec, alpha, d_inner_sol.get(), d_update_sol.get(), d_output.get()); + + GKO_ASSERT_MTX_NEAR(d_update_sol, d_inner_sol, 0); + GKO_ASSERT_MTX_NEAR(d_update_sol, update_sol, 0); + GKO_ASSERT_MTX_NEAR(d_inner_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR(d_output, output, r::value); +} + + +TEST_F(Chebyshev, KernelUpdate) +{ + value_type alpha(0.5); + value_type beta(0.25); + auto inner_sol = gen_mtx(30, 3, 3); + auto update_sol = gen_mtx(30, 3, 3); + auto output = gen_mtx(30, 3, 3); + auto d_inner_sol = gko::clone(exec, inner_sol); + auto d_update_sol = gko::clone(exec, update_sol); + auto d_output = gko::clone(exec, output); + + gko::kernels::reference::chebyshev::update( + ref, alpha, beta, inner_sol.get(), update_sol.get(), output.get()); + gko::kernels::GKO_DEVICE_NAMESPACE::chebyshev::update( + exec, alpha, beta, d_inner_sol.get(), d_update_sol.get(), + d_output.get()); + + GKO_ASSERT_MTX_NEAR(d_update_sol, d_inner_sol, 0); + GKO_ASSERT_MTX_NEAR(d_inner_sol, inner_sol, r::value); + GKO_ASSERT_MTX_NEAR(d_output, output, r::value); +} + + +#ifdef GINKGO_MIXED_PRECISION + + +TEST_F(Chebyshev, MixedKernelInitUpdate) +{ + using scalar_type = gko::next_precision; + scalar_type alpha(0.5); + auto inner_sol = gen_mtx(30, 3, 3); + auto update_sol = gen_mtx(30, 3, 3); + auto output = gen_mtx(30, 3, 3); + auto d_inner_sol = gko::clone(exec, inner_sol); + auto d_update_sol = gko::clone(exec, update_sol); + auto d_output = gko::clone(exec, output); + + gko::kernels::reference::chebyshev::init_update( + ref, alpha, inner_sol.get(), update_sol.get(), output.get()); + gko::kernels::GKO_DEVICE_NAMESPACE::chebyshev::init_update( + exec, alpha, d_inner_sol.get(), d_update_sol.get(), d_output.get()); + + GKO_ASSERT_MTX_NEAR(d_update_sol, d_inner_sol, 0); + GKO_ASSERT_MTX_NEAR(d_update_sol, update_sol, 0); + GKO_ASSERT_MTX_NEAR(d_inner_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR(d_output, output, (r_mixed())); +} + + +TEST_F(Chebyshev, MixedKernelUpdate) +{ + using scalar_type = gko::next_precision; + value_type alpha(0.5); + value_type beta(0.25); + auto inner_sol = gen_mtx(30, 3, 3); + auto update_sol = gen_mtx(30, 3, 3); + auto output = gen_mtx(30, 3, 3); + auto d_inner_sol = gko::clone(exec, inner_sol); + auto d_update_sol = gko::clone(exec, update_sol); + auto d_output = gko::clone(exec, output); + + gko::kernels::reference::chebyshev::update( + ref, alpha, beta, inner_sol.get(), update_sol.get(), output.get()); + gko::kernels::GKO_DEVICE_NAMESPACE::chebyshev::update( + exec, alpha, beta, d_inner_sol.get(), d_update_sol.get(), + d_output.get()); + + GKO_ASSERT_MTX_NEAR(d_update_sol, d_inner_sol, 0); + GKO_ASSERT_MTX_NEAR(d_inner_sol, inner_sol, r::value); + GKO_ASSERT_MTX_NEAR(d_output, output, r::value); +} + + +#endif + + +TEST_F(Chebyshev, ApplyIsEquivalentToRef) +{ + auto mtx = gen_mtx(50, 50, 52); + auto x = gen_mtx(50, 3, 8); + auto b = gen_mtx(50, 3, 5); + auto d_mtx = gko::clone(exec, mtx); + auto d_x = gko::clone(exec, x); + auto d_b = gko::clone(exec, b); + // Forget about accuracy - Chebyshev is not going to converge for a random + // matrix, just check that a couple of iterations gives the same result on + // both executors + auto chebyshev_factory = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .on(ref); + auto d_chebyshev_factory = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .on(exec); + auto solver = chebyshev_factory->generate(std::move(mtx)); + auto d_solver = d_chebyshev_factory->generate(std::move(d_mtx)); + + solver->apply(b, x); + d_solver->apply(d_b, d_x); + + GKO_ASSERT_MTX_NEAR(d_x, x, r::value); +} + + +TEST_F(Chebyshev, ApplyWithIterativeInnerSolverIsEquivalentToRef) +{ + auto mtx = gen_mtx(50, 50, 54); + auto x = gen_mtx(50, 3, 6); + auto b = gen_mtx(50, 3, 10); + auto d_mtx = gko::clone(exec, mtx); + auto d_x = gko::clone(exec, x); + auto d_b = gko::clone(exec, b); + auto chebyshev_factory = + gko::solver::Chebyshev::build() + .with_preconditioner( + gko::solver::Gmres::build().with_criteria( + gko::stop::Iteration::build().with_max_iters(1u))) + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .on(ref); + auto d_chebyshev_factory = + gko::solver::Chebyshev::build() + .with_preconditioner( + gko::solver::Gmres::build().with_criteria( + gko::stop::Iteration::build().with_max_iters(1u))) + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .on(exec); + auto solver = chebyshev_factory->generate(std::move(mtx)); + auto d_solver = d_chebyshev_factory->generate(std::move(d_mtx)); + + solver->apply(b, x); + d_solver->apply(d_b, d_x); + + // Note: r::value * 300 instead of r::value, as + // the difference in the inner gmres iteration gets amplified by the + // difference in IR. + GKO_ASSERT_MTX_NEAR(d_x, x, r::value * 300); +} + + +TEST_F(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) +{ + using initial_guess_mode = gko::solver::initial_guess_mode; + auto mtx = gko::share(gen_mtx(50, 50, 52)); + auto b = gen_mtx(50, 3, 7); + auto d_mtx = gko::share(clone(exec, mtx)); + auto d_b = gko::clone(exec, b); + for (auto guess : {initial_guess_mode::provided, initial_guess_mode::rhs, + initial_guess_mode::zero}) { + auto x = gen_mtx(50, 3, 4); + auto d_x = gko::clone(exec, x); + auto chebyshev_factory = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .with_default_initial_guess(guess) + .on(ref); + auto d_chebyshev_factory = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) + .with_default_initial_guess(guess) + .on(exec); + auto solver = chebyshev_factory->generate(mtx); + auto d_solver = d_chebyshev_factory->generate(d_mtx); + + solver->apply(b, x); + d_solver->apply(d_b, d_x); + + GKO_ASSERT_MTX_NEAR(d_x, x, r::value); + } +} diff --git a/test/solver/solver.cpp b/test/solver/solver.cpp index 57e93295940..5a3d230e2f1 100644 --- a/test/solver/solver.cpp +++ b/test/solver/solver.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -176,6 +177,21 @@ struct Ir : SimpleSolverTest> { }; +struct Chebyshev : SimpleSolverTest> { + static double tolerance() { return 1e7 * r::value; } + + static typename solver_type::parameters_type build_preconditioned( + std::shared_ptr exec, + gko::size_type iteration_count, bool check_residual = true) + { + return SimpleSolverTest>:: + build(exec, iteration_count, check_residual) + .with_preconditioner( + precond_type::build().with_max_block_size(1u).on(exec)); + } +}; + + template struct CbGmres : SimpleSolverTest> { static constexpr bool will_not_allocate() { return false; } @@ -887,9 +903,9 @@ using SolverTypes = ::testing::Types, Idr<4>,*/ - Ir, CbGmres<2>, CbGmres<10>, Gmres<2>, Gmres<10>, - FGmres<2>, FGmres<10>, Gcr<2>, Gcr<10>, LowerTrs, UpperTrs, - LowerTrsUnitdiag, UpperTrsUnitdiag + Ir, Chebyshev, CbGmres<2>, CbGmres<10>, Gmres<2>, + Gmres<10>, FGmres<2>, FGmres<10>, Gcr<2>, Gcr<10>, + LowerTrs, UpperTrs, LowerTrsUnitdiag, UpperTrsUnitdiag #ifdef GKO_COMPILING_CUDA , LowerTrsSyncfree, UpperTrsSyncfree, diff --git a/test/test_install/test_install.cpp b/test/test_install/test_install.cpp index 2f4cdeda6e4..2d4ba7a5b77 100644 --- a/test/test_install/test_install.cpp +++ b/test/test_install/test_install.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -462,6 +462,16 @@ int main() check_solver(exec, A_raw, b, x); } + // core/solver/chebyshev.hpp + { + using Solver = gko::solver::Chebyshev<>; + auto test = + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(1u).on(exec)) + .on(exec); + } + // core/solver/fcg.hpp { using Solver = gko::solver::Fcg<>;