Add possibility to split sparse matrix into block matrix

block_structures/block_matrix_t.cpp

@@ -198,4 +198,36 @@ namespace vtint {
     block_matrix_t operator*(const double &s, const block_matrix_t &v) {
         return v*s;
     }
+
+    block_matrix_t split_to_block_matrix(const Eigen::SparseMatrix<double> &A, Eigen::Index new_rows, Eigen::Index new_cols) {
+        Eigen::Index const rows = A.rows();
+        Eigen::Index const cols = A.cols();
+        if(rows % new_rows != 0 || cols % new_cols != 0){
+            throw std::runtime_error("Incompatible dimensions");
+        }
+        Eigen::Index rows_per_block = rows / new_rows;
+        Eigen::Index cols_per_block = cols / new_cols;
+        block_matrix_t A_block(new_rows, new_cols);
+
+        auto const A_triplets = sparseMatrixToTriplets(A);
+        //create a vector, with the triplets for each block
+        std::vector<std::vector<std::vector<Eigen::Triplet<double>>>> block_triplets(new_rows, std::vector<std::vector<Eigen::Triplet<double>>>(new_cols));
+        //then we go over all triplets and fill the block accordingly
+        for(auto const& t : A_triplets){
+            Eigen::Index const block_row = t.row() / rows_per_block;
+            Eigen::Index const block_col = t.col() / cols_per_block;
+            block_triplets[block_row][block_col].push_back(Eigen::Triplet<double>(t.row() % rows_per_block, t.col() % cols_per_block, t.value()));
+        }
+
+        //compress all blocks
+        for(Eigen::Index i=0; i<new_rows; ++i){
+            for(Eigen::Index j=0; j<new_cols; ++j){
+                A_block.get_block(i, j) = Eigen::SparseMatrix<double>(rows / new_rows, cols / new_cols);
+                A_block.get_block(i, j).setFromTriplets(block_triplets[i][j].begin(), block_triplets[i][j].end());
+            }
+        }
+
+        return A_block;
+
+    }
 } // vtint
\ No newline at end of file

block_structures/block_matrix_t.hpp

@@ -50,6 +50,8 @@ namespace vtint {
     block_matrix_t operator*(block_matrix_t const& v, double const& s);
     block_matrix_t operator*(double const& s, block_matrix_t const& v);
 
+    block_matrix_t split_to_block_matrix(Eigen::SparseMatrix<double> const& A, Eigen::Index new_rows, Eigen::Index new_cols);
+
 } // vtint
 
 #endif //VTINT_BLOCK_MATRIX_T_HPP

block_structures/tests/block_matrix_test.cpp

@@ -70,6 +70,13 @@ TEST_CASE("Block matrix", "[block_structures][matrix]") {
                         4, 4, 4, 5, 5, 5, 5, 6, 6,
                         4, 4, 4, 5, 5, 5, 5, 6, 6;
         REQUIRE(V_dense == V_control);
+        SECTION("Split into Blockmatrix"){
+            auto const v2 = vtint::split_to_block_matrix(V, 2, 3);
+            REQUIRE(v2.get_n_rows() == 2);
+            REQUIRE(v2.get_n_cols() == 3);
+            auto const V2_dense = v2.to_sparse_matrix().toDense();
+            REQUIRE(V2_dense == V_control);
+        }
     }SECTION("Arithmetic"){
         Eigen::SparseMatrix<double> const V = v.to_sparse_matrix();
         SECTION("Addition"){