diff --git a/src/algebra/CRSMatrix.hpp b/src/algebra/CRSMatrix.hpp
index bf650afe2784035023fbb5311c1c3a12c68b2a3b..1f7fe37e571447b104b81ff61cf38dcf2144d897 100644
--- a/src/algebra/CRSMatrix.hpp
+++ b/src/algebra/CRSMatrix.hpp
@@ -1,11 +1,11 @@
 #ifndef CRS_MATRIX_HPP
 #define CRS_MATRIX_HPP
 
+#include <algebra/Vector.hpp>
 #include <utils/Array.hpp>
+#include <utils/NaNHelper.hpp>
 #include <utils/PugsAssert.hpp>
 
-#include <algebra/Vector.hpp>
-
 #include <iostream>
 
 template <typename DataType, typename IndexType>
@@ -100,7 +100,7 @@ class CRSMatrix
       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 << ' ' << A.m_column_indices[j] << ':' << NaNHelper(A.m_values[j]);
       }
       os << '\n';
     }
diff --git a/src/algebra/CRSMatrixDescriptor.hpp b/src/algebra/CRSMatrixDescriptor.hpp
index e367db354e356035d18a8742b2410a8781dde09a..4be6a323630a9845ecb142b32f73ceca5d6dc7f6 100644
--- a/src/algebra/CRSMatrixDescriptor.hpp
+++ b/src/algebra/CRSMatrixDescriptor.hpp
@@ -4,6 +4,7 @@
 #include <algebra/CRSMatrix.hpp>
 #include <utils/Array.hpp>
 #include <utils/Exceptions.hpp>
+#include <utils/NaNHelper.hpp>
 
 #include <map>
 
@@ -105,10 +106,10 @@ class CRSMatrixDescriptor
           }();
 
           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];
+            os << ' ' << A.m_column_indices[A.m_row_map[i] + j] << ':' << NaNHelper(A.m_values[A.m_row_map[i] + j]);
             ++j;
           } else {
-            os << ' ' << i_overflow->first << ':' << i_overflow->second;
+            os << ' ' << i_overflow->first << ':' << NaNHelper(i_overflow->second);
             ++i_overflow;
           }
         }
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/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 6bf810cde020681da88afedb6a815e54a563504d..98030dd447ee6cac14b6de0f3c45d55b68cdf7ba 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -3,6 +3,7 @@
 
 #include <algebra/TinyVector.hpp>
 #include <utils/InvalidData.hpp>
+#include <utils/NaNHelper.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/Types.hpp>
@@ -134,9 +135,9 @@ class [[nodiscard]] TinyMatrix
   {
     os << '[';
     for (size_t i = 0; i < N; ++i) {
-      os << '(' << A(i, 0);
+      os << '(' << NaNHelper(A(i, 0));
       for (size_t j = 1; j < N; ++j) {
-        os << ',' << A(i, j);
+        os << ',' << NaNHelper(A(i, j));
       }
       os << ')';
     }
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index 69ecb14e39d626b1ed0e2709b6fd018af92714da..382741093f351cea80ff1e261399df9583dc4d0c 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -2,6 +2,7 @@
 #define TINY_VECTOR_HPP
 
 #include <utils/InvalidData.hpp>
+#include <utils/NaNHelper.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/Types.hpp>
@@ -98,9 +99,9 @@ class [[nodiscard]] TinyVector
   PUGS_INLINE
   constexpr friend std::ostream& operator<<(std::ostream& os, const TinyVector& v)
   {
-    os << '(' << v.m_values[0];
+    os << '(' << NaNHelper(v.m_values[0]);
     for (size_t i = 1; i < N; ++i) {
-      os << ',' << v.m_values[i];
+      os << ',' << NaNHelper(v.m_values[i]);
     }
     os << ')';
     return os;
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index 2a3e1bd296a4b07e0dd32e9f3017892a2cea290f..c918ed409746a3c212720ad8ae6769cfb0882ba9 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -28,10 +28,7 @@ class Vector   // LCOV_EXCL_LINE
   friend std::ostream&
   operator<<(std::ostream& os, const Vector& x)
   {
-    for (size_t i = 0; i < x.size(); ++i) {
-      os << ' ' << x[i];
-    }
-    return os;
+    return os << x.m_values;
   }
 
   friend Vector<std::remove_const_t<DataType>>
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index e2a4e5c71c6938e96a2d2cf63318b1f7ed9c8f4c..54786db88592d7371e7b7fa5e9acd0aaec4fffe5 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -2,6 +2,7 @@
 #define ARRAY_HPP
 
 #include <utils/InvalidData.hpp>
+#include <utils/NaNHelper.hpp>
 #include <utils/PugsAssert.hpp>
 #include <utils/PugsMacros.hpp>
 #include <utils/PugsUtils.hpp>
@@ -109,6 +110,17 @@ class [[nodiscard]] Array
 #endif   // NDEBUG
   }
 
+  friend std::ostream& operator<<(std::ostream& os, const Array& x)
+  {
+    if (x.size() > 0) {
+      os << 0 << ':' << NaNHelper(x[0]);
+    }
+    for (size_t i = 1; i < x.size(); ++i) {
+      os << ' ' << i << ':' << NaNHelper(x[i]);
+    }
+    return os;
+  }
+
   PUGS_INLINE
   Array() = default;
 
diff --git a/src/utils/NaNHelper.hpp b/src/utils/NaNHelper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c93ed333dc5faf8940ab2b5ccfcac88e6720d6b8
--- /dev/null
+++ b/src/utils/NaNHelper.hpp
@@ -0,0 +1,43 @@
+#ifndef NAN_HELPER_HPP
+#define NAN_HELPER_HPP
+
+#include <utils/PugsMacros.hpp>
+
+#include <cmath>
+#include <iostream>
+#include <type_traits>
+
+template <typename T>
+class NaNHelper
+{
+ private:
+  const T& m_t;
+
+ public:
+  PUGS_INLINE
+  friend std::ostream&
+  operator<<(std::ostream& os, const NaNHelper<T>& helper)
+  {
+    if constexpr (std::is_arithmetic_v<T>) {
+      if (std::isnan(helper.m_t)) {
+        os << "nan";
+      } else {
+        os << helper.m_t;
+      }
+    } else {
+      os << helper.m_t;
+    }
+    return os;
+  }
+
+  PUGS_INLINE
+  explicit NaNHelper(const T& t) : m_t{t} {}
+
+  NaNHelper(NaNHelper&&)      = delete;
+  NaNHelper(const NaNHelper&) = delete;
+
+  NaNHelper()  = delete;
+  ~NaNHelper() = default;
+};
+
+#endif   // NAN_HELPER_HPP
diff --git a/src/utils/Table.hpp b/src/utils/Table.hpp
index 5f621ab2719e93ea63d3be610818aab1cb164665..cd380661b5b59d33f20812f969f9df7b93cd5a73 100644
--- a/src/utils/Table.hpp
+++ b/src/utils/Table.hpp
@@ -9,6 +9,8 @@
 
 #include <Kokkos_CopyViews.hpp>
 
+#include <iostream>
+
 template <typename DataType>
 class [[nodiscard]] Table
 {
@@ -126,6 +128,18 @@ class [[nodiscard]] Table
 #endif   // NDEBUG
   }
 
+  friend std::ostream& operator<<(std::ostream& os, const Table& t)
+  {
+    for (size_t i = 0; i < t.numberOfRows(); ++i) {
+      os << i << '|';
+      for (size_t j = 0; j < t.numberOfColumns(); ++j) {
+        os << ' ' << j << ':' << NaNHelper(t(i, j));
+      }
+      os << '\n';
+    }
+    return os;
+  }
+
   PUGS_INLINE
   Table() = default;
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8cea83f7294e107b28db9f223cd989c33a5a6c2d..62a70aad353dbb9d9deb7fcc7688865ab1987024 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -85,6 +85,7 @@ add_executable (unit_tests
   test_ListAffectationProcessor.cpp
   test_MathModule.cpp
   test_NameProcessor.cpp
+  test_NaNHelper.cpp
   test_OStreamProcessor.cpp
   test_ParseError.cpp
   test_PETScUtils.cpp
diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
index dd6f57ea61bd52632477cb16fc7f6556757036f5..5b9889f396f6413b4c704e48ffa1497bf6edb8fe 100644
--- a/tests/test_Array.cpp
+++ b/tests/test_Array.cpp
@@ -232,7 +232,44 @@ TEST_CASE("Array", "[utils]")
     }
   }
 
+  SECTION("output")
+  {
+    Array<int> x{5};
+    x[0] = 2;
+    x[1] = 6;
+    x[2] = 2;
+    x[3] = 3;
+    x[4] = 7;
+
+    std::ostringstream array_ost;
+    array_ost << x;
+    std::ostringstream ref_ost;
+    ref_ost << 0 << ':' << x[0];
+    for (size_t i = 1; i < x.size(); ++i) {
+      ref_ost << ' ' << i << ':' << x[i];
+    }
+    REQUIRE(array_ost.str() == ref_ost.str());
+  }
+
 #ifndef NDEBUG
+
+  SECTION("output with signaling NaN")
+  {
+    Array<double> x{5};
+    x[0] = 2;
+    x[2] = 3;
+
+    std::ostringstream array_ost;
+    array_ost << x;
+    std::ostringstream ref_ost;
+    ref_ost << 0 << ':' << 2 << ' ';
+    ref_ost << 1 << ":nan ";
+    ref_ost << 2 << ':' << 3 << ' ';
+    ref_ost << 3 << ":nan ";
+    ref_ost << 4 << ":nan";
+    REQUIRE(array_ost.str() == ref_ost.str());
+  }
+
   SECTION("checking for bounds violation")
   {
     REQUIRE_THROWS_AS(a[10], AssertError);
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);
diff --git a/tests/test_NaNHelper.cpp b/tests/test_NaNHelper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1adbf1650f73e5592748eebad8d3534211664be7
--- /dev/null
+++ b/tests/test_NaNHelper.cpp
@@ -0,0 +1,38 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <utils/NaNHelper.hpp>
+
+#include <sstream>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("NaNHelper", "[utils]")
+{
+  SECTION("signaling NaN")
+  {
+    std::ostringstream ost_nan_helper;
+    ost_nan_helper << NaNHelper(std::numeric_limits<double>::signaling_NaN());
+
+    REQUIRE(ost_nan_helper.str() == "nan");
+  }
+
+  SECTION("valid double")
+  {
+    std::ostringstream ost_nan_helper;
+    ost_nan_helper << NaNHelper(double{3.2});
+
+    std::ostringstream ost;
+    ost << double{3.2};
+    REQUIRE(ost_nan_helper.str() == ost.str());
+  }
+
+  SECTION("non arithmetic type")
+  {
+    std::string s = "foo";
+    std::ostringstream ost_nan_helper;
+    ost_nan_helper << NaNHelper(s);
+
+    REQUIRE(ost_nan_helper.str() == s);
+  }
+}
diff --git a/tests/test_Table.cpp b/tests/test_Table.cpp
index 5372c697dc8d8bde1e2039141a8d8e226e7a68d3..a89344b75262a2eb3322589b613d6add8ae5764e 100644
--- a/tests/test_Table.cpp
+++ b/tests/test_Table.cpp
@@ -158,27 +158,62 @@ TEST_CASE("Table", "[utils]")
 
   SECTION("checking for Kokkos::View encaspulation")
   {
-    {
-      Kokkos::View<double**> kokkos_view("anonymous", 10, 3);
-      for (size_t i = 0; i < 10; ++i) {
-        for (size_t j = 0; j < 3; ++j) {
-          kokkos_view(i, j) = 3 * i + j;
-        }
+    Kokkos::View<double**> kokkos_view("anonymous", 10, 3);
+    for (size_t i = 0; i < 10; ++i) {
+      for (size_t j = 0; j < 3; ++j) {
+        kokkos_view(i, j) = 3 * i + j;
       }
+    }
 
-      Table table = encapsulate(kokkos_view);
+    Table table = encapsulate(kokkos_view);
 
-      REQUIRE(table.numberOfRows() == kokkos_view.extent(0));
-      REQUIRE(table.numberOfColumns() == kokkos_view.extent(1));
-      for (size_t i = 0; i < table.numberOfRows(); ++i) {
-        for (size_t j = 0; j < table.numberOfColumns(); ++j) {
-          REQUIRE(&table(i, j) == &kokkos_view(i, j));
-        }
+    REQUIRE(table.numberOfRows() == kokkos_view.extent(0));
+    REQUIRE(table.numberOfColumns() == kokkos_view.extent(1));
+    for (size_t i = 0; i < table.numberOfRows(); ++i) {
+      for (size_t j = 0; j < table.numberOfColumns(); ++j) {
+        REQUIRE(&table(i, j) == &kokkos_view(i, j));
       }
     }
   }
 
+  SECTION("output")
+  {
+    Table<double> table(3, 2);
+    table(0, 0) = 1;
+    table(0, 1) = 3;
+    table(1, 0) = 7;
+    table(1, 1) = -2;
+    table(2, 0) = 4;
+    table(2, 1) = -5;
+    std::ostringstream table_ost;
+    table_ost << table;
+
+    std::ostringstream ref_ost;
+    ref_ost << "0| " << 0 << ':' << 1. << ' ' << 1 << ':' << 3. << '\n';
+    ref_ost << "1| " << 0 << ':' << 7. << ' ' << 1 << ':' << -2. << '\n';
+    ref_ost << "2| " << 0 << ':' << 4. << ' ' << 1 << ':' << -5. << '\n';
+
+    REQUIRE(table_ost.str() == ref_ost.str());
+  }
+
 #ifndef NDEBUG
+  SECTION("output with signaling NaN")
+  {
+    Table<double> table(3, 2);
+    table(0, 0) = 1;
+    table(1, 1) = -2;
+    table(2, 1) = -5;
+    std::ostringstream table_ost;
+    table_ost << table;
+
+    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(table_ost.str() == ref_ost.str());
+  }
+
   SECTION("checking for bounds violation")
   {
     REQUIRE_THROWS_AS(a(4, 0), AssertError);
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 31ded1e245bc026a0a31405ddccca0413717c402..f50d2abc61f26bdbebbac4ba14ce884b352ec0aa 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -9,6 +9,8 @@
 
 #include <algebra/TinyMatrix.hpp>
 
+#include <sstream>
+
 // Instantiate to ensure full coverage is performed
 template class TinyMatrix<1, int>;
 template class TinyMatrix<2, int>;
@@ -238,6 +240,20 @@ TEST_CASE("TinyMatrix", "[algebra]")
   }
 
 #ifndef NDEBUG
+  SECTION("output with signaling NaN")
+  {
+    TinyMatrix<2> A;
+    A(0, 0) = 1;
+    A(1, 0) = 2;
+    std::ostringstream A_ost;
+    A_ost << A;
+
+    std::ostringstream ref_ost;
+    ref_ost << "[(1,nan)(2,nan)]";
+
+    REQUIRE(A_ost.str() == ref_ost.str());
+  }
+
   SECTION("checking for bounds violation")
   {
     REQUIRE_THROWS_AS(A(3, 0), AssertError);
diff --git a/tests/test_TinyVector.cpp b/tests/test_TinyVector.cpp
index c4fb06998505d0eb5e5697835958d9a970d5e01a..cf77a22325ee19137367ab68eff4f2f5a817ef6a 100644
--- a/tests/test_TinyVector.cpp
+++ b/tests/test_TinyVector.cpp
@@ -5,6 +5,8 @@
 
 #include <algebra/TinyVector.hpp>
 
+#include <sstream>
+
 // Instantiate to ensure full coverage is performed
 template class TinyVector<1, int>;
 template class TinyVector<3, int>;
@@ -81,6 +83,19 @@ TEST_CASE("TinyVector", "[algebra]")
   }
 
 #ifndef NDEBUG
+  SECTION("output with signaling NaN")
+  {
+    TinyVector<3> x;
+    x[1] = 1;
+    std::ostringstream x_ost;
+    x_ost << x;
+
+    std::ostringstream ref_ost;
+    ref_ost << "(nan,1,nan)";
+
+    REQUIRE(x_ost.str() == ref_ost.str());
+  }
+
   SECTION("checking for bounds validation")
   {
     REQUIRE_THROWS_AS(x[4] = 0, AssertError);
diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp
index 5a8f2f77ced534e9f1e2029e228dd2ed6fc372fd..e31bf166713b70f69ca108d7dace6dc07179dc5b 100644
--- a/tests/test_Vector.cpp
+++ b/tests/test_Vector.cpp
@@ -362,17 +362,18 @@ TEST_CASE("Vector", "[algebra]")
   SECTION("output")
   {
     Vector<int> x{5};
-    x[0] = 0;
-    x[1] = 1;
+    x[0] = 3;
+    x[1] = 7;
     x[2] = 2;
-    x[3] = 3;
-    x[4] = 4;
+    x[3] = 1;
+    x[4] = -4;
 
     std::ostringstream vector_ost;
     vector_ost << x;
     std::ostringstream ref_ost;
-    for (size_t i = 0; i < x.size(); ++i) {
-      ref_ost << ' ' << x[i];
+    ref_ost << 0 << ':' << x[0];
+    for (size_t i = 1; i < x.size(); ++i) {
+      ref_ost << ' ' << i << ':' << x[i];
     }
     REQUIRE(vector_ost.str() == ref_ost.str());
   }
@@ -431,6 +432,24 @@ TEST_CASE("Vector", "[algebra]")
   }
 
 #ifndef NDEBUG
+
+  SECTION("output with signaling NaN")
+  {
+    Vector<double> x{5};
+    x[0] = 3;
+    x[2] = 2;
+
+    std::ostringstream vector_ost;
+    vector_ost << x;
+    std::ostringstream ref_ost;
+    ref_ost << 0 << ':' << 3 << ' ';
+    ref_ost << 1 << ":nan ";
+    ref_ost << 2 << ':' << 2 << ' ';
+    ref_ost << 3 << ":nan ";
+    ref_ost << 4 << ":nan";
+    REQUIRE(vector_ost.str() == ref_ost.str());
+  }
+
   SECTION("invalid dot product")
   {
     Vector<int> x{5};