diff --git a/include/ginkgo/core/base/composition.hpp b/include/ginkgo/core/base/composition.hpp index 0397211d873..c5b574de641 100644 --- a/include/ginkgo/core/base/composition.hpp +++ b/include/ginkgo/core/base/composition.hpp @@ -44,6 +44,21 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. namespace gko { +/** + * CompositeionBase provides the interface + */ +class CompositionBase { +public: + /** + * Returns a list of operators of the composition. + * + * @return a list of operators + */ + virtual const std::vector>& get_operators() + const noexcept = 0; +}; + + /** * The Composition class can be used to compose linear operators `op1, op2, ..., * opn` and obtain the operator `op1 * op2 * ... * opn`. @@ -62,7 +77,8 @@ namespace gko { template class Composition : public EnableLinOp>, public EnableCreateMethod>, - public Transposable { + public Transposable, + public CompositionBase { friend class EnablePolymorphicObject; friend class EnableCreateMethod; @@ -76,7 +92,7 @@ class Composition : public EnableLinOp>, * @return a list of operators */ const std::vector>& get_operators() const - noexcept + noexcept override { return operators_; } diff --git a/include/ginkgo/core/preconditioner/ic_wrapper.hpp b/include/ginkgo/core/preconditioner/ic_wrapper.hpp new file mode 100644 index 00000000000..3e716bd5dd1 --- /dev/null +++ b/include/ginkgo/core/preconditioner/ic_wrapper.hpp @@ -0,0 +1,311 @@ +/************************************************************* +Copyright (c) 2017-2022, 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. +*************************************************************/ + +#ifndef GKO_PUBLIC_CORE_PRECONDITIONER_IC_WRAPPER_HPP_ +#define GKO_PUBLIC_CORE_PRECONDITIONER_IC_WRAPPER_HPP_ + + +#include +#include + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace gko { +namespace preconditioner { + + +/** + * The Incomplete Cholesky (IC) preconditioner solves the equation $LL^H*x = b$ + * for a given lower triangular matrix L and the right hand side b (can contain + * multiple right hand sides). + * + * It allows to set both the solver for L defaulting to solver::LowerTrs, which + * is a direct triangular solvers. The solver for L^H is the + * conjugate-transposed solver for L, ensuring that the preconditioner is + * symmetric and positive-definite. For this L solver, a factory can be provided + * (using `with_l_solver_factory`) to have more control over their behavior. In + * particular, it is possible to use an iterative method for solving the + * triangular systems. The default parameters for an iterative triangluar solver + * are: + * - reduction factor = 1e-4 + * - max iteration = + * Solvers without such criteria can also be used, in which case none are set. + * + * An object of this class can be created with a matrix or a gko::Composition + * containing two matrices. If created with a matrix, it is factorized before + * creating the solver. If a gko::Composition (containing two matrices) is + * used, the first operand will be taken as the L matrix, the second will be + * considered the L^H matrix, which helps to avoid the otherwise necessary + * transposition of L inside the solver. ParIc can be directly used, since it + * orders the factors in the correct way. + * + * @note When providing a gko::Composition, the first matrix must be the lower + * matrix ($L$), and the second matrix must be its conjugate-transpose + * ($L^H$). If they are swapped, solving might crash or return the wrong result. + * + * @note Do not use symmetric solvers (like CG) for the L solver since both + * matrices (L and L^H) are, by design, not symmetric. + * + * @note This class is not thread safe (even a const object is not) because it + * uses an internal cache to accelerate multiple (sequential) applies. + * Using it in parallel can lead to segmentation faults, wrong results + * and other unwanted behavior. + * + * + * @ingroup precond + * @ingroup LinOp + */ +class IcWrapper : public EnableLinOp, public Transposable { + friend class EnableLinOp; + friend class EnablePolymorphicObject; + +public: + using transposed_type = IcWrapper; + + GKO_CREATE_FACTORY_PARAMETERS(parameters, Factory) + { + /** + * Factory for the L solver + */ + std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( + l_solver_factory, nullptr); + + /** + * Factory for the factorization + */ + std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( + factorization_factory, nullptr); + }; + + GKO_ENABLE_LIN_OP_FACTORY(IcWrapper, parameters, Factory); + GKO_ENABLE_BUILD_METHOD(Factory); + + /** + * Returns the solver which is used for the provided L matrix. + * + * @returns the solver which is used for the provided L matrix + */ + std::shared_ptr get_l_solver() const { return l_solver_; } + + /** + * Returns the solver which is used for the L^H matrix. + * + * @returns the solver which is used for the L^H matrix + */ + std::shared_ptr get_lh_solver() const { return lh_solver_; } + + std::unique_ptr transpose() const override + { + std::unique_ptr transposed{ + new transposed_type{this->get_executor()}}; + transposed->set_size(gko::transpose(this->get_size())); + transposed->l_solver_ = + share(as(this->get_lh_solver())->transpose()); + transposed->lh_solver_ = + share(as(this->get_l_solver())->transpose()); + + return std::move(transposed); + } + + std::unique_ptr conj_transpose() const override + { + std::unique_ptr transposed{ + new transposed_type{this->get_executor()}}; + transposed->set_size(gko::transpose(this->get_size())); + transposed->l_solver_ = + share(as(this->get_lh_solver())->conj_transpose()); + transposed->lh_solver_ = + share(as(this->get_l_solver())->conj_transpose()); + + return std::move(transposed); + } + + /** + * Generates a default solver of type SolverType. + * + * Also checks whether SolverType can be assigned a criteria, and if it + * can, it is assigned default values which should be well suited for a + * preconditioner. + */ + template + static std::enable_if_t::value, + std::unique_ptr> + generate_default_solver(const std::shared_ptr& exec) + { + using value_type = typename SolverType::value_type; + constexpr gko::remove_complex default_reduce_residual{1e-4}; + + return SolverType::build() + .with_criteria(gko::stop::ResidualNorm::build() + .with_reduction_factor(default_reduce_residual) + .on(exec)) + .on(exec); + } + + /** + * @copydoc generate_default_solver + */ + template + static std::enable_if_t::value, + std::unique_ptr> + generate_default_solver(const std::shared_ptr& exec) + { + return SolverType::build().on(exec); + } + +protected: + void apply_impl(const LinOp* b, LinOp* x) const override + { + this->set_cache_to(b); + l_solver_->apply(b, cache_.intermediate.get()); + if (lh_solver_->apply_uses_initial_guess()) { + x->copy_from(cache_.intermediate.get()); + } + lh_solver_->apply(cache_.intermediate.get(), x); + } + + void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, + LinOp* x) const override + { + this->set_cache_to(b); + l_solver_->apply(b, cache_.intermediate.get()); + lh_solver_->apply(alpha, cache_.intermediate.get(), beta, x); + } + + explicit IcWrapper(std::shared_ptr exec) + : EnableLinOp(std::move(exec)) + {} + + explicit IcWrapper(const Factory* factory, + std::shared_ptr lin_op) + : EnableLinOp(factory->get_executor(), lin_op->get_size()), + parameters_{factory->get_parameters()} + { + // TODO: Composition should not need + auto comp = std::dynamic_pointer_cast(lin_op); + std::shared_ptr l_factor; + + // build factorization if we weren't passed a composition + if (!comp) { + auto exec = lin_op->get_executor(); + if (!parameters_.factorization_factory) { + GKO_NOT_SUPPORTED(lin_op); + } + auto fact = std::shared_ptr( + parameters_.factorization_factory->generate(lin_op)); + // ensure that the result is a composition + comp = std::dynamic_pointer_cast(fact); + if (!comp) { + GKO_NOT_SUPPORTED(comp); + } + } + // comp must contain one or two factors + if (comp->get_operators().size() > 2 || comp->get_operators().empty()) { + GKO_NOT_SUPPORTED(comp); + } + l_factor = comp->get_operators()[0]; + GKO_ASSERT_IS_SQUARE_MATRIX(l_factor); + + auto exec = this->get_executor(); + if (comp->get_operators().size() == 2) { + // If comp contains both factors: use the transposed factor to avoid + // transposing twice + auto lh_factor = comp->get_operators()[1]; + GKO_ASSERT_EQUAL_DIMENSIONS(l_factor, lh_factor); + } + + if (!parameters_.l_solver_factory) { + GKO_NOT_SUPPORTED(parameters_.l_solver_factory); + } else { + l_solver_ = parameters_.l_solver_factory->generate(l_factor); + lh_solver_ = as(l_solver_)->conj_transpose(); + } + } + + /** + * Prepares the intermediate vector for the solve by creating it and + * by copying the values from `b`, so `b` acts as the initial guess. + * + * @param b Right hand side of the first solve. Also acts as the + * initial guess, meaning the intermediate value will be a copy of b + */ + void set_cache_to(const LinOp* b) const + { + if (cache_.intermediate == nullptr) { + // TODO: create a create_empty from b + cache_.intermediate = b->clone(this->get_executor()); + } + // Use b as the initial guess for the first triangular solve + cache_.intermediate->copy_from(b); + } + +private: + std::shared_ptr l_solver_{}; + std::shared_ptr lh_solver_{}; + /** + * Manages a vector as a cache, so there is no need to allocate one + * every time an intermediate vector is required. Copying an instance + * will only yield an empty object since copying the cached vector would + * not make sense. + * + * @internal The struct is present so the whole class can be copyable + * (could also be done with writing `operator=` and copy + * constructor of the enclosing class by hand) + */ + mutable struct cache_struct { + cache_struct() = default; + ~cache_struct() = default; + cache_struct(const cache_struct&) {} + cache_struct(cache_struct&&) {} + cache_struct& operator=(const cache_struct&) { return *this; } + cache_struct& operator=(cache_struct&&) { return *this; } + std::unique_ptr intermediate{}; + } cache_; +}; + + +} // namespace preconditioner +} // namespace gko + + +#endif // GKO_PUBLIC_CORE_PRECONDITIONER_IC_WRAPPER_HPP_ diff --git a/include/ginkgo/ginkgo.hpp b/include/ginkgo/ginkgo.hpp index cb4eeebe268..31d267dc4ba 100644 --- a/include/ginkgo/ginkgo.hpp +++ b/include/ginkgo/ginkgo.hpp @@ -103,6 +103,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include +#include #include #include #include diff --git a/reference/test/preconditioner/CMakeLists.txt b/reference/test/preconditioner/CMakeLists.txt index dc6f56cfc49..2281dee70ad 100644 --- a/reference/test/preconditioner/CMakeLists.txt +++ b/reference/test/preconditioner/CMakeLists.txt @@ -1,5 +1,6 @@ ginkgo_create_test(ilu) ginkgo_create_test(ic) +ginkgo_create_test(ic_wrapper) ginkgo_create_test(isai_kernels) ginkgo_create_test(jacobi) ginkgo_create_test(jacobi_kernels) diff --git a/reference/test/preconditioner/ic_wrapper.cpp b/reference/test/preconditioner/ic_wrapper.cpp new file mode 100644 index 00000000000..0358f2af39d --- /dev/null +++ b/reference/test/preconditioner/ic_wrapper.cpp @@ -0,0 +1,448 @@ +/************************************************************* +Copyright (c) 2017-2022, 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 + + +#include "core/test/utils.hpp" + + +namespace { + + +template +class IcWrapper : public ::testing::Test { +protected: + using value_type = + typename std::tuple_element<0, decltype(ValueIndexType())>::type; + using index_type = + typename std::tuple_element<1, decltype(ValueIndexType())>::type; + using Mtx = gko::matrix::Csr; + using Vec = gko::matrix::Dense; + using ltr = gko::solver::LowerTrs; + using isai = gko::preconditioner::LowerIsai; + using ic_type = gko::preconditioner::IcWrapper; + using Composition = gko::Composition; + using par_ic = gko::factorization::ParIc; + + IcWrapper() + : exec(gko::ReferenceExecutor::create()), + mtx(gko::initialize({{1., 1., 0.}, {1., 2., 1.}, {0., 1., 2.}}, + exec)), + l_factor(gko::initialize( + {{1., 0., 0.}, {1., 1., 0.}, {0., 1., 1.}}, exec)), + lh_factor(gko::as(l_factor->conj_transpose())), + l_isai_factor(gko::initialize( + {{1., 0., 0.}, {-1., 1., 0.}, {0., -1., 1.}}, exec)), + lh_isai_factor(gko::as(l_isai_factor->conj_transpose())), + l_composition(Composition::create(l_factor)), + l_lh_composition(Composition::create(l_factor, lh_factor)), + ic_pre_factory( + ic_type::build() + .with_l_solver_factory( + ic_type::generate_default_solver(exec)) + .with_factorization_factory( + par_ic::build().with_both_factors(false).on(exec)) + .on(exec)), + ic_isai_pre_factory( + ic_type::build() + .with_l_solver_factory( + ic_type::generate_default_solver(exec)) + .with_factorization_factory( + par_ic::build().with_both_factors(false).on(exec)) + .on(exec)), + tol{r::value} + {} + + std::shared_ptr exec; + std::shared_ptr mtx; + std::shared_ptr l_factor; + std::shared_ptr lh_factor; + std::shared_ptr l_isai_factor; + std::shared_ptr lh_isai_factor; + std::shared_ptr l_composition; + std::shared_ptr l_lh_composition; + std::shared_ptr ic_pre_factory; + std::shared_ptr ic_isai_pre_factory; + gko::remove_complex tol; +}; + +TYPED_TEST_SUITE(IcWrapper, gko::test::ValueIndexTypes, + PairTypenameNameGenerator); + + +TYPED_TEST(IcWrapper, BuildsTwoFactorComposition) +{ + using ltr = typename TestFixture::ltr; + using utr = typename ltr::transposed_type; + auto precond = this->ic_pre_factory->generate(this->l_lh_composition); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_system_matrix(), + this->l_factor, this->tol); + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_lh_solver())->get_system_matrix(), + this->lh_factor, this->tol); +} + + +TYPED_TEST(IcWrapper, BuildsOneFactorComposition) +{ + using ltr = typename TestFixture::ltr; + using utr = typename ltr::transposed_type; + auto precond = this->ic_pre_factory->generate(this->l_composition); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_system_matrix(), + this->l_factor, this->tol); + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_lh_solver())->get_system_matrix(), + this->lh_factor, this->tol); +} + + +TYPED_TEST(IcWrapper, BuildsSymmetricTwoFactor) +{ + using ltr = typename TestFixture::ltr; + using utr = typename ltr::transposed_type; + auto precond = this->ic_pre_factory->generate(this->l_lh_composition); + + // the first factor should be identical, the second not as it was transposed + ASSERT_EQ(gko::as(precond->get_l_solver())->get_system_matrix().get(), + this->l_factor.get()); + ASSERT_NE(gko::as(precond->get_lh_solver())->get_system_matrix().get(), + this->lh_factor.get()); +} + + +TYPED_TEST(IcWrapper, BuildsFromMatrix) +{ + using ltr = typename TestFixture::ltr; + using utr = typename ltr::transposed_type; + auto precond = this->ic_pre_factory->generate(this->mtx); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_system_matrix(), + this->l_factor, this->tol); + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_lh_solver())->get_system_matrix(), + this->lh_factor, this->tol); +} + + +TYPED_TEST(IcWrapper, BuildsIsaiFromMatrix) +{ + using isai = typename TestFixture::isai; + using isai_trans = typename isai::transposed_type; + auto precond = this->ic_isai_pre_factory->generate(this->mtx); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_approximate_inverse(), + this->l_isai_factor, this->tol); + GKO_ASSERT_MTX_NEAR(gko::as(precond->get_lh_solver()) + ->get_approximate_inverse(), + this->lh_isai_factor, this->tol); +} + + +TYPED_TEST(IcWrapper, ThrowOnWrongCompositionInput) +{ + using Composition = typename TestFixture::Composition; + std::shared_ptr composition = + Composition::create(this->l_factor, this->l_factor, this->l_factor); + + ASSERT_THROW(this->ic_pre_factory->generate(composition), + gko::NotSupported); +} + + +TYPED_TEST(IcWrapper, CanBeCopied) +{ + using Mtx = typename TestFixture::Mtx; + using Composition = typename TestFixture::Composition; + auto ic = this->ic_pre_factory->generate(this->l_lh_composition); + auto before_l_solver = ic->get_l_solver(); + auto before_lh_solver = ic->get_lh_solver(); + // The switch up of matrices is intentional, to make sure they are distinct! + auto lh_l_composition = + gko::share(Composition::create(this->lh_factor, this->l_factor)); + auto copied = this->ic_pre_factory->generate(lh_l_composition); + + copied->copy_from(ic.get()); + + ASSERT_EQ(before_l_solver, copied->get_l_solver()); + ASSERT_EQ(before_lh_solver, copied->get_lh_solver()); +} + + +TYPED_TEST(IcWrapper, CanBeMoved) +{ + using Mtx = typename TestFixture::Mtx; + using Composition = typename TestFixture::Composition; + auto ic = this->ic_pre_factory->generate(this->l_lh_composition); + auto before_l_solver = ic->get_l_solver(); + auto before_lh_solver = ic->get_lh_solver(); + // The switch up of matrices is intentional, to make sure they are distinct! + auto lh_l_composition = + gko::share(Composition::create(this->lh_factor, this->l_factor)); + auto moved = this->ic_pre_factory->generate(lh_l_composition); + + moved->copy_from(std::move(ic)); + + ASSERT_EQ(before_l_solver, moved->get_l_solver()); + ASSERT_EQ(before_lh_solver, moved->get_lh_solver()); +} + + +TYPED_TEST(IcWrapper, CanBeCloned) +{ + auto ic = this->ic_pre_factory->generate(this->l_lh_composition); + auto before_l_solver = ic->get_l_solver(); + auto before_lh_solver = ic->get_lh_solver(); + + auto clone = ic->clone(); + + ASSERT_EQ(before_l_solver, clone->get_l_solver()); + ASSERT_EQ(before_lh_solver, clone->get_lh_solver()); +} + + +TYPED_TEST(IcWrapper, CanBeTransposed) +{ + using IcWrapper = typename TestFixture::ic_type; + using Mtx = typename TestFixture::Mtx; + using ltr = typename TestFixture::ltr; + using utr = typename ltr::transposed_type; + auto ic = this->ic_pre_factory->generate(this->l_lh_composition); + auto l_ref = + gko::as(gko::as(ic->get_l_solver())->get_system_matrix()); + auto lh_ref = + gko::as(gko::as(ic->get_lh_solver())->get_system_matrix()); + + auto transp = gko::as(ic->transpose()); + + auto l_transp = + gko::as(gko::as(transp->get_l_solver())->get_system_matrix()); + auto lh_transp = gko::as( + gko::as(transp->get_lh_solver())->get_system_matrix()); + GKO_ASSERT_MTX_NEAR(l_transp, l_ref, 0); + GKO_ASSERT_MTX_NEAR(lh_transp, lh_ref, 0); +} + + +TYPED_TEST(IcWrapper, CanBeConjTransposed) +{ + using IcWrapper = typename TestFixture::ic_type; + using Mtx = typename TestFixture::Mtx; + using ltr = typename TestFixture::ltr; + using utr = typename ltr::transposed_type; + auto ic = this->ic_pre_factory->generate(this->l_lh_composition); + auto l_ref = + gko::as(gko::as(ic->get_l_solver())->get_system_matrix()); + auto lh_ref = + gko::as(gko::as(ic->get_lh_solver())->get_system_matrix()); + + auto transp = gko::as(ic->conj_transpose()); + + auto l_transp = + gko::as(gko::as(transp->get_l_solver())->get_system_matrix()); + auto lh_transp = gko::as( + gko::as(transp->get_lh_solver())->get_system_matrix()); + GKO_ASSERT_MTX_NEAR(l_transp, l_ref, 0); + GKO_ASSERT_MTX_NEAR(lh_transp, lh_ref, 0); +} + + +TYPED_TEST(IcWrapper, SolvesSingleRhs) +{ + using Vec = typename TestFixture::Vec; + const auto b = gko::initialize({1.0, 3.0, 6.0}, this->exec); + auto x = Vec::create(this->exec, gko::dim<2>{3, 1}); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({3.0, -2.0, 4.0}), this->tol); +} + + +TYPED_TEST(IcWrapper, SolvesSingleRhsMixed) +{ + using T = typename TestFixture::value_type; + using Vec = gko::matrix::Dense>; + const auto b = gko::initialize({1.0, 3.0, 6.0}, this->exec); + auto x = Vec::create(this->exec, gko::dim<2>{3, 1}); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({3.0, -2.0, 4.0}), this->tol); +} + + +TYPED_TEST(IcWrapper, SolvesSingleRhsComplex) +{ + using Vec = gko::to_complex; + using T = typename Vec::value_type; + const auto b = gko::initialize( + {T{1.0, 2.0}, T{3.0, 6.0}, T{6.0, 12.0}}, this->exec); + auto x = Vec::create(this->exec, gko::dim<2>{3, 1}); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({T{3.0, 6.0}, T{-2.0, -4.0}, T{4.0, 8.0}}), + this->tol); +} + + +TYPED_TEST(IcWrapper, SolvesSingleRhsComplexMixed) +{ + using Vec = gko::matrix::Dense< + gko::next_precision>>; + using T = typename Vec::value_type; + const auto b = gko::initialize( + {T{1.0, 2.0}, T{3.0, 6.0}, T{6.0, 12.0}}, this->exec); + auto x = Vec::create(this->exec, gko::dim<2>{3, 1}); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({T{3.0, 6.0}, T{-2.0, -4.0}, T{4.0, 8.0}}), + this->tol); +} + + +TYPED_TEST(IcWrapper, AdvancedSolvesSingleRhs) +{ + using Vec = typename TestFixture::Vec; + const auto b = gko::initialize({1.0, 3.0, 6.0}, this->exec); + const auto alpha = gko::initialize({2.0}, this->exec); + const auto beta = gko::initialize({-1.0}, this->exec); + auto x = gko::initialize({1.0, 2.0, 3.0}, this->exec); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({5.0, -6.0, 5.0}), this->tol); +} + + +TYPED_TEST(IcWrapper, AdvancedSolvesSingleRhsMixed) +{ + using T = typename TestFixture::value_type; + using Vec = gko::matrix::Dense>; + const auto b = gko::initialize({1.0, 3.0, 6.0}, this->exec); + const auto alpha = gko::initialize({2.0}, this->exec); + const auto beta = gko::initialize({-1.0}, this->exec); + auto x = gko::initialize({1.0, 2.0, 3.0}, this->exec); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({5.0, -6.0, 5.0}), this->tol); +} + + +TYPED_TEST(IcWrapper, AdvancedSolvesSingleRhsComplex) +{ + using Dense = typename TestFixture::Vec; + using DenseComplex = gko::to_complex; + using T = typename DenseComplex::value_type; + const auto b = gko::initialize( + {T{1.0, 2.0}, T{3.0, 6.0}, T{6.0, 12.0}}, this->exec); + const auto alpha = gko::initialize({2.0}, this->exec); + const auto beta = gko::initialize({-1.0}, this->exec); + auto x = gko::initialize( + {T{1.0, 2.0}, T{2.0, 4.0}, T{3.0, 6.0}}, this->exec); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({T{5.0, 10.0}, T{-6.0, -12.0}, T{5.0, 10.0}}), + this->tol); +} + + +TYPED_TEST(IcWrapper, AdvancedSolvesSingleRhsComplexMixed) +{ + using MixedDense = gko::matrix::Dense< + gko::next_precision>; + using MixedDenseComplex = gko::to_complex; + using T = typename MixedDenseComplex::value_type; + const auto b = gko::initialize( + {T{1.0, 2.0}, T{3.0, 6.0}, T{6.0, 12.0}}, this->exec); + const auto alpha = gko::initialize({2.0}, this->exec); + const auto beta = gko::initialize({-1.0}, this->exec); + auto x = gko::initialize( + {T{1.0, 2.0}, T{2.0, 4.0}, T{3.0, 6.0}}, this->exec); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(alpha.get(), b.get(), beta.get(), x.get()); + + GKO_ASSERT_MTX_NEAR(x, l({T{5.0, 10.0}, T{-6.0, -12.0}, T{5.0, 10.0}}), + this->tol); +} + + +TYPED_TEST(IcWrapper, SolvesMultipleRhs) +{ + using Vec = typename TestFixture::Vec; + const auto b = gko::initialize( + {{1.0, 2.0, 3.0}, {3.0, 6.0, 9.0}, {6.0, 12.0, 18.0}}, this->exec); + auto x = Vec::create(this->exec, gko::dim<2>{3, 3}); + auto preconditioner = this->ic_pre_factory->generate(this->mtx); + + preconditioner->apply(b.get(), x.get()); + + GKO_ASSERT_MTX_NEAR( + x, l({{3.0, 6.0, 9.0}, {-2.0, -4.0, -6.0}, {4.0, 8.0, 12.0}}), + this->tol); +} + + +} // namespace diff --git a/test/test_install/test_install.cpp b/test/test_install/test_install.cpp index b3a5ec31ff4..c9840bdb9ec 100644 --- a/test/test_install/test_install.cpp +++ b/test/test_install/test_install.cpp @@ -413,6 +413,16 @@ int main() auto test = gko::multigrid::AmgxPgm<>::build().on(exec); } + // core/preconditioner/ic_wrapper.hpp + { + auto test = gko::preconditioner::IcWrapper<>::build().on(exec); + } + + // core/preconditioner/ic.hpp + { + auto test = gko::preconditioner::Ic<>::build().on(exec); + } + // core/preconditioner/ilu.hpp { auto test = gko::preconditioner::Ilu<>::build().on(exec);