diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index 2d03a9cf5e7cced7f99b6707c042395f48054199..27d1322fa25a35ebcf7845911280006dfc4af8b7 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -111,23 +111,21 @@ class LinearSolver
   void
   solveLocalSystem(const TinyMatrix<N>& A, TinyVector<N>& x, const TinyVector<N>& b)
   {
-    Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
+    Assert((x.dimension() == b.dimension()) and (x.dimension() - A.numberOfColumns() == 0) and A.isSquare());
     switch (m_options.library()) {
     case LSLibrary::builtin: {
       // Probably expensive but builtin is not the preferred way to
       // solve linear systems
-      SmallMatrix A_dense{A};
+      SmallMatrix<double> A_dense(A);
 
-      SmallVector x_vector{x};
-      SmallVector<const double> b_vector{b};
+      SmallVector<double> x_vector(x);
+      SmallVector<double> b_vector(b);
 
-      this->solveLocalSystem(A_dense, x_vector, b_vector);
+      _builtinSolveLocalSystem(A_dense, x_vector, b_vector);
 
       for (size_t i = 0; i < N; ++i) {
         x[i] = x_vector[i];
       }
-
-      _builtinSolveLocalSystem(A, x, b);
       break;
     }
       // LCOV_EXCL_START
diff --git a/src/algebra/SmallMatrix.hpp b/src/algebra/SmallMatrix.hpp
index 03d69f329dbfc9592d5eed02e6fad66885df7de6..718000c2799067e0084d5bb3460acd688a8642ea 100644
--- a/src/algebra/SmallMatrix.hpp
+++ b/src/algebra/SmallMatrix.hpp
@@ -115,8 +115,7 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
   PUGS_INLINE SmallMatrix&
   operator*=(const DataType2& a)
   {
-    parallel_for(
-      m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] *= a; });
+    parallel_for(m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] *= a; });
     return *this;
   }
 
@@ -129,8 +128,7 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
     Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
            "cannot substract matrix: incompatible sizes");
 
-    parallel_for(
-      m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] -= B.m_values[i]; });
+    parallel_for(m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] -= B.m_values[i]; });
     return *this;
   }
 
@@ -143,8 +141,7 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
     Assert((m_nb_rows == B.numberOfRows()) and (m_nb_columns == B.numberOfColumns()),
            "cannot add matrix: incompatible sizes");
 
-    parallel_for(
-      m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] += B.m_values[i]; });
+    parallel_for(m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] += B.m_values[i]; });
     return *this;
   }
 
@@ -159,8 +156,7 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
 
     SmallMatrix<std::remove_const_t<DataType>> sum{B.numberOfRows(), B.numberOfColumns()};
 
-    parallel_for(
-      m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + B.m_values[i]; });
+    parallel_for(m_values.size(), PUGS_CLASS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + B.m_values[i]; });
 
     return sum;
   }
@@ -223,8 +219,7 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
     Assert(m_nb_rows == m_nb_columns, "identity must be a square matrix");
 
     m_values.fill(0);
-    parallel_for(
-      m_nb_rows, PUGS_CLASS_LAMBDA(const index_type i) { m_values[i * m_nb_rows + i] = 1; });
+    parallel_for(m_nb_rows, PUGS_CLASS_LAMBDA(const index_type i) { m_values[i * m_nb_rows + i] = 1; });
     return *this;
   }
 
@@ -295,7 +290,8 @@ class [[nodiscard]] SmallMatrix   // LCOV_EXCL_LINE
   {}
 
   template <size_t M, size_t N>
-  explicit SmallMatrix(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, std::remove_const_t<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) {
diff --git a/src/algebra/SmallVector.hpp b/src/algebra/SmallVector.hpp
index 296b8be71ce8b20346174ce40f41edeb78f1b1c2..91093bc67203f29cad05c475e88203b99b743b82 100644
--- a/src/algebra/SmallVector.hpp
+++ b/src/algebra/SmallVector.hpp
@@ -57,8 +57,7 @@ class SmallVector   // LCOV_EXCL_LINE
   operator-() const
   {
     SmallVector<std::remove_const_t<DataType>> opposite(this->size());
-    parallel_for(
-      opposite.size(), PUGS_CLASS_LAMBDA(index_type i) { opposite.m_values[i] = -m_values[i]; });
+    parallel_for(opposite.size(), PUGS_CLASS_LAMBDA(index_type i) { opposite.m_values[i] = -m_values[i]; });
     return opposite;
   }
 
@@ -97,8 +96,7 @@ class SmallVector   // LCOV_EXCL_LINE
   PUGS_INLINE SmallVector&
   operator*=(const DataType2& a)
   {
-    parallel_for(
-      this->size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] *= a; });
+    parallel_for(this->size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] *= a; });
     return *this;
   }
 
@@ -108,8 +106,7 @@ class SmallVector   // LCOV_EXCL_LINE
   {
     Assert(this->size() == y.size(), "cannot substract vector: incompatible sizes");
 
-    parallel_for(
-      this->size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] -= y[i]; });
+    parallel_for(this->size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] -= y[i]; });
 
     return *this;
   }
@@ -120,8 +117,7 @@ class SmallVector   // LCOV_EXCL_LINE
   {
     Assert(this->size() == y.size(), "cannot add vector: incompatible sizes");
 
-    parallel_for(
-      this->size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] += y[i]; });
+    parallel_for(this->size(), PUGS_CLASS_LAMBDA(index_type i) { m_values[i] += y[i]; });
 
     return *this;
   }
@@ -133,8 +129,7 @@ class SmallVector   // LCOV_EXCL_LINE
     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_CLASS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + y[i]; });
+    parallel_for(this->size(), PUGS_CLASS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] + y[i]; });
 
     return sum;
   }
@@ -146,8 +141,7 @@ class SmallVector   // LCOV_EXCL_LINE
     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_CLASS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] - y[i]; });
+    parallel_for(this->size(), PUGS_CLASS_LAMBDA(index_type i) { sum.m_values[i] = m_values[i] - y[i]; });
 
     return sum;
   }
@@ -217,7 +211,7 @@ class SmallVector   // LCOV_EXCL_LINE
   }
 
   template <size_t N>
-  explicit SmallVector(const TinyVector<N, DataType>& x) : m_values{N}
+  explicit SmallVector(const TinyVector<N, std::remove_const_t<DataType>>& x) : m_values{N}
   {
     std::copy(&x[0], &x[0] + N, &m_values[0]);
   }
diff --git a/src/algebra/VectorWrapper.hpp b/src/algebra/VectorWrapper.hpp
index 0bdeefbfd3f10407358613e7f1d2487bba22b26d..ac20be67d4a45b4253585b892bb2fcd6f9e56fa6 100644
--- a/src/algebra/VectorWrapper.hpp
+++ b/src/algebra/VectorWrapper.hpp
@@ -33,7 +33,8 @@ class VectorWrapper
   VectorWrapper(const SmallVector<DataType>& x) : m_size{x.size()}, m_vector_ptr{(m_size > 0) ? (&x[0]) : nullptr} {}
 
   template <size_t N>
-  VectorWrapper(const TinyVector<N, DataType>& x) : m_size{N}, m_vector_ptr{(N > 0) ? (&x[0]) : nullptr}
+  VectorWrapper(const TinyVector<N, std::remove_const_t<DataType>>& x)
+    : m_size{N}, m_vector_ptr{(N > 0) ? (&x[0]) : nullptr}
   {}
 
   template <size_t N>
diff --git a/src/utils/SmallArray.hpp b/src/utils/SmallArray.hpp
index 81d70a848187216e82cf80c5208f96cb64e05dfa..df6fae320f12ee258b04da3803d686d773989b9a 100644
--- a/src/utils/SmallArray.hpp
+++ b/src/utils/SmallArray.hpp
@@ -85,7 +85,8 @@ class [[nodiscard]] SmallArray
   SmallArray& operator=(SmallArray&&) = default;
 
   PUGS_INLINE
-  explicit SmallArray(size_t size) : m_size{size}, m_values{std::shared_ptr<DataType[]>(new DataType[size])}
+  explicit SmallArray(size_t size)
+    : m_size{size}, m_values{std::shared_ptr<std::remove_const_t<DataType>[]>(new DataType[size])}
   {
     static_assert(not std::is_const<DataType>(), "Cannot allocate SmallArray of const data: only view is "
                                                  "supported");
diff --git a/tests/test_SmallMatrix.cpp b/tests/test_SmallMatrix.cpp
index 131b1e1f2fc4ae64622d410ac8e38ea051e9c9a7..3d5349c9c97cdb71aeaaf9f644af691ce42a7a33 100644
--- a/tests/test_SmallMatrix.cpp
+++ b/tests/test_SmallMatrix.cpp
@@ -113,7 +113,7 @@ TEST_CASE("SmallMatrix", "[algebra]")
   SECTION("build from TinyMatrix")
   {
     TinyMatrix<3, 4> T{1, 2, 3, 4, 0.5, 1, 1.5, 2, 0.3, 0.6, 0.9, 1.2};
-    SmallMatrix S{T};
+    SmallMatrix<double> S{T};
 
     REQUIRE(S.numberOfRows() == T.numberOfRows());
     REQUIRE(S.numberOfColumns() == T.numberOfColumns());
diff --git a/tests/test_SmallVector.cpp b/tests/test_SmallVector.cpp
index b5ffc6819736b2e55d33706aeb18f1e1d28dabca..05defc3197af10a54d746afb892b6a57dcc83960 100644
--- a/tests/test_SmallVector.cpp
+++ b/tests/test_SmallVector.cpp
@@ -437,7 +437,7 @@ TEST_CASE("SmallVector", "[algebra]")
   {
     TinyVector<5> tiny_vector{1, 3, 5, 7, 9};
 
-    SmallVector vector{tiny_vector};
+    SmallVector<double> vector{tiny_vector};
 
     REQUIRE(vector[0] == 1);
     REQUIRE(vector[1] == 3);