moved vector index counting back because it contained an error and was not faster and more ugly

src/Interpolation.hpp

@@ -77,11 +77,15 @@ void multiply_lower_triangular_inplace(It it, std::vector<boost::numeric::ublas:
     auto &Lk = L[k];
     it.reset();
 
+    size_t first_v_index = 0;
+    size_t second_v_index = 0;
+
     double *data_pointer = &v.data[0][0];
+    size_t data_size = v.data[0].size();
 
     while (not it.done()) {
       size_t last_dim_count = it.lastDimensionCount();
-      double *offset_data_pointer = data_pointer + it.getMiddleDimensionsCounter();
+      double *offset_data_pointer = data_pointer + second_v_index;
       for (size_t i = 0; i < last_dim_count; ++i) {
         double sum = 0.0;
         for (size_t j = 0; j <= i; ++j) {
@@ -89,11 +93,15 @@ void multiply_lower_triangular_inplace(It it, std::vector<boost::numeric::ublas:
         }
         w.data[i].push_back(sum);
       }
-
+      second_v_index += last_dim_count;
-      bool first_dimension_changed = it.next();
+      if (second_v_index >= data_size) {
-      if (first_dimension_changed) {
+        second_v_index = 0;
-        data_pointer = &v.data[it.firstIndex()][0];
+        first_v_index += 1;
+        data_pointer = &v.data[first_v_index][0];
+        data_size = v.data[first_v_index].size();
       }
+
+      it.next();
     }
 
     v.swap(w);
@@ -124,11 +132,15 @@ void multiply_upper_triangular_inplace(It it, std::vector<boost::numeric::ublas:
     auto &Uk = U[k];
     it.reset();
 
+    size_t first_v_index = 0;
+    size_t second_v_index = 0;
+
     double *data_pointer = &v.data[0][0];
+    size_t data_size = v.data[0].size();
 
     while (not it.done()) {
       size_t last_dim_count = it.lastDimensionCount();
-      double *offset_data_pointer = data_pointer + it.getMiddleDimensionsCounter();
+      double *offset_data_pointer = data_pointer + second_v_index;
       for (size_t i = 0; i < last_dim_count; ++i) {
         double sum = 0.0;
         for (size_t j = i; j < last_dim_count; ++j) {
@@ -136,11 +148,15 @@ void multiply_upper_triangular_inplace(It it, std::vector<boost::numeric::ublas:
         }
         w.data[i].push_back(sum);
       }
-
+      second_v_index += last_dim_count;
-      bool first_dimension_changed = it.next();
+      if (second_v_index >= data_size) {
-      if (first_dimension_changed) {
+        second_v_index = 0;
-        data_pointer = &v.data[it.firstIndex()][0];
+        first_v_index += 1;
+        data_pointer = &v.data[first_v_index][0];
+        data_size = v.data[first_v_index].size();
       }
+
+      it.next();
     }
 
     it = it.cycle();

src/Iterators.hpp

@@ -31,16 +31,11 @@ class TemplateBoundedSumIterator {
   size_t bound;
   std::vector<size_t> index_head;  // contains all entries of the index except the last one
   size_t index_head_sum;
-  size_t middle_dimensions_counter;
   bool is_done;
 
  public:
   TemplateBoundedSumIterator(size_t bound)
-      : bound(bound),
+      : bound(bound), index_head(d - 1, 0), index_head_sum(0), is_done(false){};
-        index_head(d - 1, 0),
-        index_head_sum(0),
-        middle_dimensions_counter(0),
-        is_done(false){};
 
   /**
    * At the current multi-index (i_1, ..., i_{d-1}, 0), return how many multi-indices starting with
@@ -49,9 +44,7 @@ class TemplateBoundedSumIterator {
    */
   size_t lastDimensionCount() { return bound - index_head_sum + 1; }
 
-  bool next() {
+  void next() {
-    middle_dimensions_counter += lastDimensionCount();
-
     if (bound > index_head_sum) {
       index_head_sum += 1;
       index_head[d - 2] += 1;
@@ -62,23 +55,16 @@ class TemplateBoundedSumIterator {
         // reduce dimension until entry is nonzero
       }
 
-      if (dim > 1) {
+      if (dim > 0) {
         index_head_sum -= (index_head[dim] - 1);
         index_head[dim] = 0;
         index_head[dim - 1] += 1;
-      } else if (dim == 1) {
-        index_head_sum -= (index_head[1] - 1);
-        index_head[1] = 0;
-        index_head[0] += 1;
-        middle_dimensions_counter = 0;
-        return true;
       } else if (dim == 0) {
         index_head[dim] = 0;
         index_head_sum = 0;
         is_done = true;
       }
     }
-    return false;
   }
 
   size_t firstIndex() const { return index_head[0]; }
@@ -87,12 +73,9 @@ class TemplateBoundedSumIterator {
 
   bool done() const { return is_done; }
 
-  size_t getMiddleDimensionsCounter() const { return middle_dimensions_counter; }
-
   void reset() {
     index_head = std::vector<size_t>(d - 1, 0);
     index_head_sum = 0;
-    middle_dimensions_counter = 0;
     is_done = false;
   }
 

test/main.cpp

@@ -112,8 +112,8 @@ class Timer {
 
 // using namespace fsi;
 int main() {
-  constexpr size_t d = 8;
+  constexpr size_t d = 2;
-  size_t bound = 24;
+  size_t bound = 2;
   fsi::TemplateBoundedSumIterator<d> it(bound);
   // fsi::BoundedSumIterator it(d, bound);
   std::vector<MonomialFunctions> phi(d);
@@ -122,7 +122,7 @@ int main() {
   auto rhs = evaluateFunction(it, f, x);
   auto op = createInterpolationOperator(it, phi, x);
 
-  // std::cout << rhs << "\n";
+  std::cout << rhs << "\n";
 
   Timer timer;
   op.prepareSolve();
@@ -133,12 +133,12 @@ int main() {
   auto c = op.solve(rhs);
   // auto c = fsi::interpolate(f, it, phi, x);
   std::cout << "Time for solve(): " << timer.elapsed() << " s\n";
-  // std::cout << c << "\n";
+  std::cout << c << "\n";
 
   timer.reset();
   auto b = op.apply(c);
   std::cout << "Time for apply(): " << timer.elapsed() << " s\n";
-  // std::cout << b << "\n";
+  std::cout << b << "\n";
 
   std::cout << "Number of points: " << it.numValues() << "\n";