Add more annotations, modifying LOR solvers miniapp for benchmarking

MPI-synchronizations at the start and at the end of annotation regions. These
synchronizations can be enabled or disabled (default) in code via the new
macros: MFEM_PERF_SYNC, MFEM_PERF_SYNC_STREAM, and MFEM_PERF_SYNC_MPI; the
environment variables with the same names can be set to 0/1 to control the
synchronization as well.

.gitignore

@@ -409,6 +409,8 @@ miniapps/tribol/contact-patch-test
 miniapps/diag-smoothers/abs-l1-jacobi
 miniapps/diag-smoothers/mg-abs-l1-jacobi
 
+miniapps/benchmarks/ceed-solver-bps/solver-bp
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests
@@ -429,6 +431,7 @@ tests/benchmarks/bench_ceed
 tests/benchmarks/bench_tmop
 tests/benchmarks/bench_vector
 tests/benchmarks/bench_virtuals
+tests/benchmarks/pbench_ceed
 
 # Test script output
 tests/scripts/*.err

CHANGELOG

@@ -70,6 +70,13 @@ Miscellaneous
   variable is an alternative to calling 'Device::SetGPUAwareMPI(true)'.
 - Added parallel Address Sanitizer, serial and parallel Undefined Behavior
   Sanitizer and serial Memory Sanitizer GitHub actions tests on Ubuntu.
+- MFEM_PERF_* annotations: added options to enable GPU-stream- and
+  MPI-synchronizations at the start and at the end of annotation regions. These
+  synchronizations can be enabled or disabled (default) in code via the new
+  macros: MFEM_PERF_SYNC, MFEM_PERF_SYNC_STREAM, and MFEM_PERF_SYNC_MPI; the
+  environment variables with the same names can be set to 0/1 to control the
+  synchronization as well.
+
 
 Version 4.8, released on Apr 9, 2025
 ====================================

examples/ex26.cpp

@@ -78,6 +78,7 @@ private:
       opr.SetOperatorOwner(false);
 
       CGSolver* pcg = new CGSolver();
+      // pcg->iterative_mode = false; // the multigrid algorithm does this
       pcg->SetPrintLevel(-1);
       pcg->SetMaxIter(200);
       pcg->SetRelTol(sqrt(1e-4));
@@ -100,7 +101,8 @@ private:
       Vector diag(fespace.GetTrueVSize());
       bfs[level]->AssembleDiagonal(diag);
 
-      Solver* smoother = new OperatorChebyshevSmoother(*opr, diag, ess_tdof_list, 2);
+      Solver *smoother = new OperatorChebyshevSmoother(
+         *opr, diag, ess_tdof_list, 2);
       AddLevel(opr.Ptr(), smoother, true, true);
    }
 };

examples/ex26p.cpp

@@ -88,6 +88,7 @@ private:
       amg->SetPrintLevel(-1);
 
       CGSolver* pcg = new CGSolver(MPI_COMM_WORLD);
+      // pcg->iterative_mode = false; // the multigrid algorithm does this
       pcg->SetPrintLevel(-1);
       pcg->SetMaxIter(10);
       pcg->SetRelTol(sqrt(1e-4));

fem/bilinearform_ext.cpp

@@ -255,6 +255,8 @@ PABilinearFormExtension::PABilinearFormExtension(BilinearForm *form)
 
 void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 {
+   MFEM_PERF_FUNCTION;
+
    if ( Device::Allows(Backend::CEED_MASK) ) { return; }
    ElementDofOrdering ordering = GetEVectorOrdering(*a->FESpace());
    elem_restrict = trial_fes->GetElementRestriction(ordering);
@@ -373,6 +375,8 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 
 void PABilinearFormExtension::Assemble()
 {
+   MFEM_PERF_FUNCTION;
+
    SetupRestrictionOperators(L2FaceValues::DoubleValued);
 
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
@@ -529,6 +533,8 @@ void PABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
 void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
                                            const bool useAbs) const
 {
+   MFEM_PERF_FUNCTION;
+
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
 
    const int iSz = integrators.Size();

fem/integ/bilininteg_diffusion_kernels.hpp

@@ -1212,22 +1212,18 @@ inline void SmemPADiffusionApply3D(const int NE,
 
 } // namespace internal
 
-namespace
-{
-using ApplyKernelType = DiffusionIntegrator::ApplyKernelType;
-using DiagonalKernelType = DiffusionIntegrator::DiagonalKernelType;
-}
 
 template<int DIM, int T_D1D, int T_Q1D>
-ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Kernel()
+DiffusionIntegrator::ApplyKernelType
+DiffusionIntegrator::ApplyPAKernels::Kernel()
 {
    if (DIM == 2) { return internal::SmemPADiffusionApply2D<T_D1D,T_Q1D>; }
    else if (DIM == 3) { return internal::SmemPADiffusionApply3D<T_D1D, T_Q1D>; }
    else { MFEM_ABORT(""); }
 }
 
-inline
-ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
+inline DiffusionIntegrator::ApplyKernelType
+DiffusionIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
 {
    if (DIM == 2) { return internal::PADiffusionApply2D; }
    else if (DIM == 3) { return internal::PADiffusionApply3D; }
@@ -1235,14 +1231,15 @@ ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
 }
 
 template<int DIM, int D1D, int Q1D>
-DiagonalKernelType DiffusionIntegrator::DiagonalPAKernels::Kernel()
+DiffusionIntegrator::DiagonalKernelType
+DiffusionIntegrator::DiagonalPAKernels::Kernel()
 {
    if (DIM == 2) { return internal::SmemPADiffusionDiagonal2D<D1D,Q1D>; }
    else if (DIM == 3) { return internal::SmemPADiffusionDiagonal3D<D1D, Q1D>; }
    else { MFEM_ABORT(""); }
 }
 
-inline DiagonalKernelType
+inline DiffusionIntegrator::DiagonalKernelType
 DiffusionIntegrator::DiagonalPAKernels::Fallback(int DIM, int, int)
 {
    if (DIM == 2) { return internal::PADiffusionDiagonal2D; }

fem/integ/bilininteg_diffusion_pa.cpp

@@ -39,6 +39,8 @@ void DiffusionIntegrator::AssembleDiagonalPA(Vector &diag)
 // PA Diffusion Apply kernel
 void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
    if (DeviceCanUseCeed())
    {
       ceedOp->AddMult(x, y);
@@ -88,6 +90,8 @@ void DiffusionIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
 
 void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
+   MFEM_PERF_FUNCTION;
+
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
    // Assuming the same element type

fem/integ/bilininteg_mass_kernels.hpp

@@ -1383,14 +1383,9 @@ inline void EAMassAssemble3D(const int NE,
 
 } // namespace internal
 
-namespace
-{
-using ApplyKernelType = MassIntegrator::ApplyKernelType;
-using DiagonalKernelType = MassIntegrator::DiagonalKernelType;
-}
 
 template<int DIM, int T_D1D, int T_Q1D>
-ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
+MassIntegrator::ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
 {
    if (DIM == 1) { return internal::PAMassApply1D; }
    else if (DIM == 2) { return internal::SmemPAMassApply2D<T_D1D,T_Q1D>; }
@@ -1398,8 +1393,8 @@ ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
    else { MFEM_ABORT(""); }
 }
 
-inline ApplyKernelType MassIntegrator::ApplyPAKernels::Fallback(
-   int DIM, int, int)
+inline MassIntegrator::ApplyKernelType
+MassIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
 {
    if (DIM == 1) { return internal::PAMassApply1D; }
    else if (DIM == 2) { return internal::PAMassApply2D; }
@@ -1408,7 +1403,7 @@ inline ApplyKernelType MassIntegrator::ApplyPAKernels::Fallback(
 }
 
 template<int DIM, int T_D1D, int T_Q1D>
-DiagonalKernelType MassIntegrator::DiagonalPAKernels::Kernel()
+MassIntegrator::DiagonalKernelType MassIntegrator::DiagonalPAKernels::Kernel()
 {
    if (DIM == 1) { return internal::PAMassAssembleDiagonal1D; }
    else if (DIM == 2) { return internal::SmemPAMassAssembleDiagonal2D<T_D1D,T_Q1D>; }
@@ -1416,8 +1411,8 @@ DiagonalKernelType MassIntegrator::DiagonalPAKernels::Kernel()
    else { MFEM_ABORT(""); }
 }
 
-inline DiagonalKernelType MassIntegrator::DiagonalPAKernels::Fallback(
-   int DIM, int, int)
+inline MassIntegrator::DiagonalKernelType
+MassIntegrator::DiagonalPAKernels::Fallback(int DIM, int, int)
 {
    if (DIM == 1) { return internal::PAMassAssembleDiagonal1D; }
    else if (DIM == 2) { return internal::PAMassAssembleDiagonal2D; }

fem/integ/bilininteg_mass_pa.cpp

@@ -23,6 +23,8 @@ namespace mfem
 
 void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
+   MFEM_PERF_FUNCTION;
+
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
 
@@ -170,6 +172,8 @@ void MassIntegrator::AssembleDiagonalPA(Vector &diag)
 
 void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
    if (DeviceCanUseCeed())
    {
       ceedOp->AddMult(x, y);

fem/linearform_ext.hpp

@@ -14,6 +14,7 @@
 
 #include "../general/array.hpp"
 #include "../linalg/vector.hpp"
+#include "fespace.hpp"
 
 namespace mfem
 {

fem/lor/lor.cpp

@@ -365,6 +365,8 @@ FiniteElementSpace &LORBase::GetFESpace() const
 
 void LORBase::AssembleSystem(BilinearForm &a_ho, const Array<int> &ess_dofs)
 {
+   MFEM_PERF_FUNCTION;
+
    A.Clear();
    delete a;
    if (BatchedLORAssembly::FormIsSupported(a_ho))

fem/lor/lor_batched.cpp

@@ -360,6 +360,8 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
 
 void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const
 {
+   MFEM_PERF_FUNCTION;
+
    const int nvdof = fes_ho.GetVSize();
 
    // If A contains an existing SparseMatrix, reuse it (and try to reuse its
@@ -417,6 +419,8 @@ static void Assemble_(LOR_KERNEL &kernel, int dim, int sdim, int order)
 template <typename LOR_KERNEL>
 void BatchedLORAssembly::AssemblyKernel(BilinearForm &a)
 {
+   MFEM_PERF_FUNCTION;
+
    LOR_KERNEL kernel(a, fes_ho, X_vert, sparse_ij, sparse_mapping);
 
    const int dim = fes_ho.GetMesh()->Dimension();

fem/lor/lor_h1_impl.hpp

@@ -184,6 +184,8 @@ void BatchedLOR_H1::Assemble2D()
 template <int ORDER>
 void BatchedLOR_H1::Assemble3D()
 {
+   MFEM_PERF_FUNCTION;
+
    const int nel_ho = fes_ho.GetNE();
    static constexpr int nv = 8;
    static constexpr int dim = 3;

fem/multigrid.cpp

@@ -17,7 +17,10 @@ namespace mfem
 MultigridBase::MultigridBase()
    : cycleType(CycleType::VCYCLE), preSmoothingSteps(1), postSmoothingSteps(1),
      nrhs(0)
-{}
+{
+   coarse_solver = nullptr;
+   own_coarse_solver = false;
+}
 
 MultigridBase::MultigridBase(const Array<Operator*>& operators_,
                              const Array<Solver*>& smoothers_,
@@ -29,12 +32,18 @@ MultigridBase::MultigridBase(const Array<Operator*>& operators_,
 {
    operators_.Copy(operators);
    smoothers_.Copy(smoothers);
+   coarse_solver = nullptr;
    ownedOperators_.Copy(ownedOperators);
    ownedSmoothers_.Copy(ownedSmoothers);
+   own_coarse_solver = false;
 }
 
 MultigridBase::~MultigridBase()
 {
+   if (own_coarse_solver)
+   {
+      delete coarse_solver;
+   }
    for (int i = 0; i < operators.Size(); ++i)
    {
       if (ownedOperators[i])
@@ -56,16 +65,17 @@ void MultigridBase::InitVectors() const
    X.SetSize(M, nrhs);
    Y.SetSize(M, nrhs);
    R.SetSize(M, nrhs);
-   Z.SetSize(M, nrhs);
-   for (int i = 0; i < X.NumRows(); ++i)
+   for (int i = 0; i < M; ++i)
    {
       const int n = operators[i]->Height();
-      for (int j = 0; j < X.NumCols(); ++j)
+      for (int j = 0; j < nrhs; ++j)
       {
-         X(i, j) = new Vector(n);
-         Y(i, j) = new Vector(n);
+         if (i < M - 1)
+         {
+            X(i, j) = new Vector(n);
+            Y(i, j) = new Vector(n);
+         }
          R(i, j) = new Vector(n);
-         Z(i, j) = new Vector(n);
       }
    }
 }
@@ -76,10 +86,12 @@ void MultigridBase::EraseVectors() const
    {
       for (int j = 0; j < X.NumCols(); ++j)
       {
-         delete X(i, j);
-         delete Y(i, j);
+         if (i < X.NumRows() - 1)
+         {
+            delete X(i, j);
+            delete Y(i, j);
+         }
          delete R(i, j);
-         delete Z(i, j);
       }
    }
 }
@@ -95,6 +107,12 @@ void MultigridBase::AddLevel(Operator* op, Solver* smoother,
    ownedSmoothers.Append(ownSmoother);
 }
 
+void MultigridBase::AddCoarseSolver(Solver *c_solver, bool own_c_solver)
+{
+   coarse_solver = c_solver;
+   own_coarse_solver = own_c_solver;
+}
+
 void MultigridBase::SetCycleType(CycleType cycleType_, int preSmoothingSteps_,
                                  int postSmoothingSteps_)
 {
@@ -105,10 +123,12 @@ void MultigridBase::SetCycleType(CycleType cycleType_, int preSmoothingSteps_,
 
 void MultigridBase::Mult(const Vector& x, Vector& y) const
 {
-   Array<const Vector*> X_(1);
-   Array<Vector*> Y_(1);
-   X_[0] = &x;
-   Y_[0] = &y;
+   const Vector *x_array[1] = { &x };
+   Array<const Vector*> X_(x_array, 1); // no heap allocation
+
+   Vector *y_array[1] = { &y };
+   Array<Vector*> Y_(y_array, 1); // no heap allocation
+
    ArrayMult(X_, Y_);
 }
 
@@ -119,11 +139,6 @@ void MultigridBase::ArrayMult(const Array<const Vector*>& X_,
                "Multigrid solver does not have operators set!");
    MFEM_ASSERT(X_.Size() == Y_.Size(),
                "Number of columns mismatch in MultigridBase::Mult!");
-   if (iterative_mode)
-   {
-      MFEM_WARNING("Multigrid solver does not use iterative_mode and ignores "
-                   "the initial guess!");
-   }
 
    // Add capacity as necessary
    nrhs = X_.Size();
@@ -134,96 +149,159 @@ void MultigridBase::ArrayMult(const Array<const Vector*>& X_,
    for (int j = 0; j < nrhs; ++j)
    {
       MFEM_ASSERT(X_[j] && Y_[j], "Missing Vector in MultigridBase::Mult!");
-      *X(M - 1, j) = *X_[j];
-      *Y(M - 1, j) = 0.0;
-   }
-   Cycle(M - 1);
-   for (int j = 0; j < nrhs; ++j)
-   {
-      *Y_[j] = *Y(M - 1, j);
+      X(M - 1, j) = const_cast<Vector*>(X_[j]);
+      Y(M - 1, j) = Y_[j];
    }
+   const bool zero = !iterative_mode;
+   Cycle(M - 1, zero);
 }
 
 void MultigridBase::SmoothingStep(int level, bool zero, bool transpose) const
 {
    // y = y + S (x - A y) or y = y + S^T (x - A y)
+
+   // Note: 'zero' == true means that Y(level,*) are not initialized and we
+   // should assume that the input they typically provide to this call is zeros.
+
+   // We can't use the smoothers' iterative mode since we don't know if they
+   // actually support it, so we always turn the iterative mode off to properly
+   // use smoothers that do support it.
+   smoothers[level]->iterative_mode = false;
+
    if (zero)
    {
-      Array<Vector *> X_(X[level], nrhs), Y_(Y[level], nrhs);
-      GetSmootherAtLevel(level)->ArrayMult(X_, Y_);
+      MFEM_ASSERT(!transpose, "internal error!");
+      const Array<const Vector *> cX_((const Vector **)(X[level]), nrhs);
+      Array<Vector *> Y_(Y[level], nrhs);
+
+      GetSmootherAtLevel(level)->ArrayMult(cX_, Y_);
    }
    else
    {
-      Array<Vector *> Y_(Y[level], nrhs), R_(R[level], nrhs),
-            Z_(Z[level], nrhs);
+      const Array<const Vector *> cY_((const Vector **)(Y[level]), nrhs),
+            cR_((const Vector **)(R[level]), nrhs);
+      Array<Vector *> Y_(Y[level], nrhs), R_(R[level], nrhs);
+
+      GetOperatorAtLevel(level)->ArrayMult(cY_, R_);
       for (int j = 0; j < nrhs; ++j)
       {
-         *R_[j] = *X(level, j);
+         // *R_[j] = *X(level, j) - *R_[j]
+         subtract(*X(level, j), *R_[j], *R_[j]);
       }
-      GetOperatorAtLevel(level)->ArrayAddMult(Y_, R_, -1.0);
       if (transpose)
       {
-         GetSmootherAtLevel(level)->ArrayMultTranspose(R_, Z_);
+         GetSmootherAtLevel(level)->ArrayAddMultTranspose(cR_, Y_);
       }
       else
       {
-         GetSmootherAtLevel(level)->ArrayMult(R_, Z_);
-      }
-      for (int j = 0; j < nrhs; ++j)
-      {
-         *Y_[j] += *Z_[j];
+         GetSmootherAtLevel(level)->ArrayAddMult(cR_, Y_);
       }
    }
 }
 
-void MultigridBase::Cycle(int level) const
+void MultigridBase::CoarseSolve(bool zero) const
 {
-   // Coarse solve
-   if (level == 0)
+   // See the comment about iterative mode in SmoothingStep()
+   coarse_solver->iterative_mode = false;
+
+   if (zero)
    {
-      SmoothingStep(0, true, false);
+      const Array<const Vector *> cX_((const Vector **)(X[0]), nrhs);
+      Array<Vector *> Y_(Y[0], nrhs);
+
+      coarse_solver->ArrayMult(cX_, Y_);
+   }
+   else
+   {
+      const Array<const Vector *> cY_((const Vector **)(Y[0]), nrhs),
+            cR_((const Vector **)(R[0]), nrhs);
+      Array<Vector *> Y_(Y[0], nrhs), R_(R[0], nrhs);
+
+      GetOperatorAtLevel(0)->ArrayMult(cY_, R_);
+      for (int j = 0; j < nrhs; ++j)
+      {
+         // *R_[j] = *X(0, j) - *R_[j]
+         subtract(*X(0, j), *R_[j], *R_[j]);
+      }
+      coarse_solver->ArrayAddMult(cR_, Y_);
+   }
+}
+
+void MultigridBase::Cycle(int level, bool zero) const
+{
+   // Note: 'zero' == true means that Y(level,*) are not initialized and we
+   // should assume that the input they typically provide to this call is zeros.
+
+   // Coarse solve
+   if (level == 0 && !coarse_solver)
+   {
+      SmoothingStep(0, zero, false);
       return;
    }
 
    // Pre-smooth
    for (int i = 0; i < preSmoothingSteps; ++i)
    {
-      SmoothingStep(level, (cycleType == CycleType::VCYCLE && i == 0), false);
+      SmoothingStep(level, zero && (i == 0), false);
+   }
+
+   // Coarse solve with 'coarse_solver'
+   if (level == 0)
+   {
+      CoarseSolve(preSmoothingSteps == 0 && zero);
+      goto mg_post_smooth;
    }
 
    // Compute residual and restrict
+   if (preSmoothingSteps == 0 && zero)
    {
-      Array<Vector *> Y_(Y[level], nrhs), R_(R[level], nrhs),
-            X_(X[level - 1], nrhs);
+      const Array<const Vector *> cX_l((const Vector **)(X[level]), nrhs);
+      Array<Vector *> X_lm1(X[level - 1], nrhs);
+
+      GetProlongationAtLevel(level - 1)->ArrayMultTranspose(cX_l, X_lm1);
+   }
+   else
+   {
+      const Array<const Vector *> cY_((const Vector **)(Y[level]), nrhs),
+            cR_((const Vector **)(R[level]), nrhs);
+      Array<Vector *> R_(R[level], nrhs), X_(X[level - 1], nrhs);
+
+      GetOperatorAtLevel(level)->ArrayMult(cY_, R_);
       for (int j = 0; j < nrhs; ++j)
       {
-         *R_[j] = *X(level, j);
-      }
-      GetOperatorAtLevel(level)->ArrayAddMult(Y_, R_, -1.0);
-      GetProlongationAtLevel(level - 1)->ArrayMultTranspose(R_, X_);
-      for (int j = 0; j < nrhs; ++j)
-      {
-         *Y(level - 1, j) = 0.0;
+         // *R_[j] = *X(level, j) - *R_[j]
+         subtract(*X(level, j), *R_[j], *R_[j]);
       }
+      GetProlongationAtLevel(level - 1)->ArrayMultTranspose(cR_, X_);
    }
 
    // Corrections
-   Cycle(level - 1);
+   Cycle(level - 1, true);
    if (cycleType == CycleType::WCYCLE)
    {
-      Cycle(level - 1);
+      // If the coarse solve at level 0 is "exact" solve, then we don't want to
+      // repeat it.
+      // To support multiple level 0 coarse-grid corrections, one can wrap that
+      // smoother in an SLI solver and use that instead.
+      if (level > 1) { Cycle(level - 1, false); }
    }
 
    // Prolongate and add
    {
-      Array<Vector *> Y_(Y[level - 1], nrhs), Z_(Z[level], nrhs);
-      GetProlongationAtLevel(level - 1)->ArrayMult(Y_, Z_);
-      for (int j = 0; j < nrhs; ++j)
+      const Array<const Vector *> cY_lm1((const Vector **)(Y[level - 1]), nrhs);
+      Array<Vector *> Y_l(Y[level], nrhs);
+
+      if (preSmoothingSteps == 0 && zero)
       {
-         *Y(level, j) += *Z_[j];
+         GetProlongationAtLevel(level - 1)->ArrayMult(cY_lm1, Y_l);
+      }
+      else
+      {
+         GetProlongationAtLevel(level - 1)->ArrayAddMult(cY_lm1, Y_l);
       }
    }
 
+mg_post_smooth:
    // Post-smooth
    for (int i = 0; i < postSmoothingSteps; ++i)
    {

fem/multigrid.hpp

@@ -36,12 +36,14 @@ protected:
    Array<Solver*> smoothers;
    Array<bool> ownedOperators;
    Array<bool> ownedSmoothers;
+   Solver *coarse_solver; /// can be NULL, see AddCoarseSolver()
+   bool own_coarse_solver;
 
    CycleType cycleType;
    int preSmoothingSteps;
    int postSmoothingSteps;
 
-   mutable Array2D<Vector*> X, Y, R, Z;
+   mutable Array2D<Vector*> X, Y, R;
    mutable int nrhs;
 
 public:
@@ -65,6 +67,16 @@ public:
    void AddLevel(Operator* op, Solver* smoother, bool ownOperator,
                  bool ownSmoother);
 
+   /// Adds a coarse solver for level 0 to work in tandem with the smoother
+   /** If this coarse solver is not given, the smoother at level 0 is used as
+       the coarse solver. When this coarse solver is given, the smoother at
+       level 0 is used similar to the smoothers at other levels. Thus, the
+       action at level 0 consists of:
+       - pre-smoothing steps with smoother 0,
+       - solve step with @a c_solver,
+       - post-smoothing steps with smoother 0. */
+   void AddCoarseSolver(Solver *c_solver, bool own_c_solver);
+
    /// Returns the number of levels
    int NumLevels() const { return operators.Size(); }
 
@@ -118,11 +130,14 @@ public:
 
 private:
    /// Application of a multigrid cycle at particular level
-   void Cycle(int level) const;
+   void Cycle(int level, bool zero) const;
 
    /// Application of a pre-/post-smoothing step at particular level
    void SmoothingStep(int level, bool zero, bool transpose) const;
 
+   /// Perform a coarse solve with 'coarse_solve' (must be non-NULL)
+   void CoarseSolve(bool zero) const;
+
    /// Allocate or destroy temporary storage
    void InitVectors() const;
    void EraseVectors() const;
@@ -202,6 +217,9 @@ public:
 
    /// Recover the solution of a linear system formed with FormFineLinearSystem()
    void RecoverFineFEMSolution(const Vector& X, const Vector& b, Vector& x);
+
+   const Array<int> &GetFineEssentialTrueDofs() const
+   { return *essentialTrueDofs.Last(); }
 };
 
 } // namespace mfem

fem/pbilinearform.cpp

@@ -124,6 +124,8 @@ void ParBilinearForm::pAllocMat()
 void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
                                   bool steal_loc_A)
 {
+   MFEM_PERF_FUNCTION;
+
    ParFiniteElementSpace &pfespace = *ParFESpace();
 
    // Create a block diagonal parallel matrix

fem/pfespace.cpp

@@ -5309,6 +5309,18 @@ DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
       if (recv_size > 0) { req_counter++; }
    }
    requests = new MPI_Request[req_counter];
+   if (mpi_gpu_aware)
+   {
+#if defined(MFEM_USE_HIP)
+      MFEM_GPU_CHECK(
+         hipEventCreateWithFlags(&gpu_event, hipEventDisableTiming));
+#elif defined(MFEM_USE_CUDA)
+      MFEM_GPU_CHECK(
+         cudaEventCreateWithFlags(&gpu_event, cudaEventDisableTiming));
+#else
+      MFEM_ABORT("not implemented");
+#endif
+   }
 }
 
 DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
@@ -5322,16 +5334,27 @@ DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
 }
 
 static void ExtractSubVector(const Array<int> &indices,
-                             const Vector &vin, Vector &vout)
+                             const Vector &vin, Vector &vout,
+                             real_t a, real_t b)
 {
    MFEM_ASSERT(indices.Size() == vout.Size(), "incompatible sizes!");
-   auto y = vout.Write();
+   auto y = (a == 0) ? vout.Write() : vout.ReadWrite();
    const auto x = vin.Read();
    const auto I = indices.Read();
-   mfem::forall(indices.Size(), [=] MFEM_HOST_DEVICE (int i)
+   if (a == 0)
    {
-      y[i] = x[I[i]];
-   }); // indices can be repeated
+      mfem::forall(indices.Size(), [=] MFEM_HOST_DEVICE (int i)
+      {
+         y[i] = b*x[I[i]];
+      }); // indices can be repeated
+   }
+   else
+   {
+      mfem::forall(indices.Size(), [=] MFEM_HOST_DEVICE (int i)
+      {
+         y[i] = a*y[i] + b*x[I[i]];
+      }); // indices can be repeated
+   }
 }
 
 void DeviceConformingProlongationOperator::BcastBeginCopy(
@@ -5339,10 +5362,7 @@ void DeviceConformingProlongationOperator::BcastBeginCopy(
 {
    // shr_buf[i] = src[shr_ltdof[i]]
    if (shr_ltdof.Size() == 0) { return; }
-   ExtractSubVector(shr_ltdof, x, shr_buf);
-   // If the above kernel is executed asynchronously, we should wait for it to
-   // complete
-   if (mpi_gpu_aware) { MFEM_STREAM_SYNC; }
+   ExtractSubVector(shr_ltdof, x, shr_buf, 0, 1);
 }
 
 static void SetSubVector(const Array<int> &indices,
@@ -5379,7 +5399,7 @@ void DeviceConformingProlongationOperator::Mult(const Vector &x,
                                                 Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
-   int req_counter = 0;
+   int req_counter = 0, num_recv_req = 0;
    // Make sure 'y' is marked as valid on device and for use on device.
    // This ensures that there is no unnecessary host to device copy when the
    // input 'y' is valid on host (in 'y.SetSubVector(ext_ldof, 0.0)' when local
@@ -5389,42 +5409,89 @@ void DeviceConformingProlongationOperator::Mult(const Vector &x,
    {
       // done on device since we've marked ext_ldof for use on device:
       y.SetSubVector(ext_ldof, 0.0);
+      BcastLocalCopy(x, y);
+      return;
    }
-   else
+
+   BcastBeginCopy(x); // copy to 'shr_buf'
+   if (mpi_gpu_aware && shr_ltdof.Size() != 0)
    {
-      BcastBeginCopy(x); // copy to 'shr_buf'
+      /* record a stream event to wait for later */
+#if defined(MFEM_USE_HIP)
+      MFEM_GPU_CHECK(hipEventRecord(gpu_event, 0));
+#elif defined(MFEM_USE_CUDA)
+      MFEM_GPU_CHECK(cudaEventRecord(gpu_event, 0));
+#endif
+   }
+   BcastLocalCopy(x, y);
+   // Queue all receive communications
+   if (ext_ldof.Size() != 0) // ext_ldof.Size() == ext_buf.Size()
+   {
+      auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
+      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+      {
+         const int recv_offset = ext_buf_offsets[nbr];
+         const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
+         if (recv_size > 0)
+         {
+            MPI_Irecv(recv_buf + recv_offset, recv_size,
+                      MPITypeMap<real_t>::mpi_type,
+                      gtopo.GetNeighborRank(nbr), 41822,
+                      gtopo.GetComm(), &requests[req_counter++]);
+         }
+      }
+      num_recv_req = req_counter;
+   }
+   // Queue all send communications
+   if (shr_ltdof.Size() != 0) // shr_ltdof.Size() == shr_buf.Size()
+   {
+      // The BcastBeginCopy kernel is executed asynchronously, we should wait
+      // for it to complete:
+      // - when mpi_gpu_aware == false, this is done implicily when we call
+      //   shr_buf.HostRead()
+      // - when mpi_gpu_aware == true, we need to wait for BcastBeginCopy to
+      //   complete by waiting for gpu_event.
+      if (mpi_gpu_aware)
+      {
+         /* wait for the stream event recorded above */
+#if defined(MFEM_USE_HIP)
+         MFEM_GPU_CHECK(hipEventSynchronize(gpu_event));
+#elif defined(MFEM_USE_CUDA)
+         MFEM_GPU_CHECK(cudaEventSynchronize(gpu_event));
+#endif
+      }
+      auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
       for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
       {
          const int send_offset = shr_buf_offsets[nbr];
          const int send_size = shr_buf_offsets[nbr+1] - send_offset;
          if (send_size > 0)
          {
-            auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
-            MPI_Isend(send_buf + send_offset, send_size, MPITypeMap<real_t>::mpi_type,
-                      gtopo.GetNeighborRank(nbr), 41822,
-                      gtopo.GetComm(), &requests[req_counter++]);
-         }
-         const int recv_offset = ext_buf_offsets[nbr];
-         const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
-         if (recv_size > 0)
-         {
-            auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
-            MPI_Irecv(recv_buf + recv_offset, recv_size, MPITypeMap<real_t>::mpi_type,
+            MPI_Isend(send_buf + send_offset, send_size,
+                      MPITypeMap<real_t>::mpi_type,
                       gtopo.GetNeighborRank(nbr), 41822,
                       gtopo.GetComm(), &requests[req_counter++]);
          }
       }
    }
-   BcastLocalCopy(x, y);
-   if (!local)
-   {
-      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-      BcastEndCopy(y); // copy from 'ext_buf'
-   }
+   // Wait for all receive requests
+   MPI_Waitall(num_recv_req, requests, MPI_STATUSES_IGNORE);
+   BcastEndCopy(y); // copy from 'ext_buf'
+   // Wait for all send requests
+   MPI_Waitall(req_counter - num_recv_req, requests + num_recv_req,
+               MPI_STATUSES_IGNORE);
 }
 
 DeviceConformingProlongationOperator::~DeviceConformingProlongationOperator()
 {
+   if (mpi_gpu_aware)
+   {
+#if defined(MFEM_USE_HIP)
+      MFEM_GPU_CHECK(hipEventDestroy(gpu_event));
+#elif defined(MFEM_USE_CUDA)
+      MFEM_GPU_CHECK(cudaEventDestroy(gpu_event));
+#endif
+   }
    delete [] requests;
    ext_buf_offsets.Delete();
    shr_buf_offsets.Delete();
@@ -5435,25 +5502,23 @@ void DeviceConformingProlongationOperator::ReduceBeginCopy(
 {
    // ext_buf[i] = src[ext_ldof[i]]
    if (ext_ldof.Size() == 0) { return; }
-   ExtractSubVector(ext_ldof, x, ext_buf);
-   // If the above kernel is executed asynchronously, we should wait for it to
-   // complete
-   if (mpi_gpu_aware) { MFEM_STREAM_SYNC; }
+   ExtractSubVector(ext_ldof, x, ext_buf, 0, 1);
 }
 
 void DeviceConformingProlongationOperator::ReduceLocalCopy(
-   const Vector &x, Vector &y) const
+   const Vector &x, Vector &y, real_t a, real_t b) const
 {
    // dst[i] = src[ltdof_ldof[i]]
    if (ltdof_ldof.Size() == 0) { return; }
-   ExtractSubVector(ltdof_ldof, x, y);
+   ExtractSubVector(ltdof_ldof, x, y, a, b);
 }
 
 static void AddSubVector(const Array<int> &unique_dst_indices,
                          const Array<int> &unique_to_src_offsets,
                          const Array<int> &unique_to_src_indices,
                          const Vector &src,
-                         Vector &dst)
+                         Vector &dst,
+                         real_t b)
 {
    auto y = dst.ReadWrite();
    const auto x = src.Read();
@@ -5463,56 +5528,94 @@ static void AddSubVector(const Array<int> &unique_dst_indices,
    mfem::forall(unique_dst_indices.Size(), [=] MFEM_HOST_DEVICE (int i)
    {
       const int dst_idx = DST_I[i];
-      real_t sum = y[dst_idx];
+      real_t sum = 0;
       const int end = SRC_O[i+1];
       for (int j = SRC_O[i]; j != end; ++j) { sum += x[SRC_I[j]]; }
-      y[dst_idx] = sum;
+      y[dst_idx] += b*sum;
    });
 }
 
-void DeviceConformingProlongationOperator::ReduceEndAssemble(Vector &y) const
+void DeviceConformingProlongationOperator::ReduceEndAssemble(
+   Vector &y, real_t b) const
 {
    // dst[shr_ltdof[i]] += shr_buf[i]
    if (unq_ltdof.Size() == 0) { return; }
-   AddSubVector(unq_ltdof, unq_shr_i, unq_shr_j, shr_buf, y);
+   AddSubVector(unq_ltdof, unq_shr_i, unq_shr_j, shr_buf, y, b);
 }
 
-void DeviceConformingProlongationOperator::MultTranspose(const Vector &x,
-                                                         Vector &y) const
+void DeviceConformingProlongationOperator::ApplyTranspose(
+   const Vector &x, Vector &y, real_t a, real_t b) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
-   int req_counter = 0;
-   if (!local)
+   int req_counter = 0, num_recv_req = 0;
+
+   if (local)
    {
-      ReduceBeginCopy(x); // copy to 'ext_buf'
+      ReduceLocalCopy(x, y, a, b);
+      return;
+   }
+
+   ReduceBeginCopy(x); // copy to 'ext_buf'
+   if (mpi_gpu_aware && ext_ldof.Size() != 0)
+   {
+      /* record a stream event to wait for later */
+#if defined(MFEM_USE_HIP)
+      MFEM_GPU_CHECK(hipEventRecord(gpu_event, 0));
+#elif defined(MFEM_USE_CUDA)
+      MFEM_GPU_CHECK(cudaEventRecord(gpu_event, 0));
+#endif
+   }
+   ReduceLocalCopy(x, y, a, b);
+   // Queue all receive communications
+   if (unq_ltdof.Size() != 0)
+   {
+      auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
+      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+      {
+         const int recv_offset = shr_buf_offsets[nbr];
+         const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
+         if (recv_size > 0)
+         {
+            MPI_Irecv(recv_buf + recv_offset, recv_size,
+                      MPITypeMap<real_t>::mpi_type,
+                      gtopo.GetNeighborRank(nbr), 41823,
+                      gtopo.GetComm(), &requests[req_counter++]);
+         }
+      }
+      num_recv_req = req_counter;
+   }
+   // Queue all send communications
+   if (ext_ldof.Size() != 0)
+   {
+      if (mpi_gpu_aware)
+      {
+         /* wait for the stream event recorded above */
+#if defined(MFEM_USE_HIP)
+         MFEM_GPU_CHECK(hipEventSynchronize(gpu_event));
+#elif defined(MFEM_USE_CUDA)
+         MFEM_GPU_CHECK(cudaEventSynchronize(gpu_event));
+#endif
+      }
+      auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
       for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
       {
          const int send_offset = ext_buf_offsets[nbr];
          const int send_size = ext_buf_offsets[nbr+1] - send_offset;
          if (send_size > 0)
          {
-            auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
-            MPI_Isend(send_buf + send_offset, send_size, MPITypeMap<real_t>::mpi_type,
-                      gtopo.GetNeighborRank(nbr), 41823,
-                      gtopo.GetComm(), &requests[req_counter++]);
-         }
-         const int recv_offset = shr_buf_offsets[nbr];
-         const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
-         if (recv_size > 0)
-         {
-            auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
-            MPI_Irecv(recv_buf + recv_offset, recv_size, MPITypeMap<real_t>::mpi_type,
+            MPI_Isend(send_buf + send_offset, send_size,
+                      MPITypeMap<real_t>::mpi_type,
                       gtopo.GetNeighborRank(nbr), 41823,
                       gtopo.GetComm(), &requests[req_counter++]);
          }
       }
    }
-   ReduceLocalCopy(x, y);
-   if (!local)
-   {
-      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-      ReduceEndAssemble(y); // assemble from 'shr_buf'
-   }
+   // Wait for all receive requests
+   MPI_Waitall(num_recv_req, requests, MPI_STATUSES_IGNORE);
+   ReduceEndAssemble(y, b); // assemble from 'shr_buf'
+   // Wait for all send requests
+   MPI_Waitall(req_counter - num_recv_req, requests + num_recv_req,
+               MPI_STATUSES_IGNORE);
 }
 
 } // namespace mfem

fem/pfespace.hpp

@@ -586,6 +586,7 @@ public:
    { MultTranspose(x,y); }
 };
 
+
 /// Auxiliary device class used by ParFiniteElementSpace.
 class DeviceConformingProlongationOperator: public
    ConformingProlongationOperator
@@ -598,6 +599,11 @@ protected:
    Array<int> ltdof_ldof, unq_ltdof;
    Array<int> unq_shr_i, unq_shr_j;
    MPI_Request *requests;
+#if defined(MFEM_USE_HIP)
+   hipEvent_t gpu_event;
+#elif defined(MFEM_USE_CUDA)
+   cudaEvent_t gpu_event;
+#endif
 
    // Kernel: copy ltdofs from 'src' to 'shr_buf' - prepare for send.
    //         shr_buf[i] = src[shr_ltdof[i]]
@@ -617,11 +623,12 @@ protected:
 
    // Kernel: copy owned ldofs from 'src' to ltdofs in 'dst'.
    //         dst[i] = src[ltdof_ldof[i]]
-   void ReduceLocalCopy(const Vector &src, Vector &dst) const;
+   void ReduceLocalCopy(const Vector &src, Vector &dst,
+                        real_t a, real_t b) const;
 
    // Kernel: assemble dofs from 'shr_buf' into to 'dst' - after recv.
    //         dst[shr_ltdof[i]] += shr_buf[i]
-   void ReduceEndAssemble(Vector &dst) const;
+   void ReduceEndAssemble(Vector &dst, real_t b) const;
 
 public:
    DeviceConformingProlongationOperator(
@@ -632,12 +639,20 @@ public:
 
    virtual ~DeviceConformingProlongationOperator();
 
+   void ApplyTranspose(const Vector &x, Vector &y,
+                       real_t a, real_t b) const;
+
    void Mult(const Vector &x, Vector &y) const override;
 
    void AbsMult(const Vector &x, Vector &y) const override
    { Mult(x,y); }
 
-   void MultTranspose(const Vector &x, Vector &y) const override;
+   void MultTranspose(const Vector &x, Vector &y) const override
+   { ApplyTranspose(x, y, 0, 1); }
+
+   void AddMultTranspose(const Vector &x, Vector &y,
+                         real_t a = 1) const override
+   { ApplyTranspose(x, y, 1, a); }
 
    void AbsMultTranspose(const Vector &x, Vector &y) const override
    { MultTranspose(x,y); }

fem/quadinterpolator.cpp

@@ -503,6 +503,8 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
                                   Vector &q_der,
                                   Vector &q_det) const
 {
+   MFEM_PERF_FUNCTION;
+
    using namespace internal::quadrature_interpolator;
 
    const int ne = fespace->GetNE();

fem/transfer.cpp

@@ -13,6 +13,7 @@
 #include "bilinearform.hpp"
 #include "pbilinearform.hpp"
 #include "../general/forall.hpp"
+#include "kernels.hpp"
 
 namespace mfem
 {
@@ -2322,6 +2323,74 @@ void Prolongation2D(const int NE, const int D1D, const int Q1D,
    });
 }
 
+template <int DLO, int DHI>
+static void SmemProlongation3D(const int NE,
+                               const Vector& localL, Vector& localH,
+                               const Array<real_t> &b, const Vector& mask)
+{
+   auto u_lo = Reshape(localL.Read(), DLO, DLO, DLO, NE);
+   auto u_hi = Reshape(localH.Write(), DHI, DHI, DHI, NE);
+   auto d_b = Reshape(b.Read(), DHI, DLO);
+   auto m_ = Reshape(mask.Read(), DHI, DHI, DHI, NE);
+
+   mfem::forall_2D(NE, DHI, DHI, [=] MFEM_HOST_DEVICE (int e)
+   {
+      // Load B into shared memory
+      MFEM_SHARED real_t s_B[DHI*DLO];
+      kernels::internal::LoadBt<DLO,DHI>(DLO,DHI,d_b,s_B);
+      const DeviceMatrix B(s_B, DHI, DLO);
+
+      MFEM_SHARED real_t s_u[DHI*DHI*DLO];
+      const DeviceCube u(s_u, DHI, DHI, DLO);
+      real_t v[DHI];
+
+      MFEM_FOREACH_THREAD(lx,x,DLO)
+      {
+         MFEM_FOREACH_THREAD(ly,y,DLO)
+         {
+            for (int hz = 0; hz < DHI; ++hz) { v[hz] = 0.0; }
+            for (int lz = 0; lz < DLO; ++lz)
+            {
+               const real_t XYZ = u_lo(lx,ly,lz,e);
+               for (int hz = 0; hz < DHI; ++hz) { v[hz] += XYZ * B(hz,lz); }
+            }
+            for (int hz = 0; hz < DHI; ++hz) { u(hz,ly,lx) = v[hz]; }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(hz,y,DHI)
+      {
+         MFEM_FOREACH_THREAD(lx,x,DLO)
+         {
+            for (int hy = 0; hy < DHI; ++hy) { v[hy] = 0.0; }
+            for (int ly = 0; ly < DLO; ++ly)
+            {
+               const real_t zYX = u(hz,ly,lx);
+               for (int hy = 0; hy < DHI; ++hy) { v[hy] += zYX * B(hy,ly); }
+            }
+            for (int hy = 0; hy < DHI; ++hy) { u(hz,hy,lx) = v[hy]; }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(hz,y,DHI)
+      {
+         MFEM_FOREACH_THREAD(hy,x,DHI)
+         {
+            for (int hx = 0; hx < DHI; ++hx) { v[hx] = 0.0; }
+            for (int lx = 0; lx < DLO; ++lx)
+            {
+               const real_t zyX = u(hz,hy,lx);
+               for (int hx = 0; hx < DHI; ++hx) { v[hx] += zyX * B(hx,lx); }
+            }
+            for (int hx = 0; hx < DHI; ++hx)
+            {
+               u_hi(hx,hy,hz,e) = m_(hx,hy,hz,e)*v[hx];
+            }
+         }
+      }
+   });
+}
+
 void Prolongation3D(const int NE, const int D1D, const int Q1D,
                     const Vector& localL, Vector& localH,
                     const Array<real_t>& B, const Vector& mask)
@@ -2403,9 +2472,9 @@ void Prolongation3D(const int NE, const int D1D, const int Q1D,
    });
 }
 
-void Restriction2D(const int NE, const int D1D, const int Q1D,
-                   const Vector& localH, Vector& localL,
-                   const Array<real_t>& Bt, const Vector& mask)
+void ProlongationTranspose2D(const int NE, const int D1D, const int Q1D,
+                             const Vector& localH, Vector& localL,
+                             const Array<real_t>& Bt, const Vector& mask)
 {
    auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE);
    auto y_ = Reshape(localL.Write(), D1D, D1D, NE);
@@ -2448,9 +2517,78 @@ void Restriction2D(const int NE, const int D1D, const int Q1D,
       }
    });
 }
-void Restriction3D(const int NE, const int D1D, const int Q1D,
-                   const Vector& localH, Vector& localL,
-                   const Array<real_t>& Bt, const Vector& mask)
+
+template <int DLO, int DHI>
+static void SmemProlongationTranspose3D(
+   const int NE, const Vector& localH, Vector& localL,
+   const Array<real_t>& bt, const Vector& mask)
+{
+   auto u_h = Reshape(localH.Read(), DHI, DHI, DHI, NE);
+   auto u_l = Reshape(localL.Write(), DLO, DLO, DLO, NE);
+   auto d_bt = Reshape(bt.Read(), DLO, DHI);
+   auto m_ = Reshape(mask.Read(), DHI, DHI, DHI, NE);
+
+   mfem::forall_2D(NE, DHI, DHI, [=] MFEM_HOST_DEVICE (int e)
+   {
+      // Load Bt into shared memory
+      MFEM_SHARED real_t s_Bt[DHI*DLO];
+      kernels::internal::LoadBt<DHI,DLO>(DHI,DLO,d_bt,s_Bt);
+      const DeviceMatrix Bt(s_Bt, DLO, DHI);
+
+      MFEM_SHARED real_t s_u[DLO*DHI*DHI];
+      const DeviceCube u(s_u, DLO, DHI, DHI);
+      real_t v[DLO];
+
+      MFEM_FOREACH_THREAD(hx,x,DHI)
+      {
+         MFEM_FOREACH_THREAD(hy,y,DHI)
+         {
+            for (int lz = 0; lz < DLO; ++lz) { v[lz] = 0.0; }
+            for (int hz = 0; hz < DHI; ++hz)
+            {
+               const real_t XYZ = m_(hx,hy,hz,e)*u_h(hx,hy,hz,e);
+               for (int lz = 0; lz < DLO; ++lz) { v[lz] += XYZ * Bt(lz,hz); }
+            }
+            for (int lz = 0; lz < DLO; ++lz) { u(lz,hy,hx) = v[lz]; }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(lz,y,DLO)
+      {
+         MFEM_FOREACH_THREAD(hx,x,DHI)
+         {
+            for (int ly = 0; ly < DLO; ++ly) { v[ly] = 0.0; }
+            for (int hy = 0; hy < DHI; ++hy)
+            {
+               const real_t zYX = u(lz,hy,hx);
+               for (int ly = 0; ly < DLO; ++ly) { v[ly] += zYX * Bt(ly,hy); }
+            }
+            for (int ly = 0; ly < DLO; ++ly) { u(lz,ly,hx) = v[ly]; }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(lz,y,DLO)
+      {
+         MFEM_FOREACH_THREAD(ly,x,DLO)
+         {
+            for (int lx = 0; lx < DLO; ++lx) { v[lx] = 0.0; }
+            for (int hx = 0; hx < DHI; ++hx)
+            {
+               const real_t zyX = u(lz,ly,hx);
+               for (int lx = 0; lx < DLO; ++lx) { v[lx] += zyX * Bt(lx,hx); }
+            }
+            for (int lx = 0; lx < DLO; ++lx)
+            {
+               u_l(lx,ly,lz,e) = v[lx];
+            }
+         }
+      }
+   });
+}
+
+void ProlongationTranspose3D(const int NE, const int D1D, const int Q1D,
+                             const Vector& localH, Vector& localL,
+                             const Array<real_t>& Bt, const Vector& mask)
 {
    auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE);
    auto y_ = Reshape(localL.Write(), D1D, D1D, D1D, NE);
@@ -2518,11 +2656,14 @@ void Restriction3D(const int NE, const int D1D, const int Q1D,
       }
    });
 }
+
 } // namespace TransferKernels
 
 void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
                                                     Vector& y) const
 {
+   using namespace TransferKernels;
+
    if (lFESpace.GetMesh()->GetNE() == 0)
    {
       return;
@@ -2531,11 +2672,25 @@ void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
    elem_restrict_lex_l->Mult(x, localL);
    if (dim == 2)
    {
-      TransferKernels::Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
+      Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
    }
    else if (dim == 3)
    {
-      TransferKernels::Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask);
+      switch ((D1D << 4 ) | Q1D)
+      {
+         case 0x23:
+            SmemProlongation3D<2,3>(NE, localL, localH, B, mask); break;
+         case 0x24:
+            SmemProlongation3D<2,4>(NE, localL, localH, B, mask); break;
+         case 0x35:
+            SmemProlongation3D<3,5>(NE, localL, localH, B, mask); break;
+         case 0x46:
+            SmemProlongation3D<4,6>(NE, localL, localH, B, mask); break;
+         case 0x47:
+            SmemProlongation3D<4,7>(NE, localL, localH, B, mask); break;
+         default:
+            Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask); break;
+      }
    }
    else
    {
@@ -2549,6 +2704,8 @@ void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
 void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
                                                              Vector& y) const
 {
+   using namespace TransferKernels;
+
    if (lFESpace.GetMesh()->GetNE() == 0)
    {
       return;
@@ -2557,11 +2714,25 @@ void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
    elem_restrict_lex_h->Mult(x, localH);
    if (dim == 2)
    {
-      TransferKernels::Restriction2D(NE, D1D, Q1D, localH, localL, Bt, mask);
+      ProlongationTranspose2D(NE, D1D, Q1D, localH, localL, Bt, mask);
    }
    else if (dim == 3)
    {
-      TransferKernels::Restriction3D(NE, D1D, Q1D, localH, localL, Bt, mask);
+      switch ((D1D << 4 ) | Q1D)
+      {
+         case 0x23:
+            SmemProlongationTranspose3D<2,3>(NE, localH, localL, Bt, mask); break;
+         case 0x24:
+            SmemProlongationTranspose3D<2,4>(NE, localH, localL, Bt, mask); break;
+         case 0x35:
+            SmemProlongationTranspose3D<3,5>(NE, localH, localL, Bt, mask); break;
+         case 0x46:
+            SmemProlongationTranspose3D<4,6>(NE, localH, localL, Bt, mask); break;
+         case 0x47:
+            SmemProlongationTranspose3D<4,7>(NE, localH, localL, Bt, mask); break;
+         default:
+            ProlongationTranspose3D(NE, D1D, Q1D, localH, localL, Bt, mask); break;
+      }
    }
    else
    {

fem/transfer.hpp

@@ -621,9 +621,6 @@ public:
       const FiniteElementSpace& lFESpace_,
       const FiniteElementSpace& hFESpace_);
 
-   /// Destructor
-   virtual ~TensorProductPRefinementTransferOperator() { }
-
    /// @brief Interpolation or prolongation of a vector \p x corresponding to
    /// the coarse space to the vector \p y corresponding to the fine space.
    void Mult(const Vector& x, Vector& y) const override;

general/annotation.hpp

@@ -14,24 +14,111 @@
 
 #include "../config/config.hpp"
 
-#ifdef MFEM_USE_CALIPER
+#define MFEM_CONCAT_(X,Y)  X##Y
+#define MFEM_CONCAT(X,Y)  MFEM_CONCAT_(X,Y)
 
+#ifdef MFEM_USE_CALIPER
+#include "device.hpp"
+#include "backends.hpp"
+#ifdef MFEM_USE_MPI
+#include "communication.hpp"
+#endif
+#include <optional>
 #include <caliper/cali.h>
 #include <caliper/cali-manager.h>
-#define MFEM_PERF_FUNCTION CALI_CXX_MARK_FUNCTION
-#define MFEM_PERF_BEGIN(s) CALI_MARK_BEGIN(s)
-#define MFEM_PERF_END(s) CALI_MARK_END(s)
+#endif
+
+namespace mfem
+{
+
+namespace internal
+{
+
+extern int annotation_sync_stream; // defined in globals.cpp
+extern int annotation_sync_mpi; // defined in globals.cpp
+extern int annotation_enabled; // defined in globals.cpp
+
+#ifdef MFEM_USE_CALIPER
+
+inline void AnnotationSync()
+{
+   if (!annotation_enabled) { return; }
+   if (annotation_sync_stream && Device::Allows(Backend::DEVICE_MASK))
+   {
+      MFEM_STREAM_SYNC;
+   }
+#ifdef MFEM_USE_MPI
+   if (annotation_sync_mpi && Mpi::IsInitialized() && !Mpi::IsFinalized())
+   {
+      MPI_Barrier(GetGlobalMPI_Comm());
+   }
+#endif
+}
+
+struct FunctionAnnotation
+{
+   std::optional<cali::Function> cali_func;
+
+   FunctionAnnotation(const char *fname)
+   {
+      AnnotationSync();
+      if (annotation_enabled) { cali_func.emplace(fname); }
+   }
+
+   ~FunctionAnnotation() { AnnotationSync(); }
+};
+
+struct ScopeAnnotation
+{
+   std::optional<cali::ScopeAnnotation> cali_scope;
+
+   ScopeAnnotation(const char *name)
+   {
+      AnnotationSync();
+      if (annotation_enabled) { cali_scope.emplace(name); }
+   }
+
+   ~ScopeAnnotation() { AnnotationSync(); }
+};
+
+#endif // #ifdef MFEM_USE_CALIPER
+
+} // namespace internal
+
+} // namespace mfem
+
+
+#ifdef MFEM_USE_CALIPER
+
+#define MFEM_PERF_FUNCTION \
+   mfem::internal::FunctionAnnotation mfem_func_annotation_(_MFEM_FUNC_NAME)
+#define MFEM_PERF_BEGIN(s) \
+   (mfem::internal::AnnotationSync(), CALI_MARK_BEGIN(s))
+#define MFEM_PERF_END(s) \
+   (mfem::internal::AnnotationSync(), CALI_MARK_END(s))
 #define MFEM_PERF_SCOPE(name) \
-  cali::Annotation::Guard cali_autogenerated_guard_name(cali::Annotation("function").begin(std::string(name).c_str()))
+   mfem::internal::ScopeAnnotation \
+      MFEM_CONCAT(mfem_scope_annotation_,__LINE__)(name)
 
+#define MFEM_PERF_SYNC_STREAM(b) (mfem::internal::annotation_sync_stream = (b))
+#define MFEM_PERF_SYNC_MPI(b) (mfem::internal::annotation_sync_mpi = (b))
+#define MFEM_PERF_SYNC(b) (MFEM_PERF_SYNC_STREAM(b), MFEM_PERF_SYNC_MPI(b))
+#define MFEM_PERF_ENABLE (mfem::internal::annotation_enabled = 1)
+#define MFEM_PERF_DISABLE (mfem::internal::annotation_enabled = 0)
 
-#else
+#else // #ifdef MFEM_USE_CALIPER
 
 #define MFEM_PERF_FUNCTION
 #define MFEM_PERF_BEGIN(s)
 #define MFEM_PERF_END(s)
 #define MFEM_PERF_SCOPE(name)
 
-#endif
+#define MFEM_PERF_SYNC_STREAM(b)
+#define MFEM_PERF_SYNC_MPI(b)
+#define MFEM_PERF_SYNC(b)
+#define MFEM_PERF_ENABLE
+#define MFEM_PERF_DISABLE
 
-#endif
+#endif // #ifdef MFEM_USE_CALIPER
+
+#endif // MFEM_ANNOTATION_HPP

general/communication.hpp

@@ -23,13 +23,9 @@
 #include <mpi.h>
 #include <cstdint>
 
-// can't directly use MPI_CXX_BOOL because Microsoft's MPI implementation
-// doesn't include MPI_CXX_BOOL. Fallback to MPI_C_BOOL if unavailable.
-#ifdef MPI_CXX_BOOL
-#define MFEM_MPI_CXX_BOOL MPI_CXX_BOOL
-#else
-#define MFEM_MPI_CXX_BOOL MPI_C_BOOL
-#endif
+// Some MPI implementations do not have MPI_CXX_BOOL or do not handle it
+// correctly, so we use MPI_UNSIGNED_CHAR as the MPI type for 'bool':
+#define MFEM_MPI_CXX_BOOL MPI_UNSIGNED_CHAR
 
 namespace mfem
 {

general/device.cpp

@@ -151,6 +151,22 @@ Device::Device()
    {
       SetGPUAwareMPI(true);
    }
+
+   if (const char *mfem_perf_sync = GetEnv("MFEM_PERF_SYNC"))
+   {
+      MFEM_PERF_SYNC(std::atoi(mfem_perf_sync));
+      MFEM_CONTRACT_VAR(mfem_perf_sync);
+   }
+   if (const char *mfem_perf_sync_stream = GetEnv("MFEM_PERF_SYNC_STREAM"))
+   {
+      MFEM_PERF_SYNC_STREAM(std::atoi(mfem_perf_sync_stream));
+      MFEM_CONTRACT_VAR(mfem_perf_sync_stream);
+   }
+   if (const char *mfem_perf_sync_mpi = GetEnv("MFEM_PERF_SYNC_MPI"))
+   {
+      MFEM_PERF_SYNC_MPI(std::atoi(mfem_perf_sync_mpi));
+      MFEM_CONTRACT_VAR(mfem_perf_sync_mpi);
+   }
 }
 
 Device::~Device()

general/error.cpp

@@ -193,4 +193,4 @@ void mfem_warning(const char *msg)
    }
 }
 
-}
+} // namespace mfem

general/error.hpp

@@ -208,4 +208,4 @@ __device__ void abort_msg(T & msg)
 #define MFEM_ASSERT_KERNEL(x,...)
 #endif
 
-#endif
+#endif // MFEM_ERROR_HPP

general/globals.cpp

@@ -31,6 +31,10 @@ namespace internal
 {
 bool mfem_out_initialized = false;
 bool mfem_err_initialized = false;
+
+int annotation_sync_stream = 0; // declared in annotation.hpp
+int annotation_sync_mpi = 0; // declared in annotation.hpp
+int annotation_enabled = 1; // declared in annotation.hpp
 }
 
 void OutStream::Init()

general/mem_manager.hpp

@@ -657,7 +657,8 @@ private: // Static methods used by the Memory<T> class
    /// Return the host pointer.
    MFEM_ENZYME_INACTIVE static void *Register_(void *ptr, void *h_ptr,
                                                size_t bytes, MemoryType mt,
-                                               bool own, bool alias, unsigned &flags);
+                                               bool own, bool alias,
+                                               unsigned &flags);
 
    /// Register a pair of external host and device pointers
    static void Register2_(void *h_ptr, void *d_ptr, size_t bytes,
@@ -741,7 +742,7 @@ private:
 
    /// Insert a host address @a h_ptr and size *a bytes in the memory map to be
    /// managed.
-   void Insert(void *h_ptr, size_t bytes, MemoryType h_mt,  MemoryType d_mt);
+   void Insert(void *h_ptr, size_t bytes, MemoryType h_mt, MemoryType d_mt);
 
    /// Insert a device and the host addresses in the memory map
    void InsertDevice(void *d_ptr, void *h_ptr, size_t bytes,
@@ -980,7 +981,7 @@ inline void Memory<T>::Wrap(T *ptr, int size, bool own)
 #ifdef MFEM_DEBUG
    if (own && MemoryManager::Exists())
    {
-      MemoryType h_ptr_mt = MemoryManager::GetHostMemoryType_(h_ptr);
+      MemoryType h_ptr_mt = MemoryManager::GetHostMemoryType_((void*)h_ptr);
       MFEM_VERIFY(h_mt == h_ptr_mt,
                   "h_mt = " << (int)h_mt << ", h_ptr_mt = " << (int)h_ptr_mt);
    }
@@ -988,7 +989,8 @@ inline void Memory<T>::Wrap(T *ptr, int size, bool own)
    if (own && h_mt != MemoryType::HOST)
    {
       const size_t bytes = size*sizeof(T);
-      MemoryManager::Register_(ptr, ptr, bytes, h_mt, own, false, flags);
+      MemoryManager::Register_((void*)ptr, (void*)ptr, bytes, h_mt, own, false,
+                               flags);
    }
 }
 

linalg/batched/gpu_blas.cpp

@@ -20,13 +20,13 @@
 #define MFEM_CU_or_HIP(stub) HIP##stub
 #endif
 
-#define MFEM_CONCAT(x, y, z) MFEM_CONCAT_(x, y, z)
-#define MFEM_CONCAT_(x, y, z) x ## y ## z
+#define MFEM_CONCAT3(x, y, z) MFEM_CONCAT3_(x, y, z)
+#define MFEM_CONCAT3_(x, y, z) x ## y ## z
 
 #ifdef MFEM_USE_SINGLE
-#define MFEM_GPUBLAS_PREFIX(stub) MFEM_CONCAT(MFEM_cu_or_hip(blas), S, stub)
+#define MFEM_GPUBLAS_PREFIX(stub) MFEM_CONCAT3(MFEM_cu_or_hip(blas), S, stub)
 #elif defined(MFEM_USE_DOUBLE)
-#define MFEM_GPUBLAS_PREFIX(stub) MFEM_CONCAT(MFEM_cu_or_hip(blas), D, stub)
+#define MFEM_GPUBLAS_PREFIX(stub) MFEM_CONCAT3(MFEM_cu_or_hip(blas), D, stub)
 #endif
 
 #define MFEM_BLAS_SUCCESS MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS)

linalg/hypre.cpp

@@ -4091,6 +4091,8 @@ void HypreSolver::Setup(const HypreParVector &b, HypreParVector &x) const
 {
    if (setup_called) { return; }
 
+   MFEM_PERF_FUNCTION;
+
    MFEM_VERIFY(A != NULL, "HypreParMatrix A is missing");
 
    HYPRE_Int err_flag = SetupFcn()(*this, *A, b, x);
@@ -4116,6 +4118,8 @@ void HypreSolver::Setup(const Vector &b, Vector &x) const
 
 void HypreSolver::Mult(const HypreParVector &b, HypreParVector &x) const
 {
+   MFEM_PERF_FUNCTION;
+
    HYPRE_Int err_flag;
    if (A == NULL)
    {

linalg/operator.cpp

@@ -50,17 +50,19 @@ void Operator::InitTVectors(const Operator *Po, const Operator *Ri,
 
 void Operator::AddMult(const Vector &x, Vector &y, const real_t a) const
 {
-   mfem::Vector z(y.Size());
-   Mult(x, z);
-   y.Add(a, z);
+   z_am.SetSize(y.Size());
+   z_am.UseDevice(true);
+   Mult(x, z_am);
+   y.Add(a, z_am);
 }
 
 void Operator::AddMultTranspose(const Vector &x, Vector &y,
                                 const real_t a) const
 {
-   mfem::Vector z(y.Size());
-   MultTranspose(x, z);
-   y.Add(a, z);
+   z_am.SetSize(y.Size());
+   z_am.UseDevice(true);
+   MultTranspose(x, z_am);
+   y.Add(a, z_am);
 }
 
 void Operator::ArrayMult(const Array<const Vector *> &X,

linalg/operator.hpp

@@ -23,6 +23,13 @@ class RectangularConstrainedOperator;
 /// Abstract operator
 class Operator
 {
+private:
+   /// Auxiliary Vector used by the methods AddMult() and AddMultTranspose().
+   /** @note This Vector is private to prevent derived classes from accidentaly
+       using it in their implementation of Mult() or MultTranspose() which may
+       lead to hard-to-find bugs. */
+   mutable Vector z_am;
+
 protected:
    int height; ///< Dimension of the output / number of rows in the matrix.
    int width;  ///< Dimension of the input / number of columns in the matrix.

linalg/solvers.cpp

@@ -314,25 +314,29 @@ void OperatorJacobiSmoother::Mult(const Vector &x, Vector &y) const
    MFEM_VERIFY(x.Size() == Width(), "invalid input vector");
    MFEM_VERIFY(y.Size() == Height(), "invalid output vector");
 
+   auto DI = dinv.Read();
+   auto X = x.Read();
    if (iterative_mode)
    {
       MFEM_VERIFY(oper, "iterative_mode == true requires the forward operator");
       oper->Mult(y, residual);  // r = A y
-      subtract(x, residual, residual); // r = x - A y
+      auto R = residual.Read();
+      auto Y = y.ReadWrite();
+      // y += D^{-1} (x - A y)
+      mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
+      {
+         Y[i] += DI[i] * (X[i] - R[i]);
+      });
    }
    else
    {
-      residual = x;
-      y.UseDevice(true);
-      y = 0.0;
+      auto Y = y.Write();
+      // y = D^{-1} x
+      mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
+      {
+         Y[i] = DI[i] * X[i];
+      });
    }
-   auto DI = dinv.Read();
-   auto R = residual.Read();
-   auto Y = y.ReadWrite();
-   mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
-   {
-      Y[i] += DI[i] * R[i];
-   });
 }
 
 OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
@@ -348,7 +352,8 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
    diag(d),
    coeffs(order),
    ess_tdof_list(ess_tdofs),
-   residual(N),
+   residual(order > 1 ? N : 0),
+   z(order > 1 ? N : 0),
    oper(&oper_) { Setup(); }
 
 #ifdef MFEM_USE_MPI
@@ -375,7 +380,8 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
      diag(d),
      coeffs(order),
      ess_tdof_list(ess_tdofs),
-     residual(N),
+     residual(order > 1 ? N : 0),
+     z(order > 1 ? N : 0),
      oper(&oper_)
 {
    OperatorJacobiSmoother invDiagOperator(diag, ess_tdofs, 1.0);
@@ -422,7 +428,7 @@ void OperatorChebyshevSmoother::Setup()
 {
    // Invert diagonal
    residual.UseDevice(true);
-   helperVector.UseDevice(true);
+   z.UseDevice(true);
    auto D = diag.Read();
    auto X = dinv.Write();
    mfem::forall(N, [=] MFEM_HOST_DEVICE (int i) { X[i] = 1.0 / D[i]; });
@@ -432,6 +438,20 @@ void OperatorChebyshevSmoother::Setup()
       X[I[i]] = 1.0;
    });
 
+   const int order_save = order;
+   order = -1; // avoid early exit in SetOrder() when 'new_order' == 'order'
+   SetOrder(order_save);
+}
+
+void OperatorChebyshevSmoother::SetOrder(int new_order)
+{
+   if (new_order == order) { return; }
+
+   order = new_order;
+   coeffs.SetSize(order);
+   residual.SetSize(order > 1 ? N : 0);
+   z.SetSize(order > 1 ? N : 0);
+
    // Set up Chebyshev coefficients
    // For reference, see e.g., Parallel multigrid smoothing: polynomial versus
    // Gauss-Seidel by Adams et al.
@@ -511,32 +531,55 @@ void OperatorChebyshevSmoother::Mult(const Vector& x, Vector &y) const
       MFEM_ABORT("Chebyshev smoother requires operator");
    }
 
-   residual = x;
-   helperVector.SetSize(x.Size());
-   helperVector.UseDevice(true);
-
-   y.UseDevice(true);
-   y = 0.0;
-
-   for (int k = 0; k < order; ++k)
+   // for k = 0, perform:
+   //    r = D^{-1} x
+   //    y = C_0 r
+   const real_t C_0 = coeffs[0];
+   auto Dinv = dinv.Read();
+   auto X = x.Read();
+   auto Y0 = y.Write();
+   if (order == 1)
    {
-      // Apply
-      if (k > 0)
+      mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
       {
-         oper->Mult(residual, helperVector);
-         residual = helperVector;
-      }
+         Y0[i] = C_0 * Dinv[i] * X[i];
+      });
+   }
+   else
+   {
+      auto R0 = residual.Write();
+      mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
+      {
+         Y0[i] = C_0 * (R0[i] = Dinv[i] * X[i]);
+      });
+   }
 
-      // Scale residual by inverse diagonal
-      const int n = N;
-      auto Dinv = dinv.Read();
-      auto R = residual.ReadWrite();
-      mfem::forall(n, [=] MFEM_HOST_DEVICE (int i) { R[i] *= Dinv[i]; });
+   for (int k = 1; k < order; ++k)
+   {
+      // Apply: z = A r
+      oper->Mult(residual, z);
 
-      // Add weighted contribution to y
+      // Scale residual by inverse diagonal and add weighted contribution to y:
+      //   r = D^{-1} z
+      //   y += C_k r
+      const real_t C_k = coeffs[k];
+      auto Z = z.Read();
       auto Y = y.ReadWrite();
-      auto C = coeffs.Read();
-      mfem::forall(n, [=] MFEM_HOST_DEVICE (int i) { Y[i] += C[k] * R[i]; });
+      if (k < order-1)
+      {
+         auto R = residual.Write();
+         mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
+         {
+            Y[i] += C_k * (R[i] = Dinv[i] * Z[i]);
+         });
+      }
+      else
+      {
+         mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
+         {
+            Y[i] += C_k * Dinv[i] * Z[i];
+         });
+      }
    }
 }
 
@@ -3213,7 +3256,7 @@ void ResidualBCMonitor::MonitorResidual(
    MPI_Comm comm = iter_solver->GetComm();
    if (comm != MPI_COMM_NULL)
    {
-      double glob_bc_norm_squared = 0.0;
+      real_t glob_bc_norm_squared = 0.0;
       MPI_Reduce(&bc_norm_squared, &glob_bc_norm_squared, 1,
                  MPITypeMap<real_t>::mpi_type,
                  MPI_SUM, 0, comm);

linalg/solvers.hpp

@@ -380,11 +380,11 @@ public:
    void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
 
    /// Approach the solution of the linear system by applying Jacobi smoothing.
-   void Mult(const Vector &x, Vector &y) const;
+   void Mult(const Vector &x, Vector &y) const override;
 
    /** @brief Approach the solution of the transposed linear system by applying
        Jacobi smoothing. */
-   void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
+   void MultTranspose(const Vector &x, Vector &y) const override { Mult(x, y); }
 
    /** @brief Recompute the diagonal using the method AssembleDiagonal of the
        given new Operator, @a op. */
@@ -397,7 +397,7 @@ public:
        When the new Operator, @a op, is not a (Par)BilinearForm, any previously
        set array of essential true-dofs will be thrown away because in this case
        any essential b.c. will be handled by the AssembleDiagonal method. */
-   void SetOperator(const Operator &op);
+   void SetOperator(const Operator &op) override;
 
 private:
    Vector dinv;
@@ -481,21 +481,22 @@ public:
 
    /** @brief Approach the solution of the linear system by applying Chebyshev
        smoothing. */
-   void Mult(const Vector &x, Vector &y) const;
+   void Mult(const Vector &x, Vector &y) const override;
 
    /** @brief Approach the solution of the transposed linear system by applying
        Chebyshev smoothing. */
-   void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
+   void MultTranspose(const Vector &x, Vector &y) const override { Mult(x, y); }
 
-   void SetOperator(const Operator &op_)
+   void SetOperator(const Operator &op_) override
    {
       oper = &op_;
    }
 
    void Setup();
+   void SetOrder(int new_order);
 
 private:
-   const int order;
+   int order;
    real_t max_eig_estimate;
    const int N;
    Vector dinv;
@@ -503,7 +504,7 @@ private:
    Array<real_t> coeffs;
    const Array<int>& ess_tdof_list;
    mutable Vector residual;
-   mutable Vector helperVector;
+   mutable Vector z;
    const Operator* oper;
 };
 

makefile

@@ -125,7 +125,8 @@ EXAMPLE_TEST_DIRS := examples
 
 MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
- hooke multidomain dpg hdiv-linear-solver spde diag-smoothers
+ hooke multidomain dpg hdiv-linear-solver spde diag-smoothers \
+ benchmarks/ceed-solver-bps
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
@@ -145,7 +146,7 @@ MFEM_BUILD_DIR ?= .
 BUILD_DIR := $(MFEM_BUILD_DIR)
 BUILD_REAL_DIR := $(abspath $(BUILD_DIR))
 ifneq ($(BUILD_REAL_DIR),$(MFEM_REAL_DIR))
-   BUILD_SUBDIRS = $(DIRS) config $(EM_DIRS) doc $(TEST_DIRS)
+   BUILD_SUBDIRS = $(DIRS) config $(EM_DIRS) doc $(TEST_DIRS) tests/benchmarks
    CONFIG_FILE_DEF = -DMFEM_CONFIG_FILE='"$(BUILD_REAL_DIR)/config/_config.hpp"'
    BLD := $(if $(BUILD_REAL_DIR:$(CURDIR)=),$(BUILD_DIR)/,)
    $(if $(word 2,$(BLD)),$(error Spaces in BLD = "$(BLD)" are not supported))
@@ -582,13 +583,14 @@ test-print:
 	@for dir in $(ALL_TEST_DIRS); do \
 	   $(MAKE) -j1 -C $(BLD)$${dir} test-print; done
 
-ALL_CLEAN_SUBDIRS = $(addsuffix /clean,config $(EM_DIRS) doc $(TEST_DIRS))
+ALL_CLEAN_SUBDIRS = $(addsuffix /clean,config $(EM_DIRS) doc $(TEST_DIRS) \
+ tests/benchmarks)
 .PHONY: $(ALL_CLEAN_SUBDIRS) miniapps/clean
 miniapps/clean: $(addsuffix /clean,$(MINIAPP_DIRS))
 $(ALL_CLEAN_SUBDIRS):
 	$(MAKE) -C $(BLD)$(@D) $(@F)
 
-clean: $(addsuffix /clean,$(EM_DIRS) $(TEST_DIRS))
+clean: $(addsuffix /clean,$(EM_DIRS) $(TEST_DIRS) tests/benchmarks)
 	rm -f $(addprefix $(BLD),$(foreach d,$(DIRS),$(d)/*.o))
 	rm -f $(addprefix $(BLD),$(foreach d,$(DIRS),$(d)/*~))
 	rm -rf $(addprefix $(BLD),*~ libmfem.* deps.mk)
@@ -691,7 +693,8 @@ local-config:
 .PHONY: build-config
 build-config:
 	for d in $(BUILD_SUBDIRS); do mkdir -p $(BLD)$${d}; done
-	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS)); do \
+	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS) \
+                        tests/benchmarks); do \
 	   printf "# Auto-generated file.\n%s\n%s\n" \
 	      "MFEM_DIR = $(MFEM_REAL_DIR)" \
 	      "include \$$(MFEM_DIR)/$${dir}makefile" \

mesh/mesh.cpp

@@ -883,6 +883,8 @@ const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
                                                   const int flags,
                                                   MemoryType d_mt)
 {
+   MFEM_PERF_FUNCTION;
+
    for (int i = 0; i < geom_factors.Size(); i++)
    {
       GeometricFactors *gf = geom_factors[i];
@@ -14649,6 +14651,7 @@ GeometricFactors::GeometricFactors(const GridFunction &nodes,
 void GeometricFactors::Compute(const GridFunction &nodes,
                                MemoryType d_mt)
 {
+   MFEM_PERF_FUNCTION;
 
    const FiniteElementSpace *fespace = nodes.FESpace();
    const FiniteElement *fe = fespace->GetTypicalFE();

miniapps/benchmarks/ceed-solver-bps/kershaw.hpp

@@ -0,0 +1,156 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef __KERSHAW_HPP__
+#define __KERSHAW_HPP__
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+// 1D transformation at the right boundary.
+real_t right(const real_t eps, const real_t x)
+{
+   return (x <= 0.5) ? (2-eps) * x : 1 + eps*(x-1);
+}
+
+// 1D transformation at the left boundary
+real_t left(const real_t eps, const real_t x)
+{
+   return 1-right(eps,1-x);
+}
+
+// Transition from a value of "a" for x=0, to a value of "b" for x=1. Smoothness
+// is controlled by the parameter "s", taking values 0, 1, or 2.
+real_t step(const real_t a, const real_t b, real_t x, int s)
+{
+   if (x <= 0) { return a; }
+   if (x >= 1) { return b; }
+   switch (s)
+   {
+      case 0:
+      default:
+         return a + (b-a) * (x);
+      case 1: return a + (b-a) * (x*x*(3-2*x));
+      case 2: return a + (b-a) * (x*x*x*(x*(6*x-15)+10));
+   }
+}
+
+// 3D version of a generalized Kershaw mesh transformation, see D. Kershaw,
+// "Differencing of the diffusion equation in Lagrangian hydrodynamic codes",
+// JCP, 39:375–395, 1981.
+//
+// The input mesh should be Cartesian nx x ny x nz with nx divisible by 6 and
+// ny, nz divisible by 2.
+//
+// The eps parameters are in (0, 1]. Uniform mesh is recovered for epsy=epsz=1.
+void kershaw(const real_t epsy, const real_t epsz, const int smoothness,
+             const real_t x, const real_t y, const real_t z,
+             real_t &X, real_t &Y, real_t &Z)
+{
+   X = x;
+
+   int layer = x*6.0;
+   real_t lambda = (x-layer/6.0)*6;
+
+   // The x-range is split in 6 layers going from left-to-left, left-to-right,
+   // right-to-left (2 layers), left-to-right and right-to-right yz-faces.
+   switch (layer)
+   {
+      case 0:
+         Y = left(epsy, y);
+         Z = left(epsz, z);
+         break;
+      case 1:
+      case 4:
+         Y = step(left(epsy, y), right(epsy, y), lambda, smoothness);
+         Z = step(left(epsz, z), right(epsz, z), lambda, smoothness);
+         break;
+      case 2:
+         Y = step(right(epsy, y), left(epsy, y), lambda/2, smoothness);
+         Z = step(right(epsz, z), left(epsz, z), lambda/2, smoothness);
+         break;
+      case 3:
+         Y = step(right(epsy, y), left(epsy, y), (1+lambda)/2, smoothness);
+         Z = step(right(epsz, z), left(epsz, z), (1+lambda)/2, smoothness);
+         break;
+      default:
+         Y = right(epsy, y);
+         Z = right(epsz, z);
+         break;
+   }
+}
+
+struct KershawTransformation : VectorCoefficient
+{
+   real_t epsy, epsz;
+   int dim, s;
+   KershawTransformation(int dim_, real_t epsy_, real_t epsz_, int s_=0)
+      : VectorCoefficient(dim_), epsy(epsy_), epsz(epsz_), dim(dim_), s(s_) { }
+   using VectorCoefficient::Eval;
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override
+   {
+      real_t xyz[3];
+      Vector transip(xyz, 3);
+      T.Transform(ip, transip);
+      if (dim == 1)
+      {
+         V[0] = xyz[0]; // no transformation in 1D
+      }
+      else if (dim == 2)
+      {
+         real_t z=0, zt;
+         kershaw(epsy, epsz, s, xyz[0], xyz[1], z, V[0], V[1], zt);
+      }
+      else // dim == 3
+      {
+         kershaw(epsy, epsz, s, xyz[0], xyz[1], xyz[2], V[0], V[1], V[2]);
+      }
+   }
+};
+
+ParMesh CreateKershawMesh(int nx, int ny, int nz, real_t epsy, real_t epsz)
+{
+   const bool sfc_order = true;
+   Mesh serial_mesh;
+   if (nx > 0 && ny == 0 && nz == 0)
+   {
+      serial_mesh = Mesh::MakeCartesian1D(nx, 1.0);
+   }
+   else if (nx > 0 && ny > 0 && nz == 0)
+   {
+      serial_mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL,
+                                          false, 1, 1, sfc_order);
+   }
+   else if (nx > 0 && ny > 0 && nz > 0)
+   {
+      serial_mesh = Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON,
+                                          1, 1, 1, sfc_order);
+   }
+   else
+   {
+      MFEM_ABORT("Bad grid size");
+   }
+   KershawTransformation kt(serial_mesh.Dimension(), epsy, epsz);
+   serial_mesh.Transform(kt);
+   return ParMesh(MPI_COMM_WORLD, serial_mesh);
+}
+
+ParMesh CreateKershawMesh(int n, real_t eps)
+{
+   return CreateKershawMesh(n, n, n, eps, eps);
+}
+
+}
+
+#endif

miniapps/benchmarks/ceed-solver-bps/makefile

@@ -0,0 +1,77 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../../..
+MFEM_BUILD_DIR ?= ../../..
+MFEM_INSTALL_DIR ?= ../../../mfem
+SRC = $(if $(MFEM_DIR:../../..=),$(MFEM_DIR)/miniapps/benchmarks/ceed-solver-bps/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS =
+PAR_MINIAPPS = solver-bp
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+EXTRA_SOURCES = preconditioners.cpp
+EXTRA_HEADERS = kershaw.hpp rhs.hpp preconditioners.hpp
+EXTRA_OBJECTS = $(EXTRA_SOURCES:.cpp=.o)
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS)
+
+# Rule for building solver-bp
+solver-bp: solver-bp.o $(addprefix $(SRC),$(EXTRA_HEADERS)) \
+ $(EXTRA_OBJECTS) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(EXTRA_OBJECTS) $(MFEM_LIBS)
+
+# Rules for compiling *.o files
+# -I$(MFEM_DIR) is needed for "general/forall.hpp" for out-of-source builds
+%.o: $(SRC)%.cpp $(wildcard $(SRC)%.hpp) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -I$(MFEM_DIR) -c $(<) -o $(@)
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+solver-bp-test-par: solver-bp
+	@$(call mfem-test,$<, $(RUN_MPI), CEED Solver BP,,SKIP-NO-VIS)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS) $(EXTRA_OBJECTS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@true

miniapps/benchmarks/ceed-solver-bps/plot_csv.py

@@ -0,0 +1,129 @@
+import csv
+from pylab import *
+
+fields=[
+    ['code ID', 'str'],
+    ['preconditioner ID', 'str'],
+    ['machine ID', 'str'],
+    ['number of nodes', 'int'],
+    ['number of MPI ranks', 'int'],
+    ['n_x', 'int'], ['n_y', 'int'], ['n_z', 'int'],
+    ['solution polynomial degree', 'int'],
+    ['number of 1D quadrature points', 'float'],
+    ['eps_y', 'float'], ['eps_z', 'float'],
+    ['ndofs (including Dirichlet boundary)', 'int'],
+    ['niter', 'int'],
+    ['initial residual', 'float'], ['final residual', 'float'],
+    ['error', 'float'],
+    ['t_setup (preconditioner setup)', 'float'],
+    ['t_solve (total iter time)', 'float']]
+fields_dict=dict(fields)
+
+def convert(obj, type_str):
+    ctor=getattr(__builtins__, type_str)
+    return ctor(obj)
+
+input_csv='run-001.csv'
+print('reading %s ...' % input_csv)
+runs = []
+with open(input_csv) as csvfile:
+    csvreader = csv.DictReader(csvfile, fieldnames=[f[0] for f in fields],
+                               restkey='additional notes')
+    for row in csvreader:
+        for i in fields_dict:
+            row[i]=convert(row[i], fields_dict[i])
+        runs.append(row)
+
+orders=[r['solution polynomial degree'] for r in runs]
+orders=unique(orders) # numpy function
+# orders=[1]
+
+nps=[r['number of MPI ranks'] for r in runs]
+nps=unique(nps)
+if len(nps) > 1:
+    print('multiple num-ranks present: %s' % nps)
+    quit()
+np=nps[0]
+
+# plot fx (or fx/fn) vs fy, (or fx/fn/fy, etc) for all orders
+fn='number of MPI ranks'
+fx='ndofs (including Dirichlet boundary)'
+fy='t_solve (total iter time)'
+# fy='niter'
+# fy='error'
+fz='niter'
+
+figure()
+for p in orders:
+    rr=[r for r in runs if (r['solution polynomial degree']==p and
+                            r['niter']>0)]
+    if len(rr)==0:
+        continue
+
+    # pl_data=asarray([[r[fx],r[fx]/r[fy]] for r in rr])
+    # pl_data=asarray([[r[fx],r[fy]] for r in rr])
+    # pl_data=asarray([[r[fx],r[fx]/(r[fy]/r[fz])] for r in rr])
+
+    pl_data=asarray([[r[fx]/r[fn],r[fx]/r[fn]/r[fy]] for r in rr])
+    # pl_data=asarray([[r[fx]/r[fn],r[fy]] for r in rr])
+
+    plot(pl_data[:,0],pl_data[:,1], 'o-', label='p=%i'%p)
+    rnx=asarray([r['n_x'] for r in rr])
+    rerr=asarray([r['error'] for r in rr])
+    rate=arange(1.0,len(rnx))
+    for l in range(1,len(rnx)):
+        rate[l-1]=log(rerr[l-1]/rerr[l])/log(rnx[l]/rnx[l-1])
+    set_printoptions(formatter={'float':"{:6.2f}".format},linewidth=120)
+    print(f"p={p} rate:{rate}")
+
+# xscale('log', basex=10) # older matplotlib
+xscale('log', base=10)
+# xlim(4e4,3.1e7)
+xlim(4e4,5e6)
+# yscale('log', basey=10) # older matplotlib
+# yscale('log', base=10)
+# ylim(1e5,2e7)
+# ylim(0,2.55e7)
+# ylim(0,3.25e7)
+# ylim(0,5e6)
+ymin,ymax=ylim()
+ylim(0,ymax)
+# ylim(1e-2,2e1)
+# ylim(3e-3,6e-2)
+# xlabel(fx)
+# xlabel('# DOFs')
+xlabel('# DOFs / # Ranks')
+# ylabel(fx + ' / ' + fy)
+# ylabel(fy)
+# ylabel('# DOFs / t_solve')
+ylabel('# DOFs / # Ranks / t_solve')
+# ylabel('t_solve')
+# ylabel('# DOFs / (t_solve / # Iter)')
+# ylabel('# Iter')
+# ylabel('L2 error')
+# ylabel('Grad L2 error')
+grid('on', color='gray', ls='dotted')
+grid('on', axis='both', which='minor', color='gray', ls='dotted')
+legend(ncol=2, loc='best')
+ranks='1 MPI rank'
+if np > 1:
+   ranks='%s MPI ranks' % (np,np)
+hypre='hypre CPU'
+# hypre='hypre HIP'
+# prec=hypre+', p-MG(1,1)'
+prec=hypre+', LOR'
+# prec='Jacobi'
+# eps='1'
+eps='0.3'
+mfem='MFEM CPU'
+# mfem='MFEM HIP'
+title(mfem + ', ' + prec + ', $\\varepsilon = ' + eps + '$, ' + ranks)
+
+if 1: # write .pdf file?
+    pdf_file='plot.pdf'
+    print('saving figure --> %s'%pdf_file)
+    savefig(pdf_file, format='pdf', bbox_inches='tight')
+
+if 0: # show the figures?
+    print('\nshowing figures ...')
+    show()

miniapps/benchmarks/ceed-solver-bps/preconditioners.cpp

@@ -0,0 +1,241 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "preconditioners.hpp"
+
+namespace mfem
+{
+
+AssemblyLevel GetCoarseAssemblyLevel(SolverConfig config)
+{
+   switch (config.type)
+   {
+      case SolverConfig::JACOBI:
+      case SolverConfig::LOR_HYPRE:
+      case SolverConfig::LOR_AMGX:
+         return AssemblyLevel::PARTIAL;
+      default:
+         return AssemblyLevel::FULL;
+         // return AssemblyLevel::LEGACYFULL;
+   }
+}
+
+bool NeedsLOR(SolverConfig config)
+{
+   switch (config.type)
+   {
+      case SolverConfig::LOR_HYPRE:
+      case SolverConfig::LOR_AMGX:
+         return true;
+      default:
+         return false;
+   }
+}
+
+DiffusionMultigrid::DiffusionMultigrid(
+   ParFiniteElementSpaceHierarchy& hierarchy,
+   Coefficient &coeff_,
+   Array<int>& ess_bdr,
+   SolverConfig coarse_solver_config,
+   int q1d_inc_,
+   int smoothers_cheby_order_)
+   : GeometricMultigrid(hierarchy, ess_bdr),
+     coeff(coeff_),
+     q1d_inc(q1d_inc_),
+     irs(0, Quadrature1D::GaussLegendre),
+     smoothers_cheby_order(smoothers_cheby_order_)
+{
+   ConstructCoarseOperatorAndSolver(
+      coarse_solver_config, hierarchy.GetFESpaceAtLevel(0), ess_bdr);
+   int nlevels = hierarchy.GetNumLevels();
+   for (int i=1; i<nlevels; ++i)
+   {
+      ConstructOperatorAndSmoother(hierarchy.GetFESpaceAtLevel(i), ess_bdr);
+   }
+}
+
+void DiffusionMultigrid::ConstructBilinearForm(
+   ParFiniteElementSpace &fespace, Array<int> &ess_bdr, AssemblyLevel asm_lvl)
+{
+   ParBilinearForm *form = new ParBilinearForm(&fespace);
+   form->SetAssemblyLevel(asm_lvl);
+
+   DiffusionIntegrator *integ = new DiffusionIntegrator(coeff);
+
+   int p = fespace.GetOrder(0);
+   int dim = fespace.GetMesh()->Dimension();
+   // Integration rule for high-order problem: (p+1+q1d_inc)^d Gauss-Legendre
+   // points
+   int int_order = 2*(p+1+q1d_inc) - 1;
+   Geometry::Type geom = fespace.GetMesh()->GetElementBaseGeometry(0);
+   const IntegrationRule &ir = irs.Get(geom, int_order);
+   MFEM_VERIFY(ir.Size() == pow(p+1+q1d_inc,dim), "Wrong quadrature");
+   integ->SetIntegrationRule(ir);
+
+   form->AddDomainIntegrator(integ);
+   form->Assemble();
+   bfs.Append(form);
+
+   essentialTrueDofs.Append(new Array<int>());
+   fespace.GetEssentialTrueDofs(ess_bdr, *essentialTrueDofs.Last());
+}
+
+void DiffusionMultigrid::ConstructOperatorAndSmoother(
+   ParFiniteElementSpace& fespace, Array<int>& ess_bdr)
+{
+   ConstructBilinearForm(fespace, ess_bdr, AssemblyLevel::PARTIAL);
+
+   OperatorPtr opr;
+   bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
+   opr.SetOperatorOwner(false);
+
+   Vector diag(fespace.GetTrueVSize());
+   bfs.Last()->AssembleDiagonal(diag);
+
+   Solver* smoother = new OperatorChebyshevSmoother(
+      *opr, diag, *essentialTrueDofs.Last(), smoothers_cheby_order,
+      fespace.GetParMesh()->GetComm());
+
+   AddLevel(opr.Ptr(), smoother, true, true);
+}
+
+void DiffusionMultigrid::ConstructCoarseOperatorAndSolver(
+   SolverConfig config, ParFiniteElementSpace& fespace, Array<int>& ess_bdr)
+{
+   ConstructBilinearForm(fespace, ess_bdr, GetCoarseAssemblyLevel(config));
+   ParBilinearForm &a = static_cast<ParBilinearForm&>(*bfs.Last());
+   Array<int> &ess_dofs = *essentialTrueDofs.Last();
+
+   a.FormSystemMatrix(ess_dofs, A_coarse);
+
+   OperatorPtr A_prec;
+   if (NeedsLOR(config))
+   {
+      if (Mpi::Root())
+      {
+         std::cout << "Forming LOR discretization..." << std::endl;
+      }
+      lor.reset(new ParLORDiscretization(a, ess_dofs));
+      A_prec = lor->GetAssembledSystem();
+      if (Mpi::Root())
+      {
+         std::cout << "Forming LOR discretization... Done." << std::endl;
+      }
+   }
+   else
+   {
+      A_prec = A_coarse;
+   }
+
+   if (Mpi::Root()) { std::cout << "Forming preconditioner... " << std::endl; }
+   switch (config.type)
+   {
+      case SolverConfig::JACOBI:
+         coarse_precond.reset(new OperatorJacobiSmoother(a, ess_dofs));
+         break;
+      case SolverConfig::FA_HYPRE:
+      case SolverConfig::LOR_HYPRE:
+      {
+         HypreBoomerAMG *amg = new HypreBoomerAMG(*A_prec.As<HypreParMatrix>());
+         amg->SetPrintLevel(1);
+         Vector b(amg->Height());
+         Vector x(amg->Height());
+         b = 0.0;
+         x = 0.0;
+         amg->Setup(b, x); // Force setup;
+         coarse_precond.reset(amg);
+         break;
+      }
+#ifdef MFEM_USE_AMGX
+      case SolverConfig::FA_AMGX:
+      case SolverConfig::LOR_AMGX:
+      {
+         AmgXSolver *amg = new AmgXSolver;
+         amg->ReadParameters(config.amgx_config_file, AmgXSolver::EXTERNAL);
+         amg->InitExclusiveGPU(MPI_COMM_WORLD);
+         amg->SetOperator(*A_prec.As<HypreParMatrix>());
+         coarse_precond.reset(amg);
+         break;
+      }
+#endif
+      default:
+         MFEM_ABORT("Not available.")
+   }
+
+   if (config.inner_sli) // coarse_solver = SLI
+   {
+      SLISolver *sli = new SLISolver(fespace.GetComm());
+      sli->SetPrintLevel(0);
+      sli->SetAbsTol(0.0);
+      sli->SetRelTol(0.0);
+      sli->SetMaxIter(config.inner_sli_iter);
+      sli->SetOperator(*A_coarse);
+      sli->SetPreconditioner(*coarse_precond);
+      coarse_solver.reset(sli);
+   }
+   else if (config.inner_cg)
+   {
+      CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
+      cg->SetPrintLevel(2);
+      cg->SetMaxIter(100);
+      cg->SetRelTol(1e-8);
+      cg->SetAbsTol(0.0);
+      cg->SetOperator(*A_coarse);
+      cg->SetPreconditioner(*coarse_precond);
+      cg->iterative_mode = false;
+      coarse_solver.reset(cg);
+   }
+   else
+   {
+      coarse_solver = coarse_precond;
+   }
+   if (Mpi::Root())
+   {
+      std::cout << "Forming preconditioner... Done.\n" << std::endl;
+   }
+
+   if (config.coarse_smooth)
+   {
+      Vector diag(fespace.GetTrueVSize());
+      a.AssembleDiagonal(diag);
+
+      Solver *smoother = new OperatorChebyshevSmoother(
+         *A_coarse, diag, ess_dofs, smoothers_cheby_order,
+         fespace.GetParMesh()->GetComm());
+
+      AddLevel(A_coarse.Ptr(), smoother, false, true);
+      AddCoarseSolver(coarse_solver.get(), false);
+   }
+   else
+   {
+      AddLevel(A_coarse.Ptr(), coarse_solver.get(), false, false);
+   }
+}
+
+void DiffusionMultigrid::SetSmoothersChebyshevOrder(int new_cheby_order)
+{
+   for (int level = MultigridBase::coarse_solver ? 0 : 1;
+        level < NumLevels(); level++)
+   {
+      OperatorChebyshevSmoother *cheby =
+         dynamic_cast<OperatorChebyshevSmoother*>(GetSmootherAtLevel(level));
+      if (cheby) { cheby->SetOrder(new_cheby_order); }
+   }
+   smoothers_cheby_order = new_cheby_order;
+}
+
+void DiffusionMultigrid::SetInnerSLINumIter(int inner_sli_iter)
+{
+   SLISolver *sli = dynamic_cast<SLISolver*>(coarse_solver.get());
+   if (sli) { sli->SetMaxIter(inner_sli_iter); }
+}
+
+} // namespace mfem

miniapps/benchmarks/ceed-solver-bps/preconditioners.hpp

@@ -0,0 +1,100 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef __SOLVER_BP_HPP__
+#define __SOLVER_BP_HPP__
+
+#include "mfem.hpp"
+#include <memory>
+
+namespace mfem
+{
+
+struct SolverConfig
+{
+   enum SolverType
+   {
+      JACOBI = 0,
+      FA_HYPRE = 1,
+      LOR_HYPRE = 2,
+      FA_AMGX = 3,
+      LOR_AMGX = 4
+   };
+   SolverType type;
+   const char *amgx_config_file = "amgx/amgx.json";
+   bool inner_cg = false; //<-- use inner CG iteration for coarse solver
+   bool inner_sli = false; //<-- use inner SLI iteration for coarse solver
+   int inner_sli_iter = 1; //<- number of iterations for the inner SLI solver
+   bool coarse_smooth = false; //<- enable level 0 smoothing
+   SolverConfig(SolverType type_) : type(type_) { }
+
+   void Print()
+   {
+      mfem::out << "Coarse solver: ";
+      switch (type)
+      {
+         case JACOBI: mfem::out << "Jacobi"; break;
+         case FA_HYPRE: mfem::out << "Hypre (full)"; break;
+         case LOR_HYPRE: mfem::out << "Hypre (LOR)"; break;
+         case FA_AMGX: mfem::out << "AmgX (full)"; break;
+         case LOR_AMGX: mfem::out << "AmgX (LOR)"; break;
+      }
+      mfem::out << std::endl;
+      // If inner_sli is true inner_cg is not used, see
+      // DiffusionMultigrid::ConstructCoarseOperatorAndSolver():
+      if (inner_sli) { inner_cg = false; }
+      mfem::out << "Inner CG:      "
+                << (inner_cg ? "On" : "Off")
+                << std::endl;
+      mfem::out << "Inner SLI:     " << (inner_sli ? "On" : "Off") << '\n';
+      mfem::out << "Coarse smooth: " << (coarse_smooth ? "On" : "Off") << '\n';
+   }
+};
+
+struct DiffusionMultigrid : GeometricMultigrid
+{
+   Coefficient &coeff;
+   int q1d_inc;
+   IntegrationRules irs;
+   std::unique_ptr<ParLORDiscretization> lor;
+   OperatorPtr A_coarse;
+   std::shared_ptr<Solver> coarse_solver, coarse_precond;
+   int smoothers_cheby_order;
+
+   DiffusionMultigrid(
+      ParFiniteElementSpaceHierarchy& hierarchy,
+      Coefficient &coeff_,
+      Array<int>& ess_bdr,
+      SolverConfig coarse_solver_config,
+      int q1d_inc_ = 0,
+      int smoothers_cheby_order_ = 1);
+
+   void ConstructBilinearForm(
+      ParFiniteElementSpace &fespace,
+      Array<int> &ess_bdr,
+      AssemblyLevel asm_lvl);
+
+   void ConstructOperatorAndSmoother(
+      ParFiniteElementSpace &fespace,
+      Array<int> &ess_bdr);
+
+   void ConstructCoarseOperatorAndSolver(
+      SolverConfig config,
+      ParFiniteElementSpace &fespace,
+      Array<int> &ess_bdr);
+
+   void SetSmoothersChebyshevOrder(int new_cheby_order);
+   void SetInnerSLINumIter(int inner_sli_iter);
+};
+
+} // namespace mfem
+
+#endif

miniapps/benchmarks/ceed-solver-bps/rhs.hpp

@@ -0,0 +1,334 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef __RHS_HPP__
+#define __RHS_HPP__
+
+#include "mfem.hpp"
+#include "general/forall.hpp"
+
+// 0 - Solution described in the CEED MS 36 report
+// 1 - Solution from the "ecp_special_2023" paper (option with cosine):
+//       w(n,x) = \sum_{k=0}^n a^k \cos(b^k \pi (x - 1/2)),   x \in [0,1]
+//     with a = 1/2, b = 3.
+// 2 - Solution from the "ecp_special_2023" paper (option with sine):
+//       w(n,x) = \sum_{k=0}^n a^k \sin(b^k \pi x),   x \in [0,1]
+//     with a = 1/2, b = 3.
+#define CEED_SOLVER_BP_SOLUTION_OPTION 1
+
+namespace mfem
+{
+
+constexpr real_t pi = real_t(M_PI);
+
+#if (CEED_SOLVER_BP_SOLUTION_OPTION == 0)
+
+MFEM_HOST_DEVICE inline
+real_t s(int k, real_t x)
+{
+   return sin(2*pi*k*x);
+}
+
+MFEM_HOST_DEVICE inline
+real_t u(int k, real_t x)
+{
+   real_t skx = s(k,x);
+   real_t sgn = skx < 0 ? -1.0 : 1.0;
+   return exp(-1/skx/skx)*sgn;
+}
+
+MFEM_HOST_DEVICE inline
+real_t u_xx(int k, real_t x)
+{
+   real_t kpix = k*pi*x;
+   real_t csc_2kpix = 1.0/sin(2*kpix);
+   real_t sgn = sin(2*kpix) < 0 ? -1.0 : 1.0;
+   return 2*exp(-csc_2kpix*csc_2kpix)*k*k*pi*pi
+          *(1 + 6*cos(4*kpix) + cos(8*kpix))
+          *pow(csc_2kpix,6)
+          *sgn;
+}
+
+MFEM_HOST_DEVICE inline
+real_t w(int n, real_t x)
+{
+   real_t wkx = 0.0;
+   real_t xx = 2*x - 1; // transform from [0,1] to [-1,1]
+   for (int j=0; j<n; ++j)
+   {
+      int k = pow(3, j);
+      wkx += u(k, xx);
+   }
+   return wkx;
+}
+
+MFEM_HOST_DEVICE inline
+real_t w_xx(int n, real_t x)
+{
+   real_t wkx = 0.0;
+   real_t xx = 2*x - 1; // transform from [0,1] to [-1,1]
+   if (xx == 0.0) { return 0.0; }
+   for (int j=0; j<n; ++j)
+   {
+      int k = pow(3, j);
+      wkx += 4*u_xx(k, xx); // factor of four from reference interval transf.
+   }
+   return wkx;
+}
+
+#elif (CEED_SOLVER_BP_SOLUTION_OPTION == 1)
+
+MFEM_HOST_DEVICE inline
+real_t w(int n, real_t x)
+{
+   //  w(n,x) = \sum_{k=0}^n a^k \cos(b^k \pi (x - 1/2))
+   const real_t a = 0.5, b = 3.;
+   real_t ak = 1.0;
+   real_t xk = pi * (x - 0.5);
+   real_t w_ = ak * cos(xk);
+   for (int k = 1; k <= n; k++)
+   {
+      ak *= a;
+      xk *= b;
+      w_ += ak * cos(xk);
+   }
+   return w_;
+}
+
+MFEM_HOST_DEVICE inline
+real_t w_x(int n, real_t x)
+{
+   //  w'(n,x) = -\pi \sum_{k=0}^n a^k b^k \sin(b^k \pi (x - 1/2))
+   const real_t a = 0.5, b = 3.;
+   real_t ck = -pi;
+   real_t xk = pi * (x - 0.5);
+   real_t w_x_ = ck * sin(xk);
+   for (int k = 1; k <= n; k++)
+   {
+      ck *= a * b;
+      xk *= b;
+      w_x_ += ck * sin(xk);
+   }
+   return w_x_;
+}
+
+MFEM_HOST_DEVICE inline
+real_t w_xx(int n, real_t x)
+{
+   // w''(n,x) = -\pi^2 \sum_{k=0}^n a^k b^{2 k} \cos(b^k \pi (x - 1/2))
+   const real_t a = 0.5, b = 3.;
+   real_t ck = -(pi * pi);
+   real_t xk = pi * (x - 0.5);
+   real_t w_xx_ = ck * cos(xk);
+   for (int k = 1; k <= n; k++)
+   {
+      ck *= a * b*b;
+      xk *= b;
+      w_xx_ += ck * cos(xk);
+   }
+   return w_xx_;
+}
+
+#elif (CEED_SOLVER_BP_SOLUTION_OPTION == 2)
+
+MFEM_HOST_DEVICE inline
+real_t w(int n, real_t x)
+{
+   //  w(n,x) = \sum_{k=0}^n a^k \sin(b^k \pi x)
+   const real_t a = 0.5, b = 3.;
+   real_t ak = 1.0;
+   real_t xk = pi * x;
+   real_t w_ = ak * sin(xk);
+   for (int k = 1; k <= n; k++)
+   {
+      ak *= a;
+      xk *= b;
+      w_ += ak * sin(xk);
+   }
+   return w_;
+}
+
+MFEM_HOST_DEVICE inline
+real_t w_xx(int n, real_t x)
+{
+   // w''(n,x) = -\pi^2 \sum_{k=0}^n a^k b^{2 k} \sin(b^k \pi x)
+   const real_t a = 0.5, b = 3.;
+   real_t ck = -(pi * pi);
+   real_t xk = pi * x;
+   real_t w_xx_ = ck * sin(xk);
+   for (int k = 1; k <= n; k++)
+   {
+      ck *= a * b*b;
+      xk *= b;
+      w_xx_ += ck * sin(xk);
+   }
+   return w_xx_;
+}
+
+#endif // CEED_SOLVER_BP_SOLUTION_OPTION
+
+struct ExactSolution : Coefficient
+{
+   int dim, n;
+   ExactSolution(int dim_, int n_=0) : dim(dim_), n(n_) { }
+   using Coefficient::Eval;
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      real_t xyz[3];
+      Vector transip(xyz, 3);
+      T.Transform(ip, transip);
+      if (dim == 1)
+      {
+         return w(n, xyz[0]);
+      }
+      if (dim == 2)
+      {
+         return w(n, xyz[0])*w(n, xyz[1]);
+      }
+      else // dim == 3
+      {
+         return w(n, xyz[0])*w(n, xyz[1])*w(n, xyz[2]);
+      }
+   }
+};
+
+struct ExactGrad : VectorCoefficient
+{
+   int dim, n;
+   ExactGrad(int dim_, int n_)
+      : VectorCoefficient(dim_), dim(dim_), n(n_) { }
+   using VectorCoefficient::Eval;
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override
+   {
+      real_t xyz[3];
+      Vector transip(xyz, 3);
+      T.Transform(ip, transip);
+      V.SetSize(dim);
+      if (dim == 1)
+      {
+         V(0) = w_x(n, xyz[0]);
+      }
+      if (dim == 2)
+      {
+         V(0) = w_x(n, xyz[0])*  w(n, xyz[1]);
+         V(1) =   w(n, xyz[0])*w_x(n, xyz[1]);
+      }
+      else // dim == 3
+      {
+         const real_t wnx = w(n, xyz[0]);
+         const real_t wny = w(n, xyz[1]);
+         const real_t wnz = w(n, xyz[2]);
+         V(0) = w_x(n, xyz[0])*wny           *wnz;
+         V(1) = wnx           *w_x(n, xyz[1])*wnz;
+         V(2) = wnx           *wny           *w_x(n, xyz[2]);
+      }
+   }
+};
+
+MFEM_HOST_DEVICE inline
+real_t rhs_1d(const int n, const real_t *xyz)
+{
+   return -w_xx(n, xyz[0]);
+}
+
+MFEM_HOST_DEVICE inline
+real_t rhs_2d(const int n, const real_t *xyz)
+{
+   return -w_xx(n, xyz[0])*w(n, xyz[1]) - w(n, xyz[0])*w_xx(n, xyz[1]);
+}
+
+MFEM_HOST_DEVICE inline
+real_t rhs_3d(const int n, const real_t *xyz)
+{
+   return -w_xx(n, xyz[0])*w(n, xyz[1])*w(n, xyz[2])
+          - w(n, xyz[0])*w_xx(n, xyz[1])*w(n, xyz[2])
+          - w(n, xyz[0])*w(n, xyz[1])*w_xx(n, xyz[2]);
+}
+
+using RHSFunctionType = real_t(*)(int dim, const real_t *xyz);
+
+template <RHSFunctionType F>
+void ProjectRHS_(int n, QuadratureFunction &qf)
+{
+   QuadratureSpaceBase &qs = *qf.GetSpace();
+   Mesh &mesh = *qs.GetMesh();
+   const IntegrationRule &ir = qs.GetIntRule(0);
+
+   auto *geom = mesh.GetGeometricFactors(ir, GeometricFactors::COORDINATES);
+
+   const int dim = qs.GetMesh()->Dimension();
+   const int nq = ir.Size();
+   const int N = qf.Size();
+
+   const real_t *d_x = geom->X.Read();
+   real_t *d_q = qf.Write();
+
+   mfem::forall(N, [=] MFEM_HOST_DEVICE (int ii)
+   {
+      const int i = ii / nq;
+      const int j = ii % nq;
+      real_t xvec[3];
+      for (int d = 0; d < dim; ++d)
+      {
+         xvec[d] = d_x[j + d*nq + i*dim*nq];
+      }
+      d_q[ii] = F(n, xvec);
+   });
+}
+
+void ProjectRHS(int n, QuadratureFunction &qf)
+{
+   const int dim = qf.GetSpace()->GetMesh()->Dimension();
+   switch (dim)
+   {
+      case 1: ProjectRHS_<rhs_1d>(n, qf); break;
+      case 2: ProjectRHS_<rhs_2d>(n, qf); break;
+      case 3: ProjectRHS_<rhs_3d>(n, qf); break;
+      default: MFEM_ABORT("Unsupported dimension.");
+   }
+}
+
+struct RHS : Coefficient
+{
+   int dim, n;
+   RHS(int dim_, int n_=0) : dim(dim_), n(n_) { }
+   using Coefficient::Eval;
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      real_t xyz[3];
+      Vector transip(xyz, 3);
+      T.Transform(ip, transip);
+      if (dim == 1)
+      {
+         return -w_xx(n, xyz[0]);
+      }
+      if (dim == 2)
+      {
+         return -w_xx(n, xyz[0])*w(n, xyz[1]) - w(n, xyz[0])*w_xx(n, xyz[1]);
+      }
+      else // dim == 3
+      {
+         return -w_xx(n, xyz[0])*w(n, xyz[1])*w(n, xyz[2])
+                - w(n, xyz[0])*w_xx(n, xyz[1])*w(n, xyz[2])
+                - w(n, xyz[0])*w(n, xyz[1])*w_xx(n, xyz[2]);
+      }
+   }
+
+   void Project(QuadratureFunction &qf) override
+   {
+      ProjectRHS(n,qf);
+   }
+};
+
+}
+
+#endif

miniapps/benchmarks/ceed-solver-bps/run.sh

@@ -0,0 +1,124 @@
+bsep="============================================================"
+ssep="----------------------------------------"
+# Enable GPU-aware MPI:
+# gpu_aware_mpi_env_cmd="env MPICH_GPU_SUPPORT_ENABLED=1"
+# gpu_aware_mpi="-g"
+# number of nodes, number of MPI ranks:
+nnodes=1
+np=1
+# dev="-d hip ${gpu_aware_mpi}"
+eps="0.3"
+# mpirun_np="mpirun -np"
+mpirun_np="env MFEM_REPORT_KERNELS=1 mpirun -np"
+# mpirun_np="${gpu_aware_mpi_env_cmd} flux run --exclusive -N ${nnodes} -n"
+# dry run:
+# mpirun_np="echo ${mpirun_np}"
+#    p-MG/LOR + FA-hypre, or diagonal (Jacobi smoother)
+#    prec_type: "p-mg", "lor", or "diag"
+prec_type="lor"
+p_mg_opts="-cb 1"
+# p_mg_opts="-cb 5 -sli -sli-it 6"
+# lor_opts="-cls -cb 5 -sli -sli-it 6"
+# lor_opts="-cls -cb 2 -sli -sli-it 2"
+lor_opts="-cb 2 -sli -sli-it 2"
+mg_set=("1" "1 2" "1 3" "1 2 4" "1 3 5" "1 3 6")
+# mg_set=("1 2")
+#    p=7 and p=8 fail at the moment: "1 3 5 7" "1 3 5 8"
+# per-rank limits on the number of LOR elements for different p, in 2^20 units:
+# (bigger sizes run out of GPU memory, at least with LOR prec.)
+lor_ne_max_all=(4 4 4 4 4 4 4 4)
+# lor_ne_max_all=(18 22 24 24 27 24 8 8) # MI250X
+((lor_ne_min = 40*2**10))
+((np_ = np))
+((mm = 1))
+while ((np_ > 8)); do
+   ((mm++))
+   ((np_ = (np_-1)/8+1))
+done
+((mf = 2**mm))
+((mff = 3*mf))
+echo " *** mf = ${mf}, mff = ${mff}"
+for mg in "${mg_set[@]}"; do
+   echo "${bsep}"
+   p=(${mg})
+   # p=${p[-1]}
+   p="${p[$((${#p[@]}-1))]}"
+   lor_ne_max="${lor_ne_max_all[$((p-1))]}"
+   ((lor_ne_max *= 2**20))
+   # n_max = floor(lor_ne_max^(1/3))
+   n_max=$(echo "a=e((1/3)*l(${np}*${lor_ne_max}));scale=0;a/1" | bc -l)
+   # for np*lor_ne_max=256^3, the above gives 255, so we adjust the result:
+   while (( (n_max+1)**3 <= np*lor_ne_max )); do
+      ((n_max++))
+   done
+   echo " *** p = ${p}, n_max = ${n_max}"
+   if (( n_max**3 > np*lor_ne_max )); then
+      echo "error: n_max^3 > np*lor_ne_max"
+      exit 1
+   fi
+   echo "${bsep}"
+   nx_set=()
+   for ((nx = (n_max/p/mff)*mff, last_nx = 2*nx; nx >= 6; nx -= mff)); do
+      ((last_ne = last_nx**3))
+      ((ne = nx**3))
+      ((lor_ne = (p*nx)**3))
+      if ((np*lor_ne_min > lor_ne)); then break; fi
+      if ((last_ne < ne*4/3)); then continue; fi
+      nx_set=("${nx}" "${nx_set[@]}")
+      ((ndofs = (p*nx+1)**3))
+      ((rhs_n=0))
+      while ((2*3**(rhs_n+1) <= p*nx)); do
+         ((rhs_n++))
+      done
+      # 2*3**rhs_n <= p*nx < 2*3**(rhs_n+1)
+      printf "np = ${np}, p = ${p}, nx = ${nx}, ndofs = ${ndofs}"
+      # rhs_n for eps = 1:
+      # printf ", rhs_n = ${rhs_n}"
+      printf "\n"
+      ((last_nx = nx))
+   done
+   for nx in "${nx_set[@]}"; do
+      # break;
+      if ((nx % mf != 0)); then
+         echo " *** internal error!"
+         exit 1
+      fi
+      ((rp = mm))
+      ((nx /= mf))
+      if false; then
+         # 0, 1, or 2 additional parallel refinements for 1, 8, or 64 ranks
+         ((np_=np))
+         while ((np_%8 == 0)); do
+            ((np_=np_/8))
+            ((rp++))
+         done
+      fi
+      ((ndofs = (p*nx*2**rp+1)**3))
+      echo "${bsep}"
+      echo "np = ${np}, p = ${p}, ndofs = ${ndofs}"
+      if [[ "$prec_type" == "p-mg" ]]; then
+         # p-MG
+         printf "$mpirun_np ${np} ./solver-bp ${dev}"
+         printf " -ey ${eps} -mg \"${mg}\" -cs 1 ${p_mg_opts}"
+         printf " -nx ${nx} -rp ${rp}\n"
+         echo "${ssep}"
+         $mpirun_np "${np}" ./solver-bp ${dev} \
+            -ey ${eps} -mg "${mg}" -cs 1 ${p_mg_opts} -nx "${nx}" -rp "${rp}"
+      elif [[ "$prec_type" == "lor" ]]; then
+         # LOR
+         printf "$mpirun_np ${np} ./solver-bp ${dev}"
+         printf " -ey ${eps} -mg \"${p}\" -cs 2 ${lor_opts}"
+         printf " -nx ${nx} -rp ${rp}\n"
+         echo "${ssep}"
+         $mpirun_np "${np}" ./solver-bp ${dev} \
+            -ey ${eps} -mg "${p}" -cs 2 ${lor_opts} -nx "${nx}" -rp "${rp}"
+      elif [[ "$prec_type" == "diag" ]]; then
+         # Diag
+         printf "$mpirun_np ${np} ./solver-bp ${dev}"
+         printf " -ey ${eps} -mg \"${p}\" -cs 0 -nx ${nx} -rp ${rp}\n"
+         echo "${ssep}"
+         $mpirun_np "${np}" ./solver-bp ${dev} \
+            -ey ${eps} -mg "${p}" -cs 0 -nx "${nx}" -rp "${rp}"
+      fi
+   done
+done

miniapps/benchmarks/ceed-solver-bps/solver-bp.cpp

@@ -0,0 +1,811 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// --------------------------------------------------------------
+//    MFEM Implementation of the CEED Solver Bake-off Problems
+// --------------------------------------------------------------
+//
+// Run a suite of benchmarks and view the results:
+//
+//   1. Edit 'run.sh' to adjust machine and size parameters.
+//   2. Run 'run.sh' redirecting output to a file, e.g.:
+//        bash run.sh > run-001.out
+//   3. Extract the CSV output:
+//        sed -n -e 's/^= CSV:\(.*\)$/\1/p' run-001.out > run-001.csv
+//   3. Edit the script 'plot_csv.py' set the name of your CSV file and,
+//      optionally, customize the plot it generates.
+//   4. Process the CSV file:
+//        python3 plot_csv.py
+//
+// Sample runs:
+//
+//   solver-bp -nx 6
+//   solver-bp -nx 6 -mg "1 2 3"
+//   solver-bp -nx 6 -mg "1 r r 2 3"
+//   solver-bp -nx 6 -rp 2 -mg 3 -cs 1
+//   solver-bp -nx 6 -rp 2 -mg 3 -cs 2
+//
+// Device sample runs:
+//
+//   solver-bp -d cuda -nx 6 -mg "1 r r 2 3" -cs 0
+//   solver-bp -d cuda -nx 6 -rp 2 -mg 3 -cs 3
+//   solver-bp -d cuda -nx 6 -rp 2 -mg 3 -cs 4
+//
+
+#include "mfem.hpp"
+#include "kershaw.hpp"
+#include "rhs.hpp"
+#include "preconditioners.hpp"
+#include <regex>
+#include <fem/integ/bilininteg_diffusion_kernels.hpp>
+
+using namespace std;
+using namespace mfem;
+
+struct MGRefinement
+{
+   enum Type { P_MG, H_MG };
+   const Type type;
+   const int order;
+   MGRefinement(Type type_, int order_) : type(type_), order(order_) { }
+   static MGRefinement p(int order_) { return MGRefinement(P_MG, order_); }
+   static MGRefinement h() { return MGRefinement(H_MG, 0); }
+};
+
+struct CGMonitor : IterativeSolverMonitor
+{
+   const real_t tol;
+   real_t initial_nrm, final_nrm, saved_nrm;
+   int final_it, saved_it;
+
+   CGMonitor(real_t tol_) : tol(tol_) { }
+
+   void MonitorResidual(int it, real_t norm, const Vector &r, bool final)
+   override
+   {
+      MFEM_CONTRACT_VAR(norm);
+      // Avoid recomputing the norm if it was already computed -- this method
+      // is called two times for the final iteration: once with final = false
+      // (possibly triggering the monitor convergence criterion) and a second
+      // time with final = true.
+      bool init_call = (it == 0 && !final);
+      const real_t nrm =
+         (!init_call && it == saved_it) ?
+         saved_nrm :
+         sqrt(InnerProduct(iter_solver->GetComm(), r, r));
+      if ((it == 0 || final) && Mpi::Root())
+      {
+         mfem::out << (final ? "Final" : "   Initial")
+                   << " l2 norm of residual: " << nrm << '\n';
+      }
+      if (init_call)
+      {
+         initial_nrm = nrm;
+         converged = false;
+         final_nrm = -1.0;
+         final_it = -1;
+      }
+      saved_nrm = nrm;
+      saved_it = it;
+      // Check for monitor-triggered convergence
+      converged = (nrm <= tol*initial_nrm);
+      if (final)
+      {
+         final_nrm = nrm;
+         final_it = it;
+      }
+      if (final && Mpi::Root())
+      {
+         mfem::out << "Final relative l2 residual: ";
+         if (initial_nrm == 0.0)
+         {
+            mfem::out << "N/A (initial norm is 0)" << endl;
+         }
+         else
+         {
+            const real_t rel_nrm = nrm/initial_nrm;
+            mfem::out << rel_nrm << '\n';
+            mfem::out << "Average l2 reduction factor: ";
+            if (it == 0) { mfem::out << "N/A"; }
+            else { mfem::out << pow(rel_nrm, 1.0/it); }
+            mfem::out << " [" << it << " iterations]" << endl;
+         }
+      }
+   }
+};
+
+void report_hypre_gpu_status(bool gpu_aware_mpi_requested);
+void report_env_vars();
+real_t verify_ess_bdr(const Vector &b, const Vector &x,
+                      const Array<int> &ess_tdof_list);
+
+template <typename T> void PrintPair(const string &name, T val)
+{
+   cout << setw(14) << left << name << val << '\n';
+}
+
+
+int main(int argc, char *argv[])
+{
+   DiffusionIntegrator::AddSpecialization<3,3,3>();
+   DiffusionIntegrator::AddSpecialization<3,4,4>();
+   DiffusionIntegrator::AddSpecialization<3,5,5>();
+   DiffusionIntegrator::AddSpecialization<3,6,6>();
+
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *device_config = "cpu";
+   bool gpu_aware_mpi = false;
+   int nx = 6, ny = -1, nz = -1;
+   int rhs_n = -1;
+   const char *mg_spec = "1";
+   int q1d_inc = 0; // num 1D qpts = p + 1 + q1d_inc
+   int smoothers_cheby_order = 1;
+   real_t epsy = 1.0, epsz = -1;
+   int ref_par = 0;
+   bool glvis = false;
+   bool paraview = false;
+   SolverConfig coarse_solver(SolverConfig::JACOBI);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&gpu_aware_mpi, "-g", "--gpu-aware-mpi", "-no-g",
+                  "--no-gpu-aware-mpi", "Enable GPU-aware MPI.");
+   args.AddOption(&mg_spec, "-mg", "--multigrid-spec",
+                  "Multigrid specification. See README for description.");
+   args.AddOption(&q1d_inc, "-qi", "--quadrature-points-increment",
+                  "Increment for the 1D quadrature points relative to p + 1");
+   args.AddOption(&smoothers_cheby_order, "-cb",
+                  "--smoothers-chebyshev-order",
+                  "Order of the Chebyshev smoothers for the multigrid.");
+   args.AddOption((int*)&coarse_solver.type, "-cs", "--coarse-solver-config",
+                  "Coarse solver configuration. 0: Jacobi, 1: FA-HYPRE, "
+                  "2: LOR-HYPRE, 3: FA-AMGX, 4: LOR-AMGX.");
+   args.AddOption(&coarse_solver.inner_cg, "-cg", "--inner-cg",
+                  "-no-cg", "--no-inner-cg",
+                  "Use inner CG iteration for the coarse solver.");
+   args.AddOption(&coarse_solver.inner_sli, "-sli", "--inner-sli",
+                  "-no-sli", "--no-inner-sli",
+                  "Use inner SLI iteration for the coarse solver.");
+   args.AddOption(&coarse_solver.inner_sli_iter, "-sli-it",
+                  "--inner-sli-iterations",
+                  "Number of iterations for the inner SLI solver.");
+   args.AddOption(&coarse_solver.coarse_smooth, "-cls", "--coarse-level-smooth",
+                  "-no-cls", "--no-coarse-level-smooth",
+                  "Use coarse smoothing in addition to the coarse solver.");
+   args.AddOption(&coarse_solver.amgx_config_file, "-amgx", "--amgx-config",
+                  "AmgX config JSON file.");
+   args.AddOption(&nx, "-nx", "--nx", "Number of elements in x direction.");
+   args.AddOption(&ny, "-ny", "--ny", "Number of elements in y direction.");
+   args.AddOption(&nz, "-nz", "--nz", "Number of elements in z direction.");
+   args.AddOption(&epsy, "-ey", "--epsy", "Kershaw parameter epsilon y.");
+   args.AddOption(&epsz, "-ez", "--epsz", "Kershaw parameter epsilon z.");
+   args.AddOption(&rhs_n, "-rn", "--rhs-n",
+                  "Parameter n in the RHS function; -1 for default.");
+   args.AddOption(&ref_par, "-rp", "--ref-par",
+                  "Number of uniform parallel refinements to perform.");
+   args.AddOption(&glvis, "-gv", "--glvis", "-no-gv", "--no-glvis",
+                  "Save the mesh and solution for GLVis visualization.");
+   args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Save data files for ParaView visualization.");
+   args.ParseCheck();
+
+   if (ny < 0) { ny = nx; }
+   if (nz < 0) { nz = nx; }
+   if (epsz < 0) { epsz = epsy; }
+   // rhs_n default is handled later
+
+   Device device(device_config);
+   device.SetGPUAwareMPI(gpu_aware_mpi);
+   if (Mpi::Root()) { device.Print(); }
+   // Report HYPRE's GPU config and GPU-aware MPI config. Terminates if
+   // GPU-aware MPI is requested but HYPRE's GPU-aware MPI support is disabled.
+   report_hypre_gpu_status(gpu_aware_mpi);
+   // Report environment variables like {CUDA,ROCR}_VISIBLE_DEVICES:
+   report_env_vars();
+
+   // Generate mesh
+   ParMesh mesh_coarse = CreateKershawMesh(nx, ny, nz, epsy, epsz);
+   const int dim = mesh_coarse.Dimension();
+   for (int i=0; i<ref_par; ++i) { mesh_coarse.UniformRefinement(); }
+
+   int coarse_order = 0, order = 0, h_ref = ref_par;
+   // Parse order specification
+   vector<MGRefinement> mg_refinements;
+   {
+      istringstream mg_stream(mg_spec);
+      string ref;
+      mg_stream >> coarse_order;
+      int prev_order = order = coarse_order;
+      if (Mpi::Root()) { cout << "\nCoarse order " << coarse_order << '\n'; }
+      while (mg_stream >> ref)
+      {
+         if (ref == "r")
+         {
+            if (Mpi::Root()) { cout << "h-MG uniform refinement\n"; }
+            mg_refinements.push_back(MGRefinement::h());
+            ++h_ref;
+         }
+         else
+         {
+            try { order = stoi(ref); }
+            catch (...)
+            {
+               MFEM_ABORT("Multigrid refinement must either be an integer or "
+                          "the character `r`");
+            }
+            if (Mpi::Root()) { cout << "p-MG order   " << order << '\n'; }
+            MFEM_VERIFY(order > 0, "Orders must be positive");
+            MFEM_VERIFY(order > prev_order, "Orders must be increasing");
+            mg_refinements.push_back(MGRefinement::p(order));
+            prev_order = order;
+         }
+      }
+   }
+
+   if (order == 1 && coarse_solver.type == SolverConfig::LOR_HYPRE)
+   {
+      // Using ~10^7 elements with p=1 overflows a Vector in the LOR setup.
+      // The Vector has size (3D): (p+1)^3 * 27 * num_elem_ho.
+      // In 3D, for p > 1, the overflow will happen around:
+      // - p=2: ~23.6 million dofs or 2,945,794 elements
+      // - p=3: ~33.6 million dofs or 1,242,757 elements
+      // - p=4: ~40.7 million dofs or   636,292 elements
+      // - p=5: ~46.0 million dofs or   368,225 elements
+      // - p=6: ~50.1 million dofs or   231,885 elements
+      //
+      // Note: the size of the Jacobians at quadrature points (with q1d=p+1) in
+      // 3D is: (p+1)^3 * 9 * num_elem, so 3x smaller than the above Vector.
+      //
+      // For q1d=p+2, the overflow happens around:
+      // - p=1: 8,837,382 elements or ~8.8 million dofs
+      // - p=2: 3,728,271 elements or ~29.8 million dofs
+      // - p=3: 1,908,875 elements or ~51.5 million dofs
+      // - p=4: 1,104,673 elements or ~70.7 million dofs
+      // - p=5:   695,654 elements or ~87.0 million dofs
+      // - p=6:   466,034 elements or ~100.7 million dofs
+      coarse_solver.type = SolverConfig::FA_HYPRE;
+      if (Mpi::Root())
+      {
+         cout << "\nOrder is 1: switching from LOR-HYPRE to FA-HYPRE.\n";
+      }
+   }
+#if 0
+   if (order == 1 && coarse_solver.type == SolverConfig::FA_HYPRE &&
+       coarse_solver.inner_sli)
+   {
+      coarse_solver.inner_sli = false;
+      if (Mpi::Root())
+      {
+         cout << "\nOrder is 1: turning off the inner SLI.\n";
+      }
+   }
+#endif
+
+   vector<unique_ptr<FiniteElementCollection>> fe_collections;
+   fe_collections.emplace_back(new H1_FECollection(coarse_order, dim));
+   ParFiniteElementSpace fes_coarse(&mesh_coarse, fe_collections.back().get());
+   ParFiniteElementSpaceHierarchy hierarchy(&mesh_coarse, &fes_coarse,
+                                            false, false);
+
+   for (MGRefinement ref : mg_refinements)
+   {
+      if (ref.type == MGRefinement::H_MG)
+      {
+         hierarchy.AddUniformlyRefinedLevel();
+      }
+      else // P_MG
+      {
+         fe_collections.emplace_back(new H1_FECollection(ref.order, dim));
+         hierarchy.AddOrderRefinedLevel(fe_collections.back().get());
+      }
+   }
+
+   const int nlevels = hierarchy.GetNumLevels();
+   if (Mpi::Root())
+   {
+      if (nlevels == 1)
+      {
+         cout << "1 level in MG hierarchy. Using coarse solver only." << endl;
+      }
+      else
+      {
+         cout << nlevels << " levels in MG hierarchy." << endl;
+      }
+      coarse_solver.Print();
+      cout << endl;
+   }
+
+   // Determine final nx, ny, nz and use them to determine the default rhs_n.
+   const int ref_factor = pow(2, h_ref);
+   nx *= ref_factor;
+   ny *= ref_factor;
+   nz *= ref_factor;
+   if (rhs_n < 0)
+   {
+      int n_min = min(nx, ny);
+      if (nz > 0) { n_min = min(n_min, nz); }
+      // Find rhs_n such that 2*3^rhs_n <= (order*n_min) < 2*3^{rhs_n+1}
+      rhs_n = 0;
+      for (int l = 2*3; l <= order*n_min; l *= 3) { rhs_n++; }
+      if (epsy < 0.8) { rhs_n--; }
+      if (Mpi::Root()) { cout << "Using rhs_n = " << rhs_n << '\n' << endl; }
+   }
+
+   ParFiniteElementSpace &fes = hierarchy.GetFinestFESpace();
+   ParMesh &mesh = *fes.GetParMesh();
+   mesh.PrintInfo(cout);
+   HYPRE_Int ndof = fes.GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      cout << "\nTotal number of DOFs: " << ndof << endl << endl;
+   }
+
+   // All Dirichlet boundaries
+   Array<int> ess_bdr;
+   if (mesh.bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient coeff(1.0); // Diffusion coefficient
+   // Set up RHS
+   if (Mpi::Root()) { cout << "Assembling right-hand side..." << endl; }
+   RHS rhs_coeff(dim, rhs_n);
+   ParLinearForm b(&fes);
+   const int rhs_ir_inc = 2*q1d_inc+1;
+   // --> ir_order = 2*(p+1+q1d_inc)-1 --> q1d = p+1+q1d_inc
+   b.AddDomainIntegrator(new DomainLFIntegrator(rhs_coeff, 2, rhs_ir_inc));
+   b.UseFastAssembly(true);
+   b.Assemble();
+   if (Mpi::Root()) { cout << "Assembling right-hand side... Done." << endl; }
+
+   // make sure the GPU is done with any previous tasks:
+   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
+   // make sure all ranks are done with any previous tasks:
+   MPI_Barrier(MPI_COMM_WORLD);
+   tic();
+   // Set up operators in the multigrid hierarchy
+   DiffusionMultigrid MG(hierarchy, coeff, ess_bdr, coarse_solver, q1d_inc,
+                         smoothers_cheby_order);
+   MG.SetCycleType(Multigrid::CycleType::VCYCLE, 1, 1);
+   // make sure the GPU is done with all setup tasks:
+   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
+   // make sure all ranks are done with all setup tasks:
+   MPI_Barrier(MPI_COMM_WORLD);
+   const real_t t_setup = tic_toc.RealTime();
+
+   ParGridFunction x(&fes);
+   x = 0.0;
+
+   OperatorPtr A;
+   Vector X, B;
+   MG.FormFineLinearSystem(x, b, A, X, B);
+
+   const real_t l2_tol = 1e-8;
+   CGMonitor monitor(l2_tol);
+
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(0.0); // use the 'monitor' for convergence
+   cg.SetPrintLevel(1);
+   cg.SetOperator(*A);
+   cg.SetPreconditioner(MG);
+   cg.SetMonitor(monitor);
+   // Run 2 CG iterations to ensure everything is allocated and initialized for
+   // the full CG solve:
+   if (Mpi::Root()) { cout << "Running 2 warm-up CG iterations ...\n"; }
+   cg.SetMaxIter(2);
+   {
+      Vector X_save(X);
+      cg.Mult(B, X);
+      X = X_save;
+   }
+   if (coarse_solver.inner_sli &&
+       ((coarse_solver.type == SolverConfig::FA_HYPRE /* && order > 1 */) ||
+        coarse_solver.type == SolverConfig::LOR_HYPRE))
+   {
+      // timing data: (t-solve,sli-iter,cheby-order,pcg-iter)
+      std::vector<std::tuple<double,int,int,int>> timings;
+      Vector X_save(X);
+      if (Mpi::Root()) { cout << "\nFinding optimal MG parameters ...\n"; }
+      cg.SetMaxIter(500);
+      for (int sli_it = 1; sli_it <= coarse_solver.inner_sli_iter; sli_it++)
+      {
+         MG.SetInnerSLINumIter(sli_it);
+         for (int cheby_order = 1; cheby_order <= smoothers_cheby_order;
+              cheby_order++)
+         {
+            MG.SetSmoothersChebyshevOrder(cheby_order);
+
+            if (Mpi::Root())
+            {
+               cout << "\nRunning and timing parameters (sli iter, cheby order)"
+                    << " = (" << sli_it << ',' << cheby_order << ") ...\n";
+            }
+            // make sure the GPU is done with any previous tasks:
+            if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
+            // make sure all ranks are done with any previous tasks:
+            MPI_Barrier(MPI_COMM_WORLD);
+            tic();
+            cg.Mult(B, X);
+            // make sure the GPU is done with all solve tasks:
+            if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
+            // make sure all ranks are done with all solve tasks:
+            MPI_Barrier(MPI_COMM_WORLD);
+            const double t_solve = tic_toc.RealTime();
+            if (cg.GetConverged())
+            {
+               timings.emplace_back(t_solve, sli_it, cheby_order,
+                                    cg.GetNumIterations());
+            }
+            X = X_save;
+         }
+      }
+      std::sort(timings.begin(), timings.end());
+      if (Mpi::Root())
+      {
+         cout << "\nSorted timings from rank 0:\n";
+         const auto old_prec = cout.precision(6);
+         const auto old_fmtflags = cout.flags();
+         cout << std::fixed;
+         for (size_t i = 0; i < timings.size(); i++)
+         {
+            cout << setw(2) << i << ": "
+                 << 1e3*std::get<0>(timings[i]) << " ms: ("
+                 << std::get<1>(timings[i]) << ','
+                 << std::get<2>(timings[i]) << "): "
+                 << setw(3) << std::get<3>(timings[i]) << " iter\n";
+         }
+         cout.flags(old_fmtflags);
+         cout.precision(old_prec);
+      }
+      if (timings.size() > 0)
+      {
+         // Use the fastest parameters (as timed on rank 0) for the full solve:
+         int si = std::get<1>(timings[0]);
+         int co = std::get<2>(timings[0]);
+         MPI_Bcast(&si, 1, MPI_INT, 0, MPI_COMM_WORLD);
+         MPI_Bcast(&co, 1, MPI_INT, 0, MPI_COMM_WORLD);
+         MG.SetInnerSLINumIter(si);
+         MG.SetSmoothersChebyshevOrder(co);
+         coarse_solver.inner_sli_iter = si;
+         smoothers_cheby_order = co;
+         if (Mpi::Root())
+         {
+            cout << "\nUsing the fastest option (sli iter, cheby order) = ("
+                 << si << ',' << co << ")\n";
+         }
+      }
+      else
+      {
+         MG.SetInnerSLINumIter(1);
+         MG.SetSmoothersChebyshevOrder(1);
+         coarse_solver.inner_sli_iter = 1;
+         smoothers_cheby_order = 1;
+         if (Mpi::Root())
+         {
+            cout << "\nAll options failed to converge!"
+                 << " Using (sli iter, cheby order) = (1,1)\n";
+         }
+      }
+   }
+   if (Mpi::Root()) { cout << "\nRunning and timing the full CG solve ...\n"; }
+   cg.SetMaxIter(500);
+   // make sure the GPU is done with any previous tasks:
+   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
+   // make sure all ranks are done with any previous tasks:
+   MPI_Barrier(MPI_COMM_WORLD);
+   tic();
+   cg.Mult(B, X);
+   // make sure the GPU is done with all solve tasks:
+   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
+   // make sure all ranks are done with all solve tasks:
+   MPI_Barrier(MPI_COMM_WORLD);
+   const double t_solve = tic_toc.RealTime();
+
+   const int niter = cg.GetConverged() ? cg.GetNumIterations() : -1;
+
+   const real_t bdr_err = verify_ess_bdr(B, X, MG.GetFineEssentialTrueDofs());
+   if (Mpi::Root())
+   {
+      MFEM_VERIFY(bdr_err == 0.0, "Incorrect boundary values in solution!"
+                  " bdr_err = " << bdr_err);
+   }
+
+   MG.RecoverFineFEMSolution(X, b, x);
+
+   ExactSolution exact_coeff(dim, rhs_n);
+   // ExactGrad exact_grad_coeff(dim, rhs_n);
+   real_t L2_err = x.ComputeL2Error(exact_coeff);
+   // real_t grad_err = x.ComputeGradError(&exact_grad_coeff);
+   if (Mpi::Root())
+   {
+      cout << "\nL2 Error: " << setprecision(10) << scientific
+           << L2_err << '\n';
+      // cout << "\nGrad Error: " << setprecision(10) << scientific
+      //      << grad_err << '\n';
+   }
+
+   if (glvis)
+   {
+      ofstream mesh_ofs(MakeParFilename("mesh.", Mpi::WorldRank()));
+      mesh_ofs.precision(8);
+      mesh.Print(mesh_ofs);
+
+      ofstream sol_ofs(MakeParFilename("sol.", Mpi::WorldRank()));
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   if (paraview)
+   {
+      ParGridFunction rhs_gf(&fes), exact_gf(&fes), error_gf(&fes);
+      rhs_gf.ProjectCoefficient(rhs_coeff);
+      exact_gf.ProjectCoefficient(exact_coeff);
+
+      subtract(exact_gf, x, error_gf);
+
+      ParaViewDataCollection dc("SolverBP", &mesh);
+      dc.RegisterField("u", &x);
+      dc.RegisterField("rhs", &rhs_gf);
+      dc.RegisterField("exact", &exact_gf);
+      dc.RegisterField("error", &error_gf);
+      dc.SetPrefixPath("ParaView");
+      dc.SetLevelsOfDetail(order);
+      dc.SetHighOrderOutput(true);
+      dc.SetCycle(0);
+      dc.SetTime(0.0);
+      dc.Save();
+   }
+
+   const long long nel = mesh.GetGlobalNE();
+   if (nz == 0) { MFEM_VERIFY(nel == nx*ny, "Wrong number of elements"); }
+   else { MFEM_VERIFY(nel == nx*ny*nz, "Wrong number of elements"); }
+
+   if (Mpi::Root())
+   {
+      cout << "\n= Results\n";
+      PrintPair("nranks", Mpi::WorldSize());
+      PrintPair("nx", nx);
+      PrintPair("ny", ny);
+      PrintPair("nz", nz);
+      PrintPair("degree", order);
+      PrintPair("rhs_n", rhs_n);
+      PrintPair("epsy", epsy);
+      PrintPair("epsz", epsz);
+      PrintPair("ndof", ndof);
+      PrintPair("niter", niter);
+
+      // Should also output:
+      // code id
+      // prec id
+      // machine id
+      // number of supercomputer nodes
+      // number of 1d quadrature points
+      // initial and final residuals
+      // error
+
+      // Timings
+      PrintPair("t_setup", t_setup);
+      PrintPair("t_solve", t_solve);
+
+      cout << "\nSolve MDOFs/rank/sec:   "
+           << ndof/1e6/Mpi::WorldSize()/t_solve << '\n';
+
+      // CSV fields:
+      // 1. code ID
+      // 2. preconditioner ID
+      // 3. machine ID
+      // 4. number of nodes
+      // 5. number of MPI ranks
+      // 6,7,8. n_x, n_y, n_z
+      // 9. solution polynomial degree
+      // 10. number of 1D quadrature points
+      // 11,12. eps_y, eps_z
+      // 13. ndofs (including Dirichlet boundary)
+      // 14. niter
+      // 15,16. initial and final residuals
+      // 17. error
+      // 18. t_setup (preconditioner setup)
+      // 19. t_solve (total iter time)
+      //
+      // extract the CSV lines from the output with:
+      //    grep "= CSV:" out.txt | sed -e 's/^= CSV://' > out.csv
+      cout << "\n= CSV:"
+           << "MFEM-" + string(device_config); // 1
+      string hypre_str =
+#if defined(HYPRE_USING_HIP)
+         "hypre-hip"
+#elif defined(HYPRE_USING_CUDA)
+         "hypre-cuda"
+#else
+         "hypre-cpu"
+#endif
+         ;
+      auto cs = coarse_solver.type;
+      string prec_id;
+      if (cs == SolverConfig::FA_HYPRE) // p-MG, add (sli-iter,cheby-order)
+      {
+         prec_id = hypre_str + "-pMG(";
+      }
+      else if (cs == SolverConfig::LOR_HYPRE) // LOR, add (sli-iter,cheby-order)
+      {
+         prec_id = hypre_str + "-LOR(";
+      }
+      else if (cs == SolverConfig::JACOBI)
+      {
+         prec_id = "diag(";
+      }
+      else
+      {
+         prec_id = "(unknown)(";
+      }
+      if (coarse_solver.inner_cg)
+      {
+         prec_id += "cg;";
+      }
+      if (coarse_solver.inner_sli)
+      {
+         prec_id += to_string(coarse_solver.inner_sli_iter) + ";";
+      }
+      prec_id += to_string(smoothers_cheby_order) +
+                 (coarse_solver.coarse_smooth ? "c" : "") + ")";
+      prec_id += "-" + regex_replace(mg_spec, regex(" "), "-");
+      cout << ',' << prec_id; // 2
+      const char *hostname = getenv("HOSTNAME");
+      if (!hostname) { hostname = getenv("HOST"); }
+      string host_id = regex_replace(hostname ? hostname : "(unknown)",
+                                     regex("[0-9]*$"), "");
+      cout << ',' << host_id; // 3
+      cout << ',' << (fes.GetNRanks() + 7)/8; // 4 (assuming 8 ranks/node !!)
+      cout << ',' << fes.GetNRanks(); // 5
+      cout << ',' << nx << ',' << ny << ',' << nz; // 6,7,8
+      cout << ',' << order; // 9
+      // DiffusionMultigrid::ConstructBilinearForm p+1+q1d_inc 1D points
+      real_t Q1D = order + 1 + q1d_inc;
+      cout << ',' << defaultfloat << Q1D; // 10 (note: written as real_t)
+      cout << ',' << scientific << epsy << ',' << epsz; // 11,12
+      cout << ',' << ndof; // 13
+      cout << ',' << niter; // 14
+      cout << ',' << monitor.initial_nrm << ',' << monitor.final_nrm; // 15,16
+      cout << ',' << L2_err; // 17
+      // cout << ',' << grad_err; // 17 *** for testing ***
+      cout << ',' << t_setup << ',' << t_solve; // 18,19
+      cout << endl;
+   }
+
+   return 0;
+}
+
+void report_hypre_gpu_status(bool gpu_aware_mpi_requested)
+{
+#ifdef HYPRE_WITH_GPU_AWARE_MPI
+   bool hypre_gpu_aware_mpi = true;
+#else
+   bool hypre_gpu_aware_mpi = false;
+#endif
+#if (MFEM_HYPRE_VERSION > 23000)
+   hypre_gpu_aware_mpi = hypre_gpu_aware_mpi && hypre_GetGpuAwareMPI();
+#endif
+   if (Mpi::Root())
+   {
+      MFEM_VERIFY(!gpu_aware_mpi_requested || hypre_gpu_aware_mpi,
+                  "GPU-aware MPI requested but HYPRE's GPU-aware MPI support"
+                  " is not enabled");
+      cout << "\nHYPRE GPU support: "
+#ifdef HYPRE_USING_GPU
+           << "enabled";
+#else
+           << "disabled";
+#endif
+      cout << "\nHYPRE GPU-aware MPI support: "
+           << (hypre_gpu_aware_mpi ? "enabled" : "disabled") << endl;
+   }
+}
+
+void report_env_vars()
+{
+   const int myid = Mpi::WorldRank();
+   // const int lastid = min(Mpi::WorldSize(),4)-1; // show up to 4 ranks
+   const int lastid = Mpi::WorldSize()-1;
+   if (myid > lastid) { return; }
+   Array<char> recv_buf;
+   int buflen = -1, tag = 42;
+   const char *env_vars[] =
+   {
+      "HOST", "HOSTNAME", "MPICH_GPU_SUPPORT_ENABLED", "CUDA_VISIBLE_DEVICES",
+      "ROCR_VISIBLE_DEVICES"
+   };
+   const int num_env_vars = sizeof(env_vars)/sizeof(env_vars[0]);
+   // Send strings to rank 0, so that they can be printed in order, guaranteed.
+   // Every rank > 0 sends to rank 0:
+   if (myid > 0)
+   {
+      for (int ev = 0; ev < num_env_vars; ev++)
+      {
+         const char *env_var_val = getenv(env_vars[ev]);
+         buflen = env_var_val ? int(strlen(env_var_val)+1) : -1;
+         MPI_Send(&buflen, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
+         if (env_var_val)
+         {
+            MPI_Send(env_var_val, buflen, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
+         }
+      }
+   }
+   else // myid == 0
+   {
+      cout << "\nDefined environment variables:\n";
+      for (int id = 0; id <= lastid; id++)
+      {
+         cout << "[rank " << id << "]:";
+         for (int ev = 0, vars_shown = 0; ev < num_env_vars; ev++)
+         {
+            const char *env_var_val = nullptr;
+            if (id == 0)
+            {
+               env_var_val = getenv(env_vars[ev]);
+               buflen = env_var_val ? 0 : -1;
+            }
+            else
+            {
+               MPI_Recv(&buflen, 1, MPI_INT, id, tag, MPI_COMM_WORLD,
+                        MPI_STATUS_IGNORE);
+            }
+            if (buflen != -1)
+            {
+               if (id > 0)
+               {
+                  recv_buf.SetSize(buflen);
+                  MPI_Recv(recv_buf.begin(), buflen, MPI_CHAR, id, tag,
+                           MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                  env_var_val = recv_buf.begin();
+               }
+               if (vars_shown)
+               {
+                  cout << "\n[rank " << id << "]:";
+               }
+               cout << ' ' << env_vars[ev] << '=' << env_var_val;
+               vars_shown++;
+            }
+         }
+         cout << '\n';
+      }
+      if (lastid < Mpi::WorldSize()-1)
+      {
+         cout << "... [only " << lastid+1 << '/' << Mpi::WorldSize()
+              << " ranks shown]\n";
+      }
+      cout << flush;
+   }
+}
+
+real_t verify_ess_bdr(const Vector &b, const Vector &x,
+                      const Array<int> &ess_tdof_list)
+{
+   Vector d(ess_tdof_list.Size());
+   auto d_b = b.Read();
+   auto d_x = x.Read();
+   auto d_d = d.Write();
+   auto d_ess_ind = ess_tdof_list.Read();
+   mfem::forall(ess_tdof_list.Size(), [=] MFEM_HOST_DEVICE (int i)
+   {
+      const int ind = d_ess_ind[i];
+      d_d[i] = -fabs(d_b[ind] - d_x[ind]);
+   });
+   real_t d_max = -d.Min(); // max is not implemented on device
+   MPI_Allreduce(MPI_IN_PLACE, &d_max, 1, MFEM_MPI_REAL_T, MPI_MAX,
+                 MPI_COMM_WORLD);
+   return d_max;
+}

miniapps/solvers/block-solvers.cpp

@@ -365,6 +365,7 @@ int main(int argc, char *argv[])
    std::map<const DarcySolver*, real_t> setup_time;
    chrono.Restart();
    BDPMinresSolver bdp(M, B, param);
+   bdp.iterative_mode = true;
    setup_time[&bdp] = chrono.RealTime();
 
    chrono.Restart();

miniapps/solvers/bramble_pasciak.cpp

@@ -32,6 +32,7 @@ BramblePasciakSolver::BramblePasciakSolver(ParBilinearForm &mVarf,
    std::unique_ptr<HypreParMatrix> invDBt(B_->Transpose());
    invDBt->InvScaleRows(diagM);
    S_.reset(ParMult(B_.get(), invDBt.get(), true));
+   invDBt.reset();
    M0_.Reset(new HypreDiagScale(*M_));
    M1_.Reset(new HypreBoomerAMG(*S_));
    M1_.As<HypreBoomerAMG>()->SetPrintLevel(0);

miniapps/solvers/darcy_solver.cpp

@@ -57,6 +57,7 @@ BDPMinresSolver::BDPMinresSolver(const HypreParMatrix& M,
 
 void BDPMinresSolver::Mult(const Vector & x, Vector & y) const
 {
+   solver_.iterative_mode = this->iterative_mode;
    solver_.Mult(x, y);
    for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }
 }

miniapps/solvers/darcy_solver.hpp

@@ -52,7 +52,7 @@ class BDPMinresSolver : public DarcySolver
    BlockDiagonalPreconditioner prec_;
    OperatorPtr BT_;
    OperatorPtr S_;   // S_ = B diag(M)^{-1} B^T
-   MINRESSolver solver_;
+   mutable MINRESSolver solver_;
    Array<int> ess_zero_dofs_;
 public:
    BDPMinresSolver(const HypreParMatrix& M,

miniapps/solvers/div_free_solver.cpp

@@ -84,7 +84,6 @@ DFSSpaces::DFSSpaces(int order, int num_refine, ParMesh *mesh,
    data_.Q_l2.resize(num_refine);
    hdiv_fes_->GetEssentialTrueDofs(ess_attr, data_.coarsest_ess_hdivdofs);
    data_.C.resize(num_refine+1);
-   data_.Ae.resize(num_refine+1);
 
    hcurl_fes_ = std::make_unique<ParFiniteElementSpace>(mesh, hcurl_fec_.get());
    coarse_hcurl_fes_ = std::make_unique<ParFiniteElementSpace>(*hcurl_fes_);
@@ -174,9 +173,9 @@ void DFSSpaces::CollectDFSData()
    data_.C[level_+1].Reset(curl.ParallelAssemble());
    mfem::Array<int> ess_hcurl_tdof;
    hcurl_fes_->GetEssentialTrueDofs(ess_bdr_attr_, ess_hcurl_tdof);
-   data_.Ae[level_+1].reset(
-      data_.C[level_+1].As<HypreParMatrix>()
-      ->EliminateCols(ess_hcurl_tdof));
+   HypreParMatrix *res =
+      data_.C[level_+1].As<HypreParMatrix>()->EliminateCols(ess_hcurl_tdof);
+   delete res;
 
    ++level_;
 

miniapps/solvers/div_free_solver.hpp

@@ -36,7 +36,6 @@ struct DFSParameters : IterSolveParameters
 struct DFSData
 {
    using UniqueOperatorPtr = std::unique_ptr<OperatorPtr>;
-   using UniqueHypreParMatrix = std::unique_ptr<HypreParMatrix>;
 
    std::vector<OperatorPtr> agg_hdivdof; // agglomerates to H(div) dofs table
    std::vector<OperatorPtr> agg_l2dof; // agglomerates to L2 dofs table
@@ -46,7 +45,6 @@ struct DFSData
    std::vector<OperatorPtr> Q_l2; // Q_l2[l] = (W_{l+1})^{-1} P_l2[l]^T W_l
    Array<int> coarsest_ess_hdivdofs; // coarsest level essential H(div) dofs
    std::vector<OperatorPtr> C; // discrete curl: ND -> RT, map to Null(B)
-   std::vector<UniqueHypreParMatrix> Ae;
    DFSParameters param;
 };
 

miniapps/solvers/plor_solvers.cpp

@@ -1,69 +1,3 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-//
-//                 ------------------------------------------
-//                 Parallel Low-Order Refined Solvers Miniapp
-//                 ------------------------------------------
-//
-// This miniapp illustrates the use of low-order refined preconditioners for
-// finite element problems defined using H1, H(curl), H(div), or L2 finite
-// element spaces. The following problems are solved, depending on the chosen
-// finite element space:
-//
-// H1 and L2: definite Helmholtz problem, u - Delta u = f
-//            (in L2 discretized using the symmetric interior penalty DG method)
-//
-// H(curl):   definite Maxwell problem, u + curl curl u = f
-//
-// H(div):    grad-div problem, u - grad(div u) = f
-//
-// In each case, the high-order finite element problem is preconditioned using a
-// low-order finite element discretization defined on a Gauss-Lobatto refined
-// mesh. The low-order problem is solved using hypre's AMG preconditioners:
-// BoomerAMG is used for H1 and L2 problems, AMS is used for H(curl) and 2D
-// H(div) problems, and ADS is used for 3D H(div) problems.
-//
-// For vector finite element spaces, the special "Integrated" basis type is used
-// to obtain spectral equivalence between the high-order and low-order refined
-// discretizations. This basis is defined in reference [1] and spectral
-// equivalence is shown in [2]:
-//
-// [1]. M. Gerritsma. Edge functions for spectral element methods. Spectral and
-//      High Order Methods for Partial Differential Equations. (2010)
-// [2]. C. Dohrmann. Spectral equivalence properties of higher-order tensor
-//      product finite elements and applications to preconditioning. (2021)
-//
-// The action of the high-order operator is computed using MFEM's partial
-// assembly/matrix-free algorithms (except in the case of L2, which remains
-// future work).
-//
-// Compile with: make plor_solvers
-//
-// Sample runs:
-//
-//    mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe h
-//    mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe n
-//    mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe r
-//    mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe l
-//    mpirun -np 4 plor_solvers -m ../../data/amr-hex.mesh -fe h -rs 0 -o 2
-//    mpirun -np 4 plor_solvers -m ../../data/star-surf.mesh -fe h
-//    mpirun -np 4 plor_solvers -m ../../data/star-surf.mesh -fe n
-//    mpirun -np 4 plor_solvers -m ../../data/star-surf.mesh -fe r
-//
-// Device sample runs:
-//  * mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe h -d cuda
-//  * mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe n -d cuda
-//  * mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe r -d cuda
-//  * mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe l -d cuda
-
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
@@ -74,54 +8,35 @@
 using namespace std;
 using namespace mfem;
 
-int main(int argc, char *argv[])
+struct Opts
 {
-   Mpi::Init();
-   Hypre::Init();
-
    const char *mesh_file = "../../data/star.mesh";
-   int ser_ref_levels = 1, par_ref_levels = 1;
+   int ser_ref_levels = 1;
+   int par_ref_levels = 1;
    int order = 3;
    const char *fe = "h";
-   const char *device_config = "cpu";
-   bool visualization = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
-   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(&order, "-o", "--order", "Polynomial degree.");
-   args.AddOption(&fe, "-fe", "--fe-type",
-                  "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.ParseCheck();
-
-   Device device(device_config);
-   if (Mpi::Root()) { device.Print(); }
+};
 
+int Run(const Opts &opts)
+{
    bool H1 = false, ND = false, RT = false, L2 = false;
-   if (string(fe) == "h") { H1 = true; }
-   else if (string(fe) == "n") { ND = true; }
-   else if (string(fe) == "r") { RT = true; }
-   else if (string(fe) == "l") { L2 = true; }
+   if (string(opts.fe) == "h") { H1 = true; }
+   else if (string(opts.fe) == "n") { ND = true; }
+   else if (string(opts.fe) == "r") { RT = true; }
+   else if (string(opts.fe) == "l") { L2 = true; }
    else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
 
-   real_t kappa = (order+1)*(order+1); // Penalty used for DG discretizations
+   const int order = opts.order;
+   const real_t kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
-   Mesh serial_mesh(mesh_file, 1, 1);
+   Mesh serial_mesh(opts.mesh_file, 1, 1);
    const int dim = serial_mesh.Dimension();
    const int sdim = serial_mesh.SpaceDimension();
    MFEM_VERIFY(dim == 2 || dim == 3, "Mesh dimension must be 2 or 3.");
    MFEM_VERIFY(!L2 || dim == sdim, "DG surface meshes not supported.");
-   for (int l = 0; l < ser_ref_levels; l++) { serial_mesh.UniformRefinement(); }
+   for (int l = 0; l < opts.ser_ref_levels; l++) { serial_mesh.UniformRefinement(); }
    ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
-   for (int l = 0; l < par_ref_levels; l++) { mesh.UniformRefinement(); }
+   for (int l = 0; l < opts.par_ref_levels; l++) { mesh.UniformRefinement(); }
    serial_mesh.Clear();
 
    if (mesh.ncmesh && (RT || ND))
@@ -139,6 +54,18 @@ int main(int argc, char *argv[])
    else { fec.reset(new L2_FECollection(order, dim, b1)); }
 
    ParFiniteElementSpace fes(&mesh, fec.get());
+
+   // fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+
+   {
+      MFEM_PERF_SCOPE("Ensure Nodes");
+      mesh.EnsureNodes();
+   }
+   // {
+   //    auto &ir = DiffusionIntegrator::GetRule(*fes.GetFE(0), *fes.GetFE(0));
+   //    mesh.GetGeometricFactors(ir, GeometricFactors::JACOBIANS);
+   // }
+
    HYPRE_Int ndofs = fes.GlobalTrueVSize();
    if (Mpi::Root()) { cout << "Number of DOFs: " << ndofs << endl; }
 
@@ -149,7 +76,7 @@ int main(int argc, char *argv[])
    ParBilinearForm a(&fes);
    if (H1 || L2)
    {
-      a.AddDomainIntegrator(new MassIntegrator);
+      // a.AddDomainIntegrator(new MassIntegrator);
       a.AddDomainIntegrator(new DiffusionIntegrator);
    }
    else
@@ -167,7 +94,9 @@ int main(int argc, char *argv[])
    // Partial assembly not currently supported for DG or for surface meshes with
    // vector finite elements (ND or RT).
    if (!L2 && (H1 || sdim == dim)) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+
    a.Assemble();
+   // a.Assemble();
 
    ParLinearForm b(&fes);
    if (H1 || L2) { b.AddDomainIntegrator(new DomainLFIntegrator(f_coeff)); }
@@ -188,9 +117,13 @@ int main(int argc, char *argv[])
    a.FormLinearSystem(ess_dofs, x, b, A, X, B);
 
    unique_ptr<Solver> solv_lor;
+
    if (H1 || L2)
    {
-      solv_lor.reset(new LORSolver<HypreBoomerAMG>(a, ess_dofs));
+      auto solv = new LORSolver<HypreBoomerAMG>(a, ess_dofs);
+      solv->GetSolver().SetPrintLevel(0);
+      solv->GetSolver().Setup(B, X);
+      solv_lor.reset(solv);
    }
    else if (RT && dim == 3)
    {
@@ -212,30 +145,53 @@ int main(int argc, char *argv[])
 
    a.RecoverFEMSolution(X, b, x);
 
-   if (sdim == dim)
-   {
-      real_t er =
-         (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
-      if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
-   }
-
-   if (visualization)
-   {
-      // Save the solution and mesh to disk. The output can be viewed using
-      // GLVis as follows: "glvis -np <np> -m mesh -g sol"
-      x.Save("sol");
-      mesh.Save("mesh");
-
-      // Also save the solution for visualization using ParaView
-      ParaViewDataCollection dc("PLOR", &mesh);
-      dc.SetPrefixPath("ParaView");
-      dc.SetHighOrderOutput(true);
-      dc.SetLevelsOfDetail(order);
-      dc.RegisterField("u", &x);
-      dc.SetCycle(0);
-      dc.SetTime(0.0);
-      dc.Save();
-   }
-
    return 0;
 }
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+   Hypre::Init();
+
+   Opts opts;
+   const char *device_config = "cpu";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&opts.mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&opts.ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&opts.par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&opts.order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&opts.fe, "-fe", "--fe-type",
+                  "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.ParseCheck();
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   MFEM_PERF_SYNC(true);
+
+   {
+      MFEM_PERF_SCOPE("Temporary allocations");
+      Vector tmp1(1024 * 1024 * 1024);
+      tmp1.ReadWrite();
+      Vector tmp2(1024 * 1024 * 1024);
+      tmp2.ReadWrite();
+   }
+   {
+      MFEM_PERF_SCOPE("Hypre allocations");
+      double *tmp1 = mfem_hypre_CTAlloc(double, 1024 * 1024 * 1024);
+      double *tmp2 = mfem_hypre_CTAlloc(double, 1024 * 1024 * 1024);
+      mfem_hypre_TFree(tmp2);
+      mfem_hypre_TFree(tmp1);
+   }
+
+   MFEM_PERF_DISABLE;
+   Run(opts);
+
+   MFEM_PERF_ENABLE;
+   Run(opts);
+}

tests/benchmarks/CMakeLists.txt

@@ -31,11 +31,6 @@ function(add_benchmark name)
        set_property(SOURCE ${${NAME}_BENCH_SRCS} PROPERTY LANGUAGE CUDA)
     endif(MFEM_USE_CUDA)
 
-    if (MFEM_USE_HIP)
-        set_property(SOURCE ${${NAME}_BENCH_SRCS} PROPERTY LANGUAGE 
-                     HIP_SOURCE_PROPERTY_FORMAT TRUE)
-    endif(MFEM_USE_HIP)
-
     add_executable(bench_${name} ${${NAME}_BENCH_SRCS})
     target_link_libraries(bench_${name} mfem pthread)
     add_dependencies(${MFEM_ALL_BENCHMARKS_TARGET_NAME} bench_${name})
@@ -61,3 +56,8 @@ add_benchmark(elasticity)
 add_benchmark(tmop)
 add_benchmark(vector)
 add_benchmark(virtuals)
+
+if (MFEM_USE_MPI)
+    mfem_add_executable(pbench_ceed pbench_ceed.cpp)
+    target_link_libraries(pbench_ceed mfem pthread)
+endif()

tests/benchmarks/bench.hpp

@@ -57,6 +57,14 @@ public:
    }
 };
 
+// Reporter with no output, used with MPI on non-root ranks.
+struct NoReporter : public benchmark::BenchmarkReporter
+{
+   explicit NoReporter() = default;
+   bool ReportContext(const Context &) override { return true; }
+   void ReportRuns(const std::vector<Run> &) override {}
+};
+
 } // namespace mfem
 
 #endif // MFEM_USE_BENCHMARK

tests/benchmarks/makefile

@@ -22,7 +22,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_TESTS = bench_assembly_levels bench_ceed bench_dg_amr bench_elasticity \
             bench_tmop bench_vector bench_virtuals
-PAR_TESTS = 
+PAR_TESTS = pbench_ceed
 ifeq ($(MFEM_USE_MPI),NO)
    TESTS = $(SEQ_TESTS)
 else
@@ -43,7 +43,7 @@ all: $(TESTS)
 # Rules for building the TESTS
 
 %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+	$(MFEM_CXX) $(MFEM_FLAGS) -I$(MFEM_DIR) $< -o $@ $(MFEM_LIBS)
 
 # Rules for compiling miniapp dependencies
 $($(TESTS)): \

tests/benchmarks/pbench_ceed.cpp

@@ -0,0 +1,606 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// Compile with: make pbench_ceed
+//
+// Sample runs:
+//    mpirun -np 4 pbench_ceed
+//    mpirun -np 4 pbench_ceed --benchmark_filter=BP3
+//    mpirun -np 4 pbench_ceed --benchmark_filter=BP3 --benchmark_context=local_size=5e4
+//    mpirun -np 6 pbench_ceed --benchmark_filter=BP3 --benchmark_context=proc_grid=3x2x1,local_size=5e4
+//
+// Device sample runs:
+//    mpirun -np 4 pbench_ceed --benchmark_context=device=cuda,local_size=1e6
+//    mpirun -np 4 pbench_ceed --benchmark_filter=BP3 --benchmark_context=device=cuda,local_size=1e7
+//
+// Description:
+//    This benchmark contains the implementation of the CEED's bake-off
+//    problems, BP1-BP6, and bake-off kernels, BK1-BK6: high-order benchmarks
+//    designed to test and compare the performance of high-order codes.
+//
+//    See: ceed.exascaleproject.org/bps and github.com/CEED/benchmarks
+
+
+#include "bench.hpp"
+
+#ifdef MFEM_USE_BENCHMARK
+
+// for instantiating more kernels:
+#include <fem/integ/bilininteg_mass_kernels.hpp>
+#include <fem/integ/bilininteg_diffusion_kernels.hpp>
+
+// Global parameters
+
+// local_size: desired approximate MPI-local problem size; this local size and
+// the polynomial order determine the local mesh size so that the resulting
+// problem size is (approximately) equal to local_size for all polynomial
+// orders, see MakeParMesh().
+double local_size = 61*61*61; // exact size match for all p=1,...,6
+
+// proc_grid: use processor grid given by proc_grid[0,1,2].
+int proc_grid[3] = {0,0,0};
+
+// q_gl_inc: increment for the number of GL points: q = p + 1 + q_gl_inc
+int q_gl_inc = 0;
+
+// q_gll_inc: increment for the number of GLL points: q = p + 1 + q_gll_inc
+int q_gll_inc = 0;
+
+// verbose: verbosity level: 0, 1, 2
+int verbose = 0;
+
+
+// If running on GPU, wait for GPU tasks to finish:
+inline void DeviceSync()
+{
+   if (Device::Allows(Backend::DEVICE_MASK & ~Backend::DEBUG_DEVICE))
+   {
+      MFEM_STREAM_SYNC;
+      // MFEM_DEVICE_SYNC;
+   }
+}
+
+void MakeExp2ProcGrid(int np)
+{
+   proc_grid[0] = proc_grid[1] = proc_grid[2] = 1;
+   for ( ; np >= 8; np /= 8)
+   {
+      proc_grid[0] *= 2; proc_grid[1] *= 2; proc_grid[2] *= 2;
+   }
+   if (np == 4) { proc_grid[0] *= 2; proc_grid[1] *= 2; }
+   else if (np == 2) { proc_grid[0] *= 2; }
+}
+
+// Construct the parallel mesh based on the polynomial order, p, and the
+// local_size:
+ParMesh MakeParMesh(int p)
+{
+   int nx = 0, ny = 0, nz = 0;
+   int par_ref = 0;
+   if (verbose && Mpi::Root()) { std::cout << _MFEM_FUNC_NAME << std::endl; }
+
+   const double s = local_size;
+
+   int m = floor((pow(s, 1./3)-1)/p);
+   double s_l, s_u, s_c;
+   while ((s_l=(  m  *p+1)*(  m  *p+1)*(  m  *p+1), s_l > s))  { m--; }
+   m = std::max(m, 1);
+   while ((s_u=((m+1)*p+1)*((m+1)*p+1)*((m+1)*p+1), s_u <= s)) { m++; }
+   s_l = (m*p+1)*(m*p+1)*(m*p+1);
+   if ((s_c=((m+1)*p+1)*(m*p+1)*(m*p+1), s_c > s))
+   {
+      if (s/s_l <= s_c/s) { nx = m; ny = m; nz = m; }
+      else                { nx = m; ny = m; nz = m + 1; }
+   }
+   else if ((s_l=s_c, s_c=((m+1)*p+1)*((m+1)*p+1)*(m*p+1), s_c > s))
+   {
+      if (s/s_l <= s_c/s) { nx = m; ny = m;     nz = m + 1; }
+      else                { nx = m; ny = m + 1; nz = m + 1; }
+   }
+   else
+   {
+      s_l=s_c, s_c=s_u;
+      if (s/s_l <= s_c/s) { nx = m;     ny = m + 1; nz = m + 1; }
+      else                { nx = m + 1; ny = m + 1; nz = m + 1; }
+   }
+   while (nx%2 == 0 && ny%2 == 0 && nz%2 == 0)
+   {
+      par_ref++;
+      nx /= 2; ny /= 2; nz /= 2;
+   }
+
+   nx *= proc_grid[0];
+   ny *= proc_grid[1];
+   nz *= proc_grid[2];
+
+   if (verbose && Mpi::Root())
+   {
+      std::cout
+            << '\n'
+            << "   order: " << p << '\n'
+            << "   nx: " << nx << ", ny: " << ny << ", nz: " << nz << '\n'
+            << "   px: " << proc_grid[0] << ", py: " << proc_grid[1]
+            << ", pz: " << proc_grid[2] << '\n'
+            << "   par_ref: " << par_ref << '\n'
+            << std::endl;
+   }
+
+   StopWatch timer;
+   timer.Start();
+   double t_start = timer.RealTime();
+
+   Mesh smesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON));
+   if (verbose && Mpi::Root())
+   {
+      double t_elapsed = timer.RealTime() - t_start;
+      std::cout << "   Mesh: " << 1e3*t_elapsed << " ms" << std::endl;
+   }
+   t_start = timer.RealTime();
+   Array<int> partitioning;
+   partitioning.MakeRef(smesh.CartesianPartitioning(proc_grid), smesh.GetNE(),
+                        MemoryType::HOST, true);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh, partitioning.HostRead());
+   smesh.Clear();
+   for (int i = 0; i < par_ref; i++)
+   {
+      pmesh.UniformRefinement();
+   }
+   if (verbose && Mpi::Root())
+   {
+      double t_elapsed = timer.RealTime() - t_start;
+      std::cout << "   ParMesh: " << 1e3*t_elapsed << " ms" << std::endl;
+   }
+   return pmesh;
+}
+
+template <int VDIM, bool GLL>
+struct BakeOff
+{
+   static constexpr int DIM = 3;
+   const int p, q, q_order;
+   ParMesh mesh;
+   H1_FECollection fec;
+   ParFiniteElementSpace fes;
+   const Geometry::Type geom_type;
+   IntegrationRules irs;
+   const IntegrationRule *ir;
+   ConstantCoefficient one;
+   Vector uvec;
+   VectorConstantCoefficient unit_vec;
+   const HYPRE_BigInt dofs;
+   double mdofs{};
+   StopWatch timer;
+
+   BakeOff(int p):
+      p(p),
+      q(GLL ? p + 1 + q_gll_inc : p + 1 + q_gl_inc),
+      q_order(2 * q + (GLL ? -3 : -1)),
+      mesh(MakeParMesh(p)),
+      fec(p, DIM, BasisType::GaussLobatto),
+      fes(&mesh, &fec, VDIM, VDIM == 3 ? Ordering::byVDIM : Ordering::byNODES),
+      geom_type(mesh.GetTypicalElementGeometry()),
+      irs(0, GLL ? Quadrature1D::GaussLobatto : Quadrature1D::GaussLegendre),
+      ir(&irs.Get(geom_type, q_order)),
+      one(1.0),
+      uvec(DIM),
+      unit_vec((uvec = 1.0, uvec /= uvec.Norml2(), uvec)),
+      dofs(fes.GlobalTrueVSize())
+   {
+      if (verbose && Mpi::Root())
+      {
+         std::cout << "q: " << q << ", dofs: " << dofs << std::endl;
+         // std::cout << _MFEM_FUNC_NAME << std::endl;
+      }
+      timer.Start();
+   }
+
+   virtual void benchmark(benchmark::State &state) = 0;
+
+   double SumMdofs() const { return mdofs; }
+
+   double MDofs() const { return 1e-6 * dofs; }
+};
+
+/// Bake-off Problems (BPs)
+template <typename BFI, int VDIM, bool GLL>
+struct Problem : public BakeOff<VDIM, GLL>
+{
+   const double rtol = 1e-16;
+   const int max_it = 20;
+   const int print_lvl = -1;
+
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   ParGridFunction x;
+   ParBilinearForm a;
+   LinearForm b;
+   OperatorPtr A;
+   Vector B, X;
+   CGSolver cg;
+   int bench_call_counter = 0;
+
+   using base = BakeOff<VDIM, GLL>;
+   using base::ir;
+   using base::one;
+   using base::mesh;
+   using base::fes;
+   using base::mdofs;
+   using base::timer;
+
+   Problem(int order):
+      BakeOff<VDIM, GLL>(order),
+      ess_bdr(mesh.bdr_attributes.Max()),
+      x(&fes),
+      a(&fes),
+      b(&fes),
+      cg(fes.GetComm())
+   {
+      if (verbose && Mpi::Root()) { std::cout << _MFEM_FUNC_NAME << std::endl; }
+      ess_bdr = 1;
+      x = 0.0;
+      fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      LinearFormIntegrator *integ;
+      if (VDIM == 1)
+      {
+         integ = new DomainLFIntegrator(this->one);
+      }
+      else
+      {
+         integ = new VectorDomainLFIntegrator(this->unit_vec);
+      }
+      integ->SetIntRule(ir);
+      b.AddDomainIntegrator(integ); // b takes ownership of integ
+      b.UseFastAssembly(true);
+      b.Assemble();
+
+      double t_start = timer.RealTime();
+      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      a.AddDomainIntegrator(new BFI(one, ir));
+      a.Assemble();
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+      if (verbose && Mpi::Root())
+      {
+         double t_elapsed = timer.RealTime() - t_start;
+         std::cout << "   assemble a: " << 1e3*t_elapsed << " ms" << std::endl;
+      }
+
+      cg.SetRelTol(rtol);
+      cg.SetOperator(*A);
+      cg.SetMaxIter(max_it);
+      cg.SetPrintLevel(print_lvl);
+      cg.iterative_mode = false;
+      // warmup cg
+      cg.SetMaxIter(2);
+      cg.Mult(B, X);
+      cg.SetMaxIter(max_it);
+   }
+
+   virtual ~Problem()
+   {
+      if (verbose && Mpi::Root())
+      {
+         std::cout << _MFEM_FUNC_NAME << '\n'
+                   << "   call counter: " << bench_call_counter << '\n'
+                   << "   MDofs: " << mdofs << std::endl;
+      }
+   }
+
+   void benchmark(benchmark::State &state) override
+   {
+      if (verbose > 1 && Mpi::Root())
+      {
+         std::cout << _MFEM_FUNC_NAME << std::endl;
+      }
+
+      DeviceSync();
+      MPI_Barrier(cg.GetComm());
+      double t_start = timer.RealTime();
+
+      cg.Mult(B, X);
+
+      DeviceSync();
+      MPI_Barrier(cg.GetComm());
+      double t_elapsed = timer.RealTime() - t_start;
+      // Ensure every rank gets the same time, otherwise google-benchmark may
+      // behave differently on different ranks.
+      MPI_Bcast(&t_elapsed, 1, MPI_DOUBLE, 0, cg.GetComm());
+
+      state.SetIterationTime(t_elapsed);
+      if (verbose > 1 && Mpi::Root())
+      {
+         std::cout << "   bench time: " << 1e3*t_elapsed << " ms" << std::endl;
+      }
+
+      mdofs += this->MDofs() * cg.GetNumIterations();
+      bench_call_counter++;
+   }
+};
+
+/// Bake-off Problems (BPs)
+#define BakeOff_Problem(i, Kernel, VDIM, GLL)                        \
+   static void BP##i(bm::State &state)                               \
+   {                                                                 \
+      Problem<Kernel##Integrator, VDIM, GLL> ker(state.range(0));    \
+      for (auto z : state) { ker.benchmark(state); }                 \
+      state.counters["Num Dofs"] = ker.dofs;                         \
+      state.counters["|      Dof/s"] =                               \
+         bm::Counter(1e6*ker.SumMdofs(), bm::Counter::kIsRate);      \
+      state.counters["|   Dof/s/NP"] =                               \
+         bm::Counter(1e6*ker.SumMdofs()/ker.fes.GetNRanks(),         \
+                     bm::Counter::kIsRate);                          \
+   }                                                                 \
+   BENCHMARK(BP##i)->DenseRange(1, 6)->Unit(bm::kMillisecond)->UseManualTime();
+
+// state.counters[" Q1D"] = ker.q;
+
+
+/// BP1: scalar PCG with mass matrix, GL
+BakeOff_Problem(1, Mass, 1, false)
+
+/// BP2: vector PCG with mass matrix, GL
+BakeOff_Problem(2, VectorMass, 3, false)
+
+/// BP3: scalar PCG with stiffness matrix, GL
+BakeOff_Problem(3, Diffusion, 1, false)
+
+/// BP4: vector PCG with stiffness matrix, GL
+BakeOff_Problem(4, VectorDiffusion, 3, false)
+
+/// BP5: scalar PCG with stiffness matrix, GLL
+BakeOff_Problem(5, Diffusion, 1, true)
+
+/// BP6: vector PCG with stiffness matrix, GLL
+BakeOff_Problem(6, VectorDiffusion, 3, true)
+
+/// Bake-off Kernels (BKs)
+template <typename BFI, int VDIM, bool GLL>
+struct Kernel : public BakeOff<VDIM, GLL>
+{
+   BFI bfi;
+   Vector x, y; // input and output E-vectors
+
+   using base = BakeOff<VDIM, GLL>;
+   using base::ir;
+   using base::one;
+   using base::fes;
+   using base::mdofs;
+   using base::timer;
+
+   Kernel(int order)
+      : base(order),
+        bfi(one, ir)
+   {
+      bfi.AssemblePA(fes);
+      const Table &el2dof = fes.GetElementToDofTable();
+      const int e_size = el2dof.Size_of_connections()*fes.GetVDim();
+      x.SetSize(e_size); x.UseDevice(true);
+      y.SetSize(e_size); y.UseDevice(true);
+      x.Randomize(1);
+      x.Read();
+      // warmup
+      for (int i = 0; i < 2; i++)
+      {
+         y = 0.0;
+         bfi.AddMultPA(x, y);
+      }
+   }
+
+   void benchmark(benchmark::State &state) override
+   {
+      y = 0.0;
+
+      DeviceSync();
+      MPI_Barrier(fes.GetComm());
+      double t_start = timer.RealTime();
+
+      bfi.AddMultPA(x, y);
+
+      DeviceSync();
+      MPI_Barrier(fes.GetComm());
+      double t_elapsed = timer.RealTime() - t_start;
+      // Ensure every rank gets the same time, otherwise google-benchmark may
+      // behave differently on different ranks.
+      MPI_Bcast(&t_elapsed, 1, MPI_DOUBLE, 0, fes.GetComm());
+
+      state.SetIterationTime(t_elapsed);
+      mdofs += this->MDofs();
+   }
+};
+
+/// Generic CEED BKi
+#define BakeOff_Kernel(i, KER, VDIM, GLL)                       \
+   static void BK##i(bm::State &state)                          \
+   {                                                            \
+      Kernel<KER##Integrator, VDIM, GLL> ker(state.range(0));   \
+      for (auto z : state) { ker.benchmark(state); }            \
+      state.counters["Num Dofs"] = ker.dofs;                    \
+      state.counters["|      Dof/s"] =                          \
+         bm::Counter(1e6*ker.SumMdofs(), bm::Counter::kIsRate); \
+      state.counters["|   Dof/s/NP"] =                          \
+         bm::Counter(1e6*ker.SumMdofs()/ker.fes.GetNRanks(),    \
+                     bm::Counter::kIsRate);                     \
+   }                                                            \
+   BENCHMARK(BK##i)->DenseRange(1, 6)->Unit(bm::kMillisecond)->UseManualTime();
+
+// state.counters[" Q1D"] = ker.q;
+
+
+/// BK1: scalar E-vector-to-E-vector evaluation of mass matrix, GL
+BakeOff_Kernel(1, Mass, 1, false)
+
+/// BK2: vector E-vector-to-E-vector evaluation of mass matrix, GL
+BakeOff_Kernel(2, VectorMass, 3, false)
+
+/// BK3: scalar E-vector-to-E-vector evaluation of stiffness matrix, GL
+BakeOff_Kernel(3, Diffusion, 1, false)
+
+/// BK4: vector E-vector-to-E-vector evaluation of stiffness matrix, GL
+BakeOff_Kernel(4, VectorDiffusion, 3, false)
+
+/// BK5: scalar E-vector-to-E-vector evaluation of stiffness matrix, GLL
+BakeOff_Kernel(5, Diffusion, 1, true)
+
+/// BK6: vector E-vector-to-E-vector evaluation of stiffness matrix, GLL
+BakeOff_Kernel(6, VectorDiffusion, 3, true)
+
+
+int main(int argc, char *argv[])
+{
+   // MassIntegrator specializations by <DIM, D1D, Q1D>
+   MassIntegrator::AddSpecialization<3, 3, 3>();
+   MassIntegrator::AddSpecialization<3, 4, 4>();
+   MassIntegrator::AddSpecialization<3, 5, 5>();
+   MassIntegrator::AddSpecialization<3, 6, 6>();
+   MassIntegrator::AddSpecialization<3, 7, 7>();
+
+   // DiffusionIntegrator specializations by <DIM, D1D, Q1D>
+   DiffusionIntegrator::AddSpecialization<3, 3, 3>();
+   DiffusionIntegrator::AddSpecialization<3, 4, 4>();
+   DiffusionIntegrator::AddSpecialization<3, 5, 5>();
+   DiffusionIntegrator::AddSpecialization<3, 6, 6>();
+   DiffusionIntegrator::AddSpecialization<3, 7, 7>();
+
+   Mpi::Init();
+   Hypre::Init();
+
+   bm::Initialize(&argc, argv);
+
+   // Device setup, cpu by default
+   std::string device_config = "cpu";
+   bool gpu_aware_mpi = false;
+   std::string proc_grid_str = "";
+
+   auto global_context = bmi::GetGlobalContext();
+   if (global_context != nullptr)
+   {
+      const auto device = global_context->find("device");
+      if (device != global_context->end())
+      {
+         device_config = device->second;
+      }
+      const auto ctx_gpu_aware_mpi = global_context->find("gpu_aware_mpi");
+      if (ctx_gpu_aware_mpi != global_context->end())
+      {
+         gpu_aware_mpi = std::atoi(ctx_gpu_aware_mpi->second.c_str());
+      }
+      const auto ctx_local_size = global_context->find("local_size");
+      if (ctx_local_size != global_context->end())
+      {
+         std::size_t pos;
+         local_size = std::stof(ctx_local_size->second, &pos);
+         if (ctx_local_size->second.size() != pos)
+         {
+            if (Mpi::Root())
+            {
+               std::cout << "\nerror reading local_size: "
+                         << ctx_local_size->second << '\n' << std::endl;
+            }
+            return 1;
+         }
+         if (local_size < 64.0 || local_size > std::exp2(30.0))
+         {
+            if (Mpi::Root())
+            {
+               std::cout << "\nlocal_size must be in [2^6,2^30]! local_size: "
+                         << local_size << '\n' << std::endl;
+            }
+            return 1;
+         }
+      }
+      const auto ctx_proc_grid = global_context->find("proc_grid");
+      if (ctx_proc_grid != global_context->end())
+      {
+         proc_grid_str = ctx_proc_grid->second;
+      }
+      const auto ctx_verbose = global_context->find("verbose");
+      if (ctx_verbose != global_context->end())
+      {
+         verbose = std::atoi(ctx_verbose->second.c_str());
+      }
+   }
+   const int num_procs = Mpi::WorldSize();
+   if (proc_grid_str == "" || proc_grid_str == "2^n")
+   {
+      if (((num_procs-1)&num_procs) != 0)
+      {
+         if (Mpi::Root())
+         {
+            std::cout << "\nthe number of processors is not a power of 2!"
+                      << " num_procs: " << num_procs
+                      << "\nuse a processor grid, e.g. "
+                      << "--benchmark_context=proc_grid=3x5x7\n"
+                      << std::endl;
+         }
+         return 1;
+      }
+      MakeExp2ProcGrid(num_procs);
+   }
+   else
+   {
+      int n = std::sscanf(proc_grid_str.c_str(), "%d x %d x %d",
+                          &proc_grid[0], &proc_grid[1], &proc_grid[2]);
+      if (n != 3)
+      {
+         if (Mpi::Root())
+         {
+            std::cout << "\ninvalid processor grid input: "
+                      << proc_grid_str << "\n" << std::endl;
+         }
+         return 1;
+      }
+      if (proc_grid[0]*proc_grid[1]*proc_grid[2] != num_procs ||
+          proc_grid[0] < 1 || proc_grid[1] < 1 || proc_grid[2] < 1)
+      {
+         if (Mpi::Root())
+         {
+            std::cout << "\ninvalid processor grid: " << proc_grid[0] << " x "
+                      << proc_grid[1] << " x " << proc_grid[2] << " != "
+                      << num_procs << '\n' << std::endl;
+         }
+         return 1;
+      }
+   }
+
+   if (bm::ReportUnrecognizedArguments(argc, argv)) { return 1; }
+
+   Device device(device_config.c_str());
+   device.SetGPUAwareMPI(gpu_aware_mpi);
+   if (Mpi::Root())
+   {
+      device.Print();
+      std::cout << "    num procs: " << num_procs << std::endl;
+      std::cout << "gpu aware mpi: " << gpu_aware_mpi << std::endl;
+      std::cout << "   local size: " << local_size << std::endl;
+      std::cout << "    proc grid: " << proc_grid[0] << 'x'
+                /**/                 << proc_grid[1] << 'x'
+                /**/                 << proc_grid[2] << std::endl;
+      std::cout << "     GL q_inc: " << q_gl_inc << std::endl;
+      std::cout << "    GLL q_inc: " << q_gll_inc << std::endl;
+   }
+
+   DeviceSync();
+   MPI_Barrier(MPI_COMM_WORLD);
+
+   if (Mpi::Root())
+   {
+      bm::ConsoleReporter CR;
+      bm::RunSpecifiedBenchmarks(&CR);
+   }
+   else
+   {
+      NoReporter NR;
+      bm::RunSpecifiedBenchmarks(&NR);
+   }
+
+   return 0;
+}
+
+#endif // MFEM_USE_BENCHMARK

tests/unit/linalg/test_chebyshev.cpp

@@ -14,53 +14,47 @@
 
 using namespace mfem;
 
-TEST_CASE("OperatorChebyshevSmoother", "[Chebyshev symmetry]")
+TEST_CASE("Chebyshev symmetry", "[OperatorChebyshevSmoother]")
 {
-   for (int order = 2; order < 5; ++order)
-   {
-      const int cheb_order = 2;
+   const int order = GENERATE(2, 3, 4);
+   const int cheb_order = GENERATE(2, 3);
 
-      Mesh mesh = Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
-      FiniteElementCollection *fec = new H1_FECollection(order, 3);
-      FiniteElementSpace fespace(&mesh, fec);
-      Array<int> ess_bdr(mesh.bdr_attributes.Max());
-      ess_bdr = 1;
-      Array<int> ess_tdof_list;
-      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   Mesh mesh = Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
+   H1_FECollection fec(order, 3);
+   FiniteElementSpace fespace(&mesh, &fec);
 
-      BilinearForm aform(&fespace);
-      aform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-      aform.AddDomainIntegrator(new DiffusionIntegrator);
-      aform.Assemble();
-      OperatorPtr opr;
-      opr.SetType(Operator::ANY_TYPE);
-      aform.FormSystemMatrix(ess_tdof_list, opr);
-      Vector diag(fespace.GetTrueVSize());
-      aform.AssembleDiagonal(diag);
+   Array<int> ess_tdof_list;
+   fespace.GetBoundaryTrueDofs(ess_tdof_list);
 
-      Solver* smoother = new OperatorChebyshevSmoother(*opr, diag, ess_tdof_list,
-                                                       cheb_order);
+   BilinearForm aform(&fespace);
+   aform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   aform.AddDomainIntegrator(new DiffusionIntegrator);
+   aform.Assemble();
 
-      int n = smoother->Width();
-      Vector left(n);
-      Vector right(n);
-      int seed = (int) time(0);
-      left.Randomize(seed);
-      right.Randomize(seed + 2);
+   OperatorPtr opr;
+   opr.SetType(Operator::ANY_TYPE);
+   aform.FormSystemMatrix(ess_tdof_list, opr);
 
-      // test that x^T S y = y^T S x
-      Vector smooth(n);
-      smooth = 0.0;
-      smoother->Mult(right, smooth);
-      double forward_val = left * smooth;
-      smoother->Mult(left, smooth);
-      double transpose_val = right * smooth;
+   Vector diag(fespace.GetTrueVSize());
+   aform.AssembleDiagonal(diag);
 
-      double error = fabs(forward_val - transpose_val) / fabs(forward_val);
-      CAPTURE(order, error);
-      REQUIRE(error < 1.e-13);
+   OperatorChebyshevSmoother smoother(*opr, diag, ess_tdof_list, cheb_order);
 
-      delete smoother;
-      delete fec;
-   }
+   const int n = smoother.Width();
+   Vector left(n);
+   Vector right(n);
+   left.Randomize(1);
+   right.Randomize(2);
+
+   // test that x^T S y = y^T S x
+   Vector smooth(n);
+   smoother.Mult(right, smooth);
+   real_t forward_val = left * smooth;
+
+   smoother.Mult(left, smooth);
+   real_t transpose_val = right * smooth;
+
+   real_t error = std::abs(forward_val - transpose_val) / std::abs(forward_val);
+   CAPTURE(order, error);
+   REQUIRE(error == MFEM_Approx(0.0));
 }