From 1b3671b542586e18eb1bc4f91d85af29930cb358 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Wed, 22 Feb 2023 16:42:58 +0100 Subject: [PATCH 01/19] init chebyshev --- core/CMakeLists.txt | 3 +- core/solver/chebyshev.cpp | 358 ++++++++++++++ core/test/solver/CMakeLists.txt | 3 +- core/test/solver/chebyshev.cpp | 489 ++++++++++++++++++++ include/ginkgo/core/solver/chebyshev.hpp | 266 +++++++++++ reference/CMakeLists.txt | 4 +- reference/test/solver/CMakeLists.txt | 3 +- reference/test/solver/chebyshev_kernels.cpp | 350 ++++++++++++++ 8 files changed, 1471 insertions(+), 5 deletions(-) create mode 100644 core/solver/chebyshev.cpp create mode 100644 core/test/solver/chebyshev.cpp create mode 100644 include/ginkgo/core/solver/chebyshev.hpp create mode 100644 reference/test/solver/chebyshev_kernels.cpp 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/solver/chebyshev.cpp b/core/solver/chebyshev.cpp new file mode 100644 index 00000000000..5218e488379 --- /dev/null +++ b/core/solver/chebyshev.cpp @@ -0,0 +1,358 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include +#include +#include +#include + +#include "core/distributed/helpers.hpp" +#include "core/solver/ir_kernels.hpp" +#include "core/solver/solver_base.hpp" +#include "core/solver/solver_boilerplate.hpp" + + +namespace gko { +namespace solver { +namespace chebyshev { +namespace { + + +GKO_REGISTER_OPERATION(initialize, ir::initialize); + + +} // anonymous namespace +} // namespace chebyshev + + +template +void Chebyshev::set_solver(std::shared_ptr new_solver) +{ + auto exec = this->get_executor(); + if (new_solver) { + GKO_ASSERT_EQUAL_DIMENSIONS(new_solver, this); + GKO_ASSERT_IS_SQUARE_MATRIX(new_solver); + if (new_solver->get_executor() != exec) { + new_solver = gko::clone(exec, new_solver); + } + } + solver_ = new_solver; +} + + +template +Chebyshev& Chebyshev::operator=(const Chebyshev& other) +{ + if (&other != this) { + EnableLinOp::operator=(other); + EnableSolverBase::operator=(other); + EnableIterativeBase::operator=(other); + this->parameters_ = other.parameters_; + this->set_solver(other.get_solver()); + } + return *this; +} + + +template +Chebyshev& Chebyshev::operator=(Chebyshev&& other) +{ + if (&other != this) { + EnableLinOp::operator=(std::move(other)); + EnableSolverBase::operator=(std::move(other)); + EnableIterativeBase::operator=(std::move(other)); + this->parameters_ = std::exchange(other.parameters_, parameters_type{}); + this->set_solver(other.get_solver()); + other.set_solver(nullptr); + } + 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_solver( + share(as(this->get_solver())->transpose())) + .with_criteria(this->get_stop_criterion_factory()) + .with_upper_eigval(parameters_.upper_eigval) + .with_lower_eigval(parameters_.lower_eigval) + .on(this->get_executor()) + ->generate( + share(as(this->get_system_matrix())->transpose())); +} + + +template +std::unique_ptr Chebyshev::conj_transpose() const +{ + return build() + .with_generated_solver( + share(as(this->get_solver())->conj_transpose())) + .with_criteria(this->get_stop_criterion_factory()) + .with_upper_eigval(conj(parameters_.upper_eigval)) + .with_lower_eigval(conj(parameters_.lower_eigval)) + .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; + constexpr uint8 relative_stopping_id{1}; + + 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); + // Use the scalar first + auto num_keep = this->get_parameters().num_keep; + auto alpha = this->template create_workspace_scalar( + GKO_SOLVER_TRAITS::alpha, num_keep + 1); + auto beta = this->template create_workspace_scalar( + GKO_SOLVER_TRAITS::beta, num_keep + 1); + + GKO_SOLVER_ONE_MINUS_ONE(); + + auto alpha_ref = ValueType{1} / center_eig_; + auto beta_ref = + ValueType{0.5} * (radius_eig_ * alpha_ref) * (radius_eig_ * alpha_ref); + + bool one_changed{}; + 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; + this->template log( + this, iter, residual_ptr, dense_x); + + if (iter == 0) { + // In iter 0, the iteration and residual are updated. + if (stop_criterion->update() + .num_iterations(iter) + .residual(residual_ptr) + .solution(dense_x) + .check(relative_stopping_id, true, &stop_status, + &one_changed)) { + break; + } + } else { + // In the other iterations, the residual can be updated separately. + if (stop_criterion->update() + .num_iterations(iter) + .solution(dense_x) + .check(relative_stopping_id, false, &stop_status, + &one_changed)) { + break; + } + residual_ptr = residual; + // residual = b - A * x + residual->copy_from(dense_b); + this->get_system_matrix()->apply(lend(neg_one_op), dense_x, + lend(one_op), lend(residual)); + if (stop_criterion->update() + .num_iterations(iter) + .residual(residual_ptr) + .solution(dense_x) + .check(relative_stopping_id, true, &stop_status, + &one_changed)) { + break; + } + } + + if (solver_->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); + } + solver_->apply(residual_ptr, inner_solution); + auto index = (iter >= num_keep) ? num_keep : iter; + auto alpha_scalar = + alpha->create_submatrix(span{0, 1}, span{index, index + 1}); + auto beta_scalar = + beta->create_submatrix(span{0, 1}, span{index, index + 1}); + if (iter == 0) { + if (num_generated_ < num_keep) { + alpha_scalar->fill(alpha_ref); + num_generated_++; + } + // x = x + alpha * inner_solution + dense_x->add_scaled(alpha_scalar.get(), inner_solution); + update_solution->copy_from(inner_solution); + continue; + } + // beta_ref for iter == 1 is initialized in the beginning + if (iter > 1) { + beta_ref = (radius_eig_ * alpha_ref / ValueType{2.0}) * + (radius_eig_ * alpha_ref / ValueType{2.0}); + } + alpha_ref = ValueType{1.0} / (center_eig_ - beta_ref / alpha_ref); + // The last one is always the updated one + if (num_generated_ < num_keep || iter >= num_keep) { + alpha_scalar->fill(alpha_ref); + beta_scalar->fill(beta_ref); + } + if (num_generated_ < num_keep) { + num_generated_++; + } + // z = z + beta * p + // p = z + inner_solution->add_scaled(beta_scalar.get(), update_solution); + update_solution->copy_from(inner_solution); + // x + alpha * p + dense_x->add_scaled(alpha_scalar.get(), update_solution); + } +} + + +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 7; +} + + +template +std::vector workspace_traits>::op_names( + const Solver&) +{ + return { + "residual", "inner_solution", "update_solution", "alpha", "beta", + "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/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..dbc9a65679e --- /dev/null +++ b/core/test/solver/chebyshev.cpp @@ -0,0 +1,489 @@ +// 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" + + +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).on(exec), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(exec)) + .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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// ASSERT_EQ(this->chebyshev_factory->get_executor(), this->exec); +//} + + +TYPED_TEST(Chebyshev, ChebyshevFactoryCreatesCorrectSolver) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using Solver = typename TestFixture::Solver; +// ASSERT_EQ(this->solver->get_size(), gko::dim<2>(3, 3)); +// auto cg_solver = static_cast(this->solver.get()); +// ASSERT_NE(cg_solver->get_system_matrix(), nullptr); +// ASSERT_EQ(cg_solver->get_system_matrix(), this->mtx); +//} + + +TYPED_TEST(Chebyshev, CanBeCopied) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using Mtx = typename TestFixture::Mtx; +// using Solver = typename TestFixture::Solver; +// auto copy = this->chebyshev_factory->generate(Mtx::create(this->exec)); +// +// copy->copy_from(std::move(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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// ASSERT_TRUE(this->solver->apply_uses_initial_guess()); +//} + + +TYPED_TEST(Chebyshev, CanSetInnerSolverInFactory) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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).on(this->exec), +// gko::stop::ResidualNorm::build() +// .with_reduction_factor(r::value) +// .on(this->exec)) +// .with_solver( +// Solver::build() +// .with_criteria( +// gko::stop::Iteration::build().with_max_iters(3u).on( +// this->exec)) +// .on(this->exec)) +// .on(this->exec); +// auto solver = chebyshev_factory->generate(this->mtx); +// auto inner_solver = dynamic_cast( +// static_cast(solver.get())->get_solver().get()); +// +// ASSERT_NE(inner_solver, nullptr); +// ASSERT_EQ(inner_solver->get_size(), gko::dim<2>(3, 3)); +// ASSERT_EQ(inner_solver->get_system_matrix(), this->mtx); +//} + + +TYPED_TEST(Chebyshev, CanSetGeneratedInnerSolverInFactory) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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)) +// .on(this->exec) +// ->generate(this->mtx); +// +// auto chebyshev_factory = +// Solver::build() +// .with_criteria( +// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) +// .with_generated_solver(chebyshev_solver) +// .on(this->exec); +// auto solver = chebyshev_factory->generate(this->mtx); +// auto inner_solver = solver->get_solver(); +// +// ASSERT_NE(inner_solver.get(), nullptr); +// ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); +//} + + +TYPED_TEST(Chebyshev, CanSetCriteriaAgain) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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)) +// .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)) +// .with_generated_solver(chebyshev_solver) +// .on(this->exec); +// +// ASSERT_THROW(chebyshev_factory->generate(this->mtx), +// gko::DimensionMismatch); +//} + + +TYPED_TEST(Chebyshev, CanSetInnerSolver) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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)) +// .on(this->exec) +// ->generate(this->mtx); +// +// auto chebyshev_factory = +// Solver::build() +// .with_criteria( +// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) +// .on(this->exec); +// auto solver = chebyshev_factory->generate(this->mtx); +// solver->set_solver(chebyshev_solver); +// auto inner_solver = solver->get_solver(); +// +// ASSERT_NE(inner_solver.get(), nullptr); +// ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); +//} + + +TYPED_TEST(Chebyshev, CanSetApplyWithInitialGuessMode) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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).on( +// this->exec)) +// .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) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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)) +// .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)) +// .on(this->exec); +// auto solver = chebyshev_factory->generate(this->mtx); +// +// ASSERT_THROW(solver->set_solver(chebyshev_solver), +// gko::DimensionMismatch); +//} + + +TYPED_TEST(Chebyshev, ThrowsOnRectangularMatrixInFactory) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// 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); +//} + + +TYPED_TEST(Chebyshev, DefaultRelaxationFactor) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using value_type = typename TestFixture::value_type; +// const value_type relaxation_factor{0.5}; +// +// auto richardson = +// gko::solver::Richardson::build() +// .with_criteria( +// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec), +// gko::stop::ResidualNorm::build() +// .with_reduction_factor(r::value) +// .on(this->exec)) +// .on(this->exec) +// ->generate(this->mtx); +// +// ASSERT_EQ(richardson->get_parameters().relaxation_factor, value_type{1}); +//} + + +TYPED_TEST(Chebyshev, UseAsRichardson) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using value_type = typename TestFixture::value_type; +// const value_type relaxation_factor{0.5}; +// +// auto richardson = +// gko::solver::Richardson::build() +// .with_criteria( +// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec), +// gko::stop::ResidualNorm::build() +// .with_reduction_factor(r::value) +// .on(this->exec)) +// .with_relaxation_factor(relaxation_factor) +// .on(this->exec) +// ->generate(this->mtx); +// +// ASSERT_EQ(richardson->get_parameters().relaxation_factor, +// value_type{0.5}); +//} + + +TYPED_TEST(Chebyshev, DefaultSmootherBuildWithSolver) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using value_type = typename TestFixture::value_type; +// using Solver = typename TestFixture::Solver; +// auto solver = gko::as(share(std::move(this->solver))); +// +// auto smoother_factory = gko::solver::build_smoother(solver); +// auto criteria = +// std::dynamic_pointer_cast( +// smoother_factory->get_parameters().criteria.at(0)); +// +// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, +// value_type{0.9}); +// ASSERT_NE(criteria.get(), nullptr); +// ASSERT_EQ(criteria->get_parameters().max_iters, 1); +// ASSERT_EQ(smoother_factory->get_parameters().generated_solver.get(), +// solver.get()); +//} + + +TYPED_TEST(Chebyshev, DefaultSmootherBuildWithFactory) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using value_type = typename TestFixture::value_type; +// using Solver = typename TestFixture::Solver; +// auto factory = this->chebyshev_factory; +// +// auto smoother_factory = gko::solver::build_smoother(factory); +// auto criteria = +// std::dynamic_pointer_cast( +// smoother_factory->get_parameters().criteria.at(0)); +// +// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, +// value_type{0.9}); +// ASSERT_NE(criteria.get(), nullptr); +// ASSERT_EQ(criteria->get_parameters().max_iters, 1); +// ASSERT_EQ(smoother_factory->get_parameters().solver.get(), factory.get()); +//} + + +TYPED_TEST(Chebyshev, SmootherBuildWithSolver) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using value_type = typename TestFixture::value_type; +// using Solver = typename TestFixture::Solver; +// auto solver = gko::as(gko::share(std::move(this->solver))); +// +// auto smoother_factory = +// gko::solver::build_smoother(solver, 3, value_type{0.5}); +// auto criteria = +// std::dynamic_pointer_cast( +// smoother_factory->get_parameters().criteria.at(0)); +// +// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, +// value_type{0.5}); +// ASSERT_NE(criteria.get(), nullptr); +// ASSERT_EQ(criteria->get_parameters().max_iters, 3); +// ASSERT_EQ(smoother_factory->get_parameters().generated_solver.get(), +// solver.get()); +//} + + +TYPED_TEST(Chebyshev, SmootherBuildWithFactory) +GKO_NOT_IMPLEMENTED; +//{ +// TODO (script:chebyshev): change the code imported from solver/ir if needed +// using value_type = typename TestFixture::value_type; +// using Solver = typename TestFixture::Solver; +// auto factory = this->chebyshev_factory; +// +// auto smoother_factory = +// gko::solver::build_smoother(factory, 3, value_type{0.5}); +// auto criteria = +// std::dynamic_pointer_cast( +// smoother_factory->get_parameters().criteria.at(0)); +// +// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, +// value_type{0.5}); +// ASSERT_NE(criteria.get(), nullptr); +// ASSERT_EQ(criteria->get_parameters().max_iters, 3); +// ASSERT_EQ(smoother_factory->get_parameters().solver.get(), factory.get()); +//} + + +} // namespace diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp new file mode 100644 index 00000000000..fc9ba88339c --- /dev/null +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -0,0 +1,266 @@ +// 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 + + +namespace gko { +namespace solver { + + +/** + * Chebyshev iteration is an iterative method that uses another coarse + * method to approximate the error of the current solution via the current + * residual. It has another term for the difference of solution. Moreover, this + * method requires knowledge about the spectrum of the matrix. + * + * ``` + * solution = initial_guess + * while not converged: + * residual = b - A solution + * error = solver(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 : public EnableLinOp>, + public EnableSolverBase>, + public EnableIterativeBase>, + 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; + } + + /** + * Returns the solver operator used as the inner solver. + * + * @return the solver operator used as the inner solver + */ + std::shared_ptr get_solver() const { return solver_; } + + /** + * Sets the solver operator used as the inner solver. + * + * @param new_solver the new inner solver + */ + void set_solver(std::shared_ptr new_solver); + + /** + * 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&&); + + GKO_CREATE_FACTORY_PARAMETERS(parameters, Factory) + { + /** + * Criterion factories. + */ + std::vector> + GKO_FACTORY_PARAMETER_VECTOR(criteria, nullptr); + + /** + * Inner solver factory. + */ + std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( + solver, nullptr); + + /** + * Already generated solver. If one is provided, the factory `solver` + * will be ignored. + */ + std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( + generated_solver, nullptr); + + /** + * Upper estimated eigenvalue + */ + ValueType GKO_FACTORY_PARAMETER_SCALAR(upper_eigval, value_type{1}); + + /** + * Lower estimated eigenvalue + */ + ValueType GKO_FACTORY_PARAMETER_SCALAR(lower_eigval, value_type{0}); + + /** + * 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); + + /** + * The number of scalar to keep + */ + int GKO_FACTORY_PARAMETER_SCALAR(num_keep, 2); + }; + GKO_ENABLE_LIN_OP_FACTORY(Chebyshev, parameters, Factory); + GKO_ENABLE_BUILD_METHOD(Factory); + +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; + + void set_relaxation_factor( + std::shared_ptr> new_factor); + + explicit Chebyshev(std::shared_ptr exec) + : EnableLinOp(std::move(exec)) + {} + + explicit Chebyshev(const Factory* factory, + std::shared_ptr system_matrix) + : EnableLinOp(factory->get_executor(), + gko::transpose(system_matrix->get_size())), + EnableSolverBase{std::move(system_matrix)}, + EnableIterativeBase{ + stop::combine(factory->get_parameters().criteria)}, + parameters_{factory->get_parameters()} + { + if (parameters_.generated_solver) { + this->set_solver(parameters_.generated_solver); + } else if (parameters_.solver) { + this->set_solver( + parameters_.solver->generate(this->get_system_matrix())); + } else { + this->set_solver(matrix::Identity::create( + this->get_executor(), this->get_size())); + } + this->set_default_initial_guess(parameters_.default_initial_guess); + center_eig_ = (parameters_.upper_eigval + parameters_.lower_eigval) / + ValueType{2}; + radius_eig_ = (parameters_.upper_eigval - parameters_.lower_eigval) / + ValueType{2}; + // if changing the lower/upper eig, need to reset it to zero + num_generated_ = 0; + } + +private: + std::shared_ptr solver_{}; + mutable int num_generated_; + ValueType center_eig_; + ValueType radius_eig_; +}; + + +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; + // alpha + constexpr static int alpha = 3; + // beta + constexpr static int beta = 4; + // constant 1.0 scalar + constexpr static int one = 5; + // constant -1.0 scalar + constexpr static int minus_one = 6; + + // 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..82d99105a6e 100644 --- a/reference/CMakeLists.txt +++ b/reference/CMakeLists.txt @@ -55,13 +55,13 @@ 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/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/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..ec1b9d4211c --- /dev/null +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -0,0 +1,350 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#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( + {{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).on( + exec), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(exec)) + .with_upper_eigval(value_type{1.1}) + .with_lower_eigval(value_type{0.9}) + .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, 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 value_type = gko::next_precision; + using Mtx = 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 value_type = + gko::to_complex>; + using Mtx = gko::matrix::Dense; + 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_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 inner_solver_factory = gko::share( + gko::solver::Gmres::build() + .with_criteria(gko::stop::ResidualNorm::build() + .with_reduction_factor(inner_reduction_factor) + .on(this->exec)) + .on(this->exec)); + + auto solver_factory = + gko::solver::Chebyshev::build() + .with_criteria(gko::stop::Iteration::build().with_max_iters(30u).on( + this->exec), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(this->exec)) + .with_solver(inner_solver_factory) + .with_upper_eigval(value_type{1.1}) + .with_lower_eigval(value_type{0.9}) + .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 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, 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 Scalar = gko::matrix::Dense< + gko::next_precision>; + 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::value * 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).on(this->exec)) + .with_upper_eigval(value_type{1.1}) + .with_lower_eigval(value_type{0.9}) + .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).on( + this->exec)) + .with_upper_eigval(value_type{1.1}) + .with_lower_eigval(value_type{0.9}) + .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(lend(b), lend(x)); + ref_solver->apply(lend(b), lend(ref_x)); + + GKO_ASSERT_MTX_NEAR(x, ref_x, 0.0); + } +} + + +} // namespace From 53bd8acf4d8150fe0eff72c3d32f3024680fe903 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 21 Mar 2023 13:34:02 +0100 Subject: [PATCH 02/19] add the check for kept data and remove lend --- core/solver/chebyshev.cpp | 8 ++-- reference/test/solver/chebyshev_kernels.cpp | 46 ++++++++++++++++++++- 2 files changed, 49 insertions(+), 5 deletions(-) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 5218e488379..81e068bd2a6 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -213,8 +213,8 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, residual_ptr = residual; // residual = b - A * x residual->copy_from(dense_b); - this->get_system_matrix()->apply(lend(neg_one_op), dense_x, - lend(one_op), lend(residual)); + this->get_system_matrix()->apply(neg_one_op, dense_x, one_op, + residual); if (stop_criterion->update() .num_iterations(iter) .residual(residual_ptr) @@ -232,7 +232,7 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, inner_solution->copy_from(residual_ptr); } solver_->apply(residual_ptr, inner_solution); - auto index = (iter >= num_keep) ? num_keep : iter; + size_type index = (iter >= num_keep) ? num_keep : iter; auto alpha_scalar = alpha->create_submatrix(span{0, 1}, span{index, index + 1}); auto beta_scalar = @@ -240,6 +240,8 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, if (iter == 0) { if (num_generated_ < num_keep) { alpha_scalar->fill(alpha_ref); + // unused beta for first iteration, but fill zero + beta_scalar->fill(zero()); num_generated_++; } // x = x + alpha * inner_solution diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index ec1b9d4211c..e1224a29885 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -50,6 +50,48 @@ class Chebyshev : public ::testing::Test { TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); +TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + auto upper = value_type{1.1}; + auto lower = value_type{0.9}; + auto factory = + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) + .with_upper_eigval(upper) + .with_lower_eigval(lower) + .with_num_keep(3) + .on(this->exec); + auto solver = 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()); + + auto alpha = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::alpha)); + auto beta = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::beta)); + ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); + ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 4})); + // check the num_keep alpha, beta + auto d = (upper + lower) / value_type{2}; + auto c = (upper - lower) / value_type{2}; + EXPECT_EQ(alpha->at(0, 0), value_type{1} / d); + EXPECT_EQ(beta->at(0, 0), value_type{0}); + EXPECT_EQ(beta->at(0, 1), + value_type{0.5} * (c * alpha->at(0, 0)) * (c * alpha->at(0, 0))); + EXPECT_EQ(alpha->at(0, 1), + value_type{1} / (d - beta->at(0, 1) / alpha->at(0, 0))); + EXPECT_EQ(beta->at(0, 2), (c * alpha->at(0, 1) / value_type{2}) * + (c * alpha->at(0, 1) / value_type{2})); + EXPECT_EQ(alpha->at(0, 2), + value_type{1} / (d - beta->at(0, 2) / alpha->at(0, 1))); +} + TYPED_TEST(Chebyshev, SolvesTriangularSystem) { @@ -339,8 +381,8 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) } else { ref_x = gko::initialize({0.0, 0.0, 0.0}, this->exec); } - solver->apply(lend(b), lend(x)); - ref_solver->apply(lend(b), lend(ref_x)); + solver->apply(b, x); + ref_solver->apply(b, ref_x); GKO_ASSERT_MTX_NEAR(x, ref_x, 0.0); } From 9cbb9ba5504ce38724d8b0f3519eaa348d107eb5 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 21 Mar 2023 16:02:33 +0100 Subject: [PATCH 03/19] core test --- core/test/solver/chebyshev.cpp | 614 +++++++++++++-------------------- 1 file changed, 234 insertions(+), 380 deletions(-) diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp index dbc9a65679e..f27d6adc283 100644 --- a/core/test/solver/chebyshev.cpp +++ b/core/test/solver/chebyshev.cpp @@ -62,428 +62,282 @@ TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); TYPED_TEST(Chebyshev, ChebyshevFactoryKnowsItsExecutor) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// ASSERT_EQ(this->chebyshev_factory->get_executor(), this->exec); -//} +{ + ASSERT_EQ(this->chebyshev_factory->get_executor(), this->exec); +} TYPED_TEST(Chebyshev, ChebyshevFactoryCreatesCorrectSolver) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using Solver = typename TestFixture::Solver; -// ASSERT_EQ(this->solver->get_size(), gko::dim<2>(3, 3)); -// auto cg_solver = static_cast(this->solver.get()); -// ASSERT_NE(cg_solver->get_system_matrix(), nullptr); -// ASSERT_EQ(cg_solver->get_system_matrix(), this->mtx); -//} +{ + using Solver = typename TestFixture::Solver; + ASSERT_EQ(this->solver->get_size(), gko::dim<2>(3, 3)); + auto cg_solver = static_cast(this->solver.get()); + ASSERT_NE(cg_solver->get_system_matrix(), nullptr); + ASSERT_EQ(cg_solver->get_system_matrix(), this->mtx); +} TYPED_TEST(Chebyshev, CanBeCopied) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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()); -//} +{ + 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) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using Mtx = typename TestFixture::Mtx; -// using Solver = typename TestFixture::Solver; -// auto copy = this->chebyshev_factory->generate(Mtx::create(this->exec)); -// -// copy->copy_from(std::move(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()); -//} +{ + 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) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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()); -//} +{ + 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) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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); -//} +{ + 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) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// ASSERT_TRUE(this->solver->apply_uses_initial_guess()); -//} +{ + 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).on(this->exec)) + .with_upper_eigval(value_type{1.2}) + .with_lower_eigval(value_type{0.2}) + .on(this->exec) + ->generate(this->mtx); + + ASSERT_EQ(chebyshev_solver->get_parameters().lower_eigval, value_type{0.2}); + ASSERT_EQ(chebyshev_solver->get_parameters().upper_eigval, value_type{1.2}); +} TYPED_TEST(Chebyshev, CanSetInnerSolverInFactory) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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).on(this->exec), -// gko::stop::ResidualNorm::build() -// .with_reduction_factor(r::value) -// .on(this->exec)) -// .with_solver( -// Solver::build() -// .with_criteria( -// gko::stop::Iteration::build().with_max_iters(3u).on( -// this->exec)) -// .on(this->exec)) -// .on(this->exec); -// auto solver = chebyshev_factory->generate(this->mtx); -// auto inner_solver = dynamic_cast( -// static_cast(solver.get())->get_solver().get()); -// -// ASSERT_NE(inner_solver, nullptr); -// ASSERT_EQ(inner_solver->get_size(), gko::dim<2>(3, 3)); -// ASSERT_EQ(inner_solver->get_system_matrix(), this->mtx); -//} +{ + 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).on(this->exec), + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(this->exec)) + .with_solver( + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(3u).on( + this->exec)) + .on(this->exec)) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + auto inner_solver = dynamic_cast( + static_cast(solver.get())->get_solver().get()); + + ASSERT_NE(inner_solver, nullptr); + ASSERT_EQ(inner_solver->get_size(), gko::dim<2>(3, 3)); + ASSERT_EQ(inner_solver->get_system_matrix(), this->mtx); +} TYPED_TEST(Chebyshev, CanSetGeneratedInnerSolverInFactory) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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)) -// .on(this->exec) -// ->generate(this->mtx); -// -// auto chebyshev_factory = -// Solver::build() -// .with_criteria( -// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) -// .with_generated_solver(chebyshev_solver) -// .on(this->exec); -// auto solver = chebyshev_factory->generate(this->mtx); -// auto inner_solver = solver->get_solver(); -// -// ASSERT_NE(inner_solver.get(), nullptr); -// ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); -//} +{ + 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)) + .on(this->exec) + ->generate(this->mtx); + + auto chebyshev_factory = + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) + .with_generated_solver(chebyshev_solver) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + auto inner_solver = solver->get_solver(); + + ASSERT_NE(inner_solver.get(), nullptr); + ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); +} TYPED_TEST(Chebyshev, CanSetCriteriaAgain) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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); -//} +{ + 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) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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)) -// .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)) -// .with_generated_solver(chebyshev_solver) -// .on(this->exec); -// -// ASSERT_THROW(chebyshev_factory->generate(this->mtx), -// gko::DimensionMismatch); -//} +{ + 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)) + .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)) + .with_generated_solver(chebyshev_solver) + .on(this->exec); + + ASSERT_THROW(chebyshev_factory->generate(this->mtx), + gko::DimensionMismatch); +} TYPED_TEST(Chebyshev, CanSetInnerSolver) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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)) -// .on(this->exec) -// ->generate(this->mtx); -// -// auto chebyshev_factory = -// Solver::build() -// .with_criteria( -// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) -// .on(this->exec); -// auto solver = chebyshev_factory->generate(this->mtx); -// solver->set_solver(chebyshev_solver); -// auto inner_solver = solver->get_solver(); -// -// ASSERT_NE(inner_solver.get(), nullptr); -// ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); -//} +{ + 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)) + .on(this->exec) + ->generate(this->mtx); + + auto chebyshev_factory = + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + solver->set_solver(chebyshev_solver); + auto inner_solver = solver->get_solver(); + + ASSERT_NE(inner_solver.get(), nullptr); + ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); +} TYPED_TEST(Chebyshev, CanSetApplyWithInitialGuessMode) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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).on( -// this->exec)) -// .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); -// } -//} +{ + 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).on( + this->exec)) + .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) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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)) -// .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)) -// .on(this->exec); -// auto solver = chebyshev_factory->generate(this->mtx); -// -// ASSERT_THROW(solver->set_solver(chebyshev_solver), -// gko::DimensionMismatch); -//} +{ + 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)) + .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)) + .on(this->exec); + auto solver = chebyshev_factory->generate(this->mtx); + + ASSERT_THROW(solver->set_solver(chebyshev_solver), gko::DimensionMismatch); +} TYPED_TEST(Chebyshev, ThrowsOnRectangularMatrixInFactory) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// 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); -//} - - -TYPED_TEST(Chebyshev, DefaultRelaxationFactor) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using value_type = typename TestFixture::value_type; -// const value_type relaxation_factor{0.5}; -// -// auto richardson = -// gko::solver::Richardson::build() -// .with_criteria( -// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec), -// gko::stop::ResidualNorm::build() -// .with_reduction_factor(r::value) -// .on(this->exec)) -// .on(this->exec) -// ->generate(this->mtx); -// -// ASSERT_EQ(richardson->get_parameters().relaxation_factor, value_type{1}); -//} - - -TYPED_TEST(Chebyshev, UseAsRichardson) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using value_type = typename TestFixture::value_type; -// const value_type relaxation_factor{0.5}; -// -// auto richardson = -// gko::solver::Richardson::build() -// .with_criteria( -// gko::stop::Iteration::build().with_max_iters(3u).on(this->exec), -// gko::stop::ResidualNorm::build() -// .with_reduction_factor(r::value) -// .on(this->exec)) -// .with_relaxation_factor(relaxation_factor) -// .on(this->exec) -// ->generate(this->mtx); -// -// ASSERT_EQ(richardson->get_parameters().relaxation_factor, -// value_type{0.5}); -//} - - -TYPED_TEST(Chebyshev, DefaultSmootherBuildWithSolver) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using value_type = typename TestFixture::value_type; -// using Solver = typename TestFixture::Solver; -// auto solver = gko::as(share(std::move(this->solver))); -// -// auto smoother_factory = gko::solver::build_smoother(solver); -// auto criteria = -// std::dynamic_pointer_cast( -// smoother_factory->get_parameters().criteria.at(0)); -// -// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, -// value_type{0.9}); -// ASSERT_NE(criteria.get(), nullptr); -// ASSERT_EQ(criteria->get_parameters().max_iters, 1); -// ASSERT_EQ(smoother_factory->get_parameters().generated_solver.get(), -// solver.get()); -//} - - -TYPED_TEST(Chebyshev, DefaultSmootherBuildWithFactory) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using value_type = typename TestFixture::value_type; -// using Solver = typename TestFixture::Solver; -// auto factory = this->chebyshev_factory; -// -// auto smoother_factory = gko::solver::build_smoother(factory); -// auto criteria = -// std::dynamic_pointer_cast( -// smoother_factory->get_parameters().criteria.at(0)); -// -// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, -// value_type{0.9}); -// ASSERT_NE(criteria.get(), nullptr); -// ASSERT_EQ(criteria->get_parameters().max_iters, 1); -// ASSERT_EQ(smoother_factory->get_parameters().solver.get(), factory.get()); -//} - - -TYPED_TEST(Chebyshev, SmootherBuildWithSolver) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using value_type = typename TestFixture::value_type; -// using Solver = typename TestFixture::Solver; -// auto solver = gko::as(gko::share(std::move(this->solver))); -// -// auto smoother_factory = -// gko::solver::build_smoother(solver, 3, value_type{0.5}); -// auto criteria = -// std::dynamic_pointer_cast( -// smoother_factory->get_parameters().criteria.at(0)); -// -// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, -// value_type{0.5}); -// ASSERT_NE(criteria.get(), nullptr); -// ASSERT_EQ(criteria->get_parameters().max_iters, 3); -// ASSERT_EQ(smoother_factory->get_parameters().generated_solver.get(), -// solver.get()); -//} - - -TYPED_TEST(Chebyshev, SmootherBuildWithFactory) -GKO_NOT_IMPLEMENTED; -//{ -// TODO (script:chebyshev): change the code imported from solver/ir if needed -// using value_type = typename TestFixture::value_type; -// using Solver = typename TestFixture::Solver; -// auto factory = this->chebyshev_factory; -// -// auto smoother_factory = -// gko::solver::build_smoother(factory, 3, value_type{0.5}); -// auto criteria = -// std::dynamic_pointer_cast( -// smoother_factory->get_parameters().criteria.at(0)); -// -// ASSERT_EQ(smoother_factory->get_parameters().relaxation_factor, -// value_type{0.5}); -// ASSERT_NE(criteria.get(), nullptr); -// ASSERT_EQ(criteria->get_parameters().max_iters, 3); -// ASSERT_EQ(smoother_factory->get_parameters().solver.get(), factory.get()); -//} +{ + 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 From 9bc798c6faa1d38cab36248c46a653f98c2e0b90 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 21 Mar 2023 16:03:09 +0100 Subject: [PATCH 04/19] change to foci --- core/solver/chebyshev.cpp | 19 +++++++-------- core/test/solver/chebyshev.cpp | 7 +++--- include/ginkgo/core/solver/chebyshev.hpp | 26 ++++++++++----------- reference/test/solver/chebyshev_kernels.cpp | 15 ++++-------- 4 files changed, 30 insertions(+), 37 deletions(-) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 81e068bd2a6..e888a035d65 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -93,8 +93,7 @@ std::unique_ptr Chebyshev::transpose() const .with_generated_solver( share(as(this->get_solver())->transpose())) .with_criteria(this->get_stop_criterion_factory()) - .with_upper_eigval(parameters_.upper_eigval) - .with_lower_eigval(parameters_.lower_eigval) + .with_foci(parameters_.foci) .on(this->get_executor()) ->generate( share(as(this->get_system_matrix())->transpose())); @@ -108,8 +107,8 @@ std::unique_ptr Chebyshev::conj_transpose() const .with_generated_solver( share(as(this->get_solver())->conj_transpose())) .with_criteria(this->get_stop_criterion_factory()) - .with_upper_eigval(conj(parameters_.upper_eigval)) - .with_lower_eigval(conj(parameters_.lower_eigval)) + .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())); @@ -164,9 +163,9 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, GKO_SOLVER_ONE_MINUS_ONE(); - auto alpha_ref = ValueType{1} / center_eig_; - auto beta_ref = - ValueType{0.5} * (radius_eig_ * alpha_ref) * (radius_eig_ * alpha_ref); + auto alpha_ref = ValueType{1} / center_; + auto beta_ref = ValueType{0.5} * (foci_direction_ * alpha_ref) * + (foci_direction_ * alpha_ref); bool one_changed{}; auto& stop_status = this->template create_workspace_array( @@ -251,10 +250,10 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, } // beta_ref for iter == 1 is initialized in the beginning if (iter > 1) { - beta_ref = (radius_eig_ * alpha_ref / ValueType{2.0}) * - (radius_eig_ * alpha_ref / ValueType{2.0}); + beta_ref = (foci_direction_ * alpha_ref / ValueType{2.0}) * + (foci_direction_ * alpha_ref / ValueType{2.0}); } - alpha_ref = ValueType{1.0} / (center_eig_ - beta_ref / alpha_ref); + alpha_ref = ValueType{1.0} / (center_ - beta_ref / alpha_ref); // The last one is always the updated one if (num_generated_ < num_keep || iter >= num_keep) { alpha_scalar->fill(alpha_ref); diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp index f27d6adc283..40b495390ec 100644 --- a/core/test/solver/chebyshev.cpp +++ b/core/test/solver/chebyshev.cpp @@ -146,13 +146,12 @@ TYPED_TEST(Chebyshev, CanSetEigenRegion) Solver::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) - .with_upper_eigval(value_type{1.2}) - .with_lower_eigval(value_type{0.2}) + .with_foci(value_type{0.2}, value_type{1.2}) .on(this->exec) ->generate(this->mtx); - ASSERT_EQ(chebyshev_solver->get_parameters().lower_eigval, value_type{0.2}); - ASSERT_EQ(chebyshev_solver->get_parameters().upper_eigval, value_type{1.2}); + ASSERT_EQ(chebyshev_solver->get_parameters().foci, + std::make_pair(value_type{0.2}, value_type{1.2})); } diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index fc9ba88339c..23625039db9 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -140,14 +140,11 @@ class Chebyshev : public EnableLinOp>, generated_solver, nullptr); /** - * Upper estimated eigenvalue + * The pair of foci of ellipse. It is usually be {lower bound of eigval, + * upper bound of eigval} for real matrices. */ - ValueType GKO_FACTORY_PARAMETER_SCALAR(upper_eigval, value_type{1}); - - /** - * Lower estimated eigenvalue - */ - ValueType GKO_FACTORY_PARAMETER_SCALAR(lower_eigval, value_type{0}); + std::pair GKO_FACTORY_PARAMETER_VECTOR( + foci, value_type{0}, value_type{1}); /** * Default initial guess mode. The available options are under @@ -207,10 +204,13 @@ class Chebyshev : public EnableLinOp>, this->get_executor(), this->get_size())); } this->set_default_initial_guess(parameters_.default_initial_guess); - center_eig_ = (parameters_.upper_eigval + parameters_.lower_eigval) / - ValueType{2}; - radius_eig_ = (parameters_.upper_eigval - parameters_.lower_eigval) / - ValueType{2}; + center_ = + (std::get<0>(parameters_.foci) + std::get<1>(parameters_.foci)) / + ValueType{2}; + // the absolute value of foci_direction is the focal direction + foci_direction_ = + (std::get<1>(parameters_.foci) - std::get<0>(parameters_.foci)) / + ValueType{2}; // if changing the lower/upper eig, need to reset it to zero num_generated_ = 0; } @@ -218,8 +218,8 @@ class Chebyshev : public EnableLinOp>, private: std::shared_ptr solver_{}; mutable int num_generated_; - ValueType center_eig_; - ValueType radius_eig_; + ValueType center_; + ValueType foci_direction_; }; diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index e1224a29885..1be6a85c318 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -38,8 +38,7 @@ class Chebyshev : public ::testing::Test { gko::stop::ResidualNorm::build() .with_reduction_factor(r::value) .on(exec)) - .with_upper_eigval(value_type{1.1}) - .with_lower_eigval(value_type{0.9}) + .with_foci(value_type{0.9}, value_type{1.1}) .on(exec)) {} @@ -61,8 +60,7 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) Solver::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) - .with_upper_eigval(upper) - .with_lower_eigval(lower) + .with_foci(lower, upper) .with_num_keep(3) .on(this->exec); auto solver = factory->generate(this->mtx); @@ -186,8 +184,7 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemWithIterativeInnerSolver) .with_reduction_factor(r::value) .on(this->exec)) .with_solver(inner_solver_factory) - .with_upper_eigval(value_type{1.1}) - .with_lower_eigval(value_type{0.9}) + .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); @@ -355,8 +352,7 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) gko::solver::Chebyshev::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(1u).on(this->exec)) - .with_upper_eigval(value_type{1.1}) - .with_lower_eigval(value_type{0.9}) + .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); @@ -367,8 +363,7 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) .with_criteria( gko::stop::Iteration::build().with_max_iters(1u).on( this->exec)) - .with_upper_eigval(value_type{1.1}) - .with_lower_eigval(value_type{0.9}) + .with_foci(value_type{0.9}, value_type{1.1}) .with_default_initial_guess(guess) .on(this->exec) ->generate(this->mtx); From cdabda40e98e0ee5e74f035d7d8c1d497310d5e1 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Thu, 3 Aug 2023 16:02:36 +0200 Subject: [PATCH 05/19] use iteration from stop criterion and update doc --- core/solver/chebyshev.cpp | 119 ++++++++++++---- include/ginkgo/core/solver/chebyshev.hpp | 52 ++----- reference/test/solver/chebyshev_kernels.cpp | 145 +++++++++++++++++++- 3 files changed, 246 insertions(+), 70 deletions(-) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index e888a035d65..cb5f72b969e 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -26,6 +26,37 @@ GKO_REGISTER_OPERATION(initialize, ir::initialize); } // namespace chebyshev +template +Chebyshev::Chebyshev(const Factory* factory, + std::shared_ptr system_matrix) + : EnableLinOp(factory->get_executor(), + gko::transpose(system_matrix->get_size())), + EnableSolverBase{std::move(system_matrix)}, + EnableIterativeBase{ + stop::combine(factory->get_parameters().criteria)}, + parameters_{factory->get_parameters()} +{ + if (parameters_.generated_solver) { + this->set_solver(parameters_.generated_solver); + } else if (parameters_.solver) { + this->set_solver( + parameters_.solver->generate(this->get_system_matrix())); + } else { + this->set_solver(matrix::Identity::create( + this->get_executor(), this->get_size())); + } + 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}; + // if changing the lower/upper eig, need to reset it to zero + num_generated_scalar_ = 0; + num_max_generation_ = 3; +} + + template void Chebyshev::set_solver(std::shared_ptr new_solver) { @@ -154,12 +185,29 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, GKO_SOLVER_VECTOR(residual, dense_b); GKO_SOLVER_VECTOR(inner_solution, dense_b); GKO_SOLVER_VECTOR(update_solution, dense_b); + // Use the scalar first - auto num_keep = this->get_parameters().num_keep; + // get the iteration information from stopping criterion. + if (auto combined = + std::dynamic_pointer_cast( + this->get_stop_criterion_factory())) { + for (const auto& factory : combined->get_parameters().criteria) { + if (auto iter_stop = std::dynamic_pointer_cast< + const gko::stop::Iteration::Factory>(factory)) { + num_max_generation_ = std::max( + num_max_generation_, iter_stop->get_parameters().max_iters); + } + } + } else if (auto iter_stop = std::dynamic_pointer_cast< + const gko::stop::Iteration::Factory>( + this->get_stop_criterion_factory())) { + num_max_generation_ = std::max(num_max_generation_, + iter_stop->get_parameters().max_iters); + } auto alpha = this->template create_workspace_scalar( - GKO_SOLVER_TRAITS::alpha, num_keep + 1); + GKO_SOLVER_TRAITS::alpha, num_max_generation_ + 1); auto beta = this->template create_workspace_scalar( - GKO_SOLVER_TRAITS::beta, num_keep + 1); + GKO_SOLVER_TRAITS::beta, num_max_generation_ + 1); GKO_SOLVER_ONE_MINUS_ONE(); @@ -187,26 +235,33 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, int iter = -1; while (true) { ++iter; - this->template log( - this, iter, residual_ptr, dense_x); - if (iter == 0) { // In iter 0, the iteration and residual are updated. - if (stop_criterion->update() - .num_iterations(iter) - .residual(residual_ptr) - .solution(dense_x) - .check(relative_stopping_id, true, &stop_status, - &one_changed)) { + bool 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; } } else { // In the other iterations, the residual can be updated separately. - if (stop_criterion->update() - .num_iterations(iter) - .solution(dense_x) - .check(relative_stopping_id, false, &stop_status, - &one_changed)) { + 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; @@ -214,12 +269,16 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, residual->copy_from(dense_b); this->get_system_matrix()->apply(neg_one_op, dense_x, one_op, residual); - if (stop_criterion->update() - .num_iterations(iter) - .residual(residual_ptr) - .solution(dense_x) - .check(relative_stopping_id, true, &stop_status, - &one_changed)) { + 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; } } @@ -231,17 +290,18 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, inner_solution->copy_from(residual_ptr); } solver_->apply(residual_ptr, inner_solution); - size_type index = (iter >= num_keep) ? num_keep : iter; + size_type index = + (iter >= num_max_generation_) ? num_max_generation_ : iter; auto alpha_scalar = alpha->create_submatrix(span{0, 1}, span{index, index + 1}); auto beta_scalar = beta->create_submatrix(span{0, 1}, span{index, index + 1}); if (iter == 0) { - if (num_generated_ < num_keep) { + if (num_generated_scalar_ < num_max_generation_) { alpha_scalar->fill(alpha_ref); // unused beta for first iteration, but fill zero beta_scalar->fill(zero()); - num_generated_++; + num_generated_scalar_++; } // x = x + alpha * inner_solution dense_x->add_scaled(alpha_scalar.get(), inner_solution); @@ -255,12 +315,13 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, } alpha_ref = ValueType{1.0} / (center_ - beta_ref / alpha_ref); // The last one is always the updated one - if (num_generated_ < num_keep || iter >= num_keep) { + if (num_generated_scalar_ < num_max_generation_ || + iter >= num_max_generation_) { alpha_scalar->fill(alpha_ref); beta_scalar->fill(beta_ref); } - if (num_generated_ < num_keep) { - num_generated_++; + if (num_generated_scalar_ < num_max_generation_) { + num_generated_scalar_++; } // z = z + beta * p // p = z diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 23625039db9..2dc75169e44 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -24,10 +24,12 @@ namespace solver { /** - * Chebyshev iteration is an iterative method that uses another coarse - * method to approximate the error of the current solution via the current + * Chebyshev iteration is an iterative method that uses another inner + * solver to approximate the error of the current solution via the current * residual. It has another term for the difference of solution. Moreover, this - * method requires knowledge about the spectrum of the matrix. + * method requires knowledge about the spectrum of the matrix. This + * implementation follows the algorithm in "Templates for the Solution of Linear + * Systems: Building Blocks for Iterative Methods, 2nd Edition". * * ``` * solution = initial_guess @@ -127,7 +129,8 @@ class Chebyshev : public EnableLinOp>, GKO_FACTORY_PARAMETER_VECTOR(criteria, nullptr); /** - * Inner solver factory. + * Inner solver factory. If not provided this will result in a + * non-preconditioned Chebyshev iteration. */ std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( solver, nullptr); @@ -152,11 +155,6 @@ class Chebyshev : public EnableLinOp>, */ initial_guess_mode GKO_FACTORY_PARAMETER_SCALAR( default_initial_guess, initial_guess_mode::provided); - - /** - * The number of scalar to keep - */ - int GKO_FACTORY_PARAMETER_SCALAR(num_keep, 2); }; GKO_ENABLE_LIN_OP_FACTORY(Chebyshev, parameters, Factory); GKO_ENABLE_BUILD_METHOD(Factory); @@ -186,38 +184,16 @@ class Chebyshev : public EnableLinOp>, {} explicit Chebyshev(const Factory* factory, - std::shared_ptr system_matrix) - : EnableLinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnableSolverBase{std::move(system_matrix)}, - EnableIterativeBase{ - stop::combine(factory->get_parameters().criteria)}, - parameters_{factory->get_parameters()} - { - if (parameters_.generated_solver) { - this->set_solver(parameters_.generated_solver); - } else if (parameters_.solver) { - this->set_solver( - parameters_.solver->generate(this->get_system_matrix())); - } else { - this->set_solver(matrix::Identity::create( - this->get_executor(), this->get_size())); - } - this->set_default_initial_guess(parameters_.default_initial_guess); - center_ = - (std::get<0>(parameters_.foci) + std::get<1>(parameters_.foci)) / - ValueType{2}; - // the absolute value of foci_direction is the focal direction - foci_direction_ = - (std::get<1>(parameters_.foci) - std::get<0>(parameters_.foci)) / - ValueType{2}; - // if changing the lower/upper eig, need to reset it to zero - num_generated_ = 0; - } + std::shared_ptr system_matrix); private: std::shared_ptr solver_{}; - mutable int num_generated_; + // num_generated_scalar_ is to track the number of generated scalar alpha + // and beta. + mutable size_type num_generated_scalar_; + // num_max_generation_ is the number of keeping the generated scalar in + // workspace. + mutable size_type num_max_generation_; ValueType center_; ValueType foci_direction_; }; diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index 1be6a85c318..800cd802a0a 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -49,7 +49,39 @@ class Chebyshev : public ::testing::Test { TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); -TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) + +TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithoutIteration) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + auto upper = value_type{1.1}; + auto lower = value_type{0.9}; + auto factory = + Solver::build() + .with_criteria(gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(this->exec)) + .with_foci(lower, upper) + .on(this->exec); + auto solver = 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()); + + auto alpha = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::alpha)); + auto beta = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::beta)); + // if the stop criterion does not contain iteration limit, it will use the + // default value. + ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); + ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 4})); +} + + +TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithLessIteration) { using Mtx = typename TestFixture::Mtx; using Solver = typename TestFixture::Solver; @@ -59,9 +91,11 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) auto factory = Solver::build() .with_criteria( - gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(this->exec), + gko::stop::Iteration::build().with_max_iters(1u).on(this->exec)) .with_foci(lower, upper) - .with_num_keep(3) .on(this->exec); auto solver = factory->generate(this->mtx); auto b = gko::initialize({3.9, 9.0, 2.2}, this->exec); @@ -73,8 +107,39 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) solver->get_workspace_op(gko::solver::workspace_traits::alpha)); auto beta = gko::as>( solver->get_workspace_op(gko::solver::workspace_traits::beta)); + // if the iteration limit less than the default value, it will use the + // default value. ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 4})); +} + + +TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + auto upper = value_type{1.1}; + auto lower = value_type{0.9}; + auto factory = + Solver::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) + .with_foci(lower, upper) + .on(this->exec); + auto solver = 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()); + + auto alpha = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::alpha)); + auto beta = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::beta)); + // the iteration is more than default + ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 7})); + ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 7})); // check the num_keep alpha, beta auto d = (upper + lower) / value_type{2}; auto c = (upper - lower) / value_type{2}; @@ -91,6 +156,80 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) } +TYPED_TEST(Chebyshev, NumAlphaBetaFromChangingCriterion) +{ + using Mtx = typename TestFixture::Mtx; + using Solver = typename TestFixture::Solver; + using value_type = typename TestFixture::value_type; + auto upper = value_type{1.1}; + auto lower = value_type{0.9}; + auto factory = + Solver::build() + .with_criteria( + gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value) + .on(this->exec), + gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) + .with_foci(lower, upper) + .on(this->exec); + auto solver = 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); + + // same as previous test, but it works with combined factory + solver->apply(b.get(), x.get()); + + auto alpha = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::alpha)); + auto beta = gko::as>( + solver->get_workspace_op(gko::solver::workspace_traits::beta)); + // if the iteration limit is less than the default value, it will use the + // default value. + ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 7})); + ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 7})); + { + // Set less iteration limit + solver->set_stop_criterion_factory( + gko::stop::Iteration::build().with_max_iters(4u).on(this->exec)); + + solver->apply(b.get(), x.get()); + + auto alpha_tmp = + gko::as>(solver->get_workspace_op( + gko::solver::workspace_traits::alpha)); + auto beta_tmp = + gko::as>(solver->get_workspace_op( + gko::solver::workspace_traits::beta)); + // if the iteration limit is less than the previous one, it keeps the + // storage. + ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 7})); + ASSERT_EQ(beta_tmp->get_size(), (gko::dim<2>{1, 7})); + ASSERT_EQ(alpha_tmp->get_const_values(), alpha->get_const_values()); + ASSERT_EQ(beta_tmp->get_const_values(), beta->get_const_values()); + } + { + // Set more iteration limit + solver->set_stop_criterion_factory( + gko::stop::Iteration::build().with_max_iters(10u).on(this->exec)); + + solver->apply(b.get(), x.get()); + + auto alpha_tmp = + gko::as>(solver->get_workspace_op( + gko::solver::workspace_traits::alpha)); + auto beta_tmp = + gko::as>(solver->get_workspace_op( + gko::solver::workspace_traits::beta)); + // if the iteration limit is less than the previous one, it keeps the + // storage. + ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 11})); + ASSERT_EQ(beta_tmp->get_size(), (gko::dim<2>{1, 11})); + ASSERT_NE(alpha_tmp->get_const_values(), alpha->get_const_values()); + ASSERT_NE(beta_tmp->get_const_values(), beta->get_const_values()); + } +} + + TYPED_TEST(Chebyshev, SolvesTriangularSystem) { using Mtx = typename TestFixture::Mtx; From 8bf1344ad895e611f43f1191c878d32426c41844 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 8 Aug 2023 12:05:56 +0200 Subject: [PATCH 06/19] extract the split residual update and update test Co-authored-by: Marcel Koch --- core/solver/chebyshev.cpp | 57 ++++-------------- core/solver/ir.cpp | 58 +++++------------- core/solver/residual_update.hpp | 75 ++++++++++++++++++++++++ core/test/solver/chebyshev.cpp | 5 +- include/ginkgo/core/solver/chebyshev.hpp | 9 +-- 5 files changed, 108 insertions(+), 96 deletions(-) create mode 100644 core/solver/residual_update.hpp diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index cb5f72b969e..15b25f1d20d 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -9,6 +9,7 @@ #include "core/distributed/helpers.hpp" #include "core/solver/ir_kernels.hpp" +#include "core/solver/residual_update.hpp" #include "core/solver/solver_base.hpp" #include "core/solver/solver_boilerplate.hpp" @@ -235,52 +236,20 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, int iter = -1; 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); - this->template log( - this, 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); + 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; - } + }; + bool all_stopped = residual_update( + this, iter, one_op, neg_one_op, dense_b, dense_x, residual, + residual_ptr, stop_criterion, relative_stopping_id, stop_status, + one_changed, log_func); + if (all_stopped) { + break; } if (solver_->apply_uses_initial_guess()) { diff --git a/core/solver/ir.cpp b/core/solver/ir.cpp index 074c4a7d572..de60a540249 100644 --- a/core/solver/ir.cpp +++ b/core/solver/ir.cpp @@ -11,6 +11,7 @@ #include "core/config/config_helper.hpp" #include "core/distributed/helpers.hpp" #include "core/solver/ir_kernels.hpp" +#include "core/solver/residual_update.hpp" #include "core/solver/solver_base.hpp" #include "core/solver/solver_boilerplate.hpp" @@ -223,54 +224,23 @@ 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 = residual_update( + this, iter, one_op, neg_one_op, dense_b, dense_x, residual, + residual_ptr, stop_criterion, relative_stopping_id, stop_status, + one_changed, log_func); + if (all_stopped) { + break; } + if (solver_->apply_uses_initial_guess()) { // Use the inner solver to solve // A * inner_solution = residual diff --git a/core/solver/residual_update.hpp b/core/solver/residual_update.hpp new file mode 100644 index 00000000000..8c2e43c7e6e --- /dev/null +++ b/core/solver/residual_update.hpp @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_CORE_SOLVER_RESIDUAL_UPDATE_HPP_ +#define GKO_CORE_SOLVER_RESIDUAL_UPDATE_HPP_ + + +#include +#include +#include + + +namespace gko { +namespace solver { + + +template +bool residual_update(SolverType* solver, int iter, const ScalarType* one_op, + const ScalarType* neg_one_op, const VectorType* dense_b, + VectorType* dense_x, VectorType* residual, + const VectorType*& residual_ptr, + std::unique_ptr& stop_criterion, + uint8 relative_stopping_id, + array& stop_status, bool& one_changed, + LogFunc log) +{ + 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_RESIDUAL_UPDATE_HPP_ diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp index 40b495390ec..547242a56ee 100644 --- a/core/test/solver/chebyshev.cpp +++ b/core/test/solver/chebyshev.cpp @@ -162,10 +162,7 @@ TYPED_TEST(Chebyshev, CanSetInnerSolverInFactory) auto chebyshev_factory = Solver::build() .with_criteria( - gko::stop::Iteration::build().with_max_iters(3u).on(this->exec), - gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(this->exec)) + gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) .with_solver( Solver::build() .with_criteria( diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 2dc75169e44..2a7594e6c16 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -129,8 +129,8 @@ class Chebyshev : public EnableLinOp>, GKO_FACTORY_PARAMETER_VECTOR(criteria, nullptr); /** - * Inner solver factory. If not provided this will result in a - * non-preconditioned Chebyshev iteration. + * Inner solver (preconditioner) factory. If not provided this will + * result in a non-preconditioned Chebyshev iteration. */ std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( solver, nullptr); @@ -143,8 +143,9 @@ class Chebyshev : public EnableLinOp>, generated_solver, nullptr); /** - * The pair of foci of ellipse. It is usually be {lower bound of eigval, - * upper bound of eigval} for real matrices. + * 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 GKO_FACTORY_PARAMETER_VECTOR( foci, value_type{0}, value_type{1}); From c99b9bff325dcbb7634dc477e0c2f2773327a673 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 8 Aug 2023 13:54:40 +0200 Subject: [PATCH 07/19] use preconditioner not inner_solver for chebyshev Co-authored-by: Marcel Koch --- core/solver/chebyshev.cpp | 59 +++++++-------------- core/test/solver/chebyshev.cpp | 39 +++++++------- include/ginkgo/core/solver/chebyshev.hpp | 37 +++++-------- reference/test/solver/chebyshev_kernels.cpp | 4 +- test/test_install/test_install.cpp | 12 ++++- 5 files changed, 63 insertions(+), 88 deletions(-) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 15b25f1d20d..dd9cd2de87f 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -32,20 +32,10 @@ Chebyshev::Chebyshev(const Factory* factory, std::shared_ptr system_matrix) : EnableLinOp(factory->get_executor(), gko::transpose(system_matrix->get_size())), - EnableSolverBase{std::move(system_matrix)}, - EnableIterativeBase{ - stop::combine(factory->get_parameters().criteria)}, + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, parameters_{factory->get_parameters()} { - if (parameters_.generated_solver) { - this->set_solver(parameters_.generated_solver); - } else if (parameters_.solver) { - this->set_solver( - parameters_.solver->generate(this->get_system_matrix())); - } else { - this->set_solver(matrix::Identity::create( - this->get_executor(), this->get_size())); - } this->set_default_initial_guess(parameters_.default_initial_guess); center_ = (std::get<0>(parameters_.foci) + std::get<1>(parameters_.foci)) / ValueType{2}; @@ -58,30 +48,17 @@ Chebyshev::Chebyshev(const Factory* factory, } -template -void Chebyshev::set_solver(std::shared_ptr new_solver) -{ - auto exec = this->get_executor(); - if (new_solver) { - GKO_ASSERT_EQUAL_DIMENSIONS(new_solver, this); - GKO_ASSERT_IS_SQUARE_MATRIX(new_solver); - if (new_solver->get_executor() != exec) { - new_solver = gko::clone(exec, new_solver); - } - } - solver_ = new_solver; -} - - template Chebyshev& Chebyshev::operator=(const Chebyshev& other) { if (&other != this) { EnableLinOp::operator=(other); - EnableSolverBase::operator=(other); - EnableIterativeBase::operator=(other); + EnablePreconditionedIterativeSolver< + ValueType, Chebyshev>::operator=(other); this->parameters_ = other.parameters_; - this->set_solver(other.get_solver()); + // the workspace is not copied. + this->num_generated_scalar_ = 0; + this->num_max_generation_ = 3; } return *this; } @@ -92,11 +69,11 @@ Chebyshev& Chebyshev::operator=(Chebyshev&& other) { if (&other != this) { EnableLinOp::operator=(std::move(other)); - EnableSolverBase::operator=(std::move(other)); - EnableIterativeBase::operator=(std::move(other)); - this->parameters_ = std::exchange(other.parameters_, parameters_type{}); - this->set_solver(other.get_solver()); - other.set_solver(nullptr); + EnablePreconditionedIterativeSolver< + ValueType, Chebyshev>::operator=(std::move(other)); + // the workspace is not moved. + this->num_generated_scalar_ = 0; + this->num_max_generation_ = 3; } return *this; } @@ -122,8 +99,8 @@ template std::unique_ptr Chebyshev::transpose() const { return build() - .with_generated_solver( - share(as(this->get_solver())->transpose())) + .with_generated_preconditioner( + share(as(this->get_preconditioner())->transpose())) .with_criteria(this->get_stop_criterion_factory()) .with_foci(parameters_.foci) .on(this->get_executor()) @@ -136,8 +113,8 @@ template std::unique_ptr Chebyshev::conj_transpose() const { return build() - .with_generated_solver( - share(as(this->get_solver())->conj_transpose())) + .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))) @@ -252,13 +229,13 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, break; } - if (solver_->apply_uses_initial_guess()) { + 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); } - solver_->apply(residual_ptr, inner_solution); + this->get_preconditioner()->apply(residual_ptr, inner_solution); size_type index = (iter >= num_max_generation_) ? num_max_generation_ : iter; auto alpha_scalar = diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp index 547242a56ee..d88381364f5 100644 --- a/core/test/solver/chebyshev.cpp +++ b/core/test/solver/chebyshev.cpp @@ -71,9 +71,9 @@ TYPED_TEST(Chebyshev, ChebyshevFactoryCreatesCorrectSolver) { using Solver = typename TestFixture::Solver; ASSERT_EQ(this->solver->get_size(), gko::dim<2>(3, 3)); - auto cg_solver = static_cast(this->solver.get()); - ASSERT_NE(cg_solver->get_system_matrix(), nullptr); - ASSERT_EQ(cg_solver->get_system_matrix(), this->mtx); + auto solver = static_cast(this->solver.get()); + ASSERT_NE(solver->get_system_matrix(), nullptr); + ASSERT_EQ(solver->get_system_matrix(), this->mtx); } @@ -163,7 +163,7 @@ TYPED_TEST(Chebyshev, CanSetInnerSolverInFactory) Solver::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) - .with_solver( + .with_preconditioner( Solver::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(3u).on( @@ -171,12 +171,12 @@ TYPED_TEST(Chebyshev, CanSetInnerSolverInFactory) .on(this->exec)) .on(this->exec); auto solver = chebyshev_factory->generate(this->mtx); - auto inner_solver = dynamic_cast( - static_cast(solver.get())->get_solver().get()); + auto preconditioner = dynamic_cast( + static_cast(solver.get())->get_preconditioner().get()); - ASSERT_NE(inner_solver, nullptr); - ASSERT_EQ(inner_solver->get_size(), gko::dim<2>(3, 3)); - ASSERT_EQ(inner_solver->get_system_matrix(), this->mtx); + ASSERT_NE(preconditioner, nullptr); + ASSERT_EQ(preconditioner->get_size(), gko::dim<2>(3, 3)); + ASSERT_EQ(preconditioner->get_system_matrix(), this->mtx); } @@ -194,13 +194,13 @@ TYPED_TEST(Chebyshev, CanSetGeneratedInnerSolverInFactory) Solver::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) - .with_generated_solver(chebyshev_solver) + .with_generated_preconditioner(chebyshev_solver) .on(this->exec); auto solver = chebyshev_factory->generate(this->mtx); - auto inner_solver = solver->get_solver(); + auto preconditioner = solver->get_preconditioner(); - ASSERT_NE(inner_solver.get(), nullptr); - ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); + ASSERT_NE(preconditioner.get(), nullptr); + ASSERT_EQ(preconditioner.get(), chebyshev_solver.get()); } @@ -246,7 +246,7 @@ TYPED_TEST(Chebyshev, ThrowsOnWrongInnerSolverInFactory) Solver::build() .with_criteria( gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) - .with_generated_solver(chebyshev_solver) + .with_generated_preconditioner(chebyshev_solver) .on(this->exec); ASSERT_THROW(chebyshev_factory->generate(this->mtx), @@ -270,11 +270,11 @@ TYPED_TEST(Chebyshev, CanSetInnerSolver) gko::stop::Iteration::build().with_max_iters(3u).on(this->exec)) .on(this->exec); auto solver = chebyshev_factory->generate(this->mtx); - solver->set_solver(chebyshev_solver); - auto inner_solver = solver->get_solver(); + solver->set_preconditioner(chebyshev_solver); + auto preconditioner = solver->get_preconditioner(); - ASSERT_NE(inner_solver.get(), nullptr); - ASSERT_EQ(inner_solver.get(), chebyshev_solver.get()); + ASSERT_NE(preconditioner.get(), nullptr); + ASSERT_EQ(preconditioner.get(), chebyshev_solver.get()); } @@ -320,7 +320,8 @@ TYPED_TEST(Chebyshev, ThrowOnWrongInnerSolverSet) .on(this->exec); auto solver = chebyshev_factory->generate(this->mtx); - ASSERT_THROW(solver->set_solver(chebyshev_solver), gko::DimensionMismatch); + ASSERT_THROW(solver->set_preconditioner(chebyshev_solver), + gko::DimensionMismatch); } diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 2a7594e6c16..2112408ede9 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -35,7 +35,7 @@ namespace solver { * solution = initial_guess * while not converged: * residual = b - A solution - * error = solver(A, residual) + * error = preconditioner(A) * residual * solution = solution + alpha_i * error + beta_i * (solution_i - * solution_{i-1}) * ``` @@ -47,11 +47,12 @@ namespace solver { * @ingroup LinOp */ template -class Chebyshev : public EnableLinOp>, - public EnableSolverBase>, - public EnableIterativeBase>, - public EnableApplyWithInitialGuess>, - public Transposable { +class Chebyshev + : public EnableLinOp>, + public EnablePreconditionedIterativeSolver>, + public EnableApplyWithInitialGuess>, + public Transposable { friend class EnableLinOp; friend class EnablePolymorphicObject; friend class EnableApplyWithInitialGuess; @@ -75,20 +76,6 @@ class Chebyshev : public EnableLinOp>, initial_guess_mode::provided; } - /** - * Returns the solver operator used as the inner solver. - * - * @return the solver operator used as the inner solver - */ - std::shared_ptr get_solver() const { return solver_; } - - /** - * Sets the solver operator used as the inner solver. - * - * @param new_solver the new inner solver - */ - void set_solver(std::shared_ptr new_solver); - /** * Copy-assigns a Chebyshev solver. Preserves the executor, shallow-copies * inner solver, stopping criterion and system matrix. If the executors @@ -129,18 +116,18 @@ class Chebyshev : public EnableLinOp>, GKO_FACTORY_PARAMETER_VECTOR(criteria, nullptr); /** - * Inner solver (preconditioner) factory. If not provided this will + * Preconditioner factory. If not provided this will * result in a non-preconditioned Chebyshev iteration. */ std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( - solver, nullptr); + preconditioner, nullptr); /** - * Already generated solver. If one is provided, the factory `solver` - * will be ignored. + * Already generated preconditioner. If one is provided, the factory + * `preconditioner` will be ignored. */ std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( - generated_solver, nullptr); + generated_preconditioner, nullptr); /** * The pair of foci of ellipse, which covers the eigenvalues of diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index 800cd802a0a..21750e0f81f 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -308,7 +308,7 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemWithIterativeInnerSolver) using value_type = typename TestFixture::value_type; const gko::remove_complex inner_reduction_factor = 1e-2; - auto inner_solver_factory = gko::share( + auto precond_factory = gko::share( gko::solver::Gmres::build() .with_criteria(gko::stop::ResidualNorm::build() .with_reduction_factor(inner_reduction_factor) @@ -322,7 +322,7 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemWithIterativeInnerSolver) gko::stop::ResidualNorm::build() .with_reduction_factor(r::value) .on(this->exec)) - .with_solver(inner_solver_factory) + .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); 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<>; From 10c939f96c819d956f73801b8ada1dd978139cc2 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Wed, 9 Aug 2023 09:21:21 +0200 Subject: [PATCH 08/19] add chebyshev test and fix the number of generated --- core/solver/chebyshev.cpp | 5 + test/solver/CMakeLists.txt | 1 + test/solver/chebyshev_kernels.cpp | 180 ++++++++++++++++++++++++++++++ test/solver/solver.cpp | 22 +++- 4 files changed, 205 insertions(+), 3 deletions(-) create mode 100644 test/solver/chebyshev_kernels.cpp diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index dd9cd2de87f..278df722f55 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -164,6 +164,7 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, GKO_SOLVER_VECTOR(inner_solution, dense_b); GKO_SOLVER_VECTOR(update_solution, dense_b); + auto old_num_max_generation = num_max_generation_; // Use the scalar first // get the iteration information from stopping criterion. if (auto combined = @@ -182,6 +183,10 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, num_max_generation_ = std::max(num_max_generation_, iter_stop->get_parameters().max_iters); } + // Regenerate the vector if we realloc the memory. + if (old_num_max_generation != num_max_generation_) { + num_generated_scalar_ = 0; + } auto alpha = this->template create_workspace_scalar( GKO_SOLVER_TRAITS::alpha, num_max_generation_ + 1); auto beta = this->template create_workspace_scalar( 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..5a23e52da20 --- /dev/null +++ b/test/solver/chebyshev_kernels.cpp @@ -0,0 +1,180 @@ +/************************************************************* +Copyright (c) 2017-2023, the Ginkgo authors +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its +contributors may be used to endorse or promote products derived from +this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*************************************************************/ + +#include + + +#include + + +#include +#include +#include +#include +#include +#include +#include + + +#include "core/test/utils.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, 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 = clone(exec, mtx); + auto d_x = clone(exec, x); + auto d_b = 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)) + .on(ref); + auto d_chebyshev_factory = + gko::solver::Chebyshev::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(2u).on(exec)) + .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 = clone(exec, mtx); + auto d_x = clone(exec, x); + auto d_b = 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).on( + ref)) + .on(ref)) + .with_criteria( + gko::stop::Iteration::build().with_max_iters(2u).on(ref)) + .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).on( + exec)) + .on(exec)) + .with_criteria( + gko::stop::Iteration::build().with_max_iters(2u).on(exec)) + .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 = 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 = clone(exec, x); + auto chebyshev_factory = + gko::solver::Chebyshev::build() + .with_criteria( + gko::stop::Iteration::build().with_max_iters(2u).on(ref)) + .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).on(exec)) + .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..c975d016380 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).on(exec)); + } +}; + + 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, From a63380aa70f571df1bfd9fd19bbfe3702008db30 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 12 Sep 2023 11:43:48 +0200 Subject: [PATCH 09/19] use scalar from ws and move some local variables Co-authored-by: Marcel Koch --- core/solver/chebyshev.cpp | 11 ++++----- core/solver/ir.cpp | 11 ++++----- ...esidual_update.hpp => update_residual.hpp} | 24 +++++++++++-------- 3 files changed, 22 insertions(+), 24 deletions(-) rename core/solver/{residual_update.hpp => update_residual.hpp} (77%) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 278df722f55..2875067251f 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -9,9 +9,9 @@ #include "core/distributed/helpers.hpp" #include "core/solver/ir_kernels.hpp" -#include "core/solver/residual_update.hpp" #include "core/solver/solver_base.hpp" #include "core/solver/solver_boilerplate.hpp" +#include "core/solver/update_residual.hpp" namespace gko { @@ -155,7 +155,6 @@ void Chebyshev::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(); @@ -198,7 +197,6 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, auto beta_ref = ValueType{0.5} * (foci_direction_ * alpha_ref) * (foci_direction_ * alpha_ref); - bool one_changed{}; auto& stop_status = this->template create_workspace_array( ws::stop, dense_b->get_size()[1]); exec->run(chebyshev::make_initialize(&stop_status)); @@ -226,10 +224,9 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, solver, dense_b, dense_x, iter, residual_ptr, nullptr, nullptr, &stop_status, all_stopped); }; - bool all_stopped = residual_update( - this, iter, one_op, neg_one_op, dense_b, dense_x, residual, - residual_ptr, stop_criterion, relative_stopping_id, stop_status, - one_changed, log_func); + bool all_stopped = update_residual( + this, iter, dense_b, dense_x, residual, residual_ptr, + stop_criterion, stop_status, log_func); if (all_stopped) { break; } diff --git a/core/solver/ir.cpp b/core/solver/ir.cpp index de60a540249..cbf1bcacce6 100644 --- a/core/solver/ir.cpp +++ b/core/solver/ir.cpp @@ -11,9 +11,9 @@ #include "core/config/config_helper.hpp" #include "core/distributed/helpers.hpp" #include "core/solver/ir_kernels.hpp" -#include "core/solver/residual_update.hpp" #include "core/solver/solver_base.hpp" #include "core/solver/solver_boilerplate.hpp" +#include "core/solver/update_residual.hpp" namespace gko { @@ -193,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(); @@ -203,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)); @@ -232,10 +230,9 @@ void Ir::apply_dense_impl(const VectorType* dense_b, solver, dense_b, dense_x, iter, residual_ptr, nullptr, nullptr, &stop_status, all_stopped); }; - bool all_stopped = residual_update( - this, iter, one_op, neg_one_op, dense_b, dense_x, residual, - residual_ptr, stop_criterion, relative_stopping_id, stop_status, - one_changed, log_func); + bool all_stopped = update_residual( + this, iter, dense_b, dense_x, residual, residual_ptr, + stop_criterion, stop_status, log_func); if (all_stopped) { break; } diff --git a/core/solver/residual_update.hpp b/core/solver/update_residual.hpp similarity index 77% rename from core/solver/residual_update.hpp rename to core/solver/update_residual.hpp index 8c2e43c7e6e..8105694d899 100644 --- a/core/solver/residual_update.hpp +++ b/core/solver/update_residual.hpp @@ -2,8 +2,8 @@ // // SPDX-License-Identifier: BSD-3-Clause -#ifndef GKO_CORE_SOLVER_RESIDUAL_UPDATE_HPP_ -#define GKO_CORE_SOLVER_RESIDUAL_UPDATE_HPP_ +#ifndef GKO_CORE_SOLVER_UPDATE_RESIDUAL_HPP_ +#define GKO_CORE_SOLVER_UPDATE_RESIDUAL_HPP_ #include @@ -15,17 +15,21 @@ namespace gko { namespace solver { -template -bool residual_update(SolverType* solver, int iter, const ScalarType* one_op, - const ScalarType* neg_one_op, const VectorType* dense_b, +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, - uint8 relative_stopping_id, - array& stop_status, bool& one_changed, - LogFunc log) + 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 = @@ -72,4 +76,4 @@ bool residual_update(SolverType* solver, int iter, const ScalarType* one_op, } // namespace solver } // namespace gko -#endif // GKO_CORE_SOLVER_RESIDUAL_UPDATE_HPP_ +#endif // GKO_CORE_SOLVER_UPDATE_RESIDUAL_HPP_ From 882fc5e27b65e4f3b4ce866e9c35108afb34d2fa Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Thu, 11 Jan 2024 13:04:35 +0800 Subject: [PATCH 10/19] review update Co-authored-by: Marcel Koch --- core/solver/chebyshev.cpp | 42 ++++++++-------- core/solver/ir.cpp | 1 - core/test/solver/chebyshev.cpp | 55 ++++++++------------ include/ginkgo/core/solver/chebyshev.hpp | 48 +++++++----------- reference/test/solver/chebyshev_kernels.cpp | 56 ++++++++------------- test/solver/chebyshev_kernels.cpp | 1 - 6 files changed, 81 insertions(+), 122 deletions(-) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 2875067251f..d0f0b3e2d2b 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -42,9 +42,6 @@ Chebyshev::Chebyshev(const Factory* factory, foci_direction_ = (std::get<1>(parameters_.foci) - std::get<0>(parameters_.foci)) / ValueType{2}; - // if changing the lower/upper eig, need to reset it to zero - num_generated_scalar_ = 0; - num_max_generation_ = 3; } @@ -56,9 +53,6 @@ Chebyshev& Chebyshev::operator=(const Chebyshev& other) EnablePreconditionedIterativeSolver< ValueType, Chebyshev>::operator=(other); this->parameters_ = other.parameters_; - // the workspace is not copied. - this->num_generated_scalar_ = 0; - this->num_max_generation_ = 3; } return *this; } @@ -71,9 +65,6 @@ Chebyshev& Chebyshev::operator=(Chebyshev&& other) EnableLinOp::operator=(std::move(other)); EnablePreconditionedIterativeSolver< ValueType, Chebyshev>::operator=(std::move(other)); - // the workspace is not moved. - this->num_generated_scalar_ = 0; - this->num_max_generation_ = 3; } return *this; } @@ -147,6 +138,20 @@ void Chebyshev::apply_with_initial_guess_impl( } +template +void visit_criteria(Fn&& fn, + std::shared_ptr c) +{ + fn(c); + if (auto combined = + std::dynamic_pointer_cast(c)) { + for (const auto& factory : combined->get_parameters().criteria) { + visit_criteria(std::forward(fn), factory); + } + } +} + + template template void Chebyshev::apply_dense_impl(const VectorType* dense_b, @@ -166,22 +171,15 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, auto old_num_max_generation = num_max_generation_; // Use the scalar first // get the iteration information from stopping criterion. - if (auto combined = - std::dynamic_pointer_cast( - this->get_stop_criterion_factory())) { - for (const auto& factory : combined->get_parameters().criteria) { - if (auto iter_stop = std::dynamic_pointer_cast< + visit_criteria( + [&](auto factory) { + if (auto iter = std::dynamic_pointer_cast< const gko::stop::Iteration::Factory>(factory)) { num_max_generation_ = std::max( - num_max_generation_, iter_stop->get_parameters().max_iters); + num_max_generation_, iter->get_parameters().max_iters); } - } - } else if (auto iter_stop = std::dynamic_pointer_cast< - const gko::stop::Iteration::Factory>( - this->get_stop_criterion_factory())) { - num_max_generation_ = std::max(num_max_generation_, - iter_stop->get_parameters().max_iters); - } + }, + this->get_stop_criterion_factory()); // Regenerate the vector if we realloc the memory. if (old_num_max_generation != num_max_generation_) { num_generated_scalar_ = 0; diff --git a/core/solver/ir.cpp b/core/solver/ir.cpp index cbf1bcacce6..c51365f0f0a 100644 --- a/core/solver/ir.cpp +++ b/core/solver/ir.cpp @@ -237,7 +237,6 @@ void Ir::apply_dense_impl(const VectorType* dense_b, break; } - if (solver_->apply_uses_initial_guess()) { // Use the inner solver to solve // A * inner_solution = residual diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp index d88381364f5..2bfa667f0df 100644 --- a/core/test/solver/chebyshev.cpp +++ b/core/test/solver/chebyshev.cpp @@ -33,10 +33,9 @@ class Chebyshev : public ::testing::Test { chebyshev_factory( Solver::build() .with_criteria( - gko::stop::Iteration::build().with_max_iters(3u).on(exec), + gko::stop::Iteration::build().with_max_iters(3u), gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(exec)) + .with_reduction_factor(r::value)) .on(exec)), solver(chebyshev_factory->generate(mtx)) {} @@ -70,8 +69,9 @@ TYPED_TEST(Chebyshev, ChebyshevFactoryKnowsItsExecutor) TYPED_TEST(Chebyshev, ChebyshevFactoryCreatesCorrectSolver) { using Solver = typename TestFixture::Solver; - ASSERT_EQ(this->solver->get_size(), gko::dim<2>(3, 3)); 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); } @@ -111,6 +111,7 @@ 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)); @@ -123,6 +124,7 @@ TYPED_TEST(Chebyshev, CanBeCloned) TYPED_TEST(Chebyshev, CanBeCleared) { using Solver = typename TestFixture::Solver; + this->solver->clear(); ASSERT_EQ(this->solver->get_size(), gko::dim<2>(0, 0)); @@ -142,10 +144,10 @@ 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).on(this->exec)) + .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); @@ -159,16 +161,12 @@ 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).on(this->exec)) - .with_preconditioner( - Solver::build() - .with_criteria( - gko::stop::Iteration::build().with_max_iters(3u).on( - this->exec)) - .on(this->exec)) + .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( @@ -185,15 +183,13 @@ 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)) + .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)) + .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); @@ -237,15 +233,13 @@ TYPED_TEST(Chebyshev, ThrowsOnWrongInnerSolverInFactory) 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)) + .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)) + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) .with_generated_preconditioner(chebyshev_solver) .on(this->exec); @@ -259,15 +253,13 @@ 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)) + .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)) + .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); @@ -283,13 +275,12 @@ 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).on( - this->exec)) + .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); @@ -308,15 +299,13 @@ TYPED_TEST(Chebyshev, ThrowOnWrongInnerSolverSet) 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)) + .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)) + .with_criteria(gko::stop::Iteration::build().with_max_iters(3u)) .on(this->exec); auto solver = chebyshev_factory->generate(this->mtx); diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 2112408ede9..0456792cdd2 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -24,12 +24,14 @@ namespace solver { /** - * Chebyshev iteration is an iterative method that uses another inner - * solver to approximate the error of the current solution via the current - * residual. It has another term for the difference of solution. Moreover, this - * method requires knowledge about the spectrum of the matrix. This - * implementation follows the algorithm in "Templates for the Solution of Linear - * Systems: Building Blocks for Iterative Methods, 2nd Edition". + * Chebyshev iteration is an iterative method that can solving nonsymeetric + * problems. Chebyshev Iterations avoids the inner products for computation + * which may be the bottleneck for distributed system. Chebyshev Iteration is + * developed based on the Chebyshev polynomials of the first kind. Moreover, + * this method requires knowledge about the spectrum of the preconditioned + * matrix. This implementation follows the algorithm in "Templates for the + * Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd + * Edition". * * ``` * solution = initial_guess @@ -107,28 +109,11 @@ class Chebyshev */ Chebyshev(Chebyshev&&); - GKO_CREATE_FACTORY_PARAMETERS(parameters, Factory) - { - /** - * Criterion factories. - */ - std::vector> - GKO_FACTORY_PARAMETER_VECTOR(criteria, nullptr); - - /** - * Preconditioner factory. If not provided this will - * result in a non-preconditioned Chebyshev iteration. - */ - std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( - preconditioner, nullptr); - - /** - * Already generated preconditioner. If one is provided, the factory - * `preconditioner` will be ignored. - */ - std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( - generated_preconditioner, nullptr); + 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 @@ -144,6 +129,7 @@ class Chebyshev 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); @@ -176,12 +162,12 @@ class Chebyshev private: std::shared_ptr solver_{}; - // num_generated_scalar_ is to track the number of generated scalar alpha + // num_generated_scalar_ tracks the number of generated scalar alpha // and beta. - mutable size_type num_generated_scalar_; - // num_max_generation_ is the number of keeping the generated scalar in + mutable size_type num_generated_scalar_ = 0; + // num_max_generation_ is the number of generated scalar kept in the // workspace. - mutable size_type num_max_generation_; + mutable size_type num_max_generation_ = 3; ValueType center_; ValueType foci_direction_; }; diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index 21750e0f81f..8bcbef18a68 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -33,11 +33,9 @@ class Chebyshev : public ::testing::Test { chebyshev_factory( Solver::build() .with_criteria( - gko::stop::Iteration::build().with_max_iters(30u).on( - exec), + gko::stop::Iteration::build().with_max_iters(30u), gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(exec)) + .with_reduction_factor(r::value)) .with_foci(value_type{0.9}, value_type{1.1}) .on(exec)) {} @@ -60,8 +58,7 @@ TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithoutIteration) auto factory = Solver::build() .with_criteria(gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(this->exec)) + .with_reduction_factor(r::value)) .with_foci(lower, upper) .on(this->exec); auto solver = factory->generate(this->mtx); @@ -90,11 +87,9 @@ TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithLessIteration) auto lower = value_type{0.9}; auto factory = Solver::build() - .with_criteria( - gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(this->exec), - gko::stop::Iteration::build().with_max_iters(1u).on(this->exec)) + .with_criteria(gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value), + gko::stop::Iteration::build().with_max_iters(1u)) .with_foci(lower, upper) .on(this->exec); auto solver = factory->generate(this->mtx); @@ -123,8 +118,7 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) auto lower = value_type{0.9}; auto factory = Solver::build() - .with_criteria( - gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) + .with_criteria(gko::stop::Iteration::build().with_max_iters(6u)) .with_foci(lower, upper) .on(this->exec); auto solver = factory->generate(this->mtx); @@ -156,7 +150,7 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) } -TYPED_TEST(Chebyshev, NumAlphaBetaFromChangingCriterion) +TYPED_TEST(Chebyshev, AlphaBetaFromChangingCriterion) { using Mtx = typename TestFixture::Mtx; using Solver = typename TestFixture::Solver; @@ -165,11 +159,9 @@ TYPED_TEST(Chebyshev, NumAlphaBetaFromChangingCriterion) auto lower = value_type{0.9}; auto factory = Solver::build() - .with_criteria( - gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(this->exec), - gko::stop::Iteration::build().with_max_iters(6u).on(this->exec)) + .with_criteria(gko::stop::ResidualNorm::build() + .with_reduction_factor(r::value), + gko::stop::Iteration::build().with_max_iters(6u)) .with_foci(lower, upper) .on(this->exec); auto solver = factory->generate(this->mtx); @@ -183,6 +175,8 @@ TYPED_TEST(Chebyshev, NumAlphaBetaFromChangingCriterion) solver->get_workspace_op(gko::solver::workspace_traits::alpha)); auto beta = gko::as>( solver->get_workspace_op(gko::solver::workspace_traits::beta)); + auto alpha_ref = alpha->clone(); + auto beta_ref = beta->clone(); // if the iteration limit is less than the default value, it will use the // default value. ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 7})); @@ -206,6 +200,8 @@ TYPED_TEST(Chebyshev, NumAlphaBetaFromChangingCriterion) ASSERT_EQ(beta_tmp->get_size(), (gko::dim<2>{1, 7})); ASSERT_EQ(alpha_tmp->get_const_values(), alpha->get_const_values()); ASSERT_EQ(beta_tmp->get_const_values(), beta->get_const_values()); + GKO_ASSERT_MTX_NEAR(alpha_tmp, alpha_ref, 0.0); + GKO_ASSERT_MTX_NEAR(beta_tmp, beta_ref, 0.0); } { // Set more iteration limit @@ -220,8 +216,8 @@ TYPED_TEST(Chebyshev, NumAlphaBetaFromChangingCriterion) auto beta_tmp = gko::as>(solver->get_workspace_op( gko::solver::workspace_traits::beta)); - // if the iteration limit is less than the previous one, it keeps the - // storage. + // if the iteration limit is more than the previous one, it regenerates + // workspace ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 11})); ASSERT_EQ(beta_tmp->get_size(), (gko::dim<2>{1, 11})); ASSERT_NE(alpha_tmp->get_const_values(), alpha->get_const_values()); @@ -306,22 +302,17 @@ 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)) + .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).on( - this->exec), + .with_criteria(gko::stop::Iteration::build().with_max_iters(30u), gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value) - .on(this->exec)) + .with_reduction_factor(r::value)) .with_preconditioner(precond_factory) .with_foci(value_type{0.9}, value_type{1.1}) .on(this->exec); @@ -489,8 +480,7 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) 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).on(this->exec)) + .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); @@ -499,9 +489,7 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) initial_guess_mode::zero}) { auto solver = gko::solver::Chebyshev::build() - .with_criteria( - gko::stop::Iteration::build().with_max_iters(1u).on( - this->exec)) + .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) diff --git a/test/solver/chebyshev_kernels.cpp b/test/solver/chebyshev_kernels.cpp index 5a23e52da20..ab62b3dcfb8 100644 --- a/test/solver/chebyshev_kernels.cpp +++ b/test/solver/chebyshev_kernels.cpp @@ -110,7 +110,6 @@ TEST_F(Chebyshev, ApplyWithIterativeInnerSolverIsEquivalentToRef) auto d_mtx = clone(exec, mtx); auto d_x = clone(exec, x); auto d_b = clone(exec, b); - auto chebyshev_factory = gko::solver::Chebyshev::build() .with_preconditioner( From 677cf241c7b2816b2a92fd38d9720ddcebc343d7 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Thu, 11 Jan 2024 16:29:55 +0800 Subject: [PATCH 11/19] update with new format --- core/solver/chebyshev.cpp | 3 ++- test/solver/chebyshev_kernels.cpp | 37 +++---------------------------- 2 files changed, 5 insertions(+), 35 deletions(-) diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index d0f0b3e2d2b..f31b294df6e 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -2,9 +2,10 @@ // // SPDX-License-Identifier: BSD-3-Clause +#include "ginkgo/core/solver/chebyshev.hpp" + #include #include -#include #include #include "core/distributed/helpers.hpp" diff --git a/test/solver/chebyshev_kernels.cpp b/test/solver/chebyshev_kernels.cpp index ab62b3dcfb8..e0e1e2df9aa 100644 --- a/test/solver/chebyshev_kernels.cpp +++ b/test/solver/chebyshev_kernels.cpp @@ -1,41 +1,11 @@ -/************************************************************* -Copyright (c) 2017-2023, the Ginkgo authors -All rights reserved. - -Redistribution and use in source and binary forms, with or without -modification, are permitted provided that the following conditions -are met: - -1. Redistributions of source code must retain the above copyright -notice, this list of conditions and the following disclaimer. - -2. Redistributions in binary form must reproduce the above copyright -notice, this list of conditions and the following disclaimer in the -documentation and/or other materials provided with the distribution. - -3. Neither the name of the copyright holder nor the names of its -contributors may be used to endorse or promote products derived from -this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS -IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A -PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT -HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, -SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT -LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, -DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY -THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE -OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*************************************************************/ +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause #include - #include - #include #include #include @@ -44,7 +14,6 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include - #include "core/test/utils.hpp" #include "test/utils/executor.hpp" From 9a8663e081d4ea6f3b8ab569400047cd381a394b Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Fri, 28 Jun 2024 18:35:22 +0200 Subject: [PATCH 12/19] add config for chebyshev --- core/config/config_helper.hpp | 8 +++--- core/config/registry.cpp | 1 + core/config/solver_config.cpp | 2 ++ core/solver/chebyshev.cpp | 23 +++++++++++++++ core/test/config/solver.cpp | 36 ++++++++++++++++++++++++ include/ginkgo/core/solver/chebyshev.hpp | 22 ++++++++++++++- 6 files changed, 87 insertions(+), 5 deletions(-) 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/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index f31b294df6e..c2c639eb32b 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -8,6 +8,7 @@ #include #include +#include "core/config/solver_config.hpp" #include "core/distributed/helpers.hpp" #include "core/solver/ir_kernels.hpp" #include "core/solver/solver_base.hpp" @@ -27,6 +28,28 @@ GKO_REGISTER_OPERATION(initialize, ir::initialize); } // 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, 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/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 0456792cdd2..4c895304438 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -11,6 +11,8 @@ #include #include #include +#include +#include #include #include #include @@ -49,7 +51,7 @@ namespace solver { * @ingroup LinOp */ template -class Chebyshev +class Chebyshev final : public EnableLinOp>, public EnablePreconditionedIterativeSolver>, @@ -133,6 +135,24 @@ class Chebyshev 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; From 4d44a74667c92b1c6dda710b2d03e1e453c349b3 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Fri, 3 Jan 2025 08:43:44 +0100 Subject: [PATCH 13/19] use mixed_type in test properly and update documentation Co-authored-by: Tobias Ribizel --- include/ginkgo/core/solver/chebyshev.hpp | 12 +- reference/test/solver/chebyshev_kernels.cpp | 135 ++++++++++---------- reference/test/solver/ir_kernels.cpp | 46 ++++--- test/solver/chebyshev_kernels.cpp | 32 ++--- 4 files changed, 108 insertions(+), 117 deletions(-) diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 4c895304438..611ab0fdb30 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -26,12 +26,12 @@ namespace solver { /** - * Chebyshev iteration is an iterative method that can solving nonsymeetric - * problems. Chebyshev Iterations avoids the inner products for computation - * which may be the bottleneck for distributed system. Chebyshev Iteration is - * developed based on the Chebyshev polynomials of the first kind. Moreover, - * this method requires knowledge about the spectrum of the preconditioned - * matrix. This implementation follows the algorithm in "Templates for the + * 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". * diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index 8bcbef18a68..c308ebfe364 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -15,9 +15,8 @@ #include "core/test/utils.hpp" - -namespace { - +template +using workspace_traits = gko::solver::workspace_traits; template class Chebyshev : public ::testing::Test { @@ -67,10 +66,10 @@ TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithoutIteration) solver->apply(b.get(), x.get()); - auto alpha = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::alpha)); - auto beta = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::beta)); + auto alpha = + gko::as(solver->get_workspace_op(workspace_traits::alpha)); + auto beta = + gko::as(solver->get_workspace_op(workspace_traits::beta)); // if the stop criterion does not contain iteration limit, it will use the // default value. ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); @@ -98,10 +97,10 @@ TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithLessIteration) solver->apply(b.get(), x.get()); - auto alpha = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::alpha)); - auto beta = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::beta)); + auto alpha = + gko::as(solver->get_workspace_op(workspace_traits::alpha)); + auto beta = + gko::as(solver->get_workspace_op(workspace_traits::beta)); // if the iteration limit less than the default value, it will use the // default value. ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); @@ -127,10 +126,10 @@ TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) solver->apply(b.get(), x.get()); - auto alpha = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::alpha)); - auto beta = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::beta)); + auto alpha = + gko::as(solver->get_workspace_op(workspace_traits::alpha)); + auto beta = + gko::as(solver->get_workspace_op(workspace_traits::beta)); // the iteration is more than default ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 7})); ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 7})); @@ -171,10 +170,10 @@ TYPED_TEST(Chebyshev, AlphaBetaFromChangingCriterion) // same as previous test, but it works with combined factory solver->apply(b.get(), x.get()); - auto alpha = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::alpha)); - auto beta = gko::as>( - solver->get_workspace_op(gko::solver::workspace_traits::beta)); + auto alpha = + gko::as(solver->get_workspace_op(workspace_traits::alpha)); + auto beta = + gko::as(solver->get_workspace_op(workspace_traits::beta)); auto alpha_ref = alpha->clone(); auto beta_ref = beta->clone(); // if the iteration limit is less than the default value, it will use the @@ -188,12 +187,10 @@ TYPED_TEST(Chebyshev, AlphaBetaFromChangingCriterion) solver->apply(b.get(), x.get()); - auto alpha_tmp = - gko::as>(solver->get_workspace_op( - gko::solver::workspace_traits::alpha)); - auto beta_tmp = - gko::as>(solver->get_workspace_op( - gko::solver::workspace_traits::beta)); + auto alpha_tmp = gko::as( + solver->get_workspace_op(workspace_traits::alpha)); + auto beta_tmp = gko::as( + solver->get_workspace_op(workspace_traits::beta)); // if the iteration limit is less than the previous one, it keeps the // storage. ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 7})); @@ -210,12 +207,10 @@ TYPED_TEST(Chebyshev, AlphaBetaFromChangingCriterion) solver->apply(b.get(), x.get()); - auto alpha_tmp = - gko::as>(solver->get_workspace_op( - gko::solver::workspace_traits::alpha)); - auto beta_tmp = - gko::as>(solver->get_workspace_op( - gko::solver::workspace_traits::beta)); + auto alpha_tmp = gko::as( + solver->get_workspace_op(workspace_traits::alpha)); + auto beta_tmp = gko::as( + solver->get_workspace_op(workspace_traits::beta)); // if the iteration limit is more than the previous one, it regenerates // workspace ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 11})); @@ -242,16 +237,16 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystem) TYPED_TEST(Chebyshev, SolvesTriangularSystemMixed) { - using value_type = gko::next_precision; - using Mtx = gko::matrix::Dense; + 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); + 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); + (r_mixed()) * 1e1); } @@ -278,23 +273,26 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemComplex) TYPED_TEST(Chebyshev, SolvesTriangularSystemMixedComplex) { - using value_type = + using mixed_complex_type = gko::to_complex>; - using Mtx = gko::matrix::Dense; + using MixedMtx = gko::matrix::Dense; 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}}, + 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.0, 0.0}, value_type{0.0, 0.0}, value_type{0.0, 0.0}}, + 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({value_type{1.0, -2.0}, value_type{3.0, -6.0}, - value_type{2.0, -4.0}}), - (r_mixed()) * 1e1); + 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); } @@ -361,17 +359,18 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemUsingAdvancedApply) TYPED_TEST(Chebyshev, 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->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); + 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); + GKO_ASSERT_MTX_NEAR(x, l({1.5, 5.0, 2.0}), + (r_mixed()) * 1e1); } @@ -401,26 +400,29 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemUsingAdvancedApplyComplex) TYPED_TEST(Chebyshev, 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->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}}, + 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.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); } @@ -509,6 +511,3 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) GKO_ASSERT_MTX_NEAR(x, ref_x, 0.0); } } - - -} // namespace 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/chebyshev_kernels.cpp b/test/solver/chebyshev_kernels.cpp index e0e1e2df9aa..4c56d8f8301 100644 --- a/test/solver/chebyshev_kernels.cpp +++ b/test/solver/chebyshev_kernels.cpp @@ -53,13 +53,11 @@ TEST_F(Chebyshev, ApplyIsEquivalentToRef) // both executors auto chebyshev_factory = gko::solver::Chebyshev::build() - .with_criteria( - gko::stop::Iteration::build().with_max_iters(2u).on(ref)) + .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)) + .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)); @@ -82,24 +80,16 @@ TEST_F(Chebyshev, ApplyWithIterativeInnerSolverIsEquivalentToRef) auto chebyshev_factory = gko::solver::Chebyshev::build() .with_preconditioner( - gko::solver::Gmres::build() - .with_criteria( - gko::stop::Iteration::build().with_max_iters(1u).on( - ref)) - .on(ref)) - .with_criteria( - gko::stop::Iteration::build().with_max_iters(2u).on(ref)) + 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).on( - exec)) - .on(exec)) - .with_criteria( - gko::stop::Iteration::build().with_max_iters(2u).on(exec)) + 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)); @@ -127,14 +117,12 @@ TEST_F(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) auto d_x = clone(exec, x); auto chebyshev_factory = gko::solver::Chebyshev::build() - .with_criteria( - gko::stop::Iteration::build().with_max_iters(2u).on(ref)) + .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).on(exec)) + .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) .with_default_initial_guess(guess) .on(exec); auto solver = chebyshev_factory->generate(mtx); From 9d12c658df3071ec33df24e4bd0e2b387adb1746 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Fri, 3 Jan 2025 10:41:25 +0100 Subject: [PATCH 14/19] add fused kernel for chebyshev --- common/unified/solver/chebyshev_kernels.cpp | 67 +++++++++++++++++++++ core/device_hooks/common_kernels.inc.cpp | 11 ++++ core/solver/chebyshev.cpp | 21 +++++-- core/solver/chebyshev_kernels.hpp | 58 ++++++++++++++++++ reference/CMakeLists.txt | 1 + reference/solver/chebyshev_kernels.cpp | 63 +++++++++++++++++++ 6 files changed, 215 insertions(+), 6 deletions(-) create mode 100644 common/unified/solver/chebyshev_kernels.cpp create mode 100644 core/solver/chebyshev_kernels.hpp create mode 100644 reference/solver/chebyshev_kernels.cpp diff --git a/common/unified/solver/chebyshev_kernels.cpp b/common/unified/solver/chebyshev_kernels.cpp new file mode 100644 index 00000000000..01b49273fab --- /dev/null +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -0,0 +1,67 @@ +// SPDX-FileCopyrightText: 2025 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/solver/chebyshev_kernels.hpp" + +#include + +#include "common/unified/base/kernel_launch.hpp" + + +namespace gko { +namespace kernels { +namespace GKO_DEVICE_NAMESPACE { +namespace chebyshev { + + +template +void init_update(std::shared_ptr exec, + const ScalarType* alpha, + const matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto col, auto alpha, auto inner_sol, + auto update_sol, auto output) { + const auto inner_val = inner_sol(row, col); + update_sol(row, col) = val; + output(row, col) += alpha_val * inner_val; + }, + output->get_size(), alpha, inner_sol, update_sol, output); +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_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) +{ + run_kernel( + exec, + [] GKO_KERNEL(auto row, auto col, auto alpha, auto beta, auto inner_sol, + auto update_sol, auto output) { + const auto val = + inner_sol(row, col) + beta[0] * update_sol(row, col); + inner_sol(row, col) = val; + update_sol(row, col) = val; + output(row, col) += alpha[0] * val; + }, + output->get_size(), alpha, beta, inner_sol, update_sol, output); +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_TYPE( + GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); + + +} // namespace chebyshev +} // namespace GKO_DEVICE_NAMESPACE +} // namespace kernels +} // namespace gko diff --git a/core/device_hooks/common_kernels.inc.cpp b/core/device_hooks/common_kernels.inc.cpp index f5b41ff19c9..35da5329ace 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_AND_SCALAR_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); +GKO_STUB_VALUE_AND_SCALAR_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); + + +} // namespace chebyshev + + namespace ir { diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index c2c639eb32b..00d9aecbf9e 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -10,6 +10,7 @@ #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" @@ -23,6 +24,8 @@ namespace { GKO_REGISTER_OPERATION(initialize, ir::initialize); +GKO_REGISTER_OPERATION(init_update, chebyshev::init_update); +GKO_REGISTER_OPERATION(update, chebyshev::update); } // anonymous namespace @@ -274,8 +277,12 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, num_generated_scalar_++; } // x = x + alpha * inner_solution - dense_x->add_scaled(alpha_scalar.get(), inner_solution); - update_solution->copy_from(inner_solution); + // update_solultion = inner_solution + exec->run(chebyshev::make_init_update( + alpha_scalar->get_const_values(), + 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 @@ -295,10 +302,12 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, } // z = z + beta * p // p = z - inner_solution->add_scaled(beta_scalar.get(), update_solution); - update_solution->copy_from(inner_solution); - // x + alpha * p - dense_x->add_scaled(alpha_scalar.get(), update_solution); + // x += alpha * p + exec->run(chebyshev::make_update( + alpha_scalar->get_const_values(), beta_scalar->get_const_values(), + gko::detail::get_local(inner_solution), + gko::detail::get_local(update_solution), + gko::detail::get_local(dense_x))); } } diff --git a/core/solver/chebyshev_kernels.hpp b/core/solver/chebyshev_kernels.hpp new file mode 100644 index 00000000000..d8d9b15b708 --- /dev/null +++ b/core/solver/chebyshev_kernels.hpp @@ -0,0 +1,58 @@ +// 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 "core/base/kernel_declaration.hpp" + + +namespace gko { +namespace kernels { +namespace chebyshev { + + +#define GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType, ScalarType) \ + void init_update(std::shared_ptr exec, \ + const ScalarType* alpha, \ + const matrix::Dense* inner_sol, \ + matrix::Dense* update_sol, \ + matrix::Dense* output) + +#define GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType, ScalarType) \ + void update(std::shared_ptr exec, \ + const ScalarType* alpha, const ScalarType* 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, ScalarType); \ + template \ + GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType, ScalarType) + + +} // 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/reference/CMakeLists.txt b/reference/CMakeLists.txt index 82d99105a6e..87858d18812 100644 --- a/reference/CMakeLists.txt +++ b/reference/CMakeLists.txt @@ -58,6 +58,7 @@ target_sources( 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 diff --git a/reference/solver/chebyshev_kernels.cpp b/reference/solver/chebyshev_kernels.cpp new file mode 100644 index 00000000000..dc319f0694b --- /dev/null +++ b/reference/solver/chebyshev_kernels.cpp @@ -0,0 +1,63 @@ +// SPDX-FileCopyrightText: 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 +void init_update(std::shared_ptr exec, + const ScalarType* alpha, + const matrix::Dense* inner_sol, + matrix::Dense* update_sol, + matrix::Dense* output) +{ + const auto alpha_val = alpha[0]; + 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 = inner_sol->at(row, col); + update_sol->at(row, col) = inner_val; + output->at(row, col) += alpha_val * inner_val; + } + } +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_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) +{ + const auto alpha_val = alpha[0]; + const auto beta_val = beta[0]; + 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 = + inner_sol->at(row, col) + beta[0] * update_sol->at(row, col); + inner_sol->at(row, col) = val; + update_sol->at(row, col) = val; + output->at(row, col) += alpha_val * val; + } + } +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_TYPE( + GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); + + +} // namespace chebyshev +} // namespace reference +} // namespace kernels +} // namespace gko From 83df7758284e23a5bd0c0d5c414cf55e6e6d97f7 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 7 Jan 2025 08:36:49 +0100 Subject: [PATCH 15/19] pass alpha/beta by value to kernel and add test --- common/unified/CMakeLists.txt | 1 + common/unified/solver/chebyshev_kernels.cpp | 16 +- core/solver/chebyshev.cpp | 67 +---- core/solver/chebyshev_kernels.hpp | 4 +- include/ginkgo/core/solver/chebyshev.hpp | 14 +- reference/solver/chebyshev_kernels.cpp | 16 +- reference/test/solver/chebyshev_kernels.cpp | 257 ++++++++------------ test/solver/chebyshev_kernels.cpp | 117 ++++++++- 8 files changed, 239 insertions(+), 253 deletions(-) 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 index 01b49273fab..358d62b11fa 100644 --- a/common/unified/solver/chebyshev_kernels.cpp +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -17,7 +17,7 @@ namespace chebyshev { template void init_update(std::shared_ptr exec, - const ScalarType* alpha, + const ScalarType alpha, const matrix::Dense* inner_sol, matrix::Dense* update_sol, matrix::Dense* output) @@ -27,8 +27,8 @@ void init_update(std::shared_ptr exec, [] GKO_KERNEL(auto row, auto col, auto alpha, auto inner_sol, auto update_sol, auto output) { const auto inner_val = inner_sol(row, col); - update_sol(row, col) = val; - output(row, col) += alpha_val * inner_val; + update_sol(row, col) = inner_val; + output(row, col) += alpha * inner_val; }, output->get_size(), alpha, inner_sol, update_sol, output); } @@ -38,9 +38,8 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_TYPE( template -void update(std::shared_ptr exec, - const ScalarType* alpha, const ScalarType* beta, - matrix::Dense* inner_sol, +void update(std::shared_ptr exec, const ScalarType alpha, + const ScalarType beta, matrix::Dense* inner_sol, matrix::Dense* update_sol, matrix::Dense* output) { @@ -48,11 +47,10 @@ void update(std::shared_ptr exec, exec, [] GKO_KERNEL(auto row, auto col, auto alpha, auto beta, auto inner_sol, auto update_sol, auto output) { - const auto val = - inner_sol(row, col) + beta[0] * update_sol(row, col); + const auto val = inner_sol(row, col) + beta * update_sol(row, col); inner_sol(row, col) = val; update_sol(row, col) = val; - output(row, col) += alpha[0] * val; + output(row, col) += alpha * val; }, output->get_size(), alpha, beta, inner_sol, update_sol, output); } diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 00d9aecbf9e..864203e5a69 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -165,20 +165,6 @@ void Chebyshev::apply_with_initial_guess_impl( } -template -void visit_criteria(Fn&& fn, - std::shared_ptr c) -{ - fn(c); - if (auto combined = - std::dynamic_pointer_cast(c)) { - for (const auto& factory : combined->get_parameters().criteria) { - visit_criteria(std::forward(fn), factory); - } - } -} - - template template void Chebyshev::apply_dense_impl(const VectorType* dense_b, @@ -195,27 +181,6 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, GKO_SOLVER_VECTOR(inner_solution, dense_b); GKO_SOLVER_VECTOR(update_solution, dense_b); - auto old_num_max_generation = num_max_generation_; - // Use the scalar first - // get the iteration information from stopping criterion. - visit_criteria( - [&](auto factory) { - if (auto iter = std::dynamic_pointer_cast< - const gko::stop::Iteration::Factory>(factory)) { - num_max_generation_ = std::max( - num_max_generation_, iter->get_parameters().max_iters); - } - }, - this->get_stop_criterion_factory()); - // Regenerate the vector if we realloc the memory. - if (old_num_max_generation != num_max_generation_) { - num_generated_scalar_ = 0; - } - auto alpha = this->template create_workspace_scalar( - GKO_SOLVER_TRAITS::alpha, num_max_generation_ + 1); - auto beta = this->template create_workspace_scalar( - GKO_SOLVER_TRAITS::beta, num_max_generation_ + 1); - GKO_SOLVER_ONE_MINUS_ONE(); auto alpha_ref = ValueType{1} / center_; @@ -263,24 +228,11 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, inner_solution->copy_from(residual_ptr); } this->get_preconditioner()->apply(residual_ptr, inner_solution); - size_type index = - (iter >= num_max_generation_) ? num_max_generation_ : iter; - auto alpha_scalar = - alpha->create_submatrix(span{0, 1}, span{index, index + 1}); - auto beta_scalar = - beta->create_submatrix(span{0, 1}, span{index, index + 1}); if (iter == 0) { - if (num_generated_scalar_ < num_max_generation_) { - alpha_scalar->fill(alpha_ref); - // unused beta for first iteration, but fill zero - beta_scalar->fill(zero()); - num_generated_scalar_++; - } // x = x + alpha * inner_solution // update_solultion = inner_solution exec->run(chebyshev::make_init_update( - alpha_scalar->get_const_values(), - gko::detail::get_local(inner_solution), + alpha_ref, gko::detail::get_local(inner_solution), gko::detail::get_local(update_solution), gko::detail::get_local(dense_x))); continue; @@ -291,21 +243,11 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, (foci_direction_ * alpha_ref / ValueType{2.0}); } alpha_ref = ValueType{1.0} / (center_ - beta_ref / alpha_ref); - // The last one is always the updated one - if (num_generated_scalar_ < num_max_generation_ || - iter >= num_max_generation_) { - alpha_scalar->fill(alpha_ref); - beta_scalar->fill(beta_ref); - } - if (num_generated_scalar_ < num_max_generation_) { - num_generated_scalar_++; - } // z = z + beta * p // p = z // x += alpha * p exec->run(chebyshev::make_update( - alpha_scalar->get_const_values(), beta_scalar->get_const_values(), - gko::detail::get_local(inner_solution), + alpha_ref, beta_ref, gko::detail::get_local(inner_solution), gko::detail::get_local(update_solution), gko::detail::get_local(dense_x))); } @@ -351,7 +293,7 @@ int workspace_traits>::num_arrays(const Solver&) template int workspace_traits>::num_vectors(const Solver&) { - return 7; + return 5; } @@ -360,8 +302,7 @@ std::vector workspace_traits>::op_names( const Solver&) { return { - "residual", "inner_solution", "update_solution", "alpha", "beta", - "one", "minus_one", + "residual", "inner_solution", "update_solution", "one", "minus_one", }; } diff --git a/core/solver/chebyshev_kernels.hpp b/core/solver/chebyshev_kernels.hpp index d8d9b15b708..28a7dd03c75 100644 --- a/core/solver/chebyshev_kernels.hpp +++ b/core/solver/chebyshev_kernels.hpp @@ -22,14 +22,14 @@ namespace chebyshev { #define GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType, ScalarType) \ void init_update(std::shared_ptr exec, \ - const ScalarType* alpha, \ + const ScalarType alpha, \ const matrix::Dense* inner_sol, \ matrix::Dense* update_sol, \ matrix::Dense* output) #define GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType, ScalarType) \ void update(std::shared_ptr exec, \ - const ScalarType* alpha, const ScalarType* beta, \ + const ScalarType alpha, const ScalarType beta, \ matrix::Dense* inner_sol, \ matrix::Dense* update_sol, \ matrix::Dense* output) diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index 611ab0fdb30..fcad32b1bf4 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -182,12 +182,6 @@ class Chebyshev final private: std::shared_ptr solver_{}; - // num_generated_scalar_ tracks the number of generated scalar alpha - // and beta. - mutable size_type num_generated_scalar_ = 0; - // num_max_generation_ is the number of generated scalar kept in the - // workspace. - mutable size_type num_max_generation_ = 3; ValueType center_; ValueType foci_direction_; }; @@ -215,14 +209,10 @@ struct workspace_traits> { constexpr static int inner_solution = 1; // update solution constexpr static int update_solution = 2; - // alpha - constexpr static int alpha = 3; - // beta - constexpr static int beta = 4; // constant 1.0 scalar - constexpr static int one = 5; + constexpr static int one = 3; // constant -1.0 scalar - constexpr static int minus_one = 6; + constexpr static int minus_one = 4; // stopping status array constexpr static int stop = 0; diff --git a/reference/solver/chebyshev_kernels.cpp b/reference/solver/chebyshev_kernels.cpp index dc319f0694b..8150c59c22d 100644 --- a/reference/solver/chebyshev_kernels.cpp +++ b/reference/solver/chebyshev_kernels.cpp @@ -14,17 +14,16 @@ namespace chebyshev { template void init_update(std::shared_ptr exec, - const ScalarType* alpha, + const ScalarType alpha, const matrix::Dense* inner_sol, matrix::Dense* update_sol, matrix::Dense* output) { - const auto alpha_val = alpha[0]; 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 = inner_sol->at(row, col); update_sol->at(row, col) = inner_val; - output->at(row, col) += alpha_val * inner_val; + output->at(row, col) += alpha * inner_val; } } } @@ -34,21 +33,18 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_TYPE( template -void update(std::shared_ptr exec, - const ScalarType* alpha, const ScalarType* beta, - matrix::Dense* inner_sol, +void update(std::shared_ptr exec, const ScalarType alpha, + const ScalarType beta, matrix::Dense* inner_sol, matrix::Dense* update_sol, matrix::Dense* output) { - const auto alpha_val = alpha[0]; - const auto beta_val = beta[0]; 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 = - inner_sol->at(row, col) + beta[0] * update_sol->at(row, col); + inner_sol->at(row, col) + beta * update_sol->at(row, col); inner_sol->at(row, col) = val; update_sol->at(row, col) = val; - output->at(row, col) += alpha_val * val; + output->at(row, col) += alpha * val; } } } diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index c308ebfe364..07d5ff76287 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -2,6 +2,8 @@ // // SPDX-License-Identifier: BSD-3-Clause +#include "core/solver/chebyshev_kernels.hpp" + #include #include @@ -15,9 +17,6 @@ #include "core/test/utils.hpp" -template -using workspace_traits = gko::solver::workspace_traits; - template class Chebyshev : public ::testing::Test { protected: @@ -47,180 +46,140 @@ class Chebyshev : public ::testing::Test { TYPED_TEST_SUITE(Chebyshev, gko::test::ValueTypes, TypenameNameGenerator); -TYPED_TEST(Chebyshev, CheckDefaultNumAlphaBetaWithoutIteration) +TYPED_TEST(Chebyshev, KernelInitUpdate) { - using Mtx = typename TestFixture::Mtx; - using Solver = typename TestFixture::Solver; using value_type = typename TestFixture::value_type; - auto upper = value_type{1.1}; - auto lower = value_type{0.9}; - auto factory = - Solver::build() - .with_criteria(gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value)) - .with_foci(lower, upper) - .on(this->exec); - auto solver = 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); + 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); - solver->apply(b.get(), x.get()); + gko::kernels::reference::chebyshev::init_update( + this->exec, alpha, inner_sol.get(), update_sol.get(), output.get()); - auto alpha = - gko::as(solver->get_workspace_op(workspace_traits::alpha)); - auto beta = - gko::as(solver->get_workspace_op(workspace_traits::beta)); - // if the stop criterion does not contain iteration limit, it will use the - // default value. - ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); - ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 4})); + 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, CheckDefaultNumAlphaBetaWithLessIteration) +TYPED_TEST(Chebyshev, KernelUpdate) { - using Mtx = typename TestFixture::Mtx; - using Solver = typename TestFixture::Solver; using value_type = typename TestFixture::value_type; - auto upper = value_type{1.1}; - auto lower = value_type{0.9}; - auto factory = - Solver::build() - .with_criteria(gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value), - gko::stop::Iteration::build().with_max_iters(1u)) - .with_foci(lower, upper) - .on(this->exec); - auto solver = 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); + 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); - solver->apply(b.get(), x.get()); + gko::kernels::reference::chebyshev::update(this->exec, alpha, beta, + inner_sol.get(), + update_sol.get(), output.get()); - auto alpha = - gko::as(solver->get_workspace_op(workspace_traits::alpha)); - auto beta = - gko::as(solver->get_workspace_op(workspace_traits::beta)); - // if the iteration limit less than the default value, it will use the - // default value. - ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 4})); - ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 4})); + 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); } -TYPED_TEST(Chebyshev, CheckStoredAlphaBeta) +#ifdef GINKGO_MIXED_PRECISION + + +TYPED_TEST(Chebyshev, MixedKernelInitUpdate) { - using Mtx = typename TestFixture::Mtx; - using Solver = typename TestFixture::Solver; using value_type = typename TestFixture::value_type; - auto upper = value_type{1.1}; - auto lower = value_type{0.9}; - auto factory = - Solver::build() - .with_criteria(gko::stop::Iteration::build().with_max_iters(6u)) - .with_foci(lower, upper) - .on(this->exec); - auto solver = 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); + 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); - solver->apply(b.get(), x.get()); + gko::kernels::reference::chebyshev::init_update( + this->exec, alpha, inner_sol.get(), update_sol.get(), output.get()); - auto alpha = - gko::as(solver->get_workspace_op(workspace_traits::alpha)); - auto beta = - gko::as(solver->get_workspace_op(workspace_traits::beta)); - // the iteration is more than default - ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 7})); - ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 7})); - // check the num_keep alpha, beta - auto d = (upper + lower) / value_type{2}; - auto c = (upper - lower) / value_type{2}; - EXPECT_EQ(alpha->at(0, 0), value_type{1} / d); - EXPECT_EQ(beta->at(0, 0), value_type{0}); - EXPECT_EQ(beta->at(0, 1), - value_type{0.5} * (c * alpha->at(0, 0)) * (c * alpha->at(0, 0))); - EXPECT_EQ(alpha->at(0, 1), - value_type{1} / (d - beta->at(0, 1) / alpha->at(0, 0))); - EXPECT_EQ(beta->at(0, 2), (c * alpha->at(0, 1) / value_type{2}) * - (c * alpha->at(0, 1) / value_type{2})); - EXPECT_EQ(alpha->at(0, 2), - value_type{1} / (d - beta->at(0, 2) / alpha->at(0, 1))); + 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, AlphaBetaFromChangingCriterion) +TYPED_TEST(Chebyshev, MixedKernelUpdate) { - using Mtx = typename TestFixture::Mtx; - using Solver = typename TestFixture::Solver; using value_type = typename TestFixture::value_type; - auto upper = value_type{1.1}; - auto lower = value_type{0.9}; - auto factory = - Solver::build() - .with_criteria(gko::stop::ResidualNorm::build() - .with_reduction_factor(r::value), - gko::stop::Iteration::build().with_max_iters(6u)) - .with_foci(lower, upper) - .on(this->exec); - auto solver = 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); + 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); - // same as previous test, but it works with combined factory - solver->apply(b.get(), x.get()); + gko::kernels::reference::chebyshev::update(this->exec, alpha, beta, + inner_sol.get(), + update_sol.get(), output.get()); - auto alpha = - gko::as(solver->get_workspace_op(workspace_traits::alpha)); - auto beta = - gko::as(solver->get_workspace_op(workspace_traits::beta)); - auto alpha_ref = alpha->clone(); - auto beta_ref = beta->clone(); - // if the iteration limit is less than the default value, it will use the - // default value. - ASSERT_EQ(alpha->get_size(), (gko::dim<2>{1, 7})); - ASSERT_EQ(beta->get_size(), (gko::dim<2>{1, 7})); - { - // Set less iteration limit - solver->set_stop_criterion_factory( - gko::stop::Iteration::build().with_max_iters(4u).on(this->exec)); - - solver->apply(b.get(), x.get()); - - auto alpha_tmp = gko::as( - solver->get_workspace_op(workspace_traits::alpha)); - auto beta_tmp = gko::as( - solver->get_workspace_op(workspace_traits::beta)); - // if the iteration limit is less than the previous one, it keeps the - // storage. - ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 7})); - ASSERT_EQ(beta_tmp->get_size(), (gko::dim<2>{1, 7})); - ASSERT_EQ(alpha_tmp->get_const_values(), alpha->get_const_values()); - ASSERT_EQ(beta_tmp->get_const_values(), beta->get_const_values()); - GKO_ASSERT_MTX_NEAR(alpha_tmp, alpha_ref, 0.0); - GKO_ASSERT_MTX_NEAR(beta_tmp, beta_ref, 0.0); - } - { - // Set more iteration limit - solver->set_stop_criterion_factory( - gko::stop::Iteration::build().with_max_iters(10u).on(this->exec)); - - solver->apply(b.get(), x.get()); - - auto alpha_tmp = gko::as( - solver->get_workspace_op(workspace_traits::alpha)); - auto beta_tmp = gko::as( - solver->get_workspace_op(workspace_traits::beta)); - // if the iteration limit is more than the previous one, it regenerates - // workspace - ASSERT_EQ(alpha_tmp->get_size(), (gko::dim<2>{1, 11})); - ASSERT_EQ(beta_tmp->get_size(), (gko::dim<2>{1, 11})); - ASSERT_NE(alpha_tmp->get_const_values(), alpha->get_const_values()); - ASSERT_NE(beta_tmp->get_const_values(), beta->get_const_values()); - } + 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; diff --git a/test/solver/chebyshev_kernels.cpp b/test/solver/chebyshev_kernels.cpp index 4c56d8f8301..55c71bdc68b 100644 --- a/test/solver/chebyshev_kernels.cpp +++ b/test/solver/chebyshev_kernels.cpp @@ -2,6 +2,8 @@ // // SPDX-License-Identifier: BSD-3-Clause +#include "core/solver/chebyshev_kernels.hpp" + #include #include @@ -15,6 +17,7 @@ #include #include "core/test/utils.hpp" +#include "test/utils/common_fixture.hpp" #include "test/utils/executor.hpp" @@ -40,14 +43,112 @@ class Chebyshev : public CommonTestFixture { }; +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 = clone(exec, mtx); - auto d_x = clone(exec, x); - auto d_b = clone(exec, b); + 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 @@ -74,9 +175,9 @@ 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 = clone(exec, mtx); - auto d_x = clone(exec, x); - auto d_b = clone(exec, b); + 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( @@ -110,11 +211,11 @@ TEST_F(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) 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 = clone(exec, b); + 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 = clone(exec, x); + auto d_x = gko::clone(exec, x); auto chebyshev_factory = gko::solver::Chebyshev::build() .with_criteria(gko::stop::Iteration::build().with_max_iters(2u)) From c3728788b2eb7e27a5d6c2ac34eff8d11e7efc36 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Fri, 14 Feb 2025 15:57:59 +0100 Subject: [PATCH 16/19] use double/complex for fosi Co-authored-by: Tobias Ribizel --- common/unified/solver/chebyshev_kernels.cpp | 81 +++++++++++++++++---- core/device_hooks/common_kernels.inc.cpp | 4 +- core/solver/chebyshev.cpp | 11 +-- core/solver/chebyshev_kernels.hpp | 36 +++++---- include/ginkgo/core/solver/chebyshev.hpp | 22 ++++-- reference/solver/chebyshev_kernels.cpp | 44 ++++++----- 6 files changed, 136 insertions(+), 62 deletions(-) diff --git a/common/unified/solver/chebyshev_kernels.cpp b/common/unified/solver/chebyshev_kernels.cpp index 358d62b11fa..dea4ca1cb9a 100644 --- a/common/unified/solver/chebyshev_kernels.cpp +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -4,6 +4,8 @@ #include "core/solver/chebyshev_kernels.hpp" +#include + #include #include "common/unified/base/kernel_launch.hpp" @@ -15,48 +17,95 @@ namespace GKO_DEVICE_NAMESPACE { namespace chebyshev { -template +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 ScalarType alpha, + 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 = inner_sol(row, col); - update_sol(row, col) = inner_val; - output(row, col) += alpha * inner_val; + 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, inner_sol, update_sol, output); + output->get_size(), alpha_val, inner_sol, update_sol, output); } -GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_TYPE( - GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); -template +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 = inner_sol(row, col) + beta * update_sol(row, col); - inner_sol(row, col) = val; - update_sol(row, col) = val; - output(row, col) += alpha * val; + 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, beta, inner_sol, update_sol, output); + output->get_size(), alpha_val, beta_val, inner_sol, update_sol, output); } -GKO_INSTANTIATE_FOR_EACH_VALUE_AND_SCALAR_TYPE( - GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); } // namespace chebyshev diff --git a/core/device_hooks/common_kernels.inc.cpp b/core/device_hooks/common_kernels.inc.cpp index 35da5329ace..1abe27e9558 100644 --- a/core/device_hooks/common_kernels.inc.cpp +++ b/core/device_hooks/common_kernels.inc.cpp @@ -681,8 +681,8 @@ GKO_STUB_CB_GMRES_CONST(GKO_DECLARE_CB_GMRES_SOLVE_KRYLOV_KERNEL); namespace chebyshev { -GKO_STUB_VALUE_AND_SCALAR_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); -GKO_STUB_VALUE_AND_SCALAR_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); +GKO_STUB_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); +GKO_STUB_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); } // namespace chebyshev diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 864203e5a69..93098cdc124 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -173,6 +173,7 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, { using Vector = matrix::Dense; using ws = workspace_traits; + using coeff_type = detail::coeff_type; auto exec = this->get_executor(); this->setup_workspace(); @@ -183,8 +184,8 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, GKO_SOLVER_ONE_MINUS_ONE(); - auto alpha_ref = ValueType{1} / center_; - auto beta_ref = ValueType{0.5} * (foci_direction_ * alpha_ref) * + 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( @@ -239,10 +240,10 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, } // beta_ref for iter == 1 is initialized in the beginning if (iter > 1) { - beta_ref = (foci_direction_ * alpha_ref / ValueType{2.0}) * - (foci_direction_ * alpha_ref / ValueType{2.0}); + beta_ref = (foci_direction_ * alpha_ref / coeff_type{2.0}) * + (foci_direction_ * alpha_ref / coeff_type{2.0}); } - alpha_ref = ValueType{1.0} / (center_ - beta_ref / alpha_ref); + alpha_ref = coeff_type{1.0} / (center_ - beta_ref / alpha_ref); // z = z + beta * p // p = z // x += alpha * p diff --git a/core/solver/chebyshev_kernels.hpp b/core/solver/chebyshev_kernels.hpp index 28a7dd03c75..59c0db1a477 100644 --- a/core/solver/chebyshev_kernels.hpp +++ b/core/solver/chebyshev_kernels.hpp @@ -20,25 +20,31 @@ namespace kernels { namespace chebyshev { -#define GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType, ScalarType) \ - void init_update(std::shared_ptr exec, \ - const ScalarType alpha, \ - const matrix::Dense* inner_sol, \ - matrix::Dense* update_sol, \ +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, ScalarType) \ - void update(std::shared_ptr exec, \ - const ScalarType alpha, const ScalarType beta, \ - matrix::Dense* inner_sol, \ - matrix::Dense* update_sol, \ +#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, ScalarType); \ - template \ - GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType, ScalarType) +#define GKO_DECLARE_ALL_AS_TEMPLATES \ + template \ + GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL(ValueType); \ + template \ + GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL(ValueType) } // namespace chebyshev diff --git a/include/ginkgo/core/solver/chebyshev.hpp b/include/ginkgo/core/solver/chebyshev.hpp index fcad32b1bf4..9bd8fbf8b4d 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -23,8 +23,16 @@ 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 @@ -121,8 +129,11 @@ class Chebyshev final * preconditioned system. It is usually be {lower bound of eigval, upper * bound of eigval} of preconditioned real matrices. */ - std::pair GKO_FACTORY_PARAMETER_VECTOR( - foci, value_type{0}, value_type{1}); + 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 @@ -170,9 +181,6 @@ class Chebyshev final const LinOp* beta, LinOp* x, initial_guess_mode guess) const override; - void set_relaxation_factor( - std::shared_ptr> new_factor); - explicit Chebyshev(std::shared_ptr exec) : EnableLinOp(std::move(exec)) {} @@ -182,8 +190,8 @@ class Chebyshev final private: std::shared_ptr solver_{}; - ValueType center_; - ValueType foci_direction_; + detail::coeff_type center_; + detail::coeff_type foci_direction_; }; diff --git a/reference/solver/chebyshev_kernels.cpp b/reference/solver/chebyshev_kernels.cpp index 8150c59c22d..9bd0347446f 100644 --- a/reference/solver/chebyshev_kernels.cpp +++ b/reference/solver/chebyshev_kernels.cpp @@ -12,45 +12,55 @@ namespace reference { namespace chebyshev { -template +template +using coeff_type = gko::kernels::chebyshev::coeff_type; + + +template void init_update(std::shared_ptr exec, - const ScalarType alpha, + 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 = inner_sol->at(row, col); - update_sol->at(row, col) = inner_val; - output->at(row, col) += alpha * inner_val; + 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_AND_SCALAR_TYPE( - GKO_DECLARE_CHEBYSHEV_INIT_UPDATE_KERNEL); +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, +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 = - inner_sol->at(row, col) + beta * update_sol->at(row, col); - inner_sol->at(row, col) = val; - update_sol->at(row, col) = val; - output->at(row, col) += alpha * val; + 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_AND_SCALAR_TYPE( - GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); +GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_CHEBYSHEV_UPDATE_KERNEL); } // namespace chebyshev From e8cba616587ccb35eb2abd059dced9c8e915f5e6 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Thu, 20 Mar 2025 17:05:47 +0100 Subject: [PATCH 17/19] rename type and variable, update documentation and format --- common/unified/solver/chebyshev_kernels.cpp | 58 ++++---- core/solver/chebyshev.cpp | 63 ++++++--- core/solver/chebyshev_kernels.hpp | 28 ++-- core/test/solver/chebyshev.cpp | 35 ++--- include/ginkgo/core/base/std_extensions.hpp | 12 +- include/ginkgo/core/solver/chebyshev.hpp | 12 +- reference/solver/chebyshev_kernels.cpp | 39 +++--- reference/test/solver/chebyshev_kernels.cpp | 139 +++++--------------- test/solver/chebyshev_kernels.cpp | 60 +-------- test/solver/solver.cpp | 2 +- 10 files changed, 173 insertions(+), 275 deletions(-) diff --git a/common/unified/solver/chebyshev_kernels.cpp b/common/unified/solver/chebyshev_kernels.cpp index dea4ca1cb9a..47a032f2478 100644 --- a/common/unified/solver/chebyshev_kernels.cpp +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -6,7 +6,9 @@ #include +#include #include +#include #include "common/unified/base/kernel_launch.hpp" @@ -17,10 +19,6 @@ namespace GKO_DEVICE_NAMESPACE { namespace chebyshev { -template -using coeff_type = gko::kernels::chebyshev::coeff_type; - - #if GINKGO_DPCPP_SINGLE_MODE @@ -35,14 +33,8 @@ using if_single_only_type = #else -template -struct type_identity { - using type = T; -}; - - template -using if_single_only_type = typename type_identity::type; +using if_single_only_type = xstd::type_identity_t; #endif @@ -50,26 +42,30 @@ using if_single_only_type = typename type_identity::type; template void init_update(std::shared_ptr exec, - const coeff_type alpha, + const solver::detail::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; + using backend_coeff_type = + if_single_only_type>; + // the backend_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 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)); + 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); + static_cast(output(row, col)) + + alpha * inner_val); }, output->get_size(), alpha_val, inner_sol, update_sol, output); } @@ -78,29 +74,33 @@ 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, +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 actual_coeff_type = if_single_only_type>; - using type = device_type; + using backend_coeff_type = + if_single_only_type>; + // the backend_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); + 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)); + 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)) + - static_cast(alpha) * val); + static_cast(output(row, col)) + alpha * val); }, output->get_size(), alpha_val, beta_val, inner_sol, update_sol, output); } diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 93098cdc124..7b4fd54f3e7 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -19,11 +19,21 @@ namespace gko { namespace solver { -namespace chebyshev { +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); @@ -31,6 +41,7 @@ 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, @@ -43,8 +54,9 @@ typename Chebyshev::parameters_type Chebyshev::parse( 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))); + 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( @@ -54,6 +66,12 @@ typename Chebyshev::parameters_type Chebyshev::parse( } +template +Chebyshev::Chebyshev(std::shared_ptr exec) + : EnableLinOp(std::move(exec)) +{} + + template Chebyshev::Chebyshev(const Factory* factory, std::shared_ptr system_matrix) @@ -64,11 +82,13 @@ Chebyshev::Chebyshev(const Factory* factory, 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}; + auto left_foci = std::get<0>(parameters_.foci); + auto right_foci = std::get<1>(parameters_.foci); + GKO_ASSERT(real(left_foci) <= real(right_foci)); + center_ = + (left_foci + right_foci) / solver::detail::coeff_type{2}; foci_direction_ = - (std::get<1>(parameters_.foci) - std::get<0>(parameters_.foci)) / - ValueType{2}; + (right_foci - left_foci) / solver::detail::coeff_type{2}; } @@ -80,6 +100,8 @@ Chebyshev& Chebyshev::operator=(const Chebyshev& other) EnablePreconditionedIterativeSolver< ValueType, Chebyshev>::operator=(other); this->parameters_ = other.parameters_; + this->center_ = other.center_; + this->foci_direction_ = other.foci_direction_; } return *this; } @@ -92,6 +114,11 @@ Chebyshev& Chebyshev::operator=(Chebyshev&& other) 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; } @@ -173,7 +200,7 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, { using Vector = matrix::Dense; using ws = workspace_traits; - using coeff_type = detail::coeff_type; + using coeff_type = solver::detail::coeff_type; auto exec = this->get_executor(); this->setup_workspace(); @@ -184,13 +211,13 @@ void Chebyshev::apply_dense_impl(const VectorType* 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 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(chebyshev::make_initialize(&stop_status)); + 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); @@ -233,22 +260,22 @@ void Chebyshev::apply_dense_impl(const VectorType* dense_b, // x = x + alpha * inner_solution // update_solultion = inner_solution exec->run(chebyshev::make_init_update( - alpha_ref, gko::detail::get_local(inner_solution), + alpha_host, 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 + // beta_host 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}); + beta_host = (foci_direction_ * alpha_host / coeff_type{2.0}) * + (foci_direction_ * alpha_host / coeff_type{2.0}); } - alpha_ref = coeff_type{1.0} / (center_ - beta_ref / alpha_ref); + 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_ref, beta_ref, gko::detail::get_local(inner_solution), + alpha_host, beta_host, gko::detail::get_local(inner_solution), gko::detail::get_local(update_solution), gko::detail::get_local(dense_x))); } diff --git a/core/solver/chebyshev_kernels.hpp b/core/solver/chebyshev_kernels.hpp index 59c0db1a477..9e2af545a6e 100644 --- a/core/solver/chebyshev_kernels.hpp +++ b/core/solver/chebyshev_kernels.hpp @@ -11,6 +11,7 @@ #include #include #include +#include #include "core/base/kernel_declaration.hpp" @@ -20,24 +21,19 @@ 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, \ +#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 coeff_type alpha, \ - const coeff_type beta, \ - matrix::Dense* inner_sol, \ - matrix::Dense* update_sol, \ +#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 \ diff --git a/core/test/solver/chebyshev.cpp b/core/test/solver/chebyshev.cpp index 2bfa667f0df..5c2152d128e 100644 --- a/core/test/solver/chebyshev.cpp +++ b/core/test/solver/chebyshev.cpp @@ -14,9 +14,7 @@ #include #include "core/test/utils.hpp" - - -namespace { +#include "core/test/utils/assertions.hpp" template @@ -44,17 +42,6 @@ class Chebyshev : public ::testing::Test { 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); @@ -87,8 +74,8 @@ TYPED_TEST(Chebyshev, CanBeCopied) 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()); + GKO_ASSERT_MTX_NEAR(static_cast(copy_mtx.get()), + this->mtx.get(), 0.0); } @@ -102,8 +89,8 @@ TYPED_TEST(Chebyshev, CanBeMoved) 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()); + GKO_ASSERT_MTX_NEAR(static_cast(copy_mtx.get()), + this->mtx.get(), 0.0); } @@ -116,8 +103,8 @@ TYPED_TEST(Chebyshev, CanBeCloned) 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()); + GKO_ASSERT_MTX_NEAR(static_cast(clone_mtx.get()), + this->mtx.get(), 0.0); } @@ -144,16 +131,17 @@ 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(value_type{0.2}, value_type{1.2}) + .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(value_type{0.2}, value_type{1.2})); + std::make_pair(coeff_type{0.2}, coeff_type{1.2})); } @@ -324,6 +312,3 @@ TYPED_TEST(Chebyshev, ThrowsOnRectangularMatrixInFactory) ASSERT_THROW(this->chebyshev_factory->generate(rectangular_mtx), gko::DimensionMismatch); } - - -} // namespace 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 index 9bd8fbf8b4d..cc04282e624 100644 --- a/include/ginkgo/core/solver/chebyshev.hpp +++ b/include/ginkgo/core/solver/chebyshev.hpp @@ -28,7 +28,7 @@ namespace detail { template using coeff_type = - std::conditional_t, std::complex, double>; + std::conditional_t(), std::complex, double>; } @@ -126,8 +126,10 @@ class Chebyshev final 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. + * 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> @@ -181,9 +183,7 @@ class Chebyshev final const LinOp* beta, LinOp* x, initial_guess_mode guess) const override; - explicit Chebyshev(std::shared_ptr exec) - : EnableLinOp(std::move(exec)) - {} + explicit Chebyshev(std::shared_ptr exec); explicit Chebyshev(const Factory* factory, std::shared_ptr system_matrix); diff --git a/reference/solver/chebyshev_kernels.cpp b/reference/solver/chebyshev_kernels.cpp index 9bd0347446f..c02406eed36 100644 --- a/reference/solver/chebyshev_kernels.cpp +++ b/reference/solver/chebyshev_kernels.cpp @@ -5,6 +5,7 @@ #include "core/solver/chebyshev_kernels.hpp" #include +#include namespace gko { namespace kernels { @@ -12,25 +13,24 @@ 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 solver::detail::coeff_type alpha, const matrix::Dense* inner_sol, matrix::Dense* update_sol, matrix::Dense* output) { - using type = coeff_type; + // the backend_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)); + 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); + output->at(row, col) = static_cast( + static_cast(output->at(row, col)) + + alpha * inner_val); } } } @@ -40,22 +40,25 @@ 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, + const solver::detail::coeff_type alpha, + const solver::detail::coeff_type beta, matrix::Dense* inner_sol, matrix::Dense* update_sol, matrix::Dense* output) { - using type = coeff_type; + // the backend_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)) + - static_cast(beta) * - static_cast(update_sol->at(row, 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)) + - static_cast(alpha) * val); + output->at(row, col) = static_cast( + static_cast(output->at(row, col)) + + alpha * val); } } } diff --git a/reference/test/solver/chebyshev_kernels.cpp b/reference/test/solver/chebyshev_kernels.cpp index 07d5ff76287..457555bb640 100644 --- a/reference/test/solver/chebyshev_kernels.cpp +++ b/reference/test/solver/chebyshev_kernels.cpp @@ -23,6 +23,7 @@ class Chebyshev : public ::testing::Test { 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( @@ -34,13 +35,29 @@ class Chebyshev : public ::testing::Test { 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)) + .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); @@ -50,22 +67,13 @@ 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()); + this->exec, this->alpha, this->inner_sol.get(), this->update_sol.get(), + this->output.get()); - GKO_ASSERT_MTX_NEAR(update_sol, inner_sol, 0); - GKO_ASSERT_MTX_NEAR(output, + 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}}, @@ -78,29 +86,19 @@ 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::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(update_sol, inner_sol, 0); + GKO_ASSERT_MTX_NEAR(this->update_sol, this->inner_sol, 0); GKO_ASSERT_MTX_NEAR( - inner_sol, + 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(output, + 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}}, @@ -109,77 +107,6 @@ TYPED_TEST(Chebyshev, KernelUpdate) } -#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; @@ -259,6 +186,7 @@ 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() @@ -271,7 +199,7 @@ TYPED_TEST(Chebyshev, SolvesTriangularSystemWithIterativeInnerSolver) gko::stop::ResidualNorm::build() .with_reduction_factor(r::value)) .with_preconditioner(precond_factory) - .with_foci(value_type{0.9}, value_type{1.1}) + .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); @@ -439,10 +367,11 @@ 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(value_type{0.9}, value_type{1.1}) + .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); @@ -451,7 +380,7 @@ TYPED_TEST(Chebyshev, ApplyWithGivenInitialGuessModeIsEquivalentToRef) 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_foci(coeff_type{0.9}, coeff_type{1.1}) .with_default_initial_guess(guess) .on(this->exec) ->generate(this->mtx); diff --git a/test/solver/chebyshev_kernels.cpp b/test/solver/chebyshev_kernels.cpp index 55c71bdc68b..1447ad8cf4a 100644 --- a/test/solver/chebyshev_kernels.cpp +++ b/test/solver/chebyshev_kernels.cpp @@ -24,6 +24,7 @@ class Chebyshev : public CommonTestFixture { protected: using Mtx = gko::matrix::Dense; + using coeff_type = gko::solver::detail::coeff_type; Chebyshev() : rand_engine(30) {} @@ -45,7 +46,7 @@ class Chebyshev : public CommonTestFixture { TEST_F(Chebyshev, KernelInitUpdate) { - value_type alpha(0.5); + 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); @@ -67,8 +68,8 @@ TEST_F(Chebyshev, KernelInitUpdate) TEST_F(Chebyshev, KernelUpdate) { - value_type alpha(0.5); - value_type beta(0.25); + 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); @@ -88,59 +89,6 @@ TEST_F(Chebyshev, KernelUpdate) } -#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); diff --git a/test/solver/solver.cpp b/test/solver/solver.cpp index c975d016380..d81781a433d 100644 --- a/test/solver/solver.cpp +++ b/test/solver/solver.cpp @@ -188,7 +188,7 @@ struct Chebyshev : SimpleSolverTest> { return SimpleSolverTest>:: build(exec, iteration_count, check_residual) .with_preconditioner( - precond_type::build().with_max_block_size(1u).on(exec)); + precond_type::build().with_max_block_size(1u)); } }; From 0351f48dd9887d6eb0e7451d7625ad8b2ec8f5d7 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Tue, 25 Mar 2025 13:57:35 +0100 Subject: [PATCH 18/19] update the type name and add more assert to check foci Co-authored-by: Marcel Koch --- common/unified/solver/chebyshev_kernels.cpp | 18 +++++++++--------- core/solver/chebyshev.cpp | 3 +++ reference/solver/chebyshev_kernels.cpp | 4 ++-- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/common/unified/solver/chebyshev_kernels.cpp b/common/unified/solver/chebyshev_kernels.cpp index 47a032f2478..d78b1cc888b 100644 --- a/common/unified/solver/chebyshev_kernels.cpp +++ b/common/unified/solver/chebyshev_kernels.cpp @@ -47,13 +47,13 @@ void init_update(std::shared_ptr exec, matrix::Dense* update_sol, matrix::Dense* output) { - using backend_coeff_type = + using coeff_type = if_single_only_type>; - // the backend_coeff_type always be the highest precision, so we need + // 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; + using arithmetic_type = device_type; - auto alpha_val = static_cast(alpha); + auto alpha_val = static_cast(alpha); run_kernel( exec, @@ -81,14 +81,14 @@ void update(std::shared_ptr exec, matrix::Dense* update_sol, matrix::Dense* output) { - using backend_coeff_type = + using coeff_type = if_single_only_type>; - // the backend_coeff_type always be the highest precision, so we need + // 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; + using arithmetic_type = device_type; - auto alpha_val = static_cast(alpha); - auto beta_val = static_cast(beta); + auto alpha_val = static_cast(alpha); + auto beta_val = static_cast(beta); run_kernel( exec, diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 7b4fd54f3e7..2021b2a99be 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -85,10 +85,13 @@ Chebyshev::Chebyshev(const Factory* factory, 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_)); } diff --git a/reference/solver/chebyshev_kernels.cpp b/reference/solver/chebyshev_kernels.cpp index c02406eed36..b4f17309d86 100644 --- a/reference/solver/chebyshev_kernels.cpp +++ b/reference/solver/chebyshev_kernels.cpp @@ -20,7 +20,7 @@ void init_update(std::shared_ptr exec, matrix::Dense* update_sol, matrix::Dense* output) { - // the backend_coeff_type always be the highest precision, so we need + // 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++) { @@ -46,7 +46,7 @@ void update(std::shared_ptr exec, matrix::Dense* update_sol, matrix::Dense* output) { - // the backend_coeff_type always be the highest precision, so we need + // 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++) { From 79045621a8dfdd7f531f1b3ee624ed6c0de43c7d Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Wed, 16 Apr 2025 10:10:02 +0200 Subject: [PATCH 19/19] fix the circular dep --- core/solver/update_residual.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/core/solver/update_residual.hpp b/core/solver/update_residual.hpp index 8105694d899..688dcdb4b18 100644 --- a/core/solver/update_residual.hpp +++ b/core/solver/update_residual.hpp @@ -8,6 +8,7 @@ #include #include +#include #include