Merge branch 'feature/solvers-verbosity' into 'develop'

Add a static switch to choose verbosity of PCG/BiCGStab

See merge request !48

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)
   {
+    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()};
 
+      if constexpr (verbose) {
         std::cout << "   initial residu: " << resid0 << '\n';
+      }
       for (size_t i = 1; i <= max_iter; ++i) {
+        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";
   }
 };
 

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)
   {
+    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;
+        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';
+      }
 
       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';

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)));
   }

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)));
   }