Add MC64

.github/workflows/windows-msvc-ref.yml

@@ -27,8 +27,10 @@ jobs:
       fail-fast: false
       matrix:
         config:
-        - {shared: "ON", build_type: "Debug", name: "reference/debug/shared"}
+        # Debug shared exceeds symbol limit
+        # - {shared: "ON", build_type: "Debug", name: "reference/debug/shared"}
         - {shared: "OFF", build_type: "Release", name: "reference/release/static"}
+        - {shared: "ON", build_type: "Release", name: "reference/release/shared"}
         # Debug static needs too much storage
         # - {shared: "OFF", build_type: "Debug", name: "reference/debug/static"}
     name: msvc/${{ matrix.config.name }}

core/reorder/mc64.cpp

@@ -580,25 +580,29 @@ std::unique_ptr<LinOp> Mc64<ValueType, IndexType>::generate_impl(
     const auto row_ptrs = mtx->get_const_row_ptrs();
     const auto col_idxs = mtx->get_const_col_idxs();
 
-    exec->run(make_initialize_weights(mtx.get(), weights, dual_u, row_maxima,
-                                      parameters_.strategy));
-
-    // Compute an initial maximum matching from the nonzero entries for which
-    // the reduced weight (W(i, j) - u(j) - v(i)) is zero. Here, W is the
-    // weight matrix and u and v are the dual vectors. Note that v initially
-    // only contains zeros and hence can still be ignored here.
-    exec->run(make_initial_matching(
-        num_rows, row_ptrs, col_idxs, weights, dual_u, permutation,
-        inv_permutation, matched_idxs, unmatched_rows, parameters_.tolerance));
-
-    exec->run(make_augment_matching(
-        mtx.get(), weights, dual_u, distance, permutation, inv_permutation,
-        unmatched_rows, parents, generation, marked_cols, matched_idxs,
-        this->get_parameters().tolerance));
-
-    exec->run(make_compute_scaling(
-        mtx.get(), weights, dual_u, row_maxima, permutation, matched_idxs,
-        parameters_.strategy, row_scaling.get_data(), col_scaling.get_data()));
+    if (num_rows > 0) {
+        exec->run(make_initialize_weights(mtx.get(), weights, dual_u,
+                                          row_maxima, parameters_.strategy));
+
+        // Compute an initial maximum matching from the nonzero entries for
+        // which the reduced weight (W(i, j) - u(j) - v(i)) is zero. Here, W is
+        // the weight matrix and u and v are the dual vectors. Note that v
+        // initially only contains zeros and hence can still be ignored here.
+        exec->run(make_initial_matching(num_rows, row_ptrs, col_idxs, weights,
+                                        dual_u, permutation, inv_permutation,
+                                        matched_idxs, unmatched_rows,
+                                        parameters_.tolerance));
+
+        exec->run(make_augment_matching(
+            mtx.get(), weights, dual_u, distance, permutation, inv_permutation,
+            unmatched_rows, parents, generation, marked_cols, matched_idxs,
+            this->get_parameters().tolerance));
+
+        exec->run(make_compute_scaling(
+            mtx.get(), weights, dual_u, row_maxima, permutation, matched_idxs,
+            parameters_.strategy, row_scaling.get_data(),
+            col_scaling.get_data()));
+    }
 
     array<index_type> identity_permutation{exec, num_rows};
     exec->run(make_fill_seq_array(identity_permutation.get_data(), num_rows));

reference/test/reorder/mc64.cpp

@@ -134,6 +134,22 @@ TYPED_TEST(Mc64, CanBeCreatedWithReorderingStrategy)
 }
 
 
+TYPED_TEST(Mc64, CanBeCreatedEmpty)
+{
+    using reorder_type = typename TestFixture::reorder_type;
+    using matrix_type = typename TestFixture::CsrMtx;
+
+    auto mc64 =
+        reorder_type::build()
+            .with_strategy(
+                gko::experimental::reorder::mc64_strategy::max_diagonal_sum)
+            .on(this->exec)
+            ->generate(matrix_type::create(this->exec));
+
+    ASSERT_FALSE(mc64->get_size());
+}
+
+
 TYPED_TEST(Mc64, CanBeCreatedWithTolerance)
 {
     using reorder_type = typename TestFixture::reorder_type;

reference/test/reorder/mc64_kernels.cpp

@@ -68,8 +68,12 @@ class Mc64 : public ::testing::Test {
     using real_type = gko::remove_complex<value_type>;
     using matrix_type = gko::matrix::Csr<value_type, index_type>;
     using perm_type = gko::matrix::ScaledPermutation<value_type, index_type>;
-    static constexpr auto inf = std::numeric_limits<real_type>::infinity();
-    static constexpr real_type tol = 1e-14;
+    // this is a constexpr member functions to avoid having to export the symbol
+    // for the constant variable
+    static constexpr real_type inf()
+    {
+        return std::numeric_limits<real_type>::infinity();
+    }
 
     Mc64()
         : ref(gko::ReferenceExecutor::create()),
@@ -122,7 +126,8 @@ class Mc64 : public ::testing::Test {
                                 static_cast<real_type>(std::log2(4.)),
                                 static_cast<real_type>(std::log2(4.)),
                                 static_cast<real_type>(std::log2(8.))}},
-          initialized_distance{ref, I<real_type>{inf, inf, inf, inf, inf, inf}},
+          initialized_distance{
+              ref, I<real_type>{inf(), inf(), inf(), inf(), inf(), inf()}},
           empty_permutation{ref, I<index_type>{-1, -1, -1, -1, -1, -1}},
           empty_inverse_permutation{ref, I<index_type>{-1, -1, -1, -1, -1, -1}},
           empty_matched_idxs{ref, I<index_type>{0, 0, 0, 0, 0, 0}},
@@ -144,8 +149,8 @@ class Mc64 : public ::testing::Test {
           final_weights{ref, I<real_type>{2., 1., 0., 0., 4., 0., 2., 0., 1.,
                                           0., 2., 3., 0.}},
           final_dual_u{ref, I<real_type>{0., 1., -1., -2., 0., 0.}},
-          final_distance{ref, I<real_type>{inf, inf, 1., 0., inf, 1.}},
-          tolerance{10 * std::numeric_limits<real_type>::epsilon()}
+          final_distance{ref, I<real_type>{inf(), inf(), 1., 0., inf(), 1.}},
+          zero_tol{1e-14}
     {}
 
     std::pair<std::shared_ptr<const perm_type>,
@@ -156,6 +161,20 @@ class Mc64 : public ::testing::Test {
                               gko::as<perm_type>(result->get_operators()[1]));
     }
 
+    void assert_array_near(const gko::array<real_type>& a,
+                           const gko::array<real_type>& b, std::string name)
+    {
+        ASSERT_EQ(a.get_num_elems(), b.get_num_elems());
+        for (gko::size_type i = 0; i < a.get_num_elems(); i++) {
+            if (std::isfinite(a.get_const_data()[i]) ||
+                std::isfinite(b.get_const_data()[i])) {
+                ASSERT_NEAR(a.get_const_data()[i], b.get_const_data()[i],
+                            r<value_type>::value)
+                    << name << '[' << i << ']';
+            }
+        }
+    }
+
     std::shared_ptr<const gko::ReferenceExecutor> ref;
     gko::array<real_type> tmp;
     gko::array<real_type> weights;
@@ -190,7 +209,7 @@ class Mc64 : public ::testing::Test {
     gko::array<index_type> final_marked_cols;
     gko::array<index_type> final_matched_idxs;
     std::shared_ptr<matrix_type> mtx;
-    const real_type tolerance;
+    const real_type zero_tol;
 };
 
 TYPED_TEST_SUITE(Mc64, gko::test::ValueIndexTypes, PairTypenameNameGenerator);
@@ -205,9 +224,12 @@ TYPED_TEST(Mc64, InitializeWeightsSum)
         this->mtx.get(), this->weights, this->dual_u, this->row_maxima,
         gko::experimental::reorder::mc64_strategy::max_diagonal_sum);
 
-    GKO_ASSERT_ARRAY_EQ(this->weights, this->initialized_weights_sum);
-    GKO_ASSERT_ARRAY_EQ(this->dual_u, this->initialized_dual_u_sum);
-    GKO_ASSERT_ARRAY_EQ(this->row_maxima, this->initialized_row_maxima_sum);
+    this->assert_array_near(this->weights, this->initialized_weights_sum,
+                            "weights");
+    this->assert_array_near(this->dual_u, this->initialized_dual_u_sum,
+                            "dual_u");
+    this->assert_array_near(this->row_maxima, this->initialized_row_maxima_sum,
+                            "row_maxima");
 }
 
 
@@ -220,9 +242,12 @@ TYPED_TEST(Mc64, InitializeWeightsProduct)
         this->mtx.get(), this->weights, this->dual_u, this->row_maxima,
         gko::experimental::reorder::mc64_strategy::max_diagonal_product);
 
-    GKO_ASSERT_ARRAY_EQ(this->weights, this->initialized_weights_product);
-    GKO_ASSERT_ARRAY_EQ(this->dual_u, this->initialized_dual_u_product);
-    GKO_ASSERT_ARRAY_EQ(this->row_maxima, this->initialized_row_maxima_product);
+    this->assert_array_near(this->weights, this->initialized_weights_product,
+                            "weights");
+    this->assert_array_near(this->dual_u, this->initialized_dual_u_product,
+                            "dual_u");
+    this->assert_array_near(this->row_maxima,
+                            this->initialized_row_maxima_product, "row_maxima");
 }
 
 
@@ -235,7 +260,7 @@ TYPED_TEST(Mc64, InitialMatching)
         this->mtx->get_const_col_idxs(), this->initialized_weights_sum,
         this->initialized_dual_u_sum, this->empty_permutation,
         this->empty_inverse_permutation, this->empty_matched_idxs,
-        this->empty_unmatched_rows, this->tol);
+        this->empty_unmatched_rows, this->zero_tol);
 
     GKO_ASSERT_ARRAY_EQ(this->empty_permutation,
                         this->initial_matching_permutation);
@@ -263,7 +288,7 @@ TYPED_TEST(Mc64, ShortestAugmentingPath)
         this->initial_matching_inverse_permutation, 4 * gko::one<index_type>(),
         this->initial_parents, this->initial_generation,
         this->initial_marked_cols, this->initial_matched_idxs, Q, q_j,
-        this->tol);
+        this->zero_tol);
 
     GKO_ASSERT_ARRAY_EQ(this->initial_matching_permutation,
                         this->final_permutation);
@@ -273,9 +298,12 @@ TYPED_TEST(Mc64, ShortestAugmentingPath)
     GKO_ASSERT_ARRAY_EQ(this->initial_generation, this->final_generation);
     GKO_ASSERT_ARRAY_EQ(this->initial_marked_cols, this->final_marked_cols);
     GKO_ASSERT_ARRAY_EQ(this->initial_matched_idxs, this->final_matched_idxs);
-    GKO_ASSERT_ARRAY_EQ(this->initialized_weights_sum, this->final_weights);
-    GKO_ASSERT_ARRAY_EQ(this->initialized_dual_u_sum, this->final_dual_u);
-    GKO_ASSERT_ARRAY_EQ(this->initialized_distance, this->final_distance);
+    this->assert_array_near(this->initialized_weights_sum, this->final_weights,
+                            "weights");
+    this->assert_array_near(this->initialized_dual_u_sum, this->final_dual_u,
+                            "dual_u");
+    this->assert_array_near(this->initialized_distance, this->final_distance,
+                            "distance");
 }
 
 
@@ -326,18 +354,18 @@ TYPED_TEST(Mc64, CreatesCorrectPermutationAndScalingExampleProduct)
     ASSERT_EQ(perm[3], 4);
     ASSERT_EQ(perm[4], 0);
     ASSERT_EQ(perm[5], 2);
-    GKO_ASSERT_NEAR(row_scaling[0], value_type{1. / 3.}, this->tolerance);
-    GKO_ASSERT_NEAR(row_scaling[1], value_type{0.2}, this->tolerance);
-    GKO_ASSERT_NEAR(row_scaling[2], value_type{0.2}, this->tolerance);
-    GKO_ASSERT_NEAR(row_scaling[3], value_type{4. / 15.}, this->tolerance);
-    GKO_ASSERT_NEAR(row_scaling[4], value_type{0.3}, this->tolerance);
-    GKO_ASSERT_NEAR(row_scaling[5], value_type{2. / 15.}, this->tolerance);
-    GKO_ASSERT_NEAR(col_scaling[0], value_type{1.}, this->tolerance);
-    GKO_ASSERT_NEAR(col_scaling[1], value_type{1.5}, this->tolerance);
-    GKO_ASSERT_NEAR(col_scaling[2], value_type{0.9375}, this->tolerance);
-    GKO_ASSERT_NEAR(col_scaling[3], value_type{5. / 6.}, this->tolerance);
-    GKO_ASSERT_NEAR(col_scaling[4], value_type{1.}, this->tolerance);
-    GKO_ASSERT_NEAR(col_scaling[5], value_type{1.25}, this->tolerance);
+    GKO_ASSERT_NEAR(row_scaling[0], value_type{1. / 3.}, r<value_type>::value);
+    GKO_ASSERT_NEAR(row_scaling[1], value_type{0.2}, r<value_type>::value);
+    GKO_ASSERT_NEAR(row_scaling[2], value_type{0.2}, r<value_type>::value);
+    GKO_ASSERT_NEAR(row_scaling[3], value_type{4. / 15.}, r<value_type>::value);
+    GKO_ASSERT_NEAR(row_scaling[4], value_type{0.3}, r<value_type>::value);
+    GKO_ASSERT_NEAR(row_scaling[5], value_type{2. / 15.}, r<value_type>::value);
+    GKO_ASSERT_NEAR(col_scaling[0], value_type{1.}, r<value_type>::value);
+    GKO_ASSERT_NEAR(col_scaling[1], value_type{1.5}, r<value_type>::value);
+    GKO_ASSERT_NEAR(col_scaling[2], value_type{0.9375}, r<value_type>::value);
+    GKO_ASSERT_NEAR(col_scaling[3], value_type{5. / 6.}, r<value_type>::value);
+    GKO_ASSERT_NEAR(col_scaling[4], value_type{1.}, r<value_type>::value);
+    GKO_ASSERT_NEAR(col_scaling[5], value_type{1.25}, r<value_type>::value);
 }
 
 
@@ -366,7 +394,7 @@ TYPED_TEST(Mc64, CreatesCorrectPermutationAndScalingLargeTrivialExampleProduct)
 
     mtx = mtx->scale_permute(row_perm, col_perm);
 
-    GKO_ASSERT_MTX_NEAR(mtx, expected_result, this->tolerance);
+    GKO_ASSERT_MTX_NEAR(mtx, expected_result, r<value_type>::value);
 }