Add reordering function.

ypodlesov · May 19, 2024 · 2e1f4e9 · 2e1f4e9
1 parent f8ddea6
commit 2e1f4e9
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 31 deletions.
diff --git a/basic_la/vector.h b/basic_la/vector.h
@@ -40,7 +40,7 @@ struct Vector: public CommonContainer<Vector<T>, T, Hold> {
     ~Vector() = default;
 };
 
-template <typename T, typename std::enable_if<std::is_floating_point<T>::value>::type* = nullptr>
+template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type* = nullptr>
 bool operator ==(const Vector<T>& a, const Vector<T>& b) {
     assert(a.mem_size_ == b.mem_size_);
     T diff = 0;

diff --git a/matrix_powers_mv/matrix_powers_mv.cpp b/matrix_powers_mv/matrix_powers_mv.cpp
@@ -1,15 +1,15 @@
 #include "matrix_powers_mv.h"
+#include <cstdint>
+#include <cstring>
 
-int ReorderMatrix(const SparseMatrix<double>& sp_matrix) {
+bool ReorderMatrix(SparseMatrix<double>& sp_matrix) {
     SparseMatrix<double> a_no_diag;
     sp_matrix.RemoveDiag(a_no_diag);
     idx_t* perm = new idx_t[a_no_diag.row_cnt_];
     idx_t* iperm = new idx_t[a_no_diag.row_cnt_];
     idx_t options[METIS_NOPTIONS];
     METIS_SetDefaultOptions(options);
     options[METIS_OPTION_NUMBERING] = 0;
-    // int METIS_NodeND(idx_t *nvtxs, idx_t *xadj, idx_t *adjncy, idx_t *vwgt,
-    //   idx_t *options, idx_t *perm, idx_t *iperm);
 
     int result = METIS_NodeND(
         &a_no_diag.row_cnt_
@@ -20,9 +20,22 @@ int ReorderMatrix(const SparseMatrix<double>& sp_matrix) {
         , perm
         , iperm
     );
-    return result;
+
+    if (result != METIS_OK) {
+        return false;
+    }
+    for (int64_t i = 0; i < sp_matrix.row_cnt_; ++i) {
+        int64_t cur_row_cnt = sp_matrix.i_a_[i + 1] - sp_matrix.i_a_[i];
+        std::memcpy(&sp_matrix.data_[sp_matrix.i_a_[i]], &sp_matrix.data_[sp_matrix.i_a_[perm[i]]], cur_row_cnt * sizeof(double));
+        std::memcpy(&sp_matrix.j_a_[sp_matrix.i_a_[i]], &sp_matrix.j_a_[sp_matrix.i_a_[perm[i]]], cur_row_cnt * sizeof(int64_t));
+        sp_matrix.i_a_[i] = sp_matrix.i_a_[perm[i]];
+    }
+    return true;
 }
 
-void MatrixPowersMV(const SparseMatrix<double>& /* sp_matrix */, const Vector<double>& /* x */, const Vector<Vector<double>*>& /* res */) {
-    return;
+bool MatrixPowersMV(SparseMatrix<double>& sp_matrix, const Vector<double>& /* x */, const std::vector<Vector<double>*>& /* res */) {
+    if (!ReorderMatrix(sp_matrix)) {
+        return false;
+    }
+    return true;
 }
diff --git a/matrix_powers_mv/matrix_powers_mv.h b/matrix_powers_mv/matrix_powers_mv.h
@@ -4,5 +4,5 @@
 #include <sparse_matrix.h>
 #include <vector.h>
 
-int ReorderMatrix(const SparseMatrix<double>& sp_matrix);
-void MatrixPowersMV(const SparseMatrix<double>& sp_matrix, const Vector<double>& x, const Vector<Vector<double>*>& res);
+bool ReorderMatrix(SparseMatrix<double>& sp_matrix);
+bool MatrixPowersMV(SparseMatrix<double>& sp_matrix, const Vector<double>& x, const std::vector<Vector<double>*>& res);
diff --git a/matrix_powers_mv/test.cpp b/matrix_powers_mv/test.cpp
@@ -21,29 +21,12 @@ static void Test(const uint32_t n) {
         fstream >> a;
         fstream.close();
         REQUIRE(a.data_);
-        SparseMatrix<double> a_no_diag;
-        a.RemoveDiag(a_no_diag);
-        idx_t* perm = new idx_t[a_no_diag.row_cnt_];
-        idx_t* iperm = new idx_t[a_no_diag.row_cnt_];
-        idx_t options[METIS_NOPTIONS];
-        METIS_SetDefaultOptions(options);
-        options[METIS_OPTION_NUMBERING] = 0;
-        // int METIS_NodeND(idx_t *nvtxs, idx_t *xadj, idx_t *adjncy, idx_t *vwgt,
-        //   idx_t *options, idx_t *perm, idx_t *iperm);
-
-        int result = METIS_NodeND(
-            &a_no_diag.row_cnt_
-            , a_no_diag.i_a_
-            , a_no_diag.j_a_
-            , nullptr
-            , options
-            , perm
-            , iperm
-        );
-        REQUIRE(result == METIS_OK);
     }
-
-    // REQUIRE(ReorderMatrix(a_no_diag) == METIS_OK);
+    std::vector<Vector<double>*> res;
+    Vector<double> x(n);
+    NHelpers::GenRandomVector(x, n, true);
+
+    REQUIRE(MatrixPowersMV(a, x, res));
 }
 
 TEST_CASE("Size 128") {