diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index c931385b214c91ca2807db4f07fca027ebd79b25..bc822668a7f63c4d5d7cc4c45f4d1b1d70084fc0 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -5,14 +5,19 @@
 #include <iomanip>
 #include <iostream>
 
+#include <rang.hpp>
+
+template <bool verbose = true>
 struct BiCGStab
 {
   template <typename VectorType, typename MatrixType>
   BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
   {
-    std::cout << "- bi-conjugate gradient stabilized\n";
-    std::cout << "  epsilon = " << epsilon << '\n';
-    std::cout << "  maximum number of iterations: " << max_iter << '\n';
+    if constexpr (verbose) {
+      std::cout << "- bi-conjugate gradient stabilized\n";
+      std::cout << "  epsilon = " << epsilon << '\n';
+      std::cout << "  maximum number of iterations: " << max_iter << '\n';
+    }
 
     VectorType r_k_1{b.size()};
 
@@ -33,10 +38,14 @@ struct BiCGStab
 
       VectorType r_k{x.size()};
 
-      std::cout << "   initial residu: " << resid0 << '\n';
+      if constexpr (verbose) {
+        std::cout << "   initial residu: " << resid0 << '\n';
+      }
       for (size_t i = 1; i <= max_iter; ++i) {
-        std::cout << "  - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
-                  << "\tabsolute: " << residu << '\n';
+        if constexpr (verbose) {
+          std::cout << "  - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
+                    << "\tabsolute: " << residu << '\n';
+        }
 
         Ap_k = A * p_k;
 
@@ -64,14 +73,14 @@ struct BiCGStab
       }
 
       if (residu / resid0 > epsilon) {
-        std::cout << "  bi_conjugate gradient stabilized: *NOT CONVERGED*\n";
+        std::cout << "  bi-conjugate gradient stabilized: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset
+                  << '\n';
+        ;
         std::cout << "  - epsilon:          " << epsilon << '\n';
         std::cout << "  - relative residu : " << residu / resid0 << '\n';
         std::cout << "  - absolute residu : " << residu << '\n';
       }
     }
-
-    std::cout << "- bi-conjugate gradient stabilized: done\n";
   }
 };
 
diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp
index 74c17ee40971b96c512c196d630f1edfadfa6ddc..36a4635b313d2bba13d97e1b7a4ea1882e2a66a2 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/PCG.hpp
@@ -4,6 +4,9 @@
 #include <iomanip>
 #include <iostream>
 
+#include <rang.hpp>
+
+template <bool verbose = true>
 struct PCG
 {
   template <typename VectorType, typename MatrixType, typename PreconditionerType>
@@ -14,14 +17,16 @@ struct PCG
       const size_t maxiter,
       const double epsilon = 1e-6)
   {
-    std::cout << "- conjugate gradient\n";
-    std::cout << "  epsilon = " << epsilon << '\n';
-    std::cout << "  maximum number of iterations: " << maxiter << '\n';
+    if constexpr (verbose) {
+      std::cout << "- conjugate gradient\n";
+      std::cout << "  epsilon = " << epsilon << '\n';
+      std::cout << "  maximum number of iterations: " << maxiter << '\n';
+    }
 
     VectorType h{f.size()};
     VectorType b = copy(f);
 
-    if (true) {
+    if constexpr (verbose) {
       h = A * x;
       h -= f;
       std::cout << "- initial *real* residu :   " << (h, h) << '\n';
@@ -69,10 +74,14 @@ struct PCG
       if ((i == 1) && (gcg != 0)) {
         relativeEpsilon = epsilon * gcg;
         gcg0            = gcg;
-        std::cout << "  initial residu: " << gcg << '\n';
+        if constexpr (verbose) {
+          std::cout << "  initial residu: " << gcg << '\n';
+        }
+      }
+      if constexpr (verbose) {
+        std::cout << "  - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
+        std::cout << "\tabsolute: " << gcg << '\n';
       }
-      std::cout << "  - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
-      std::cout << "\tabsolute: " << gcg << '\n';
 
       if (gcg < relativeEpsilon) {
         break;
@@ -85,7 +94,7 @@ struct PCG
     }
 
     if (gcg > relativeEpsilon) {
-      std::cout << "  conjugate gradient: *NOT CONVERGED*\n";
+      std::cout << "  conjugate gradient: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset << '\n';
       std::cout << "  - epsilon:          " << epsilon << '\n';
       std::cout << "  - relative residu : " << gcg / gcg0 << '\n';
       std::cout << "  - absolute residu : " << gcg << '\n';
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index b7cb2d3bd8afe280a7d8c416a662016851e0d820..99b94515ba7795903434db761d509b9052fcaf71 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    BiCGStab{b, A, x, 10, 1e-12};
+    BiCGStab<false>{b, A, x, 10, 1e-12};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
   }
diff --git a/tests/test_PCG.cpp b/tests/test_PCG.cpp
index 4b2f2df9bed4bec0ea31502be3c5480e1eec6989..d4564761d8c564082e4668360e5a2e2a8ecce008 100644
--- a/tests/test_PCG.cpp
+++ b/tests/test_PCG.cpp
@@ -62,7 +62,7 @@ TEST_CASE("PCG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    PCG{b, A, A, x, 10, 1e-12};
+    PCG<false>{b, A, A, x, 10, 1e-12};
     REQUIRE(std::sqrt((x, x)) == 0);
   }
 
@@ -104,7 +104,7 @@ TEST_CASE("PCG", "[algebra]")
     Vector<double> x{5};
     x = 0;
 
-    PCG{b, A, A, x, 1, 1e-12};
+    PCG<false>{b, A, A, x, 1, 1e-12};
     Vector error = x - x_exact;
     REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
   }