add new interface for RCM

core/reorder/rcm.cpp

@@ -66,22 +66,85 @@ GKO_REGISTER_OPERATION(get_degree_of_nodes, rcm::get_degree_of_nodes);
 
 
 template <typename ValueType, typename IndexType>
-void Rcm<ValueType, IndexType>::generate(
-    std::shared_ptr<const Executor>& exec,
-    std::unique_ptr<SparsityMatrix> adjacency_matrix) const
+void rcm_reorder(matrix::SparsityCsr<ValueType, IndexType>* mtx,
+                 IndexType* permutation, IndexType* inv_permutation,
+                 starting_strategy strategy)
 {
-    const IndexType num_rows = adjacency_matrix->get_size()[0];
-    const auto mtx = adjacency_matrix.get();
-    auto degrees = array<IndexType>(exec, num_rows);
-    // RCM is only valid for symmetric matrices. Need to add an expensive check
-    // for symmetricity here ?
+    const auto exec = mtx->get_executor();
+    const IndexType num_rows = mtx->get_size()[0];
+    array<IndexType> degrees{exec, mtx->get_size()[0]};
     exec->run(rcm::make_get_degree_of_nodes(num_rows, mtx->get_const_row_ptrs(),
                                             degrees.get_data()));
     exec->run(rcm::make_get_permutation(
         num_rows, mtx->get_const_row_ptrs(), mtx->get_const_col_idxs(),
-        degrees.get_const_data(), permutation_->get_permutation(),
-        inv_permutation_.get() ? inv_permutation_->get_permutation() : nullptr,
-        parameters_.strategy));
+        degrees.get_const_data(), permutation, inv_permutation, strategy));
+}
+
+
+template <typename ValueType, typename IndexType>
+Rcm<ValueType, IndexType>::Rcm(std::shared_ptr<const Executor> exec)
+    : EnablePolymorphicObject<Rcm, ReorderingBase<IndexType>>(std::move(exec))
+{}
+
+
+template <typename ValueType, typename IndexType>
+Rcm<ValueType, IndexType>::Rcm(const Factory* factory,
+                               const ReorderingBaseArgs& args)
+    : EnablePolymorphicObject<Rcm, ReorderingBase<IndexType>>(
+          factory->get_executor()),
+      parameters_{factory->get_parameters()}
+{
+    // Always execute the reordering on the cpu.
+    const auto is_gpu_executor =
+        this->get_executor() != this->get_executor()->get_master();
+    auto cpu_exec = is_gpu_executor ? this->get_executor()->get_master()
+                                    : this->get_executor();
+
+    auto adjacency_matrix = SparsityMatrix::create(cpu_exec);
+    array<IndexType> degrees;
+
+    // The adjacency matrix has to be square.
+    GKO_ASSERT_IS_SQUARE_MATRIX(args.system_matrix);
+    // This is needed because it does not make sense to call the copy and
+    // convert if the existing matrix is empty.
+    if (args.system_matrix->get_size()) {
+        auto tmp =
+            copy_and_convert_to<SparsityMatrix>(cpu_exec, args.system_matrix);
+        // This function provided within the Sparsity matrix format removes
+        // the diagonal elements and outputs an adjacency matrix.
+        adjacency_matrix = tmp->to_adjacency_matrix();
+    }
+
+    auto const dim = adjacency_matrix->get_size();
+    permutation_ = PermutationMatrix::create(cpu_exec, dim);
+
+    // To make it explicit.
+    inv_permutation_ = nullptr;
+    if (parameters_.construct_inverse_permutation) {
+        inv_permutation_ = PermutationMatrix::create(cpu_exec, dim);
+    }
+
+    rcm_reorder(
+        adjacency_matrix.get(), permutation_->get_permutation(),
+        inv_permutation_ ? inv_permutation_->get_permutation() : nullptr,
+        parameters_.strategy);
+
+    // Copy back results to gpu if necessary.
+    if (is_gpu_executor) {
+        const auto gpu_exec = this->get_executor();
+        auto gpu_perm = share(PermutationMatrix::create(gpu_exec, dim));
+        gpu_perm->copy_from(permutation_);
+        permutation_ = gpu_perm;
+        if (inv_permutation_) {
+            auto gpu_inv_perm = share(PermutationMatrix::create(gpu_exec, dim));
+            gpu_inv_perm->copy_from(inv_permutation_);
+            inv_permutation_ = gpu_inv_perm;
+        }
+    }
+    auto permutation_array =
+        make_array_view(this->get_executor(), permutation_->get_size()[0],
+                        permutation_->get_permutation());
+    this->set_permutation_array(permutation_array);
 }
 
 
@@ -90,4 +153,103 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_RCM);
 
 
 }  // namespace reorder
+
+
+namespace experimental {
+namespace reorder {
+
+
+template <typename IndexType>
+Rcm<IndexType>::Rcm(std::shared_ptr<const Executor> exec,
+                    const parameters_type& params)
+    : EnablePolymorphicObject<Rcm, LinOpFactory>(std::move(exec)),
+      parameters_{params}
+{}
+
+
+template <typename IndexType>
+std::unique_ptr<matrix::Permutation<IndexType>> Rcm<IndexType>::generate(
+    std::shared_ptr<const LinOp> system_matrix) const
+{
+    auto product =
+        std::unique_ptr<permutation_type>(static_cast<permutation_type*>(
+            this->LinOpFactory::generate(std::move(system_matrix)).release()));
+    return product;
+}
+
+
+template <typename IndexType>
+std::unique_ptr<LinOp> Rcm<IndexType>::generate_impl(
+    std::shared_ptr<const LinOp> system_matrix) const
+{
+    GKO_ASSERT_IS_SQUARE_MATRIX(system_matrix);
+    const auto exec = this->get_executor();
+    const auto host_exec = exec->get_master();
+    const auto num_rows = system_matrix->get_size()[0];
+    using complex_scalar = matrix::Dense<std::complex<float>>;
+    using real_scalar = matrix::Dense<float>;
+    using complex_identity = matrix::Identity<std::complex<float>>;
+    using real_identity = matrix::Identity<float>;
+    using complex_mtx = matrix::Csr<std::complex<float>, IndexType>;
+    using real_mtx = matrix::Csr<float, IndexType>;
+    using sparsity_mtx = matrix::SparsityCsr<float, IndexType>;
+    std::unique_ptr<LinOp> converted;
+    // extract row pointers and column indices
+    IndexType* d_row_ptrs{};
+    IndexType* d_col_idxs{};
+    size_type d_nnz{};
+    if (auto convertible = dynamic_cast<const ConvertibleTo<complex_mtx>*>(
+            system_matrix.get())) {
+        auto conv_csr = complex_mtx::create(exec);
+        convertible->convert_to(conv_csr);
+        if (!parameters_.skip_symmetrize) {
+            auto scalar =
+                initialize<complex_scalar>({one<std::complex<float>>()}, exec);
+            auto id = complex_identity::create(exec, conv_csr->get_size()[0]);
+            // compute A^T + A
+            conv_csr->transpose()->apply(scalar, id, scalar, conv_csr);
+        }
+        d_nnz = conv_csr->get_num_stored_elements();
+        d_row_ptrs = conv_csr->get_row_ptrs();
+        d_col_idxs = conv_csr->get_col_idxs();
+        converted = std::move(conv_csr);
+    } else {
+        auto conv_csr = real_mtx::create(exec);
+        as<ConvertibleTo<real_mtx>>(system_matrix)->convert_to(conv_csr);
+        if (!parameters_.skip_symmetrize) {
+            auto scalar = initialize<real_scalar>({one<float>()}, exec);
+            auto id = real_identity::create(exec, conv_csr->get_size()[0]);
+            // compute A^T + A
+            conv_csr->transpose()->apply(scalar, id, scalar, conv_csr);
+        }
+        d_nnz = conv_csr->get_num_stored_elements();
+        d_row_ptrs = conv_csr->get_row_ptrs();
+        d_col_idxs = conv_csr->get_col_idxs();
+        converted = std::move(conv_csr);
+    }
+
+    array<IndexType> permutation(host_exec, num_rows);
+
+    // remove diagonal entries
+    auto pattern =
+        sparsity_mtx::create(exec, gko::dim<2>{num_rows, num_rows},
+                             make_array_view(exec, d_nnz, d_col_idxs),
+                             make_array_view(exec, num_rows + 1, d_row_ptrs));
+    pattern = pattern->to_adjacency_matrix();
+    rcm_reorder(pattern.get(), permutation.get_data(),
+                static_cast<IndexType*>(nullptr), parameters_.strategy);
+
+    // permutation gets copied to device via gko::array constructor
+    return permutation_type::create(exec, dim<2>{num_rows, num_rows},
+                                    std::move(permutation));
+}
+
+
+#undef GKO_DECLARE_RCM
+#define GKO_DECLARE_RCM(IndexType) class Rcm<IndexType>
+GKO_INSTANTIATE_FOR_EACH_INDEX_TYPE(GKO_DECLARE_RCM);
+
+
+}  // namespace reorder
+}  // namespace experimental
 }  // namespace gko

core/test/reorder/rcm.cpp

@@ -52,6 +52,7 @@ class Rcm : public ::testing::Test {
     using v_type = double;
     using i_type = int;
     using reorder_type = gko::reorder::Rcm<v_type, i_type>;
+    using new_reorder_type = gko::experimental::reorder::Rcm<i_type>;
 
     Rcm()
         : exec(gko::ReferenceExecutor::create()),
@@ -67,4 +68,25 @@ TEST_F(Rcm, RcmFactoryKnowsItsExecutor)
     ASSERT_EQ(this->rcm_factory->get_executor(), this->exec);
 }
 
+
+TEST_F(Rcm, NewInterfaceDefaults)
+{
+    auto param = new_reorder_type::build();
+
+    ASSERT_EQ(param.skip_symmetrize, false);
+    ASSERT_EQ(param.strategy,
+              gko::reorder::starting_strategy::pseudo_peripheral);
+}
+
+
+TEST_F(Rcm, NewInterfaceSetParameters)
+{
+    auto param =
+        new_reorder_type::build().with_skip_symmetrize(true).with_strategy(
+            gko::reorder::starting_strategy::minimum_degree);
+
+    ASSERT_EQ(param.skip_symmetrize, false);
+    ASSERT_EQ(param.strategy, gko::reorder::starting_strategy::minimum_degree);
+}
+
 }  // namespace

cuda/test/reorder/rcm_kernels.cpp

@@ -49,34 +49,39 @@ class Rcm : public CudaTestFixture {
     using i_type = int;
     using CsrMtx = gko::matrix::Csr<v_type, i_type>;
     using reorder_type = gko::reorder::Rcm<v_type, i_type>;
+    using new_reorder_type = gko::experimental::reorder::Rcm<i_type>;
     using perm_type = gko::matrix::Permutation<i_type>;
 
 
     Rcm()
-        :  // clang-format off
-          p_mtx(gko::initialize<CsrMtx>({{1.0, 2.0, 0.0, -1.3, 2.1},
+        : p_mtx(gko::initialize<CsrMtx>({{1.0, 2.0, 0.0, -1.3, 2.1},
                                          {2.0, 5.0, 1.5, 0.0, 0.0},
                                          {0.0, 1.5, 1.5, 1.1, 0.0},
                                          {-1.3, 0.0, 1.1, 2.0, 0.0},
                                          {2.1, 0.0, 0.0, 0.0, 1.0}},
-                                        exec)),
-          // clang-format on
-          rcm_factory(reorder_type::build().on(exec)),
-          reorder_op(rcm_factory->generate(p_mtx))
+                                        exec))
     {}
 
-    std::unique_ptr<reorder_type::Factory> rcm_factory;
     std::shared_ptr<CsrMtx> p_mtx;
-    std::unique_ptr<reorder_type> reorder_op;
 };
 
 
-TEST_F(Rcm, IsExecutedOnCpuExecutor)
+TEST_F(Rcm, IsEquivalentToRef)
 {
-    // This only executes successfully if computed on cpu executor.
-    auto p = reorder_op->get_permutation();
+    auto reorder_op = reorder_type::build().on(ref)->generate(p_mtx);
+    auto dreorder_op = reorder_type::build().on(exec)->generate(p_mtx);
 
-    ASSERT_TRUE(true);
+    GKO_ASSERT_ARRAY_EQ(dreorder_op->get_permutation_array(),
+                        reorder_op->get_permutation_array());
+}
+
+
+TEST_F(Rcm, IsEquivalentToRefNewInterface)
+{
+    auto reorder_op = new_reorder_type::build().on(ref)->generate(p_mtx);
+    auto dreorder_op = new_reorder_type::build().on(exec)->generate(p_mtx);
+
+    GKO_ASSERT_MTX_EQ_SPARSITY(dreorder_op, reorder_op);
 }