Repository: SparseLinearAlgebra/spla (language: python)

Issue:
Algorithms: BFS using library API



Discussion / hints:

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index cf987b9fa..15e8248d0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -208,7 +208,9 @@ include(include/CMakeLists.txt)
 include(sources/CMakeLists.txt)
 
 add_library(spla SHARED
+        ${SPLA_ALGO_HEADERS}
         ${SPLA_CPP_HEADERS}
+        ${SPLA_ALGO_SOURCES}
         ${SPLA_CPP_SOURCES}
         ${SPLA_ALGORITHM_SOURCES}
         ${SPLA_COMPUTE_SOURCES}
diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt
index ed5bfc94c..71c884f96 100644
--- a/include/CMakeLists.txt
+++ b/include/CMakeLists.txt
@@ -16,4 +16,9 @@ set(SPLA_CPP_HEADERS
         include/spla-cpp/SplaScalar.hpp
         include/spla-cpp/SplaType.hpp
         include/spla-cpp/SplaTypes.hpp
-        include/spla-cpp/SplaVector.hpp)
\ No newline at end of file
+        include/spla-cpp/SplaVector.hpp)
+
+set(SPLA_ALGO_HEADERS
+        include/spla-algo/SplaAlgo.hpp
+        include/spla-algo/SplaAlgoBfs.hpp
+        include/spla-algo/SplaAlgoCommon.hpp)
\ No newline at end of file
diff --git a/include/spla-algo/SplaAlgo.hpp b/include/spla-algo/SplaAlgo.hpp
new file mode 100644
index 000000000..66aeadf97
--- /dev/null
+++ b/include/spla-algo/SplaAlgo.hpp
@@ -0,0 +1,54 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLAALGO_HPP
+#define SPLA_SPLAALGO_HPP
+
+/**
+ * @defgroup Algorithm
+ *
+ * @brief Public library algorithms interface
+ *
+ * @details Algorithm module provides interface to comment graph
+ * algorithms, such as breadth-first search (bfs), single source
+ * shortest paths (ssps), triangles counting (tc), page rank,
+ * connected components and etc. implemented both using
+ * spla library API (matrix, vectors, expressions) for multi-GPU evaluation
+ * and using standard C++ primitives (for reference and conformance checks only).
+ *
+ * This algorithms can be safely used in user applications.
+ *
+ * Implementation details are hidden in private Internal sources module.
+ * Header files has no other dependencies, except standard c++ library files and core SPLA API interface.
+ *
+ * File SplaAlgo.hpp provides access to all algorithm module components.
+ */
+
+#include <spla-algo/SplaAlgoBfs.hpp>
+#include <spla-algo/SplaAlgoCommon.hpp>
+
+#endif//SPLA_SPLAALGO_HPP
diff --git a/include/spla-algo/SplaAlgoBfs.hpp b/include/spla-algo/SplaAlgoBfs.hpp
new file mode 100644
index 000000000..c98d3b4fb
--- /dev/null
+++ b/include/spla-algo/SplaAlgoBfs.hpp
@@ -0,0 +1,68 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLAALGOBFS_HPP
+#define SPLA_SPLAALGOBFS_HPP
+
+#include <spla-algo/SplaAlgoCommon.hpp>
+#include <spla-cpp/SplaMatrix.hpp>
+#include <spla-cpp/SplaVector.hpp>
+
+namespace spla {
+
+    /**
+     * @addtogroup Algorithm
+     * @{
+     */
+
+    /**
+     * @brief Breadth-first search algorithm
+     *
+     * @param[out] v Vector where to store levels of the reached vertices
+     * @param A Input adjacency matrix of the graph; must be n x n
+     * @param s Index of the source vertex to begin bfs
+     */
+    SPLA_API void Bfs(RefPtr<Vector> &v, const RefPtr<Matrix> &A, Index s);
+
+    /**
+     * @brief Breadth-first search algorithm
+     *
+     * @note Naive cpu reference algo implementation for correctness only
+     *
+     * @param[out] v Vector where to store levels of the reached vertices
+     * @param A Input adjacency matrix of the graph; must be n x n
+     * @param s Index of the source vertex to begin bfs
+     */
+    SPLA_API void Bfs(RefPtr<HostVector> &v, const RefPtr<HostMatrix> &A, Index s);
+
+    /**
+     * @}
+     */
+
+}// namespace spla
+
+#endif//SPLA_SPLAALGOBFS_HPP
diff --git a/include/spla-algo/SplaAlgoCommon.hpp b/include/spla-algo/SplaAlgoCommon.hpp
new file mode 100644
index 000000000..307207278
--- /dev/null
+++ b/include/spla-algo/SplaAlgoCommon.hpp
@@ -0,0 +1,115 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLAALGOCOMMON_HPP
+#define SPLA_SPLAALGOCOMMON_HPP
+
+#include <spla-cpp/SplaConfig.hpp>
+#include <spla-cpp/SplaData.hpp>
+
+#include <vector>
+
+namespace spla {
+
+    /**
+     * @addtogroup Algorithm
+     * @{
+     */
+
+    /**
+     * @class HostVector
+     * @brief Represents host vector
+     *
+     * Host data wrapper for the vector primitive.
+     * Used in host (cpu) reference algorithms implementations.
+     */
+    class HostVector final : public RefCnt {
+    public:
+        SPLA_API HostVector(Size nrows, std::vector<Index> rows, std::vector<unsigned char> vals);
+        SPLA_API ~HostVector() override = default;
+
+        SPLA_API RefPtr<DataVector> GetData(Library &library);
+
+        [[nodiscard]] bool HasValues() const { return !mValues.empty(); }
+
+        [[nodiscard]] Size GetNrows() const { return mNrows; }
+        [[nodiscard]] Size GetNnvals() const { return mNnvals; }
+        [[nodiscard]] Size GetElementSize() const { return mElementSize; }
+
+        [[nodiscard]] const std::vector<Index> &GetRowIndices() { return mRowIndices; }
+        [[nodiscard]] const std::vector<unsigned char> &GetValues() { return mValues; }
+
+    private:
+        Size mNrows;
+        Size mNnvals;
+        Size mElementSize;
+        std::vector<Index> mRowIndices;
+        std::vector<unsigned char> mValues;
+    };
+
+    /**
+     * @class HostMatrix
+     * @brief Represents host matrix
+     *
+     * Host data wrapper for the matrix primitive.
+     * Used in host (cpu) reference algorithms implementations.
+     */
+    class HostMatrix final : public RefCnt {
+    public:
+        SPLA_API HostMatrix(Size nrows, Size ncols, std::vector<Index> rows, std::vector<Index> cols, std::vector<unsigned char> vals);
+        SPLA_API ~HostMatrix() override = default;
+
+        SPLA_API RefPtr<DataMatrix> GetData(Library &library);
+
+        [[nodiscard]] bool HasValues() const { return !mValues.empty(); }
+
+        [[nodiscard]] Size GetNrows() const { return mNrows; }
+        [[nodiscard]] Size GetNcols() const { return mNcols; }
+        [[nodiscard]] Size GetNnvals() const { return mNnvals; }
+        [[nodiscard]] Size GetElementSize() const { return mElementSize; }
+
+        [[nodiscard]] const std::vector<Index> &GetRowIndices() { return mRowIndices; }
+        [[nodiscard]] const std::vector<Index> &GetColIndices() { return mColIndices; }
+        [[nodiscard]] const std::vector<unsigned char> &GetValues() { return mValues; }
+
+    private:
+        Size mNrows;
+        Size mNcols;
+        Size mNnvals;
+        Size mElementSize;
+        std::vector<Index> mRowIndices;
+        std::vector<Index> mColIndices;
+        std::vector<unsigned char> mValues;
+    };
+
+    /**
+     * @}
+     */
+
+}// namespace spla
+
+#endif//SPLA_SPLAALGOCOMMON_HPP
diff --git a/include/spla-cpp/Spla.hpp b/include/spla-cpp/Spla.hpp
index d0ff0f1c3..f79703086 100644
--- a/include/spla-cpp/Spla.hpp
+++ b/include/spla-cpp/Spla.hpp
@@ -28,6 +28,22 @@
 #ifndef SPLA_SPLA_HPP
 #define SPLA_SPLA_HPP
 
+/**
+ * @defgroup API
+ *
+ * @brief Public library interface
+ *
+ * @details API module provides access to all common library primitives.
+ * Primitives are exposed as a number of Object reference-counted classes.
+ * Reference-counting is used to automate objects life-time checks and
+ * to pass around this objects as expressions params in a safe way.
+ *
+ * Implementation details are hidden in private Internal sources module.
+ * Header files has no other dependencies, except standard c++ library files.
+ *
+ * File Spla.hpp provides access to all API components.
+ */
+
 #include <spla-cpp/SplaConfig.hpp>
 #include <spla-cpp/SplaData.hpp>
 #include <spla-cpp/SplaDescriptor.hpp>
diff --git a/include/spla-cpp/SplaConfig.hpp b/include/spla-cpp/SplaConfig.hpp
index 9cb6bccdf..5ce4afdc4 100644
--- a/include/spla-cpp/SplaConfig.hpp
+++ b/include/spla-cpp/SplaConfig.hpp
@@ -45,21 +45,6 @@
 namespace spla {
 
     /**
-     * @defgroup API
-     *
-     * @brief Public library interface
-     *
-     * @details API module provides access to all common library primitives.
-     * Primitives are exposed as a number of Object reference-counted classes.
-     * Reference-counting is used to automate objects life-time checks and
-     * to pass around this objects as expressions params in a safe way.
-     *
-     * Implementation details are hidden in private Internal sources module.
-     * Header files has no other dependencies, except standard c++ library files.
-     *
-     * File Spla.hpp provides access to all API components.
-     *
-     *
      * @defgroup Internal
      *
      * @brief Implementation details
@@ -83,6 +68,12 @@ namespace spla {
     #error Unsupported platfrom
 #endif
 
+    /** @brief Index of matrix/vector values */
+    using Index = unsigned int;
+
+    /** @brief Size (count) type */
+    using Size = std::size_t;
+
     /**
      * @}
      */
diff --git a/include/spla-cpp/SplaData.hpp b/include/spla-cpp/SplaData.hpp
index 147242f23..640cb11bf 100644
--- a/include/spla-cpp/SplaData.hpp
+++ b/include/spla-cpp/SplaData.hpp
@@ -115,8 +115,12 @@ namespace spla {
         /** @return Number of entries in data */
         std::size_t GetNvals() const;
 
+        /** @return Emtpy data object */
         static RefPtr<DataMatrix> Make(Library &library);
 
+        /** @return Transient data object; passed pointers must stay valid util data is consumed; */
+        static RefPtr<DataMatrix> Make(unsigned int *rows, unsigned int *cols, void *values, std::size_t nvals, Library &library);
+
     private:
         explicit DataMatrix(class Library &library);
 
@@ -164,8 +168,12 @@ namespace spla {
         /** @return Number of entries in data */
         std::size_t GetNvals() const;
 
+        /** @return Emtpy data object */
         static RefPtr<DataVector> Make(Library &library);
 
+        /** @return Transient data object; passed pointers must stay valid util data is consumed; */
+        static RefPtr<DataVector> Make(unsigned int *rows, void *values, std::size_t nvals, Library &library);
+
     private:
         explicit DataVector(class Library &library);
 
@@ -193,8 +201,12 @@ namespace spla {
         /** @return Scalar values */
         void *GetValue() const;
 
+        /** @return Emtpy data object */
         static RefPtr<DataScalar> Make(Library &library);
 
+        /** @return Transient data object; passed pointers must stay valid util data is consumed; */
+        static RefPtr<DataScalar> Make(void *value, Library &library);
+
     private:
         void *mValue = nullptr;
     };
diff --git a/include/spla-cpp/SplaExpression.hpp b/include/spla-cpp/SplaExpression.hpp
index b03deac69..4999db289 100644
--- a/include/spla-cpp/SplaExpression.hpp
+++ b/include/spla-cpp/SplaExpression.hpp
@@ -103,6 +103,19 @@ namespace spla {
          */
         void Wait();
 
+        /**
+         * Submit expression for the execution and wait until completion.
+         *
+         * @note Expression asynchronously submitted for execution
+         * @note Expression cannot be modified after submission
+         * @note Blocks current thread until evaluated or aborted.
+         * @note Check expression state after completion to ensure correctness.
+         *
+         * @see Expression::Submit();
+         * @see Expression::Wait();
+         */
+        void SubmitWait();
+
         /**
          * Makes dependency between provided expression nodes.
          * Next `succ` node will be evaluated only after `pred` node evaluation is finished.
diff --git a/sources/CMakeLists.txt b/sources/CMakeLists.txt
index ad4d75630..a7282e0ac 100644
--- a/sources/CMakeLists.txt
+++ b/sources/CMakeLists.txt
@@ -15,6 +15,10 @@ set(SPLA_CPP_SOURCES
         sources/SplaTypes.cpp
         sources/SplaVector.cpp)
 
+set(SPLA_ALGO_SOURCES
+        sources/SplaAlgoBfs.cpp
+        sources/SplaAlgoCommon.cpp)
+
 set(SPLA_ALGORITHM_SOURCES
         sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp
         sources/algo/matrix/SplaMatrixEWiseAddCOO.hpp
diff --git a/sources/SplaAlgoBfs.cpp b/sources/SplaAlgoBfs.cpp
new file mode 100644
index 000000000..4bb3a5004
--- /dev/null
+++ b/sources/SplaAlgoBfs.cpp
@@ -0,0 +1,147 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#include <spla-algo/SplaAlgoBfs.hpp>
+#include <spla-cpp/Spla.hpp>
+
+#include <compute/SplaIndicesToRowOffsets.hpp>
+#include <core/SplaError.hpp>
+
+#include <limits>
+#include <queue>
+#include <vector>
+
+void spla::Bfs(RefPtr<Vector> &sp_v, const RefPtr<Matrix> &sp_A, Index s) {
+    CHECK_RAISE_ERROR(sp_A.IsNotNull(), NullPointer, "Passed null argument");
+    CHECK_RAISE_ERROR(sp_A->GetNrows() == sp_A->GetNcols(), DimensionMismatch, "Matrix must be nxn");
+    CHECK_RAISE_ERROR(s < sp_A->GetNrows(), InvalidArgument, "Start index must be withing A bounds");
+
+    auto &library = sp_A->GetLibrary();
+    auto n = sp_A->GetNrows();
+
+    // Vector with reached levels
+    sp_v = Vector::Make(n, Types::Int32(library), library);
+    // Vector-front of the bfs
+    auto sp_q = Vector::Make(n, Types::Void(library), library);
+    // Scalar to update depth of the v
+    auto sp_depth = Scalar::Make(Types::Int32(library), library);
+
+    // Used to assign to v new reached level and preserve v values
+    auto sp_desc_accum = Descriptor::Make(library);
+    sp_desc_accum->SetParam(Descriptor::Param::AccumResult);
+
+    // Used to apply !v mask when try to discover new values and ignore previously found
+    auto sp_desc_comp = Descriptor::Make(library);
+    sp_desc_comp->SetParam(Descriptor::Param::MaskComplement);
+
+    // Set q[s]: start vertex
+    auto sp_setup = Expression::Make(library);
+    sp_setup->MakeDataWrite(sp_q, DataVector::Make(&s, nullptr, 1, library));
+    sp_setup->SubmitWait();
+
+    // Start for depth 1: v[s]=1
+    std::int32_t depth = 1;
+
+    while (sp_q->GetNvals() != 0) {
+        auto sp_iter = Expression::Make(library);
+
+        auto t1 = sp_iter->MakeDataWrite(sp_depth, DataScalar::Make(&depth, library));     // Update depth scalar
+        auto t2 = sp_iter->MakeAssign(sp_v, sp_q, nullptr, sp_depth, sp_desc_accum);       // New reached vertices v[q] = depth
+        auto t3 = sp_iter->MakeVxM(sp_q, sp_v, nullptr, nullptr, sp_q, sp_A, sp_desc_comp);// Discover new front q[!v] = q x A
+
+        sp_iter->Dependency(t1, t2);
+        sp_iter->Dependency(t2, t3);
+        sp_iter->SubmitWait();
+
+        depth += 1;
+    }
+}
+
+void spla::Bfs(RefPtr<HostVector> &v, const RefPtr<HostMatrix> &A, Index s) {
+    CHECK_RAISE_ERROR(A.IsNotNull(), NullPointer, "Passed null argument");
+    CHECK_RAISE_ERROR(A->GetNrows() == A->GetNcols(), DimensionMismatch, "Matrix must be nxn");
+    CHECK_RAISE_ERROR(s < A->GetNrows(), InvalidArgument, "Start index must be withing A bounds");
+
+    auto n = A->GetNrows();
+
+    // Offset to fetch rows by index
+    std::vector<Index> offsets;
+    // Depths of reached vertices, initially - +inf
+    std::vector<std::int32_t> depths(n, std::numeric_limits<std::int32_t>::max());
+    // Track already reached vertices
+    std::vector<bool> visited(n, false);
+    // Queue of the bfs front
+    std::queue<Index> front;
+
+    // Compute row offsets
+    IndicesToRowOffsets(A->GetRowIndices(), offsets, n);
+
+    // Start from s
+    front.push(s);
+    depths[s] = 1;
+    visited[s] = true;
+
+    while (!front.empty()) {
+        // Get next vertex and remove it from the queue
+        auto i = front.front();
+        front.pop();
+
+        // Traverse all adjacent neighbors
+        for (Index k = offsets[i]; k < offsets[i + 1]; k++) {
+            auto j = A->GetColIndices()[k];
+
+            // Reached new vertex
+            if (!visited[j]) {
+                // Update depth
+                depths[j] = depths[i] + 1;
+                // Mark as visited
+                visited[j] = true;
+                // Add to the front
+                front.push(j);
+            }
+        }
+    }
+
+    // Define reached vertices
+    std::vector<Index> rows;
+    std::vector<std::int32_t> values;
+
+    // Add in order only reached vertices
+    for (Index k = 0; k < n; k++) {
+        if (visited[k]) {
+            rows.push_back(k);
+            values.push_back(depths[k]);
+        }
+    }
+
+    // Convert to raw data
+    std::vector<unsigned char> data(values.size() * sizeof(std::int32_t));
+    std::memcpy(data.data(), values.data(), values.size() * sizeof(std::int32_t));
+
+    // Build result vector
+    v = RefPtr<HostVector>(new HostVector(n, std::move(rows), std::move(data)));
+}
\ No newline at end of file
diff --git a/sources/SplaAlgoCommon.cpp b/sources/SplaAlgoCommon.cpp
new file mode 100644
index 000000000..7133d6da0
--- /dev/null
+++ b/sources/SplaAlgoCommon.cpp
@@ -0,0 +1,44 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#include <spla-algo/SplaAlgoCommon.hpp>
+
+spla::HostVector::HostVector(spla::Size nrows, std::vector<Index> rows, std::vector<unsigned char> vals)
+    : mNrows(nrows), mNnvals(rows.size()), mElementSize(vals.size() / rows.size()), mRowIndices(std::move(rows)), mValues(std::move(vals)) {
+}
+
+spla::RefPtr<spla::DataVector> spla::HostVector::GetData(Library &library) {
+    return DataVector::Make(mRowIndices.data(), HasValues() ? mValues.data() : nullptr, GetNnvals(), library);
+}
+
+spla::HostMatrix::HostMatrix(spla::Size nrows, spla::Size ncols, std::vector<Index> rows, std::vector<Index> cols, std::vector<unsigned char> vals)
+    : mNrows(nrows), mNcols(ncols), mNnvals(rows.size()), mElementSize(vals.size() / rows.size()), mRowIndices(std::move(rows)), mColIndices(std::move(cols)), mValues(std::move(vals)) {
+}
+
+spla::RefPtr<spla::DataMatrix> spla::HostMatrix::GetData(Library &library) {
+    return DataMatrix::Make(mRowIndices.data(), mColIndices.data(), HasValues() ? mValues.data() : nullptr, GetNnvals(), library);
+}
\ No newline at end of file
diff --git a/sources/SplaData.cpp b/sources/SplaData.cpp
index 0241777e2..cec6e5cd6 100644
--- a/sources/SplaData.cpp
+++ b/sources/SplaData.cpp
@@ -89,6 +89,15 @@ spla::RefPtr<spla::DataMatrix> spla::DataMatrix::Make(spla::Library &library) {
     return spla::RefPtr<spla::DataMatrix>(new DataMatrix(library));
 }
 
+spla::RefPtr<spla::DataMatrix> spla::DataMatrix::Make(unsigned int *rows, unsigned int *cols, void *values, std::size_t nvals, Library &library) {
+    auto data = Make(library);
+    data->SetRows(rows);
+    data->SetCols(cols);
+    data->SetVals(values);
+    data->SetNvals(nvals);
+    return data;
+}
+
 spla::DataVector::DataVector(spla::Library &library) : Data(Object::TypeName::DataVector, library) {
 }
 
@@ -133,6 +142,14 @@ spla::RefPtr<spla::DataVector> spla::DataVector::Make(spla::Library &library) {
     return spla::RefPtr<spla::DataVector>(new DataVector(library));
 }
 
+spla::RefPtr<spla::DataVector> spla::DataVector::Make(unsigned int *rows, void *values, std::size_t nvals, Library &library) {
+    auto data = Make(library);
+    data->SetRows(rows);
+    data->SetVals(values);
+    data->SetNvals(nvals);
+    return data;
+}
+
 spla::DataScalar::DataScalar(spla::Library &library) : Data(Object::TypeName::DataScalar, library) {
 }
 
@@ -156,3 +173,9 @@ void *spla::DataScalar::GetValue() const {
 spla::RefPtr<spla::DataScalar> spla::DataScalar::Make(spla::Library &library) {
     return spla::RefPtr<spla::DataScalar>(new DataScalar(library));
 }
+
+spla::RefPtr<spla::DataScalar> spla::DataScalar::Make(void *value, spla::Library &library) {
+    auto data = Make(library);
+    data->SetValue(value);
+    return data;
+}
\ No newline at end of file
diff --git a/sources/SplaExpression.cpp b/sources/SplaExpression.cpp
index 15c68331b..4669add48 100644
--- a/sources/SplaExpression.cpp
+++ b/sources/SplaExpression.cpp
@@ -63,6 +63,11 @@ void spla::Expression::Wait() {
         mFuture->Get().wait();
 }
 
+void spla::Expression::SubmitWait() {
+    Submit();
+    Wait();
+}
+
 void spla::Expression::Dependency(const spla::RefPtr<spla::ExpressionNode> &pred,
                                   const spla::RefPtr<spla::ExpressionNode> &succ) {
     CHECK_RAISE_ERROR(GetState() == State::Default, InvalidState, "Expression must be in default state");
diff --git a/sources/compute/SplaIndicesToRowOffsets.hpp b/sources/compute/SplaIndicesToRowOffsets.hpp
index 0ef3ba7ec..e7e8c8455 100644
--- a/sources/compute/SplaIndicesToRowOffsets.hpp
+++ b/sources/compute/SplaIndicesToRowOffsets.hpp
@@ -31,6 +31,10 @@
 #include <boost/compute/algorithm.hpp>
 #include <boost/compute/command_queue.hpp>
 
+#include <algorithm>
+#include <numeric>
+#include <vector>
+
 namespace spla {
 
     /**
@@ -90,6 +94,47 @@ namespace spla {
         IndicesToRowOffsets(indices, offsets, lengths, n, queue);
     }
 
+    /**
+     * @brief Compute row offsets from coo row indices buffer.
+     *
+     * Result offsets array has size n + 1.
+     * Offsets[i] stores first index of row i in indices buffer.
+     * Offsets[i + 1] - Offsets[i] equals number of nnz values in row i in indices buffer.
+     * Offsets[n] equals number off values in indices buffer.
+     * Lengths[i] equals number of nnz values in row i in indices buffer.
+     *
+     * @param indices Array of rows indices
+     * @param[out] offsets Output array with offsets
+     * @param[out] lengths Output array of row lengths
+     * @param n Number of rows in matrix
+     */
+    inline void IndicesToRowOffsets(const std::vector<unsigned int> &indices,
+                                    std::vector<Index> &offsets,
+                                    std::vector<Index> &lengths,
+                                    std::size_t n) {
+        offsets.resize(n + 1);
+        lengths.resize(n + 1);
+        std::fill(lengths.begin(), lengths.end(), 0u);
+
+        if (indices.empty()) {
+            std::fill(offsets.begin(), offsets.end(), 0u);
+            return;
+        }
+
+        std::for_each(indices.begin(), indices.end(), [&](auto i) { lengths[i] += 1; });
+        std::exclusive_scan(lengths.begin(), lengths.end(), offsets.begin(), 0u);
+    }
+
+    /**
+     * Overload
+     */
+    inline void IndicesToRowOffsets(const std::vector<unsigned int> &indices,
+                                    std::vector<Index> &offsets,
+                                    std::size_t n) {
+        std::vector<unsigned int> lengths;
+        IndicesToRowOffsets(indices, offsets, lengths, n);
+    }
+
     /**
      * @}
      */
diff --git a/sources/expression/vector/SplaVectorAssign.cpp b/sources/expression/vector/SplaVectorAssign.cpp
index 2c4725adf..1afee182b 100644
--- a/sources/expression/vector/SplaVectorAssign.cpp
+++ b/sources/expression/vector/SplaVectorAssign.cpp
@@ -94,7 +94,7 @@ void spla::VectorAssign::Process(std::size_t nodeIdx, const spla::Expression &ex
 
         if (applyAccum) {
             auto accumTask = builder.Emplace([=]() {
-                auto tmpBlock = w->GetStorage()->GetBlock(i);
+                auto tmpBlock = tmp->GetStorage()->GetBlock(i);
 
                 if (tmpBlock.IsNotNull()) {
                     ParamsVectorEWiseAdd params;


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index bc5736620..fb72fed2e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -6,6 +6,7 @@ function(spla_test_target target)
     target_include_directories(${target} PRIVATE ../sources)
 endfunction()
 
+spla_test_target(TestAlgoBfs)
 spla_test_target(TestBasic)
 spla_test_target(TestDataMatrix)
 spla_test_target(TestDataScalar)
diff --git a/tests/TestAlgoBfs.cpp b/tests/TestAlgoBfs.cpp
new file mode 100644
index 000000000..5591e668e
--- /dev/null
+++ b/tests/TestAlgoBfs.cpp
@@ -0,0 +1,84 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#include <Testing.hpp>
+
+void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::size_t seed = 0) {
+    auto rnd = utils::UniformIntGenerator<spla::Index>(seed, 0, M - 1);
+    auto sp_Int32 = spla::Types::Int32(library);
+    auto sp_Void = spla::Types::Void(library);
+    auto sp_s = rnd();
+
+    utils::Matrix A = utils::Matrix<char>::Generate(M, M, nvals).SortReduceDuplicates();
+
+    auto sp_v = spla::Vector::Make(M, sp_Int32, library);
+    auto sp_A = spla::Matrix::Make(M, M, sp_Void, library);
+
+    auto sp_setup = spla::Expression::Make(library);
+    sp_setup->MakeDataWrite(sp_A, A.GetDataIndices(library));
+    sp_setup->SubmitWait();
+    ASSERT_EQ(sp_setup->GetState(), spla::Expression::State::Evaluated);
+
+    spla::Bfs(sp_v, sp_A, sp_s);
+
+    auto host_A = A.ToHostMatrix();
+    auto host_v = spla::RefPtr<spla::HostVector>();
+
+    spla::Bfs(host_v, host_A, sp_s);
+
+    auto result = utils::Vector<std::int32_t>::FromHostVector(host_v);
+    ASSERT_TRUE(result.Equals(sp_v));
+}
+
+void test(std::size_t M, std::size_t base, std::size_t step, std::size_t iter, const std::vector<std::size_t> &blocksSizes) {
+    utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCase(library, M, nvals, i);
+        }
+    });
+}
+
+TEST(BFS, Small) {
+    std::vector<std::size_t> blockSizes = {100, 1000};
+    std::size_t M = 120;
+    test(M, M, M, 10, blockSizes);
+}
+
+TEST(BFS, Medium) {
+    std::vector<std::size_t> blockSizes = {1000, 10000};
+    std::size_t M = 1220;
+    test(M, M, M, 10, blockSizes);
+}
+
+TEST(BFS, Large) {
+    std::vector<std::size_t> blockSizes = {10000, 100000};
+    std::size_t M = 12400;
+    test(M, M, M, 5, blockSizes);
+}
+
+SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/TestVxM.cpp b/tests/TestVxM.cpp
index 18e837b87..a45239b9b 100644
--- a/tests/TestVxM.cpp
+++ b/tests/TestVxM.cpp
@@ -209,13 +209,13 @@ TEST(VxM, Small) {
     test(M, N, M, M, 10, blockSizes);
 }
 
-TEST(MxM, Medium) {
+TEST(VxM, Medium) {
     std::vector<std::size_t> blockSizes = {1000, 10000};
     std::size_t M = 1220, N = 880;
     test(M, N, M, M, 10, blockSizes);
 }
 
-TEST(MxM, Large) {
+TEST(VxM, Large) {
     std::vector<std::size_t> blockSizes = {10000, 100000};
     std::size_t M = 12400, N = 8080;
     test(M, N, M, M, 5, blockSizes);
diff --git a/tests/utils/Matrix.hpp b/tests/utils/Matrix.hpp
index 4fe7430ea..0deeb4164 100644
--- a/tests/utils/Matrix.hpp
+++ b/tests/utils/Matrix.hpp
@@ -33,6 +33,7 @@
 #include <functional>
 #include <iostream>
 #include <random>
+#include <spla-algo/SplaAlgo.hpp>
 #include <spla-cpp/Spla.hpp>
 #include <utility>
 #include <utils/Compute.hpp>
@@ -495,6 +496,22 @@ namespace utils {
             return MxM<W>(b, multOp, addOp).Mask(mask, complement);
         }
 
+        [[nodiscard]] spla::RefPtr<spla::HostMatrix> ToHostMatrix() const {
+            std::vector<Index> rows = mRows;
+            std::vector<Index> cols = mCols;
+            std::vector<unsigned char> data(mVals.size() * sizeof(T));
+            std::memcpy(data.data(), mVals.data(), mVals.size() * sizeof(T));
+            return spla::RefPtr<spla::HostMatrix>(new spla::HostMatrix(mNrows, mNcols, std::move(rows), std::move(cols), std::move(data)));
+        }
+
+        static Matrix FromHostMatrix(const spla::RefPtr<spla::HostMatrix> &m) {
+            std::vector<Index> rows = m->GetRowIndices();
+            std::vector<Index> cols = m->GetColIndices();
+            std::vector<T> values(m->GetNnvals());
+            std::memcpy(values.data(), m->GetValues().data(), values.size() * sizeof(T));
+            return Matrix(m->GetNrows(), m->GetNcols(), std::move(rows), std::move(cols), std::move(values));
+        }
+
         static Matrix Empty(std::size_t nrows, std::size_t ncols) {
             return Matrix(nrows, ncols, std::vector<Index>{}, std::vector<Index>{}, std::vector<T>{});
         }
diff --git a/tests/utils/Vector.hpp b/tests/utils/Vector.hpp
index 6b5e30eca..28e6f6159 100644
--- a/tests/utils/Vector.hpp
+++ b/tests/utils/Vector.hpp
@@ -202,8 +202,10 @@ namespace utils {
                 auto a = GetVals()[i];
                 auto b = spVals[i];
 
-                if ((useError && std::abs(a - b) > error) || (!useError && a != b))
+                if ((useError && std::abs(a - b) > error) || (!useError && a != b)) {
+                    std::cout << "Values not equal at " << i << " a=" << a << " b=" << b << std::endl;
                     return false;
+                }
             }
 
             return true;
@@ -372,6 +374,20 @@ namespace utils {
             return Vector<T>(nrows, std::move(rows), std::move(vals));
         }
 
+        [[nodiscard]] spla::RefPtr<spla::HostVector> ToHostVector() const {
+            std::vector<Index> rows = mRows;
+            std::vector<unsigned char> data(mVals.size() * sizeof(T));
+            std::memcpy(data.data(), mVals.data(), mVals.size() * sizeof(T));
+            return spla::RefPtr<spla::HostVector>(new spla::HostVector(mNrows, std::move(rows), std::move(data)));
+        }
+
+        static Vector FromHostVector(const spla::RefPtr<spla::HostVector> &m) {
+            std::vector<Index> rows = m->GetRowIndices();
+            std::vector<T> values(m->GetNnvals());
+            std::memcpy(values.data(), m->GetValues().data(), values.size() * sizeof(T));
+            return Vector(m->GetNrows(), std::move(rows), std::move(values));
+        }
+
         static Vector Empty(std::size_t nrows) {
             return Vector(nrows, std::vector<Index>{}, std::vector<T>{});
         }