results

lab07/results/aaron_benchmark_failed/jacobiWaveSplitWork.cpp

@@ -0,0 +1,217 @@
+/**
+ *
+ * Schritt 1
+ * -------------
+ * |T0|--|--|--|
+ * -------------
+ * |--|--|--|--|
+ * -------------
+ * |--|--|--|--|
+ * -------------
+ * |--|--|--|--|
+ * -------------
+ *
+ * Schritt 2
+ * -------------
+ * |xx|T0|--|--|
+ * -------------
+ * |T1|--|--|--|
+ * -------------
+ * |--|--|--|--|
+ * -------------
+ * |--|--|--|--|
+ * -------------
+ *
+ * Schritt 3
+ * -------------
+ * |xx|xx|T0|--|
+ * -------------
+ * |xx|T1|--|--|
+ * -------------
+ * |T2|--|--|--|
+ * -------------
+ * |--|--|--|--|
+ * -------------
+ *
+ * Schritt 4
+ * -------------
+ * |xx|xx|xx|T0|
+ * -------------
+ * |xx|xx|T1|--|
+ * -------------
+ * |xx|T2|--|--|
+ * -------------
+ * |T3|--|--|--|
+ * -------------
+ */
+
+#include <benchmark/benchmark.h>
+#include "matrix.h"
+#include <omp.h>
+#include <atomic>
+#include <vector>
+
+void gauss_seidel(Matrix &phi, int maxNumIter)
+{
+  const int m = phi.dim1();
+  const int n = phi.dim2();
+  const double osth = 1. / 4;
+  for (int iIter = 0; iIter < maxNumIter; ++iIter)
+  {
+    // Shared data per iteration
+    std::vector<std::atomic<int>> column(m);
+    for (int i = 0; i < m; ++i)
+    {
+      column[i].store(1);
+    }
+    std::atomic<int> threadsCount(0);
+
+#pragma omp parallel for schedule(static)
+    for (int rowToCalculate = 1; rowToCalculate < (n - 1); rowToCalculate++)
+    {
+      int row = rowToCalculate;
+      while (column[row] < n - 1)
+      {
+        // All threads beneath T_row1 have to check specific circumstances, where
+        // T_row1 has no condition to wait for
+        if (row != 1)
+        {
+          // T_rowx has to wait at least until T_row(x-1) has calculated values from
+          // its last row
+          // T_rowx will calculate value x at row,column if T_row(x-1) has calculated
+          // its value at row-1, column-1
+          while (column[row] == column[row - 1])
+          {
+            // Wait for wave (thread) above
+          }
+        }
+
+        // Central jacobi calculation
+        phi(row, column[row]) = osth * (phi(row + 1, column[row]) + phi(row - 1, column[row]) + phi(row, column[row] + 1) +
+                                        phi(row, column[row] - 1));
+
+        // Increment column index
+        column[row]++;
+      }
+    }
+  }
+}
+
+void gauss_seidel_lin(Matrix &phi, int maxNumIter)
+{
+  const int m = phi.dim1();
+  const int n = phi.dim2();
+
+  const double osth = 1. / 4;
+
+  for (int iter = 0; iter < maxNumIter; ++iter)
+  {
+    for (int i = 1; i < m - 1; ++i)
+    {
+      for (int j = 1; j < n - 1; ++j)
+      {
+        phi(i, j) = osth * (phi(i + 1, j) + phi(i - 1, j) + phi(i, j + 1) +
+                            phi(i, j - 1));
+      }
+    }
+  }
+}
+
+void fill_matrix(Matrix &matrix, const int filler)
+{
+  const int m = matrix.dim1();
+  const int n = matrix.dim2();
+
+  for (int i = 0; i < m; ++i)
+  {
+    for (int j = 0; j < n; ++j)
+    {
+      if (i == 0 || i == m - 1 || j == 0 || j == n - 1)
+      {
+        matrix(i, j) = filler;
+      }
+      else
+        matrix(i, j) = 0;
+    }
+  }
+}
+
+void check(const Matrix &a, const Matrix &b)
+{
+  const int m = a.dim1();
+  const int n = a.dim2();
+
+  for (int i = 0; i < m; ++i)
+  {
+    for (int j = 0; j < n; ++j)
+    {
+      if (a(i, j) != b(i, j))
+      {
+        std::cout << "Not equal at (" << i << ", " << j << "), a: " << a(i, j)
+                  << " != " << b(i, j) << " :b" << std::endl;
+      }
+    }
+  }
+}
+
+void print_matrix(const Matrix &matrix)
+{
+  const int m = matrix.dim1();
+  const int n = matrix.dim2();
+
+  for (int i = 0; i < m; ++i)
+  {
+    for (int j = 0; j < n; ++j)
+    {
+      std::cout << matrix(i, j) << " ";
+    }
+    std::cout << std::endl;
+  }
+}
+
+void benchmarkComputeResult(benchmark::State& state) {
+  int iterations = state.range(0);
+  Matrix sec = Matrix(1000, 1000);
+  fill_matrix(sec, 10);
+
+  for (auto _ : state) {
+    gauss_seidel(sec, iterations);
+    benchmark::DoNotOptimize(sec);
+  }
+}
+
+int main(int argc, char** argv) {
+  ::benchmark::Initialize(&argc, argv);
+
+  for (int iterations = 0; iterations < 10000; iterations++) {
+    benchmark::RegisterBenchmark("idk", benchmarkComputeResult)
+        ->Arg(iterations)
+        ->Unit(benchmark::kMillisecond);
+  }
+
+  ::benchmark::RunSpecifiedBenchmarks();
+
+  return 0;
+}
+
+// int main()
+// {
+//   Matrix first = Matrix(8, 8);
+//   Matrix sec = Matrix(8, 8);
+
+//   fill_matrix(first, 10);
+//   fill_matrix(sec, 10);
+
+//   print_matrix(first);
+//   print_matrix(sec);
+
+//   gauss_seidel_lin(first, 1000);
+//   gauss_seidel(sec, 1000);
+
+//   check(first, sec);
+
+//   print_matrix(first);
+//   print_matrix(sec);
+
+//   return 0;
+// }

lab07/results/aaron_benchmark_failed/jacobiWaveSplitWork.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+#SBATCH --job-name=jacobiWaveSplitWorkTest
+#SBATCH --output=benchmark_jacobiWaveSplitWork.out
+#SBATCH --ntasks=1
+#SBATCH --cpus-per-task=12
+
+# Load modules or activate environment if needed
+# module load gcc cmake ...
+
+# Run the benchmark binary
+./jacobiWaveSplitWork --benchmark_out=results.json --benchmark_out_format=json

lab07/results/aaron_benchmark_failed/matrix.h

@@ -0,0 +1,343 @@
+/**
+ * matrix.h a very simplistic class for m times n matrices.
+ */
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+// A very simplistic vector class for vectors of size n
+class Vector {
+ public:
+  // constructors
+  Vector(int n) : n_(n), data_(n_, 0) {}
+  Vector(const Vector& other) = default;
+  Vector(Vector&& other) = default;
+  ~Vector() = default;
+
+  // assignment operators
+  Vector& operator=(const Vector& other) = default;
+  Vector& operator=(Vector&& other) = default;
+
+  // element access
+  double& operator()(int i) { return data_[i]; }
+  const double& operator()(int i) const { return data_[i]; }
+
+  // getter functions for the dimensions
+  int dim() const { return n_; }
+
+  // comparison operators
+  bool operator==(const Vector& b) { return (data_ == b.data_); }
+  bool operator!=(const Vector& b) { return (data_ != b.data_); }
+
+  // addition
+  Vector& operator+=(const Vector& b) {
+    for (int i = 0; i < n_; ++i) {
+      operator()(i) += b(i);
+    }
+    return *this;
+  }
+
+  // subtraction
+  Vector& operator-=(const Vector& b) {
+    for (int i = 0; i < n_; ++i) {
+      operator()(i) -= b(i);
+    }
+    return *this;
+  }
+
+  // scalar multiplication
+  Vector& operator*=(double x) {
+    for (int i = 0; i < n_; ++i) {
+      operator()(i) *= x;
+    }
+    return *this;
+  }
+
+  // dot product between two vectors
+  double dot(const Vector& other) const {
+    double sum = 0;
+    for (int i = 0; i < n_; ++i) {
+      sum += operator()(i) * other(i);
+    }
+    return sum;
+  }
+
+ private:
+  int n_;                     // vector dimension
+  std::vector<double> data_;  // the vectors entries
+};
+
+inline double dot(const Vector& v1, const Vector& v2) {
+  return v1.dot(v2);
+}
+
+// Print the vector as a table
+inline std::ostream& operator<<(std::ostream& os, const Vector& a) {
+  const int width = 10;
+  const int precision = 4;
+
+  const auto originalPrecision = os.precision();
+  os << std::setprecision(precision);
+
+  for (int i = 0; i < a.dim(); ++i) {
+    os << std::setw(width) << a(i) << " ";
+  }
+
+  os << "\n";
+
+  os << std::setprecision(originalPrecision);
+  return os;
+}
+
+// A very simple class for m times n matrices
+class Matrix {
+ public:
+  // constructors
+  Matrix() : Matrix(0, 0) {}
+  Matrix(int m, int n) : m_(m), n_(n), data_(m_ * n_, 0) {}
+  Matrix(std::pair<int, int> dim) : Matrix(dim.first, dim.second) {}
+  Matrix(int n) : Matrix(n, n) {}
+  Matrix(const Matrix& other) = default;
+  Matrix(Matrix&& other) = default;
+  ~Matrix() = default;
+
+  // assignment operators
+  Matrix& operator=(const Matrix& other) = default;
+  Matrix& operator=(Matrix&& other) = default;
+
+  // element access
+  double& operator()(int i, int j) { return data_[i * n_ + j]; }
+  const double& operator()(int i, int j) const { return data_[i * n_ + j]; }
+
+  // getter functions for the dimensions
+  std::pair<int, int> dim() const { return std::pair<int, int>(m_, n_); }
+  int dim1() const { return m_; }
+  int dim2() const { return n_; }
+  int numEntries() const { return data_.size(); }
+
+  // comparison operators
+  bool operator==(const Matrix& b) { return (data_ == b.data_); }
+  bool operator!=(const Matrix& b) { return (data_ != b.data_); }
+
+  // addition
+  Matrix& operator+=(const Matrix& b) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) += b(i, j);
+      }
+    }
+    return *this;
+  }
+
+  // subtraction
+  Matrix& operator-=(const Matrix& b) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) -= b(i, j);
+      }
+    }
+    return *this;
+  }
+
+  // scalar multiplication
+  Matrix& operator*=(double x) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) *= x;
+      }
+    }
+    return *this;
+  }
+
+  // scalar division
+  Matrix& operator/=(double x) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) /= x;
+      }
+    }
+    return *this;
+  }
+
+  // matrix product (only for square matrices of equal dimension)
+  Matrix& operator*=(const Matrix& b) {
+    if (dim1() != dim2()) {
+      std::cout << "Error in matrix multiplication: no square matrix\n";
+    } else if (dim1() != b.dim1() || dim2() != b.dim2()) {
+      std::cout << "Error in matrix multiplication: dimensions do not match\n";
+    } else {
+      Matrix a = *this;
+      Matrix& c = *this;
+      const int m = dim1();
+      for (int i = 0; i < m; ++i) {
+        for (int j = 0; j < m; ++j) {
+          for (int k = 0; k < m; ++k) {
+            c(i, j) += a(i, k) * b(k, j);
+          }
+        }
+      }
+    }
+
+    return *this;
+  }
+
+ public:
+  int m_;                     // first dimension
+  int n_;                     // second dimension
+  std::vector<double> data_;  // the matrix' entries
+};
+
+// Print the matrix as a table
+inline std::ostream& operator<<(std::ostream& os, const Matrix& a) {
+  const int width = 10;
+  const int precision = 4;
+
+  const auto originalPrecision = os.precision();
+  os << std::setprecision(precision);
+
+  for (int i = 0; i < a.dim1(); ++i) {
+    for (int j = 0; j < a.dim2(); ++j) {
+      os << std::setw(width) << a(i, j) << " ";
+    }
+    if (i != a.dim1() - 1)
+      os << "\n";
+  }
+
+  os << std::setprecision(originalPrecision);
+  return os;
+}
+
+// matrix product
+inline Matrix operator*(const Matrix& a, const Matrix& b) {
+  if (a.dim2() == b.dim1()) {
+    int m = a.dim1();
+    int n = a.dim2();
+    int p = b.dim2();
+    Matrix c(m, p);
+    for (int i = 0; i < m; ++i) {
+      for (int j = 0; j < p; ++j) {
+        for (int k = 0; k < n; ++k) {
+          c(i, j) += a(i, k) * b(k, j);
+        }
+      }
+    }
+    return c;
+  } else {
+    return Matrix(0, 0);
+  }
+}
+
+inline bool equalWithinRange(const Matrix& a,
+                             const Matrix& b,
+                             double eps = 1e-12) {
+  if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
+    return false;
+
+  int m = a.dim1();
+  int n = a.dim2();
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (fabs(a(i, j) - b(i, j)) > eps) {
+        return false;
+      }
+    }
+  }
+
+  return true;
+}
+
+// A very simple class for "3D-Matrices" (tensors) with dimension l x m x n
+class Matrix3D {
+ public:
+  // constructors
+  Matrix3D(int l, int m, int n) : l_(l), m_(m), n_(n), data_(l) {
+    for (int i = 0; i < l_; ++i) {
+      data_[i] = std::vector<std::vector<double>>(m_);
+      for (int j = 0; j < m_; ++j) {
+        data_[i][j] = std::vector<double>(n_, 0);
+      }
+    }
+  }
+  Matrix3D(int n) : Matrix3D(n, n, n) {}
+  Matrix3D(const Matrix3D& other) = default;
+  Matrix3D(Matrix3D&& other) = default;
+  ~Matrix3D() = default;
+
+  // assignment operators
+  Matrix3D& operator=(const Matrix3D& other) = default;
+  Matrix3D& operator=(Matrix3D&& other) = default;
+
+  // element access
+  double& operator()(int i, int j, int k) { return data_[i][j][k]; }
+  const double& operator()(int i, int j, int k) const { return data_[i][j][k]; }
+
+  // getter functions for the dimensions
+  int dim1() const { return l_; }
+  int dim2() const { return m_; }
+  int dim3() const { return n_; }
+
+  // comparison operators
+  bool operator==(const Matrix3D& b) { return (data_ == b.data_); }
+  bool operator!=(const Matrix3D& b) { return (data_ != b.data_); }
+
+  // addition
+  Matrix3D& operator+=(const Matrix3D& b) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) += b(i, j, k);
+        }
+      }
+    }
+    return *this;
+  }
+
+  // substraction
+  Matrix3D& operator-=(const Matrix3D& b) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) -= b(i, j, k);
+        }
+      }
+    }
+    return *this;
+  }
+
+  // scalar multiplication
+  Matrix3D& operator*=(double x) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) *= x;
+        }
+      }
+    }
+    return *this;
+  }
+
+  // scalar division
+  Matrix3D& operator/=(double x) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) /= x;
+        }
+      }
+    }
+    return *this;
+  }
+
+ private:
+  int l_;                                               // first dimension
+  int m_;                                               // second dimension
+  int n_;                                               // third dimension
+  std::vector<std::vector<std::vector<double>>> data_;  // the tensors' entries
+};
+
+#endif  // MATRIX_H

lab07/results/kai_benchmark/kaiJacobiWave.cpp

@@ -0,0 +1,123 @@
+#include "matrix.h"
+#include <benchmark/benchmark.h>
+#include <iostream>
+#include <omp.h>
+
+void gauss_seidel(Matrix &phi, int maxNumIter) {
+  const int m = phi.dim1();
+  const int n = phi.dim2();
+
+  const double osth = 1. / 4;
+
+  for (int iter = 0; iter < maxNumIter; ++iter) {
+    for (int i = 1; i < m - 1; ++i) {
+      for (int j = 1; j < n - 1; ++j) {
+        phi(i, j) = osth * (phi(i + 1, j) + phi(i - 1, j) + phi(i, j + 1) +
+                            phi(i, j - 1));
+      }
+    }
+  }
+}
+
+void gauss_seidel_par(Matrix &phi, int maxNumIter) {
+  const int m = phi.dim1();
+  const int n = phi.dim2();
+
+  const double osth = 1. / 4;
+
+  for (int iter = 0; iter < maxNumIter; ++iter) {
+#pragma omp parallel num_threads(10)
+    {
+      int num_theads = omp_get_num_threads();
+      int chunk = (m - 2) / num_theads;
+      int start = 1 + chunk * omp_get_thread_num();
+      int end = chunk * (omp_get_thread_num() + 1);
+
+      // printf("thread %d, start: %d, end: %d\n", omp_get_thread_num(), start,
+      //        end);
+
+      for (int j = 1; j < n + omp_get_num_threads() - 1; ++j) {
+        for (int i = start; i <= end; ++i) {
+
+          int k = j - omp_get_thread_num();
+
+          // printf("thread %d, i: %d, j: %d, k: %d\n", omp_get_thread_num(), i,
+          // j,
+          //        k);
+
+          if (k > 0 && k < n - 1) {
+            phi(i, k) = osth * (phi(i + 1, k) + phi(i - 1, k) + phi(i, k + 1) +
+                                phi(i, k - 1));
+          }
+#pragma omp barrier
+        }
+      }
+    }
+  }
+}
+
+void fill_matrix(Matrix &matrix, const int filler) {
+  const int m = matrix.dim1();
+  const int n = matrix.dim2();
+
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (i == 0 || i == m - 1 || j == 0 || j == n - 1) {
+        matrix(i, j) = filler;
+      } else
+        matrix(i, j) = 0;
+    }
+  }
+}
+
+void check(const Matrix &a, const Matrix &b) {
+  const int m = a.dim1();
+  const int n = a.dim2();
+
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (a(i, j) != b(i, j)) {
+        std::cout << "Not equal at (" << i << ", " << j << "), a: " << a(i, j)
+                  << " != " << b(i, j) << " :b" << std::endl;
+      }
+    }
+  }
+}
+
+void print_matrix(const Matrix &matrix) {
+  const int m = matrix.dim1();
+  const int n = matrix.dim2();
+
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      std::cout << matrix(i, j) << " ";
+    }
+    std::cout << std::endl;
+  }
+}
+
+void benchmarkComputeResult(benchmark::State &state) {
+  int iterations = state.range(0);
+  Matrix matrix = Matrix(1000, 1000);
+
+  fill_matrix(matrix, 10);
+
+  for (auto _ : state) {
+    gauss_seidel(matrix, iterations);
+    benchmark::DoNotOptimize(matrix);
+  }
+}
+
+int main(int argc, char **argv) {
+  ::benchmark::Initialize(&argc, argv);
+
+  for (int iterations = 0; iterations < 10000; iterations++) {
+    benchmark::RegisterBenchmark("bench_gaus_seidel", benchmarkComputeResult)
+        ->Arg(iterations)
+        ->Unit(benchmark::kMillisecond);
+  }
+
+  ::benchmark::RunSpecifiedBenchmarks();
+
+  return 0;
+}

lab07/results/kai_benchmark/kai_slurm.sh

@@ -0,0 +1,11 @@
+#!/bin/bash
+#SBATCH --job-name=kaiJacobiWave
+#SBATCH --output=kaiJacobiWaveSplitWork.out
+#SBATCH --ntasks=1
+#SBATCH --cpus-per-task=12
+
+# Load modules or activate environment if needed
+# module load gcc cmake ...
+
+# Run the benchmark binary
+./kaiJacobiWave --benchmark_out=kai_results.json --benchmark_out_format=json

lab07/results/kai_benchmark/matrix.h

@@ -0,0 +1,343 @@
+/**
+ * matrix.h a very simplistic class for m times n matrices.
+ */
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+// A very simplistic vector class for vectors of size n
+class Vector {
+ public:
+  // constructors
+  Vector(int n) : n_(n), data_(n_, 0) {}
+  Vector(const Vector& other) = default;
+  Vector(Vector&& other) = default;
+  ~Vector() = default;
+
+  // assignment operators
+  Vector& operator=(const Vector& other) = default;
+  Vector& operator=(Vector&& other) = default;
+
+  // element access
+  double& operator()(int i) { return data_[i]; }
+  const double& operator()(int i) const { return data_[i]; }
+
+  // getter functions for the dimensions
+  int dim() const { return n_; }
+
+  // comparison operators
+  bool operator==(const Vector& b) { return (data_ == b.data_); }
+  bool operator!=(const Vector& b) { return (data_ != b.data_); }
+
+  // addition
+  Vector& operator+=(const Vector& b) {
+    for (int i = 0; i < n_; ++i) {
+      operator()(i) += b(i);
+    }
+    return *this;
+  }
+
+  // subtraction
+  Vector& operator-=(const Vector& b) {
+    for (int i = 0; i < n_; ++i) {
+      operator()(i) -= b(i);
+    }
+    return *this;
+  }
+
+  // scalar multiplication
+  Vector& operator*=(double x) {
+    for (int i = 0; i < n_; ++i) {
+      operator()(i) *= x;
+    }
+    return *this;
+  }
+
+  // dot product between two vectors
+  double dot(const Vector& other) const {
+    double sum = 0;
+    for (int i = 0; i < n_; ++i) {
+      sum += operator()(i) * other(i);
+    }
+    return sum;
+  }
+
+ private:
+  int n_;                     // vector dimension
+  std::vector<double> data_;  // the vectors entries
+};
+
+inline double dot(const Vector& v1, const Vector& v2) {
+  return v1.dot(v2);
+}
+
+// Print the vector as a table
+inline std::ostream& operator<<(std::ostream& os, const Vector& a) {
+  const int width = 10;
+  const int precision = 4;
+
+  const auto originalPrecision = os.precision();
+  os << std::setprecision(precision);
+
+  for (int i = 0; i < a.dim(); ++i) {
+    os << std::setw(width) << a(i) << " ";
+  }
+
+  os << "\n";
+
+  os << std::setprecision(originalPrecision);
+  return os;
+}
+
+// A very simple class for m times n matrices
+class Matrix {
+ public:
+  // constructors
+  Matrix() : Matrix(0, 0) {}
+  Matrix(int m, int n) : m_(m), n_(n), data_(m_ * n_, 0) {}
+  Matrix(std::pair<int, int> dim) : Matrix(dim.first, dim.second) {}
+  Matrix(int n) : Matrix(n, n) {}
+  Matrix(const Matrix& other) = default;
+  Matrix(Matrix&& other) = default;
+  ~Matrix() = default;
+
+  // assignment operators
+  Matrix& operator=(const Matrix& other) = default;
+  Matrix& operator=(Matrix&& other) = default;
+
+  // element access
+  double& operator()(int i, int j) { return data_[i * n_ + j]; }
+  const double& operator()(int i, int j) const { return data_[i * n_ + j]; }
+
+  // getter functions for the dimensions
+  std::pair<int, int> dim() const { return std::pair<int, int>(m_, n_); }
+  int dim1() const { return m_; }
+  int dim2() const { return n_; }
+  int numEntries() const { return data_.size(); }
+
+  // comparison operators
+  bool operator==(const Matrix& b) { return (data_ == b.data_); }
+  bool operator!=(const Matrix& b) { return (data_ != b.data_); }
+
+  // addition
+  Matrix& operator+=(const Matrix& b) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) += b(i, j);
+      }
+    }
+    return *this;
+  }
+
+  // subtraction
+  Matrix& operator-=(const Matrix& b) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) -= b(i, j);
+      }
+    }
+    return *this;
+  }
+
+  // scalar multiplication
+  Matrix& operator*=(double x) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) *= x;
+      }
+    }
+    return *this;
+  }
+
+  // scalar division
+  Matrix& operator/=(double x) {
+    for (int i = 0; i < m_; ++i) {
+      for (int j = 0; j < n_; ++j) {
+        operator()(i, j) /= x;
+      }
+    }
+    return *this;
+  }
+
+  // matrix product (only for square matrices of equal dimension)
+  Matrix& operator*=(const Matrix& b) {
+    if (dim1() != dim2()) {
+      std::cout << "Error in matrix multiplication: no square matrix\n";
+    } else if (dim1() != b.dim1() || dim2() != b.dim2()) {
+      std::cout << "Error in matrix multiplication: dimensions do not match\n";
+    } else {
+      Matrix a = *this;
+      Matrix& c = *this;
+      const int m = dim1();
+      for (int i = 0; i < m; ++i) {
+        for (int j = 0; j < m; ++j) {
+          for (int k = 0; k < m; ++k) {
+            c(i, j) += a(i, k) * b(k, j);
+          }
+        }
+      }
+    }
+
+    return *this;
+  }
+
+ public:
+  int m_;                     // first dimension
+  int n_;                     // second dimension
+  std::vector<double> data_;  // the matrix' entries
+};
+
+// Print the matrix as a table
+inline std::ostream& operator<<(std::ostream& os, const Matrix& a) {
+  const int width = 10;
+  const int precision = 4;
+
+  const auto originalPrecision = os.precision();
+  os << std::setprecision(precision);
+
+  for (int i = 0; i < a.dim1(); ++i) {
+    for (int j = 0; j < a.dim2(); ++j) {
+      os << std::setw(width) << a(i, j) << " ";
+    }
+    if (i != a.dim1() - 1)
+      os << "\n";
+  }
+
+  os << std::setprecision(originalPrecision);
+  return os;
+}
+
+// matrix product
+inline Matrix operator*(const Matrix& a, const Matrix& b) {
+  if (a.dim2() == b.dim1()) {
+    int m = a.dim1();
+    int n = a.dim2();
+    int p = b.dim2();
+    Matrix c(m, p);
+    for (int i = 0; i < m; ++i) {
+      for (int j = 0; j < p; ++j) {
+        for (int k = 0; k < n; ++k) {
+          c(i, j) += a(i, k) * b(k, j);
+        }
+      }
+    }
+    return c;
+  } else {
+    return Matrix(0, 0);
+  }
+}
+
+inline bool equalWithinRange(const Matrix& a,
+                             const Matrix& b,
+                             double eps = 1e-12) {
+  if (a.dim1() != b.dim1() || a.dim2() != b.dim2())
+    return false;
+
+  int m = a.dim1();
+  int n = a.dim2();
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (fabs(a(i, j) - b(i, j)) > eps) {
+        return false;
+      }
+    }
+  }
+
+  return true;
+}
+
+// A very simple class for "3D-Matrices" (tensors) with dimension l x m x n
+class Matrix3D {
+ public:
+  // constructors
+  Matrix3D(int l, int m, int n) : l_(l), m_(m), n_(n), data_(l) {
+    for (int i = 0; i < l_; ++i) {
+      data_[i] = std::vector<std::vector<double>>(m_);
+      for (int j = 0; j < m_; ++j) {
+        data_[i][j] = std::vector<double>(n_, 0);
+      }
+    }
+  }
+  Matrix3D(int n) : Matrix3D(n, n, n) {}
+  Matrix3D(const Matrix3D& other) = default;
+  Matrix3D(Matrix3D&& other) = default;
+  ~Matrix3D() = default;
+
+  // assignment operators
+  Matrix3D& operator=(const Matrix3D& other) = default;
+  Matrix3D& operator=(Matrix3D&& other) = default;
+
+  // element access
+  double& operator()(int i, int j, int k) { return data_[i][j][k]; }
+  const double& operator()(int i, int j, int k) const { return data_[i][j][k]; }
+
+  // getter functions for the dimensions
+  int dim1() const { return l_; }
+  int dim2() const { return m_; }
+  int dim3() const { return n_; }
+
+  // comparison operators
+  bool operator==(const Matrix3D& b) { return (data_ == b.data_); }
+  bool operator!=(const Matrix3D& b) { return (data_ != b.data_); }
+
+  // addition
+  Matrix3D& operator+=(const Matrix3D& b) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) += b(i, j, k);
+        }
+      }
+    }
+    return *this;
+  }
+
+  // substraction
+  Matrix3D& operator-=(const Matrix3D& b) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) -= b(i, j, k);
+        }
+      }
+    }
+    return *this;
+  }
+
+  // scalar multiplication
+  Matrix3D& operator*=(double x) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) *= x;
+        }
+      }
+    }
+    return *this;
+  }
+
+  // scalar division
+  Matrix3D& operator/=(double x) {
+    for (int i = 0; i < l_; ++i) {
+      for (int j = 0; j < m_; ++j) {
+        for (int k = 0; k < n_; ++k) {
+          operator()(i, j, k) /= x;
+        }
+      }
+    }
+    return *this;
+  }
+
+ private:
+  int l_;                                               // first dimension
+  int m_;                                               // second dimension
+  int n_;                                               // third dimension
+  std::vector<std::vector<std::vector<double>>> data_;  // the tensors' entries
+};
+
+#endif  // MATRIX_H