From 52b0077922c855ca17e2348ae04fcea2ba0e76a6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?David=20Holzm=C3=BCller?= <david.holzmueller@online.de>
Date: Thu, 11 Apr 2019 17:38:33 +0200
Subject: [PATCH] resize instead of reserve gained some efficiency

---
 src/Interpolation.hpp | 16 ++++++++++------
 src/Iterators.hpp     |  2 ++
 test/main.cpp         |  4 ++--
 3 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/src/Interpolation.hpp b/src/Interpolation.hpp
index 9fee70d..956b478 100644
--- a/src/Interpolation.hpp
+++ b/src/Interpolation.hpp
@@ -69,9 +69,10 @@ void multiply_lower_triangular_inplace(It it, std::vector<boost::numeric::ublas:
 
   for (int k = d - 1; k >= 0; --k) {
     MultiDimVector w(n[k]);
+    std::vector<size_t> indexes(n[k], 0);
+    auto sizes = it.numValuesPerFirstIndex();
     for (size_t idx = 0; idx < n[k]; ++idx) {
-      // TODO: make compatible with other iterators
-      w.data[idx].reserve(v.data[idx].size());
+      w.data[idx].resize(sizes[idx]);
     }
 
     auto &Lk = L[k];
@@ -91,7 +92,8 @@ void multiply_lower_triangular_inplace(It it, std::vector<boost::numeric::ublas:
         for (size_t j = 0; j <= i; ++j) {
           sum += Lk(i, j) * (*(offset_data_pointer + j));
         }
-        w.data[i].push_back(sum);
+        w.data[i][indexes[i]] = sum;
+        ++indexes[i];
       }
       second_v_index += last_dim_count;
       if (second_v_index >= data_size) {
@@ -124,9 +126,10 @@ void multiply_upper_triangular_inplace(It it, std::vector<boost::numeric::ublas:
 
   for (int k = d - 1; k >= 0; --k) {
     MultiDimVector w(n[k]);
+    std::vector<size_t> indexes(n[k], 0);
+    auto sizes = it.numValuesPerFirstIndex();
     for (size_t idx = 0; idx < n[k]; ++idx) {
-      // TODO: make compatible with other iterators
-      w.data[idx].reserve(v.data[idx].size());
+      w.data[idx].resize(sizes[idx]);
     }
 
     auto &Uk = U[k];
@@ -146,7 +149,8 @@ void multiply_upper_triangular_inplace(It it, std::vector<boost::numeric::ublas:
         for (size_t j = i; j < last_dim_count; ++j) {
           sum += Uk(i, j) * (*(offset_data_pointer + j));
         }
-        w.data[i].push_back(sum);
+        w.data[i][indexes[i]] = sum;
+        ++indexes[i];
       }
       second_v_index += last_dim_count;
       if (second_v_index >= data_size) {
diff --git a/src/Iterators.hpp b/src/Iterators.hpp
index 2a840ae..4f7fc20 100644
--- a/src/Iterators.hpp
+++ b/src/Iterators.hpp
@@ -192,6 +192,8 @@ class BoundedSumIterator {
 
   size_t dim() const { return d; }
 
+  size_t firstIndexBound() const { return bound + 1; }
+
   std::vector<size_t> indexBounds() const { return std::vector<size_t>(d, bound + 1); }
 
   size_t numValues() const { return binom(bound + d, d); }
diff --git a/test/main.cpp b/test/main.cpp
index 512ab68..5e4b4ae 100644
--- a/test/main.cpp
+++ b/test/main.cpp
@@ -26,8 +26,8 @@ limitations under the License.
  */
 
 void runFunctions() {
-  constexpr size_t d = 3;
-  size_t bound = 100;
+  constexpr size_t d = 30;
+  size_t bound = 8;
   fsi::TemplateBoundedSumIterator<d> it(bound);
   // fsi::BoundedSumIterator it(d, bound);
   std::vector<MonomialFunctions> phi(d);