[Code] Add seq reduce impl && fix cuda reduce impl

Author: EgorOrachyov

README.md

@@ -30,7 +30,7 @@ as well as python high-level wrapper with automated resources management.
 - [X] CSR submatrix
 - [X] CSR matrix reduce
 - [X] CSR slicing
-- [ ] Sequential fallback backend for CPU
+- [X] Sequential fallback backend for CPU
 - [ ] IO matrix loading from mtx file
 - [ ] IO matrix saving into mtx file
 - [ ] IO matrix saving into gviz format

cubool/sources/cuda/kernels/spreduce.cuh

@@ -50,22 +50,16 @@ namespace cubool {
             MatrixType operator()(const MatrixType& a) {
                 auto nrows = a.m_rows;
                 auto ncols = a.m_cols;
-                auto nvals = a.m_vals;
 
                 ContainerType<IndexType> rowIndex(nrows + 1);
                 ContainerType<IndexType> tmpRowIndex(nrows + 1);
 
-                // Assign for each row 0 (empty)
-                thrust::fill(tmpRowIndex.begin(), tmpRowIndex.end(), (IndexType) 0);
-
-                // For each value find row and set this row in 1
-                thrust::for_each(thrust::counting_iterator<IndexType>(0), thrust::counting_iterator<IndexType>(nvals),
-                    [aRowIndex = a.m_row_index.data(),
-                     tmpRowIndex = tmpRowIndex.data(),
-                     nrows]
-                    __device__ (IndexType valuedIdx) {
-                        auto rowId = findNearestRowIdx(valuedIdx, nrows, aRowIndex);
-                        atomicOr((tmpRowIndex + rowId).get(), (IndexType) 1);
+                // For each row if not empty then row will have 1 value
+                thrust::for_each(thrust::counting_iterator<IndexType>(0), thrust::counting_iterator<IndexType>(nrows),
+                    [aRowIndex = a.m_row_index.data(), tmpRowIndex = tmpRowIndex.data()]
+                    __device__ (IndexType i) {
+                        auto nnz = aRowIndex[i + 1] - aRowIndex[i];
+                        tmpRowIndex[i] = nnz > 0? 1: 0;
                     }
                 );
 

cubool/sources/sequential/sq_matrix.cpp

@@ -28,6 +28,7 @@
 #include <sequential/sq_kronecker.hpp>
 #include <sequential/sq_ewiseadd.hpp>
 #include <sequential/sq_spgemm.hpp>
+#include <sequential/sq_reduce.hpp>
 #include <sequential/sq_exclusive_scan.hpp>
 #include <core/error.hpp>
 #include <algorithm>
@@ -173,7 +174,17 @@ namespace cubool {
         CHECK_RAISE_ERROR(other != nullptr, InvalidArgument, "Provided matrix does not belongs to sequential matrix class");
 
         auto M = other->getNrows();
-        auto N = other->getNcols();
+        auto N = 1;
+
+        assert(M == this->getNrows());
+        assert(N == this->getNcols());
+
+        CsrData out;
+        out.nrows = this->getNrows();
+        out.ncols = this->getNcols();
+
+        sq_reduce(other->mData, out);
+        this->mData = std::move(out);
     }
 
     void SqMatrix::multiply(const MatrixBase &aBase, const MatrixBase &bBase, bool accumulate) {

cubool/sources/sequential/sq_reduce.cpp

@@ -23,7 +23,26 @@
 /**********************************************************************************/
 
 #include <sequential/sq_reduce.hpp>
+#include <sequential/sq_exclusive_scan.hpp>
 
 namespace cubool {
 
+    void sq_reduce(const CsrData& a, CsrData& out) {
+        out.rowOffsets.resize(a.nrows + 1);
+
+        for (index i = 0; i < a.nrows; i++) {
+            index nnz = a.rowOffsets[i + 1] - a.rowOffsets[i];
+
+            if (nnz > 0) {
+                out.rowOffsets[i] = 1;
+            }
+        }
+
+        sq_exclusive_scan(out.rowOffsets.begin(), out.rowOffsets.end(), 0);
+
+        out.nvals = out.rowOffsets.back();
+        out.colIndices.clear();
+        out.colIndices.resize(out.nvals, 0);
+    }
+
 }

cubool/sources/sequential/sq_reduce.hpp

@@ -25,8 +25,18 @@
 #ifndef CUBOOL_SP_REDUCE_HPP
 #define CUBOOL_SP_REDUCE_HPP
 
+#include <sequential/sq_csr_data.hpp>
+
 namespace cubool {
 
+    /**
+     * Reduce matrix `a` to column vector `out`
+     *
+     * @param a Input matrix
+     * @param[out] out Where to store result
+     */
+    void sq_reduce(const CsrData& a, CsrData& out);
+
 }
 
 #endif //CUBOOL_SP_REDUCE_HPP

cubool/tests/test_matrix_reduce.cpp

@@ -74,4 +74,22 @@ TEST(cuBool_Matrix, ReduceLarge) {
     testRun(m, n, step, CUBOOL_HINT_NO);
 }
 
+TEST(cuBool_Matrix, ReduceSmallFallback) {
+    cuBool_Index m = 100, n = 200;
+    float step = 0.05f;
+    testRun(m, n, step, CUBOOL_HINT_CPU_BACKEND);
+}
+
+TEST(cuBool_Matrix, ReduceMediumFallback) {
+    cuBool_Index m = 400, n = 700;
+    float step = 0.05f;
+    testRun(m, n, step, CUBOOL_HINT_CPU_BACKEND);
+}
+
+TEST(cuBool_Matrix, ReduceLargeFallback) {
+    cuBool_Index m = 2000, n = 4000;
+    float step = 0.01f;
+    testRun(m, n, step, CUBOOL_HINT_CPU_BACKEND);
+}
+
 CUBOOL_GTEST_MAIN