Iterative refinement (IR)

core/CMakeLists.txt

@@ -26,6 +26,7 @@ target_sources(ginkgo
         solver/cgs.cpp
         solver/fcg.cpp
         solver/gmres.cpp
+        solver/ir.cpp
         stop/combined.cpp
         stop/criterion.cpp
         stop/iteration.cpp

core/device_hooks/common_kernels.inc.cpp

@@ -46,6 +46,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "core/solver/cgs_kernels.hpp"
 #include "core/solver/fcg_kernels.hpp"
 #include "core/solver/gmres_kernels.hpp"
+#include "core/solver/ir_kernels.hpp"
 #include "core/stop/criterion_kernels.hpp"
 #include "core/stop/residual_norm_reduction_kernels.hpp"
 
@@ -313,6 +314,16 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_GMRES_STEP_2_KERNEL);
 }  // namespace gmres
 
 
+namespace ir {
+
+
+GKO_DECLARE_IR_INITIALIZE_KERNEL
+GKO_NOT_COMPILED(GKO_HOOK_MODULE);
+
+
+}  // namespace ir
+
+
 namespace csr {
 
 

core/solver/ir.cpp

@@ -0,0 +1,121 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright 2017-2019
+
+Karlsruhe Institute of Technology
+Universitat Jaume I
+University of Tennessee
+
+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/solver/ir.hpp>
+
+
+#include <ginkgo/core/matrix/dense.hpp>
+
+
+#include "core/solver/ir_kernels.hpp"
+
+
+namespace gko {
+namespace solver {
+namespace ir {
+
+
+GKO_REGISTER_OPERATION(initialize, ir::initialize);
+
+
+}  // namespace ir
+
+
+template <typename ValueType>
+void Ir<ValueType>::apply_impl(const LinOp *b, LinOp *x) const
+{
+    using Vector = matrix::Dense<ValueType>;
+    constexpr uint8 relative_stopping_id{1};
+
+    auto exec = this->get_executor();
+    auto one_op = initialize<Vector>({one<ValueType>()}, exec);
+    auto neg_one_op = initialize<Vector>({-one<ValueType>()}, exec);
+
+    auto dense_b = as<const Vector>(b);
+    auto dense_x = as<Vector>(x);
+    auto residual = Vector::create_with_config_of(dense_b);
+
+    bool one_changed{};
+    Array<stopping_status> stop_status(exec, dense_b->get_size()[1]);
+    exec->run(ir::make_initialize(&stop_status));
+
+    residual->copy_from(dense_b);
+    system_matrix_->apply(lend(neg_one_op), dense_x, lend(one_op),
+                          lend(residual));
+
+    auto stop_criterion = stop_criterion_factory_->generate(
+        system_matrix_, std::shared_ptr<const LinOp>(b, [](const LinOp *) {}),
+        x, lend(residual));
+
+    int iter = -1;
+    while (true) {
+        ++iter;
+        this->template log<log::Logger::iteration_complete>(
+            this, iter, lend(residual), dense_x);
+
+        if (stop_criterion->update()
+                .num_iterations(iter)
+                .residual(lend(residual))
+                .solution(dense_x)
+                .check(relative_stopping_id, true, &stop_status,
+                       &one_changed)) {
+            break;
+        }
+
+        solver_->apply(lend(one_op), lend(residual), lend(one_op), dense_x);
+        residual->copy_from(dense_b);
+        system_matrix_->apply(lend(neg_one_op), dense_x, lend(one_op),
+                              lend(residual));
+    }
+}
+
+
+template <typename ValueType>
+void Ir<ValueType>::apply_impl(const LinOp *alpha, const LinOp *b,
+                               const LinOp *beta, LinOp *x) const
+{
+    auto dense_x = as<matrix::Dense<ValueType>>(x);
+
+    auto x_clone = dense_x->clone();
+    this->apply(b, x_clone.get());
+    dense_x->scale(beta);
+    dense_x->add_scaled(alpha, x_clone.get());
+}
+
+
+#define GKO_DECLARE_IR(_type) class Ir<_type>
+GKO_INSTANTIATE_FOR_EACH_VALUE_TYPE(GKO_DECLARE_IR);
+
+
+}  // namespace solver
+}  // namespace gko

core/solver/ir_kernels.hpp

@@ -0,0 +1,94 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright 2017-2019
+
+Karlsruhe Institute of Technology
+Universitat Jaume I
+University of Tennessee
+
+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>*******************************/
+
+#ifndef GKO_CORE_SOLVER_IR_KERNELS_HPP_
+#define GKO_CORE_SOLVER_IR_KERNELS_HPP_
+
+
+#include <ginkgo/core/base/array.hpp>
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/base/types.hpp>
+#include <ginkgo/core/stop/stopping_status.hpp>
+
+
+namespace gko {
+namespace kernels {
+namespace ir {
+
+
+#define GKO_DECLARE_IR_INITIALIZE_KERNEL                         \
+    void initialize(std::shared_ptr<const DefaultExecutor> exec, \
+                    Array<stopping_status> *stop_status)
+
+
+#define GKO_DECLARE_ALL_AS_TEMPLATES GKO_DECLARE_IR_INITIALIZE_KERNEL
+
+
+}  // namespace ir
+
+
+namespace omp {
+namespace ir {
+
+GKO_DECLARE_ALL_AS_TEMPLATES;
+
+}  // namespace ir
+}  // namespace omp
+
+
+namespace cuda {
+namespace ir {
+
+GKO_DECLARE_ALL_AS_TEMPLATES;
+
+}  // namespace ir
+}  // namespace cuda
+
+
+namespace reference {
+namespace ir {
+
+GKO_DECLARE_ALL_AS_TEMPLATES;
+
+}  // namespace ir
+}  // namespace reference
+
+
+#undef GKO_DECLARE_ALL_AS_TEMPLATES
+
+
+}  // namespace kernels
+}  // namespace gko
+
+
+#endif  // GKO_CORE_SOLVER_IR_KERNELS_HPP

core/test/log/record.cpp

@@ -478,7 +478,7 @@ TEST(Record, CatchesCriterionCheckCompleted)
         criterion.get(), 1, nullptr, nullptr, nullptr, RelativeStoppingId, true,
         &stop_status, true, true);
 
-    stop_status.get_data()->clear();
+    stop_status.get_data()->reset();
     stop_status.get_data()->stop(RelativeStoppingId);
     auto &data = logger->get().criterion_check_completed.back();
     ASSERT_NE(data->criterion, nullptr);

core/test/log/stream.cpp

@@ -668,7 +668,7 @@ TEST(Stream, CatchesCriterionCheckCompletedWithVerbose)
     std::stringstream true_in_stream;
     true_in_stream << true;
 
-    stop_status.get_data()->clear();
+    stop_status.get_data()->reset();
     stop_status.get_data()->stop(RelativeStoppingId);
     logger->on<gko::log::Logger::criterion_check_completed>(
         criterion.get(), 1, nullptr, nullptr, nullptr, RelativeStoppingId, true,

core/test/solver/CMakeLists.txt

@@ -3,3 +3,4 @@ ginkgo_create_test(cg)
 ginkgo_create_test(cgs)
 ginkgo_create_test(fcg)
 ginkgo_create_test(gmres)
+ginkgo_create_test(ir)