added documentation/comments

README.md

@@ -1,10 +1,25 @@
-This is a header-only library providing a fast matrix-vector product and linear system solver for tensor product matrices with a downward-closed index set restriction. This can especially be applied to compute interpolation coefficients for a sparse grid basis or evaluate a sparse interpolant at sparse grid points.
+This is a header-only library providing a fast matrix-vector product and linear system solver for tensor product matrices with a downward-closed index set restriction. This can especially be applied to compute interpolation coefficients for a sparse grid basis or evaluate a sparse interpolant at sparse grid points. The library code is located in the src folder.
 
-The fast_sparse_interpolation library is published under an Apache 2.0 license. If you use this project for research purposes, please cite one of the following publications:
+The fast_sparse_interpolation library is published under an Apache 2.0 license. If you use this project for research purposes, please cite the following publication, which describes the mathematical background:
 - David Holzmüller, Dirk Pflüger: Fast Sparse Grid Interpolation Coefficients Via LU Decomposition (2019).
 
-Example code for usage can be found in test/main.cpp. This file can be compiled by executing the top-level Makefile (make debug / make release). This Makefile executes cmake on CMakeLists.txt to generate a Makefile in the build folder, which is then automatically executed using make to compile the library. 
+The algorithm was first proposed in:
+- Gustavo Avila and Tucker Carrington Jr.: A multi-dimensional Smolyak collocation method in curvilinear coordinates for computing vibrational spectra (2015).
 
+Example code for usage can be found in test/main.cpp. This file can be compiled by executing the top-level Makefile (make debug / make release). This Makefile executes cmake on CMakeLists.txt to generate a Makefile in the build folder, which is then automatically executed using make to compile the executable test code. The test folder also contains code for performance measurement and plotting of performance data, which we used for our 2019 paper.
+
+All functions and classes lie inside the namespace fsi (short for fast_sparse_interpolation). The basic structure is as follows:
+- The class MultiDimVector<IteratorType> stores vectors indexed by multi-indices.
+- The class SparseTPOperator<IteratorType> implicitly stores matrices of a certain type (restrictions of tensor product matrices to a downward closed index set), as described in the paper. It implements a matrix-vector product (apply()) and an inverse-matrix-vector product (solve()).
+- The class BoundedSumIterator provides an iterator that allows to iterate over a downward closed index set, specifically a set of the form {(i_1, ..., i_d) | i_1, ..., i_d >= 0, i_1 + ... + i_d <= bound}. It is possible to provide other iterator classes with the same interface in order to use different downward closed index sets.
+- The function evaluateFunction() evaluates a function at grid points and returns a MultiDimVector with the function values.
+- The function createInterpolationOperator() creates a SparseTPOperator whose matrix corresponds to a sparse grid interpolation problem with given basis functions and grid points.
+- For custom extensions (such as evaluating derivatives), helper functions for multiplication with matrices of the form ($I \otimes \hdots \otimes I \otimes M$) are available (for arbitrary M, for lower triangular M, for upper triangular M and for M = I). Each single multiplication permutes the indices in the MultiDimVector, after d multiplications (where d is the number of dimensions), the indices are ordered as in the beginning.
+
+The basis functions and grid points used in the examples are only easy-to-implement examples. For practical purposes, you should probably use other basis functions and grid points, for example ones provided in the SG++ library. Recommendations:
+- Legendre polynomials and Leja points (or L2-Leja points, see SG++)
+- Splines and uniform points
+- Maybe also kernel functions and uniform points
 
 Note for developers:
 To run the script ./prepend_header.sh, modify it to only prepend the header in files that do not contain one and run

src/Interpolation.hpp

@@ -23,32 +23,21 @@ limitations under the License.
 
 namespace fsi {
 
-inline std::ostream &operator<<(std::ostream &os, boost::numeric::ublas::matrix<double> matrix) {
+/**
-  os << "[";
+ * Represents a vector indexed by multi-indices which are specified by an iterator of type It.
-  for (size_t i = 0; i < matrix.size1(); ++i) {
+ */
-    for (size_t j = 0; j < matrix.size2(); ++j) {
-      os << matrix(i, j) << "   ";
-    }
-    os << "\n";
-  }
-  os << " ]";
-  return os;
-}
-
-template <class T>
-inline std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
-  os << "[";
-  for (auto ii = v.begin(); ii != v.end(); ++ii) {
-    os << " " << *ii;
-  }
-  os << " ]";
-  return os;
-}
-
 template <typename It>
 class MultiDimVector {
   It it;
 
+ public:
+  // put the first dimension into an outer vector for processing reasons
+  // data[i] contains all entries whose multi-index starts with i, ordered lexicographically.
+  std::vector<std::vector<double>> data;
+
+  /**
+   * Iterator class used for begin() / end()
+   */
   class Iterator {
     MultiDimVector<It> &v;
     StepIterator<It> stepIt;
@@ -90,10 +79,6 @@ class MultiDimVector {
     bool operator!=(const self_type &rhs) { return !(*this == rhs); }
   };
 
- public:
-  // put the first dimension into an outer vector for processing reasons
-  std::vector<std::vector<double>> data;
-
   MultiDimVector(It it) : it(it), data(it.firstIndexBound()) {
     auto sizes = it.numValuesPerFirstIndex();
     for (size_t i = 0; i < sizes.size(); ++i) {
@@ -106,14 +91,15 @@ class MultiDimVector {
     std::swap(it, other.it);
   }
 
-  void clear() {
+  /**
-    for (size_t dim = 0; dim < data.size(); ++dim) {
+   * Returns the associated iterator object.
-      data[dim].clear();
+   */
-    }
-  }
-
   It getJumpIterator() const { return it; }
 
+  /**
+   * Changes the associated iterator object and adjusts the structure of the data vectors (without
+   * initializing them).
+   */
   void resetWithJumpIterator(It const &other_it) {
     it = other_it;
     size_t n_1 = it.firstIndexBound();
@@ -151,6 +137,34 @@ class MultiDimVector {
   }
 };
 
+/**
+ * Helper for printing a matrix
+ */
+inline std::ostream &operator<<(std::ostream &os, boost::numeric::ublas::matrix<double> matrix) {
+  os << "[";
+  for (size_t i = 0; i < matrix.size1(); ++i) {
+    for (size_t j = 0; j < matrix.size2(); ++j) {
+      os << matrix(i, j) << "   ";
+    }
+    os << "\n";
+  }
+  os << " ]";
+  return os;
+}
+
+/**
+ * Helper for printing a std::vector
+ */
+template <class T>
+inline std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
+  os << "[";
+  for (auto ii = v.begin(); ii != v.end(); ++ii) {
+    os << " " << *ii;
+  }
+  os << " ]";
+  return os;
+}
+
 template <typename It>
 MultiDimVector<It> operator+(MultiDimVector<It> first, MultiDimVector<It> const &second) {
   first += second;
@@ -175,6 +189,21 @@ double squared_l2_norm(MultiDimVector<It> const &v) {
   return sum;
 }
 
+// ----- The next block of methods contains methods for efficient multi-indexed matrix-vector
+// products.
+
+/**
+ * Multiplies the multi-indexed vector v with the matrix \widehat{I \otimes \hdots \otimes I \otimes
+ * L}, according to the notation from the paper. Saves the result in v. The vector buffer is used
+ * for storing intermediate results. It should be provided externally and can be used for multiple
+ * matrix-vector products to omit allocating new memory each time.
+ * Note that this method permutes indices, specifically it moves the last index to the front. This
+ * means e.g. that an element at index (0, 1, 2, 3) of the resulting vector would have been at index
+ * (1, 2, 3, 0) if the indices hadn't been permuted.
+ * We cycle indices because it allows for almost continuous memory access during each
+ * multiplication, hence likely accelerating the method due to better cache performance and better
+ * iterator performance.
+ */
 template <typename It>
 void multiply_single_lower_triangular_inplace(boost::numeric::ublas::matrix<double> const &L,
                                               MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
@@ -216,6 +245,9 @@ void multiply_single_lower_triangular_inplace(boost::numeric::ublas::matrix<doub
   v.swap(buffer);
 }
 
+/**
+ * Same as above, but with an upper triangular matrix U instead of L.
+ */
 template <typename It>
 void multiply_single_upper_triangular_inplace(boost::numeric::ublas::matrix<double> const &U,
                                               MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
@@ -257,6 +289,9 @@ void multiply_single_upper_triangular_inplace(boost::numeric::ublas::matrix<doub
   v.swap(buffer);
 }
 
+/**
+ * Same as above, but with an arbitrary Matrix M instead of L.
+ */
 template <typename It>
 void multiply_single_arbitrary_inplace(boost::numeric::ublas::matrix<double> const &M,
                                        MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
@@ -298,6 +333,11 @@ void multiply_single_arbitrary_inplace(boost::numeric::ublas::matrix<double> con
   v.swap(buffer);
 }
 
+/**
+ * Same as above, but multiplies with the identity matrix, i.e. it only cycles indices. This should
+ * be somewhat faster than explicitly multiplying with the identity matrix using one of the other
+ * three methods.
+ */
 template <typename It>
 void multiply_single_identity_inplace(MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
   It it = v.getJumpIterator();
@@ -334,6 +374,9 @@ void multiply_single_identity_inplace(MultiDimVector<It> &v, MultiDimVector<It>
   v.swap(buffer);
 }
 
+/**
+ * Repeats calls to the identity multiplication to cycle the indices several times.
+ */
 template <typename It>
 void cycle_vector_inplace(size_t num_times, MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
   for (size_t i = 0; i < num_times; ++i) {
@@ -341,6 +384,10 @@ void cycle_vector_inplace(size_t num_times, MultiDimVector<It> &v, MultiDimVecto
   }
 }
 
+/**
+ * Repeats calls to the identity multiplication to cycle the indices several times. Uses its own
+ * buffer (potentially slower than reusing an already created buffer).
+ */
 template <typename It>
 void cycle_vector_inplace(size_t num_times, MultiDimVector<It> &v) {
   MultiDimVector<It> buffer(v.getJumpIterator());
@@ -349,6 +396,10 @@ void cycle_vector_inplace(size_t num_times, MultiDimVector<It> &v) {
   }
 }
 
+/**
+ * Performs a multiplication with a tensor product of lower triangular matrices \widehat{L[1]
+ * \otimes \hdots \otimes L[k]}.
+ */
 template <typename It>
 void multiply_lower_triangular_inplace(std::vector<boost::numeric::ublas::matrix<double>> L,
                                        MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
@@ -362,6 +413,10 @@ void multiply_lower_triangular_inplace(std::vector<boost::numeric::ublas::matrix
   }
 }
 
+/**
+ * Performs a multiplication with a tensor product of upper triangular matrices \widehat{U[1]
+ * \otimes \hdots \otimes U[k]}.
+ */
 template <typename It>
 void multiply_upper_triangular_inplace(std::vector<boost::numeric::ublas::matrix<double>> U,
                                        MultiDimVector<It> &v, MultiDimVector<It> &buffer) {
@@ -376,24 +431,46 @@ void multiply_upper_triangular_inplace(std::vector<boost::numeric::ublas::matrix
 }
 
 /**
- * Represents a sparse linear tensor product operator defined by a matrix for each dimension.
+ * Represents a sparse linear tensor product operator defined by a matrix for each dimension and an
+ * iterator that defines the multi-index set.
  */
 template <typename It>
 class SparseTPOperator {
+  // Iterator that defines the multi-index set
   It it;
+
+  // Matrices associated with the operator.
   std::vector<boost::numeric::ublas::matrix<double>> M;
+
+  // Matrices containing the L and U parts of the LU decompositions
   std::vector<boost::numeric::ublas::matrix<double>> LU;
+
+  // Matrices containing the L parts of the LU decompositions
   std::vector<boost::numeric::ublas::matrix<double>> L;
+
+  // Matrices containing the U parts of the LU decomposition
   std::vector<boost::numeric::ublas::matrix<double>> U;
+
+  // Inverses of the L matrices
   std::vector<boost::numeric::ublas::matrix<double>> Linv;
+
+  // Inverses of the U matrices
   std::vector<boost::numeric::ublas::matrix<double>> Uinv;
+
+  // is used to store intermediate data
   MultiDimVector<It> buffer;
+
+  // dimension of the indices
   size_t d;
 
  public:
   SparseTPOperator(It it, std::vector<boost::numeric::ublas::matrix<double>> matrices)
       : it(it), M(matrices), buffer(it), d(matrices.size()){};
 
+  /**
+   * This method performs preparations common to apply() and solve() if they haven't already been
+   * performed. The preparations are the LU decompositions of the M matrices.
+   */
   void prepareCommon() {
     if (LU.size() > 0) {
       return;  // already prepared
@@ -412,6 +489,10 @@ class SparseTPOperator {
     }
   }
 
+  /**
+   * This method performs preparations specific to apply() if they haven't already been performed.
+   * These are the (cheap) computation of L and U from LU.
+   */
   void prepareApply() {
     if (L.size() > 0) {
       return;  // already prepared
@@ -445,6 +526,10 @@ class SparseTPOperator {
     }
   }
 
+  /**
+   * This method performs preparations specific to solve() if they haven't already been performed.
+   * This is the inversion of the L and U matrices.
+   */
   void prepareSolve() {
     if (Linv.size() > 0) {
       return;  // already prepared
@@ -467,6 +552,9 @@ class SparseTPOperator {
     }
   }
 
+  /**
+   * Performs a matrix-vector product with the matrix that is implicitly defined by this operator.
+   */
   MultiDimVector<It> apply(MultiDimVector<It> input) {
     prepareApply();
     multiply_upper_triangular_inplace(U, input, buffer);
@@ -474,6 +562,10 @@ class SparseTPOperator {
     return input;
   }
 
+  /**
+   * Performs a matrix-vector product with the inverse of the matrix that is implicitly defined by
+   * this operator.
+   */
   MultiDimVector<It> solve(MultiDimVector<It> rhs) {
     prepareSolve();
     multiply_lower_triangular_inplace(Linv, rhs, buffer);
@@ -482,6 +574,11 @@ class SparseTPOperator {
   }
 };
 
+/**
+ * Creates a SparseTPOperator corresponding to the interpolation problem. That means its
+ * one-dimensional matrices contain the entries phi[k](j)(x[k](i)). The type of phi and x must be
+ * such that this expression is defined.
+ */
 template <typename It, typename X, typename Phi>
 SparseTPOperator<It> createInterpolationOperator(It it, Phi phi, X x) {
   auto n = it.indexBounds();
@@ -508,6 +605,12 @@ SparseTPOperator<It> createInterpolationOperator(It it, Phi phi, X x) {
   return SparseTPOperator<It>(it, matrices);
 }
 
+/**
+ * Creates a MultiDimVector containing the values of the function f at the grid points specified by
+ * x and the multi-index iterator it. The function f should take a std::vector<double> as an
+ * argument. The point object x should have a type such that the expression x[k](i) yields a double
+ * specifying the i-th point (starting from i=0) in dimension k (also starting from k=0).
+ */
 template <typename It, typename Func, typename X>
 MultiDimVector<It> evaluateFunction(It it, Func f, X x) {
   size_t d = it.dim();
@@ -538,6 +641,11 @@ MultiDimVector<It> evaluateFunction(It it, Func f, X x) {
   return v;
 }
 
+/**
+ * Convenience function that computes a vector of interpolation coefficients for the sparse grid
+ * basis given by phi and it on the points given by phi and it. The arguments f, it, phi, x need to
+ * satisfy the same requirements as explained in the documentation of the functions above.
+ */
 template <typename Func, typename It, typename Phi, typename X>
 MultiDimVector<It> interpolate(Func f, It it, Phi phi, X x) {
   auto rhs = evaluateFunction(it, f, x);

src/Iterators.hpp

@@ -9,6 +9,9 @@
 
 namespace fsi {
 
+/**
+ * Computes a binomial coefficient.
+ */
 inline size_t binom(size_t n, size_t k) {
   if (2 * k > n) {
     k = n - k;
@@ -26,74 +29,55 @@ inline size_t binom(size_t n, size_t k) {
   return prod;
 }
 
-template <typename JumpIt>
+/**
-class StepIterator {
+ * An iterator class for iterating lexicographically over a multi-index set of the form
-  JumpIt jump_it;
+ * {(i_1, ..., i_d) | i_1, ..., i_d >= 0, i_1 + ... + i_d <= bound}. Other iterators for other index
-  size_t first_dim_value;
+ * sets may be implemented following the interface of this iterator class.
-  size_t last_dim_value;
+ *
-  size_t last_dim_count;
+ * Since the fast matrix-vector product sums over the last dimension, this iterator does not iterate
-  size_t tail_dims_counter;
+ * over the last dimension. For example, if d = 3 and bound = 2, the index set is (lexicographically
-
+ * ordered)
- public:
+ * {(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 1, 0), (0, 1, 1), (0, 2, 0), (1, 0, 0), (1, 0, 1), (1, 1,
-  StepIterator(JumpIt jump_it)
+ * 0), (2, 0, 0)}. The iterator iterates over the indices
-      : jump_it(jump_it),
+ * {(0, 0, 0), (0, 1, 0), (0, 2, 0), (1, 0, 0), (1, 1, 0), (2, 0, 0)} and for each index provides a
-        first_dim_value(0),
+ * number n >= 1 that indicates how many values the last dimension can take. For example, when the
-        last_dim_value(0),
+ * iterator is at the index (0, 1, 0), we would have n = 2 since (0, 1, 0) and (0, 1, 1) belong to
-        last_dim_count(jump_it.lastDimensionCount()),
+ * the index set. In contrast, we would have n = 3 for (0, 0, 0) and n = 1 for (2, 0, 0).
-        tail_dims_counter(0){};
+ */
-
-  void next() {
-    tail_dims_counter += 1;
-    last_dim_value += 1;
-    if (last_dim_value >= last_dim_count) {
-      last_dim_value = 0;
-      jump_it.next();
-      if (jump_it.firstIndex() != first_dim_value) {
-        first_dim_value = jump_it.firstIndex();
-        tail_dims_counter = 0;
-      }
-      last_dim_count = jump_it.lastDimensionCount();
-    }
-  }
-
-  bool valid() const { return jump_it.valid(); }
-
-  void reset() {
-    jump_it.reset();
-    first_dim_value = 0;
-    last_dim_value = 0;
-    last_dim_count = jump_it.lastDimensionCount();
-    tail_dims_counter = 0;
-  }
-
-  size_t firstIndex() const { return first_dim_value; }
-  size_t tailDimsCounter() const { return tail_dims_counter; }
-
-  std::vector<size_t> index() const {
-    std::vector<size_t> result = jump_it.getIndexHead();
-    result.push_back(last_dim_value);
-    return result;
-  }
-};
-
 class BoundedSumIterator {
   size_t d;
   size_t bound;
-  std::vector<size_t> index_head;  // contains all entries of the index except the last one
+
+  // contains all entries of the index except the last one
+  std::vector<size_t> index_head;
+
+  // contains the sum of the entries of index_head
   size_t index_head_sum;
+
+  // stores whether the iterator is still at a valid position.
+  // This is set to false whenever the iterator is in the state (bound, 0, ..., 0) and next() is
+  // called.
   bool is_valid;
 
  public:
+  /**
+   * Constructor with dimension d of the indices (d includes the dimension that is not iterated
+   * over) and the bound for the sum of indices. Initially, the iterator points to the multi-index
+   * (0, ..., 0).
+   */
   BoundedSumIterator(size_t d, size_t bound)
       : d(d), bound(bound), index_head(d - 1, 0), index_head_sum(0), is_valid(true){};
 
   /**
    * At the current multi-index (i_1, ..., i_{d-1}, 0), return how many multi-indices starting with
    * (i_1, ..., i_{d-1}) are contained in the multi-index set, then advance to the next multi-index
-   * that ends with a zero.
+   * that ends with a zero. This corresponds to the value n in the class description.
    */
   size_t lastDimensionCount() { return bound - index_head_sum + 1; }
 
+  /**
+   * Advances to the next multi-index ending with zero.
+   */
   void next() {
     if (bound > index_head_sum) {
       index_head_sum += 1;
@@ -117,28 +101,62 @@ class BoundedSumIterator {
     }
   }
 
+  /**
+   * Returns the component i_1 of the current multi-index (i_1, ..., i_{d-1}, 0).
+   */
   size_t firstIndex() const { return index_head[0]; }
 
+  /**
+   * Returns the component of the multi-index at dimension dim, where dim < d-1.
+   */
   size_t indexAt(size_t dim) const { return index_head[dim]; }
 
+  /**
+   * Returns the vector (i_1, ..., i_{d-1}) where (i_1, \hdots, i_{d-1}, 0) is the current
+   * multi-index.
+   */
   std::vector<size_t> getIndexHead() const { return index_head; }
 
+  /**
+   * Returns true if the iterator points to a valid multi-index (i.e. the iteration is not
+   * finished).
+   */
   bool valid() const { return is_valid; }
 
+  /**
+   * Resets the iterator to (0, ..., 0) and resets valid() to true.
+   */
   void reset() {
     index_head = std::vector<size_t>(d - 1, 0);
     index_head_sum = 0;
     is_valid = true;
   }
 
+  /**
+   * Returns the dimension of the multi-indices (including the last-dimension index).
+   */
   size_t dim() const { return d; }
 
+  /**
+   * Returns the maximum value that the first dimension of the index can take, plus one.
+   */
   size_t firstIndexBound() const { return bound + 1; }
 
+  /**
+   * Returns the maximum values that the dimensions can take, plus one.
+   */
   std::vector<size_t> indexBounds() const { return std::vector<size_t>(d, bound + 1); }
 
+  /**
+   * Returns the total number of multi-indices that are valid (including those with last dimension
+   * != 0).
+   */
   size_t numValues() const { return binom(bound + d, d); }
 
+  /**
+   * For any value 0 <= i_1, i_1 < firstIndexBound(), returns how many valid multi-indices there are
+   * that start with i_1.
+   */
   std::vector<size_t> numValuesPerFirstIndex() const {
     std::vector<size_t> result;
     for (size_t firstIndex = 0; firstIndex <= bound; ++firstIndex) {
@@ -147,6 +165,10 @@ class BoundedSumIterator {
     return result;
   }
 
+  /**
+   * Moves the iterator directly to the position that is reached after the iteration is finished
+   * (with valid() == false). This is helpful for implementing MultiDimVector.end().
+   */
   void goToEnd() {
     index_head = std::vector<size_t>(d - 1, 0);
     index_head_sum = 0;
@@ -163,4 +185,83 @@ class BoundedSumIterator {
     return it;
   }
 };
+
+/**
+ * Takes an iterator that does not iterate over the last index (such as a BoundedSumIterator) and
+ * provides an interface to also iterate over the last index. This is used internally in
+ * MultiDimVector, but can also used for other purposes.
+ */
+template <typename JumpIt>
+class StepIterator {
+  JumpIt jump_it;
+  size_t first_dim_value;
+  size_t last_dim_value;
+  size_t last_dim_count;
+  size_t tail_dims_counter;
+
+ public:
+  /**
+   * Initializes to the multi-index (0, ..., 0).
+   * JumpIt is the iterator object that does not iterate over the last index.
+   */
+  StepIterator(JumpIt jump_it)
+      : jump_it(jump_it),
+        first_dim_value(0),
+        last_dim_value(0),
+        last_dim_count(jump_it.lastDimensionCount()),
+        tail_dims_counter(0){};
+
+  /**
+   * Moves to the next multi-index in lexicographical order.
+   */
+  void next() {
+    tail_dims_counter += 1;
+    last_dim_value += 1;
+    if (last_dim_value >= last_dim_count) {
+      last_dim_value = 0;
+      jump_it.next();
+      if (jump_it.firstIndex() != first_dim_value) {
+        first_dim_value = jump_it.firstIndex();
+        tail_dims_counter = 0;
+      }
+      last_dim_count = jump_it.lastDimensionCount();
+    }
+  }
+
+  /**
+   * Returns true if the iterator points to a valid multi-index.
+   */
+  bool valid() const { return jump_it.valid(); }
+
+  /**
+   * Resets the iterator to its initial state.
+   */
+  void reset() {
+    jump_it.reset();
+    first_dim_value = 0;
+    last_dim_value = 0;
+    last_dim_count = jump_it.lastDimensionCount();
+    tail_dims_counter = 0;
+  }
+
+  /**
+   * Returns the value i_1 of the current multi-index (i_1, ..., i_d).
+   */
+  size_t firstIndex() const { return first_dim_value; }
+
+  /**
+   * Returns n, such that the current multi-index is the (n+1)-th multi-index that starts with
+   * firstIndex().
+   */
+  size_t tailDimsCounter() const { return tail_dims_counter; }
+
+  /**
+   * Returns a std::vector<size_t> containing the multi-index.
+   */
+  std::vector<size_t> index() const {
+    std::vector<size_t> result = jump_it.getIndexHead();
+    result.push_back(last_dim_value);
+    return result;
+  }
+};
 }

test/main.cpp

@@ -15,18 +15,104 @@ limitations under the License.
 
 #include "common.hpp"
 
-// TODO: refactoring:
 /**
- * - Tests?
+ * Possible additions:
- * - More point distributions / Basis functions?
+ * - Tests
- * - Computation of derivatives?
+ * - More point distributions / Basis functions
- * - Examples?
+ * - Computation of derivatives
  */
 
+/**
+ * Simple interpolation code. The function to evaluate is f(x, y, z) = x*y*z, so it can be
+ * perfectly interpolated using the tensor product of the linear monomials at index (1, 1, 1).
+ */
+void example1() {
+  // Choose a dimensionality of the problem.
+  size_t d = 3;
+
+  // As multi-indices indexing our points and basis functions, this implementation provides an
+  // iterator for the standard index set
+  // {(i_1, ..., i_d) | i_1, ..., i_d >= 0, i_1 + ... + i_d <= bound}.
+  // We configure the bound here.
+  size_t bound = 3;
+
+  // This is the iterator that (implicitly) specifies the multi-index set used for the grid points
+  // and basis functions.
+  fsi::BoundedSumIterator it(d, bound);
+
+  // This is a vector of function-like objects such that we can write phi[k](i)(x) for the i-th
+  // basis function (size_t i >= 0) in dimension k, evaluated at the point x (double).
+  // As an example, we provide Monomials (phi[k](i)(x) = x^i).
+  // Note that the Vandermonde matrix (the interpolation matrix) for monomials is badly conditioned.
+  // Therefore, the result gets quite bad when we use bounds >= 30.
+  // For this purpose, it would be better to use Legendre polynomials or similar things.
+  std::vector<MonomialFunctions> phi(d);
+
+  // This is a vector of function-like objects such that we can write x[k](i) for the i-th point
+  // (size_t i >= 0) in dimension k. This distribution is only an example, in practice things like
+  // uniform distributions for
+  std::vector<GoldenPointDistribution> x(d);
+
+  // Determine a (multi-dimensional) vector with the interpolation coefficients.
+  auto coefficients = fsi::interpolate(f, it, phi, x);
+
+  // Print the coefficients and their multi-indices. These multi-indices are those that the iterator
+  // it implicitly specifies.
+  for (auto coeff_it = coefficients.begin(); coeff_it != coefficients.end(); ++coeff_it) {
+    std::cout << "Coefficient at index " << coeff_it.index() << ": " << *coeff_it << "\n";
+  }
+}
+
+/**
+ * More detailed interpolation and evaluation code.
+ */
+void example2() {
+  size_t d = 8;
+  size_t bound = 24;
+  fsi::BoundedSumIterator it(d, bound);
+  std::vector<MonomialFunctions> phi(d);
+  std::vector<GoldenPointDistribution> x(d);
+
+  // Create a vector containing the function values at the grid points
+  auto rhs = fsi::evaluateFunction(it, f, x);
+
+  // Create an object representing the linear sparse tensor product matrix corresponding to the
+  // interpolation problem
+  auto op = fsi::createInterpolationOperator(it, phi, x);
+
+  // Run the LU decomposition and invert L and U. If this function is not called, it is
+  // automatically called once you call solve(). It can be called separately for performance
+  // measurements. This preparation will only done once for an operator object, no matter how often
+  // solve() or prepareSolve() is called.
+  op.prepareSolve();
+
+  // Solve the linear system with the vector of function values as a right hand side to get a vector
+  // of coefficients.
+  auto coefficients = op.solve(rhs);
+
+  // We can also prepare the apply operation, which will do almost nothing here since the LU
+  // decomposition has already been done. As op.prepareSolve(), this is not necessary.
+  op.prepareApply();
+
+  // The apply operation performs the matrix-vector product of the operator's sparse tensor product
+  // matrix with the coefficients vector. In our case, this means that we compute the values of the
+  // interpolant at the grid points based on its coefficients.
+  auto b = op.apply(coefficients);
+
+  // Verify that solve and apply are numerically stable, i.e. that b and rhs are (almost) equal.
+  std::cout << "Reconstruction error (L2 norm): " << sqrt(fsi::squared_l2_norm(b - rhs)) << "\n";
+
+  // Print the number of grid points (which equals the number of basis functions) that are used in
+  // total.
+  std::cout << "Number of points: " << it.numValues() << "\n";
+}
+
+/**
+ * Similar to example 2, but includes timing information.
+ */
 void runFunctions() {
-  constexpr size_t d = 8;
+  size_t d = 8;
   size_t bound = 30;
-  // fsi::TemplateBoundedSumIterator<d> it(bound);
   fsi::BoundedSumIterator it(d, bound);
   std::vector<MonomialFunctions> phi(d);
   std::vector<GoldenPointDistribution> x(d);
@@ -39,12 +125,10 @@ void runFunctions() {
   Timer timer;
   timer.reset();
   op.prepareSolve();
-  // auto result = fsi::interpolate(f, it, phi, x);
   std::cout << "Time for prepareSolve(): " << timer.elapsed() << " s\n";
 
   timer.reset();
   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";
 
@@ -73,5 +157,7 @@ void runFunctions() {
 // using namespace fsi;
 int main() {
   // runFunctions();
-  measurePerformance();
+  // measurePerformance();
+  // example1();
+  example2();
 }