diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..357ab3a1d1275d6f79ceb82627e28f0003f06edd
--- /dev/null
+++ b/src/algebra/PCG.hpp
@@ -0,0 +1,94 @@
+#ifndef PCG_HPP
+#define PCG_HPP
+
+struct PCG
+{
+  template <typename VectorType, typename MatrixType, typename PreconditionerType>
+  PCG(const VectorType& f,
+      const MatrixType& A,
+      [[maybe_unused]] const PreconditionerType& C,
+      VectorType& x,
+      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';
+
+    VectorType h{f.size()};
+    VectorType b = copy(f);
+
+    if (true) {
+      h = A * x;
+      h -= f;
+      std::cout << "- initial *real* residu :   " << (h, h) << '\n';
+    }
+
+    VectorType g{b.size()};
+    VectorType cg = copy(b);
+
+    double gcg  = 0;
+    double gcg0 = 1;
+
+    double relativeEpsilon = epsilon;
+
+    for (size_t i = 1; i <= maxiter; ++i) {
+      if (i == 1) {
+        h = A * x;
+
+        cg -= h;
+
+        g = copy(cg);   // C.computes(cg, g);
+
+        gcg = (g, cg);
+
+        h = copy(g);
+      }
+
+      b = A * h;
+
+      double hAh = (h, b);
+
+      if (hAh == 0) {
+        hAh = 1.;
+      }
+      double ro = gcg / hAh;
+      cg -= ro * b;
+
+      // VectorType b2 = b;
+      // C.computes(b2, b);
+
+      x += ro * h;
+      g -= ro * b;
+
+      double gamma = gcg;
+      gcg          = (g, cg);
+
+      if ((i == 1) && (gcg != 0)) {
+        relativeEpsilon = epsilon * gcg;
+        gcg0            = gcg;
+        std::cout << "  initial residu: " << gcg << '\n';
+      }
+      std::cout << "  - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
+      std::cout << "\tabsolute: " << gcg << '\n';
+
+      if (gcg < relativeEpsilon) {
+        break;
+      }
+
+      gamma = gcg / gamma;
+
+      h *= gamma;
+      h += g;
+    }
+
+    if (gcg > relativeEpsilon) {
+      std::cout << "  conjugate gradient: *NOT CONVERGED*\n";
+      std::cout << "  - epsilon:          " << epsilon << '\n';
+      std::cout << "  - relative residu : " << gcg / gcg0 << '\n';
+      std::cout << "  - absolute residu : " << gcg << '\n';
+    }
+  }
+};
+
+#endif   // PCG_HPP
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index bed915b78d443b85bc77c7540fb2ae8b7df61534..5f6cff9415ff431123070ac5f2363adb328365df 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -20,6 +20,16 @@ class Vector
   static_assert(std::is_same_v<typename decltype(m_values)::index_type, index_type>);
 
  public:
+  friend PUGS_INLINE Vector
+  copy(const Vector& x)
+  {
+    auto values = copy(x.m_values);
+
+    Vector x_copy{0};
+    x_copy.m_values = values;
+    return x_copy;
+  }
+
   friend std::ostream&
   operator<<(std::ostream& os, const Vector<DataType>& x)
   {
@@ -29,6 +39,12 @@ class Vector
     return os;
   }
 
+  friend Vector operator*(const DataType& a, const Vector& x)
+  {
+    Vector y = copy(x);
+    return y *= a;
+  }
+
   PUGS_INLINE
   DataType
   operator,(const Vector& y)
@@ -37,7 +53,7 @@ class Vector
 
     // Does not use parallel_for to preserve sum order
     for (index_type i = 0; i < this->size(); ++i) {
-      sum += m_values[i] + y.m_values[i];
+      sum += m_values[i] * y.m_values[i];
     }
 
     return sum;
@@ -102,6 +118,7 @@ class Vector
   operator=(const DataType& value) noexcept
   {
     m_values.fill(value);
+    return *this;
   }
 
   template <typename DataType2>
diff --git a/src/main.cpp b/src/main.cpp
index bcd9161bd60990975229bbb93caaa26860f91d5b..760164fcf200da6f1f70183ba572d26fdea8b38e 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -25,6 +25,7 @@
 #include <SynchronizerManager.hpp>
 
 #include <CRSMatrix.hpp>
+#include <PCG.hpp>
 
 #include <iostream>
 
@@ -44,9 +45,10 @@ main(int argc, char* argv[])
       FlexibleMatrix F(matrix_size);
 
       for (size_t i = 0; i < matrix_size; ++i) {
-        F(i, i) = 2;
+        F(i, i) = 4;
         if (i + 1 < matrix_size) {
-          F(i, i + 1) = 1;
+          F(i, i + 1) = -1;
+          F(i + 1, i) = -1;
         }
       }
       std::cout << F << '\n';
@@ -64,6 +66,12 @@ main(int argc, char* argv[])
 
       Vector<double> Ax = A * x;
       std::cout << "Ax=\n" << Ax << '\n';
+
+      Vector<double> u{x.size()};
+      u = 0;
+
+      PCG{Ax, A, A, u, 9, 1e-4};
+      std::cout << "u=" << u << '\n';
     }
 
     Kokkos::finalize();