From c6971711c26e129e9ceb0066bbdcb0733259f28f Mon Sep 17 00:00:00 2001 From: Fritz Goebel Date: Tue, 9 Mar 2021 17:16:40 +0100 Subject: [PATCH 1/5] backup --- core/preconditioner/isai.cpp | 22 +++++++----- include/ginkgo/core/preconditioner/isai.hpp | 36 ++++++++++++++++--- .../test/preconditioner/isai_kernels.cpp | 19 ++++++++++ 3 files changed, 64 insertions(+), 13 deletions(-) diff --git a/core/preconditioner/isai.cpp b/core/preconditioner/isai.cpp index 7247b9a04e7..91a35509ea9 100644 --- a/core/preconditioner/isai.cpp +++ b/core/preconditioner/isai.cpp @@ -268,16 +268,19 @@ void Isai::generate_inverse( template std::unique_ptr Isai::transpose() const { + auto exec = this->get_executor(); auto is_spd = IsaiType == isai_type::spd; if (is_spd) { return this->clone(); } - std::unique_ptr transp{ - new transposed_type{this->get_executor()}}; + std::unique_ptr transp{new transposed_type{exec}}; transp->set_size(gko::transpose(this->get_size())); - transp->approximate_inverse_ = - share(as(this->get_approximate_inverse())->transpose()); + auto csr_transp = + share(as(as(this->get_approximate_inverse())->transpose())); + auto ell_transp = Ell::create(exec); + csr_transp->convert_to(ell_transp.get()); + transp->approximate_inverse_ = share(ell_transp); return std::move(transp); } @@ -287,16 +290,19 @@ template std::unique_ptr Isai::conj_transpose() const { + auto exec = this->get_executor(); auto is_spd = IsaiType == isai_type::spd; if (is_spd) { return this->clone(); } - std::unique_ptr transp{ - new transposed_type{this->get_executor()}}; + std::unique_ptr transp{new transposed_type{exec}}; transp->set_size(gko::transpose(this->get_size())); - transp->approximate_inverse_ = - share(as(this->get_approximate_inverse())->conj_transpose()); + auto csr_transp = share( + as(as(this->get_approximate_inverse())->conj_transpose())); + auto ell_transp = Ell::create(exec); + csr_transp->convert_to(ell_transp.get()); + transp->approximate_inverse_ = share(ell_transp); return std::move(transp); } diff --git a/include/ginkgo/core/preconditioner/isai.hpp b/include/ginkgo/core/preconditioner/isai.hpp index d3ff87d7b10..53b9eb0ac48 100644 --- a/include/ginkgo/core/preconditioner/isai.hpp +++ b/include/ginkgo/core/preconditioner/isai.hpp @@ -43,6 +43,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include #include +#include namespace gko { @@ -126,6 +127,7 @@ class Isai : public EnableLinOp>, using Comp = Composition; using Csr = matrix::Csr; using Dense = matrix::Dense; + using Ell = matrix::Ell; static constexpr isai_type type{IsaiType}; /** @@ -139,8 +141,25 @@ class Isai : public EnableLinOp>, Comp, Csr>::type> get_approximate_inverse() const { + auto exec = this->get_executor(); + std::shared_ptr return_inv; + if (IsaiType == isai_type::spd) { + auto ops = as(approximate_inverse_)->get_operators(); + auto inv = as(ops[1]); + auto inv_transp = as(ops[0]); + auto csr_inv = Csr::create(exec); + auto csr_inv_transp = Csr::create(exec); + inv->convert_to(csr_inv.get()); + inv_transp->convert_to(csr_inv_transp.get()); + return_inv = share(Composition::create( + share(csr_inv_transp), share(csr_inv))); + } else { + auto csr_inv = Csr::create(exec); + as(approximate_inverse_)->convert_to(csr_inv.get()); + return_inv = share(csr_inv); + } return as::type>(approximate_inverse_); + Csr>::type>(return_inv); } GKO_CREATE_FACTORY_PARAMETERS(parameters, Factory) @@ -211,15 +230,22 @@ class Isai : public EnableLinOp>, : EnableLinOp(factory->get_executor(), system_matrix->get_size()), parameters_{factory->get_parameters()} { + const auto exec = this->get_executor(); const auto skip_sorting = parameters_.skip_sorting; const auto power = parameters_.sparsity_power; const auto excess_limit = parameters_.excess_limit; generate_inverse(system_matrix, skip_sorting, power, excess_limit); + auto inv = share(as(approximate_inverse_)); + auto ell_inv = Ell::create(exec); + inv->move_to(ell_inv.get()); if (IsaiType == isai_type::spd) { - auto inv = share(as(approximate_inverse_)); - auto inv_transp = share(inv->conj_transpose()); - approximate_inverse_ = - Composition::create(inv_transp, inv); + auto inv_transp = share(as(inv->conj_transpose())); + auto ell_inv_transp = Ell::create(exec); + inv_transp->move_to(ell_inv_transp.get()); + approximate_inverse_ = Composition::create( + share(ell_inv_transp), share(ell_inv)); + } else { + approximate_inverse_ = share(ell_inv); } } diff --git a/reference/test/preconditioner/isai_kernels.cpp b/reference/test/preconditioner/isai_kernels.cpp index c501f956cfc..68e9968ef40 100644 --- a/reference/test/preconditioner/isai_kernels.cpp +++ b/reference/test/preconditioner/isai_kernels.cpp @@ -1110,6 +1110,25 @@ TYPED_TEST(Isai, ReturnsCorrectInverseULongrow) auto u_inv = isai->get_approximate_inverse(); + auto inv_d = gko::matrix::Dense::create(this->exec); + inv_d->copy_from(u_inv.get()); + for (auto i = 0; i < u_inv->get_size()[0]; i++) { + for (auto j = 0; j < u_inv->get_size()[1]; j++) { + std::cout << inv_d->at(i, j) << " "; + } + std::cout << std::endl; + } + + std::cout << "----------------------------------------" << std::endl; + auto inv_dd = gko::matrix::Dense::create(this->exec); + inv_dd->copy_from(this->u_csr_longrow.get()); + for (auto i = 0; i < u_inv->get_size()[0]; i++) { + for (auto j = 0; j < u_inv->get_size()[1]; j++) { + std::cout << inv_dd->at(i, j) << " "; + } + std::cout << std::endl; + } + GKO_ASSERT_MTX_EQ_SPARSITY(u_inv, this->u_csr_longrow_inv); GKO_ASSERT_MTX_NEAR(u_inv, this->u_csr_longrow_inv, r::value); } From 3960d29f7dc9de596eed9869dd86dbf9a47ecff4 Mon Sep 17 00:00:00 2001 From: Fritz Goebel Date: Wed, 10 Mar 2021 17:33:36 +0100 Subject: [PATCH 2/5] Make ISAI support reduced storage precision --- core/preconditioner/isai.cpp | 48 +++++---- include/ginkgo/core/base/types.hpp | 19 ++++ include/ginkgo/core/preconditioner/isai.hpp | 100 ++++++++++++------ .../test/preconditioner/isai_kernels.cpp | 39 +++---- 4 files changed, 127 insertions(+), 79 deletions(-) diff --git a/core/preconditioner/isai.cpp b/core/preconditioner/isai.cpp index 91a35509ea9..ab9b63fd9f4 100644 --- a/core/preconditioner/isai.cpp +++ b/core/preconditioner/isai.cpp @@ -117,8 +117,9 @@ std::shared_ptr extend_sparsity(std::shared_ptr &exec, } -template -void Isai::generate_inverse( +template +void Isai::generate_inverse( std::shared_ptr input, bool skip_sorting, int power, IndexType excess_limit) { @@ -265,8 +266,10 @@ void Isai::generate_inverse( } -template -std::unique_ptr Isai::transpose() const +template +std::unique_ptr +Isai::transpose() const { auto exec = this->get_executor(); auto is_spd = IsaiType == isai_type::spd; @@ -278,17 +281,17 @@ std::unique_ptr Isai::transpose() const transp->set_size(gko::transpose(this->get_size())); auto csr_transp = share(as(as(this->get_approximate_inverse())->transpose())); - auto ell_transp = Ell::create(exec); - csr_transp->convert_to(ell_transp.get()); + auto ell_transp = convert_matrix_formats(csr_transp); transp->approximate_inverse_ = share(ell_transp); return std::move(transp); } -template -std::unique_ptr Isai::conj_transpose() - const +template +std::unique_ptr +Isai::conj_transpose() const { auto exec = this->get_executor(); auto is_spd = IsaiType == isai_type::spd; @@ -300,29 +303,28 @@ std::unique_ptr Isai::conj_transpose() transp->set_size(gko::transpose(this->get_size())); auto csr_transp = share( as(as(this->get_approximate_inverse())->conj_transpose())); - auto ell_transp = Ell::create(exec); - csr_transp->convert_to(ell_transp.get()); + auto ell_transp = convert_matrix_formats(csr_transp); transp->approximate_inverse_ = share(ell_transp); return std::move(transp); } -#define GKO_DECLARE_LOWER_ISAI(ValueType, IndexType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI); +#define GKO_DECLARE_LOWER_ISAI(ValueType, IndexType, StorageType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_LOWER_ISAI); -#define GKO_DECLARE_UPPER_ISAI(ValueType, IndexType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UPPER_ISAI); +#define GKO_DECLARE_UPPER_ISAI(ValueType, IndexType, StorageType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_UPPER_ISAI); -#define GKO_DECLARE_GENERAL_ISAI(ValueType, IndexType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GENERAL_ISAI); +#define GKO_DECLARE_GENERAL_ISAI(ValueType, IndexType, StorageType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_GENERAL_ISAI); -#define GKO_DECLARE_SPD_ISAI(ValueType, IndexType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SPD_ISAI); +#define GKO_DECLARE_SPD_ISAI(ValueType, IndexType, StorageType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_SPD_ISAI); } // namespace preconditioner diff --git a/include/ginkgo/core/base/types.hpp b/include/ginkgo/core/base/types.hpp index 2c77aaef826..94d30596d59 100644 --- a/include/ginkgo/core/base/types.hpp +++ b/include/ginkgo/core/base/types.hpp @@ -524,6 +524,25 @@ GKO_ATTRIBUTES constexpr bool operator!=(precision_reduction x, #endif +#define GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(_macro) \ + template _macro(float, int32, float); \ + template _macro(float, int32, double); \ + template _macro(double, int32, double); \ + template _macro(double, int32, float); \ + template _macro(std::complex, int32, std::complex); \ + template _macro(std::complex, int32, std::complex); \ + template _macro(std::complex, int32, std::complex); \ + template _macro(std::complex, int32, std::complex); \ + template _macro(float, int64, float); \ + template _macro(float, int64, double); \ + template _macro(double, int64, double); \ + template _macro(double, int64, float); \ + template _macro(std::complex, int64, std::complex); \ + template _macro(std::complex, int64, std::complex); \ + template _macro(std::complex, int64, std::complex); \ + template _macro(std::complex, int64, std::complex) + + /** * Instantiates a template for each value type conversion pair compiled by * Ginkgo. diff --git a/include/ginkgo/core/preconditioner/isai.hpp b/include/ginkgo/core/preconditioner/isai.hpp index 53b9eb0ac48..321d2be1b6e 100644 --- a/include/ginkgo/core/preconditioner/isai.hpp +++ b/include/ginkgo/core/preconditioner/isai.hpp @@ -103,15 +103,17 @@ enum struct isai_type { lower, upper, general, spd }; * @ingroup precond * @ingroup LinOp */ -template -class Isai : public EnableLinOp>, - public Transposable { +template +class Isai + : public EnableLinOp>, + public Transposable { friend class EnableLinOp; friend class EnablePolymorphicObject; - friend class Isai; - friend class Isai; - friend class Isai; - friend class Isai; + friend class Isai; + friend class Isai; + friend class Isai; + friend class Isai; public: using value_type = ValueType; @@ -123,11 +125,11 @@ class Isai : public EnableLinOp>, ? isai_type::spd : IsaiType == isai_type::lower ? isai_type::upper : isai_type::lower, - ValueType, IndexType>; + ValueType, IndexType, StorageType>; using Comp = Composition; using Csr = matrix::Csr; using Dense = matrix::Dense; - using Ell = matrix::Ell; + using Ell = matrix::Ell; static constexpr isai_type type{IsaiType}; /** @@ -142,21 +144,18 @@ class Isai : public EnableLinOp>, get_approximate_inverse() const { auto exec = this->get_executor(); - std::shared_ptr return_inv; + std::shared_ptr return_inv; if (IsaiType == isai_type::spd) { auto ops = as(approximate_inverse_)->get_operators(); auto inv = as(ops[1]); auto inv_transp = as(ops[0]); - auto csr_inv = Csr::create(exec); - auto csr_inv_transp = Csr::create(exec); - inv->convert_to(csr_inv.get()); - inv_transp->convert_to(csr_inv_transp.get()); - return_inv = share(Composition::create( - share(csr_inv_transp), share(csr_inv))); + auto csr_inv = convert_matrix_formats(inv); + auto csr_inv_transp = convert_matrix_formats(inv_transp); + return_inv = + share(Composition::create(csr_inv_transp, csr_inv)); } else { - auto csr_inv = Csr::create(exec); - as(approximate_inverse_)->convert_to(csr_inv.get()); - return_inv = share(csr_inv); + return_inv = convert_matrix_formats( + as(approximate_inverse_)); } return as::type>(return_inv); @@ -205,6 +204,15 @@ class Isai : public EnableLinOp>, */ std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( excess_solver_factory, nullptr); + + /** + * @brief Parameter for reducing preconditioner storage precision. + * + * Reduces the preconditioner storage precision up to reduce_precision + * times. For value types double and complex, storage precision + * can be reduced up to two times, for float and complex once. + */ + int GKO_FACTORY_PARAMETER_SCALAR(reduce_precision, 0); }; GKO_ENABLE_LIN_OP_FACTORY(Isai, parameters, Factory); @@ -236,12 +244,10 @@ class Isai : public EnableLinOp>, const auto excess_limit = parameters_.excess_limit; generate_inverse(system_matrix, skip_sorting, power, excess_limit); auto inv = share(as(approximate_inverse_)); - auto ell_inv = Ell::create(exec); - inv->move_to(ell_inv.get()); + auto ell_inv = convert_matrix_formats(inv); if (IsaiType == isai_type::spd) { auto inv_transp = share(as(inv->conj_transpose())); - auto ell_inv_transp = Ell::create(exec); - inv_transp->move_to(ell_inv_transp.get()); + auto ell_inv_transp = convert_matrix_formats(inv_transp); approximate_inverse_ = Composition::create( share(ell_inv_transp), share(ell_inv)); } else { @@ -249,6 +255,34 @@ class Isai : public EnableLinOp>, } } + template + std::shared_ptr convert_matrix_formats( + std::shared_ptr in) const + { + if (std::is_same::value) { + return as(in); + } + + auto out = OutType::create(in->get_executor()); + if (std::is_same< + InType, gko::matrix::Csr>::value) { + auto tmp = gko::matrix::Csr< + typename OutType::value_type, + typename InType::index_type>::create(in->get_executor()); + as(in)->convert_to(tmp.get()); + tmp->convert_to(out.get()); + return share(out); + } else { + auto tmp = gko::matrix::Ell< + typename OutType::value_type, + typename InType::index_type>::create(in->get_executor()); + as(in)->convert_to(tmp.get()); + tmp->convert_to(out.get()); + return share(out); + } + } + void apply_impl(const LinOp *b, LinOp *x) const override { approximate_inverse_->apply(b, x); @@ -280,17 +314,21 @@ class Isai : public EnableLinOp>, }; -template -using LowerIsai = Isai; +template +using LowerIsai = Isai; -template -using UpperIsai = Isai; +template +using UpperIsai = Isai; -template -using GeneralIsai = Isai; +template +using GeneralIsai = Isai; -template -using SpdIsai = Isai; +template +using SpdIsai = Isai; } // namespace preconditioner diff --git a/reference/test/preconditioner/isai_kernels.cpp b/reference/test/preconditioner/isai_kernels.cpp index 68e9968ef40..d76932db9bf 100644 --- a/reference/test/preconditioner/isai_kernels.cpp +++ b/reference/test/preconditioner/isai_kernels.cpp @@ -323,7 +323,7 @@ class Isai : public ::testing::Test { std::shared_ptr spd_sparse_inv; }; -TYPED_TEST_SUITE(Isai, gko::test::ValueIndexTypes); +TYPED_TEST_SUITE(Isai, gko::test::RealValueIndexTypes); TYPED_TEST(Isai, KernelGenerateA) @@ -825,7 +825,7 @@ TYPED_TEST(Isai, KernelGenerateULongrow) this->exec, lend(this->u_csr_longrow), lend(result), a1.get_data(), a2.get_data(), false); - GKO_ASSERT_MTX_EQ_SPARSITY(result, this->u_csr_longrow_inv_partial); + // GKO_ASSERT_MTX_EQ_SPARSITY(result, this->u_csr_longrow_inv_partial); GKO_ASSERT_MTX_NEAR(result, this->u_csr_longrow_inv_partial, r::value); GKO_ASSERT_ARRAY_EQ(a1, a1_expect); @@ -1063,7 +1063,9 @@ TYPED_TEST(Isai, ReturnsCorrectInverseLLongrow) auto l_inv = isai->get_approximate_inverse(); - GKO_ASSERT_MTX_EQ_SPARSITY(l_inv, this->l_csr_longrow_inv); + // Disable sparsity check since ISAI is stored in ELL and the conversion + // removes explicitly stored zeros. + // GKO_ASSERT_MTX_EQ_SPARSITY(l_inv, this->l_csr_longrow_inv); GKO_ASSERT_MTX_NEAR(l_inv, this->l_csr_longrow_inv, r::value); } @@ -1081,7 +1083,9 @@ TYPED_TEST(Isai, ReturnsCorrectInverseLLongrowWithExcessSolver) auto l_inv = isai->get_approximate_inverse(); - GKO_ASSERT_MTX_EQ_SPARSITY(l_inv, this->l_csr_longrow_inv); + // Disable sparsity check since ISAI is stored in ELL and the conversion + // removes explicitly stored zeros. + // GKO_ASSERT_MTX_EQ_SPARSITY(l_inv, this->l_csr_longrow_inv); // need to drastically reduce precision due to using different excess solver // factory. GKO_ASSERT_MTX_NEAR(l_inv, this->l_csr_longrow_inv, @@ -1110,26 +1114,9 @@ TYPED_TEST(Isai, ReturnsCorrectInverseULongrow) auto u_inv = isai->get_approximate_inverse(); - auto inv_d = gko::matrix::Dense::create(this->exec); - inv_d->copy_from(u_inv.get()); - for (auto i = 0; i < u_inv->get_size()[0]; i++) { - for (auto j = 0; j < u_inv->get_size()[1]; j++) { - std::cout << inv_d->at(i, j) << " "; - } - std::cout << std::endl; - } - - std::cout << "----------------------------------------" << std::endl; - auto inv_dd = gko::matrix::Dense::create(this->exec); - inv_dd->copy_from(this->u_csr_longrow.get()); - for (auto i = 0; i < u_inv->get_size()[0]; i++) { - for (auto j = 0; j < u_inv->get_size()[1]; j++) { - std::cout << inv_dd->at(i, j) << " "; - } - std::cout << std::endl; - } - - GKO_ASSERT_MTX_EQ_SPARSITY(u_inv, this->u_csr_longrow_inv); + // Disable sparsity check since ISAI is stored in ELL and the conversion + // removes explicitly stored zeros. + // GKO_ASSERT_MTX_EQ_SPARSITY(u_inv, this->u_csr_longrow_inv); GKO_ASSERT_MTX_NEAR(u_inv, this->u_csr_longrow_inv, r::value); } @@ -1147,7 +1134,9 @@ TYPED_TEST(Isai, ReturnsCorrectInverseULongrowWithExcessSolver) auto u_inv = isai->get_approximate_inverse(); - GKO_ASSERT_MTX_EQ_SPARSITY(u_inv, this->u_csr_longrow_inv); + // Disable sparsity check since ISAI is stored in ELL and the conversion + // removes explicitly stored zeros. + // GKO_ASSERT_MTX_EQ_SPARSITY(u_inv, this->u_csr_longrow_inv); // need to drastically reduce precision due to using different excess solver // factory. GKO_ASSERT_MTX_NEAR(u_inv, this->u_csr_longrow_inv, From 5f336019118b7fbf08a5b79f43581cd047b42d7b Mon Sep 17 00:00:00 2001 From: Fritz Goebel Date: Thu, 11 Mar 2021 18:01:26 +0100 Subject: [PATCH 3/5] get rid of extra set of types --- core/preconditioner/isai.cpp | 49 +++++++++++++------ include/ginkgo/core/base/types.hpp | 19 ------- .../test/preconditioner/isai_kernels.cpp | 44 ++++++++++++++++- 3 files changed, 77 insertions(+), 35 deletions(-) diff --git a/core/preconditioner/isai.cpp b/core/preconditioner/isai.cpp index ab9b63fd9f4..80946fff8ae 100644 --- a/core/preconditioner/isai.cpp +++ b/core/preconditioner/isai.cpp @@ -310,21 +310,40 @@ Isai::conj_transpose() const } -#define GKO_DECLARE_LOWER_ISAI(ValueType, IndexType, StorageType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_LOWER_ISAI); - -#define GKO_DECLARE_UPPER_ISAI(ValueType, IndexType, StorageType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_UPPER_ISAI); - -#define GKO_DECLARE_GENERAL_ISAI(ValueType, IndexType, StorageType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_GENERAL_ISAI); - -#define GKO_DECLARE_SPD_ISAI(ValueType, IndexType, StorageType) \ - class Isai -GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(GKO_DECLARE_SPD_ISAI); +#define GKO_DECLARE_LOWER_ISAI1(ValueType, IndexType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI1); + +#define GKO_DECLARE_UPPER_ISAI1(ValueType, IndexType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UPPER_ISAI1); + +#define GKO_DECLARE_GENERAL_ISAI1(ValueType, IndexType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GENERAL_ISAI1); + +#define GKO_DECLARE_SPD_ISAI1(ValueType, IndexType) \ + class Isai +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SPD_ISAI1); + +#define GKO_DECLARE_LOWER_ISAI2(ValueType, IndexType) \ + class Isai> +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI2); + +#define GKO_DECLARE_UPPER_ISAI2(ValueType, IndexType) \ + class Isai> +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UPPER_ISAI2); + +#define GKO_DECLARE_GENERAL_ISAI2(ValueType, IndexType) \ + class Isai> +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GENERAL_ISAI2); + +#define GKO_DECLARE_SPD_ISAI2(ValueType, IndexType) \ + class Isai> +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SPD_ISAI2); } // namespace preconditioner diff --git a/include/ginkgo/core/base/types.hpp b/include/ginkgo/core/base/types.hpp index 94d30596d59..2c77aaef826 100644 --- a/include/ginkgo/core/base/types.hpp +++ b/include/ginkgo/core/base/types.hpp @@ -524,25 +524,6 @@ GKO_ATTRIBUTES constexpr bool operator!=(precision_reduction x, #endif -#define GKO_INSTANTIATE_FOR_EACH_VALUE_INDEX_AND_STORAGE_TYPE(_macro) \ - template _macro(float, int32, float); \ - template _macro(float, int32, double); \ - template _macro(double, int32, double); \ - template _macro(double, int32, float); \ - template _macro(std::complex, int32, std::complex); \ - template _macro(std::complex, int32, std::complex); \ - template _macro(std::complex, int32, std::complex); \ - template _macro(std::complex, int32, std::complex); \ - template _macro(float, int64, float); \ - template _macro(float, int64, double); \ - template _macro(double, int64, double); \ - template _macro(double, int64, float); \ - template _macro(std::complex, int64, std::complex); \ - template _macro(std::complex, int64, std::complex); \ - template _macro(std::complex, int64, std::complex); \ - template _macro(std::complex, int64, std::complex) - - /** * Instantiates a template for each value type conversion pair compiled by * Ginkgo. diff --git a/reference/test/preconditioner/isai_kernels.cpp b/reference/test/preconditioner/isai_kernels.cpp index d76932db9bf..198609bffe6 100644 --- a/reference/test/preconditioner/isai_kernels.cpp +++ b/reference/test/preconditioner/isai_kernels.cpp @@ -73,6 +73,12 @@ class Isai : public ::testing::Test { using GeneralIsai = gko::preconditioner::GeneralIsai; using SpdIsai = gko::preconditioner::SpdIsai; + using GeneralIsai2 = + gko::preconditioner::GeneralIsai>; + using SpdIsai2 = + gko::preconditioner::SpdIsai>; using Mtx = gko::matrix::Csr; using Dense = gko::matrix::Dense; using Csr = gko::matrix::Csr; @@ -201,6 +207,8 @@ class Isai : public ::testing::Test { upper_isai_factory = UpperIsai::build().on(exec); general_isai_factory = GeneralIsai::build().on(exec); spd_isai_factory = SpdIsai::build().on(exec); + general_isai_factory2 = GeneralIsai2::build().on(exec); + spd_isai_factory2 = SpdIsai2::build().on(exec); a_dense->convert_to(lend(a_csr)); a_dense_inv->convert_to(lend(a_csr_inv)); l_dense->convert_to(lend(l_csr)); @@ -268,6 +276,8 @@ class Isai : public ::testing::Test { std::unique_ptr upper_isai_factory; std::unique_ptr general_isai_factory; std::unique_ptr spd_isai_factory; + std::unique_ptr general_isai_factory2; + std::unique_ptr spd_isai_factory2; std::shared_ptr a_dense; std::shared_ptr a_dense_inv; std::shared_ptr l_dense; @@ -323,7 +333,7 @@ class Isai : public ::testing::Test { std::shared_ptr spd_sparse_inv; }; -TYPED_TEST_SUITE(Isai, gko::test::RealValueIndexTypes); +TYPED_TEST_SUITE(Isai, gko::test::ValueIndexTypes); TYPED_TEST(Isai, KernelGenerateA) @@ -1007,6 +1017,18 @@ TYPED_TEST(Isai, ReturnsCorrectInverseA) } +TYPED_TEST(Isai, ReturnsCorrectInverseAMixed) +{ + using value_type = typename TestFixture::value_type; + const auto isai = this->general_isai_factory2->generate(this->a_sparse); + + auto l_inv = isai->get_approximate_inverse(); + + GKO_ASSERT_MTX_EQ_SPARSITY(l_inv, this->a_sparse_inv); + GKO_ASSERT_MTX_NEAR(l_inv, this->a_sparse_inv, r::value); +} + + TYPED_TEST(Isai, ReturnsCorrectInverseALongrow) { using value_type = typename TestFixture::value_type; @@ -1224,6 +1246,26 @@ TYPED_TEST(Isai, ReturnsCorrectInverseSpd) } +TYPED_TEST(Isai, ReturnsCorrectInverseSpdMixed) +{ + using Csr = typename TestFixture::Csr; + using value_type = typename TestFixture::value_type; + const auto isai = this->spd_isai_factory2->generate(this->spd_sparse); + const auto expected_transpose = + gko::as(this->spd_sparse_inv->transpose()); + + // In the spd case, the approximate inverse is a composition of L^T and L. + const auto composition = isai->get_approximate_inverse()->get_operators(); + const auto lower_t = gko::as(composition[0]); + const auto lower = gko::as(composition[1]); + + GKO_ASSERT_MTX_EQ_SPARSITY(lower, this->spd_sparse_inv); + GKO_ASSERT_MTX_EQ_SPARSITY(lower_t, expected_transpose); + GKO_ASSERT_MTX_NEAR(lower, this->spd_sparse_inv, r::value); + GKO_ASSERT_MTX_NEAR(lower_t, expected_transpose, r::value); +} + + TYPED_TEST(Isai, ReturnsCorrectInverseSpdLongrow) { using Csr = typename TestFixture::Csr; From 414eb8fbecf6751f748ee95f5d1ef741e897f824 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Gr=C3=BCtzmacher?= Date: Tue, 23 Mar 2021 10:25:10 +0100 Subject: [PATCH 4/5] Improve format conversion code in ISAI --- core/preconditioner/isai.cpp | 28 ++--- include/ginkgo/core/preconditioner/isai.hpp | 116 ++++++++++++-------- 2 files changed, 80 insertions(+), 64 deletions(-) diff --git a/core/preconditioner/isai.cpp b/core/preconditioner/isai.cpp index 80946fff8ae..967ac10ec30 100644 --- a/core/preconditioner/isai.cpp +++ b/core/preconditioner/isai.cpp @@ -271,18 +271,16 @@ template Isai::transpose() const { - auto exec = this->get_executor(); - auto is_spd = IsaiType == isai_type::spd; - if (is_spd) { + if (IsaiType == isai_type::spd) { return this->clone(); } - std::unique_ptr transp{new transposed_type{exec}}; + std::unique_ptr transp{ + new transposed_type{this->get_executor()}}; transp->set_size(gko::transpose(this->get_size())); - auto csr_transp = - share(as(as(this->get_approximate_inverse())->transpose())); - auto ell_transp = convert_matrix_formats(csr_transp); - transp->approximate_inverse_ = share(ell_transp); + + transp->approximate_inverse_ = + convert_csr_to_ell(this->get_approximate_inverse()->transpose().get()); return std::move(transp); } @@ -293,18 +291,16 @@ template Isai::conj_transpose() const { - auto exec = this->get_executor(); - auto is_spd = IsaiType == isai_type::spd; - if (is_spd) { + if (IsaiType == isai_type::spd) { return this->clone(); } - std::unique_ptr transp{new transposed_type{exec}}; + std::unique_ptr transp{ + new transposed_type{this->get_executor()}}; transp->set_size(gko::transpose(this->get_size())); - auto csr_transp = share( - as(as(this->get_approximate_inverse())->conj_transpose())); - auto ell_transp = convert_matrix_formats(csr_transp); - transp->approximate_inverse_ = share(ell_transp); + + transp->approximate_inverse_ = convert_csr_to_ell( + this->get_approximate_inverse()->conj_transpose().get()); return std::move(transp); } diff --git a/include/ginkgo/core/preconditioner/isai.hpp b/include/ginkgo/core/preconditioner/isai.hpp index 321d2be1b6e..3ebac28cfa6 100644 --- a/include/ginkgo/core/preconditioner/isai.hpp +++ b/include/ginkgo/core/preconditioner/isai.hpp @@ -119,13 +119,11 @@ class Isai using value_type = ValueType; using index_type = IndexType; using transposed_type = - Isai; + using Comp = Composition; using Csr = matrix::Csr; using Dense = matrix::Dense; @@ -147,15 +145,12 @@ class Isai std::shared_ptr return_inv; if (IsaiType == isai_type::spd) { auto ops = as(approximate_inverse_)->get_operators(); - auto inv = as(ops[1]); - auto inv_transp = as(ops[0]); - auto csr_inv = convert_matrix_formats(inv); - auto csr_inv_transp = convert_matrix_formats(inv_transp); + auto csr_inv = convert_ell_to_csr(ops[1].get()); + auto csr_inv_transp = convert_ell_to_csr(ops[0].get()); return_inv = share(Composition::create(csr_inv_transp, csr_inv)); } else { - return_inv = convert_matrix_formats( - as(approximate_inverse_)); + return_inv = convert_ell_to_csr(approximate_inverse_.get()); } return as::type>(return_inv); @@ -243,43 +238,16 @@ class Isai const auto power = parameters_.sparsity_power; const auto excess_limit = parameters_.excess_limit; generate_inverse(system_matrix, skip_sorting, power, excess_limit); - auto inv = share(as(approximate_inverse_)); - auto ell_inv = convert_matrix_formats(inv); - if (IsaiType == isai_type::spd) { - auto inv_transp = share(as(inv->conj_transpose())); - auto ell_inv_transp = convert_matrix_formats(inv_transp); - approximate_inverse_ = Composition::create( - share(ell_inv_transp), share(ell_inv)); - } else { - approximate_inverse_ = share(ell_inv); - } - } - template - std::shared_ptr convert_matrix_formats( - std::shared_ptr in) const - { - if (std::is_same::value) { - return as(in); - } - - auto out = OutType::create(in->get_executor()); - if (std::is_same< - InType, gko::matrix::Csr>::value) { - auto tmp = gko::matrix::Csr< - typename OutType::value_type, - typename InType::index_type>::create(in->get_executor()); - as(in)->convert_to(tmp.get()); - tmp->convert_to(out.get()); - return share(out); + auto csr_inv = as(approximate_inverse_); + auto ell_inv = convert_csr_to_ell(approximate_inverse_.get()); + if (IsaiType == isai_type::spd) { + auto csr_inv_transp = csr_inv->conj_transpose(); + auto ell_inv_transp = convert_csr_to_ell(csr_inv_transp.get()); + approximate_inverse_ = + Composition::create(ell_inv_transp, ell_inv); } else { - auto tmp = gko::matrix::Ell< - typename OutType::value_type, - typename InType::index_type>::create(in->get_executor()); - as(in)->convert_to(tmp.get()); - tmp->convert_to(out.get()); - return share(out); + approximate_inverse_ = ell_inv; } } @@ -295,6 +263,59 @@ class Isai } private: + /** + * Converts a given matrix in ELL format to a matrix in CSR format. + * @note ELL format uses StorageType as the storage + * + * @param in Matrix in ELL format + * + * @returns the given matrix converted to the CSR format. + */ + static std::shared_ptr convert_ell_to_csr(const LinOp *in) + { + using EllValue = matrix::Ell; + auto ell = as(in); + auto out = Csr::create(ell->get_executor()); + if (std::is_same::value) { + // This cast is necessary for the compiler + as(ell)->convert_to(out.get()); + } + // In-between step required to convert StorageType to value_type + else { + auto tmp = EllValue::create(ell->get_executor()); + ell->convert_to(tmp.get()); + tmp->convert_to(out.get()); + } + return std::move(out); + } + + /** + * Converts a given matrix in CSR format to a matrix in ELL format + * @note ELL format uses StorageType as the storage + * + * @param in Matrix in CSR format + * + * @returns the given matrix converted to the ELL format. + */ + static std::shared_ptr convert_csr_to_ell(const LinOp *in) + { + using CsrStorage = matrix::Csr; + + auto csr = as(in); + auto out = Ell::create(in->get_executor()); + if (std::is_same::value) { + // This cast is necessary for the compiler + as(csr)->convert_to(out.get()); + } + // In-between step required to convert value_type to StorageType + else { + auto tmp = CsrStorage::create(csr->get_executor()); + csr->convert_to(tmp.get()); + tmp->convert_to(out.get()); + } + return std::move(out); + } + /** * Generates the approximate inverse for a triangular matrix and * stores the result in `approximate_inverse_`. @@ -309,7 +330,6 @@ class Isai bool skip_sorting, int power, index_type excess_limit); -private: std::shared_ptr approximate_inverse_; }; From 334b5ca24d7b4b0a743e7c81b101b4b681364d8c Mon Sep 17 00:00:00 2001 From: Fritz Goebel Date: Sun, 6 Jun 2021 13:58:43 +0200 Subject: [PATCH 5/5] Review comments --- include/ginkgo/core/preconditioner/isai.hpp | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/include/ginkgo/core/preconditioner/isai.hpp b/include/ginkgo/core/preconditioner/isai.hpp index 3ebac28cfa6..04e71065736 100644 --- a/include/ginkgo/core/preconditioner/isai.hpp +++ b/include/ginkgo/core/preconditioner/isai.hpp @@ -96,8 +96,9 @@ enum struct isai_type { lower, upper, general, spd }; * @tparam IsaiType determines if the ISAI is generated for a general square * matrix, a lower triangular matrix, an upper triangular matrix or an * spd matrix - * @tparam ValueType precision of matrix elements - * @tparam IndexType precision of matrix indexes + * @tparam ValueType arithmetic precision of matrix elements + * @tparam IndexType precision of matrix indexes + * @tparam StorageType storage precision of matrix elements * * @ingroup isai * @ingroup precond @@ -141,7 +142,6 @@ class Isai Comp, Csr>::type> get_approximate_inverse() const { - auto exec = this->get_executor(); std::shared_ptr return_inv; if (IsaiType == isai_type::spd) { auto ops = as(approximate_inverse_)->get_operators(); @@ -199,15 +199,6 @@ class Isai */ std::shared_ptr GKO_FACTORY_PARAMETER_SCALAR( excess_solver_factory, nullptr); - - /** - * @brief Parameter for reducing preconditioner storage precision. - * - * Reduces the preconditioner storage precision up to reduce_precision - * times. For value types double and complex, storage precision - * can be reduced up to two times, for float and complex once. - */ - int GKO_FACTORY_PARAMETER_SCALAR(reduce_precision, 0); }; GKO_ENABLE_LIN_OP_FACTORY(Isai, parameters, Factory); @@ -233,7 +224,6 @@ class Isai : EnableLinOp(factory->get_executor(), system_matrix->get_size()), parameters_{factory->get_parameters()} { - const auto exec = this->get_executor(); const auto skip_sorting = parameters_.skip_sorting; const auto power = parameters_.sparsity_power; const auto excess_limit = parameters_.excess_limit; @@ -325,6 +315,11 @@ class Isai * * @param skip_sorting dictates if the sorting of the input matrix should * be skipped. + * + * @param power determines which power of the input matrix should be used + * for the sparsity pattern. + * + * @param excess_limit the size limit for the excess system. */ void generate_inverse(std::shared_ptr to_invert, bool skip_sorting, int power,