diff --git a/core/factorization/factorization.cpp b/core/factorization/factorization.cpp index 1df1f49aa13..4a16e84a1e0 100644 --- a/core/factorization/factorization.cpp +++ b/core/factorization/factorization.cpp @@ -31,10 +31,25 @@ GKO_REGISTER_OPERATION(initialize_l, factorization::initialize_l); template std::unique_ptr> -Factorization::unpack() const +Factorization::unpack( + std::shared_ptr lower_factor_strategy, + std::shared_ptr upper_factor_strategy) + const { const auto exec = this->get_executor(); const auto size = this->get_size(); + auto create_matrix = [](auto exec, auto size, auto vals, auto col_idxs, + auto row_ptrs, auto strategy) { + if (strategy == nullptr) { + return matrix_type::create(exec, size, std::move(vals), + std::move(col_idxs), + std::move(row_ptrs)); + } else { + return matrix_type::create(exec, size, std::move(vals), + std::move(col_idxs), std::move(row_ptrs), + strategy); + } + }; switch (this->get_storage_type()) { case storage_type::empty: GKO_NOT_SUPPORTED(nullptr); @@ -53,12 +68,14 @@ Factorization::unpack() const const auto u_nnz = static_cast(get_element(u_row_ptrs, size[0])); // create matrices - auto l_mtx = matrix_type::create( - exec, size, array{exec, l_nnz}, - array{exec, l_nnz}, std::move(l_row_ptrs)); - auto u_mtx = matrix_type::create( - exec, size, array{exec, u_nnz}, - array{exec, u_nnz}, std::move(u_row_ptrs)); + auto l_mtx = + create_matrix(exec, size, array{exec, l_nnz}, + array{exec, l_nnz}, std::move(l_row_ptrs), + lower_factor_strategy); + auto u_mtx = + create_matrix(exec, size, array{exec, u_nnz}, + array{exec, u_nnz}, std::move(u_row_ptrs), + upper_factor_strategy); // fill matrices exec->run(make_initialize_l_u(mtx.get(), l_mtx.get(), u_mtx.get())); return create_from_composition( @@ -72,9 +89,10 @@ Factorization::unpack() const const auto l_nnz = static_cast(get_element(l_row_ptrs, size[0])); // create matrices - auto l_mtx = matrix_type::create( - exec, size, array{exec, l_nnz}, - array{exec, l_nnz}, std::move(l_row_ptrs)); + auto l_mtx = + create_matrix(exec, size, array{exec, l_nnz}, + array{exec, l_nnz}, std::move(l_row_ptrs), + lower_factor_strategy); // fill matrices exec->run(make_initialize_l(mtx.get(), l_mtx.get(), false)); auto u_mtx = l_mtx->conj_transpose(); diff --git a/core/factorization/ilu.cpp b/core/factorization/ilu.cpp index dccfabcb51a..32812f699c1 100644 --- a/core/factorization/ilu.cpp +++ b/core/factorization/ilu.cpp @@ -80,7 +80,6 @@ std::unique_ptr> Ilu::generate_l_u( // Converts the system matrix to CSR. // Throws an exception if it is not convertible. auto local_system_matrix = share(matrix_type::create(exec)); - std::shared_ptr ilu; as>(system_matrix.get()) ->convert_to(local_system_matrix); @@ -104,25 +103,26 @@ std::unique_ptr> Ilu::generate_l_u( make_const_array_view( exec, local_system_matrix->get_size()[0] + 1, local_system_matrix->get_const_row_ptrs()))); - ilu = + auto unpack = gko::experimental::factorization::Lu::build() .with_has_all_fillin(false) .with_symbolic_factorization(sparsity) .on(exec) ->generate(local_system_matrix) - ->get_combined(); - } else { - exec->run( - ilu_factorization::make_compute_ilu(local_system_matrix.get())); - ilu = local_system_matrix; + ->unpack(parameters_.l_strategy, parameters_.u_strategy); + return Composition::create(unpack->get_lower_factor(), + unpack->get_upper_factor()); } + exec->run(ilu_factorization::make_compute_ilu(local_system_matrix.get())); + // Separate L and U factors: nnz - const auto matrix_size = ilu->get_size(); + const auto matrix_size = local_system_matrix->get_size(); const auto num_rows = matrix_size[0]; array l_row_ptrs{exec, num_rows + 1}; array u_row_ptrs{exec, num_rows + 1}; exec->run(ilu_factorization::make_initialize_row_ptrs_l_u( - ilu.get(), l_row_ptrs.get_data(), u_row_ptrs.get_data())); + local_system_matrix.get(), l_row_ptrs.get_data(), + u_row_ptrs.get_data())); // Get nnz from device memory auto l_nnz = static_cast(get_element(l_row_ptrs, num_rows)); @@ -141,8 +141,8 @@ std::unique_ptr> Ilu::generate_l_u( std::move(u_row_ptrs), parameters_.u_strategy); // Separate L and U: columns and values - exec->run(ilu_factorization::make_initialize_l_u(ilu.get(), l_factor.get(), - u_factor.get())); + exec->run(ilu_factorization::make_initialize_l_u( + local_system_matrix.get(), l_factor.get(), u_factor.get())); return Composition::create(std::move(l_factor), std::move(u_factor)); diff --git a/include/ginkgo/core/factorization/factorization.hpp b/include/ginkgo/core/factorization/factorization.hpp index 39345f59a44..5b203fb7e09 100644 --- a/include/ginkgo/core/factorization/factorization.hpp +++ b/include/ginkgo/core/factorization/factorization.hpp @@ -88,10 +88,21 @@ class Factorization : public EnableLinOp> { * for triangular solves to a composition representation that can also be * used to access individual factors and multiply with the factorization. * + * @param lower_factor_strategy the Csr strategy for the lower factor and + * the transposed lower factor. + * @param upper_factor_strategy the Csr strategy for the upper factor + * * @return a new Factorization object containing this factorization * represented as storage_type::composition. + * + * @note The strategy only has effect when it is unpacked from the combined + * matrix. */ - std::unique_ptr unpack() const; + std::unique_ptr unpack( + std::shared_ptr + lower_factor_strategy = nullptr, + std::shared_ptr + upper_factor_strategy = nullptr) const; /** Returns the storage type used by this factorization. */ storage_type get_storage_type() const; diff --git a/reference/test/factorization/factorization.cpp b/reference/test/factorization/factorization.cpp index 2ded81d4867..402307e452a 100644 --- a/reference/test/factorization/factorization.cpp +++ b/reference/test/factorization/factorization.cpp @@ -252,6 +252,31 @@ TYPED_TEST(Factorization, UnpackCombinedLUWorks) } +TYPED_TEST(Factorization, UnpackCombinedLUWorksWithStrategy) +{ + using factorization_type = typename TestFixture::factorization_type; + using matrix_type = typename TestFixture::matrix_type; + auto fact = factorization_type::create_from_combined_lu( + this->combined_mtx->clone()); + + auto separated = + fact->unpack(std::make_shared(), + std::make_shared()); + + ASSERT_EQ(separated->get_storage_type(), + gko::experimental::factorization::storage_type::composition); + ASSERT_EQ(separated->get_combined(), nullptr); + ASSERT_EQ(separated->get_diagonal(), nullptr); + GKO_ASSERT_MTX_NEAR(separated->get_lower_factor(), this->lower_mtx, 0.0); + GKO_ASSERT_MTX_NEAR(separated->get_upper_factor(), this->upper_nonunit_mtx, + 0.0); + ASSERT_EQ(separated->get_lower_factor()->get_strategy()->get_name(), + "classical"); + ASSERT_EQ(separated->get_upper_factor()->get_strategy()->get_name(), + "merge_path"); +} + + TYPED_TEST(Factorization, UnpackSymmCombinedCholeskyWorks) { using matrix_type = typename TestFixture::matrix_type; @@ -273,6 +298,35 @@ TYPED_TEST(Factorization, UnpackSymmCombinedCholeskyWorks) } +TYPED_TEST(Factorization, UnpackSymmCombinedCholeskyWorksWithStrategy) +{ + using matrix_type = typename TestFixture::matrix_type; + using factorization_type = typename TestFixture::factorization_type; + auto fact = factorization_type::create_from_combined_cholesky( + this->combined_mtx->clone()); + + // second one is ignored in cholesky to keep the same behavior as + // factorization::Ic + auto separated = + fact->unpack(std::make_shared(), + std::make_shared()); + + ASSERT_EQ(separated->get_storage_type(), + gko::experimental::factorization::storage_type::symm_composition); + ASSERT_EQ(separated->get_combined(), nullptr); + ASSERT_EQ(separated->get_diagonal(), nullptr); + GKO_ASSERT_MTX_NEAR(separated->get_lower_factor(), this->lower_cholesky_mtx, + 0.0); + GKO_ASSERT_MTX_NEAR( + separated->get_upper_factor(), + gko::as(this->lower_cholesky_mtx->conj_transpose()), 0.0); + ASSERT_EQ(separated->get_lower_factor()->get_strategy()->get_name(), + "classical"); + ASSERT_EQ(separated->get_upper_factor()->get_strategy()->get_name(), + "classical"); +} + + TYPED_TEST(Factorization, UnpackCompositionWorks) { using factorization_type = typename TestFixture::factorization_type;