Add more tests

core/components/addressable_pq.hpp

@@ -34,9 +34,13 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #define GKO_CORE_COMPONENTS_ADRESSABLE_PQ_HPP_
 
 
+#include <algorithm>
 #include <vector>
 
 
+#include <ginkgo/core/base/types.hpp>
+
+
 namespace gko {
 
 
@@ -46,11 +50,13 @@ namespace gko {
  * It allows inserting key-value pairs, modifying their key as well as accessing
  * and removing the key-value pair with the minimum key.
  *
- * @tparam Degree_Log2 the binary logarithm of the heap arity, i.e.,
- *         `k = 1 << Degree_Log2`
+ * @tparam KeyType    The type of the keys
+ * @tparam ValueType  The type of the values
  */
-template <typename KeyType, typename ValueType, int Degree_Log2>
+template <typename KeyType, typename ValueType>
 struct addressable_priority_queue {
+    explicit addressable_priority_queue(int deg_log2) : degree{1 << deg_log2} {}
+
     /**
      * Inserts the given key-value pair into the PQ.
      * Duplicate keys are allowed, they may be returned in an arbitrary order.
@@ -70,12 +76,14 @@ struct addressable_priority_queue {
         return handle;
     }
 
-    /** Updates the key of the pair with the given handle. */
+    /**
+     * Updates the key of the pair with the given handle.
+     */
     void update_key(std::size_t handle, KeyType new_key)
     {
         auto pos = m_handle_pos[handle];
-        assert(pos < size());
-        assert(m_handles[pos] == handle);
+        GKO_ASSERT(pos < size());
+        GKO_ASSERT(m_handles[pos] == handle);
         auto old_key = m_keys[pos];
         m_keys[pos] = new_key;
         if (old_key < new_key) {
@@ -85,16 +93,24 @@ struct addressable_priority_queue {
         }
     }
 
-    /** Returns the minimum key from the queue. */
+    /**
+     * Returns the minimum key from the queue.
+     */
     KeyType min_key() const { return m_keys[0]; }
 
-    /** Returns the value belonging to the minimum key from the queue. */
+    /**
+     * Returns the value belonging to the minimum key from the queue.
+     */
     ValueType min_val() const { return m_values[0]; }
 
-    /** Returns the key-value pair with the minimum key from the queue. */
+    /**
+     * Returns the key-value pair with the minimum key from the queue.
+     */
     std::pair<KeyType, ValueType> min() const { return {min_key(), min_val()}; }
 
-    /** Removes the key-value pair with the minimum key from the queue. */
+    /**
+     * Removes the key-value pair with the minimum key from the queue.
+     */
     void pop_min()
     {
         swap(0, size() - 1);
@@ -105,10 +121,14 @@ struct addressable_priority_queue {
         sift_down(0);
     }
 
-    /** Returns the number of key-value pairs in the queue. */
+    /**
+     * Returns the number of key-value pairs in the queue.
+     */
     std::size_t size() const { return m_keys.size(); }
 
-    /** Returns true if and only if the queue has size 0. */
+    /**
+     * Returns true if and only if the queue has size 0.
+     */
     bool empty() const { return size() == 0; }
 
     void reset()
@@ -120,8 +140,7 @@ struct addressable_priority_queue {
     }
 
 private:
-    constexpr static int degree = 1 << 4;       // Degree_Log2;
-    constexpr static auto invalid_handle = -1;  //((std::size_t)-1);
+    // constexpr static int degree = 1 << Degree_Log2;
 
     std::size_t parent(std::size_t i) const { return (i - 1) / degree; }
 
@@ -166,6 +185,7 @@ struct addressable_priority_queue {
 
     std::size_t next_handle() const { return m_handle_pos.size(); }
 
+    const int degree;
     std::vector<KeyType> m_keys;
     std::vector<ValueType> m_values;
     std::vector<std::size_t> m_handles;

core/reorder/mc64.cpp

@@ -77,47 +77,66 @@ void Mc64<ValueType, IndexType>::generate(std::shared_ptr<const Executor>& exec,
     size_type num_rows = mtx->get_size()[0];
     size_type nnz = mtx->get_num_stored_elements();
 
-    array<remove_complex<ValueType>> workspace{exec, nnz + 3 * num_rows};
+    // A real valued array with space for:
+    //     - nnz entries for weights
+    //     - num_rows entries each for the dual vector u, distance information
+    //       and the max weight per row
+    array<remove_complex<ValueType>> value_workspace{exec, nnz + 3 * num_rows};
+    // A zero initialized index array with space for n entries each for parent
+    // information, priority queue handles, generation information, marked
+    // columns, indices corresponding to matched columns in the according row
+    // and still unmatched rows
+    array<IndexType> index_workspace{exec, 6 * num_rows};
+    index_workspace.fill(0);
+
     array<IndexType> permutation{exec, num_rows};
     array<IndexType> inv_permutation{exec, num_rows};
     permutation.fill(-one<IndexType>());
     inv_permutation.fill(-one<IndexType>());
+
     const auto row_ptrs = mtx->get_const_row_ptrs();
     const auto col_idxs = mtx->get_const_col_idxs();
 
-    exec->run(mc64::make_initialize_weights(mtx.get(), workspace,
+    exec->run(mc64::make_initialize_weights(mtx.get(), value_workspace,
                                             parameters_.strategy));
 
-    array<IndexType> parents{exec, 6 * num_rows};
-    parents.fill(0);
-    exec->run(mc64::make_initial_matching(num_rows, row_ptrs, col_idxs,
-                                          workspace, permutation,
-                                          inv_permutation, parents));
-
-    addressable_priority_queue<remove_complex<ValueType>, IndexType, 2> Q{};
+    // Compute an initial extreme 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(mc64::make_initial_matching(
+        num_rows, row_ptrs, col_idxs, value_workspace, permutation,
+        inv_permutation, index_workspace, parameters_.tolerance));
+
+    // For each row that is not contained in the initial matching, search for
+    // an augmenting path, update the matching and compute the new entries
+    // of the dual vectors.
+    addressable_priority_queue<remove_complex<ValueType>, IndexType> Q{
+        parameters_.log2_degree};
     std::vector<IndexType> q_j{};
-    const auto unmatched = parents.get_data() + 5 * num_rows;
+    const auto unmatched = index_workspace.get_data() + 5 * num_rows;
     auto um = 0;
     auto root = unmatched[um];
     while (root != 0 && um < num_rows) {
         if (root != -1)
             exec->run(mc64::make_shortest_augmenting_path(
-                num_rows, row_ptrs, col_idxs, workspace, permutation,
-                inv_permutation, root, parents, Q, q_j));
+                num_rows, row_ptrs, col_idxs, value_workspace, permutation,
+                inv_permutation, root, index_workspace, Q, q_j,
+                parameters_.tolerance));
         root = unmatched[++um];
     }
 
-    permutation_ = std::move(share(PermutationMatrix::create(
-        exec, system_matrix->get_size(), permutation,
-        gko::matrix::row_permute | matrix::inverse_permute)));
-    inv_permutation_ = std::move(share(
+    permutation_ = std::move(share(
         PermutationMatrix::create(exec, system_matrix->get_size(),
-                                  inv_permutation, matrix::column_permute)));
+                                  inv_permutation, gko::matrix::row_permute)));
+    inv_permutation_ = std::move(share(PermutationMatrix::create(
+        exec, system_matrix->get_size(), permutation, matrix::row_permute)));
     row_scaling_ = std::move(DiagonalMatrix::create(exec, num_rows));
     col_scaling_ = std::move(DiagonalMatrix::create(exec, num_rows));
+
     exec->run(mc64::make_compute_scaling(
-        mtx.get(), workspace, permutation, parents, parameters_.strategy,
-        row_scaling_.get(), col_scaling_.get()));
+        mtx.get(), value_workspace, permutation, index_workspace,
+        parameters_.strategy, row_scaling_.get(), col_scaling_.get()));
 }
 
 

core/reorder/mc64_kernels.hpp

@@ -52,41 +52,43 @@ namespace gko {
 namespace kernels {
 
 
-#define GKO_DECLARE_MC64_INITIALIZE_WEIGHTS_KERNEL(ValueType, IndexType)  \
-    void initialize_weights(std::shared_ptr<const DefaultExecutor> exec,  \
-                            const matrix::Csr<ValueType, IndexType>* mtx, \
-                            array<remove_complex<ValueType>>& workspace,  \
+#define GKO_DECLARE_MC64_INITIALIZE_WEIGHTS_KERNEL(ValueType, IndexType)       \
+    void initialize_weights(std::shared_ptr<const DefaultExecutor> exec,       \
+                            const matrix::Csr<ValueType, IndexType>* mtx,      \
+                            array<remove_complex<ValueType>>& value_workspace, \
                             gko::reorder::reordering_strategy strategy)
 
 
 #define GKO_DECLARE_MC64_INITIAL_MATCHING_KERNEL(ValueType, IndexType)    \
     void initial_matching(                                                \
         std::shared_ptr<const DefaultExecutor> exec, size_type num_rows,  \
         const IndexType* row_ptrs, const IndexType* col_idxs,             \
-        const array<ValueType>& workspace, array<IndexType>& permutation, \
-        array<IndexType>& inv_permutation, array<IndexType>& parents)
+        const array<ValueType>& value_workspace,                          \
+        array<IndexType>& permutation, array<IndexType>& inv_permutation, \
+        array<IndexType>& index_workspace, ValueType tolerance)
 
 
 #define GKO_DECLARE_MC64_SHORTEST_AUGMENTING_PATH_KERNEL(ValueType, IndexType) \
     void shortest_augmenting_path(                                             \
         std::shared_ptr<const DefaultExecutor> exec, size_type num_rows,       \
         const IndexType* row_ptrs, const IndexType* col_idxs,                  \
-        array<ValueType>& workspace, array<IndexType>& permutation,            \
+        array<ValueType>& value_workspace, array<IndexType>& permutation,      \
         array<IndexType>& inv_permutation, IndexType root,                     \
-        array<IndexType>& parents,                                             \
-        addressable_priority_queue<ValueType, IndexType, 2>& Q,                \
-        std::vector<IndexType>& q_j)
-
-
-#define GKO_DECLARE_MC64_COMPUTE_SCALING_KERNEL(ValueType, IndexType)       \
-    void compute_scaling(std::shared_ptr<const DefaultExecutor> exec,       \
-                         const matrix::Csr<ValueType, IndexType>* mtx,      \
-                         const array<remove_complex<ValueType>>& workspace, \
-                         const array<IndexType>& permutation,               \
-                         const array<IndexType>& parents,                   \
-                         gko::reorder::reordering_strategy strategy,        \
-                         gko::matrix::Diagonal<ValueType>* row_scaling,     \
-                         gko::matrix::Diagonal<ValueType>* col_scaling)
+        array<IndexType>& index_workspace,                                     \
+        addressable_priority_queue<ValueType, IndexType>& Q,                   \
+        std::vector<IndexType>& q_j, ValueType tolerance)
+
+
+#define GKO_DECLARE_MC64_COMPUTE_SCALING_KERNEL(ValueType, IndexType) \
+    void compute_scaling(                                             \
+        std::shared_ptr<const DefaultExecutor> exec,                  \
+        const matrix::Csr<ValueType, IndexType>* mtx,                 \
+        const array<remove_complex<ValueType>>& value_workspace,      \
+        const array<IndexType>& permutation,                          \
+        const array<IndexType>& index_workspace,                      \
+        gko::reorder::reordering_strategy strategy,                   \
+        gko::matrix::Diagonal<ValueType>* row_scaling,                \
+        gko::matrix::Diagonal<ValueType>* col_scaling)
 
 
 #define GKO_DECLARE_ALL_AS_TEMPLATES                                        \

core/test/components/CMakeLists.txt

@@ -1 +1,2 @@
+ginkgo_create_test(addressable_pq)
 ginkgo_create_test(disjoint_sets)