From b8bd9a8a5e4f7eca8aae2ed46ec1af9415ff5c19 Mon Sep 17 00:00:00 2001
From: Alex Bilger <alexbilger0@gmail.com>
Date: Fri, 20 Oct 2023 10:55:28 +0200
Subject: [PATCH 1/2] Measure perfs of BTDMatrix::add

---
 CMakeLists.txt                                |  1 +
 .../Sofa.LinearAlgebra/BTDMatrix.cpp          | 53 +++++++++++++++++++
 2 files changed, 54 insertions(+)
 create mode 100644 src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7f6adb8..8416a78 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -36,6 +36,7 @@ list(APPEND HEADER_FILES
     ${SOFABENCHMARK_SRC}/utils/thread_pool.hpp
 )
 list(APPEND SOURCE_FILES
+    ${SOFABENCHMARK_SRC}/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp
     ${SOFABENCHMARK_SRC}/benchmarks/Sofa.LinearAlgebra/SparseMatrixCompression.cpp
     ${SOFABENCHMARK_SRC}/benchmarks/Sofa.LinearAlgebra/SparseMatrixMulTranspose.cpp
     ${SOFABENCHMARK_SRC}/benchmarks/Sofa.LinearAlgebra/SparseMatrixProduct.cpp
diff --git a/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp b/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp
new file mode 100644
index 0000000..969d922
--- /dev/null
+++ b/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp
@@ -0,0 +1,53 @@
+#include <benchmark/benchmark.h>
+
+#include <sofa/linearalgebra/BTDMatrix.h>
+#include <utils/RandomValuePool.h>
+
+constexpr int64_t maxSubIterations = 8 << 7;
+
+template<std::size_t N, typename T>
+static void BM_BTDMatrix_add(benchmark::State& state)
+{
+    sofa::linearalgebra::BTDMatrix<N, T> matrix;
+
+    constexpr auto totalsize = maxSubIterations * N * N;
+    const std::array<T, totalsize>& values = RandomValuePool<T, totalsize>::get();
+
+    std::vector<sofa::type::Mat<N, N, T>> matricesToAdd;
+    matricesToAdd.resize(state.range(0));
+    for (int64_t i = 0; i < state.range(0); i++)
+    {
+        for (sofa::Index j = 0; j < N; ++j)
+        {
+            for (sofa::Index k = 0; k < N; ++k)
+            {
+                matricesToAdd[i][j][k] = values[i * N * N + j * N + k];
+            }
+        }
+    }
+
+    auto randIndices = RandomValuePool<sofa::Index, maxSubIterations * 2>::get();
+    for (auto& id : randIndices)
+    {
+        id *= N;
+    }
+    std::vector<std::pair<sofa::Index, sofa::Index>> indices;
+    indices.reserve(state.range(0));
+    for (int64_t i = 0; i < state.range(0) * 2; i += 2)
+    {
+        indices.emplace_back(randIndices[i], randIndices[i+1]);
+    }
+
+    const auto max = *std::max_element(randIndices.begin(), randIndices.end());
+    matrix.resize(N * ((max+N) / N), N * ((max+N) / N));
+
+    for (auto _ : state)
+    {
+        for (int64_t i = 0; i < state.range(0); ++i)
+        {
+            matrix.add(indices[i].first, indices[i].second, matricesToAdd[i]);
+        }
+    }
+}
+
+BENCHMARK_TEMPLATE(BM_BTDMatrix_add, 6, SReal)->Range(8 << 3, maxSubIterations);

From e81033dcc650093e25c5f6fafe3038ee6ccc0a8a Mon Sep 17 00:00:00 2001
From: Alex Bilger <alexbilger0@gmail.com>
Date: Fri, 20 Oct 2023 11:07:12 +0200
Subject: [PATCH 2/2] comparison with the scalar accumulation

---
 .../Sofa.LinearAlgebra/BTDMatrix.cpp          | 42 ++++++++++++++++---
 1 file changed, 36 insertions(+), 6 deletions(-)

diff --git a/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp b/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp
index 969d922..c560047 100644
--- a/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp
+++ b/src/benchmarks/Sofa.LinearAlgebra/BTDMatrix.cpp
@@ -5,15 +5,13 @@
 
 constexpr int64_t maxSubIterations = 8 << 7;
 
-template<std::size_t N, typename T>
-static void BM_BTDMatrix_add(benchmark::State& state)
-{
-    sofa::linearalgebra::BTDMatrix<N, T> matrix;
 
+template <std::size_t N, typename T>
+static void prepareData(benchmark::State& state, sofa::linearalgebra::BTDMatrix<N, T>& matrix, std::vector<sofa::type::Mat<N, N, T>>& matricesToAdd, std::vector<std::pair<sofa::Index, sofa::Index>>& indices)
+{
     constexpr auto totalsize = maxSubIterations * N * N;
     const std::array<T, totalsize>& values = RandomValuePool<T, totalsize>::get();
 
-    std::vector<sofa::type::Mat<N, N, T>> matricesToAdd;
     matricesToAdd.resize(state.range(0));
     for (int64_t i = 0; i < state.range(0); i++)
     {
@@ -31,7 +29,6 @@ static void BM_BTDMatrix_add(benchmark::State& state)
     {
         id *= N;
     }
-    std::vector<std::pair<sofa::Index, sofa::Index>> indices;
     indices.reserve(state.range(0));
     for (int64_t i = 0; i < state.range(0) * 2; i += 2)
     {
@@ -40,6 +37,38 @@ static void BM_BTDMatrix_add(benchmark::State& state)
 
     const auto max = *std::max_element(randIndices.begin(), randIndices.end());
     matrix.resize(N * ((max+N) / N), N * ((max+N) / N));
+}
+
+template<std::size_t N, typename T>
+static void BM_BTDMatrix_add(benchmark::State& state)
+{
+    sofa::linearalgebra::BTDMatrix<N, T> matrix;
+    std::vector<sofa::type::Mat<N, N, T>> matricesToAdd;
+    std::vector<std::pair<sofa::Index, sofa::Index>> indices;
+    prepareData<N, T>(state, matrix, matricesToAdd, indices);
+
+    for (auto _ : state)
+    {
+        for (int64_t i = 0; i < state.range(0); ++i)
+        {
+            for (sofa::Index j = 0; j < N; ++j)
+            {
+                for (sofa::Index k = 0; k < N; ++k)
+                {
+                    matrix.add(indices[i].first, indices[i].second, matricesToAdd[i](j, k));
+                }
+            }
+        }
+    }
+}
+
+template<std::size_t N, typename T>
+static void BM_BTDMatrix_addBlock(benchmark::State& state)
+{
+    sofa::linearalgebra::BTDMatrix<N, T> matrix;
+    std::vector<sofa::type::Mat<N, N, T>> matricesToAdd;
+    std::vector<std::pair<sofa::Index, sofa::Index>> indices;
+    prepareData<N, T>(state, matrix, matricesToAdd, indices);
 
     for (auto _ : state)
     {
@@ -51,3 +80,4 @@ static void BM_BTDMatrix_add(benchmark::State& state)
 }
 
 BENCHMARK_TEMPLATE(BM_BTDMatrix_add, 6, SReal)->Range(8 << 3, maxSubIterations);
+BENCHMARK_TEMPLATE(BM_BTDMatrix_addBlock, 6, SReal)->Range(8 << 3, maxSubIterations);