diff --git a/MathLib/LinAlg/Eigen/EigenMatrix.h b/MathLib/LinAlg/Eigen/EigenMatrix.h
index 4cbd970c03d1b4800ffd2ed887dd4e1fddc0d20c..44e15aff2d1604572bff397e5ee5845fbd5656f9 100644
--- a/MathLib/LinAlg/Eigen/EigenMatrix.h
+++ b/MathLib/LinAlg/Eigen/EigenMatrix.h
@@ -18,6 +18,7 @@
 #include <Eigen/Sparse>
 
 #include "MathLib/LinAlg/RowColumnIndices.h"
+#include "MathLib/LinAlg/SetMatrixSparsity.h"
 #include "EigenVector.h"
 
 namespace MathLib
@@ -184,6 +185,27 @@ void EigenMatrix::add(std::vector<IndexType> const& row_pos,
     }
 };
 
+/// Sets the sparsity pattern of the underlying EigenMatrix.
+template <typename SPARSITY_PATTERN>
+struct SetMatrixSparsity<EigenMatrix, SPARSITY_PATTERN>
+{
+
+/// \note This operator relies on row-major storage order of the underlying
+/// eigen matrix i.e. of the RawMatrixType.
+void operator()(EigenMatrix &matrix, SPARSITY_PATTERN const& sparsity_pattern)
+{
+    static_assert(EigenMatrix::RawMatrixType::IsRowMajor,
+                  "Set matrix sparsity relies on the EigenMatrix to be in "
+                  "row-major storage order.");
+    auto const n_rows = matrix.getNRows();
+    Eigen::VectorXi row_sizes(n_rows);
+
+    for (auto i = decltype(n_rows){0}; i < n_rows; i++)
+        row_sizes[i] = sparsity_pattern.getNodeDegree(i);
+
+    matrix.getRawMatrix().reserve(row_sizes);
+}
+};
 
 } // end namespace MathLib
 
diff --git a/MathLib/LinAlg/SetMatrixSparsity.h b/MathLib/LinAlg/SetMatrixSparsity.h
index f1ba8bb5eae15d95db878e1d8562375b5dd9a454..c13f171d863f3b6f11c516a63bb3e5d4ddc64cef 100644
--- a/MathLib/LinAlg/SetMatrixSparsity.h
+++ b/MathLib/LinAlg/SetMatrixSparsity.h
@@ -38,4 +38,8 @@ void setMatrixSparsity(MATRIX& matrix, SPARSITY_PATTERN const& sparsity_pattern)
 #include "Lis/LisMatrix.h"
 #endif  // USE_LIS
 
+#ifdef OGS_USE_EIGEN
+#include "Eigen/EigenMatrix.h"
+#endif  // OGS_USE_EIGEN
+
 #endif  // MATHLIB_SETMATRIXSPARSITY_H_