GCR solver rewrite

benchmark/solver/solver_common.hpp

@@ -71,8 +71,11 @@ DEFINE_uint32(
     nrhs, 1,
     "The number of right hand sides. Record the residual only when nrhs == 1.");
 
+DEFINE_uint32(gcr_restart, 100,
+              "Maximum dimension of the Krylov space to use in GCR");
+
 DEFINE_uint32(gmres_restart, 100,
-              "What maximum dimension of the Krylov space to use in GMRES");
+              "Maximum dimension of the Krylov space to use in GMRES");
 
 DEFINE_uint32(idr_subspace_dim, 2,
               "What dimension of the subspace to use in IDR");

common/CMakeLists.txt

@@ -23,6 +23,7 @@ set(UNIFIED_SOURCES
     solver/cgs_kernels.cpp
     solver/common_gmres_kernels.cpp
     solver/fcg_kernels.cpp
+    solver/gcr_kernels.cpp
     solver/gmres_kernels.cpp
     solver/ir_kernels.cpp
     )

common/unified/solver/gcr_kernels.cpp

@@ -0,0 +1,126 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2023, 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 "core/solver/gcr_kernels.hpp"
+
+
+#include <ginkgo/core/base/math.hpp>
+
+
+#include "common/unified/base/kernel_launch_solver.hpp"
+
+
+namespace gko {
+namespace kernels {
+namespace GKO_DEVICE_NAMESPACE {
+/**
+ * @brief The GCR solver namespace.
+ *
+ * @ingroup grc
+ */
+namespace gcr {
+
+
+template <typename ValueType>
+void initialize(std::shared_ptr<const DefaultExecutor> exec,
+                const matrix::Dense<ValueType>* b,
+                matrix::Dense<ValueType>* residual,
+                stopping_status* stop_status)
+{
+    run_kernel_solver(
+        exec,
+        [] GKO_KERNEL(auto row, auto col, auto b, auto residual, auto stop) {
+            if (row == 0) {
+                stop[col].reset();
+            }
+            residual(row, col) = b(row, col);
+        },
+        b->get_size(), b->get_stride(), default_stride(b),
+        default_stride(residual), stop_status);
+}
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_GCR_INITIALIZE_KERNEL);
+
+
+template <typename ValueType>
+void restart(std::shared_ptr<const DefaultExecutor> exec,
+             const matrix::Dense<ValueType>* residual,
+             const matrix::Dense<ValueType>* A_residual,
+             matrix::Dense<ValueType>* p_bases,
+             matrix::Dense<ValueType>* Ap_bases, size_type* final_iter_nums)
+{
+    run_kernel_solver(
+        exec,
+        [] GKO_KERNEL(auto row, auto col, auto residual, auto A_residual,
+                      auto p_bases, auto Ap_bases, auto final_iter_nums) {
+            if (row == 0) {
+                final_iter_nums[col] = 0;
+            }
+            p_bases(row, col) = residual(row, col);
+            Ap_bases(row, col) = A_residual(row, col);
+        },
+        residual->get_size(), residual->get_stride(), default_stride(residual),
+        default_stride(A_residual), p_bases, Ap_bases, final_iter_nums);
+}
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_GCR_RESTART_KERNEL);
+
+
+template <typename ValueType>
+void step_1(std::shared_ptr<const DefaultExecutor> exec,
+            matrix::Dense<ValueType>* x, matrix::Dense<ValueType>* residual,
+            const matrix::Dense<ValueType>* p,
+            const matrix::Dense<ValueType>* Ap,
+            const matrix::Dense<remove_complex<ValueType>>* Ap_norm,
+            const matrix::Dense<ValueType>* rAp,
+            const stopping_status* stop_status)
+{
+    run_kernel_solver(
+        exec,
+        [] GKO_KERNEL(auto row, auto col, auto x, auto residual, auto p,
+                      auto Ap, auto Ap_norm, auto rAp, auto stop) {
+            if (!stop[col].has_stopped()) {
+                auto tmp = rAp[col] / Ap_norm[col];
+                x(row, col) += tmp * p(row, col);
+                residual(row, col) -= tmp * Ap(row, col);
+            }
+        },
+        x->get_size(), p->get_stride(), x, residual, p, Ap, Ap_norm, rAp,
+        stop_status);
+}
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_GCR_STEP_1_KERNEL);
+
+}  // namespace gcr
+}  // namespace GKO_DEVICE_NAMESPACE
+}  // namespace kernels
+}  // namespace gko

core/CMakeLists.txt

@@ -62,6 +62,7 @@ target_sources(ginkgo
     solver/cgs.cpp
     solver/direct.cpp
     solver/fcg.cpp
+    solver/gcr.cpp
     solver/gmres.cpp
     solver/idr.cpp
     solver/ir.cpp

core/device_hooks/common_kernels.inc.cpp

@@ -76,6 +76,7 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "core/solver/cgs_kernels.hpp"
 #include "core/solver/common_gmres_kernels.hpp"
 #include "core/solver/fcg_kernels.hpp"
+#include "core/solver/gcr_kernels.hpp"
 #include "core/solver/gmres_kernels.hpp"
 #include "core/solver/idr_kernels.hpp"
 #include "core/solver/ir_kernels.hpp"
@@ -437,6 +438,15 @@ GKO_STUB_VALUE_TYPE(GKO_DECLARE_CGS_STEP_3_KERNEL);
 
 }  // namespace cgs
 
+namespace gcr {
+
+
+GKO_STUB_VALUE_TYPE(GKO_DECLARE_GCR_INITIALIZE_KERNEL);
+GKO_STUB_VALUE_TYPE(GKO_DECLARE_GCR_RESTART_KERNEL);
+GKO_STUB_VALUE_TYPE(GKO_DECLARE_GCR_STEP_1_KERNEL);
+
+
+}  // namespace gcr
 
 namespace common_gmres {
 

core/matrix/dense.cpp

@@ -226,6 +226,15 @@ void Dense<ValueType>::compute_norm1(ptr_param<LinOp> result) const
 }
 
 
+template <typename ValueType>
+void Dense<ValueType>::compute_squared_norm2(LinOp* result) const
+{
+    auto exec = this->get_executor();
+    this->compute_squared_norm2_impl(
+        make_temporary_output_clone(exec, result).get());
+}
+
+
 template <typename ValueType>
 void Dense<ValueType>::inv_scale_impl(const LinOp* alpha)
 {
@@ -470,6 +479,33 @@ void Dense<ValueType>::compute_norm1_impl(LinOp* result) const
 }
 
 
+template <typename ValueType>
+void Dense<ValueType>::compute_squared_norm2(LinOp* result,
+                                             array<char>& tmp) const
+{
+    GKO_ASSERT_EQUAL_DIMENSIONS(result, dim<2>(1, this->get_size()[1]));
+    auto exec = this->get_executor();
+    if (tmp.get_executor() != exec) {
+        tmp.clear();
+        tmp.set_executor(exec);
+    }
+    auto local_result = make_temporary_clone(exec, result);
+    auto dense_res = make_temporary_conversion<remove_complex<ValueType>>(
+        local_result.get());
+    exec->run(dense::make_compute_squared_norm2(this, dense_res.get(), tmp));
+}
+
+
+template <typename ValueType>
+void Dense<ValueType>::compute_squared_norm2_impl(LinOp* result) const
+{
+    auto exec = this->get_executor();
+    array<char> tmp{exec};
+    this->compute_squared_norm2(make_temporary_output_clone(exec, result).get(),
+                                tmp);
+}
+
+
 template <typename ValueType>
 Dense<ValueType>& Dense<ValueType>::operator=(const Dense& other)
 {