Add a uniform coarse generation algorithm

common/unified/CMakeLists.txt

@@ -18,6 +18,7 @@ set(UNIFIED_SOURCES
     matrix/sparsity_csr_kernels.cpp
     matrix/diagonal_kernels.cpp
     multigrid/pgm_kernels.cpp
+    multigrid/uniform_coarsening_kernels.cpp
     preconditioner/jacobi_kernels.cpp
     solver/bicg_kernels.cpp
     solver/bicgstab_kernels.cpp

common/unified/multigrid/uniform_coarsening_kernels.cpp

@@ -0,0 +1,65 @@
+// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include "core/multigrid/uniform_coarsening_kernels.hpp"
+
+#include <ginkgo/core/base/math.hpp>
+
+#include "common/unified/base/kernel_launch.hpp"
+#include "common/unified/base/kernel_launch_reduction.hpp"
+#include "core/components/prefix_sum_kernels.hpp"
+
+
+namespace gko {
+namespace kernels {
+namespace GKO_DEVICE_NAMESPACE {
+namespace uniform_coarsening {
+
+
+template <typename ValueType, typename IndexType>
+void fill_restrict_op(std::shared_ptr<const DefaultExecutor> exec,
+                      const array<IndexType>* coarse_rows,
+                      matrix::Csr<ValueType, IndexType>* restrict_op)
+{
+    run_kernel(
+        exec,
+        [] GKO_KERNEL(auto tidx, const auto coarse_data,
+                      auto restrict_col_idxs) {
+            if (coarse_data[tidx] >= 0) {
+                restrict_col_idxs[coarse_data[tidx]] = tidx;
+            }
+        },
+        coarse_rows->get_size(), coarse_rows->get_const_data(),
+        restrict_op->get_col_idxs());
+}
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
+    GKO_DECLARE_UNIFORM_COARSENING_FILL_RESTRICT_OP);
+
+
+template <typename IndexType>
+void fill_incremental_indices(std::shared_ptr<const DefaultExecutor> exec,
+                              size_type coarse_skip,
+                              array<IndexType>* coarse_rows)
+{
+    IndexType num_elems = coarse_rows->get_size();
+    run_kernel(
+        exec,
+        [] GKO_KERNEL(auto tidx, auto coarse_skip, auto coarse_data,
+                      auto size) {
+            if (tidx % coarse_skip == 0 && tidx < size) {
+                coarse_data[tidx] = tidx / coarse_skip;
+            }
+        },
+        num_elems, coarse_skip, coarse_rows->get_data(), num_elems);
+}
+
+GKO_INSTANTIATE_FOR_EACH_INDEX_TYPE(
+    GKO_DECLARE_UNIFORM_COARSENING_FILL_INCREMENTAL_INDICES);
+
+
+}  // namespace uniform_coarsening
+}  // namespace GKO_DEVICE_NAMESPACE
+}  // namespace kernels
+}  // namespace gko

core/CMakeLists.txt

@@ -83,6 +83,7 @@ target_sources(${ginkgo_core}
     matrix/sparsity_csr.cpp
     multigrid/pgm.cpp
     multigrid/fixed_coarsening.cpp
+    multigrid/uniform_coarsening.cpp
     preconditioner/batch_jacobi.cpp
     preconditioner/ic.cpp
     preconditioner/ilu.cpp

core/device_hooks/common_kernels.inc.cpp

@@ -45,6 +45,7 @@
 #include "core/matrix/sellp_kernels.hpp"
 #include "core/matrix/sparsity_csr_kernels.hpp"
 #include "core/multigrid/pgm_kernels.hpp"
+#include "core/multigrid/uniform_coarsening_kernels.hpp"
 #include "core/preconditioner/batch_jacobi_kernels.hpp"
 #include "core/preconditioner/isai_kernels.hpp"
 #include "core/preconditioner/jacobi_kernels.hpp"
@@ -963,6 +964,16 @@ GKO_STUB_INDEX_TYPE(GKO_DECLARE_PGM_GATHER_INDEX);
 }  // namespace pgm
 
 
+namespace uniform_coarsening {
+
+
+GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UNIFORM_COARSENING_FILL_RESTRICT_OP);
+GKO_STUB_INDEX_TYPE(GKO_DECLARE_UNIFORM_COARSENING_FILL_INCREMENTAL_INDICES);
+
+
+}  // namespace uniform_coarsening
+
+
 namespace set_all_statuses {
 
 

core/multigrid/uniform_coarsening.cpp

@@ -0,0 +1,106 @@
+// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <ginkgo/core/base/array.hpp>
+#include <ginkgo/core/base/exception_helpers.hpp>
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/base/polymorphic_object.hpp>
+#include <ginkgo/core/base/types.hpp>
+#include <ginkgo/core/base/utils.hpp>
+#include <ginkgo/core/matrix/csr.hpp>
+#include <ginkgo/core/matrix/dense.hpp>
+#include <ginkgo/core/matrix/identity.hpp>
+#include <ginkgo/core/matrix/row_gatherer.hpp>
+#include <ginkgo/core/matrix/sparsity_csr.hpp>
+#include <ginkgo/core/multigrid/uniform_coarsening.hpp>
+
+#include "core/base/utils.hpp"
+#include "core/components/fill_array_kernels.hpp"
+#include "core/matrix/csr_builder.hpp"
+#include "core/multigrid/uniform_coarsening_kernels.hpp"
+#include "ginkgo/core/base/math.hpp"
+
+
+namespace gko {
+namespace multigrid {
+namespace uniform_coarsening {
+namespace {
+
+
+GKO_REGISTER_OPERATION(fill_restrict_op, uniform_coarsening::fill_restrict_op);
+GKO_REGISTER_OPERATION(fill_incremental_indices,
+                       uniform_coarsening::fill_incremental_indices);
+GKO_REGISTER_OPERATION(fill_array, components::fill_array);
+GKO_REGISTER_OPERATION(fill_seq_array, components::fill_seq_array);
+
+
+}  // anonymous namespace
+}  // namespace uniform_coarsening
+
+
+template <typename ValueType, typename IndexType>
+void UniformCoarsening<ValueType, IndexType>::generate()
+{
+    using csr_type = matrix::Csr<ValueType, IndexType>;
+    using sparsity_type = matrix::SparsityCsr<ValueType, IndexType>;
+    using real_type = remove_complex<ValueType>;
+    using weight_csr_type = remove_complex<csr_type>;
+    auto exec = this->get_executor();
+    const auto num_rows = this->system_matrix_->get_size()[0];
+
+    // Only support csr matrix currently.
+    const csr_type* uniform_coarsening_op =
+        dynamic_cast<const csr_type*>(system_matrix_.get());
+    std::shared_ptr<const csr_type> uniform_coarsening_op_shared_ptr{};
+    // If system matrix is not csr or need sorting, generate the csr.
+    if (!parameters_.skip_sorting || !uniform_coarsening_op) {
+        uniform_coarsening_op_shared_ptr = convert_to_with_sorting<csr_type>(
+            exec, system_matrix_, parameters_.skip_sorting);
+        uniform_coarsening_op = uniform_coarsening_op_shared_ptr.get();
+        // keep the same precision data in fine_op
+        this->set_fine_op(uniform_coarsening_op_shared_ptr);
+    }
+    // Use -1 as sentinel value
+    coarse_rows_ = array<IndexType>(exec, num_rows);
+    coarse_rows_.fill(-one<IndexType>());
+
+    // Fill with incremental local indices.
+    exec->run(uniform_coarsening::make_fill_incremental_indices(
+        parameters_.coarse_skip, &coarse_rows_));
+
+    gko::dim<2>::dimension_type coarse_dim =
+        gko::ceildiv(coarse_rows_.get_size(), parameters_.coarse_skip);
+    auto fine_dim = system_matrix_->get_size()[0];
+    auto restrict_op = share(
+        csr_type::create(exec, gko::dim<2>{coarse_dim, fine_dim}, coarse_dim,
+                         uniform_coarsening_op->get_strategy()));
+    exec->run(uniform_coarsening::make_fill_restrict_op(&coarse_rows_,
+                                                        restrict_op.get()));
+    exec->run(uniform_coarsening::make_fill_array(
+        restrict_op->get_values(), coarse_dim, one<ValueType>()));
+    exec->run(uniform_coarsening::make_fill_seq_array(
+        restrict_op->get_row_ptrs(), coarse_dim + 1));
+
+    auto prolong_op = gko::as<csr_type>(share(restrict_op->transpose()));
+
+    // TODO: Can be done with submatrix index_set.
+    auto coarse_matrix =
+        share(csr_type::create(exec, gko::dim<2>{coarse_dim, coarse_dim}));
+    coarse_matrix->set_strategy(uniform_coarsening_op->get_strategy());
+    auto tmp = csr_type::create(exec, gko::dim<2>{fine_dim, coarse_dim});
+    tmp->set_strategy(uniform_coarsening_op->get_strategy());
+    uniform_coarsening_op->apply(prolong_op, tmp);
+    restrict_op->apply(tmp, coarse_matrix);
+
+    this->set_multigrid_level(prolong_op, coarse_matrix, restrict_op);
+}
+
+
+#define GKO_DECLARE_UNIFORM_COARSENING(_vtype, _itype) \
+    class UniformCoarsening<_vtype, _itype>
+GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_UNIFORM_COARSENING);
+
+
+}  // namespace multigrid
+}  // namespace gko

core/multigrid/uniform_coarsening_kernels.hpp

@@ -0,0 +1,55 @@
+// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#ifndef GKO_CORE_MULTIGRID_UNIFORM_COARSENING_KERNELS_HPP_
+#define GKO_CORE_MULTIGRID_UNIFORM_COARSENING_KERNELS_HPP_
+
+
+#include <memory>
+
+#include <ginkgo/core/base/array.hpp>
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/matrix/csr.hpp>
+#include <ginkgo/core/multigrid/uniform_coarsening.hpp>
+
+#include "core/base/kernel_declaration.hpp"
+
+
+namespace gko {
+namespace kernels {
+namespace uniform_coarsening {
+
+#define GKO_DECLARE_UNIFORM_COARSENING_FILL_RESTRICT_OP(ValueType, IndexType) \
+    void fill_restrict_op(std::shared_ptr<const DefaultExecutor> exec,        \
+                          const array<IndexType>* coarse_rows,                \
+                          matrix::Csr<ValueType, IndexType>* restrict_op)
+
+#define GKO_DECLARE_UNIFORM_COARSENING_FILL_INCREMENTAL_INDICES(IndexType)     \
+    void fill_incremental_indices(std::shared_ptr<const DefaultExecutor> exec, \
+                                  size_type coarse_skip,                       \
+                                  array<IndexType>* coarse_rows)
+
+
+#define GKO_DECLARE_ALL_AS_TEMPLATES                                       \
+    template <typename ValueType, typename IndexType>                      \
+    GKO_DECLARE_UNIFORM_COARSENING_FILL_RESTRICT_OP(ValueType, IndexType); \
+    template <typename IndexType>                                          \
+    GKO_DECLARE_UNIFORM_COARSENING_FILL_INCREMENTAL_INDICES(IndexType)
+
+
+}  // namespace uniform_coarsening
+
+
+GKO_DECLARE_FOR_ALL_EXECUTOR_NAMESPACES(uniform_coarsening,
+                                        GKO_DECLARE_ALL_AS_TEMPLATES);
+
+
+#undef GKO_DECLARE_ALL_AS_TEMPLATES
+
+
+}  // namespace kernels
+}  // namespace gko
+
+
+#endif  // GKO_CORE_MULTIGRID_UNIFORM_COARSENING_KERNELS_HPP_

core/test/multigrid/CMakeLists.txt

@@ -1,2 +1,3 @@
 ginkgo_create_test(pgm)
 ginkgo_create_test(fixed_coarsening)
+ginkgo_create_test(uniform_coarsening)

core/test/multigrid/uniform_coarsening.cpp

@@ -0,0 +1,73 @@
+// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <memory>
+
+#include <gtest/gtest.h>
+
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/multigrid/uniform_coarsening.hpp>
+
+#include "core/test/utils.hpp"
+
+
+namespace {
+
+
+template <typename ValueIndexType>
+class UniformCoarseningFactory : 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<value_type, index_type>;
+    using Vec = gko::matrix::Dense<value_type>;
+    using MgLevel = gko::multigrid::UniformCoarsening<value_type, index_type>;
+    UniformCoarseningFactory()
+        : exec(gko::ReferenceExecutor::create()),
+          uniform_coarsening1_factory(
+              MgLevel::build().with_coarse_skip(4).with_skip_sorting(true).on(
+                  exec))
+    {}
+
+    std::shared_ptr<const gko::Executor> exec;
+    std::unique_ptr<typename MgLevel::Factory> uniform_coarsening1_factory;
+};
+
+TYPED_TEST_SUITE(UniformCoarseningFactory, gko::test::ValueIndexTypes,
+                 PairTypenameNameGenerator);
+
+
+TYPED_TEST(UniformCoarseningFactory, FactoryKnowsItsExecutor)
+{
+    ASSERT_EQ(this->uniform_coarsening1_factory->get_executor(), this->exec);
+}
+
+
+TYPED_TEST(UniformCoarseningFactory, DefaultSetting)
+{
+    using MgLevel = typename TestFixture::MgLevel;
+    auto factory = MgLevel::build().on(this->exec);
+
+    ASSERT_EQ(factory->get_parameters().coarse_skip, 2);
+    ASSERT_EQ(factory->get_parameters().skip_sorting, false);
+}
+
+
+TYPED_TEST(UniformCoarseningFactory, SetNumJumps)
+{
+    ASSERT_EQ(this->uniform_coarsening1_factory->get_parameters().coarse_skip,
+              4);
+}
+
+
+TYPED_TEST(UniformCoarseningFactory, SetSkipSorting)
+{
+    ASSERT_EQ(this->uniform_coarsening1_factory->get_parameters().skip_sorting,
+              true);
+}
+
+
+}  // namespace

examples/multigrid-preconditioned-solver/doc/results.dox

@@ -3,18 +3,19 @@ This is the expected output:
 
 @code{.cpp}
 
+Problem size: (8000, 8000)
 Initial residual norm sqrt(r^T r):
 %%MatrixMarket matrix array real general
 1 1
-4.3589
+163.578
 Final residual norm sqrt(r^T r):
 %%MatrixMarket matrix array real general
 1 1
-1.69858e-09
-CG iteration count:     39
-CG generation time [ms]: 2.04293
-CG execution time [ms]: 22.3874
-CG execution time per iteration[ms]: 0.574036
+7.61038e-09
+CG iteration count:     1
+CG generation time [ms]: 14.4458
+CG execution time [ms]: 10.975
+CG execution time per iteration[ms]: 10.975
 
 @endcode