diff --git a/common/unified/CMakeLists.txt b/common/unified/CMakeLists.txt index 8ed33f1c201..d35f62b20cb 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..d78b1cc888b --- /dev/null +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -0,0 +1,114 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include + +#include +#include +#include + +#include "common/unified/base/kernel_launch.hpp" + + +namespace gko { +namespace kernels { +namespace GKO_DEVICE_NAMESPACE { +namespace chebyshev { + + +#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 +using if_single_only_type = xstd::type_identity_t; + + +#endif + + +template +void init_update(std::shared_ptr exec, + const solver::detail::coeff_type alpha, + const matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + using coeff_type = + if_single_only_type>; + // the coeff_type always be the highest precision, so we need + // to cast the others from ValueType to this precision. + using arithmetic_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)) + + 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 solver::detail::coeff_type alpha, + const solver::detail::coeff_type beta, + matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + using coeff_type = + if_single_only_type>; + // the coeff_type always be the highest precision, so we need + // to cast the others from ValueType to this precision. + using arithmetic_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)) + + 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)) + 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 57fe85df222..815545e5f8b 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -93,8 +93,8 @@ target_sources( 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 @@ -113,6 +113,7 @@ target_sources( 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 d99d9ff0ff7..61768bbce16 100644 --- a/core/config/config_helper.hpp +++ b/core/config/config_helper.hpp @@ -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) \ @@ -53,6 +52,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 6647d3f3a30..96f10119217 100644 --- a/core/config/registry.cpp +++ b/core/config/registry.cpp @@ -33,6 +33,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 7de23a49738..ef4a3cd5637 100644 --- a/core/config/solver_config.cpp +++ b/core/config/solver_config.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -45,6 +46,7 @@ GKO_PARSE_VALUE_TYPE(Minres, gko::solver::Minres); 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 f5b41ff19c9..1abe27e9558 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" @@ -677,6 +678,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..2021b2a99be --- /dev/null +++ b/core/solver/chebyshev.cpp @@ -0,0 +1,371 @@ +// SPDX-FileCopyrightText: 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 ir { +namespace { + + +GKO_REGISTER_OPERATION(initialize, ir::initialize); + + +} +} // namespace ir + + +namespace chebyshev { +namespace { + + +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(std::shared_ptr exec) + : EnableLinOp(std::move(exec)) +{} + + +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); + auto left_foci = std::get<0>(parameters_.foci); + auto right_foci = std::get<1>(parameters_.foci); + GKO_ASSERT(real(left_foci) <= real(right_foci)); + GKO_ASSERT(is_nonzero(left_foci) || is_nonzero(right_foci)); + center_ = + (left_foci + right_foci) / solver::detail::coeff_type{2}; + foci_direction_ = + (right_foci - left_foci) / solver::detail::coeff_type{2}; + // if the center is zero then the alpha will be inf + GKO_ASSERT(is_nonzero(center_)); +} + + +template +Chebyshev& Chebyshev::operator=(const Chebyshev& other) +{ + if (&other != this) { + EnableLinOp::operator=(other); + EnablePreconditionedIterativeSolver< + ValueType, Chebyshev>::operator=(other); + this->parameters_ = other.parameters_; + this->center_ = other.center_; + this->foci_direction_ = other.foci_direction_; + } + return *this; +} + + +template +Chebyshev& Chebyshev::operator=(Chebyshev&& other) +{ + if (&other != this) { + EnableLinOp::operator=(std::move(other)); + EnablePreconditionedIterativeSolver< + ValueType, Chebyshev>::operator=(std::move(other)); + this->parameters_ = std::exchange(other.parameters_, parameters_type{}); + this->center_ = std::exchange(other.center_, + solver::detail::coeff_type{}); + this->foci_direction_ = std::exchange( + other.foci_direction_, solver::detail::coeff_type{}); + } + 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 = solver::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_host = coeff_type{1} / center_; + auto beta_host = coeff_type{0.5} * (foci_direction_ * alpha_host) * + (foci_direction_ * alpha_host); + + auto& stop_status = this->template create_workspace_array( + ws::stop, dense_b->get_size()[1]); + exec->run(ir::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_host, gko::detail::get_local(inner_solution), + gko::detail::get_local(update_solution), + gko::detail::get_local(dense_x))); + continue; + } + // beta_host for iter == 1 is initialized in the beginning + if (iter > 1) { + beta_host = (foci_direction_ * alpha_host / coeff_type{2.0}) * + (foci_direction_ * alpha_host / coeff_type{2.0}); + } + alpha_host = coeff_type{1.0} / (center_ - beta_host / alpha_host); + // z = z + beta * p + // p = z + // x += alpha * p + exec->run(chebyshev::make_update( + alpha_host, beta_host, 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..9e2af545a6e --- /dev/null +++ b/core/solver/chebyshev_kernels.hpp @@ -0,0 +1,60 @@ +// SPDX-FileCopyrightText: 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 + +#include "core/base/kernel_declaration.hpp" + + +namespace gko { +namespace kernels { +namespace chebyshev { + + +#define GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType) \ + void init_update(std::shared_ptr exec, \ + const solver::detail::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 solver::detail::coeff_type alpha, \ + const solver::detail::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 074c4a7d572..c51365f0f0a 100644 --- a/core/solver/ir.cpp +++ b/core/solver/ir.cpp @@ -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 { @@ -192,7 +193,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(); @@ -202,7 +202,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)); @@ -223,52 +222,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..688dcdb4b18 --- /dev/null +++ b/core/solver/update_residual.hpp @@ -0,0 +1,80 @@ +// SPDX-FileCopyrightText: 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 +#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 e23f45fd71c..5f829ed4c75 100644 --- a/core/test/config/solver.cpp +++ b/core/test/config/solver.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -449,6 +450,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..5c2152d128e --- /dev/null +++ b/core/test/solver/chebyshev.cpp @@ -0,0 +1,314 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" +#include "core/test/utils/assertions.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( + {{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; +}; + +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(); + GKO_ASSERT_MTX_NEAR(static_cast(copy_mtx.get()), + this->mtx.get(), 0.0); +} + + +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(); + GKO_ASSERT_MTX_NEAR(static_cast(copy_mtx.get()), + this->mtx.get(), 0.0); +} + + +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(); + GKO_ASSERT_MTX_NEAR(static_cast(clone_mtx.get()), + this->mtx.get(), 0.0); +} + + +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; + using coeff_type = gko::solver::detail::coeff_type; + + std::shared_ptr chebyshev_solver = + Solver::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) + .with_foci(coeff_type{0.2}, coeff_type{1.2}) + .on(this->exec) + ->generate(this->mtx); + + ASSERT_EQ(chebyshev_solver->get_parameters().foci, + std::make_pair(coeff_type{0.2}, coeff_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); +} diff --git a/include/ginkgo/core/base/std_extensions.hpp b/include/ginkgo/core/base/std_extensions.hpp index a950fcc2003..70c9b605ae2 100644 --- a/include/ginkgo/core/base/std_extensions.hpp +++ b/include/ginkgo/core/base/std_extensions.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 @@ -105,6 +105,16 @@ template using conjunction = std::conjunction; +// Provide the type_identity from C++20 +template +struct type_identity { + using type = T; +}; + +template +using type_identity_t = typename type_identity::type; + + } // namespace xstd } // namespace gko diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp new file mode 100644 index 00000000000..cc04282e624 --- /dev/null +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -0,0 +1,234 @@ +// SPDX-FileCopyrightText: 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 a pair {lower bound of eigval, + * upper bound of eigval} of the preconditioned system if the + * preconditioned system only contains non-complex eigenvalues. The foci + * value must satisfy real(foci(1)) >= real(foci(0)). + */ + 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); + + 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/reference/CMakeLists.txt b/reference/CMakeLists.txt index cbc6e70b466..87858d18812 100644 --- a/reference/CMakeLists.txt +++ b/reference/CMakeLists.txt @@ -55,13 +55,14 @@ target_sources( 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..b4f17309d86 --- /dev/null +++ b/reference/solver/chebyshev_kernels.cpp @@ -0,0 +1,72 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include +#include + +namespace gko { +namespace kernels { +namespace reference { +namespace chebyshev { + + +template +void init_update(std::shared_ptr exec, + const solver::detail::coeff_type alpha, + const matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + // the coeff_type always be the highest precision, so we need + // to cast the others from ValueType to this precision. + using arithmetic_type = solver::detail::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)) + + alpha * inner_val); + } + } +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); + + +template +void update(std::shared_ptr exec, + const solver::detail::coeff_type alpha, + const solver::detail::coeff_type beta, + matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + // the coeff_type always be the highest precision, so we need + // to cast the others from ValueType to this precision. + using arithmetic_type = solver::detail::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)) + + 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)) + + 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 ee527587881..14a0d4cee8b 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..457555bb640 --- /dev/null +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -0,0 +1,401 @@ +// SPDX-FileCopyrightText: 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; + using coeff_type = gko::solver::detail::coeff_type; + 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(coeff_type{0.9}, coeff_type{1.1}) + .on(exec)), + alpha(0.5), + beta(0.25), + inner_sol(gko::initialize( + {{0.5, 0.125, -0.125}, {0.25, 0.5, -1.0}, {1.5, -0.25, 1.5}}, + this->exec)), + update_sol(gko::initialize( + {{1.0, 0.0, -0.5}, {0.5, 1.0, -1.0}, {-1.5, -0.5, 1.0}}, + this->exec)), + output(gko::initialize( + {{-1.0, 0.5, -0.0}, {0.75, 0.25, -1.25}, {1.0, -1.25, 3.0}}, + this->exec)) + {} + + std::shared_ptr exec; + std::shared_ptr mtx; + std::unique_ptr chebyshev_factory; + coeff_type alpha; + coeff_type beta; + std::shared_ptr inner_sol; + std::shared_ptr update_sol; + std::shared_ptr output; +}; + +TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); + + +TYPED_TEST(Chebyshev, KernelInitUpdate) +{ + using value_type = typename TestFixture::value_type; + using Mtx = typename TestFixture::Mtx; + + gko::kernels::reference::chebyshev::init_update( + this->exec, this->alpha, this->inner_sol.get(), this->update_sol.get(), + this->output.get()); + + GKO_ASSERT_MTX_NEAR(this->update_sol, this->inner_sol, 0); + GKO_ASSERT_MTX_NEAR(this->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; + + gko::kernels::reference::chebyshev::update( + this->exec, this->alpha, this->beta, this->inner_sol.get(), + this->update_sol.get(), this->output.get()); + + GKO_ASSERT_MTX_NEAR(this->update_sol, this->inner_sol, 0); + GKO_ASSERT_MTX_NEAR( + this->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(this->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); +} + + +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; + using coeff_type = typename TestFixture::coeff_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(coeff_type{0.9}, coeff_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; + using coeff_type = typename TestFixture::coeff_type; + auto ref_solver = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(1u)) + .with_foci(coeff_type{0.9}, coeff_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(coeff_type{0.9}, coeff_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 4181d382144..21c0878070f 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..1447ad8cf4a --- /dev/null +++ b/test/solver/chebyshev_kernels.cpp @@ -0,0 +1,185 @@ +// 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; + using coeff_type = gko::solver::detail::coeff_type; + + 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) +{ + coeff_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) +{ + coeff_type alpha(0.5); + coeff_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); +} + + +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 8fc89881faa..d81781a433d 100644 --- a/test/solver/solver.cpp +++ b/test/solver/solver.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -177,6 +178,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)); + } +}; + + template struct CbGmres : SimpleSolverTest> { static constexpr bool will_not_allocate() { return false; } @@ -904,9 +920,9 @@ using SolverTypes = ::testing::Types, Idr<4>,*/ - Ir, CbGmres<2>, CbGmres<10>, Gmres<2>, Gmres<10>, - FGmres<2>, FGmres<10>, Gcr<2>, Gcr<10>, Minres, LowerTrs, - UpperTrs, LowerTrsUnitdiag, UpperTrsUnitdiag + Ir, Chebyshev, CbGmres<2>, CbGmres<10>, Gmres<2>, + Gmres<10>, FGmres<2>, FGmres<10>, Gcr<2>, Gcr<10>, Minres, + 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<>;