Repository: SparseLinearAlgebra/spla (language: python)

Issue:
Algorithms: SSSP using library API



Discussion / hints:

Fix patch:
diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt
index 71c884f96..7fb638c3c 100644
--- a/include/CMakeLists.txt
+++ b/include/CMakeLists.txt
@@ -21,4 +21,5 @@ set(SPLA_CPP_HEADERS
 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
+        include/spla-algo/SplaAlgoCommon.hpp
+        include/spla-algo/SplaAlgoSssp.hpp)
\ No newline at end of file
diff --git a/include/spla-algo/SplaAlgo.hpp b/include/spla-algo/SplaAlgo.hpp
index 66aeadf97..cfb006d47 100644
--- a/include/spla-algo/SplaAlgo.hpp
+++ b/include/spla-algo/SplaAlgo.hpp
@@ -50,5 +50,6 @@
 
 #include <spla-algo/SplaAlgoBfs.hpp>
 #include <spla-algo/SplaAlgoCommon.hpp>
+#include <spla-algo/SplaAlgoSssp.hpp>
 
 #endif//SPLA_SPLAALGO_HPP
diff --git a/include/spla-algo/SplaAlgoSssp.hpp b/include/spla-algo/SplaAlgoSssp.hpp
new file mode 100644
index 000000000..37ffe7d12
--- /dev/null
+++ b/include/spla-algo/SplaAlgoSssp.hpp
@@ -0,0 +1,71 @@
+/**********************************************************************************/
+/* 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_SPLAALGOSSSP_HPP
+#define SPLA_SPLAALGOSSSP_HPP
+
+#include <spla-algo/SplaAlgoCommon.hpp>
+#include <spla-cpp/SplaMatrix.hpp>
+#include <spla-cpp/SplaVector.hpp>
+
+namespace spla {
+
+    /**
+     * @addtogroup Algorithm
+     * @{
+     */
+
+    /**
+     * @brief Single-source shortest path
+     *
+     * @note Accepts input matrix with float32 type.
+     *
+     * @param[out] v Vector where to store shortest paths 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 paths search
+     */
+    SPLA_API void Sssp(RefPtr<Vector> &v, const RefPtr<Matrix> &A, Index s);
+
+    /**
+     * @brief Single-source shortest path
+     *
+     * @note Naive cpu reference algo implementation for correctness only
+     * @note Accepts input matrix with float32 type.
+     *
+     * @param[out] v Vector where to store shortest paths 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 paths search
+     */
+    SPLA_API void Sssp(RefPtr<HostVector> &v, const RefPtr<HostMatrix> &A, Index s);
+
+    /**
+     * @}
+     */
+
+}// namespace spla
+
+#endif//SPLA_SPLAALGOSSSP_HPP
diff --git a/sources/CMakeLists.txt b/sources/CMakeLists.txt
index e86b2b6a3..6e46905d4 100644
--- a/sources/CMakeLists.txt
+++ b/sources/CMakeLists.txt
@@ -17,7 +17,8 @@ set(SPLA_CPP_SOURCES
 
 set(SPLA_ALGO_SOURCES
         sources/SplaAlgoBfs.cpp
-        sources/SplaAlgoCommon.cpp)
+        sources/SplaAlgoCommon.cpp
+        sources/SplaAlgoSssp.cpp)
 
 set(SPLA_ALGORITHM_SOURCES
         sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp
diff --git a/sources/SplaAlgoSssp.cpp b/sources/SplaAlgoSssp.cpp
new file mode 100644
index 000000000..55bdb51c0
--- /dev/null
+++ b/sources/SplaAlgoSssp.cpp
@@ -0,0 +1,162 @@
+/**********************************************************************************/
+/* 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/SplaAlgoSssp.hpp>
+#include <spla-cpp/Spla.hpp>
+
+#include <compute/SplaIndicesToRowOffsets.hpp>
+#include <core/SplaError.hpp>
+
+#include <limits>
+#include <queue>
+#include <vector>
+
+void spla::Sssp(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 type = sp_A->GetType();
+    auto n = sp_A->GetNrows();
+
+    CHECK_RAISE_ERROR(type->IsBuiltIn(), InvalidArgument, "Must be built-in type");
+    CHECK_RAISE_ERROR(type == Types::Float32(library), InvalidArgument, "Must be float32 type");
+
+    float zero = 0.0f;
+
+    // Used to check if something changed
+    auto reduceOp = Functions::PlusFloat32(library);
+    // Used to concatenate paths lengths
+    auto multOp = Functions::PlusFloat32(library);
+    // Used find shortest path almond others
+    auto addOp = Functions::MinFloat32(library);
+
+    // Vector with reached paths
+    sp_v = Vector::Make(n, type, library);
+    // Temporary to store new reached paths
+    auto sp_w = Vector::Make(n, type, library);
+    // Scalar to store reduce succ
+    auto sp_succ = Scalar::Make(type, library);
+
+    // Set v[s] = w[s] = 0.0f - so start is reached in 0.0
+    auto sp_setup = Expression::Make(library);
+    sp_setup->MakeDataWrite(sp_v, DataVector::Make(&s, &zero, 1, library));
+    sp_setup->MakeDataWrite(sp_w, DataVector::Make(&s, &zero, 1, library));
+    sp_setup->SubmitWait();
+
+    std::size_t iterations = 0;
+    float succLast = 0.0f;
+    float succ = 1.0f;
+
+    while (iterations < n && succ != succLast) {
+        succLast = succ;
+
+        auto sp_iter = Expression::Make(library);
+
+        auto t1 = sp_iter->MakeVxM(sp_w, nullptr, multOp, addOp, sp_w, sp_A);
+        auto t2 = sp_iter->MakeEWiseAdd(sp_v, nullptr, addOp, sp_v, sp_w);
+        auto t3 = sp_iter->MakeReduce(sp_succ, reduceOp, sp_v);
+        auto t4 = sp_iter->MakeDataRead(sp_succ, DataScalar::Make(&succ, library));
+
+        sp_iter->Dependency(t1, t2);
+        sp_iter->Dependency(t2, t3);
+        sp_iter->Dependency(t3, t4);
+        sp_iter->SubmitWait();
+
+        iterations += 1;
+    }
+}
+
+void spla::Sssp(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");
+    CHECK_RAISE_ERROR(A->GetElementSize() == sizeof(float), InvalidType, "Matrix A must have float values of size 4");
+
+    auto n = A->GetNrows();
+    auto inf = std::numeric_limits<float>::max();
+
+    // Note: implementation of shortest path fastest algorithms
+    // see for details https://en.wikipedia.org/wiki/Shortest_Path_Faster_Algorithm
+
+    // Offset to fetch rows by index
+    std::vector<Index> offsets;
+    // Initial distances to vertices
+    std::vector<float> distances(n, inf);
+    // Track already reached vertices
+    std::vector<bool> isInQueue(n, false);
+    // Queue of the search
+    std::queue<Index> front;
+
+    // Compute row offsets
+    IndicesToRowOffsets(A->GetRowIndices(), offsets, n);
+
+    // Start from s
+    front.push(s);
+    distances[s] = 0.0f;
+    isInQueue[s] = true;
+
+    while (!front.empty()) {
+        auto i = front.front();
+        front.pop();
+        isInQueue[i] = false;
+
+        // Traverse all adjacent neighbors
+        for (Index k = offsets[i]; k < offsets[i + 1]; k++) {
+            auto j = A->GetColIndices()[k];
+            auto w = *(reinterpret_cast<const float *>(A->GetValues().data()) + k);
+
+            if (distances[j] == inf || distances[i] + w < distances[j]) {
+                distances[j] = distances[i] + w;
+                if (!isInQueue[j]) {
+                    isInQueue[j] = true;
+                    front.push(j);
+                }
+            }
+        }
+    }
+
+    // Define reached vertices
+    std::vector<Index> rows;
+    std::vector<float> values;
+
+    // Add in order only reached vertices
+    for (Index k = 0; k < n; k++) {
+        if (distances[k] != inf) {
+            rows.push_back(k);
+            values.push_back(distances[k]);
+        }
+    }
+
+    // Convert to raw data
+    std::vector<unsigned char> data(values.size() * sizeof(float));
+    std::memcpy(data.data(), values.data(), values.size() * sizeof(float));
+
+    // 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/expression/vector/SplaVectorReduce.cpp b/sources/expression/vector/SplaVectorReduce.cpp
index 3f13aad83..0a0d70f0a 100644
--- a/sources/expression/vector/SplaVectorReduce.cpp
+++ b/sources/expression/vector/SplaVectorReduce.cpp
@@ -114,21 +114,24 @@ void spla::VectorReduce::Process(std::size_t nodeIdx, const spla::Expression &ex
 
     for (std::size_t i = 0; i < blocksInVector; ++i) {
         auto deviceId = deviceIds[i];
-
-        tf::Task reduceIthBlock = builder.Emplace([=]() {
-            ParamsVectorReduce params;
-            params.desc = desc;
-            params.deviceId = deviceId;
-            params.vec = argV->GetStorage()->GetBlock(i);
-            params.type = argV->GetType();
-            params.reduce = argOp;
-            library->GetAlgoManager()->Dispatch(Algorithm::Type::VectorReduce, params);
-            if (params.scalar.IsNotNull()) {
-                intermediateBuffer->Add(params.scalar);
-            }
-            SPDLOG_LOGGER_TRACE(logger, "Reduce block i={} nnz={}", i, params.vec->GetNvals());
-        });
-        reduceBlocksTasks.push_back(std::move(reduceIthBlock));
+        auto blockVec = argV->GetStorage()->GetBlock(i);
+
+        if (blockVec.IsNotNull()) {
+            tf::Task reduceIthBlock = builder.Emplace([=]() {
+                ParamsVectorReduce params;
+                params.desc = desc;
+                params.deviceId = deviceId;
+                params.vec = argV->GetStorage()->GetBlock(i);
+                params.type = argV->GetType();
+                params.reduce = argOp;
+                library->GetAlgoManager()->Dispatch(Algorithm::Type::VectorReduce, params);
+                if (params.scalar.IsNotNull()) {
+                    intermediateBuffer->Add(params.scalar);
+                }
+                SPDLOG_LOGGER_TRACE(logger, "Reduce block i={} nnz={}", i, params.vec->GetNvals());
+            });
+            reduceBlocksTasks.push_back(std::move(reduceIthBlock));
+        }
     }
 
     auto lastReduceDeviceId = deviceIds[blocksInVector];


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index ae8cbc541..b84ca8b70 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -7,6 +7,7 @@ function(spla_test_target target)
 endfunction()
 
 spla_test_target(TestAlgoBfs)
+spla_test_target(TestAlgoSssp)
 spla_test_target(TestBasic)
 spla_test_target(TestDataMatrix)
 spla_test_target(TestDataScalar)
diff --git a/tests/TestAlgoBfs.cpp b/tests/TestAlgoBfs.cpp
index 5f3b4764d..4b35e4b8d 100644
--- a/tests/TestAlgoBfs.cpp
+++ b/tests/TestAlgoBfs.cpp
@@ -43,12 +43,16 @@ void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::siz
     sp_setup->SubmitWait();
     ASSERT_EQ(sp_setup->GetState(), spla::Expression::State::Evaluated);
 
+    SPLA_TIME_BEGIN(bfs_spla);
     spla::Bfs(sp_v, sp_A, sp_s);
+    SPLA_TIME_END(bfs_spla, "spla");
 
     auto host_A = A.ToHostMatrix();
     auto host_v = spla::RefPtr<spla::HostVector>();
 
+    SPLA_TIME_BEGIN(bfs_cpu);
     spla::Bfs(host_v, host_A, sp_s);
+    SPLA_TIME_END(bfs_cpu, "cpu");
 
     auto result = utils::Vector<std::int32_t>::FromHostVector(host_v);
     ASSERT_TRUE(result.Equals(sp_v));
@@ -81,4 +85,16 @@ TEST(BFS, Large) {
     test(M, M, M, 5, blockSizes);
 }
 
+TEST(BFS, MegaLarge) {
+    std::vector<std::size_t> blockSizes = {1000000};
+    std::size_t M = 990990;
+    test(M, 10 * M, M, 1, blockSizes);
+}
+
+TEST(BFS, UltraLarge) {
+    std::vector<std::size_t> blockSizes = {10000000};
+    std::size_t M = 4500000;
+    test(M, 10 * M, 10 * M, 2, blockSizes);
+}
+
 SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/TestAlgoSssp.cpp b/tests/TestAlgoSssp.cpp
new file mode 100644
index 000000000..d9d491dde
--- /dev/null
+++ b/tests/TestAlgoSssp.cpp
@@ -0,0 +1,94 @@
+/**********************************************************************************/
+/* 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_Float32 = spla::Types::Float32(library);
+    auto sp_s = rnd();
+
+    utils::Matrix A = utils::Matrix<float>::Generate(M, M, nvals).SortReduceDuplicates();
+    A.Fill(utils::UniformRealGenerator<float>());
+
+    auto sp_v = spla::Vector::Make(M, sp_Float32, library);
+    auto sp_A = spla::Matrix::Make(M, M, sp_Float32, library);
+
+    auto sp_setup = spla::Expression::Make(library);
+    sp_setup->MakeDataWrite(sp_A, A.GetData(library));
+    sp_setup->SubmitWait();
+    ASSERT_EQ(sp_setup->GetState(), spla::Expression::State::Evaluated);
+
+    SPLA_TIME_BEGIN(sssp_spla);
+    spla::Sssp(sp_v, sp_A, sp_s);
+    SPLA_TIME_END(sssp_spla, "spla");
+
+    auto host_A = A.ToHostMatrix();
+    auto host_v = spla::RefPtr<spla::HostVector>();
+
+    SPLA_TIME_BEGIN(sssp_cpu);
+    spla::Sssp(host_v, host_A, sp_s);
+    SPLA_TIME_END(sssp_cpu, "cpu");
+
+    auto result = utils::Vector<float>::FromHostVector(host_v);
+    EXPECT_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(SSSP, Small) {
+    std::vector<std::size_t> blockSizes = {100, 1000};
+    std::size_t M = 120;
+    test(M, M, M, 10, blockSizes);
+}
+
+TEST(SSSP, Medium) {
+    std::vector<std::size_t> blockSizes = {1000, 10000};
+    std::size_t M = 1220;
+    test(M, M, M, 10, blockSizes);
+}
+
+TEST(SSSP, Large) {
+    std::vector<std::size_t> blockSizes = {10000, 100000};
+    std::size_t M = 12400;
+    test(M, M, M, 5, blockSizes);
+}
+
+TEST(SSSP, MegaLarge) {
+    std::vector<std::size_t> blockSizes = {1000000};
+    std::size_t M = 990990;
+    test(M, 5 * M, M, 3, blockSizes);
+}
+
+SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/Testing.hpp b/tests/Testing.hpp
index 2b21eb0c3..61c7d09f1 100644
--- a/tests/Testing.hpp
+++ b/tests/Testing.hpp
@@ -37,6 +37,16 @@
 #include <utils/Typetraits.hpp>
 #include <utils/Vector.hpp>
 
+#include <chrono>
+#include <iostream>
+
+#define SPLA_TIME_BEGIN(timer) \
+    auto timer = std::chrono::steady_clock::now();
+
+#define SPLA_TIME_END(timer, message) \
+    std::cout << message << " " << std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - timer).count() << std::endl;
+
+
 // Put in the end of the unit test file
 #define SPLA_GTEST_MAIN                                  \
     int main(int argc, char *argv[]) {                   \