Add Reference ILU(0)

reference/factorization/ic_kernels.cpp

@@ -33,6 +33,12 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "core/factorization/ic_kernels.hpp"
 
 
+#include <ginkgo/core/base/math.hpp>
+
+
+#include "core/base/allocator.hpp"
+
+
 namespace gko {
 namespace kernels {
 namespace reference {
@@ -46,7 +52,52 @@ namespace ic_factorization {
 
 template <typename ValueType, typename IndexType>
 void compute(std::shared_ptr<const DefaultExecutor> exec,
-             matrix::Csr<ValueType, IndexType> *m) GKO_NOT_IMPLEMENTED;
+             matrix::Csr<ValueType, IndexType> *m)
+{
+    vector<IndexType> diagonals{m->get_size()[0], -1, exec};
+    const auto row_ptrs = m->get_const_row_ptrs();
+    const auto col_idxs = m->get_const_col_idxs();
+    const auto values = m->get_values();
+    for (size_type row = 0; row < m->get_size()[0]; row++) {
+        const auto begin = row_ptrs[row];
+        const auto end = row_ptrs[row + 1];
+        for (auto nz = begin; nz < end; nz++) {
+            const auto col = col_idxs[nz];
+            if (col == row) {
+                diagonals[row] = nz;
+            }
+            if (col > row) {
+                continue;
+            }
+            // accumulate l(row,:) * l(col,:) without the last entry l(col, col)
+            ValueType sum{};
+            auto l_idx = begin;
+            const auto l_end = end;
+            auto lh_idx = row_ptrs[col];
+            const auto lh_end = row_ptrs[col + 1];
+            while (l_idx < l_end && lh_idx < lh_end) {
+                const auto l_col = col_idxs[l_idx];
+                const auto lh_row = col_idxs[lh_idx];
+                // only consider lower triangle of L
+                if (max(l_col, lh_row) > row) {
+                    break;
+                }
+                // ignore l(col, col)
+                if (l_col == lh_row && l_col < col) {
+                    sum += values[l_idx] * conj(values[lh_idx]);
+                }
+                l_idx += l_col <= lh_row ? 1 : 0;
+                lh_idx += lh_row <= l_col ? 1 : 0;
+            }
+            if (row == col) {
+                values[nz] = sqrt(values[nz] - sum);
+            } else {
+                GKO_ASSERT(diagonals[col] != -1);
+                values[nz] = (values[nz] - sum) / values[diagonals[col]];
+            }
+        }
+    }
+}
 
 GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_IC_COMPUTE_KERNEL);
 

reference/factorization/ilu_kernels.cpp

@@ -33,6 +33,15 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "core/factorization/ilu_kernels.hpp"
 
 
+#include <algorithm>
+
+
+#include <ginkgo/core/base/math.hpp>
+
+
+#include "core/base/allocator.hpp"
+
+
 namespace gko {
 namespace kernels {
 namespace reference {
@@ -46,7 +55,46 @@ namespace ilu_factorization {
 
 template <typename ValueType, typename IndexType>
 void compute_lu(std::shared_ptr<const DefaultExecutor> exec,
-                matrix::Csr<ValueType, IndexType> *m) GKO_NOT_IMPLEMENTED;
+                matrix::Csr<ValueType, IndexType> *m)
+{
+    vector<IndexType> diagonals{m->get_size()[0], -1, exec};
+    const auto row_ptrs = m->get_const_row_ptrs();
+    const auto col_idxs = m->get_const_col_idxs();
+    const auto values = m->get_values();
+    for (IndexType row = 0; row < static_cast<IndexType>(m->get_size()[0]);
+         row++) {
+        const auto begin = row_ptrs[row];
+        const auto end = row_ptrs[row + 1];
+        for (auto nz = begin; nz < end; nz++) {
+            const auto col = col_idxs[nz];
+            if (col == row) {
+                diagonals[row] = nz;
+            }
+            auto value = values[nz];
+            for (auto l_nz = begin; l_nz < end; l_nz++) {
+                // for each lower triangular entry l_ik
+                const auto l_col = col_idxs[l_nz];
+                if (l_col >= min(row, col)) {
+                    continue;
+                }
+                // find corresponding entry u_kj
+                const auto u_begin_it = col_idxs + row_ptrs[l_col];
+                const auto u_end_it = col_idxs + row_ptrs[l_col + 1];
+                const auto u_it = std::lower_bound(u_begin_it, u_end_it, col);
+                const auto u_nz = std::distance(col_idxs, u_it);
+                if (u_it != u_end_it && *u_it == col) {
+                    value -= values[l_nz] * values[u_nz];
+                }
+            }
+            if (row <= col) {
+                values[nz] = value;
+            } else {
+                GKO_ASSERT(diagonals[col] != -1);
+                values[nz] = value / values[diagonals[col]];
+            }
+        }
+    }
+}
 
 GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
     GKO_DECLARE_ILU_COMPUTE_LU_KERNEL);

reference/test/factorization/CMakeLists.txt

@@ -1,3 +1,5 @@
+ginkgo_create_test(ic_kernels)
+ginkgo_create_test(ilu_kernels)
 ginkgo_create_test(par_ic_kernels)
 ginkgo_create_test(par_ict_kernels)
 ginkgo_create_test(par_ilu_kernels)

reference/test/factorization/ic_kernels.cpp

@@ -0,0 +1,225 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2021, 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/factorization/ic.hpp>
+
+
+#include <algorithm>
+#include <memory>
+#include <vector>
+
+
+#include <gtest/gtest.h>
+
+
+#include <ginkgo/core/base/executor.hpp>
+#include <ginkgo/core/matrix/coo.hpp>
+#include <ginkgo/core/matrix/csr.hpp>
+#include <ginkgo/core/matrix/dense.hpp>
+
+
+#include "core/test/utils.hpp"
+
+
+namespace {
+
+
+class DummyLinOp : public gko::EnableLinOp<DummyLinOp>,
+                   public gko::EnableCreateMethod<DummyLinOp> {
+public:
+    DummyLinOp(std::shared_ptr<const gko::Executor> exec,
+               gko::dim<2> size = gko::dim<2>{})
+        : EnableLinOp<DummyLinOp>(exec, size)
+    {}
+
+protected:
+    void apply_impl(const gko::LinOp *b, gko::LinOp *x) const override {}
+
+    void apply_impl(const gko::LinOp *alpha, const gko::LinOp *b,
+                    const gko::LinOp *beta, gko::LinOp *x) const override
+    {}
+};
+
+
+template <typename ValueIndexType>
+class Ic : public ::testing::Test {
+protected:
+    using value_type =
+        typename std::tuple_element<0, decltype(ValueIndexType())>::type;
+    using index_type =
+        typename std::tuple_element<1, decltype(ValueIndexType())>::type;
+    using factorization_type = gko::factorization::Ic<value_type, index_type>;
+    using Coo = gko::matrix::Coo<value_type, index_type>;
+    using Csr = gko::matrix::Csr<value_type, index_type>;
+    using Dense = gko::matrix::Dense<value_type>;
+
+    Ic()
+        : ref(gko::ReferenceExecutor::create()),
+          exec(std::static_pointer_cast<const gko::Executor>(ref)),
+          identity(gko::initialize<Csr>(
+              {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}}, ref)),
+          banded(gko::initialize<Csr>(
+              {{1., 1., 0.}, {1., 2., 1.}, {0., 1., 2.}}, ref)),
+          banded_l_expect(gko::initialize<Csr>(
+              {{1., 0., 0.}, {1., 1., 0.}, {0., 1., 1.}}, ref)),
+          mtx_system(gko::initialize<Csr>({{9., 0., -6., 3.},
+                                           {0., 36., 18., 24.},
+                                           {-6., 18., 17., 14.},
+                                           {-3., 24., 14., 18.}},
+                                          ref)),
+          mtx_l_it_expect(gko::initialize<Csr>({{3., 0., 0., 0.},
+                                                {0., 6., 0., 0.},
+                                                {-2., 3., 2., 0.},
+                                                {-1., 4., 0., 1.}},
+                                               ref)),
+          fact_fact(factorization_type::build().on(exec)),
+          tol{r<value_type>::value}
+    {}
+
+    std::shared_ptr<const gko::ReferenceExecutor> ref;
+    std::shared_ptr<const gko::Executor> exec;
+    std::shared_ptr<Csr> identity;
+    std::shared_ptr<Csr> banded;
+    std::shared_ptr<Csr> banded_l_expect;
+    std::shared_ptr<Csr> mtx_system;
+    std::unique_ptr<Csr> mtx_l_it_expect;
+    std::unique_ptr<typename factorization_type::Factory> fact_fact;
+    gko::remove_complex<value_type> tol;
+};
+
+TYPED_TEST_SUITE(Ic, gko::test::ValueIndexTypes);
+
+
+TYPED_TEST(Ic, ThrowNotSupportedForWrongLinOp)
+{
+    auto lin_op = DummyLinOp::create(this->ref);
+
+    ASSERT_THROW(this->fact_fact->generate(gko::share(lin_op)),
+                 gko::NotSupported);
+}
+
+
+TYPED_TEST(Ic, ThrowDimensionMismatch)
+{
+    using Csr = typename TestFixture::Csr;
+    auto matrix = Csr::create(this->ref, gko::dim<2>{2, 3}, 4);
+
+    ASSERT_THROW(this->fact_fact->generate(gko::share(matrix)),
+                 gko::DimensionMismatch);
+}
+
+
+TYPED_TEST(Ic, SetStrategy)
+{
+    using Csr = typename TestFixture::Csr;
+    using factorization_type = typename TestFixture::factorization_type;
+    auto l_strategy = std::make_shared<typename Csr::merge_path>();
+
+    auto factory =
+        factorization_type::build().with_l_strategy(l_strategy).on(this->ref);
+    auto fact = factory->generate(this->mtx_system);
+
+    ASSERT_EQ(factory->get_parameters().l_strategy, l_strategy);
+    ASSERT_EQ(fact->get_l_factor()->get_strategy()->get_name(),
+              l_strategy->get_name());
+    ASSERT_EQ(fact->get_lt_factor()->get_strategy()->get_name(),
+              l_strategy->get_name());
+}
+
+
+TYPED_TEST(Ic, IsConsistentWithComposition)
+{
+    auto fact = this->fact_fact->generate(this->mtx_system);
+
+    auto lin_op_l_factor = gko::as<gko::LinOp>(fact->get_l_factor());
+    auto lin_op_lt_factor = gko::as<gko::LinOp>(fact->get_lt_factor());
+    auto first_operator = fact->get_operators()[0];
+    auto second_operator = fact->get_operators()[1];
+
+    ASSERT_EQ(lin_op_l_factor, first_operator);
+    ASSERT_EQ(lin_op_lt_factor, second_operator);
+}
+
+
+TYPED_TEST(Ic, GenerateIdentity)
+{
+    auto fact = this->fact_fact->generate(this->identity);
+
+    GKO_ASSERT_MTX_NEAR(fact->get_l_factor(), this->identity, this->tol);
+    GKO_ASSERT_MTX_NEAR(fact->get_lt_factor(), this->identity, this->tol);
+}
+
+
+TYPED_TEST(Ic, GenerateDenseIdentity)
+{
+    using Dense = typename TestFixture::Dense;
+    auto dense_id = Dense::create(this->exec, this->identity->get_size());
+    this->identity->convert_to(dense_id.get());
+
+    auto fact = this->fact_fact->generate(gko::share(dense_id));
+
+    GKO_ASSERT_MTX_NEAR(fact->get_l_factor(), this->identity, this->tol);
+    GKO_ASSERT_MTX_NEAR(fact->get_lt_factor(), this->identity, this->tol);
+}
+
+
+TYPED_TEST(Ic, GenerateBanded)
+{
+    using factorization_type = typename TestFixture::factorization_type;
+    using Csr = typename TestFixture::Csr;
+
+    auto fact =
+        factorization_type::build().on(this->exec)->generate(this->banded);
+
+    GKO_ASSERT_MTX_NEAR(fact->get_l_factor(), this->banded_l_expect, this->tol);
+    GKO_ASSERT_MTX_NEAR(fact->get_lt_factor(),
+                        gko::as<Csr>(this->banded_l_expect->conj_transpose()),
+                        this->tol);
+}
+
+
+TYPED_TEST(Ic, GenerateGeneral)
+{
+    using factorization_type = typename TestFixture::factorization_type;
+    using Csr = typename TestFixture::Csr;
+
+    auto fact =
+        factorization_type::build().on(this->exec)->generate(this->mtx_system);
+
+    GKO_ASSERT_MTX_NEAR(fact->get_l_factor(), this->mtx_l_it_expect, this->tol);
+    GKO_ASSERT_MTX_NEAR(fact->get_lt_factor(),
+                        gko::as<Csr>(this->mtx_l_it_expect->conj_transpose()),
+                        this->tol);
+}
+
+
+}  // namespace