[Code] Add sequential sub-matrix and kronecker impl

Author: EgorOrachyov

README.md

@@ -29,13 +29,16 @@ as well as python high-level wrapper with automated resources management.
 - [X] CSR transpose
 - [X] CSR submatrix
 - [X] CSR matrix reduce
-- [ ] CSR slicing
+- [X] CSR slicing
+- [ ] Sequential fallback backend for CPU
 - [ ] IO matrix loading from mtx file
 - [ ] IO matrix saving into mtx file
 - [ ] IO matrix saving into gviz format
 - [ ] IO user-defined file logging
 - [X] Wrapper for Python API
-- [ ] Wrapper tests in Python 
+- [ ] Wrapper syntax sugar
+- [ ] Tests for Python wrapper
+- [ ] Code examples
 - [ ] User guide
 - [X] Unit Tests collection
 - [X] Dummy library implementation for testing

cubool/sources/core/matrix.cpp

@@ -41,8 +41,8 @@ namespace cubool {
     }
 
     void Matrix::build(const index *rows, const index *cols, size_t nvals, bool isSorted) {
-        CHECK_RAISE_ERROR(rows != nullptr, InvalidArgument, "Null ptr rows array");
-        CHECK_RAISE_ERROR(cols != nullptr, InvalidArgument, "Null ptr cols array");
+        CHECK_RAISE_ERROR(rows != nullptr || nvals == 0, InvalidArgument, "Null ptr rows array");
+        CHECK_RAISE_ERROR(cols != nullptr || nvals == 0, InvalidArgument, "Null ptr cols array");
 
         mHnd->build(rows, cols, nvals, isSorted);
     }
@@ -69,8 +69,8 @@ namespace cubool {
         CHECK_RAISE_ERROR(bI <= other->getNrows(), InvalidArgument, "Provided sub-matrix range must be within matrix bounds");
         CHECK_RAISE_ERROR(bJ <= other->getNcols(), InvalidArgument, "Provided sub-matrix range must be within matrix bounds");
 
-        CHECK_RAISE_ERROR(nrows <= this->getNrows(), InvalidArgument, "Result matrix has incompatible size for extracted sub-matrix range");
-        CHECK_RAISE_ERROR(ncols <= this->getNcols(), InvalidArgument, "Result matrix has incompatible size for extracted sub-matrix range");
+        CHECK_RAISE_ERROR(nrows == this->getNrows(), InvalidArgument, "Result matrix has incompatible size for extracted sub-matrix range");
+        CHECK_RAISE_ERROR(ncols == this->getNcols(), InvalidArgument, "Result matrix has incompatible size for extracted sub-matrix range");
 
         mHnd->extractSubMatrix(*other->mHnd, i, j, nrows, ncols);
     }

cubool/sources/cuBool_Matrix_Build.cpp

@@ -34,8 +34,6 @@ cuBool_Status cuBool_Matrix_Build(
     CUBOOL_BEGIN_BODY
         CUBOOL_VALIDATE_LIBRARY
         CUBOOL_ARG_NOT_NULL(matrix)
-        CUBOOL_ARG_NOT_NULL(rows)
-        CUBOOL_ARG_NOT_NULL(cols)
         auto m = (cubool::Matrix *) matrix;
         m->build(rows, cols, nvals, hints & CUBOOL_HINT_VALUES_SORTED);
     CUBOOL_END_BODY

cubool/sources/sequential/sq_kronecker.cpp

@@ -26,4 +26,55 @@
 
 namespace cubool {
 
+    void scan(std::vector<index> &v) {
+        index sum = 0;
+        for (auto& val: v) {
+            index next = sum + val;
+            val = sum;
+            sum = next;
+        }
+    }
+
+    void sq_kronecker(const CooData& a, const CooData& b, CooData& out) {
+        size_t nvals = a.mNvals * b.mNvals;
+
+        out.mNvals = nvals;
+        out.mRowIndices.resize(nvals);
+        out.mColIndices.resize(nvals);
+
+        std::vector<index> aoffsets(a.mNrows + 1, 0);
+        std::vector<index> boffsets(b.mNrows + 1, 0);
+
+        for (index k = 0; k < a.mNvals; k++) {
+            aoffsets[a.mRowIndices[k]]++;
+        }
+
+        for (index k = 0; k < b.mNvals; k++) {
+            boffsets[b.mRowIndices[k]]++;
+        }
+
+        scan(aoffsets);
+        scan(boffsets);
+
+        size_t id = 0;
+
+        for (index ai = 0; ai < a.mNrows; ai++) {
+            for (index bi = 0; bi < b.mNrows; bi++) {
+                index rowId = ai * b.mNrows + bi;
+
+                for (index k = aoffsets[ai]; k < aoffsets[ai + 1]; k++) {
+                    index colIdBase = a.mColIndices[k] * b.mNcols;
+
+                    for (index l = boffsets[bi]; l < boffsets[bi + 1]; l++) {
+                        index colId = colIdBase + b.mColIndices[l];
+
+                        out.mRowIndices[id] = rowId;
+                        out.mColIndices[id] = colId;
+                        id += 1;
+                    }
+                }
+            }
+        }
+    }
+
 }

cubool/sources/sequential/sq_kronecker.hpp

@@ -25,8 +25,18 @@
 #ifndef CUBOOL_SP_KRONECKER_HPP
 #define CUBOOL_SP_KRONECKER_HPP
 
+#include <sequential/sq_coo_data.hpp>
+
 namespace cubool {
 
+    /**
+     * Kronecker product of `a` and `b` matrices.
+     * @param a Input matrix
+     * @param b Input matrix
+     * @param[out] out Result matrix
+     */
+    void sq_kronecker(const CooData& a, const CooData& b, CooData& out);
+
 }
 
 #endif //CUBOOL_SP_KRONECKER_HPP

cubool/sources/sequential/sq_matrix.cpp

@@ -24,6 +24,8 @@
 
 #include <sequential/sq_matrix.hpp>
 #include <sequential/sq_transpose.hpp>
+#include <sequential/sq_submatrix.hpp>
+#include <sequential/sq_kronecker.hpp>
 #include <core/error.hpp>
 #include <algorithm>
 #include <cassert>
@@ -109,7 +111,17 @@ namespace cubool {
     }
 
     void SqMatrix::extractSubMatrix(const MatrixBase &otherBase, index i, index j, index nrows, index ncols) {
+        auto other = dynamic_cast<const SqMatrix*>(&otherBase);
+
+        CHECK_RAISE_ERROR(other != nullptr, InvalidArgument, "Provided matrix does not belongs to sequential matrix class");
+
+        auto M = this->getNrows();
+        auto N = this->getNcols();
 
+        assert(M == nrows);
+        assert(N == ncols);
+
+        sq_submatrix(other->mData, this->mData, i, j, nrows, ncols);
     }
 
     void SqMatrix::clone(const MatrixBase &otherBase) {
@@ -154,7 +166,21 @@ namespace cubool {
     }
 
     void SqMatrix::kronecker(const MatrixBase &aBase, const MatrixBase &bBase) {
+        auto a = dynamic_cast<const SqMatrix*>(&aBase);
+        auto b = dynamic_cast<const SqMatrix*>(&bBase);
+
+        CHECK_RAISE_ERROR(a != nullptr, InvalidArgument, "Provided matrix does not belongs to sequential matrix class");
+        CHECK_RAISE_ERROR(b != nullptr, InvalidArgument, "Provided matrix does not belongs to sequential matrix class");
+
+        auto M = a->getNrows();
+        auto N = a->getNcols();
+        auto K = b->getNrows();
+        auto T = b->getNcols();
+
+        assert(M * K == this->getNrows());
+        assert(N * T == this->getNcols());
 
+        sq_kronecker(a->mData, b->mData, this->mData);
     }
 
     void SqMatrix::eWiseAdd(const MatrixBase &aBase, const MatrixBase &bBase) {

cubool/sources/sequential/sq_submatrix.cpp

@@ -26,4 +26,37 @@
 
 namespace cubool {
 
+    void sq_submatrix(const CooData& a, CooData& sub, index i, index j, index nrows, index ncols) {
+        index first = a.mNvals;
+        index last = 0;
+        size_t nvals = 0;
+
+        for (index k = 0; k < a.mNvals; k++) {
+            auto ai = a.mRowIndices[k];
+            auto aj = a.mColIndices[k];
+
+            if (i <= ai && ai < i + nrows && j <= aj && aj < j + ncols) {
+                first = std::min(first, k);
+                last = std::max(last, k);
+                nvals += 1;
+            }
+        }
+
+        sub.mNvals = nvals;
+        sub.mRowIndices.resize(nvals);
+        sub.mColIndices.resize(nvals);
+
+        size_t idx = 0;
+        for (index k = first; k <= last; k++) {
+            auto ai = a.mRowIndices[k];
+            auto aj = a.mColIndices[k];
+
+            if (i <= ai && ai < i + nrows && j <= aj && aj < j + ncols) {
+                sub.mRowIndices[idx] = ai - i;
+                sub.mColIndices[idx] = aj - j;
+                idx += 1;
+            }
+        }
+    }
+
 }

cubool/sources/sequential/sq_submatrix.hpp

@@ -25,8 +25,22 @@
 #ifndef CUBOOL_SQ_SUBMATRIX_HPP
 #define CUBOOL_SQ_SUBMATRIX_HPP
 
+#include <sequential/sq_coo_data.hpp>
+
 namespace cubool {
 
+    /**
+     * Extracts sub-matrix from matrix `a`
+     *
+     * @param a Source
+     * @param[out] sub Result
+     * @param i First sub-matrix row
+     * @param j First sub-matrix col
+     * @param nrows Sub-matrix size
+     * @param ncols Sub-matrix size
+     */
+    void sq_submatrix(const CooData& a, CooData& sub, index i, index j, index nrows, index ncols);
+
 }
 
 #endif //CUBOOL_SQ_SUBMATRIX_HPP

cubool/sources/sequential/sq_transpose.cpp

@@ -27,7 +27,7 @@
 namespace cubool {
 
     void sq_transpose(const CooData& a, CooData& at) {
-        std::vector<index> offsets(a.mNcols + 1, 0);
+        std::vector<index> offsets(a.mNcols, 0);
 
         for (size_t k = 0; k < a.mNvals; k++) {
             offsets[a.mColIndices[k]]++;

cubool/tests/test_matrix_extract_sub_matrix.cpp

@@ -30,16 +30,23 @@ void testMatrixExtractSubMatrix(cuBool_Index m, cuBool_Index n, cuBool_Index N,
     auto ta = testing::Matrix::generateSparse(m, n, density);
 
     ASSERT_EQ(cuBool_Matrix_New(&a, m, n), CUBOOL_STATUS_SUCCESS);
-
     ASSERT_EQ(cuBool_Matrix_Build(a, ta.mRowsIndex.data(), ta.mColsIndex.data(), ta.mNvals, CUBOOL_HINT_VALUES_SORTED), CUBOOL_STATUS_SUCCESS);
 
+    testing::TimeQuery query;
+
     for (size_t i = 0; i < N; i++) {
         for (size_t j = 0; j < N; j++) {
             auto k = m / N;
             auto l = n / N;
 
+            testing::Timer t;
+
             ASSERT_EQ(cuBool_Matrix_New(&r, k, l), CUBOOL_STATUS_SUCCESS);
-            ASSERT_EQ(cuBool_Matrix_ExtractSubMatrix(r, a, i * k, j * l, k, l, CUBOOL_HINT_NO), CUBOOL_STATUS_SUCCESS);
+
+            {
+                testing::TimeScope scope(query);
+                ASSERT_EQ(cuBool_Matrix_ExtractSubMatrix(r, a, i * k, j * l, k, l, CUBOOL_HINT_NO), CUBOOL_STATUS_SUCCESS);
+            }
 
             auto tr = ta.subMatrix(i * k, j * l, k, l);
 
@@ -48,6 +55,7 @@ void testMatrixExtractSubMatrix(cuBool_Index m, cuBool_Index n, cuBool_Index N,
         }
     }
 
+    std::cout << query.getAverageTimeMs() << "ms" << std::endl;
 
     ASSERT_EQ(cuBool_Matrix_Free(a), CUBOOL_STATUS_SUCCESS);
 }
@@ -84,4 +92,22 @@ TEST(cuBool_Matrix, SubMatrixExtractLarge) {
     testRun(m, n, step, CUBOOL_HINT_NO);
 }
 
+TEST(cuBool_Matrix, SubMatrixExtractSmallFallback) {
+    cuBool_Index m = 100, n = 200;
+    float step = 0.05f;
+    testRun(m, n, step, CUBOOL_HINT_CPU_BACKEND);
+}
+
+TEST(cuBool_Matrix, SubMatrixExtractMediumFallback) {
+    cuBool_Index m = 400, n = 700;
+    float step = 0.05f;
+    testRun(m, n, step, CUBOOL_HINT_CPU_BACKEND);
+}
+
+TEST(cuBool_Matrix, SubMatrixExtractLargeFallback) {
+    cuBool_Index m = 2000, n = 4000;
+    float step = 0.01f;
+    testRun(m, n, step, CUBOOL_HINT_CPU_BACKEND);
+}
+
 CUBOOL_GTEST_MAIN