Merge pull request #1059 from ginkgo-project/reordered_interface

This PR adds an interface for wrapping an inner operator (like a direct or iterative solver) into a (for now symmetrical, e.g. RCM) reordering via a reordering factory and diagonal scaling via one row and one column scaling diagonal matrix.

The ScaledReordered LinOp applies the given scaling and reordering to the system matrix and generates the inner operator on the modified problem. It manages three vectors, a right-hand side and a solution vector for the inner operator, and an intermediate vector used for the permutations.

With a permutation matrix P, a row scaling matrix R, and a column matrix C, the system matrix and vectors are scaled and permuted such that instead of the original system of equations
A*x = b
the inner operator solves
P*R*A*C*P^T*y = P*R*b.
Hence, when generating a ScaledReordered LinOp, the new system matrix P*R*A*C*P^T is computed. When it is applied to vectors b and x, the right-hand side vector is modified to P*R*b, and the solution y to the inner system gives the solution x to the original system via x = C*P^T*y.
When using an inner operator that uses an initial guess, currently the initial guess for the inner system is computed by y_0 = P*C^(-1)*x_0, which is the same as using x_0 for the system P*R*A*x = P*R*b. For the inverse diagonal scaling here I added an inverse_apply to the matrix::Diagonal.

Related PR: #1059

common/cuda_hip/matrix/diagonal_kernels.hpp.inc

@@ -37,7 +37,7 @@ template <typename ValueType, typename IndexType>
 __global__ __launch_bounds__(default_block_size) void apply_to_csr(
     size_type num_rows, const ValueType* __restrict__ diag,
     const IndexType* __restrict__ row_ptrs,
-    ValueType* __restrict__ result_values)
+    ValueType* __restrict__ result_values, bool inverse)
 {
     constexpr auto warp_size = config::warp_size;
     auto warp_tile =
@@ -49,7 +49,7 @@ __global__ __launch_bounds__(default_block_size) void apply_to_csr(
         return;
     }
 
-    const auto diag_val = diag[row];
+    const auto diag_val = inverse ? one<ValueType>() / diag[row] : diag[row];
 
     for (size_type idx = row_ptrs[row] + tid_in_warp; idx < row_ptrs[row + 1];
          idx += warp_size) {

common/unified/matrix/diagonal_kernels.cpp

@@ -54,14 +54,16 @@ template <typename ValueType>
 void apply_to_dense(std::shared_ptr<const DefaultExecutor> exec,
                     const matrix::Diagonal<ValueType>* a,
                     const matrix::Dense<ValueType>* b,
-                    matrix::Dense<ValueType>* c)
+                    matrix::Dense<ValueType>* c, bool inverse)
 {
     run_kernel(
         exec,
-        [] GKO_KERNEL(auto row, auto col, auto diag, auto source, auto result) {
-            result(row, col) = source(row, col) * diag[row];
+        [] GKO_KERNEL(auto row, auto col, auto diag, auto source, auto result,
+                      bool inverse) {
+            result(row, col) = inverse ? source(row, col) / diag[row]
+                                       : source(row, col) * diag[row];
         },
-        b->get_size(), a->get_const_values(), b, c);
+        b->get_size(), a->get_const_values(), b, c, inverse);
 }
 
 GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_DIAGONAL_APPLY_TO_DENSE_KERNEL);

core/CMakeLists.txt

@@ -45,6 +45,7 @@ target_sources(ginkgo
     preconditioner/isai.cpp
     preconditioner/jacobi.cpp
     reorder/rcm.cpp
+    reorder/scaled_reordered.cpp
     solver/bicg.cpp
     solver/bicgstab.cpp
     solver/cb_gmres.cpp

core/matrix/diagonal.cpp

@@ -73,17 +73,19 @@ void Diagonal<ValueType>::apply_impl(const LinOp* b, LinOp* x) const
 
     if (dynamic_cast<const Csr<ValueType, int32>*>(b) &&
         dynamic_cast<Csr<ValueType, int32>*>(x)) {
-        exec->run(diagonal::make_apply_to_csr(
-            this, as<Csr<ValueType, int32>>(b), as<Csr<ValueType, int32>>(x)));
+        exec->run(
+            diagonal::make_apply_to_csr(this, as<Csr<ValueType, int32>>(b),
+                                        as<Csr<ValueType, int32>>(x), false));
     } else if (dynamic_cast<const Csr<ValueType, int64>*>(b) &&
                dynamic_cast<Csr<ValueType, int64>*>(x)) {
-        exec->run(diagonal::make_apply_to_csr(
-            this, as<Csr<ValueType, int64>>(b), as<Csr<ValueType, int64>>(x)));
+        exec->run(
+            diagonal::make_apply_to_csr(this, as<Csr<ValueType, int64>>(b),
+                                        as<Csr<ValueType, int64>>(x), false));
     } else {
         precision_dispatch_real_complex<ValueType>(
             [this, &exec](auto dense_b, auto dense_x) {
-                exec->run(
-                    diagonal::make_apply_to_dense(this, dense_b, dense_x));
+                exec->run(diagonal::make_apply_to_dense(this, dense_b, dense_x,
+                                                        false));
             },
             b, x);
     }
@@ -116,6 +118,32 @@ void Diagonal<ValueType>::rapply_impl(const LinOp* b, LinOp* x) const
 }
 
 
+template <typename ValueType>
+void Diagonal<ValueType>::inverse_apply_impl(const LinOp* b, LinOp* x) const
+{
+    auto exec = this->get_executor();
+
+    if (dynamic_cast<const Csr<ValueType, int32>*>(b) &&
+        dynamic_cast<Csr<ValueType, int32>*>(x)) {
+        exec->run(
+            diagonal::make_apply_to_csr(this, as<Csr<ValueType, int32>>(b),
+                                        as<Csr<ValueType, int32>>(x), true));
+    } else if (dynamic_cast<const Csr<ValueType, int64>*>(b) &&
+               dynamic_cast<Csr<ValueType, int64>*>(x)) {
+        exec->run(
+            diagonal::make_apply_to_csr(this, as<Csr<ValueType, int64>>(b),
+                                        as<Csr<ValueType, int64>>(x), true));
+    } else {
+        precision_dispatch_real_complex<ValueType>(
+            [this, &exec](auto dense_b, auto dense_x) {
+                exec->run(diagonal::make_apply_to_dense(this, dense_b, dense_x,
+                                                        true));
+            },
+            b, x);
+    }
+}
+
+
 template <typename ValueType>
 void Diagonal<ValueType>::apply_impl(const LinOp* alpha, const LinOp* b,
                                      const LinOp* beta, LinOp* x) const

core/matrix/diagonal_kernels.hpp

@@ -53,7 +53,7 @@ namespace kernels {
     void apply_to_dense(std::shared_ptr<const DefaultExecutor> exec, \
                         const matrix::Diagonal<value_type>* a,       \
                         const matrix::Dense<value_type>* b,          \
-                        matrix::Dense<value_type>* c)
+                        matrix::Dense<value_type>* c, bool inverse)
 
 #define GKO_DECLARE_DIAGONAL_RIGHT_APPLY_TO_DENSE_KERNEL(value_type)       \
     void right_apply_to_dense(std::shared_ptr<const DefaultExecutor> exec, \
@@ -65,7 +65,7 @@ namespace kernels {
     void apply_to_csr(std::shared_ptr<const DefaultExecutor> exec,       \
                       const matrix::Diagonal<value_type>* a,             \
                       const matrix::Csr<value_type, index_type>* b,      \
-                      matrix::Csr<value_type, index_type>* c)
+                      matrix::Csr<value_type, index_type>* c, bool inverse)
 
 #define GKO_DECLARE_DIAGONAL_RIGHT_APPLY_TO_CSR_KERNEL(value_type, index_type) \
     void right_apply_to_csr(std::shared_ptr<const DefaultExecutor> exec,       \

core/reorder/scaled_reordered.cpp

@@ -0,0 +1,124 @@
+/*******************************<GINKGO LICENSE>******************************
+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.
+******************************<GINKGO LICENSE>*******************************/
+
+#include <ginkgo/core/reorder/scaled_reordered.hpp>
+
+
+#include <utility>
+
+
+#include <ginkgo/core/base/precision_dispatch.hpp>
+#include <ginkgo/core/matrix/permutation.hpp>
+
+
+namespace gko {
+namespace experimental {
+namespace reorder {
+
+
+template <typename ValueType, typename IndexType>
+void ScaledReordered<ValueType, IndexType>::apply_impl(const LinOp* b,
+                                                       LinOp* x) const
+{
+    precision_dispatch_real_complex<ValueType>(
+        [this](auto dense_b, auto dense_x) {
+            auto exec = this->get_executor();
+            this->set_cache_to(dense_b, dense_x);
+
+            // Preprocess the input vectors before applying the inner operator.
+            if (row_scaling_) {
+                row_scaling_->apply(cache_.inner_b.get(),
+                                    cache_.intermediate.get());
+                std::swap(cache_.inner_b, cache_.intermediate);
+            }
+            // Col scaling for x is only necessary if the inner operator uses an
+            // initial guess. Otherwise x is overwritten anyway.
+            if (col_scaling_ && inner_operator_->apply_uses_initial_guess()) {
+                col_scaling_->inverse_apply(cache_.inner_x.get(),
+                                            cache_.intermediate.get());
+                std::swap(cache_.inner_x, cache_.intermediate);
+            }
+            if (permutation_array_.get_num_elems() > 0) {
+                cache_.inner_b->row_permute(&permutation_array_,
+                                            cache_.intermediate.get());
+                std::swap(cache_.inner_b, cache_.intermediate);
+                if (inner_operator_->apply_uses_initial_guess()) {
+                    cache_.inner_x->row_permute(&permutation_array_,
+                                                cache_.intermediate.get());
+                    std::swap(cache_.inner_x, cache_.intermediate);
+                }
+            }
+
+            inner_operator_->apply(cache_.inner_b.get(), cache_.inner_x.get());
+
+            // Permute and scale the solution vector back.
+            if (permutation_array_.get_num_elems() > 0) {
+                cache_.inner_x->inverse_row_permute(&permutation_array_,
+                                                    cache_.intermediate.get());
+                std::swap(cache_.inner_x, cache_.intermediate);
+            }
+            if (col_scaling_) {
+                col_scaling_->apply(cache_.inner_x.get(),
+                                    cache_.intermediate.get());
+                std::swap(cache_.inner_x, cache_.intermediate);
+            }
+            dense_x->copy_from(cache_.inner_x.get());
+        },
+        b, x);
+}
+
+
+template <typename ValueType, typename IndexType>
+void ScaledReordered<ValueType, IndexType>::apply_impl(const LinOp* alpha,
+                                                       const LinOp* b,
+                                                       const LinOp* beta,
+                                                       LinOp* x) const
+{
+    precision_dispatch_real_complex<ValueType>(
+        [this](auto dense_alpha, auto dense_b, auto dense_beta, auto dense_x) {
+            auto x_clone = dense_x->clone();
+            this->apply_impl(dense_b, x_clone.get());
+            dense_x->scale(dense_beta);
+            dense_x->add_scaled(dense_alpha, x_clone.get());
+        },
+        alpha, b, beta, x);
+}
+
+
+#define GKO_DECLARE_SCALED_REORDERED(ValueType, IndexType) \
+    class ScaledReordered<ValueType, IndexType>
+GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SCALED_REORDERED);
+
+
+}  // namespace reorder
+}  // namespace experimental
+}  // namespace gko

core/test/reorder/CMakeLists.txt

@@ -1 +1,2 @@
 ginkgo_create_test(rcm)
+ginkgo_create_test(scaled_reordered)

core/test/reorder/scaled_reordered.cpp

@@ -0,0 +1,134 @@
+/*******************************<GINKGO LICENSE>******************************
+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.
+******************************<GINKGO LICENSE>*******************************/
+
+#include <ginkgo/core/reorder/scaled_reordered.hpp>
+
+
+#include <memory>
+
+
+#include <gtest/gtest.h>
+
+
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/factorization/par_ic.hpp>
+#include <ginkgo/core/matrix/diagonal.hpp>
+#include <ginkgo/core/reorder/rcm.hpp>
+#include <ginkgo/core/solver/bicgstab.hpp>
+
+
+namespace {
+
+
+class ScaledReorderedFactory : public ::testing::Test {
+protected:
+    using value_type = double;
+    using index_type = gko::int32;
+    using diag = gko::matrix::Diagonal<value_type>;
+    using solver_type = gko::solver::Bicgstab<value_type>;
+    using reorder_type = gko::reorder::Rcm<value_type, index_type>;
+    using scaled_reordered_type =
+        gko::experimental::reorder::ScaledReordered<value_type, index_type>;
+
+    ScaledReorderedFactory()
+        : exec(gko::ReferenceExecutor::create()),
+          scaled_reordered_factory(scaled_reordered_type::build().on(exec)),
+          reordering_factory(reorder_type::build().on(exec)),
+          solver_factory(solver_type::build().on(exec)),
+          diag_matrix(diag::create(exec))
+    {}
+
+    std::shared_ptr<const gko::Executor> exec;
+    std::shared_ptr<diag> diag_matrix;
+    std::shared_ptr<typename scaled_reordered_type::Factory>
+        scaled_reordered_factory;
+    std::shared_ptr<typename reorder_type::Factory> reordering_factory;
+    std::shared_ptr<typename solver_type::Factory> solver_factory;
+};
+
+
+TEST_F(ScaledReorderedFactory, KnowsItsExecutor)
+{
+    auto scaled_reordered_factory =
+        scaled_reordered_type::build().on(this->exec);
+
+    ASSERT_EQ(scaled_reordered_factory->get_executor(), this->exec);
+}
+
+
+TEST_F(ScaledReorderedFactory, CanSetReorderingFactory)
+{
+    auto scaled_reordered_factory =
+        scaled_reordered_type::build()
+            .with_reordering(this->reordering_factory)
+            .on(this->exec);
+
+    ASSERT_EQ(scaled_reordered_factory->get_parameters().reordering,
+              this->reordering_factory);
+}
+
+
+TEST_F(ScaledReorderedFactory, CanSetInnerOperatorFactory)
+{
+    auto scaled_reordered_factory =
+        scaled_reordered_type::build()
+            .with_inner_operator(this->solver_factory)
+            .on(this->exec);
+
+    ASSERT_EQ(scaled_reordered_factory->get_parameters().inner_operator,
+              this->solver_factory);
+}
+
+
+TEST_F(ScaledReorderedFactory, CanSetRowScaling)
+{
+    auto scaled_reordered_factory = scaled_reordered_type::build()
+                                        .with_row_scaling(this->diag_matrix)
+                                        .on(this->exec);
+
+    ASSERT_EQ(scaled_reordered_factory->get_parameters().row_scaling,
+              this->diag_matrix);
+}
+
+
+TEST_F(ScaledReorderedFactory, CanSetColScaling)
+{
+    auto scaled_reordered_factory = scaled_reordered_type::build()
+                                        .with_col_scaling(this->diag_matrix)
+                                        .on(this->exec);
+
+    ASSERT_EQ(scaled_reordered_factory->get_parameters().col_scaling,
+              this->diag_matrix);
+}
+
+
+}  // namespace