Add BiCGStab and make few clean-up

src/algebra/BiCGStab.hpp

+#ifndef BICG_STAB_HPP
+#define BICG_STAB_HPP
+
+#include <iostream>
+
+class BiCGStab
+{
+ public:
+  template <typename VectorType, typename MatrixType, typename PreconditionerType>
+  BiCGStab(const VectorType& b,
+           const MatrixType& A,
+           [[maybe_unused]] const PreconditionerType& P,
+           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';
+
+    VectorType r_k_1{b.size()};
+
+    r_k_1 = b - A * x;
+
+    double residu = std::sqrt((r_k_1, r_k_1));   // Norm(r_k_1);
+
+    if (residu != 0) {
+      double resid0 = residu;
+
+      VectorType rTilda_0 = copy(r_k_1);
+      VectorType p_k      = copy(r_k_1);
+
+      VectorType s_k{x.size()};
+
+      VectorType Ap_k{x.size()};
+      VectorType As_k{x.size()};
+
+      VectorType r_k{x.size()};
+
+      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';
+
+        Ap_k = A * p_k;
+
+        const double alpha_k = (r_k_1, rTilda_0) / (Ap_k, rTilda_0);
+
+        s_k  = r_k_1 - alpha_k * Ap_k;
+        As_k = A * s_k;
+
+        const double w_k = (As_k, s_k) / (As_k, As_k);
+
+        x += alpha_k * p_k + w_k * s_k;
+        r_k = s_k - w_k * As_k;
+
+        const double beta_k = (r_k, rTilda_0) / (r_k_1, rTilda_0) * (alpha_k / w_k);
+
+        p_k -= w_k * Ap_k;
+        p_k *= beta_k;
+        p_k += r_k;
+
+        if ((residu = std::sqrt((r_k, r_k))) / resid0 < epsilon) {
+          break;
+        }
+
+        r_k_1 = r_k;
+      }
+
+      if (residu / resid0 > epsilon) {
+        std::cout << "  bi_conjugate gradient stabilized: *NOT CONVERGED*\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";
+  }
+};
+
+#endif   // BICG_STAB_HPP

src/algebra/PCG.hpp

@@ -38,7 +38,7 @@ struct PCG
 
         cg -= h;
 
-        g = copy(cg);   // C.computes(cg, g);
+        g = copy(cg);   // TODO: precond: g = g/C
 
         gcg = (g, cg);
 
@@ -55,8 +55,7 @@ struct PCG
       double ro = gcg / hAh;
       cg -= ro * b;
 
-      // VectorType b2 = b;
-      // C.computes(b2, b);
+      // TODO: precond: b <- b/C
 
       x += ro * h;
       g -= ro * b;

src/algebra/Vector.hpp

@@ -103,7 +103,7 @@ class Vector
 
   PUGS_INLINE
   Vector
-  operator+(const Vector& y)
+  operator+(const Vector& y) const
   {
     Assert(this->size() == y.size());
     Vector sum{y.size()};
@@ -116,7 +116,7 @@ class Vector
 
   PUGS_INLINE
   Vector
-  operator-(const Vector& y)
+  operator-(const Vector& y) const
   {
     Assert(this->size() == y.size());
     Vector sum{y.size()};

src/main.cpp

@@ -24,6 +24,7 @@
 
 #include <SynchronizerManager.hpp>
 
+#include <BiCGStab.hpp>
 #include <CRSMatrix.hpp>
 #include <PCG.hpp>
 
@@ -41,6 +42,7 @@ main(int argc, char* argv[])
     Kokkos::initialize({4, -1, -1, true});
 
     {
+      {   // PCG simple test
         size_t matrix_size = 10;
         FlexibleMatrix F(matrix_size);
 
@@ -51,31 +53,59 @@ main(int argc, char* argv[])
             F(i + 1, i) = -1;
           }
         }
-      std::cout << F << '\n';
 
         CRSMatrix A{F};
 
-      std::cout << A << '\n';
-
         Vector<double> x{10};
         for (size_t i = 0; i < x.size(); ++i) {
           x[i] = 2 * i + 1;
         }
 
-      std::cout << "x=\n" << x << '\n';
-
         Vector<double> Ax = A * x;
-      std::cout << "Ax=\n" << Ax << '\n';
 
         Vector<double> u{x.size()};
         u = 0;
 
-      PCG{Ax, A, A, u, static_cast<size_t>(-1), 1e-12};
-      std::cout << "u=\n" << u << '\n';
+        PCG{Ax, A, A, u, 100, 1e-12};
+
+        std::cout << "== PCG report ==\n";
+        std::cout << "L2 Error = " << std::sqrt((x - u, x - u)) << '\n';
+        std::cout << "L2 Norm  = " << std::sqrt((x, x)) << '\n';
+        std::cout << "L2 RelEr = " << std::sqrt((x - u, x - u)) / std::sqrt((x, x)) << '\n';
+        std::cout << "================\n";
+      }
+
+      {   // BiCGStab simple test
+        size_t matrix_size = 10;
+        FlexibleMatrix F(matrix_size);
+
+        for (size_t i = 0; i < matrix_size; ++i) {
+          F(i, i) = 4;
+          if (i + 1 < matrix_size) {
+            F(i, i + 1) = -1;
+            F(i + 1, i) = -0.3;
+          }
+        }
+
+        CRSMatrix A{F};
+
+        Vector<double> x{10};
+        for (size_t i = 0; i < x.size(); ++i) {
+          x[i] = 2 * i + 1;
+        }
+
+        Vector<double> Ax = A * x;
+        Vector<double> u{x.size()};
+        u = 0;
+
+        BiCGStab{Ax, A, A, u, 100, 1e-12};
 
+        std::cout << "== BiCGStab report ==\n";
         std::cout << "L2 Error = " << std::sqrt((x - u, x - u)) << '\n';
         std::cout << "L2 Norm  = " << std::sqrt((x, x)) << '\n';
         std::cout << "L2 RelEr = " << std::sqrt((x - u, x - u)) / std::sqrt((x, x)) << '\n';
+        std::cout << "=====================\n";
+      }
     }
 
     Kokkos::finalize();