diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index edb85d5e87b6a347ad5bc4f82987fd9255154484..bf650afe2784035023fbb5311c1c3a12c68b2a3b 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -1,28 +1,31 @@
 #ifndef CRS_MATRIX_HPP
 #define CRS_MATRIX_HPP
 
-#include <Kokkos_StaticCrsGraph.hpp>
-
 #include <utils/Array.hpp>
 #include <utils/PugsAssert.hpp>
 
-#include <algebra/SparseMatrixDescriptor.hpp>
 #include <algebra/Vector.hpp>
 
 #include <iostream>
 
-template <typename DataType = double, typename IndexType = size_t>
+template <typename DataType, typename IndexType>
+class CRSMatrixDescriptor;
+
+template <typename DataType = double, typename IndexType = int>
 class CRSMatrix
 {
+ public:
   using MutableDataType = std::remove_const_t<DataType>;
+  using index_type      = IndexType;
+  using data_type       = DataType;
 
  private:
-  using HostMatrix = Kokkos::StaticCrsGraph<const IndexType, Kokkos::HostSpace>;
+  const IndexType m_nb_rows;
+  const IndexType m_nb_columns;
 
-  HostMatrix m_host_matrix;
+  Array<const IndexType> m_row_map;
   Array<const DataType> m_values;
-  size_t m_nb_rows;
-  size_t m_nb_columns;
+  Array<const IndexType> m_column_indices;
 
  public:
   PUGS_INLINE
@@ -33,14 +36,14 @@ class CRSMatrix
   }
 
   PUGS_INLINE
-  size_t
+  IndexType
   numberOfRows() const
   {
     return m_nb_rows;
   }
 
   PUGS_INLINE
-  size_t
+  IndexType
   numberOfColumns() const
   {
     return m_nb_columns;
@@ -53,15 +56,15 @@ class CRSMatrix
   }
 
   auto
-  rowIndices() const
+  rowMap() const
   {
-    return encapsulate(m_host_matrix.row_map);
+    return m_row_map;
   }
 
   auto
-  row(size_t i) const
+  columnIndices() const
   {
-    return m_host_matrix.rowConst(i);
+    return m_column_indices;
   }
 
   template <typename DataType2>
@@ -71,18 +74,17 @@ class CRSMatrix
     static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
                   "Cannot multiply matrix and vector of different type");
 
-    Assert(x.size() == m_host_matrix.numRows());
+    Assert(x.size() - m_nb_columns == 0);
 
-    Vector<MutableDataType> Ax{m_host_matrix.numRows()};
-    auto host_row_map = m_host_matrix.row_map;
+    Vector<MutableDataType> Ax(m_nb_rows);
 
     parallel_for(
-      m_host_matrix.numRows(), PUGS_LAMBDA(const IndexType& i_row) {
-        const auto row_begin = host_row_map(i_row);
-        const auto row_end   = host_row_map(i_row + 1);
+      m_nb_rows, PUGS_LAMBDA(const IndexType& i_row) {
+        const auto row_begin = m_row_map[i_row];
+        const auto row_end   = m_row_map[i_row + 1];
         MutableDataType sum{0};
         for (IndexType j = row_begin; j < row_end; ++j) {
-          sum += m_values[j] * x[m_host_matrix.entries(j)];
+          sum += m_values[j] * x[m_column_indices[j]];
         }
         Ax[i_row] = sum;
       });
@@ -90,21 +92,37 @@ class CRSMatrix
     return Ax;
   }
 
-  CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
-    : m_nb_rows{M.numberOfRows()}, m_nb_columns{M.numberOfColumns()}
+  friend PUGS_INLINE std::ostream&
+  operator<<(std::ostream& os, const CRSMatrix& A)
   {
-    {
-      auto host_matrix =
-        Kokkos::create_staticcrsgraph<Kokkos::StaticCrsGraph<IndexType, Kokkos::HostSpace>>("connectivity_matrix",
-                                                                                            M.graphVector());
-
-      // This is a bit crappy but it is the price to pay to avoid
-      m_host_matrix.entries           = host_matrix.entries;
-      m_host_matrix.row_map           = host_matrix.row_map;
-      m_host_matrix.row_block_offsets = host_matrix.row_block_offsets;
+    for (IndexType i = 0; i < A.m_nb_rows; ++i) {
+      const auto row_begin = A.m_row_map[i];
+      const auto row_end   = A.m_row_map[i + 1];
+      os << i << "|";
+      for (IndexType j = row_begin; j < row_end; ++j) {
+        os << ' ' << A.m_column_indices[j] << ':' << A.m_values[j];
+      }
+      os << '\n';
     }
-    m_values = M.valueArray();
+    return os;
   }
+
+ private:
+  friend class CRSMatrixDescriptor<DataType, IndexType>;
+
+  CRSMatrix(const IndexType nb_rows,
+            const IndexType nb_columns,
+            const Array<const IndexType>& row_map,
+            const Array<const DataType>& values,
+            const Array<const IndexType>& column_indices)
+    : m_nb_rows{nb_rows},
+      m_nb_columns{nb_columns},
+      m_row_map{row_map},
+      m_values{values},
+      m_column_indices{column_indices}
+  {}
+
+ public:
   ~CRSMatrix() = default;
 };
 
diff --git a/src/algebra/CRSMatrixDescriptor.hpp b/src/algebra/CRSMatrixDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cbf82a904e85dcb8aa735e978612a4baec0affa7
--- /dev/null
+++ b/src/algebra/CRSMatrixDescriptor.hpp
@@ -0,0 +1,265 @@
+#ifndef CRS_MATRIX_DESCRIPTOR_HPP
+#define CRS_MATRIX_DESCRIPTOR_HPP
+
+#include <algebra/CRSMatrix.hpp>
+#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <map>
+
+template <typename DataType = double, typename IndexType = int>
+class CRSMatrixDescriptor
+{
+ public:
+  using index_type = IndexType;
+  using data_type  = DataType;
+
+  static_assert(std::is_integral_v<IndexType>);
+
+ private:
+  const IndexType m_nb_rows;
+  const IndexType m_nb_columns;
+
+  Array<const IndexType> m_row_map;
+  Array<data_type> m_values;
+  Array<IndexType> m_column_indices;
+
+  Array<IndexType> m_nb_values_per_row;
+  Array<std::map<IndexType, DataType>> m_overflow_per_row;
+
+  Array<const IndexType>
+  _computeRowMap(const Array<IndexType>& non_zeros) const
+  {
+    Assert(non_zeros.size() - m_nb_rows == 0);
+
+    Array<IndexType> row_map(m_nb_rows + 1);
+    row_map[0] = 0;
+    for (IndexType i = 0; i < m_nb_rows; ++i) {
+      row_map[i + 1] = row_map[i] + non_zeros[i];
+    }
+
+    return row_map;
+  }
+
+ public:
+  IndexType
+  numberOfRows() const
+  {
+    return m_nb_rows;
+  }
+
+  IndexType
+  numberOfColumns() const
+  {
+    return m_nb_columns;
+  }
+
+  bool PUGS_INLINE
+  hasOverflow() const
+  {
+    for (IndexType i = 0; i < m_nb_rows; ++i) {
+      if (m_overflow_per_row[i].size() > 0) {
+        return true;
+      }
+    }
+    return false;
+  }
+
+  bool PUGS_INLINE
+  isFilled() const
+  {
+    for (IndexType i = 0; i < m_nb_rows; ++i) {
+      if (m_nb_values_per_row[i] < m_row_map[i + 1] - m_row_map[i]) {
+        return false;
+      }
+    }
+    return true;
+  }
+
+  friend PUGS_INLINE std::ostream&
+  operator<<(std::ostream& os, const CRSMatrixDescriptor& A)
+  {
+    for (IndexType i = 0; i < A.m_nb_rows; ++i) {
+      const auto& overflow_row = A.m_overflow_per_row[i];
+      os << i << "|";
+
+      if (A.m_nb_values_per_row[i] > 0) {
+        IndexType j     = 0;
+        auto i_overflow = overflow_row.begin();
+        while (j < A.m_nb_values_per_row[i] or i_overflow != overflow_row.end()) {
+          const IndexType j_index = [&] {
+            if (j < A.m_nb_values_per_row[i]) {
+              return A.m_column_indices[A.m_row_map[i] + j];
+            } else {
+              return std::numeric_limits<IndexType>::max();
+            }
+          }();
+
+          const IndexType overflow_index = [&] {
+            if (i_overflow != overflow_row.end()) {
+              return i_overflow->first;
+            } else {
+              return std::numeric_limits<IndexType>::max();
+            }
+          }();
+
+          if (j_index < overflow_index) {
+            os << ' ' << A.m_column_indices[A.m_row_map[i] + j] << ':' << A.m_values[A.m_row_map[i] + j];
+            ++j;
+          } else {
+            os << ' ' << i_overflow->first << ':' << i_overflow->second;
+            ++i_overflow;
+          }
+        }
+      }
+
+      os << '\n';
+    }
+    return os;
+  }
+
+  [[nodiscard]] PUGS_INLINE DataType&
+  operator()(const IndexType i, const IndexType j)
+  {
+    Assert(i < m_nb_rows, "invalid row index");
+    Assert(j < m_nb_columns, "invalid column index");
+
+    const IndexType row_start = m_row_map[i];
+    const IndexType row_end   = m_row_map[i + 1];
+
+    IndexType position_candidate = 0;
+    bool found                   = false;
+    for (IndexType i_row = 0; i_row < m_nb_values_per_row[i]; ++i_row) {
+      if (m_column_indices[row_start + i_row] > j) {
+        position_candidate = i_row;
+        break;
+      } else {
+        if (m_column_indices[row_start + i_row] == j) {
+          position_candidate = i_row;
+          found              = true;
+          break;
+        } else {
+          ++position_candidate;
+        }
+      }
+    }
+
+    if (not found) {
+      if (m_nb_values_per_row[i] < row_end - row_start) {
+        for (IndexType i_shift = 0; i_shift < m_nb_values_per_row[i] - position_candidate; ++i_shift) {
+          Assert(std::make_signed_t<IndexType>(row_start + m_nb_values_per_row[i]) -
+                   std::make_signed_t<IndexType>(i_shift + 1) >=
+                 0);
+          const IndexType i_destination = row_start + m_nb_values_per_row[i] - i_shift;
+          const IndexType i_source      = i_destination - 1;
+
+          m_column_indices[i_destination] = m_column_indices[i_source];
+          m_values[i_destination]         = m_values[i_source];
+        }
+
+        m_column_indices[row_start + position_candidate] = j;
+        m_values[row_start + position_candidate]         = 0;
+        ++m_nb_values_per_row[i];
+
+        return m_values[row_start + position_candidate];
+      } else {
+        auto& overflow_row = m_overflow_per_row[i];
+
+        auto iterator = overflow_row.insert(std::make_pair(j, DataType{0})).first;
+
+        return iterator->second;
+      }
+    } else {
+      return m_values[row_start + position_candidate];
+    }
+  }
+
+  CRSMatrix<DataType, IndexType>
+  getCRSMatrix() const
+  {
+    const bool is_filled    = this->isFilled();
+    const bool has_overflow = this->hasOverflow();
+    if (is_filled and not has_overflow) {
+      return CRSMatrix<DataType, IndexType>{m_nb_rows, m_nb_columns, m_row_map, m_values, m_column_indices};
+    } else {
+      std::cout << rang::fgB::yellow << "warning:" << rang::style::reset
+                << " CRSMatrix profile was not properly set:\n";
+      if (not is_filled) {
+        std::cout << "- some rows are " << rang::style::bold << "too long" << rang::style::reset << '\n';
+      }
+      if (has_overflow) {
+        std::cout << "- some rows are " << rang::style::bold << "too short" << rang::style::reset << '\n';
+      }
+
+      Array<IndexType> row_map(m_nb_rows + 1);
+
+      row_map[0] = 0;
+      for (IndexType i = 0; i < m_nb_rows; ++i) {
+        row_map[i + 1] = row_map[i] + m_nb_values_per_row[i] + m_overflow_per_row[i].size();
+      }
+
+      IndexType nb_values = row_map[row_map.size() - 1];
+
+      Array<data_type> values(nb_values);
+      Array<IndexType> column_indices(nb_values);
+
+      IndexType l = 0;
+      for (IndexType i = 0; i < m_nb_rows; ++i) {
+        const auto& overflow_row = m_overflow_per_row[i];
+
+        if (m_nb_values_per_row[i] > 0) {
+          IndexType j     = 0;
+          auto i_overflow = overflow_row.begin();
+          while (j < m_nb_values_per_row[i] or i_overflow != overflow_row.end()) {
+            const IndexType j_index = [&] {
+              if (j < m_nb_values_per_row[i]) {
+                return m_column_indices[m_row_map[i] + j];
+              } else {
+                return std::numeric_limits<IndexType>::max();
+              }
+            }();
+
+            const IndexType overflow_index = [&] {
+              if (i_overflow != overflow_row.end()) {
+                return i_overflow->first;
+              } else {
+                return std::numeric_limits<IndexType>::max();
+              }
+            }();
+
+            if (j_index < overflow_index) {
+              column_indices[l] = m_column_indices[m_row_map[i] + j];
+              values[l]         = m_values[m_row_map[i] + j];
+              ++j;
+
+            } else {
+              column_indices[l] = i_overflow->first;
+              values[l]         = i_overflow->second;
+              ++i_overflow;
+            }
+            ++l;
+          }
+        }
+      }
+
+      return CRSMatrix<DataType, IndexType>{m_nb_rows, m_nb_columns, row_map, values, column_indices};
+    }
+  }
+
+  CRSMatrixDescriptor(const IndexType nb_rows, const IndexType nb_columns, const Array<IndexType>& non_zeros)
+    : m_nb_rows{nb_rows},
+      m_nb_columns{nb_columns},
+      m_row_map{this->_computeRowMap(non_zeros)},
+      m_values(m_row_map[m_row_map.size() - 1]),
+      m_column_indices(m_row_map[m_row_map.size() - 1]),
+      m_nb_values_per_row(m_nb_rows),
+      m_overflow_per_row(m_nb_rows)
+  {
+    m_nb_values_per_row.fill(0);
+    m_values.fill(0);
+  }
+
+  ~CRSMatrixDescriptor() = default;
+};
+
+#endif   // CRS_MATRIX_DESCRIPTOR_HPP
diff --git a/src/algebra/LinearSolver.cpp b/src/algebra/LinearSolver.cpp
index 98e2a3a1c2caf8a041b6e69e71502bb6441647e7..3cffca7298b0a7377996c130195ac9e43a66ca44 100644
--- a/src/algebra/LinearSolver.cpp
+++ b/src/algebra/LinearSolver.cpp
@@ -175,7 +175,7 @@ struct LinearSolver::Internals
                         const RHSVectorType& b,
                         const LinearSolverOptions& options)
   {
-    Assert(x.size() == b.size() and x.size() == A.numberOfColumns() and A.isSquare());
+    Assert((x.size() == b.size()) and (x.size() - A.numberOfColumns() == 0) and A.isSquare());
 
     Vec petscB;
     VecCreateMPIWithArray(PETSC_COMM_SELF, 1, b.size(), b.size(), &b[0], &petscB);
@@ -329,7 +329,7 @@ LinearSolver::checkOptions(const LinearSolverOptions& options) const
 }
 
 void
-LinearSolver::solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b)
+LinearSolver::solveLocalSystem(const CRSMatrix<double, int>& A, Vector<double>& x, const Vector<const double>& b)
 {
   Internals::solveLocalSystem(A, x, b, m_options);
 }
diff --git a/src/algebra/LinearSolver.hpp b/src/algebra/LinearSolver.hpp
index 3873374002a146c79a447a815d4ef6a57669883e..1ab43a1405ee312fdb9255c6b1f4508ade4e168e 100644
--- a/src/algebra/LinearSolver.hpp
+++ b/src/algebra/LinearSolver.hpp
@@ -39,7 +39,7 @@ class LinearSolver
     }
   }
 
-  void solveLocalSystem(const CRSMatrix<double, size_t>& A, Vector<double>& x, const Vector<const double>& b);
+  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);
 
   LinearSolver();
diff --git a/src/algebra/PETScUtils.cpp b/src/algebra/PETScUtils.cpp
index 5ed261353d93673267c3f4107c35f855bf643211..63c5a68ac5f3a248ba7e5bc51bf381f5418dbe4a 100644
--- a/src/algebra/PETScUtils.cpp
+++ b/src/algebra/PETScUtils.cpp
@@ -35,32 +35,22 @@ PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const size_t nb_rows, const size_
   MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
 }
 
-PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, size_t>& A)
-  : m_nb_rows{A.numberOfRows()}, m_nb_columns{A.numberOfColumns()}
+PETScAijMatrixEmbedder::PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A)
+  : m_nb_rows(A.numberOfRows()), m_nb_columns(A.numberOfColumns())
 {
-  const Array<PetscReal>& values = copy(A.values());
+  static_assert(std::is_same_v<PetscInt, typename std::decay_t<decltype(A)>::index_type>,
+                "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());
 
-  {
-    const auto& A_row_indices = A.rowIndices();
-    Array<PetscInt> row_indices{A_row_indices.size()};
-    for (size_t i = 0; i < row_indices.size(); ++i) {
-      row_indices[i] = A_row_indices[i];
-    }
-    m_row_indices = row_indices;
-
-    Array<PetscInt> column_indices{values.size()};
-    size_t l = 0;
-    for (size_t i = 0; i < A.numberOfRows(); ++i) {
-      const auto row_i = A.row(i);
-      for (size_t j = 0; j < row_i.length; ++j) {
-        column_indices[l++] = row_i.colidx(j);
-      }
-    }
-    m_column_indices = column_indices;
-  }
-
-  MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.numberOfRows(), A.numberOfColumns(), &m_row_indices[0],
-                            &m_column_indices[0], &values[0], &m_petscMat);
+  MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, A.numberOfRows(), A.numberOfColumns(),
+                            const_cast<PetscInt*>(&m_row_indices[0]),      // These ugly const_cast
+                            const_cast<PetscInt*>(&m_column_indices[0]),   // permit CRSMatrix to
+                            const_cast<PetscScalar*>(&m_values[0]),        // be non-mutable
+                            &m_petscMat);
 
   MatAssemblyBegin(m_petscMat, MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(m_petscMat, MAT_FINAL_ASSEMBLY);
diff --git a/src/algebra/PETScUtils.hpp b/src/algebra/PETScUtils.hpp
index 68748c02f5b36092cba6dce8ba3281a6d1cbb93b..3e7cff5b87b843f84575ed9b2af4a425e6ed4f5a 100644
--- a/src/algebra/PETScUtils.hpp
+++ b/src/algebra/PETScUtils.hpp
@@ -15,8 +15,9 @@ class PETScAijMatrixEmbedder
 {
  private:
   Mat m_petscMat;
-  Array<PetscInt> m_row_indices;
-  Array<PetscInt> m_column_indices;
+  Array<const PetscScalar> m_values;
+  Array<const PetscInt> m_row_indices;
+  Array<const PetscInt> m_column_indices;
   const size_t m_nb_rows;
   const size_t m_nb_columns;
 
@@ -57,7 +58,7 @@ class PETScAijMatrixEmbedder
     : PETScAijMatrixEmbedder{A.numberOfRows(), A.numberOfColumns(), &A(0, 0)}
   {}
 
-  PETScAijMatrixEmbedder(const CRSMatrix<double, size_t>& A);
+  PETScAijMatrixEmbedder(const CRSMatrix<double, int>& A);
 
   ~PETScAijMatrixEmbedder();
 };
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
deleted file mode 100644
index d94ed468606a10ee15857a971abab20fb4be5b1a..0000000000000000000000000000000000000000
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ /dev/null
@@ -1,170 +0,0 @@
-#ifndef SPARSE_MATRIX_DESCRIPTOR_HPP
-#define SPARSE_MATRIX_DESCRIPTOR_HPP
-
-#include <utils/Array.hpp>
-
-#include <map>
-#include <type_traits>
-
-template <typename DataType = double, typename IndexType = size_t>
-class SparseMatrixDescriptor
-{
-  static_assert(std::is_integral_v<IndexType>, "Index type must be an integral type");
-  static_assert(std::is_unsigned_v<IndexType>, "Index type must be unsigned");
-
- public:
-  using data_type  = DataType;
-  using index_type = IndexType;
-
-  class SparseRowDescriptor
-  {
-    friend class SparseMatrixDescriptor;
-    std::map<IndexType, DataType> m_id_value_map;
-
-   public:
-    IndexType
-    numberOfValues() const noexcept
-    {
-      return m_id_value_map.size();
-    }
-
-    const DataType&
-    operator[](const IndexType& j) const
-    {
-      auto i_value = m_id_value_map.find(j);
-      Assert(i_value != m_id_value_map.end());
-      return i_value->second;
-    }
-
-    DataType&
-    operator[](const IndexType& j)
-    {
-      auto i_value = m_id_value_map.find(j);
-      if (i_value != m_id_value_map.end()) {
-        return i_value->second;
-      } else {
-        auto [i_inserted, success] = m_id_value_map.insert(std::make_pair(j, 0));
-        Assert(success);
-        return i_inserted->second;
-      }
-    }
-
-    friend std::ostream&
-    operator<<(std::ostream& os, const SparseRowDescriptor& row)
-    {
-      for (auto [j, value] : row.m_id_value_map) {
-        os << ' ' << static_cast<size_t>(j) << ':' << value;
-      }
-      return os;
-    }
-
-    SparseRowDescriptor() = default;
-  };
-
- private:
-  Array<SparseRowDescriptor> m_row_array;
-  const size_t m_nb_columns;
-
- public:
-  PUGS_INLINE
-  size_t
-  numberOfRows() const noexcept
-  {
-    return m_row_array.size();
-  }
-
-  PUGS_INLINE
-  size_t
-  numberOfColumns() const noexcept
-  {
-    return m_nb_columns;
-  }
-
-  SparseRowDescriptor&
-  row(const IndexType i)
-  {
-    Assert(i < m_row_array.size());
-    return m_row_array[i];
-  }
-
-  const SparseRowDescriptor&
-  row(const IndexType i) const
-  {
-    Assert(i < m_row_array.size());
-    return m_row_array[i];
-  }
-
-  DataType&
-  operator()(const IndexType& i, const IndexType& j)
-  {
-    Assert(i < m_row_array.size());
-    Assert(j < m_nb_columns);
-    return m_row_array[i][j];
-  }
-
-  const DataType&
-  operator()(const IndexType& i, const IndexType& j) const
-  {
-    Assert(i < m_row_array.size());
-    Assert(j < m_nb_columns);
-    const auto& r = m_row_array[i];   // split to ensure const-ness of call
-    return r[j];
-  }
-
-  friend std::ostream&
-  operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
-  {
-    for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
-      os << static_cast<size_t>(i) << " |" << M.m_row_array[i] << '\n';
-    }
-    return os;
-  }
-
-  auto
-  graphVector() const
-  {
-    std::vector<std::vector<IndexType>> graph_vector(m_row_array.size());
-    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
-      const SparseRowDescriptor& row = m_row_array[i_row];
-      for (auto [id, value] : row.m_id_value_map) {
-        graph_vector[i_row].push_back(id);
-      }
-    }
-    return graph_vector;
-  }
-
-  Array<DataType>
-  valueArray() const
-  {
-    IndexType size = 0;
-    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
-      size += m_row_array[i_row].m_id_value_map.size();
-    }
-
-    Array<DataType> values(size);
-
-    IndexType cpt = 0;
-    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
-      const SparseRowDescriptor& row = m_row_array[i_row];
-      for (auto [id, value] : row.m_id_value_map) {
-        values[cpt++] = value;
-      }
-    }
-    return values;
-  }
-
-  SparseMatrixDescriptor(size_t nb_rows, size_t nb_columns) : m_row_array{nb_rows}, m_nb_columns{nb_columns}
-  {
-    if (nb_rows == nb_columns) {
-      for (size_t i = 0; i < nb_rows; ++i) {
-        m_row_array[i][i] = 0;
-      }
-    }
-  }
-
-  SparseMatrixDescriptor(size_t nb_rows) : SparseMatrixDescriptor{nb_rows, nb_rows} {}
-
-  ~SparseMatrixDescriptor() = default;
-};
-
-#endif   // SPARSE_MATRIX_DESCRIPTOR_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index a56fb1f9d0b49108f1770f186cec453758cfd793..7904eedc03c7a203667e5c604096a5cded6e136b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -56,6 +56,7 @@ add_executable (unit_tests
   test_ConcatExpressionProcessor.cpp
   test_CRSGraph.cpp
   test_CRSMatrix.cpp
+  test_CRSMatrixDescriptor.cpp
   test_DataVariant.cpp
   test_DenseMatrix.cpp
   test_Demangle.cpp
@@ -89,7 +90,6 @@ add_executable (unit_tests
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp
   test_RevisionInfo.cpp
-  test_SparseMatrixDescriptor.cpp
   test_SymbolTable.cpp
   test_Table.cpp
   test_Timer.cpp
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index e1be8069885b80a3e26ede233dac3941b7525521..231783090756e52187b46da02a97d44cf2279ea9 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -3,6 +3,7 @@
 
 #include <algebra/BiCGStab.hpp>
 #include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -10,7 +11,13 @@ TEST_CASE("BiCGStab", "[algebra]")
 {
   SECTION("no preconditionner")
   {
-    SparseMatrixDescriptor S{5};
+    Array<int> non_zeros(5);
+    non_zeros[0] = 2;
+    non_zeros[1] = 3;
+    non_zeros[2] = 3;
+    non_zeros[3] = 3;
+    non_zeros[4] = 2;
+    CRSMatrixDescriptor<double> S{5, 5, non_zeros};
     S(0, 0) = 2;
     S(0, 1) = -1;
 
@@ -29,7 +36,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     S(4, 3) = 1;
     S(4, 4) = 3;
 
-    CRSMatrix A{S};
+    CRSMatrix A{S.getCRSMatrix()};
 
     Vector<const double> x_exact = [] {
       Vector<double> y{5};
@@ -53,7 +60,9 @@ TEST_CASE("BiCGStab", "[algebra]")
 
   SECTION("no preconditionner non-converged")
   {
-    SparseMatrixDescriptor S{5};
+    Array<int> non_zeros(5);
+    non_zeros.fill(2);
+    CRSMatrixDescriptor<double> S{5, 5, non_zeros};
     S(0, 0) = 2;
     S(0, 1) = -1;
 
@@ -72,7 +81,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     S(4, 3) = 1;
     S(4, 4) = 3;
 
-    CRSMatrix A{S};
+    CRSMatrix A{S.getCRSMatrix()};
 
     Vector<const double> x_exact = [] {
       Vector<double> y{5};
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index 05919858f3ada681eb371883c6fcf6fd7d15a724..b76d90a047f2aab41e0576ab6f69ed5b133b474e 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -3,6 +3,7 @@
 
 #include <algebra/CG.hpp>
 #include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -10,7 +11,11 @@ TEST_CASE("CG", "[algebra]")
 {
   SECTION("no preconditionner")
   {
-    SparseMatrixDescriptor S{5};
+    Array<int> non_zeros{5};
+    non_zeros.fill(3);
+    non_zeros[0] = 2;
+    non_zeros[4] = 2;
+    CRSMatrixDescriptor S{5, 5, non_zeros};
     S(0, 0) = 2;
     S(0, 1) = -1;
 
@@ -29,7 +34,7 @@ TEST_CASE("CG", "[algebra]")
     S(4, 3) = -1;
     S(4, 4) = 2;
 
-    CRSMatrix A{S};
+    CRSMatrix A{S.getCRSMatrix()};
 
     Vector<const double> x_exact = [] {
       Vector<double> y{5};
@@ -53,9 +58,11 @@ TEST_CASE("CG", "[algebra]")
 
   SECTION("singular matrix with RHS in its image")
   {
-    SparseMatrixDescriptor<double, uint8_t> S{5};
+    Array<int> non_zeros(5);
+    non_zeros.fill(0);
+    CRSMatrixDescriptor<double> S{5, 5, non_zeros};
 
-    CRSMatrix<double, uint8_t> A{S};
+    CRSMatrix<double> A{S.getCRSMatrix()};
 
     Vector<double> b{5};
     b = 0;
@@ -69,7 +76,11 @@ TEST_CASE("CG", "[algebra]")
 
   SECTION("no preconditionner non-converged")
   {
-    SparseMatrixDescriptor S{5};
+    Array<int> non_zeros(5);
+    non_zeros.fill(3);
+    non_zeros[0] = 2;
+    non_zeros[4] = 2;
+    CRSMatrixDescriptor<double> S{5, 5, non_zeros};
     S(0, 0) = 2;
     S(0, 1) = -1;
 
@@ -88,7 +99,7 @@ TEST_CASE("CG", "[algebra]")
     S(4, 3) = -1;
     S(4, 4) = 2;
 
-    CRSMatrix A{S};
+    CRSMatrix A{S.getCRSMatrix()};
 
     Vector<const double> x_exact = [] {
       Vector<double> y{5};
diff --git a/tests/test_CRSMatrix.cpp b/tests/test_CRSMatrix.cpp
index d6bfdaba966c0c36a5245af457954a6f53a0919d..e0aa27b57c1a9fd773c99d7d96c6d5e47d79f74e 100644
--- a/tests/test_CRSMatrix.cpp
+++ b/tests/test_CRSMatrix.cpp
@@ -2,9 +2,7 @@
 #include <catch2/matchers/catch_matchers_all.hpp>
 
 #include <algebra/CRSMatrix.hpp>
-
-// Instantiate to ensure full coverage is performed
-template class SparseMatrixDescriptor<int, uint8_t>;
+#include <algebra/CRSMatrixDescriptor.hpp>
 
 // clazy:excludeall=non-pod-global-static
 
@@ -12,7 +10,10 @@ TEST_CASE("CRSMatrix", "[algebra]")
 {
   SECTION("matrix size")
   {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
+    Array<int> non_zeros{6};
+    non_zeros.fill(3);
+
+    CRSMatrixDescriptor<int> S{6, 5, non_zeros};
     S(0, 2) = 3;
     S(1, 2) = 11;
     S(1, 1) = 1;
@@ -22,16 +23,19 @@ TEST_CASE("CRSMatrix", "[algebra]")
     S(2, 2) = 4;
     S(4, 3) = 2;
     S(4, 4) = -2;
+    S(5, 2) = 1;
 
-    CRSMatrix<int, uint8_t> A{S};
+    CRSMatrix<int> A{S.getCRSMatrix()};
 
     REQUIRE(A.numberOfRows() == S.numberOfRows());
     REQUIRE(A.numberOfColumns() == S.numberOfColumns());
   }
 
-  SECTION("matrix vector product (simple)")
+  SECTION("square matrix vector product (simple)")
   {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
+    Array<int> non_zeros{5};
+    non_zeros.fill(3);
+    CRSMatrixDescriptor<int> S{5, 5, non_zeros};
     S(0, 2) = 3;
     S(1, 2) = 11;
     S(1, 1) = 1;
@@ -42,9 +46,9 @@ TEST_CASE("CRSMatrix", "[algebra]")
     S(4, 3) = 2;
     S(4, 4) = -2;
 
-    CRSMatrix<int, uint8_t> A{S};
+    CRSMatrix<int> A{S.getCRSMatrix()};
 
-    Vector<int> x{A.numberOfColumns()};
+    Vector<int> x(A.numberOfColumns());
     x    = 0;
     x[0] = 1;
 
@@ -96,9 +100,39 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(y[4] == -2);
   }
 
+  SECTION("rectangular matrix vector product (simple)")
+  {
+    Array<int> non_zeros{3};
+    non_zeros.fill(3);
+    CRSMatrixDescriptor<int> S{3, 5, non_zeros};
+    S(0, 2) = 3;
+    S(0, 4) = 2;
+    S(1, 2) = 11;
+    S(1, 1) = 2;
+    S(1, 0) = 1;
+    S(2, 3) = 6;
+    S(2, 2) = 4;
+
+    CRSMatrix<int> A{S.getCRSMatrix()};
+
+    Vector<int> x(A.numberOfColumns());
+    x[0] = 1;
+    x[1] = 3;
+    x[2] = 2;
+    x[3] = -3;
+    x[4] = 5;
+
+    Vector<int> y = A * x;
+    REQUIRE(y[0] == 16);
+    REQUIRE(y[1] == 29);
+    REQUIRE(y[2] == -10);
+  }
+
   SECTION("matrix vector product (complete)")
   {
-    SparseMatrixDescriptor<int, uint8_t> S{4};
+    Array<int> non_zeros{4};
+    non_zeros.fill(4);
+    CRSMatrixDescriptor<int> S{4, 4, non_zeros};
     S(0, 0) = 1;
     S(0, 1) = 2;
     S(0, 2) = 3;
@@ -116,9 +150,9 @@ TEST_CASE("CRSMatrix", "[algebra]")
     S(3, 2) = 15;
     S(3, 3) = 16;
 
-    CRSMatrix<int, uint8_t> A{S};
+    CRSMatrix<int> A{S.getCRSMatrix()};
 
-    Vector<int> x{A.numberOfColumns()};
+    Vector<int> x(A.numberOfColumns());
     x[0] = 1;
     x[1] = 2;
     x[2] = 3;
@@ -133,7 +167,10 @@ TEST_CASE("CRSMatrix", "[algebra]")
 
   SECTION("check values")
   {
-    SparseMatrixDescriptor<int, uint8_t> S{4};
+    Array<int> non_zeros{4};
+    non_zeros.fill(3);
+    non_zeros[2] = 4;
+    CRSMatrixDescriptor<int> S{4, 4, non_zeros};
     S(3, 0) = 13;
     S(0, 0) = 1;
     S(0, 1) = 2;
@@ -148,7 +185,7 @@ TEST_CASE("CRSMatrix", "[algebra]")
     S(1, 0) = 5;
     S(2, 1) = 10;
 
-    CRSMatrix<int, uint8_t> A{S};
+    CRSMatrix<int> A{S.getCRSMatrix()};
 
     auto values = A.values();
     REQUIRE(values.size() == 13);
@@ -166,29 +203,46 @@ TEST_CASE("CRSMatrix", "[algebra]")
     REQUIRE(values[11] == 15);
     REQUIRE(values[12] == 16);
 
-    auto row_indices = A.rowIndices();
-
-    REQUIRE(row_indices.size() == 5);
-
-    REQUIRE(A.row(0).colidx(0) == 0);
-    REQUIRE(A.row(0).colidx(1) == 1);
-    REQUIRE(A.row(0).colidx(2) == 3);
-    REQUIRE(A.row(1).colidx(0) == 0);
-    REQUIRE(A.row(1).colidx(1) == 1);
-    REQUIRE(A.row(1).colidx(2) == 2);
-    REQUIRE(A.row(2).colidx(0) == 0);
-    REQUIRE(A.row(2).colidx(1) == 1);
-    REQUIRE(A.row(2).colidx(2) == 2);
-    REQUIRE(A.row(2).colidx(3) == 3);
-    REQUIRE(A.row(3).colidx(0) == 0);
-    REQUIRE(A.row(3).colidx(1) == 2);
-    REQUIRE(A.row(3).colidx(2) == 3);
+    auto row_map = A.rowMap();
+    REQUIRE(row_map.size() == 5);
+    REQUIRE(row_map[0] == 0);
+    REQUIRE(row_map[1] == 3);
+    REQUIRE(row_map[2] == 6);
+    REQUIRE(row_map[3] == 10);
+    REQUIRE(row_map[4] == 13);
+
+    auto column_indices = A.columnIndices();
+    REQUIRE(column_indices.size() == 13);
+    REQUIRE(column_indices[0] == 0);
+    REQUIRE(column_indices[1] == 1);
+    REQUIRE(column_indices[2] == 3);
+    REQUIRE(column_indices[3] == 0);
+    REQUIRE(column_indices[4] == 1);
+    REQUIRE(column_indices[5] == 2);
+    REQUIRE(column_indices[6] == 0);
+    REQUIRE(column_indices[7] == 1);
+    REQUIRE(column_indices[8] == 2);
+    REQUIRE(column_indices[9] == 3);
+    REQUIRE(column_indices[10] == 0);
+    REQUIRE(column_indices[11] == 2);
+    REQUIRE(column_indices[12] == 3);
+
+    std::ostringstream crs_descriptor_os;
+    std::ostringstream crs_os;
+
+    crs_descriptor_os << S;
+    crs_os << A;
+
+    REQUIRE(crs_descriptor_os.str() == crs_os.str());
   }
 
 #ifndef NDEBUG
   SECTION("incompatible runtime matrix/vector product")
   {
-    CRSMatrix<int, uint8_t> A{SparseMatrixDescriptor<int, uint8_t>{4}};
+    Array<int> non_zeros(2);
+    non_zeros.fill(0);
+    CRSMatrixDescriptor<int> S{2, 4, non_zeros};
+    CRSMatrix<int> A{S.getCRSMatrix()};
     Vector<int> x{2};
     REQUIRE_THROWS_AS(A * x, AssertError);
   }
diff --git a/tests/test_CRSMatrixDescriptor.cpp b/tests/test_CRSMatrixDescriptor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2bf59af22201ee9a8f164d2d6344b16357c56303
--- /dev/null
+++ b/tests/test_CRSMatrixDescriptor.cpp
@@ -0,0 +1,166 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <algebra/CRSMatrixDescriptor.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class CRSMatrixDescriptor<int>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("CRSMatrixDescriptor", "[algebra]")
+{
+  SECTION("has overflow / not filled")
+  {
+    const size_t nb_lines   = 2;
+    const size_t nb_columns = 5;
+
+    Array<int> non_zeros{nb_lines};
+    non_zeros.fill(2);
+    CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+    S(1, 3) = 5;
+    S(1, 1) += 3;
+    S(0, 1) = 4;
+    S(1, 0) = 2;
+    S(1, 2) = 7;
+
+    REQUIRE(S.numberOfRows() == 2);
+    REQUIRE(S.numberOfColumns() == 5);
+
+    REQUIRE(S.isFilled() == false);
+    REQUIRE(S.hasOverflow() == true);
+
+    std::ostringstream os;
+    os << S;
+
+    std::string_view output =
+      R"(0| 1:4
+1| 0:2 1:3 2:7 3:5
+)";
+    REQUIRE(os.str() == output);
+  }
+
+  SECTION("no overflow / not filled")
+  {
+    const size_t nb_lines   = 2;
+    const size_t nb_columns = 5;
+
+    Array<int> non_zeros{nb_lines};
+    non_zeros.fill(2);
+    CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+    S(1, 3) = 5;
+    S(1, 1) += 3;
+    S(0, 1) = 4;
+
+    REQUIRE(S.numberOfRows() == 2);
+    REQUIRE(S.numberOfColumns() == 5);
+
+    REQUIRE(S.isFilled() == false);
+    REQUIRE(S.hasOverflow() == false);
+
+    std::ostringstream os;
+    os << S;
+
+    std::string_view output =
+      R"(0| 1:4
+1| 1:3 3:5
+)";
+    REQUIRE(os.str() == output);
+  }
+
+  SECTION("no overflow / filled")
+  {
+    const size_t nb_lines   = 2;
+    const size_t nb_columns = 5;
+
+    Array<int> non_zeros{nb_lines};
+    non_zeros[0] = 1;
+    non_zeros[1] = 5;
+    CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+    S(1, 3) = 2;
+    S(1, 3) += 3;
+    S(1, 1) += 3;
+    S(0, 1) = 4;
+    S(1, 2) = 2;
+    S(1, 0) = 1;
+    S(1, 4) = 13;
+
+    REQUIRE(S.numberOfRows() == 2);
+    REQUIRE(S.numberOfColumns() == 5);
+
+    REQUIRE(S.isFilled() == true);
+    REQUIRE(S.hasOverflow() == false);
+
+    std::ostringstream os;
+    os << S;
+
+    std::string_view output =
+      R"(0| 1:4
+1| 0:1 1:3 2:2 3:5 4:13
+)";
+    REQUIRE(os.str() == output);
+  }
+
+  SECTION("has overflow / filled")
+  {
+    const size_t nb_lines   = 2;
+    const size_t nb_columns = 5;
+
+    Array<int> non_zeros{nb_lines};
+    non_zeros.fill(2);
+    CRSMatrixDescriptor<int> S(nb_lines, nb_columns, non_zeros);
+
+    S(1, 3) = 3;
+    S(1, 0) = -2;
+    S(1, 1) += 1;
+    S(1, 1) += 1;
+    S(1, 1) = 3;
+    S(0, 1) = 4;
+    S(0, 0) = 1;
+    S(0, 3) = 5;
+    S(0, 2) = 1;
+
+    REQUIRE(S.numberOfRows() == 2);
+    REQUIRE(S.numberOfColumns() == 5);
+
+    REQUIRE(S.isFilled() == true);
+    REQUIRE(S.hasOverflow() == true);
+
+    std::ostringstream os;
+    os << S;
+
+    std::string_view output =
+      R"(0| 0:1 1:4 2:1 3:5
+1| 0:-2 1:3 3:3
+)";
+    REQUIRE(os.str() == output);
+  }
+
+#ifndef NDEBUG
+  SECTION("incompatible non zero vector size and number of columns")
+  {
+    Array<int> non_zeros(3);
+    REQUIRE_THROWS_AS((CRSMatrixDescriptor<int>{2, 4, non_zeros}), AssertError);
+  }
+
+  SECTION("bad row number")
+  {
+    Array<int> non_zeros(2);
+    CRSMatrixDescriptor<int> S{2, 4, non_zeros};
+    REQUIRE_THROWS_AS(S(2, 3) = 1, AssertError);
+  }
+
+  SECTION("bad column number")
+  {
+    Array<int> non_zeros(2);
+    CRSMatrixDescriptor<int> S{2, 4, non_zeros};
+    REQUIRE_THROWS_AS(S(0, 5) = 1, AssertError);
+  }
+
+#endif   // NDEBUG
+}
diff --git a/tests/test_EigenvalueSolver.cpp b/tests/test_EigenvalueSolver.cpp
index e779bce5600e32cf26b0fde7339693f8da69bcb1..3ca8f7f5136d72f7003357113b7b597d1b84ef54 100644
--- a/tests/test_EigenvalueSolver.cpp
+++ b/tests/test_EigenvalueSolver.cpp
@@ -2,9 +2,9 @@
 #include <catch2/matchers/catch_matchers_all.hpp>
 
 #include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
 #include <algebra/DenseMatrix.hpp>
 #include <algebra/EigenvalueSolver.hpp>
-#include <algebra/SparseMatrixDescriptor.hpp>
 
 #include <utils/pugs_config.hpp>
 
@@ -16,7 +16,10 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
   {
     SECTION("symmetric system")
     {
-      SparseMatrixDescriptor S{3};
+      Array<int> non_zeros(3);
+      non_zeros.fill(3);
+      CRSMatrixDescriptor<double> S{3, 3, non_zeros};
+
       S(0, 0) = 3;
       S(0, 1) = 2;
       S(0, 2) = 4;
@@ -29,7 +32,7 @@ TEST_CASE("EigenvalueSolver", "[algebra]")
       S(2, 1) = 2;
       S(2, 2) = 3;
 
-      CRSMatrix A{S};
+      CRSMatrix A{S.getCRSMatrix()};
 
       SECTION("eigenvalues")
       {
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 1219c48e62d28a749383e115aef8b6e51a17e5ec..a9b4cab1fe49efafa6a8ec4d38e34b7b947c6eff 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -3,6 +3,7 @@
 
 #include <utils/pugs_config.hpp>
 
+#include <algebra/CRSMatrixDescriptor.hpp>
 #include <algebra/LinearSolver.hpp>
 
 // clazy:excludeall=non-pod-global-static
@@ -147,7 +148,11 @@ TEST_CASE("LinearSolver", "[algebra]")
     {
       SECTION("symmetric system")
       {
-        SparseMatrixDescriptor S{5};
+        Array<int> non_zeros{5};
+        non_zeros.fill(3);
+        non_zeros[0] = 2;
+        non_zeros[4] = 2;
+        CRSMatrixDescriptor S{5, 5, non_zeros};
         S(0, 0) = 2;
         S(0, 1) = -1;
 
@@ -166,7 +171,7 @@ TEST_CASE("LinearSolver", "[algebra]")
         S(4, 3) = -1;
         S(4, 4) = 2;
 
-        CRSMatrix A{S};
+        CRSMatrix A{S.getCRSMatrix()};
 
         Vector<const double> x_exact = [] {
           Vector<double> y{5};
@@ -302,7 +307,12 @@ TEST_CASE("LinearSolver", "[algebra]")
 
       SECTION("none symmetric system")
       {
-        SparseMatrixDescriptor S{5};
+        Array<int> non_zeros{5};
+        non_zeros.fill(3);
+        non_zeros[0] = 2;
+        non_zeros[4] = 2;
+        CRSMatrixDescriptor S{5, 5, non_zeros};
+
         S(0, 0) = 2;
         S(0, 1) = -1;
 
@@ -321,7 +331,7 @@ TEST_CASE("LinearSolver", "[algebra]")
         S(4, 3) = 1;
         S(4, 4) = 3;
 
-        CRSMatrix A{S};
+        CRSMatrix A{S.getCRSMatrix()};
 
         Vector<const double> x_exact = [] {
           Vector<double> y{5};
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
deleted file mode 100644
index 77bc890f7c4a53ea170de08bebb83de9b7e0fbed..0000000000000000000000000000000000000000
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ /dev/null
@@ -1,201 +0,0 @@
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_all.hpp>
-
-#include <utils/PugsAssert.hpp>
-
-#include <algebra/SparseMatrixDescriptor.hpp>
-
-// Instantiate to ensure full coverage is performed
-template class SparseMatrixDescriptor<int, uint8_t>;
-
-// clazy:excludeall=non-pod-global-static
-
-TEST_CASE("SparseMatrixDescriptor", "[algebra]")
-{
-  SECTION("SparseRowDescriptor subclass")
-  {
-    SECTION("number of values")
-    {
-      SparseMatrixDescriptor<int, uint8_t>::SparseRowDescriptor r;
-      REQUIRE(r.numberOfValues() == 0);
-
-      r[0] = 3;
-      r[0] += 2;
-
-      REQUIRE(r.numberOfValues() == 1);
-
-      r[1] = -2;
-
-      REQUIRE(r.numberOfValues() == 2);
-    }
-
-    SECTION("values access")
-    {
-      SparseMatrixDescriptor<int, uint8_t>::SparseRowDescriptor r;
-
-      r[0] = 3;
-      r[0] += 2;
-
-      r[3] = 8;
-
-      REQUIRE(r[0] == 5);
-      REQUIRE(r[3] == 8);
-
-      const auto& const_r = r;
-      REQUIRE(const_r[0] == 5);
-      REQUIRE(const_r[3] == 8);
-
-#ifndef NDEBUG
-      REQUIRE_THROWS_AS(const_r[2], AssertError);
-#endif   // NDEBUG
-    }
-  }
-
-  SECTION("matrix size")
-  {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
-    REQUIRE(S.numberOfRows() == 5);
-    REQUIRE(S.numberOfColumns() == 5);
-  }
-
-  SECTION("location operators")
-  {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
-    S(0, 2) = 3;
-    S(1, 1) = 1;
-    S(1, 2) = 11;
-    S(0, 2) += 2;
-    S(2, 2) = 4;
-    S(3, 1) = 5;
-    S(4, 1) = 1;
-    S(4, 4) = -2;
-
-    REQUIRE(S.row(0).numberOfValues() == 2);
-    REQUIRE(S.row(1).numberOfValues() == 2);
-    REQUIRE(S.row(2).numberOfValues() == 1);
-    REQUIRE(S.row(3).numberOfValues() == 2);
-    REQUIRE(S.row(4).numberOfValues() == 2);
-
-    const auto& const_S = S;
-
-    REQUIRE(const_S.row(0).numberOfValues() == 2);
-    REQUIRE(const_S.row(1).numberOfValues() == 2);
-    REQUIRE(const_S.row(2).numberOfValues() == 1);
-    REQUIRE(const_S.row(3).numberOfValues() == 2);
-    REQUIRE(const_S.row(4).numberOfValues() == 2);
-
-#ifndef NDEBUG
-    REQUIRE_THROWS_AS(S.row(5), AssertError);
-    REQUIRE_THROWS_AS(const_S.row(5), AssertError);
-#endif   // NDEBUG
-
-    REQUIRE(S(0, 2) == 5);
-    REQUIRE(S(1, 1) == 1);
-    REQUIRE(S(1, 2) == 11);
-    REQUIRE(S(2, 2) == 4);
-    REQUIRE(S(3, 1) == 5);
-    REQUIRE(S(4, 1) == 1);
-    REQUIRE(S(4, 4) == -2);
-
-    REQUIRE(const_S(0, 2) == 5);
-    REQUIRE(const_S(1, 1) == 1);
-    REQUIRE(const_S(1, 2) == 11);
-    REQUIRE(const_S(2, 2) == 4);
-    REQUIRE(const_S(3, 1) == 5);
-    REQUIRE(const_S(4, 1) == 1);
-    REQUIRE(const_S(4, 4) == -2);
-
-#ifndef NDEBUG
-    REQUIRE_THROWS_AS(S(5, 0), AssertError);
-    REQUIRE_THROWS_AS(const_S(6, 1), AssertError);
-    REQUIRE_THROWS_AS(const_S(0, 1), AssertError);
-#endif   // NDEBUG
-  }
-
-  SECTION("vector-graph")
-  {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
-    S(0, 2) = 3;
-    S(1, 2) = 11;
-    S(1, 1) = 1;
-    S(0, 2) += 2;
-    S(3, 3) = 5;
-    S(3, 1) = 5;
-    S(4, 1) = 1;
-    S(2, 2) = 4;
-    S(4, 4) = -2;
-
-    const auto graph = S.graphVector();
-
-    REQUIRE(graph.size() == S.numberOfRows());
-    REQUIRE(graph[0].size() == 2);
-    REQUIRE(graph[1].size() == 2);
-    REQUIRE(graph[2].size() == 1);
-    REQUIRE(graph[3].size() == 2);
-    REQUIRE(graph[4].size() == 2);
-
-    REQUIRE(graph[0][0] == 0);
-    REQUIRE(graph[0][1] == 2);
-    REQUIRE(graph[1][0] == 1);
-    REQUIRE(graph[1][1] == 2);
-    REQUIRE(graph[2][0] == 2);
-    REQUIRE(graph[3][0] == 1);
-    REQUIRE(graph[3][1] == 3);
-    REQUIRE(graph[4][0] == 1);
-    REQUIRE(graph[4][1] == 4);
-  }
-
-  SECTION("value array")
-  {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
-    S(0, 2) = 3;
-    S(1, 2) = 11;
-    S(1, 1) = 1;
-    S(0, 2) += 2;
-    S(3, 3) = 5;
-    S(3, 1) = -3;
-    S(4, 1) = 1;
-    S(2, 2) = 4;
-    S(4, 4) = -2;
-
-    const auto value_array = S.valueArray();
-
-    REQUIRE(value_array.size() == 9);
-
-    REQUIRE(value_array[0] == 0);
-    REQUIRE(value_array[1] == 5);
-    REQUIRE(value_array[2] == 1);
-    REQUIRE(value_array[3] == 11);
-    REQUIRE(value_array[4] == 4);
-    REQUIRE(value_array[5] == -3);
-    REQUIRE(value_array[6] == 5);
-    REQUIRE(value_array[7] == 1);
-    REQUIRE(value_array[8] == -2);
-  }
-
-  SECTION("output")
-  {
-    SparseMatrixDescriptor<int, uint8_t> S{5};
-    S(0, 2) = 3;
-    S(1, 2) = 11;
-    S(1, 1) = 1;
-    S(0, 2) += 2;
-    S(3, 3) = 5;
-    S(3, 1) = -3;
-    S(4, 1) = 1;
-    S(2, 2) = 4;
-    S(4, 4) = -2;
-
-    std::ostringstream output;
-    output << '\n' << S;
-
-    std::string expected_output = R"(
-0 | 0:0 2:5
-1 | 1:1 2:11
-2 | 2:4
-3 | 1:-3 3:5
-4 | 1:1 4:-2
-)";
-    REQUIRE(output.str() == expected_output);
-  }
-}