diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index 010ebd95d037212db2d0fb55adcf85e9bc76b04c..491a35242b97d5dd32e76b9790c97fb0eb26b977 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -5,129 +5,12 @@
 #include <Kokkos_StaticCrsGraph.hpp>
 #include <PugsAssert.hpp>
 
+#include <SparseMatrixDescriptor.hpp>
+
 #include <Vector.hpp>
 
 #include <iostream>
 
-template <typename DataType = double, typename IndexType = size_t>
-class FlexibleMatrix
-{
- public:
-  using data_type  = DataType;
-  using index_type = IndexType;
-
-  class FlexibleRow
-  {
-    friend class FlexibleMatrix;
-    std::map<IndexType, DataType> m_id_value_map;
-
-   public:
-    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 FlexibleRow& row)
-    {
-      for (auto [j, value] : row.m_id_value_map) {
-        os << ' ' << j << ':' << value;
-      }
-      return os;
-    }
-
-    FlexibleRow() = default;
-  };
-
- private:
-  Array<FlexibleRow> m_row_array;
-
- public:
-  FlexibleRow&
-  row(const IndexType i)
-  {
-    return m_row_array[i];
-  }
-
-  const FlexibleRow&
-  row(const IndexType i) const
-  {
-    return m_row_array[i];
-  }
-
-  DataType&
-  operator()(const IndexType& i, const IndexType& j)
-  {
-    return m_row_array[i][j];
-  }
-
-  const DataType&
-  operator()(const IndexType& i, const IndexType& j) const
-  {
-    return m_row_array[i][j];
-  }
-
-  friend std::ostream&
-  operator<<(std::ostream& os, const FlexibleMatrix& M)
-  {
-    for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
-      os << 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 FlexibleRow& 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 FlexibleRow& row = m_row_array[i_row];
-      for (auto [id, value] : row.m_id_value_map) {
-        values[cpt++] = value;
-      }
-    }
-    return values;
-  }
-
-  FlexibleMatrix(IndexType nb_row) : m_row_array{nb_row} {}
-
-  ~FlexibleMatrix() = default;
-};
-
 template <typename DataType = double, typename IndexType = size_t>
 class CRSMatrix
 {
@@ -173,7 +56,7 @@ class CRSMatrix
     return os;
   }
 
-  CRSMatrix(const FlexibleMatrix<DataType, IndexType>& M)
+  CRSMatrix(const SparseMatrixDescriptor<DataType, IndexType>& M)
   {
     m_host_matrix = Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", M.graphVector());
     m_values      = M.valueArray();
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a36121aff6e099fe1298f7c1ef9c4bda04784dda
--- /dev/null
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -0,0 +1,126 @@
+#ifndef SPARSE_MATRIX_DESCRIPTOR_HPP
+#define SPARSE_MATRIX_DESCRIPTOR_HPP
+
+#include <Array.hpp>
+#include <map>
+
+template <typename DataType = double, typename IndexType = size_t>
+class SparseMatrixDescriptor
+{
+ public:
+  using data_type  = DataType;
+  using index_type = IndexType;
+
+  class SparseRowDescriptor
+  {
+    friend class SparseMatrixDescriptor;
+    std::map<IndexType, DataType> m_id_value_map;
+
+   public:
+    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 << ' ' << j << ':' << value;
+      }
+      return os;
+    }
+
+    SparseRowDescriptor() = default;
+  };
+
+ private:
+  Array<SparseRowDescriptor> m_row_array;
+
+ public:
+  SparseRowDescriptor&
+  row(const IndexType i)
+  {
+    return m_row_array[i];
+  }
+
+  const SparseRowDescriptor&
+  row(const IndexType i) const
+  {
+    return m_row_array[i];
+  }
+
+  DataType&
+  operator()(const IndexType& i, const IndexType& j)
+  {
+    return m_row_array[i][j];
+  }
+
+  const DataType&
+  operator()(const IndexType& i, const IndexType& j) const
+  {
+    return m_row_array[i][j];
+  }
+
+  friend std::ostream&
+  operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
+  {
+    for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
+      os << 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(IndexType nb_row) : m_row_array{nb_row} {}
+
+  ~SparseMatrixDescriptor() = default;
+};
+
+#endif   // SPARSE_MATRIX_DESCRIPTOR_HPP
diff --git a/src/main.cpp b/src/main.cpp
index c2f6b3a32d2d869611c753e69d501c2c29dbb225..e2802e4654f6024ce88dd082c81113f1441982f2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -44,17 +44,17 @@ main(int argc, char* argv[])
     {
       {   // PCG simple test
         size_t matrix_size = 10;
-        FlexibleMatrix F(matrix_size);
+        SparseMatrixDescriptor D(matrix_size);
 
         for (size_t i = 0; i < matrix_size; ++i) {
-          F(i, i) = 4;
+          D(i, i) = 4;
           if (i + 1 < matrix_size) {
-            F(i, i + 1) = -1;
-            F(i + 1, i) = -1;
+            D(i, i + 1) = -1;
+            D(i + 1, i) = -1;
           }
         }
 
-        CRSMatrix A{F};
+        CRSMatrix A{D};
 
         Vector<double> x{10};
         for (size_t i = 0; i < x.size(); ++i) {
@@ -77,17 +77,17 @@ main(int argc, char* argv[])
 
       {   // BiCGStab simple test
         size_t matrix_size = 10;
-        FlexibleMatrix F(matrix_size);
+        SparseMatrixDescriptor D(matrix_size);
 
         for (size_t i = 0; i < matrix_size; ++i) {
-          F(i, i) = 4;
+          D(i, i) = 4;
           if (i + 1 < matrix_size) {
-            F(i, i + 1) = -1;
-            F(i + 1, i) = -0.3;
+            D(i, i + 1) = -1;
+            D(i + 1, i) = -0.3;
           }
         }
 
-        CRSMatrix A{F};
+        CRSMatrix A{D};
 
         Vector<double> x{10};
         for (size_t i = 0; i < x.size(); ++i) {