comment out monitor methods

examples/main.cpp

@@ -0,0 +1,392 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../data/star.mesh
+//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
+//               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
+//               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
+//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Device sample runs:
+//               mpirun -np 4 ex1p -pa -d cuda
+//               mpirun -np 4 ex1p -fa -d cuda
+//               mpirun -np 4 ex1p -pa -d occa-cuda
+//               mpirun -np 4 ex1p -pa -d raja-omp
+//               mpirun -np 4 ex1p -pa -d ceed-cpu
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
+//             * mpirun -np 4 ex1p -pa -d ceed-cuda
+//             * mpirun -np 4 ex1p -pa -d ceed-hip
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = true;
+   bool fa = false;
+   int ser_ref_levels = 2;
+   int par_ref_levels = 0;
+   const char *device_config = "cpu";
+   bool visualization = false;
+   bool algebraic_ceed = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+#ifdef MFEM_USE_CEED
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
+                  "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      for (int lev = 0; lev < ser_ref_levels; lev++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      for (int lev = 0; lev < par_ref_levels; lev++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (pmesh.GetNodes())
+   {
+      fec = pmesh.GetNodes()->OwnFEC();
+      delete_fec = false;
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+   ParFiniteElementSpace fespace(&pmesh, fec);
+   HYPRE_BigInt size = fespace.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid
+   //     function corresponding to fespace. Initialize x with initial guess of
+   //     zero, which satisfies the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   ParGridFunction x_pp(&fespace);
+   x_pp = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the
+   //     Diffusion domain integrator.
+   ParBilinearForm a(&fespace);
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (fa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::FULL);
+      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
+      // when Device::IsEnabled() returns true). This makes the results
+      // bit-for-bit deterministic at the cost of somewhat longer run time.
+      a.EnableSparseMatrixSorting(Device::IsEnabled());
+   }
+   a.AddDomainIntegrator(new MassIntegrator(one));
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a.EnableStaticCondensation(); }
+   a.Assemble();
+
+   OperatorPtr A;
+   Vector B, X, X_pp;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // 13. Solve the linear system A X = B.
+   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
+   //     * With partial assembly, use Jacobi smoothing, for now.
+   Solver *prec = NULL;
+   if (pa)
+   {
+      prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+   }
+   else
+   {
+      prec = new HypreBoomerAMG;
+   }
+
+#if 1
+   if (myid == 0)
+   {
+      std::cout<<"\n Running CGSolver"<<std::endl;
+   }
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   //cg.SetPrintLevel(1);
+   if (prec) { cg.SetPreconditioner(*prec); }
+   cg.SetOperator(*A);
+
+   tic_toc.Clear();
+   // Start & Stop CG timing.
+   {
+      tic_toc.Start();
+      cg.Mult(B, X);
+      tic_toc.Stop();
+
+      double rt_min, rt_max, my_rt;
+      my_rt = tic_toc.RealTime();
+      MPI_Reduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, 0, pmesh.GetComm());
+      MPI_Reduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, 0, pmesh.GetComm());
+
+      if (myid == 0)
+      {
+         int cg_iter = cg.GetNumIterations();
+         std::cout << "No of iterations = "<<cg_iter <<std::endl;
+         std::cout << "Total CG time:    " << rt_max << " (" << rt_min << ") sec."
+                   << std::endl;
+      }
+   }
+   //delete prec;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a.RecoverFEMSolution(X, b, x);
+
+   //std::cout<<"\n Reference solution"<<std::endl;
+   //x.Print();
+#endif
+
+   if (myid == 0)
+   {
+      std::cout<<"\n Running PipelinedPCGSolver"<<std::endl;
+   }
+   a.FormLinearSystem(ess_tdof_list, x_pp, b, A, X_pp, B);
+   PipelinedPCGSolver ppcg(MPI_COMM_WORLD);
+   ppcg.SetRelTol(1e-12);
+   ppcg.SetMaxIter(2000);
+   //ppcg.SetPrintLevel(1);
+   if (prec) { ppcg.SetPreconditioner(*prec); }
+   ppcg.SetOperator(*A);
+
+
+   tic_toc.Clear();
+   // Start & Stop CG timing.
+   {
+      tic_toc.Start();
+      ppcg.Mult(B, X_pp);
+      tic_toc.Stop();
+
+      double rt_min, rt_max, my_rt;
+      my_rt = tic_toc.RealTime();
+      MPI_Reduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, 0, pmesh.GetComm());
+      MPI_Reduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, 0, pmesh.GetComm());
+
+      if (myid == 0)
+      {
+         int ppcg_iter = ppcg.GetNumIterations();
+         std::cout << "No of iterations = "<<ppcg_iter <<std::endl;
+         std::cout << "Total CG time:    " << rt_max << " (" << rt_min << ") sec."
+                   << std::endl;
+      }
+   }
+   delete prec;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a.RecoverFEMSolution(X, b, x_pp);
+
+   //std::cout<<"\n Pipelined solution"<<std::endl;
+   //x_pp.Print();
+
+   //x.Print();
+   x_pp -= x;
+   double norm_diff = x_pp.Normlinf();
+   MPI_Allreduce(MPI_IN_PLACE, &norm_diff, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   if (myid == 0)
+   {
+      std::cout<<"\n Diff in solution = "<<norm_diff<<std::endl;
+   }
+
+
+
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   if(false) // don't write out
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   if (delete_fec)
+   {
+      delete fec;
+   }
+
+   return 0;
+}

examples/makefile

@@ -21,7 +21,7 @@ CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
+SEQ_EXAMPLES = main ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
   ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
   ex31 ex33
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \

linalg/solvers.cpp

@@ -853,6 +853,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
          add(r, beta, d, d);
       }
       oper->Mult(d, z);       //  z = A d
+
       den = Dot(d, z);
       MFEM_ASSERT(IsFinite(den), "den = " << den);
       if (den <= 0.0)
@@ -927,6 +928,227 @@ void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
 }
 
 
+void PipelinedPCGSolver::UpdateVectors()
+{
+   MemoryType mt = GetMemoryType(oper->GetMemoryClass());
+
+   r.SetSize(width, mt); r.UseDevice(true);
+   u.SetSize(width, mt); u.UseDevice(true);
+   w.SetSize(width, mt); w.UseDevice(true);
+   m.SetSize(width, mt); m.UseDevice(true);
+   n.SetSize(width, mt); n.UseDevice(true);
+
+   z.SetSize(width, mt); z.UseDevice(true);
+   q.SetSize(width, mt); q.UseDevice(true);
+   s.SetSize(width, mt); s.UseDevice(true);
+   p.SetSize(width, mt); p.UseDevice(true);
+}
+
+void PipelinedPCGSolver::Mult(const Vector &b, Vector &x) const
+{
+
+   double r0;
+   x.UseDevice(true);
+
+   oper->Mult(x, r);
+   subtract(b, r, r); // r = b - Ax
+
+   prec->Mult(r, u); // u = inv(M) r
+
+   oper->Mult(u, w); // w = A u
+
+   double gamma = 0; //intialize so we may copy the value in the loop
+   double gamma_0, gamma_old;
+   double delta;
+
+   double beta, alpha;
+   double alpha_old;
+
+   MPI_Request request;
+
+   Vector loc_gamma_delta(2);
+   Vector glo_gamma_delta(2);
+
+   bool converged = false;
+   final_iter = max_iter;
+
+   for (int i = 0; true;) //has to start at zero
+   {
+      //Fuse these two operations below
+      gamma_old = gamma;
+
+      //gamma = Dot(r, u);
+      //delta = Dot(w, u);
+      loc_gamma_delta(0) = r * u; //local inner product
+      loc_gamma_delta(1) = w * u; //local inner product
+
+      MPI_Iallreduce(loc_gamma_delta.HostRead(), glo_gamma_delta.HostWrite(), 2, MPI_DOUBLE, MPI_SUM, GetComm(), &request);
+
+      //Do products computation may be overlapped with the following below
+      prec->Mult(w, m);
+      oper->Mult(m, n);
+
+      MPI_Wait(&request, MPI_STATUS_IGNORE);
+
+      gamma = glo_gamma_delta(0);
+      delta = glo_gamma_delta(1);
+
+
+      if (i == 0)
+      {
+         gamma_0 = initial_norm = gamma;
+         if (gamma >= 0.0) { initial_norm = sqrt(gamma);}
+         if (print_options.iterations || print_options.first_and_last)
+         {
+            mfem::out << "   Iteration : " << setw(3) << 0 << "  (B r, r) = "
+                      << gamma << (print_options.first_and_last ? " ...\n" : "\n");
+         }
+         //Monitor(0, gamma, r, x);
+
+         if (gamma < 0.0)
+         {
+            if (print_options.warnings)
+            {
+               mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
+                         << gamma << '\n';
+            }
+            converged = false;
+            final_iter = 0;
+            initial_norm = gamma;
+            final_norm = gamma;
+            return;
+         }
+         r0 = std::max(gamma*rel_tol*rel_tol, abs_tol*abs_tol);
+         if (gamma <= r0)
+         {
+            converged = true;
+            final_iter = 0;
+            final_norm = sqrt(gamma);
+            return;
+         }
+
+         MFEM_ASSERT(IsFinite(delta), "delta = " << delta);
+         if (delta <= 0.0)
+         {
+           /*
+            if (Dot(w, w) > 0.0 && print_options.warnings)
+            {
+               mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
+                         << delta << '\n';
+            }
+           */
+            if (delta == 0.0)
+            {
+               converged = false;
+               final_iter = 0;
+               final_norm = sqrt(gamma);
+               return;
+            }
+         }
+      }// i == 0
+      else
+      {
+
+         //Check if preconditioner is positive definite
+         MFEM_ASSERT(IsFinite(gamma), "gamma = " << gamma);
+         if (gamma < 0.0)
+         {
+            if (print_options.warnings)
+            {
+               mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
+                         << gamma << '\n';
+            }
+            converged = false;
+            final_iter = i;
+            break;
+         }
+
+         //Report iterations
+         if (print_options.iterations)
+         {
+            mfem::out << "   Iteration : " << setw(3) << i << "  (B r, r) = "
+                      << gamma << std::endl;
+         }
+
+         //Monitor(i, gamma, r, x);
+
+         if (gamma <= r0)
+         {
+            converged = true;
+            final_iter = i;
+            break;
+         }
+
+      } // i > 0
+
+      if ( i > 0 )
+      {
+         alpha_old = alpha;
+         beta = gamma/gamma_old;
+         alpha = gamma/(delta - beta*gamma/alpha_old);
+      }
+      else
+      {
+         beta = 0;
+         alpha = gamma/delta;
+      }
+
+      //If maxed out on iterations break
+      if (++i > max_iter)
+      {
+         break;
+      }
+
+      //z_i = n_i + beta_i * z_{i-1}
+      add(n, beta, z, z);
+
+      //q_i = m_i + beta_i q_{i-1}
+      add(m, beta, q, q);
+
+      //s_i = w_i + beta_i s_{i-1}
+      add(w, beta, s, s);
+
+      //p_i = u_i + beta_i p_{i-1}
+      add(u, beta, p, p);
+
+      //x_{i+1} = x_i + alpha_i p_i
+      add(x, alpha, p, x);
+
+      //r_{i+1} = r_i - alpha_i s_i
+      add(r, (-alpha), s, r);
+
+      //u_{i+1} = u_i  - alpha_i q_i
+      add(u, (-alpha), q, u);
+
+      //w_{i+1} = w_i - alpha_i z_i
+      add(w, (-alpha), z, w);
+   }
+
+   if (print_options.first_and_last && !print_options.iterations)
+   {
+      mfem::out << "   Iteration : " << setw(3) << final_iter << "  (B r, r) = "
+                << gamma << '\n';
+   }
+   if (print_options.summary || (print_options.warnings && !converged))
+   {
+      mfem::out << "PCG: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_options.summary || print_options.iterations ||
+       print_options.first_and_last)
+   {
+      const auto arf = pow (gamma/gamma_0, 0.5/final_iter);
+      mfem::out << "Average reduction factor = " << arf << '\n';
+   }
+   if (print_options.warnings && !converged)
+   {
+      mfem::out << "PCG: No convergence!" << '\n';
+   }
+
+   final_norm = sqrt(gamma);
+
+   ///Monitor(final_iter, final_norm, r, x, true);
+}
+
 inline void GeneratePlaneRotation(double &dx, double &dy,
                                   double &cs, double &sn)
 {

linalg/solvers.hpp

@@ -521,6 +521,30 @@ void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
          double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
 
 
+/// Preconditioned Conjugate gradient method
+/// Algorithm 4. Preconditioned pipelined CG
+/// https://www.sciencedirect.com/science/article/pii/S0167819113000719?via%3Dihub
+class PipelinedPCGSolver : public IterativeSolver
+{
+protected:
+   mutable Vector r, u, w, m, n;
+   mutable Vector z, q, s, p;
+
+   void UpdateVectors();
+
+public:
+   PipelinedPCGSolver() { }
+
+#ifdef MFEM_USE_MPI
+   PipelinedPCGSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
+#endif
+
+   virtual void SetOperator(const Operator &op)
+   { IterativeSolver::SetOperator(op); UpdateVectors(); }
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
 /// GMRES method
 class GMRESSolver : public IterativeSolver
 {