ScaledReordered interface

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