diff --git a/app/benchmark_matrix/main.cpp b/app/benchmark_matrix/main.cpp
index ba02904..95cba2f 100644
--- a/app/benchmark_matrix/main.cpp
+++ b/app/benchmark_matrix/main.cpp
@@ -8,13 +8,13 @@ using namespace pls;
 
 using namespace comparison_benchmarks::base;
 
-template<typename T, int SIZE>
-class pls_matrix : public matrix::matrix<T, SIZE> {
+template<typename T>
+class pls_matrix : public matrix::matrix<T> {
  public:
-  pls_matrix() : matrix::matrix<T, SIZE>() {}
+  explicit pls_matrix(size_t size) : matrix::matrix<T>(size) {}
 
-  void multiply(const matrix::matrix<T, SIZE> &a, const matrix::matrix<T, SIZE> &b) override {
-    pls::algorithm::for_each_range(0, SIZE, [&](int i) {
+  void multiply(const matrix::matrix<T> &a, const matrix::matrix<T> &b) override {
+    pls::algorithm::for_each_range(0, this->size_, [&](int i) {
       this->multiply_column(i, a, b);
     });
   }
@@ -34,9 +34,9 @@ int main(int argc, char **argv) {
   runner.enable_memory_stats();
   runner.pre_allocate_stats();
 
-  pls_matrix<double, matrix::MATRIX_SIZE> a;
-  pls_matrix<double, matrix::MATRIX_SIZE> b;
-  pls_matrix<double, matrix::MATRIX_SIZE> result;
+  pls_matrix<double> a{matrix::MATRIX_SIZE};
+  pls_matrix<double> b{matrix::MATRIX_SIZE};
+  pls_matrix<double> result{matrix::MATRIX_SIZE};
 
   scheduler scheduler{(unsigned) num_threads, MAX_NUM_TASKS, MAX_STACK_SIZE};
 
diff --git a/extern/benchmark_base/include/benchmark_base/heat.h b/extern/benchmark_base/include/benchmark_base/heat.h
index 007d643..65eb2e8 100644
--- a/extern/benchmark_base/include/benchmark_base/heat.h
+++ b/extern/benchmark_base/include/benchmark_base/heat.h
@@ -2,7 +2,7 @@
 #ifndef COMPARISON_BENCHMARKS_BASE_HEAT_H
 #define COMPARISON_BENCHMARKS_BASE_HEAT_H
 
-#include <array>
+#include <vector>
 #include <iostream>
 #include <memory>
 
@@ -16,37 +16,39 @@ const int DIFFUSION_STEPS = 256;
 const int NUM_ITERATIONS = 100;
 const int WARMUP_ITERATIONS = 20;
 
-template<typename T, int SIZE>
+template<typename T>
 class heat_diffusion {
-  // Center portion is SIZExSIZE, borders are fixed temperature values
-  using matrix = std::array<std::array<T, SIZE + 2>, SIZE + 2>;
+  // Center portion is SIZExSIZE, borders are fixed temperature values.
+  // Should be a matrix of (SIZE + 2) x (SIZE + 2)
+  using matrix = std::vector<T>;
+ public:
+  size_t size_;
 
- protected:
   // Sane default values for the simulation (form paper).
   // This is not about perfect simulation results but the speedup of the workload.
-  double c = 0.1;
-  double d_s = 1.0 / (SIZE + 1);
-  double d_t = (d_s * d_s) / (4 * c);
+  double c_{0.1};
+  double d_s_{1.0 / (size_ + 1)};
+  double d_t_{(d_s_ * d_s_) / (4 * c_)};
 
- public:
-  matrix *current_data;
-  matrix *next_data;
+  matrix data_1_;
+  matrix data_2_;
 
-  explicit heat_diffusion() {
-    current_data = new matrix;
-    next_data = new matrix;
-    reset_data();
-  }
+  matrix *current_data_;
+  matrix *next_data_;
+
+  explicit heat_diffusion(size_t size) : size_{size},
+                                         data_1_((size + 2) * (size + 2)),
+                                         data_2_((size + 2) * (size + 2)) {
+    current_data_ = &data_1_;
+    next_data_ = &data_2_;
 
-  ~heat_diffusion() {
-    delete current_data;
-    delete next_data;
+    reset_data();
   }
 
   virtual void run_simulation(int n) {
     for (int i = 0; i < n; i++) {
-      for (int row = 1; row <= SIZE; row++) {
-        for (int column = 1; column <= SIZE; column++) {
+      for (int row = 1; row <= size_; row++) {
+        for (int column = 1; column <= size_; column++) {
           update_element(row, column);
         }
       }
@@ -56,40 +58,45 @@ class heat_diffusion {
     }
   }
 
- protected:
+  T &current_data_at(int row, int column) {
+    return (*current_data_)[row * (size_ + 2) + column];
+  }
+
+  T &next_data_at(int row, int column) {
+    return (*next_data_)[row * (size_ + 2) + column];
+  }
+
   void update_element(int row, int column) {
-    (*next_data)[row][column] = (*current_data)[row][column] + ((c * d_t) / (d_s * d_s)) *
-        ((*current_data)[row + 1][column] + (*current_data)[row - 1][column]
-            - 4 * (*current_data)[row][column]
-            + (*current_data)[row][column + 1] + (*current_data)[row][column - 1]);
+    next_data_at(row, column) = current_data_at(row, column) + ((c_ * d_t_) / (d_s_ * d_s_)) *
+        current_data_at(row + 1, column) + current_data_at(row - 1, column)
+        - 4 * current_data_at(row, column)
+        + current_data_at(row, column + 1) + current_data_at(row, column - 1);
   }
 
   void swap_data_arrays() {
-    matrix *tmp = current_data;
-    current_data = next_data;
-    next_data = tmp;
+    std::swap(current_data_, next_data_);
   }
 
   void reset_data() {
-    for (int row = 0; row < SIZE + 2; row++) {
-      for (int column = 0; column < SIZE + 2; column++) {
-        (*current_data)[row][column] = 0.0;
-        (*next_data)[row][column] = 0.0;
+    for (int row = 0; row < size_ + 2; row++) {
+      for (int column = 0; column < size_ + 2; column++) {
+        current_data_at(row, column) = 0.0;
+        next_data_at(row, column) = 0.0;
 
         // Edges are a fixed, hot temperature
-        if (row == 0 || row == SIZE + 1) {
-          (*current_data)[row][column] = 1.0;
-          (*next_data)[row][column] = 1.0;
+        if (row == 0 || row == size_ + 1) {
+          current_data_at(row, column) = 1.0;
+          next_data_at(row, column) = 1.0;
         }
       }
     }
   }
 };
 
-template<typename T, int SIZE>
-std::ostream &operator<<(std::ostream &strm, const heat_diffusion<T, SIZE> &simulation) {
-  for (int i = 0; i < SIZE + 2; i++) {
-    for (int j = 0; j < SIZE + 2; j++) {
+template<typename T>
+std::ostream &operator<<(std::ostream &strm, const heat_diffusion<T> &simulation) {
+  for (int i = 0; i < simulation.size_ + 2; i++) {
+    for (int j = 0; j < simulation.size_ + 2; j++) {
       // 'color' our output according to temperature
       char out;
       if (simulation.current_data[i][j] < 0.1) {
diff --git a/extern/benchmark_base/include/benchmark_base/matrix.h b/extern/benchmark_base/include/benchmark_base/matrix.h
index 1d8af13..3144baa 100644
--- a/extern/benchmark_base/include/benchmark_base/matrix.h
+++ b/extern/benchmark_base/include/benchmark_base/matrix.h
@@ -2,6 +2,7 @@
 #ifndef COMPARISON_BENCHMARKS_BASE_MATRIX_H
 #define COMPARISON_BENCHMARKS_BASE_MATRIX_H
 
+#include <vector>
 #include <algorithm>
 #include <iostream>
 
@@ -14,43 +15,51 @@ const int MATRIX_SIZE = 128;
 const int NUM_ITERATIONS = 5000;
 const int WARMUP_ITERATIONS = 1000;
 
-template<typename T, int SIZE>
+template<typename T>
 class matrix {
  public:
-  T data[SIZE][SIZE];
+  size_t size_;
+  std::vector<T> data_;
+
+  T &data_at(size_t row, size_t column) {
+    return data_[row * size_ + column];
+  }
+  const T &data_at(size_t row, size_t column) const {
+    return data_[row * size_ + column];
+  }
 
-  explicit matrix() {
-    for (int i = 0; i < SIZE; i++) {
-      for (int j = 0; j < SIZE; j++) {
-        data[i][j] = i;
+  explicit matrix(size_t size) : size_{size}, data_(size * size) {
+    for (int i = 0; i < size_; i++) {
+      for (int j = 0; j < size_; j++) {
+        data_at(i, j) = i;
       }
     }
   }
 
-  virtual void multiply(const matrix<T, SIZE> &a, const matrix<T, SIZE> &b) {
-    for (int i = 0; i < SIZE; i++) {
+  virtual void multiply(const matrix<T> &a, const matrix<T> &b) {
+    for (int i = 0; i < size_; i++) {
       multiply_column(i, a, b);
     }
   }
 
  protected:
-  void multiply_column(int i, const matrix<T, SIZE> &a, const matrix<T, SIZE> &b) {
-    for (int j = 0; j < SIZE; ++j) {
-      data[i][j] = 0;
+  void multiply_column(int i, const matrix<T> &a, const matrix<T> &b) {
+    for (int j = 0; j < size_; ++j) {
+      data_at(i, j) = 0;
     }
-    for (int k = 0; k < SIZE; ++k) {
-      for (int j = 0; j < SIZE; ++j) {
-        data[i][j] += a.data[i][k] * b.data[k][j];
+    for (int k = 0; k < size_; ++k) {
+      for (int j = 0; j < size_; ++j) {
+        data_at(i, j) += a.data_at(i, k) * b.data_at(i, k);
       }
     }
   }
 };
 
 template<typename T, int SIZE>
-std::ostream &operator<<(std::ostream &strm, const matrix<T, SIZE> &matrix) {
+std::ostream &operator<<(std::ostream &strm, const matrix<T> &matrix) {
   for (int i = 0; i < SIZE; i++) {
     for (int j = 0; j < SIZE; j++) {
-      strm << matrix.data[i][j] << "\t";
+      strm << matrix.data_at(i, j) << "\t";
     }
     strm << std::endl;
   }