From 5c19e5b084a884dba2c37980422bd9aec072893a Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 26 Aug 2021 11:29:42 +0200
Subject: [PATCH] Replace DenseMatrix by SmallMatrix

Fixes #24

On the way added
- SmallArray: which is intend to deal with "small" dynamic data sets
which follow Kokkos semantic
- SmallVector which is the companion vector class for DenseMatrix

These tools are made to define "small" linear systems (defined locally
on some items).
---
 src/algebra/EigenvalueSolver.cpp              |  20 +-
 src/algebra/EigenvalueSolver.hpp              |  26 +-
 src/algebra/LinearSolver.cpp                  |   3 +-
 src/algebra/LinearSolver.hpp                  |  15 +-
 src/algebra/PETScUtils.cpp                    |  20 +-
 src/algebra/PETScUtils.hpp                    |   6 +-
 .../{DenseMatrix.hpp => SmallMatrix.hpp}      | 127 ++---
 src/algebra/SmallVector.hpp                   | 216 ++++++++
 src/utils/SmallArray.hpp                      | 144 ++++++
 tests/CMakeLists.txt                          |   2 +-
 tests/test_DenseMatrix.cpp                    | 470 ------------------
 tests/test_EigenvalueSolver.cpp               |  26 +-
 tests/test_LinearSolver.cpp                   |  80 +--
 tests/test_PETScUtils.cpp                     |   8 +-
 tests/test_SmallArray.cpp                     | 285 +++++++++++
 15 files changed, 802 insertions(+), 646 deletions(-)
 rename src/algebra/{DenseMatrix.hpp => SmallMatrix.hpp} (64%)
 create mode 100644 src/algebra/SmallVector.hpp
 create mode 100644 src/utils/SmallArray.hpp
 delete mode 100644 tests/test_DenseMatrix.cpp
 create mode 100644 tests/test_SmallArray.cpp

diff --git a/src/algebra/EigenvalueSolver.cpp b/src/algebra/EigenvalueSolver.cpp
index 1d96d01f9..6125e67e6 100644
--- a/src/algebra/EigenvalueSolver.cpp
+++ b/src/algebra/EigenvalueSolver.cpp
@@ -63,7 +63,7 @@ struct EigenvalueSolver::Internals
 };
 
 void
-EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues)
+EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues)
 {
   EPS eps;
 
@@ -76,7 +76,7 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Arr
   PetscInt nb_eigenvalues;
   EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
 
-  eigenvalues = Array<double>(nb_eigenvalues);
+  eigenvalues = SmallArray<double>(nb_eigenvalues);
   for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
     EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, nullptr, nullptr);
   }
@@ -86,8 +86,8 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Arr
 
 void
 EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
-                                            Array<double>& eigenvalues,
-                                            std::vector<Vector<double>>& eigenvector_list)
+                                            SmallArray<double>& eigenvalues,
+                                            std::vector<SmallVector<double>>& eigenvector_list)
 {
   EPS eps;
 
@@ -100,11 +100,11 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
   PetscInt nb_eigenvalues;
   EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
 
-  eigenvalues = Array<double>(nb_eigenvalues);
+  eigenvalues = SmallArray<double>(nb_eigenvalues);
   eigenvector_list.reserve(nb_eigenvalues);
   for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
     Vec Vr;
-    Vector<double> eigenvector{A.numberOfRows()};
+    SmallVector<double> eigenvector{A.numberOfRows()};
     VecCreateSeqWithArray(PETSC_COMM_SELF, 1, A.numberOfRows(), &(eigenvector[0]), &Vr);
     EPSGetEigenpair(eps, i, &(eigenvalues[i]), nullptr, Vr, nullptr);
     VecDestroy(&Vr);
@@ -116,8 +116,8 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
 
 void
 EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
-                                            Array<double>& eigenvalues,
-                                            DenseMatrix<double>& P)
+                                            SmallArray<double>& eigenvalues,
+                                            SmallMatrix<double>& P)
 {
   EPS eps;
 
@@ -130,8 +130,8 @@ EigenvalueSolver::computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
   PetscInt nb_eigenvalues;
   EPSGetDimensions(eps, &nb_eigenvalues, nullptr, nullptr);
 
-  eigenvalues = Array<double>(nb_eigenvalues);
-  P           = DenseMatrix<double>(nb_eigenvalues, nb_eigenvalues);
+  eigenvalues = SmallArray<double>(nb_eigenvalues);
+  P           = SmallMatrix<double>(nb_eigenvalues, nb_eigenvalues);
 
   Array<double> eigenvector(nb_eigenvalues);
   for (PetscInt i = 0; i < nb_eigenvalues; ++i) {
diff --git a/src/algebra/EigenvalueSolver.hpp b/src/algebra/EigenvalueSolver.hpp
index a9dd941b9..a9b462509 100644
--- a/src/algebra/EigenvalueSolver.hpp
+++ b/src/algebra/EigenvalueSolver.hpp
@@ -3,10 +3,10 @@
 
 #include <algebra/PETScUtils.hpp>
 
-#include <algebra/DenseMatrix.hpp>
-#include <algebra/Vector.hpp>
-#include <utils/Array.hpp>
+#include <algebra/SmallMatrix.hpp>
+#include <algebra/SmallVector.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/SmallArray.hpp>
 
 class EigenvalueSolver
 {
@@ -14,19 +14,21 @@ class EigenvalueSolver
   struct Internals;
 
 #ifdef PUGS_HAS_SLEPC
-  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues);
+  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues);
 
   void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
-                                 Array<double>& eigenvalues,
-                                 std::vector<Vector<double>>& eigenvectors);
+                                 SmallArray<double>& eigenvalues,
+                                 std::vector<SmallVector<double>>& eigenvectors);
 
-  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, Array<double>& eigenvalues, DenseMatrix<double>& P);
+  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
+                                 SmallArray<double>& eigenvalues,
+                                 SmallMatrix<double>& P);
 #endif   // PUGS_HAS_SLEPC
 
  public:
   template <typename MatrixType>
   void
-  computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] Array<double>& eigenvalues)
+  computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallArray<double>& eigenvalues)
   {
 #ifdef PUGS_HAS_SLEPC
     this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues);
@@ -38,8 +40,8 @@ class EigenvalueSolver
   template <typename MatrixType>
   void
   computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
-                            [[maybe_unused]] Array<double>& eigenvalues,
-                            [[maybe_unused]] std::vector<Vector<double>>& eigenvectors)
+                            [[maybe_unused]] SmallArray<double>& eigenvalues,
+                            [[maybe_unused]] std::vector<SmallVector<double>>& eigenvectors)
   {
 #ifdef PUGS_HAS_SLEPC
     this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, eigenvectors);
@@ -51,8 +53,8 @@ class EigenvalueSolver
   template <typename MatrixType>
   void
   computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
-                            [[maybe_unused]] Array<double>& eigenvalues,
-                            [[maybe_unused]] DenseMatrix<double>& P)
+                            [[maybe_unused]] SmallArray<double>& eigenvalues,
+                            [[maybe_unused]] SmallMatrix<double>& P)
   {
 #ifdef PUGS_HAS_SLEPC
     this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, P);
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index d7a1d4119..99bff927c 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -224,6 +224,7 @@ struct LinearSolver::Internals
     case LSMethod::cholesky: {
       KSPSetType(ksp, KSPPREONLY);
       PCSetType(pc, PCCHOLESKY);
+      PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
       direct_solver = true;
       break;
     }
@@ -335,7 +336,7 @@ LinearSolver::solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>&
 }
 
 void
-LinearSolver::solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b)
+LinearSolver::solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b)
 {
   Internals::solveLocalSystem(A, x, b, m_options);
 }
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index 1ab43a140..e90402c8d 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -2,8 +2,8 @@
 #define LINEAR_SOLVER_HPP
 
 #include <algebra/CRSMatrix.hpp>
-#include <algebra/DenseMatrix.hpp>
 #include <algebra/LinearSolverOptions.hpp>
+#include <algebra/SmallMatrix.hpp>
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
 #include <algebra/Vector.hpp>
@@ -24,13 +24,10 @@ class LinearSolver
   void
   solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
   {
-    DenseMatrix A_dense{A};
-    Vector<double> x_vector{N};
-    Vector<double> b_vector{N};
-    for (size_t i = 0; i < N; ++i) {
-      x_vector[i] = x[i];
-      b_vector[i] = b[i];
-    }
+    SmallMatrix A_dense{A};
+
+    SmallVector x_vector{x};
+    SmallVector b_vector{b};
 
     this->solveLocalSystem(A_dense, x_vector, b_vector);
 
@@ -40,7 +37,7 @@ class LinearSolver
   }
 
   void solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b);
-  void solveLocalSystem(const DenseMatrix<double>& A, Vector<double>& x, const Vector<const double>& b);
+  void solveLocalSystem(const SmallMatrix<double>& A, SmallVector<double>& x, const SmallVector<const double>& b);
 
   LinearSolver();
   LinearSolver(const LinearSolverOptions& options);
diff --git a/src/algebra/PETScUtils.cpp b/src/algebra/PETScUtils.cpp
index 63c5a68ac..299041035 100644
--- a/src/algebra/PETScUtils.cpp
+++ b/src/algebra/PETScUtils.cpp
@@ -12,24 +12,25 @@ PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const size_t nb_rows, const size_
   MatSetUp(m_petscMat);
 
   {
-    Array<PetscInt> row_indices(nb_rows);
+    SmallArray<PetscInt> row_indices(nb_rows);
     for (size_t i = 0; i < nb_rows; ++i) {
       row_indices[i] = i;
     }
-    m_row_indices = row_indices;
+    m_small_row_indices = row_indices;
   }
 
   if (nb_rows == nb_columns) {
-    m_column_indices = m_row_indices;
+    m_small_column_indices = m_small_row_indices;
   } else {
-    Array<PetscInt> column_indices(nb_columns);
+    SmallArray<PetscInt> column_indices(nb_columns);
     for (size_t i = 0; i < nb_columns; ++i) {
       column_indices[i] = i;
     }
-    m_column_indices = column_indices;
+    m_small_column_indices = column_indices;
   }
 
-  MatSetValuesBlocked(m_petscMat, nb_rows, &(m_row_indices[0]), nb_columns, &(m_column_indices[0]), A, INSERT_VALUES);
+  MatSetValuesBlocked(m_petscMat, nb_rows, &(m_small_row_indices[0]), nb_columns, &(m_small_column_indices[0]), A,
+                      INSERT_VALUES);
 
   MatAssemblyBegin(m_petscMat, MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
@@ -42,9 +43,10 @@ PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A)
                 "index type conversion is not implemented yet");
   static_assert(std::is_same_v<PetscScalar, typename std::decay_t<decltype(A)>::data_type>,
                 "value type conversion is not implemented yet");
-  m_values         = copy(A.values());
-  m_row_indices    = copy(A.rowMap());
-  m_column_indices = copy(A.columnIndices());
+
+  m_values         = A.values();
+  m_row_indices    = A.rowMap();
+  m_column_indices = A.columnIndices();
 
   MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.numberOfRows(), A.numberOfColumns(),
                             const_cast<PetscInt*>(&m_row_indices[0]),      // These ugly const_cast
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index 3e7cff5b8..c4db9ee59 100644
--- a/src/algebra/PETScUtils.hpp
+++ b/src/algebra/PETScUtils.hpp
@@ -6,7 +6,7 @@
 #ifdef PUGS_HAS_PETSC
 
 #include <algebra/CRSMatrix.hpp>
-#include <algebra/DenseMatrix.hpp>
+#include <algebra/SmallMatrix.hpp>
 #include <algebra/TinyMatrix.hpp>
 
 #include <petsc.h>
@@ -18,6 +18,8 @@ class PETScAijMatrixEmbedder
   Array<const PetscScalar> m_values;
   Array<const PetscInt> m_row_indices;
   Array<const PetscInt> m_column_indices;
+  SmallArray<const PetscInt> m_small_row_indices;
+  SmallArray<const PetscInt> m_small_column_indices;
   const size_t m_nb_rows;
   const size_t m_nb_columns;
 
@@ -54,7 +56,7 @@ class PETScAijMatrixEmbedder
   PETScAijMatrixEmbedder(const TinyMatrix<N>& A) : PETScAijMatrixEmbedder{N, N, &A(0, 0)}
   {}
 
-  PETScAijMatrixEmbedder(const DenseMatrix<double>& A)
+  PETScAijMatrixEmbedder(const SmallMatrix<double>& A)
     : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
   {}
 
diff --git a/src/algebra/DenseMatrix.hpp b/src/algebra/SmallMatrix.hpp
similarity index 64%
rename from src/algebra/DenseMatrix.hpp
rename to src/algebra/SmallMatrix.hpp
index 15967441b..d0d6ee329 100644
--- a/src/algebra/DenseMatrix.hpp
+++ b/src/algebra/SmallMatrix.hpp
@@ -1,18 +1,18 @@
-#ifndef DENSE_MATRIX_HPP
-#define DENSE_MATRIX_HPP
+#ifndef SMALL_MATRIX_HPP
+#define SMALL_MATRIX_HPP
 
+#include <algebra/SmallVector.hpp>
 #include <algebra/TinyMatrix.hpp>
-#include <algebra/Vector.hpp>
-#include <utils/Array.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/PugsUtils.hpp>
+#include <utils/SmallArray.hpp>
 #include <utils/Types.hpp>
 
 #include <iostream>
 
 template <typename DataType>
-class DenseMatrix   // LCOV_EXCL_LINE
+class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
 {
  public:
   using data_type  = DataType;
@@ -21,32 +21,29 @@ class DenseMatrix   // LCOV_EXCL_LINE
  private:
   size_t m_nb_rows;
   size_t m_nb_columns;
-  Array<DataType> m_values;
+  SmallArray<DataType> m_values;
 
   static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
-  static_assert(std::is_arithmetic_v<DataType>, "Dense matrices expect arithmetic data");
+  static_assert(std::is_arithmetic_v<DataType>, "Small matrices expect arithmetic data");
 
   // Allows const version to access our data
-  friend DenseMatrix<std::add_const_t<DataType>>;
+  friend SmallMatrix<std::add_const_t<DataType>>;
 
  public:
   PUGS_INLINE
-  bool
-  isSquare() const noexcept
+  bool isSquare() const noexcept
   {
     return m_nb_rows == m_nb_columns;
   }
 
-  friend PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
-  copy(const DenseMatrix& A) noexcept
+  friend PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> copy(const SmallMatrix& A) noexcept
   {
     return {A.m_nb_rows, A.m_nb_columns, copy(A.m_values)};
   }
 
-  friend PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
-  transpose(const DenseMatrix& A)
+  friend PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> transpose(const SmallMatrix& A)
   {
-    DenseMatrix<std::remove_const_t<DataType>> A_transpose{A.m_nb_columns, A.m_nb_rows};
+    SmallMatrix<std::remove_const_t<DataType>> A_transpose{A.m_nb_columns, A.m_nb_rows};
     for (size_t i = 0; i < A.m_nb_rows; ++i) {
       for (size_t j = 0; j < A.m_nb_columns; ++j) {
         A_transpose(j, i) = A(i, j);
@@ -55,20 +52,20 @@ class DenseMatrix   // LCOV_EXCL_LINE
     return A_transpose;
   }
 
-  friend PUGS_INLINE DenseMatrix operator*(const DataType& a, const DenseMatrix& A)
+  friend PUGS_INLINE SmallMatrix operator*(const DataType& a, const SmallMatrix& A)
   {
-    DenseMatrix<std::remove_const_t<DataType>> aA = copy(A);
+    SmallMatrix<std::remove_const_t<DataType>> aA = copy(A);
     return aA *= a;
   }
 
   template <typename DataType2>
-  PUGS_INLINE Vector<std::remove_const_t<DataType>> operator*(const Vector<DataType2>& x) const
+  PUGS_INLINE SmallVector<std::remove_const_t<DataType>> operator*(const SmallVector<DataType2>& x) const
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
     Assert(m_nb_columns == x.size(), "cannot compute matrix-vector product: incompatible sizes");
-    const DenseMatrix& A = *this;
-    Vector<std::remove_const_t<DataType>> Ax{m_nb_rows};
+    const SmallMatrix& A = *this;
+    SmallVector<std::remove_const_t<DataType>> Ax{m_nb_rows};
     for (size_t i = 0; i < m_nb_rows; ++i) {
       std::remove_const_t<DataType> Axi = A(i, 0) * x[0];
       for (size_t j = 1; j < m_nb_columns; ++j) {
@@ -80,13 +77,13 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>> operator*(const DenseMatrix<DataType2>& B) const
+  PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> operator*(const SmallMatrix<DataType2>& B) const
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
     Assert(m_nb_columns == B.numberOfRows(), "cannot compute matrix product: incompatible sizes");
-    const DenseMatrix& A = *this;
-    DenseMatrix<std::remove_const_t<DataType>> AB{m_nb_rows, B.numberOfColumns()};
+    const SmallMatrix& A = *this;
+    SmallMatrix<std::remove_const_t<DataType>> AB{m_nb_rows, B.numberOfColumns()};
 
     for (size_t i = 0; i < m_nb_rows; ++i) {
       for (size_t j = 0; j < B.numberOfColumns(); ++j) {
@@ -101,16 +98,14 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix&
-  operator/=(const DataType2& a)
+  PUGS_INLINE SmallMatrix& operator/=(const DataType2& a)
   {
     const auto inv_a = 1. / a;
     return (*this) *= inv_a;
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix&
-  operator*=(const DataType2& a)
+  PUGS_INLINE SmallMatrix& operator*=(const DataType2& a)
   {
     parallel_for(
       m_values.size(), PUGS_LAMBDA(index_type i) { m_values[i] *= a; });
@@ -118,8 +113,7 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix&
-  operator-=(const DenseMatrix<DataType2>& B)
+  PUGS_INLINE SmallMatrix& operator-=(const SmallMatrix<DataType2>& B)
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
@@ -132,8 +126,7 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix&
-  operator+=(const DenseMatrix<DataType2>& B)
+  PUGS_INLINE SmallMatrix& operator+=(const SmallMatrix<DataType2>& B)
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
@@ -145,15 +138,14 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
-  operator+(const DenseMatrix<DataType2>& B) const
+  PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> operator+(const SmallMatrix<DataType2>& B) const
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
     Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
            "cannot compute matrix sum: incompatible sizes");
 
-    DenseMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
+    SmallMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
 
     parallel_for(
       m_values.size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + B.m_values[i]; });
@@ -162,15 +154,14 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix<std::remove_const_t<DataType>>
-  operator-(const DenseMatrix<DataType2>& B) const
+  PUGS_INLINE SmallMatrix<std::remove_const_t<DataType>> operator-(const SmallMatrix<DataType2>& B) const
   {
     static_assert(std::is_same_v<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>,
                   "incompatible data types");
     Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
            "cannot compute matrix difference: incompatible sizes");
 
-    DenseMatrix<std::remove_const_t<DataType>> difference{B.numberOfRows(), B.numberOfColumns()};
+    SmallMatrix<std::remove_const_t<DataType>> difference{B.numberOfRows(), B.numberOfColumns()};
 
     parallel_for(
       m_values.size(), PUGS_LAMBDA(index_type i) { difference.m_values[i] = m_values[i] - B.m_values[i]; });
@@ -179,40 +170,36 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   PUGS_INLINE
-  DataType&
-  operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
+  DataType& operator()(index_type i, index_type j) const noexcept(NO_ASSERT)
   {
     Assert(i < m_nb_rows and j < m_nb_columns, "cannot access element: invalid indices");
     return m_values[i * m_nb_columns + j];
   }
 
   PUGS_INLINE
-  size_t
-  numberOfRows() const noexcept
+  size_t numberOfRows() const noexcept
   {
     return m_nb_rows;
   }
 
   PUGS_INLINE
-  size_t
-  numberOfColumns() const noexcept
+  size_t numberOfColumns() const noexcept
   {
     return m_nb_columns;
   }
 
-  PUGS_INLINE void
-  fill(const DataType& value) noexcept
+  PUGS_INLINE void fill(const DataType& value) noexcept
   {
     m_values.fill(value);
   }
 
-  PUGS_INLINE DenseMatrix& operator=(ZeroType) noexcept
+  PUGS_INLINE SmallMatrix& operator=(ZeroType) noexcept
   {
     m_values.fill(0);
     return *this;
   }
 
-  PUGS_INLINE DenseMatrix& operator=(IdentityType) noexcept(NO_ASSERT)
+  PUGS_INLINE SmallMatrix& operator=(IdentityType) noexcept(NO_ASSERT)
   {
     Assert(m_nb_rows == m_nb_columns, "identity must be a square matrix");
 
@@ -223,15 +210,14 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  PUGS_INLINE DenseMatrix&
-  operator=(const DenseMatrix<DataType2>& A) noexcept
+  PUGS_INLINE SmallMatrix& operator=(const SmallMatrix<DataType2>& A) noexcept
   {
     // ensures that DataType is the same as source DataType2
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
-                  "Cannot assign DenseMatrix of different type");
+                  "Cannot assign SmallMatrix of different type");
     // ensures that const is not lost through copy
     static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
-                  "Cannot assign DenseMatrix of const to DenseMatrix of non-const");
+                  "Cannot assign SmallMatrix of const to SmallMatrix of non-const");
 
     m_nb_rows    = A.m_nb_rows;
     m_nb_columns = A.m_nb_columns;
@@ -240,13 +226,12 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   PUGS_INLINE
-  DenseMatrix& operator=(const DenseMatrix&) = default;
+  SmallMatrix& operator=(const SmallMatrix&) = default;
 
   PUGS_INLINE
-  DenseMatrix& operator=(DenseMatrix&&) = default;
+  SmallMatrix& operator=(SmallMatrix&&) = default;
 
-  friend std::ostream&
-  operator<<(std::ostream& os, const DenseMatrix& A)
+  friend std::ostream& operator<<(std::ostream& os, const SmallMatrix& A)
   {
     for (size_t i = 0; i < A.numberOfRows(); ++i) {
       os << i << '|';
@@ -259,50 +244,42 @@ class DenseMatrix   // LCOV_EXCL_LINE
   }
 
   template <typename DataType2>
-  DenseMatrix(const DenseMatrix<DataType2>& A)
+  SmallMatrix(const SmallMatrix<DataType2>& A)
   {
     // ensures that DataType is the same as source DataType2
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
-                  "Cannot assign DenseMatrix of different type");
+                  "Cannot assign SmallMatrix of different type");
     // ensures that const is not lost through copy
     static_assert(((std::is_const<DataType2>() and std::is_const<DataType>()) or not std::is_const<DataType2>()),
-                  "Cannot assign DenseMatrix of const to DenseMatrix of non-const");
+                  "Cannot assign SmallMatrix of const to SmallMatrix of non-const");
 
     this->operator=(A);
   }
 
-  DenseMatrix(const DenseMatrix&) = default;
+  SmallMatrix(const SmallMatrix&) = default;
 
-  DenseMatrix(DenseMatrix&&) = default;
+  SmallMatrix(SmallMatrix &&) = default;
 
-  explicit DenseMatrix(size_t nb_rows, size_t nb_columns) noexcept
+  explicit SmallMatrix(size_t nb_rows, size_t nb_columns) noexcept
     : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{nb_rows * nb_columns}
   {}
 
-  explicit DenseMatrix(size_t nb_rows) noexcept : m_nb_rows{nb_rows}, m_nb_columns{nb_rows}, m_values{nb_rows * nb_rows}
+  explicit SmallMatrix(size_t nb_rows) noexcept : m_nb_rows{nb_rows}, m_nb_columns{nb_rows}, m_values{nb_rows * nb_rows}
   {}
 
   template <size_t M, size_t N>
-  explicit DenseMatrix(const TinyMatrix<M, N, DataType>& A) noexcept : m_nb_rows{M}, m_nb_columns{N}, m_values{M * N}
+  explicit SmallMatrix(const TinyMatrix<M, N, DataType>& A) noexcept : m_nb_rows{M}, m_nb_columns{N}, m_values{M * N}
   {
     for (size_t i = 0; i < M; ++i) {
       for (size_t j = 0; j < N; ++j) {
-        m_values[i * M + j] = A(i, j);
+        m_values[i * N + j] = A(i, j);
       }
     }
   }
 
-  DenseMatrix() noexcept : m_nb_rows{0}, m_nb_columns{0} {}
+  SmallMatrix() noexcept : m_nb_rows{0}, m_nb_columns{0} {}
 
- private:
-  DenseMatrix(size_t nb_rows, size_t nb_columns, const Array<DataType> values) noexcept(NO_ASSERT)
-    : m_nb_rows{nb_rows}, m_nb_columns{nb_columns}, m_values{values}
-  {
-    Assert(m_values.size() == m_nb_rows * m_nb_columns, "cannot create matrix: incorrect number of values");
-  }
-
- public:
-  ~DenseMatrix() = default;
+  ~SmallMatrix() = default;
 };
 
-#endif   // DENSE_MATRIX_HPP
+#endif   // SMALL_MATRIX_HPP
diff --git a/src/algebra/SmallVector.hpp b/src/algebra/SmallVector.hpp
new file mode 100644
index 000000000..d20bb48f4
--- /dev/null
+++ b/src/algebra/SmallVector.hpp
@@ -0,0 +1,216 @@
+#ifndef SMALL_VECTOR_HPP
+#define SMALL_VECTOR_HPP
+
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+#include <utils/SmallArray.hpp>
+#include <utils/Types.hpp>
+
+template <typename DataType>
+class SmallVector   // LCOV_EXCL_LINE
+{
+ public:
+  using data_type  = DataType;
+  using index_type = size_t;
+
+ private:
+  SmallArray<DataType> m_values;
+
+  static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
+  static_assert(std::is_arithmetic_v<DataType>, "SmallVector is only defined for arithmetic values");
+
+  // Allows const version to access our data
+  friend SmallVector<std::add_const_t<DataType>>;
+
+ public:
+  friend std::ostream&
+  operator<<(std::ostream& os, const SmallVector& x)
+  {
+    return os << x.m_values;
+  }
+
+  friend SmallVector<std::remove_const_t<DataType>>
+  copy(const SmallVector& x)
+  {
+    auto values = copy(x.m_values);
+
+    SmallVector<std::remove_const_t<DataType>> x_copy{0};
+    x_copy.m_values = values;
+    return x_copy;
+  }
+
+  friend SmallVector operator*(const DataType& a, const SmallVector& x)
+  {
+    SmallVector<std::remove_const_t<DataType>> y = copy(x);
+    return y *= a;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE friend auto
+  dot(const SmallVector& x, const SmallVector<DataType2>& y)
+  {
+    Assert(x.size() == y.size(), "cannot compute dot product: incompatible vector sizes");
+    // Quite ugly, TODO: fix me in C++20
+    auto promoted = [] {
+      DataType a{0};
+      DataType2 b{0};
+      return a * b;
+    }();
+
+    decltype(promoted) sum = 0;
+
+    // Does not use parallel_for to preserve sum order
+    for (index_type i = 0; i < x.size(); ++i) {
+      sum += x[i] * y[i];
+    }
+
+    return sum;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector&
+  operator/=(const DataType2& a)
+  {
+    const auto inv_a = 1. / a;
+    return (*this) *= inv_a;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector&
+  operator*=(const DataType2& a)
+  {
+    parallel_for(
+      this->size(), PUGS_LAMBDA(index_type i) { m_values[i] *= a; });
+    return *this;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector&
+  operator-=(const SmallVector<DataType2>& y)
+  {
+    Assert(this->size() == y.size(), "cannot substract vector: incompatible sizes");
+
+    parallel_for(
+      this->size(), PUGS_LAMBDA(index_type i) { m_values[i] -= y[i]; });
+
+    return *this;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector&
+  operator+=(const SmallVector<DataType2>& y)
+  {
+    Assert(this->size() == y.size(), "cannot add vector: incompatible sizes");
+
+    parallel_for(
+      this->size(), PUGS_LAMBDA(index_type i) { m_values[i] += y[i]; });
+
+    return *this;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector<std::remove_const_t<DataType>>
+  operator+(const SmallVector<DataType2>& y) const
+  {
+    Assert(this->size() == y.size(), "cannot compute vector sum: incompatible sizes");
+    SmallVector<std::remove_const_t<DataType>> sum{y.size()};
+
+    parallel_for(
+      this->size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + y[i]; });
+
+    return sum;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector<std::remove_const_t<DataType>>
+  operator-(const SmallVector<DataType2>& y) const
+  {
+    Assert(this->size() == y.size(), "cannot compute vector difference: incompatible sizes");
+    SmallVector<std::remove_const_t<DataType>> sum{y.size()};
+
+    parallel_for(
+      this->size(), PUGS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] - y[i]; });
+
+    return sum;
+  }
+
+  PUGS_INLINE
+  DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  {
+    return m_values[i];
+  }
+
+  PUGS_INLINE
+  size_t
+  size() const noexcept
+  {
+    return m_values.size();
+  }
+
+  PUGS_INLINE SmallVector&
+  fill(const DataType& value) noexcept
+  {
+    m_values.fill(value);
+    return *this;
+  }
+
+  PUGS_INLINE SmallVector& operator=(const DataType&) = delete;
+
+  PUGS_INLINE SmallVector&
+  operator=(const ZeroType&) noexcept
+  {
+    m_values.fill(0);
+    return *this;
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallVector&
+  operator=(const SmallVector<DataType2>& x)
+  {
+    // ensures that DataType is the same as source DataType2
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot assign SmallVector of different type");
+    // ensures that const is not lost through copy. If this point is
+    // reached data types only differ through their constness
+    static_assert((not std::is_const_v<DataType2>), "Cannot assign SmallVector of const to SmallVector of non-const");
+
+    m_values = x.m_values;
+    return *this;
+  }
+
+  PUGS_INLINE
+  SmallVector& operator=(const SmallVector&) = default;
+
+  PUGS_INLINE
+  SmallVector& operator=(SmallVector&& x) = default;
+
+  template <typename DataType2>
+  SmallVector(const SmallVector<DataType2>& x)
+  {
+    // ensures that DataType is the same as source DataType2
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot assign SmallVector of different type");
+    // ensures that const is not lost through copy
+    static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
+                  "Cannot assign SmallVector of const to SmallVector of non-const");
+
+    m_values = x.m_values;
+  }
+
+  template <size_t N>
+  explicit SmallVector(const TinyVector<N, DataType>& x) : m_values{N}
+  {
+    std::copy(&x[0], &x[0] + N, &m_values[0]);
+  }
+
+  SmallVector(const SmallVector&) = default;
+  SmallVector(SmallVector&&)      = default;
+
+  SmallVector(size_t size) : m_values{size} {}
+
+  SmallVector()  = default;
+  ~SmallVector() = default;
+};
+
+#endif   // SMALL_VECTOR_HPP
diff --git a/src/utils/SmallArray.hpp b/src/utils/SmallArray.hpp
new file mode 100644
index 000000000..b5186a7f6
--- /dev/null
+++ b/src/utils/SmallArray.hpp
@@ -0,0 +1,144 @@
+#ifndef SMALL_ARRAY_HPP
+#define SMALL_ARRAY_HPP
+
+#include <utils/InvalidData.hpp>
+#include <utils/NaNHelper.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+#include <utils/PugsUtils.hpp>
+
+#include <algorithm>
+
+template <typename DataType>
+class [[nodiscard]] SmallArray
+{
+ public:
+  using data_type  = DataType;
+  using index_type = size_t;
+
+ private:
+  index_type m_size;
+  std::shared_ptr<DataType[]> m_values;
+
+  // Allows const version to access our data
+  friend SmallArray<std::add_const_t<DataType>>;
+
+ public:
+  PUGS_INLINE size_t size() const noexcept
+  {
+    return m_size;
+  }
+
+  friend PUGS_INLINE SmallArray<std::remove_const_t<DataType>> copy(const SmallArray<DataType>& source)
+  {
+    SmallArray<std::remove_const_t<DataType>> image(source.m_size);
+    std::copy(source.m_values.get(), source.m_values.get() + source.m_size, image.m_values.get());
+
+    return image;
+  }
+
+  friend PUGS_INLINE void copy_to(const SmallArray<DataType>& source,
+                                  const SmallArray<std::remove_const_t<DataType>>& destination)
+  {
+    Assert(source.size() == destination.size());
+    std::copy(source.m_values.get(), source.m_values.get() + source.m_size, destination.m_values.get());
+  }
+
+  PUGS_INLINE DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  {
+    Assert(i < m_size);
+    return m_values[i];
+  }
+
+  PUGS_INLINE
+  void fill(const DataType& data) const
+  {
+    static_assert(not std::is_const_v<DataType>, "Cannot modify SmallArray of const");
+    std::fill(m_values.get(), m_values.get() + m_size, data);
+  }
+
+  template <typename DataType2>
+  PUGS_INLINE SmallArray& operator=(const SmallArray<DataType2>& array) noexcept
+  {
+    // ensures that DataType is the same as source DataType2
+    static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+                  "Cannot assign SmallArray of different type");
+    // ensures that const is not lost through copy
+    static_assert(((std::is_const_v<DataType2> and std::is_const_v<DataType>) or not std::is_const_v<DataType2>),
+                  "Cannot assign SmallArray of const to SmallArray of non-const");
+
+    m_size   = array.size();
+    m_values = array.m_values;
+    return *this;
+  }
+
+  PUGS_INLINE
+  SmallArray& operator=(const SmallArray&) = default;
+
+  PUGS_INLINE
+  SmallArray& operator=(SmallArray&&) = default;
+
+  PUGS_INLINE
+  explicit SmallArray(size_t size) : m_size{size}, m_values{std::shared_ptr<DataType[]>(new DataType[size])}
+  {
+    static_assert(not std::is_const<DataType>(), "Cannot allocate SmallArray of const data: only view is "
+                                                 "supported");
+
+#ifndef NDEBUG
+    if constexpr (not std::is_const_v<DataType>) {
+      using T = std::decay_t<DataType>;
+      if constexpr (std::is_arithmetic_v<T>) {
+        this->fill(invalid_data_v<T>);
+      } else if constexpr (is_tiny_vector_v<T>) {
+        this->fill(T{});
+      } else if constexpr (is_tiny_matrix_v<T>) {
+        this->fill(T{});
+      }
+    }
+#endif   // NDEBUG
+  }
+
+  friend std::ostream& operator<<(std::ostream& os, const SmallArray& x)
+  {
+    if (x.size() > 0) {
+      os << 0 << ':' << NaNHelper(x[0]);
+    }
+    for (size_t i = 1; i < x.size(); ++i) {
+      os << ' ' << i << ':' << NaNHelper(x[i]);
+    }
+    return os;
+  }
+
+  PUGS_INLINE
+  SmallArray() = default;
+
+  PUGS_INLINE
+  SmallArray(const SmallArray&) = default;
+
+  template <typename DataType2>
+  PUGS_INLINE SmallArray(const SmallArray<DataType2>& array) noexcept
+  {
+    this->operator=(array);
+  }
+
+  PUGS_INLINE
+  SmallArray(SmallArray &&) = default;
+
+  PUGS_INLINE
+  ~SmallArray() = default;
+};
+
+// map, multimap, unordered_map and stack cannot be copied this way
+template <typename Container>
+PUGS_INLINE SmallArray<std::remove_const_t<typename Container::value_type>>
+convert_to_small_array(const Container& given_container)
+{
+  using DataType = typename Container::value_type;
+  SmallArray<std::remove_const_t<DataType>> array(given_container.size());
+  if (given_container.size() > 0) {
+    std::copy(begin(given_container), end(given_container), &(array[0]));
+  }
+  return array;
+}
+
+#endif   // SMALL_ARRAY_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 714bb23bf..4bac18760 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -59,7 +59,6 @@ add_executable (unit_tests
   test_CRSMatrix.cpp
   test_CRSMatrixDescriptor.cpp
   test_DataVariant.cpp
-  test_DenseMatrix.cpp
   test_Demangle.cpp
   test_DiscreteFunctionDescriptorP0.cpp
   test_DiscreteFunctionDescriptorP0Vector.cpp
@@ -93,6 +92,7 @@ add_executable (unit_tests
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp
   test_RevisionInfo.cpp
+  test_SmallArray.cpp
   test_SymbolTable.cpp
   test_Table.cpp
   test_Timer.cpp
diff --git a/tests/test_DenseMatrix.cpp b/tests/test_DenseMatrix.cpp
deleted file mode 100644
index b794689a0..000000000
--- a/tests/test_DenseMatrix.cpp
+++ /dev/null
@@ -1,470 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <utils/PugsAssert.hpp>
-
-#include <algebra/DenseMatrix.hpp>
-#include <algebra/Vector.hpp>
-
-// Instantiate to ensure full coverage is performed
-template class DenseMatrix<int>;
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("DenseMatrix", "[algebra]")
-{
-  SECTION("size")
-  {
-    DenseMatrix<int> A{2, 3};
-    REQUIRE(A.numberOfRows() == 2);
-    REQUIRE(A.numberOfColumns() == 3);
-  }
-
-  SECTION("is square")
-  {
-    REQUIRE(not DenseMatrix<int>{2, 3}.isSquare());
-    REQUIRE(not DenseMatrix<int>{4, 3}.isSquare());
-    REQUIRE(DenseMatrix<int>{1, 1}.isSquare());
-    REQUIRE(DenseMatrix<int>{2, 2}.isSquare());
-    REQUIRE(DenseMatrix<int>{3, 3}.isSquare());
-  }
-
-  SECTION("write access")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-    REQUIRE(A(0, 0) == 0);
-    REQUIRE(A(0, 1) == 1);
-    REQUIRE(A(0, 2) == 2);
-    REQUIRE(A(1, 0) == 3);
-    REQUIRE(A(1, 1) == 4);
-    REQUIRE(A(1, 2) == 5);
-    DenseMatrix<const int> const_A = A;
-    REQUIRE(const_A(0, 0) == 0);
-    REQUIRE(const_A(0, 1) == 1);
-    REQUIRE(const_A(0, 2) == 2);
-    REQUIRE(const_A(1, 0) == 3);
-    REQUIRE(const_A(1, 1) == 4);
-    REQUIRE(const_A(1, 2) == 5);
-  }
-
-  SECTION("fill")
-  {
-    DenseMatrix<int> A{2, 3};
-    A.fill(2);
-    REQUIRE(A(0, 0) == 2);
-    REQUIRE(A(0, 1) == 2);
-    REQUIRE(A(0, 2) == 2);
-    REQUIRE(A(1, 0) == 2);
-    REQUIRE(A(1, 1) == 2);
-    REQUIRE(A(1, 2) == 2);
-  }
-
-  SECTION("zero matrix")
-  {
-    DenseMatrix<int> A{2, 3};
-    A = zero;
-    for (size_t i = 0; i < A.numberOfRows(); ++i) {
-      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
-        REQUIRE(A(i, j) == 0);
-      }
-    }
-  }
-
-  SECTION("identity matrix")
-  {
-    DenseMatrix<int> A{3, 3};
-    A = identity;
-    for (size_t i = 0; i < A.numberOfRows(); ++i) {
-      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
-        REQUIRE(A(i, j) == (i == j));
-      }
-    }
-  }
-
-  SECTION("copy constructor (shallow)")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    const DenseMatrix<int> B = A;
-    REQUIRE(B(0, 0) == 0);
-    REQUIRE(B(0, 1) == 1);
-    REQUIRE(B(0, 2) == 2);
-    REQUIRE(B(1, 0) == 3);
-    REQUIRE(B(1, 1) == 4);
-    REQUIRE(B(1, 2) == 5);
-  }
-
-  SECTION("copy (deep)")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    const DenseMatrix<int> B = copy(A);
-
-    A(0, 0) = 10;
-    A(0, 1) = 11;
-    A(0, 2) = 12;
-    A(1, 0) = 13;
-    A(1, 1) = 14;
-    A(1, 2) = 15;
-
-    REQUIRE(B(0, 0) == 0);
-    REQUIRE(B(0, 1) == 1);
-    REQUIRE(B(0, 2) == 2);
-    REQUIRE(B(1, 0) == 3);
-    REQUIRE(B(1, 1) == 4);
-    REQUIRE(B(1, 2) == 5);
-
-    REQUIRE(A(0, 0) == 10);
-    REQUIRE(A(0, 1) == 11);
-    REQUIRE(A(0, 2) == 12);
-    REQUIRE(A(1, 0) == 13);
-    REQUIRE(A(1, 1) == 14);
-    REQUIRE(A(1, 2) == 15);
-  }
-
-  SECTION("self scalar multiplication")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    A *= 2;
-    REQUIRE(A(0, 0) == 0);
-    REQUIRE(A(0, 1) == 2);
-    REQUIRE(A(0, 2) == 4);
-    REQUIRE(A(1, 0) == 6);
-    REQUIRE(A(1, 1) == 8);
-    REQUIRE(A(1, 2) == 10);
-  }
-
-  SECTION("left scalar multiplication")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    const DenseMatrix<int> B = 2 * A;
-
-    REQUIRE(B(0, 0) == 0);
-    REQUIRE(B(0, 1) == 2);
-    REQUIRE(B(0, 2) == 4);
-    REQUIRE(B(1, 0) == 6);
-    REQUIRE(B(1, 1) == 8);
-    REQUIRE(B(1, 2) == 10);
-  }
-
-  SECTION("product matrix vector")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 6;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    Vector<int> x{3};
-    x[0] = 7;
-    x[1] = 3;
-    x[2] = 4;
-
-    Vector y = A * x;
-    REQUIRE(y[0] == 53);
-    REQUIRE(y[1] == 53);
-  }
-
-  SECTION("self scalar division")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 2;
-    A(0, 2) = 4;
-    A(1, 0) = 6;
-    A(1, 1) = 8;
-    A(1, 2) = 10;
-
-    A /= 2;
-
-    REQUIRE(A(0, 0) == 0);
-    REQUIRE(A(0, 1) == 1);
-    REQUIRE(A(0, 2) == 2);
-    REQUIRE(A(1, 0) == 3);
-    REQUIRE(A(1, 1) == 4);
-    REQUIRE(A(1, 2) == 5);
-  }
-
-  SECTION("self minus")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    DenseMatrix<int> B{2, 3};
-    B(0, 0) = 5;
-    B(0, 1) = 6;
-    B(0, 2) = 4;
-    B(1, 0) = 2;
-    B(1, 1) = 1;
-    B(1, 2) = 3;
-
-    A -= B;
-
-    REQUIRE(A(0, 0) == -5);
-    REQUIRE(A(0, 1) == -5);
-    REQUIRE(A(0, 2) == -2);
-    REQUIRE(A(1, 0) == 1);
-    REQUIRE(A(1, 1) == 3);
-    REQUIRE(A(1, 2) == 2);
-  }
-
-  SECTION("self sum")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    DenseMatrix<int> B{2, 3};
-    B(0, 0) = 5;
-    B(0, 1) = 6;
-    B(0, 2) = 4;
-    B(1, 0) = 2;
-    B(1, 1) = 1;
-    B(1, 2) = 3;
-
-    A += B;
-
-    REQUIRE(A(0, 0) == 5);
-    REQUIRE(A(0, 1) == 7);
-    REQUIRE(A(0, 2) == 6);
-    REQUIRE(A(1, 0) == 5);
-    REQUIRE(A(1, 1) == 5);
-    REQUIRE(A(1, 2) == 8);
-  }
-
-  SECTION("sum")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 6;
-    A(0, 1) = 5;
-    A(0, 2) = 4;
-    A(1, 0) = 3;
-    A(1, 1) = 2;
-    A(1, 2) = 1;
-
-    DenseMatrix<int> B{2, 3};
-    B(0, 0) = 0;
-    B(0, 1) = 1;
-    B(0, 2) = 2;
-    B(1, 0) = 3;
-    B(1, 1) = 4;
-    B(1, 2) = 5;
-
-    DenseMatrix C = A + B;
-    REQUIRE(C(0, 0) == 6);
-    REQUIRE(C(0, 1) == 6);
-    REQUIRE(C(0, 2) == 6);
-    REQUIRE(C(1, 0) == 6);
-    REQUIRE(C(1, 1) == 6);
-    REQUIRE(C(1, 2) == 6);
-  }
-
-  SECTION("difference")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 6;
-    A(0, 1) = 5;
-    A(0, 2) = 4;
-    A(1, 0) = 3;
-    A(1, 1) = 2;
-    A(1, 2) = 1;
-
-    DenseMatrix<int> B{2, 3};
-    B(0, 0) = 0;
-    B(0, 1) = 1;
-    B(0, 2) = 2;
-    B(1, 0) = 3;
-    B(1, 1) = 4;
-    B(1, 2) = 5;
-
-    DenseMatrix C = A - B;
-    REQUIRE(C(0, 0) == 6);
-    REQUIRE(C(0, 1) == 4);
-    REQUIRE(C(0, 2) == 2);
-    REQUIRE(C(1, 0) == 0);
-    REQUIRE(C(1, 1) == -2);
-    REQUIRE(C(1, 2) == -4);
-  }
-
-  SECTION("transpose")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 0;
-    A(0, 1) = 1;
-    A(0, 2) = 2;
-    A(1, 0) = 3;
-    A(1, 1) = 4;
-    A(1, 2) = 5;
-
-    DenseMatrix B = transpose(A);
-    REQUIRE(B(0, 0) == 0);
-    REQUIRE(B(0, 1) == 3);
-    REQUIRE(B(1, 0) == 1);
-    REQUIRE(B(1, 1) == 4);
-    REQUIRE(B(2, 0) == 2);
-    REQUIRE(B(2, 1) == 5);
-  }
-
-  SECTION("product matrix vector")
-  {
-    DenseMatrix<int> A{2, 3};
-    A(0, 0) = 1;
-    A(0, 1) = 2;
-    A(0, 2) = 3;
-    A(1, 0) = 4;
-    A(1, 1) = 5;
-    A(1, 2) = 6;
-
-    DenseMatrix<int> B{3, 2};
-    B(0, 0) = 2;
-    B(0, 1) = 8;
-    B(1, 0) = 4;
-    B(1, 1) = 9;
-    B(2, 0) = 6;
-    B(2, 1) = 10;
-
-    DenseMatrix C = A * B;
-    REQUIRE(C(0, 0) == 28);
-    REQUIRE(C(0, 1) == 56);
-    REQUIRE(C(1, 0) == 64);
-    REQUIRE(C(1, 1) == 137);
-  }
-
-  SECTION("output")
-  {
-    DenseMatrix<double> A(3, 2);
-    A(0, 0) = 1;
-    A(0, 1) = 3;
-    A(1, 0) = -8;
-    A(1, 1) = -2;
-    A(2, 0) = 4;
-    A(2, 1) = -5;
-    std::ostringstream A_ost;
-    A_ost << A;
-
-    std::ostringstream ref_ost;
-    ref_ost << "0| " << 0 << ':' << 1. << ' ' << 1 << ':' << 3. << '\n';
-    ref_ost << "1| " << 0 << ':' << -8. << ' ' << 1 << ':' << -2. << '\n';
-    ref_ost << "2| " << 0 << ':' << 4. << ' ' << 1 << ':' << -5. << '\n';
-
-    REQUIRE(A_ost.str() == ref_ost.str());
-  }
-
-#ifndef NDEBUG
-  SECTION("output with signaling NaN")
-  {
-    DenseMatrix<double> A(3, 2);
-    A(0, 0) = 1;
-    A(1, 1) = -2;
-    A(2, 1) = -5;
-    std::ostringstream A_ost;
-    A_ost << A;
-
-    std::ostringstream ref_ost;
-    ref_ost << "0| " << 0 << ':' << 1. << ' ' << 1 << ":nan\n";
-    ref_ost << "1| " << 0 << ":nan " << 1 << ':' << -2. << '\n';
-    ref_ost << "2| " << 0 << ":nan " << 1 << ':' << -5. << '\n';
-
-    REQUIRE(A_ost.str() == ref_ost.str());
-  }
-
-  SECTION("non square identity matrix")
-  {
-    DenseMatrix<int> A(2, 3);
-
-    REQUIRE_THROWS_WITH(A = identity, "identity must be a square matrix");
-  }
-
-  SECTION("invalid matrix vector product")
-  {
-    DenseMatrix<int> A(2, 3);
-    Vector<int> u(2);
-
-    REQUIRE_THROWS_WITH(A * u, "cannot compute matrix-vector product: incompatible sizes");
-  }
-
-  SECTION("invalid matrix product")
-  {
-    DenseMatrix<int> A(2, 3);
-    DenseMatrix<int> B(2, 3);
-
-    REQUIRE_THROWS_WITH(A * B, "cannot compute matrix product: incompatible sizes");
-  }
-
-  SECTION("invalid matrix substract")
-  {
-    DenseMatrix<int> A(2, 3);
-    DenseMatrix<int> B(3, 2);
-
-    REQUIRE_THROWS_WITH(A -= B, "cannot substract matrix: incompatible sizes");
-  }
-
-  SECTION("invalid matrix sum")
-  {
-    DenseMatrix<int> A(2, 3);
-    DenseMatrix<int> B(3, 2);
-
-    REQUIRE_THROWS_WITH(A + B, "cannot compute matrix sum: incompatible sizes");
-  }
-
-  SECTION("invalid matrix difference")
-  {
-    DenseMatrix<int> A(2, 3);
-    DenseMatrix<int> B(3, 2);
-
-    REQUIRE_THROWS_WITH(A - B, "cannot compute matrix difference: incompatible sizes");
-  }
-
-  SECTION("invalid indices")
-  {
-    DenseMatrix<int> A(2, 3);
-
-    REQUIRE_THROWS_WITH(A(2, 0), "cannot access element: invalid indices");
-    REQUIRE_THROWS_WITH(A(0, 3), "cannot access element: invalid indices");
-  }
-
-#endif   // NDEBUG
-}
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index 3ca8f7f51..fb7bf0790 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -3,8 +3,8 @@
 
 #include <algebra/CRSMatrix.hpp>
 #include <algebra/CRSMatrixDescriptor.hpp>
-#include <algebra/DenseMatrix.hpp>
 #include <algebra/EigenvalueSolver.hpp>
+#include <algebra/SmallMatrix.hpp>
 
 #include <utils/pugs_config.hpp>
 
@@ -36,7 +36,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
 
       SECTION("eigenvalues")
       {
-        Array<double> eigenvalues;
+        SmallArray<double> eigenvalues;
 
 #ifdef PUGS_HAS_SLEPC
         EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalues);
@@ -51,8 +51,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
 
       SECTION("eigenvalues and eigenvectors")
       {
-        Array<double> eigenvalue_list;
-        std::vector<Vector<double>> eigenvector_list;
+        SmallArray<double> eigenvalue_list;
+        std::vector<SmallVector<double>> eigenvector_list;
 
 #ifdef PUGS_HAS_SLEPC
         EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, eigenvector_list);
@@ -61,9 +61,9 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
         REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
         REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
 
-        DenseMatrix<double> P{3};
-        DenseMatrix<double> PT{3};
-        DenseMatrix<double> D{3};
+        SmallMatrix<double> P{3};
+        SmallMatrix<double> PT{3};
+        SmallMatrix<double> D{3};
         D = zero;
 
         for (size_t i = 0; i < 3; ++i) {
@@ -74,7 +74,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
           D(i, i) = eigenvalue_list[i];
         }
 
-        DenseMatrix PDPT = P * D * PT;
+        SmallMatrix PDPT = P * D * PT;
         for (size_t i = 0; i < 3; ++i) {
           for (size_t j = 0; j < 3; ++j) {
             REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
@@ -88,8 +88,8 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
 
       SECTION("eigenvalues and passage matrix")
       {
-        Array<double> eigenvalue_list;
-        DenseMatrix<double> P{};
+        SmallArray<double> eigenvalue_list;
+        SmallMatrix<double> P{};
 
 #ifdef PUGS_HAS_SLEPC
         EigenvalueSolver{}.computeForSymmetricMatrix(A, eigenvalue_list, P);
@@ -98,14 +98,14 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
         REQUIRE(eigenvalue_list[1] == Catch::Approx(-1));
         REQUIRE(eigenvalue_list[2] == Catch::Approx(8));
 
-        DenseMatrix<double> D{3};
+        SmallMatrix<double> D{3};
         D = zero;
         for (size_t i = 0; i < 3; ++i) {
           D(i, i) = eigenvalue_list[i];
         }
-        DenseMatrix PT = transpose(P);
+        SmallMatrix PT = transpose(P);
 
-        DenseMatrix PDPT = P * D * PT;
+        SmallMatrix PDPT = P * D * PT;
         for (size_t i = 0; i < 3; ++i) {
           for (size_t j = 0; j < 3; ++j) {
             REQUIRE(PDPT(i, j) - S(i, j) == Catch::Approx(0).margin(1E-13));
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index e60f42df7..11195faab 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -545,7 +545,7 @@ TEST_CASE("LinearSolver", "[algebra]")
     {
       SECTION("symmetric system")
       {
-        DenseMatrix<double> A{5};
+        SmallMatrix<double> A{5};
         A = zero;
 
         A(0, 0) = 2;
@@ -566,8 +566,8 @@ TEST_CASE("LinearSolver", "[algebra]")
         A(4, 3) = -1;
         A(4, 4) = 2;
 
-        Vector<const double> x_exact = [] {
-          Vector<double> y{5};
+        SmallVector<const double> x_exact = [] {
+          SmallVector<double> y{5};
           y[0] = 1;
           y[1] = 3;
           y[2] = 2;
@@ -576,7 +576,7 @@ TEST_CASE("LinearSolver", "[algebra]")
           return y;
         }();
 
-        Vector<double> b = A * x_exact;
+        SmallVector<double> b = A * x_exact;
 
         SECTION("builtin")
         {
@@ -587,13 +587,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.method()  = LSMethod::cg;
             options.precond() = LSPrecond::none;
 
-            Vector<double> x{5};
+            SmallVector<double> x{5};
             x = zero;
 
             LinearSolver solver{options};
 
             solver.solveLocalSystem(A, x, b);
-            Vector error = x - x_exact;
+            SmallVector error = x - x_exact;
             REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
@@ -614,13 +614,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::none;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -628,13 +628,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::diagonal;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -642,13 +642,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::incomplete_cholesky;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -656,13 +656,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::amg;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
           }
@@ -674,13 +674,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.method()  = LSMethod::cholesky;
             options.precond() = LSPrecond::none;
 
-            Vector<double> x{5};
+            SmallVector<double> x{5};
             x = zero;
 
             LinearSolver solver{options};
 
             solver.solveLocalSystem(A, x, b);
-            Vector error = x - x_exact;
+            SmallVector error = x - x_exact;
             REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
@@ -700,7 +700,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
       SECTION("none symmetric system")
       {
-        DenseMatrix<double> A{5};
+        SmallMatrix<double> A{5};
         A = zero;
 
         A(0, 0) = 2;
@@ -721,8 +721,8 @@ TEST_CASE("LinearSolver", "[algebra]")
         A(4, 3) = 1;
         A(4, 4) = 3;
 
-        Vector<const double> x_exact = [] {
-          Vector<double> y{5};
+        SmallVector<const double> x_exact = [] {
+          SmallVector<double> y{5};
           y[0] = 1;
           y[1] = 3;
           y[2] = 2;
@@ -731,7 +731,7 @@ TEST_CASE("LinearSolver", "[algebra]")
           return y;
         }();
 
-        Vector<double> b = A * x_exact;
+        SmallVector<double> b = A * x_exact;
 
         SECTION("builtin")
         {
@@ -742,13 +742,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.method()  = LSMethod::bicgstab;
             options.precond() = LSPrecond::none;
 
-            Vector<double> x{5};
+            SmallVector<double> x{5};
             x = zero;
 
             LinearSolver solver{options};
 
             solver.solveLocalSystem(A, x, b);
-            Vector error = x - x_exact;
+            SmallVector error = x - x_exact;
             REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
@@ -769,13 +769,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::none;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -783,13 +783,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::diagonal;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -797,13 +797,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::incomplete_LU;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
           }
@@ -819,13 +819,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::none;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -833,13 +833,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::diagonal;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
           }
@@ -855,13 +855,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::none;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -869,13 +869,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::diagonal;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
 
@@ -883,13 +883,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             {
               options.precond() = LSPrecond::incomplete_LU;
 
-              Vector<double> x{5};
+              SmallVector<double> x{5};
               x = zero;
 
               LinearSolver solver{options};
 
               solver.solveLocalSystem(A, x, b);
-              Vector error = x - x_exact;
+              SmallVector error = x - x_exact;
               REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
             }
           }
@@ -901,13 +901,13 @@ TEST_CASE("LinearSolver", "[algebra]")
             options.method()  = LSMethod::lu;
             options.precond() = LSPrecond::none;
 
-            Vector<double> x{5};
+            SmallVector<double> x{5};
             x = zero;
 
             LinearSolver solver{options};
 
             solver.solveLocalSystem(A, x, b);
-            Vector error = x - x_exact;
+            SmallVector error = x - x_exact;
             REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
diff --git a/tests/test_PETScUtils.cpp b/tests/test_PETScUtils.cpp
index 2e0daa851..493d2dd44 100644
--- a/tests/test_PETScUtils.cpp
+++ b/tests/test_PETScUtils.cpp
@@ -29,9 +29,9 @@ TEST_CASE("PETScUtils", "[algebra]")
       }
     }
 
-    SECTION("from DenseMatrix")
+    SECTION("from SmallMatrix")
     {
-      DenseMatrix<double> A(3, 3);
+      SmallMatrix<double> A(3, 3);
       for (size_t i = 0; i < 3; ++i) {
         for (size_t j = 0; j < 3; ++j) {
           A(i, j) = 1 + i * 3 + j;
@@ -52,9 +52,9 @@ TEST_CASE("PETScUtils", "[algebra]")
       }
     }
 
-    SECTION("from DenseMatrix [non-square]")
+    SECTION("from SmallMatrix [non-square]")
     {
-      DenseMatrix<double> A(4, 3);
+      SmallMatrix<double> A(4, 3);
       for (size_t i = 0; i < 4; ++i) {
         for (size_t j = 0; j < 3; ++j) {
           A(i, j) = 1 + i * 3 + j;
diff --git a/tests/test_SmallArray.cpp b/tests/test_SmallArray.cpp
new file mode 100644
index 000000000..9b5d48279
--- /dev/null
+++ b/tests/test_SmallArray.cpp
@@ -0,0 +1,285 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/SmallArray.hpp>
+#include <utils/Types.hpp>
+
+#include <deque>
+#include <list>
+#include <set>
+#include <unordered_set>
+#include <valarray>
+#include <vector>
+
+// Instantiate to ensure full coverage is performed
+template class SmallArray<int>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SmallArray", "[utils]")
+{
+  SmallArray<int> a(10);
+  REQUIRE(a.size() == 10);
+
+  for (size_t i = 0; i < a.size(); ++i) {
+    a[i] = 2 * i;
+  }
+
+  REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+           (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+  SECTION("checking for copies")
+  {
+    SmallArray<const int> b{a};
+
+    REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
+             (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));
+
+    SmallArray<int> c{a};
+
+    REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
+             (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));
+
+    SmallArray<int> d = std::move(c);
+
+    REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
+             (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
+  }
+
+  SECTION("checking for fill")
+  {
+    SmallArray<int> b(10);
+    b.fill(3);
+
+    REQUIRE(((b[0] == 3) and (b[1] == 3) and (b[2] == 3) and (b[3] == 3) and (b[4] == 3) and (b[5] == 3) and
+             (b[6] == 3) and (b[7] == 3) and (b[8] == 3) and (b[9] == 3)));
+  }
+
+  SECTION("checking for affectations (shallow copy)")
+  {
+    SmallArray<const int> b;
+    b = a;
+
+    REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
+             (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));
+
+    SmallArray<int> c;
+    c = a;
+
+    REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
+             (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));
+
+    SmallArray<int> d;
+    d = std::move(c);
+
+    REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
+             (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
+  }
+
+  SECTION("checking for affectations (deep copy)")
+  {
+    SmallArray<int> b(copy(a));
+
+    REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
+             (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));
+
+    b.fill(2);
+
+    REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+             (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+    REQUIRE(((b[0] == 2) and (b[1] == 2) and (b[2] == 2) and (b[3] == 2) and (b[4] == 2) and (b[5] == 2) and
+             (b[6] == 2) and (b[7] == 2) and (b[8] == 2) and (b[9] == 2)));
+
+    SmallArray<int> c;
+    c = a;
+
+    REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
+             (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));
+
+    c = copy(b);
+
+    REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+             (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+    REQUIRE(((c[0] == 2) and (c[1] == 2) and (c[2] == 2) and (c[3] == 2) and (c[4] == 2) and (c[5] == 2) and
+             (c[6] == 2) and (c[7] == 2) and (c[8] == 2) and (c[9] == 2)));
+
+    SmallArray<int> d{a.size()};
+    copy_to(a, d);
+
+    REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
+             (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
+
+    REQUIRE(((c[0] == 2) and (c[1] == 2) and (c[2] == 2) and (c[3] == 2) and (c[4] == 2) and (c[5] == 2) and
+             (c[6] == 2) and (c[7] == 2) and (c[8] == 2) and (c[9] == 2)));
+
+    copy_to(c, d);
+
+    REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
+             (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));
+
+    REQUIRE(((d[0] == 2) and (d[1] == 2) and (d[2] == 2) and (d[3] == 2) and (d[4] == 2) and (d[5] == 2) and
+             (d[6] == 2) and (d[7] == 2) and (d[8] == 2) and (d[9] == 2)));
+  }
+
+  SECTION("checking for std container conversion")
+  {
+    {
+      std::vector<int> v{1, 2, 5, 3};
+      {
+        SmallArray<int> v_array = convert_to_small_array(v);
+
+        REQUIRE(v_array.size() == v.size());
+        REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
+      }
+
+      {
+        SmallArray<const int> v_array = convert_to_small_array(v);
+
+        REQUIRE(v_array.size() == v.size());
+        REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
+      }
+    }
+
+    {
+      std::vector<int> w;
+      {
+        SmallArray<int> w_array = convert_to_small_array(w);
+        REQUIRE(w_array.size() == 0);
+      }
+      {
+        SmallArray<const int> w_array = convert_to_small_array(w);
+        REQUIRE(w_array.size() == 0);
+      }
+    }
+
+    {
+      std::valarray<int> v{1, 2, 5, 3};
+      SmallArray<int> v_array = convert_to_small_array(v);
+
+      REQUIRE(v_array.size() == v.size());
+      REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
+    }
+
+    {
+      std::set<int> s{4, 2, 5, 3, 1, 3, 2};
+      SmallArray<int> s_array = convert_to_small_array(s);
+
+      REQUIRE(s_array.size() == s.size());
+      REQUIRE(
+        ((s_array[0] == 1) and (s_array[1] == 2) and (s_array[2] == 3) and (s_array[3] == 4) and (s_array[4] == 5)));
+    }
+
+    {
+      std::unordered_set<int> us{4, 2, 5, 3, 1, 3, 2};
+      SmallArray<int> us_array = convert_to_small_array(us);
+
+      REQUIRE(us_array.size() == us.size());
+
+      std::set<int> s;
+      for (size_t i = 0; i < us_array.size(); ++i) {
+        REQUIRE((us.find(us_array[i]) != us.end()));
+        s.insert(us_array[i]);
+      }
+      REQUIRE(s.size() == us_array.size());
+    }
+
+    {
+      std::multiset<int> ms{4, 2, 5, 3, 1, 3, 2};
+      SmallArray<int> ms_array = convert_to_small_array(ms);
+
+      REQUIRE(ms_array.size() == ms.size());
+      REQUIRE(((ms_array[0] == 1) and (ms_array[1] == 2) and (ms_array[2] == 2) and (ms_array[3] == 3) and
+               (ms_array[4] == 3) and (ms_array[5] == 4) and (ms_array[6] == 5)));
+    }
+
+    {
+      std::list<int> l{1, 3, 5, 6, 2};
+      SmallArray<int> l_array = convert_to_small_array(l);
+
+      REQUIRE(l_array.size() == l.size());
+      REQUIRE(
+        ((l_array[0] == 1) and (l_array[1] == 3) and (l_array[2] == 5) and (l_array[3] == 6) and (l_array[4] == 2)));
+    }
+
+    {
+      std::deque<int> q{1, 3, 5, 6, 2};
+      q.push_front(2);
+      SmallArray<int> q_array = convert_to_small_array(q);
+
+      REQUIRE(q_array.size() == q.size());
+      REQUIRE(((q_array[0] == 2) and (q_array[1] == 1) and (q_array[2] == 3) and (q_array[3] == 5) and
+               (q_array[4] == 6) and (q_array[5] == 2)));
+    }
+  }
+
+  SECTION("output")
+  {
+    SmallArray<int> x{5};
+    x[0] = 2;
+    x[1] = 6;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 7;
+
+    std::ostringstream array_ost;
+    array_ost << x;
+    std::ostringstream ref_ost;
+    ref_ost << 0 << ':' << x[0];
+    for (size_t i = 1; i < x.size(); ++i) {
+      ref_ost << ' ' << i << ':' << x[i];
+    }
+    REQUIRE(array_ost.str() == ref_ost.str());
+  }
+
+#ifndef NDEBUG
+
+  SECTION("output with signaling NaN")
+  {
+    SmallArray<double> x{5};
+    x[0] = 2;
+    x[2] = 3;
+
+    std::ostringstream array_ost;
+    array_ost << x;
+    std::ostringstream ref_ost;
+    ref_ost << 0 << ':' << 2 << ' ';
+    ref_ost << 1 << ":nan ";
+    ref_ost << 2 << ':' << 3 << ' ';
+    ref_ost << 3 << ":nan ";
+    ref_ost << 4 << ":nan";
+    REQUIRE(array_ost.str() == ref_ost.str());
+  }
+
+  SECTION("checking for bounds violation")
+  {
+    REQUIRE_THROWS_AS(a[10], AssertError);
+  }
+
+  SECTION("invalid copy_to")
+  {
+    SmallArray<int> b{2 * a.size()};
+    REQUIRE_THROWS_AS(copy_to(a, b), AssertError);
+  }
+
+  SECTION("checking for nan initialization")
+  {
+    SmallArray<double> array(10);
+
+    for (size_t i = 0; i < array.size(); ++i) {
+      REQUIRE(std::isnan(array[i]));
+    }
+  }
+
+  SECTION("checking for bad initialization")
+  {
+    SmallArray<int> array(10);
+
+    for (size_t i = 0; i < array.size(); ++i) {
+      REQUIRE(array[i] == std::numeric_limits<int>::max() / 2);
+    }
+  }
+#endif   // NDEBUG
+}
-- 
GitLab