diff --git a/src/algebra/DenseMatrix.hpp b/src/algebra/DenseMatrix.hpp
index 25f607e7bf0bba35d8b919655e63f6e9431724fd..952342fedbd04246879356d67425275d471af2c9 100644
--- a/src/algebra/DenseMatrix.hpp
+++ b/src/algebra/DenseMatrix.hpp
@@ -9,6 +9,8 @@
 #include <utils/PugsUtils.hpp>
 #include <utils/Types.hpp>
 
+#include <iostream>
+
 template <typename DataType>
 class DenseMatrix   // LCOV_EXCL_LINE
 {
@@ -246,6 +248,19 @@ class DenseMatrix   // LCOV_EXCL_LINE
   PUGS_INLINE
   DenseMatrix& operator=(DenseMatrix&&) = default;
 
+  friend std::ostream&
+  operator<<(std::ostream& os, const DenseMatrix& A)
+  {
+    for (size_t i = 0; i < A.numberOfRows(); ++i) {
+      os << i << '|';
+      for (size_t j = 0; j < A.numberOfColumns(); ++j) {
+        os << ' ' << j << ':' << NaNHelper(A(i, j));
+      }
+      os << '\n';
+    }
+    return os;
+  }
+
   template <typename DataType2>
   DenseMatrix(const DenseMatrix<DataType2>& A)
   {
diff --git a/tests/test_DenseMatrix.cpp b/tests/test_DenseMatrix.cpp
index 938ebefbda40aae6ff6d7a3c33f11355c51b23c1..b794689a09793f0bd844144c9b7ccdb39338f648 100644
--- a/tests/test_DenseMatrix.cpp
+++ b/tests/test_DenseMatrix.cpp
@@ -373,7 +373,44 @@ TEST_CASE("DenseMatrix", "[algebra]")
     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);