Vector::DeleteAt on device using InclusiveScan

# Conflicts:
#	mesh/mesh.cpp

.github/workflows/builds-and-tests.yml

@@ -168,10 +168,13 @@ jobs:
           env
         shell: bash
 
+      # For info on Xcode see:
+      #  - https://github.com/actions/runner-images/issues/12541
+      #  - https://github.com/actions/runner-images/blob/releases/macos-15-arm64/20250811/images/macos/macos-15-arm64-Readme.md#xcode
       - name: Xcode version setup (MacOS)
         if: matrix.os == 'macos-latest'
         run: |
-          XCODE_PATH="/Applications/Xcode_15.3.app"
+          XCODE_PATH="/Applications/Xcode_16.4.app"
           echo "> sudo xcode-select -s ${XCODE_PATH}"
           sudo xcode-select -s ${XCODE_PATH}
           echo "> g++ -v"

.gitignore

@@ -411,8 +411,6 @@ 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

CHANGELOG

@@ -80,13 +80,6 @@ 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,7 +78,6 @@ 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));
@@ -101,8 +100,7 @@ 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,7 +88,6 @@ 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.cpp

@@ -1275,6 +1275,22 @@ void BilinearForm::Update(FiniteElementSpace *nfes)
    height = width = fes->GetVSize();
 
    if (ext) { ext->Update(); }
+   for (int k = 0; k < domain_integs.Size(); ++k)
+   {
+      domain_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_integs[k]->Update();
+   }
+   for (int k = 0; k < interior_face_integs.Size(); ++k)
+   {
+      interior_face_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_face_integs[k]->Update();
+   }
 }
 
 void BilinearForm::SetDiagonalPolicy(DiagonalPolicy policy)
@@ -2337,6 +2353,31 @@ void MixedBilinearForm::Update()
    height = test_fes->GetVSize();
    width = trial_fes->GetVSize();
    if (ext) { ext->Update(); }
+
+   for (int k = 0; k < domain_integs.Size(); ++k)
+   {
+      domain_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_integs[k]->Update();
+   }
+   for (int k = 0; k < interior_face_integs.Size(); ++k)
+   {
+      interior_face_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_integs.Size(); ++k)
+   {
+      boundary_face_integs[k]->Update();
+   }
+   for (int k = 0; k < trace_face_integs.Size(); ++k)
+   {
+      trace_face_integs[k]->Update();
+   }
+   for (int k = 0; k < boundary_trace_face_integs.Size(); ++k)
+   {
+      boundary_trace_face_integs[k]->Update();
+   }
 }
 
 MixedBilinearForm::~MixedBilinearForm()

fem/bilinearform_ext.cpp

@@ -255,8 +255,6 @@ 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);
@@ -333,8 +331,6 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 
 void PABilinearFormExtension::Assemble()
 {
-   MFEM_PERF_FUNCTION;
-
    SetupRestrictionOperators(L2FaceValues::DoubleValued);
 
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
@@ -491,8 +487,6 @@ 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();
@@ -865,7 +859,6 @@ EABilinearFormExtension::EABilinearFormExtension(BilinearForm *form)
 
 void EABilinearFormExtension::Assemble()
 {
-   MFEM_PERF_FUNCTION;
    SetupRestrictionOperators(L2FaceValues::SingleValued);
 
    ne = trial_fes->GetMesh()->GetNE();
@@ -1414,7 +1407,6 @@ FABilinearFormExtension::FABilinearFormExtension(BilinearForm *form)
 
 void FABilinearFormExtension::Assemble()
 {
-   MFEM_PERF_FUNCTION;
    EABilinearFormExtension::Assemble();
    FiniteElementSpace &fes = *a->FESpace();
    int width = fes.GetVSize();

fem/bilininteg.cpp

@@ -21,6 +21,11 @@ using namespace std;
 namespace mfem
 {
 
+void BilinearFormIntegrator::Update()
+{
+   // default no-op
+}
+
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&)
 {
    MFEM_ABORT("BilinearFormIntegrator::AssemblePA(fes)\n"
@@ -3460,6 +3465,12 @@ real_t ElasticityIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
    return energy;
 }
 
+void DGTraceIntegrator::Update()
+{
+   qspace[0].reset();
+   qspace[1].reset();
+}
+
 void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
                                            const FiniteElement &el2,
                                            FaceElementTransformations &Trans,

fem/bilininteg.hpp

@@ -23,6 +23,8 @@
 
 namespace mfem
 {
+class QuadratureSpace;
+class FaceQuadratureSpace;
 
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator
@@ -44,6 +46,10 @@ public:
    // make sense for the action of the nonlinear operator (but they all make
    // sense for its Jacobian).
 
+   /// Signal this integrator that something about either the trial or test space has changed.
+   virtual void Update();
+
+
    /// Method defining partial assembly.
    /** The result of the partial assembly is stored internally so that it can be
        used later in the methods AddMultPA() and AddMultTransposePA(). */
@@ -3311,6 +3317,7 @@ protected:
    VectorCoefficient *u;
    real_t alpha, beta;
    // PA extension
+   std::unique_ptr<FaceQuadratureSpace> qspace[2];
    Vector pa_data;
    const DofToQuad *maps;             ///< Not owned
    const FaceGeometricFactors *geom;  ///< Not owned
@@ -3333,6 +3340,8 @@ public:
                      real_t a, real_t b)
    { rho = &rho_; u = &u_; alpha = a; beta = b; }
 
+   void Update() override;
+
    using BilinearFormIntegrator::AssembleFaceMatrix;
    void AssembleFaceMatrix(const FiniteElement &el1,
                            const FiniteElement &el2,

fem/coefficient.cpp

@@ -50,7 +50,6 @@ ElementTransformation *RefinedToCoarse(
 
 void Coefficient::Project(QuadratureFunction &qf)
 {
-   MFEM_PERF_FUNCTION;
    QuadratureSpaceBase &qspace = *qf.GetSpace();
    const int ne = qspace.GetNE();
    Vector values;

fem/fespace.cpp

@@ -101,10 +101,7 @@ FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
 FiniteElementSpace::FiniteElementSpace(Mesh *mesh,
                                        const FiniteElementCollection *fec,
                                        int vdim, int ordering)
-{
-   MFEM_PERF_FUNCTION;
-   Constructor(mesh, NULL, fec, vdim, ordering);
-}
+{ Constructor(mesh, NULL, fec, vdim, ordering); }
 
 FiniteElementSpace::FiniteElementSpace(Mesh *mesh, NURBSExtension *ext,
                                        const FiniteElementCollection *fec,
@@ -396,8 +393,6 @@ void FiniteElementSpace::BuildElementToDofTable() const
 {
    if (elem_dof) { return; }
 
-   MFEM_PERF_FUNCTION;
-
    // TODO: can we call GetElementDofs only once per element?
    Table *el_dof = new Table;
    Table *el_fos = (mesh->Dimension() > 2) ? (new Table) : NULL;
@@ -2753,8 +2748,6 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const
 
 void FiniteElementSpace::Construct()
 {
-   MFEM_PERF_FUNCTION;
-
    // This method should be used only for non-NURBS spaces.
    MFEM_VERIFY(!NURBSext, "internal error");
 

fem/gridfunc.cpp

@@ -19,7 +19,6 @@
 #include "../mesh/nurbs.hpp"
 #include "../mesh/vtkhdf.hpp"
 #include "../general/text.hpp"
-#include "../general/reducers.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -3327,126 +3326,8 @@ real_t GridFunction::ComputeLpError(const real_t p, Coefficient &exsol,
                                     const IntegrationRule *irs[],
                                     const Array<int> *elems) const
 {
-   MFEM_PERF_FUNCTION;
-
-   MFEM_VERIFY(fes->GetVDim() == 1, "invalid vector dimension!");
-
    real_t error = 0.0;
-
-   bool device_eval = true;
-   // TODO: check for cases that are not supported on device:
-   //       * mixed meshes
-   //       * meshes with non-tensor-product elements can have negative weights
-   //       * variable orders
-   //       * weight is not NULL
-   //       * elems is not NULL
-   //       * map type is not VALUE
-   //       * ...
-   Mesh *mesh = fes->GetMesh();
-   const FiniteElement *fe = fes->GetTypicalFE();
-   if (mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
-       (mesh->Dimension() > 1 && mesh->MeshGenerator() != 2) ||
-       fes->IsVariableOrder() ||
-       weight != nullptr ||
-       elems != nullptr ||
-       fe->GetMapType() != FiniteElement::MapType::VALUE)
-   {
-      device_eval = false;
-   }
-   if (device_eval)
-   {
-      Geometry::Type geom = mesh->GetTypicalElementGeometry();
-      const IntegrationRule *ir_p;
-      if (irs)
-      {
-         ir_p = irs[geom];
-      }
-      else
-      {
-         int intorder = 2*fe->GetOrder() + 3; // <----------
-         ir_p = &(IntRules.Get(geom, intorder));
-      }
-      const IntegrationRule &ir = *ir_p;
-      QuadratureSpace qs(*mesh, ir);
-      CoefficientVector coeff(exsol, qs, CoefficientStorage::FULL);
-
-      const QVectorLayout ql = QVectorLayout::byNODES;
-      const MemoryType d_mt = MemoryType::DEFAULT;
-      Vector q_vals;
-      // TODO: make this a method
-      {
-         // const FiniteElement *fe = fes->GetTypicalFE();
-         const int vdim = fes->GetVDim();
-         const int NE   = fes->GetNE();
-         const int ND   = fe->GetDof();
-         const int NQ   = ir.GetNPoints();
-         MemoryType my_d_mt = (d_mt != MemoryType::DEFAULT) ? d_mt :
-                              Device::GetDeviceMemoryType();
-         // byNODES : NQPT x VDIM x NE
-         // byVDIM  : VDIM x NQPT x NE
-         q_vals.SetSize(vdim*NQ*NE, my_d_mt);
-         const QuadratureInterpolator &qi = *fes->GetQuadratureInterpolator(ir);
-         qi.SetOutputLayout(ql);
-         const bool use_tensor_products = UsesTensorBasis(*fes);
-         qi.DisableTensorProducts(!use_tensor_products);
-         const ElementDofOrdering e_ordering =
-            use_tensor_products ?
-            ElementDofOrdering::LEXICOGRAPHIC :
-            ElementDofOrdering::NATIVE;
-         const Operator *elem_restr = fes->GetElementRestriction(e_ordering);
-         if (fe->GetMapType() == FiniteElement::MapType::INTEGRAL)
-         {
-            // Pre-compute the geometric factors in order to set the desired
-            // MemoryType they use:
-            fes->GetMesh()->GetGeometricFactors(
-               ir, GeometricFactors::DETERMINANTS, my_d_mt);
-         }
-         if (elem_restr)
-         {
-            Vector f_e(vdim*ND*NE, my_d_mt);
-            elem_restr->Mult(*this, f_e);
-            qi.PhysValues(f_e, q_vals);
-         }
-         else
-         {
-            qi.PhysValues(*this, q_vals);
-         }
-      }
-
-      const real_t *exact_d = coeff.Read();
-      const real_t *gridf_d = q_vals.Read();
-      // FIXME: reuse the workspace vector from vector.cpp?
-      static Array<real_t> workspace;
-      if (p < infinity())
-      {
-         MemoryType my_d_mt = (d_mt != MemoryType::DEFAULT) ? d_mt :
-                              Device::GetDeviceMemoryType();
-         const GeometricFactors *geom_factors =
-            fes->GetMesh()->GetGeometricFactors(
-               ir, GeometricFactors::DETERMINANTS, my_d_mt);
-         const real_t *detJ_d = geom_factors->detJ.Read();
-         const real_t *w_d = ir.GetWeights().Read();
-         const int NQ = ir.GetNPoints();
-         mfem::reduce(q_vals.Size(), error,
-                      [=] MFEM_HOST_DEVICE(int i, real_t &r)
-         {
-            const real_t diff = fabs(exact_d[i] - gridf_d[i]);
-            r += w_d[i%NQ] * detJ_d[i] * pow(diff, p);
-         }, SumReducer<real_t> {}, true, workspace);
-         error = pow(error, 1./p);
-      }
-      else
-      {
-         mfem::reduce(q_vals.Size(), error,
-                      [=] MFEM_HOST_DEVICE(int i, real_t &r)
-         {
-            const real_t diff = fabs(exact_d[i] - gridf_d[i]);
-            r = fmax(r, diff);
-         }, MaxReducer<real_t> {}, true, workspace);
-      }
-      return error;
-   }
-
+   const FiniteElement *fe;
    ElementTransformation *T;
    Vector vals;
 

fem/integ/bilininteg_dgtrace_pa.cpp

@@ -139,8 +139,6 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
 
-   nf = fes.GetNFbyType(type);
-   if (nf==0) { return; }
    // Assumes tensor-product elements
    Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetTypicalTraceElement();
@@ -148,6 +146,17 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
                                IntRule:
                                &GetRule(el.GetGeomType(), el.GetOrder(),
                                         *mesh->GetTypicalElementTransformation());
+
+   if (!qspace[static_cast<int>(type)])
+   {
+      qspace[static_cast<int>(type)].reset(
+         new FaceQuadratureSpace(*mesh, *ir, type));
+   }
+
+   FaceQuadratureSpace& qs = *qspace[static_cast<int>(type)];
+   nf = qs.GetNumFaces();
+   if (nf==0) { return; }
+
    const int symmDims = 4;
    nq = ir->GetNPoints();
    dim = mesh->Dimension();
@@ -159,8 +168,6 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
-
-   FaceQuadratureSpace qs(*mesh, *ir, type);
    CoefficientVector vel(*u, qs, CoefficientStorage::COMPRESSED);
 
    CoefficientVector r(qs, CoefficientStorage::COMPRESSED);

fem/integ/bilininteg_diffusion_ea.cpp

@@ -161,8 +161,7 @@ static void EADiffusionAssemble3D(const int NE,
    auto B = Reshape(b.Read(), Q1D, D1D);
    auto G = Reshape(g.Read(), Q1D, D1D);
    auto D = Reshape(padata.Read(), Q1D, Q1D, Q1D, 6, NE);
-   auto A = Reshape(add ? eadata.ReadWrite() : eadata.Write(),
-                    D1D, D1D, D1D, D1D, D1D, D1D, NE);
+   auto A = Reshape(eadata.ReadWrite(), D1D, D1D, D1D, D1D, D1D, D1D, NE);
    mfem::forall_3D(NE, D1D, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
    {
       const int D1D = T_D1D ? T_D1D : d1d;
@@ -247,60 +246,58 @@ void DiffusionIntegrator::AssembleEA(const FiniteElementSpace &fes,
                                      Vector &ea_data,
                                      const bool add)
 {
-   MFEM_PERF_FUNCTION;
    AssemblePA(fes);
    ne = fes.GetMesh()->GetNE();
    const Array<real_t> &B = maps->B;
    const Array<real_t> &G = maps->G;
-   decltype(&EADiffusionAssemble1D<>) kernel = nullptr;
    if (dim == 1)
    {
       switch ((dofs1D << 4 ) | quad1D)
       {
-         case 0x22: kernel = EADiffusionAssemble1D<2,2>;
-         case 0x33: kernel = EADiffusionAssemble1D<3,3>;
-         case 0x44: kernel = EADiffusionAssemble1D<4,4>;
-         case 0x55: kernel = EADiffusionAssemble1D<5,5>;
-         case 0x66: kernel = EADiffusionAssemble1D<6,6>;
-         case 0x77: kernel = EADiffusionAssemble1D<7,7>;
-         case 0x88: kernel = EADiffusionAssemble1D<8,8>;
-         case 0x99: kernel = EADiffusionAssemble1D<9,9>;
-         default:   kernel = EADiffusionAssemble1D<>;
+         case 0x22: return EADiffusionAssemble1D<2,2>(ne,B,G,pa_data,ea_data,add);
+         case 0x33: return EADiffusionAssemble1D<3,3>(ne,B,G,pa_data,ea_data,add);
+         case 0x44: return EADiffusionAssemble1D<4,4>(ne,B,G,pa_data,ea_data,add);
+         case 0x55: return EADiffusionAssemble1D<5,5>(ne,B,G,pa_data,ea_data,add);
+         case 0x66: return EADiffusionAssemble1D<6,6>(ne,B,G,pa_data,ea_data,add);
+         case 0x77: return EADiffusionAssemble1D<7,7>(ne,B,G,pa_data,ea_data,add);
+         case 0x88: return EADiffusionAssemble1D<8,8>(ne,B,G,pa_data,ea_data,add);
+         case 0x99: return EADiffusionAssemble1D<9,9>(ne,B,G,pa_data,ea_data,add);
+         default:   return EADiffusionAssemble1D(ne,B,G,pa_data,ea_data,add,
+                                                    dofs1D,quad1D);
       }
    }
    else if (dim == 2)
    {
       switch ((dofs1D << 4 ) | quad1D)
       {
-         case 0x22: kernel = EADiffusionAssemble2D<2,2>;
-         case 0x33: kernel = EADiffusionAssemble2D<3,3>;
-         case 0x44: kernel = EADiffusionAssemble2D<4,4>;
-         case 0x55: kernel = EADiffusionAssemble2D<5,5>;
-         case 0x66: kernel = EADiffusionAssemble2D<6,6>;
-         case 0x77: kernel = EADiffusionAssemble2D<7,7>;
-         case 0x88: kernel = EADiffusionAssemble2D<8,8>;
-         case 0x99: kernel = EADiffusionAssemble2D<9,9>;
-         default:   kernel = EADiffusionAssemble2D<>;
+         case 0x22: return EADiffusionAssemble2D<2,2>(ne,B,G,pa_data,ea_data,add);
+         case 0x33: return EADiffusionAssemble2D<3,3>(ne,B,G,pa_data,ea_data,add);
+         case 0x44: return EADiffusionAssemble2D<4,4>(ne,B,G,pa_data,ea_data,add);
+         case 0x55: return EADiffusionAssemble2D<5,5>(ne,B,G,pa_data,ea_data,add);
+         case 0x66: return EADiffusionAssemble2D<6,6>(ne,B,G,pa_data,ea_data,add);
+         case 0x77: return EADiffusionAssemble2D<7,7>(ne,B,G,pa_data,ea_data,add);
+         case 0x88: return EADiffusionAssemble2D<8,8>(ne,B,G,pa_data,ea_data,add);
+         case 0x99: return EADiffusionAssemble2D<9,9>(ne,B,G,pa_data,ea_data,add);
+         default:   return EADiffusionAssemble2D(ne,B,G,pa_data,ea_data,add,
+                                                    dofs1D,quad1D);
       }
    }
    else if (dim == 3)
    {
       switch ((dofs1D << 4 ) | quad1D)
       {
-         case 0x23: kernel = EADiffusionAssemble3D<2,3>;
-         case 0x34: kernel = EADiffusionAssemble3D<3,4>;
-         case 0x45: kernel = EADiffusionAssemble3D<4,5>;
-         case 0x56: kernel = EADiffusionAssemble3D<5,6>;
-         case 0x67: kernel = EADiffusionAssemble3D<6,7>;
-         case 0x78: kernel = EADiffusionAssemble3D<7,8>;
-         case 0x89: kernel = EADiffusionAssemble3D<8,9>;
-         default:   kernel = EADiffusionAssemble3D<>;
+         case 0x23: return EADiffusionAssemble3D<2,3>(ne,B,G,pa_data,ea_data,add);
+         case 0x34: return EADiffusionAssemble3D<3,4>(ne,B,G,pa_data,ea_data,add);
+         case 0x45: return EADiffusionAssemble3D<4,5>(ne,B,G,pa_data,ea_data,add);
+         case 0x56: return EADiffusionAssemble3D<5,6>(ne,B,G,pa_data,ea_data,add);
+         case 0x67: return EADiffusionAssemble3D<6,7>(ne,B,G,pa_data,ea_data,add);
+         case 0x78: return EADiffusionAssemble3D<7,8>(ne,B,G,pa_data,ea_data,add);
+         case 0x89: return EADiffusionAssemble3D<8,9>(ne,B,G,pa_data,ea_data,add);
+         default:   return EADiffusionAssemble3D(ne,B,G,pa_data,ea_data,add,
+                                                    dofs1D,quad1D);
       }
    }
-   MFEM_VERIFY(kernel != nullptr, "Unknown kernel.");
-   kernel(ne,B,G,pa_data,ea_data,add,dofs1D,quad1D);
-   // Free the PA data:
-   pa_data.Destroy();
+   MFEM_ABORT("Unknown kernel.");
 }
 
 }

fem/integ/bilininteg_diffusion_pa.cpp

@@ -39,8 +39,6 @@ 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);
@@ -90,8 +88,6 @@ 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_pa.cpp

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

fem/integ/lininteg_domain.cpp

@@ -242,8 +242,6 @@ void DomainLFIntegrator::AssembleDevice(const FiniteElementSpace &fes,
                                         const Array<int> &markers,
                                         Vector &b)
 {
-   MFEM_PERF_FUNCTION;
-
    const FiniteElement &fe = *fes.GetTypicalFE();
    const int qorder = oa * fe.GetOrder() + ob;
    const Geometry::Type gtype = fe.GetGeomType();

fem/linearform.cpp

@@ -161,7 +161,6 @@ bool LinearForm::SupportsDevice() const
 
 void LinearForm::UseFastAssembly(bool use_fa)
 {
-   MFEM_PERF_FUNCTION;
    fast_assembly = use_fa;
 
    if (fast_assembly && SupportsDevice() && !ext)
@@ -172,8 +171,6 @@ void LinearForm::UseFastAssembly(bool use_fa)
 
 void LinearForm::Assemble()
 {
-   MFEM_PERF_FUNCTION;
-
    Array<int> vdofs;
    ElementTransformation *eltrans;
    Vector elemvect;

fem/linearform_ext.cpp

@@ -15,16 +15,10 @@
 namespace mfem
 {
 
-LinearFormExtension::LinearFormExtension(LinearForm *lf): lf(lf)
-{
-   MFEM_PERF_FUNCTION;
-   Update();
-}
+LinearFormExtension::LinearFormExtension(LinearForm *lf): lf(lf) { Update(); }
 
 void LinearFormExtension::Assemble()
 {
-   MFEM_PERF_FUNCTION;
-
    const FiniteElementSpace &fes = *lf->FESpace();
    MFEM_VERIFY(lf->SupportsDevice(), "Not supported.");
    MFEM_VERIFY(lf->Size() == fes.GetVSize(), "LinearForm size does not "
@@ -119,8 +113,6 @@ void LinearFormExtension::Assemble()
 
 void LinearFormExtension::Update()
 {
-   MFEM_PERF_FUNCTION;
-
    const FiniteElementSpace &fes = *lf->FESpace();
    const Mesh &mesh = *fes.GetMesh();
    constexpr ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;

fem/linearform_ext.hpp

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

fem/lor/lor.cpp

@@ -365,8 +365,6 @@ 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,8 +360,6 @@ 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
@@ -419,8 +417,6 @@ 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,8 +184,6 @@ 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

@@ -10,7 +10,6 @@
 // CONTRIBUTING.md for details.
 
 #include "multigrid.hpp"
-#include "../general/annotation.hpp"
 
 namespace mfem
 {
@@ -18,10 +17,7 @@ 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_,
@@ -33,18 +29,12 @@ 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])
@@ -66,17 +56,16 @@ void MultigridBase::InitVectors() const
    X.SetSize(M, nrhs);
    Y.SetSize(M, nrhs);
    R.SetSize(M, nrhs);
-   for (int i = 0; i < M; ++i)
+   Z.SetSize(M, nrhs);
+   for (int i = 0; i < X.NumRows(); ++i)
    {
       const int n = operators[i]->Height();
-      for (int j = 0; j < nrhs; ++j)
+      for (int j = 0; j < X.NumCols(); ++j)
       {
-         if (i < M - 1)
-         {
-            X(i, j) = new Vector(n);
-            Y(i, j) = new Vector(n);
-         }
+         X(i, j) = new Vector(n);
+         Y(i, j) = new Vector(n);
          R(i, j) = new Vector(n);
+         Z(i, j) = new Vector(n);
       }
    }
 }
@@ -87,12 +76,10 @@ void MultigridBase::EraseVectors() const
    {
       for (int j = 0; j < X.NumCols(); ++j)
       {
-         if (i < X.NumRows() - 1)
-         {
-            delete X(i, j);
-            delete Y(i, j);
-         }
+         delete X(i, j);
+         delete Y(i, j);
          delete R(i, j);
+         delete Z(i, j);
       }
    }
 }
@@ -108,12 +95,6 @@ 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_)
 {
@@ -124,24 +105,25 @@ void MultigridBase::SetCycleType(CycleType cycleType_, int preSmoothingSteps_,
 
 void MultigridBase::Mult(const Vector& x, Vector& y) const
 {
-   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
-
+   Array<const Vector*> X_(1);
+   Array<Vector*> Y_(1);
+   X_[0] = &x;
+   Y_[0] = &y;
    ArrayMult(X_, Y_);
 }
 
 void MultigridBase::ArrayMult(const Array<const Vector*>& X_,
                               Array<Vector*>& Y_) const
 {
-   MFEM_PERF_FUNCTION;
-
    MFEM_ASSERT(operators.Size() > 0,
                "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();
@@ -152,163 +134,96 @@ 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) = const_cast<Vector*>(X_[j]);
-      Y(M - 1, j) = Y_[j];
+      *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);
    }
-   const bool zero = !iterative_mode;
-   Cycle(M - 1, zero);
 }
 
 void MultigridBase::SmoothingStep(int level, bool zero, bool transpose) const
 {
-   MFEM_PERF_FUNCTION;
-
    // 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)
    {
-      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_);
+      Array<Vector *> X_(X[level], nrhs), Y_(Y[level], nrhs);
+      GetSmootherAtLevel(level)->ArrayMult(X_, Y_);
    }
    else
    {
-      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_);
+      Array<Vector *> Y_(Y[level], nrhs), R_(R[level], nrhs),
+            Z_(Z[level], nrhs);
       for (int j = 0; j < nrhs; ++j)
       {
-         // *R_[j] = *X(level, j) - *R_[j]
-         subtract(*X(level, j), *R_[j], *R_[j]);
+         *R_[j] = *X(level, j);
       }
+      GetOperatorAtLevel(level)->ArrayAddMult(Y_, R_, -1.0);
       if (transpose)
       {
-         GetSmootherAtLevel(level)->ArrayAddMultTranspose(cR_, Y_);
+         GetSmootherAtLevel(level)->ArrayMultTranspose(R_, Z_);
       }
       else
       {
-         GetSmootherAtLevel(level)->ArrayAddMult(cR_, Y_);
+         GetSmootherAtLevel(level)->ArrayMult(R_, Z_);
       }
-   }
-}
-
-void MultigridBase::CoarseSolve(bool zero) const
-{
-   MFEM_PERF_FUNCTION;
-
-   // See the comment about iterative mode in SmoothingStep()
-   coarse_solver->iterative_mode = false;
-
-   if (zero)
-   {
-      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]);
+         *Y_[j] += *Z_[j];
       }
-      coarse_solver->ArrayAddMult(cR_, Y_);
    }
 }
 
-void MultigridBase::Cycle(int level, bool zero) const
+void MultigridBase::Cycle(int level) 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)
+   if (level == 0)
    {
-      SmoothingStep(0, zero, false);
+      SmoothingStep(0, true, false);
       return;
    }
 
    // Pre-smooth
    for (int i = 0; i < preSmoothingSteps; ++i)
    {
-      SmoothingStep(level, zero && (i == 0), false);
-   }
-
-   // Coarse solve with 'coarse_solver'
-   if (level == 0)
-   {
-      CoarseSolve(preSmoothingSteps == 0 && zero);
-      goto mg_post_smooth;
+      SmoothingStep(level, (cycleType == CycleType::VCYCLE && i == 0), false);
    }
 
    // Compute residual and restrict
-   if (preSmoothingSteps == 0 && zero)
    {
-      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_);
+      Array<Vector *> Y_(Y[level], nrhs), R_(R[level], nrhs),
+            X_(X[level - 1], nrhs);
       for (int j = 0; j < nrhs; ++j)
       {
-         // *R_[j] = *X(level, j) - *R_[j]
-         subtract(*X(level, j), *R_[j], *R_[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;
       }
-      GetProlongationAtLevel(level - 1)->ArrayMultTranspose(cR_, X_);
    }
 
    // Corrections
-   Cycle(level - 1, true);
+   Cycle(level - 1);
    if (cycleType == CycleType::WCYCLE)
    {
-      // 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); }
+      Cycle(level - 1);
    }
 
    // Prolongate and add
    {
-      const Array<const Vector *> cY_lm1((const Vector **)(Y[level - 1]), nrhs);
-      Array<Vector *> Y_l(Y[level], nrhs);
-
-      if (preSmoothingSteps == 0 && zero)
+      Array<Vector *> Y_(Y[level - 1], nrhs), Z_(Z[level], nrhs);
+      GetProlongationAtLevel(level - 1)->ArrayMult(Y_, Z_);
+      for (int j = 0; j < nrhs; ++j)
       {
-         GetProlongationAtLevel(level - 1)->ArrayMult(cY_lm1, Y_l);
-      }
-      else
-      {
-         GetProlongationAtLevel(level - 1)->ArrayAddMult(cY_lm1, Y_l);
+         *Y(level, j) += *Z_[j];
       }
    }
 
-mg_post_smooth:
    // Post-smooth
    for (int i = 0; i < postSmoothingSteps; ++i)
    {

fem/multigrid.hpp

@@ -36,14 +36,12 @@ 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;
+   mutable Array2D<Vector*> X, Y, R, Z;
    mutable int nrhs;
 
 public:
@@ -67,16 +65,6 @@ 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(); }
 
@@ -130,14 +118,11 @@ public:
 
 private:
    /// Application of a multigrid cycle at particular level
-   void Cycle(int level, bool zero) const;
+   void Cycle(int level) 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;
@@ -217,9 +202,6 @@ 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,8 +124,6 @@ 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

@@ -63,11 +63,9 @@ ParFiniteElementSpace::ParFiniteElementSpace(
 
 ParFiniteElementSpace::ParFiniteElementSpace(
    ParMesh *pm, const FiniteElementCollection *f, int dim, int ordering)
-   : FiniteElementSpace((MFEM_PERF_BEGIN(_MFEM_FUNC_NAME), pm),
-                        f, dim, ordering)
+   : FiniteElementSpace(pm, f, dim, ordering)
 {
    ParInit(pm);
-   MFEM_PERF_END(_MFEM_FUNC_NAME);
 }
 
 ParFiniteElementSpace::ParFiniteElementSpace(
@@ -94,7 +92,6 @@ ParNURBSExtension *ParFiniteElementSpace::MakeLocalNURBSext(
 
 void ParFiniteElementSpace::ParInit(ParMesh *pm)
 {
-   MFEM_PERF_FUNCTION;
    pmesh = pm;
    pncmesh = nullptr;
 
@@ -184,7 +181,6 @@ void ParFiniteElementSpace::CommunicateGhostOrder()
 
 void ParFiniteElementSpace::Construct()
 {
-   MFEM_PERF_FUNCTION;
    if (NURBSext)
    {
       ConstructTrueNURBSDofs();
@@ -843,8 +839,6 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
 
    if (P) { return; }
 
-   MFEM_PERF_FUNCTION;
-
    if (!nd_strias)
    {
       // Safe to assume 1-1 correspondence between shared dofs
@@ -1430,7 +1424,6 @@ const Operator *ParFiniteElementSpace::GetRestrictionOperator() const
       if (NRanks == 1)
       {
          R_transpose.reset(new IdentityOperator(GetTrueVSize()));
-         Rconf = new IdentityOperator(GetTrueVSize());
       }
       else
       {
@@ -1443,8 +1436,8 @@ const Operator *ParFiniteElementSpace::GetRestrictionOperator() const
             R_transpose.reset(
                new DeviceConformingProlongationOperator(*this, true));
          }
-         Rconf = new TransposeOperator(*R_transpose);
       }
+      Rconf = new TransposeOperator(*R_transpose);
       return Rconf;
    }
    else

fem/plinearform.cpp

@@ -45,8 +45,6 @@ void ParLinearForm::MakeRef(ParFiniteElementSpace *pf, Vector &v, int v_offset)
 
 void ParLinearForm::Assemble()
 {
-   MFEM_PERF_FUNCTION;
-
    LinearForm::Assemble();
 
    if (interior_face_integs.Size())

fem/qspace.cpp

@@ -17,8 +17,9 @@ namespace mfem
 {
 
 QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                                         const IntegrationRule &ir)
-   : mesh(mesh_), order(ir.GetOrder())
+                                         const IntegrationRule &ir,
+                                         QSpaceStorage storage)
+   : mesh(mesh_), order(ir.GetOrder()), storage(storage)
 {
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
@@ -96,11 +97,10 @@ void QuadratureSpaceBase::Integrate(VectorCoefficient &coeff,
 
 void QuadratureSpace::ConstructOffsets()
 {
-   MFEM_PERF_FUNCTION;
-   const int num_elem = mesh.GetNE();
-   ne = num_elem;
+   const int num_elem = ne;
 
-   if (mesh.GetNumGeometries(mesh.Dimension()) == 1)
+   if (storage == QSpaceStorage::COMPRESSED &&
+       mesh.GetNumGeometries(mesh.Dimension()) == 1)
    {
       Array<Geometry::Type> geoms;
       mesh.GetGeometries(mesh.Dimension(), geoms);
@@ -125,14 +125,9 @@ void QuadratureSpace::ConstructOffsets()
    }
 }
 
-void QuadratureSpace::Construct()
-{
-   ConstructIntRules(mesh.Dimension());
-   ConstructOffsets();
-}
-
-QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
-   : QuadratureSpaceBase(*mesh_)
+QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in,
+                                 QSpaceStorage storage)
+   : QuadratureSpaceBase(*mesh_, 0, storage)
 {
    const char *msg = "invalid input stream";
    std::string ident;
@@ -151,15 +146,24 @@ QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
       return;
    }
 
-   Construct();
+   ne = mesh.GetNE();
+   ConstructIntRules(mesh.Dimension());
 }
 
-QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir)
+QuadratureSpace::QuadratureSpace(Mesh *mesh_, int order_, QSpaceStorage storage)
+   : QuadratureSpaceBase(*mesh_, order_, storage)
+{
+   ne = mesh.GetNE();
+   ConstructIntRules(mesh.Dimension());
+}
+
+QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
+                                 QSpaceStorage storage)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir, storage)
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) <= 1,
                "Constructor not valid for mixed meshes");
-   ConstructOffsets();
+   ne = mesh.GetNE();
 }
 
 void QuadratureSpace::Save(std::ostream &os) const
@@ -181,55 +185,53 @@ const Vector &QuadratureSpace::GetGeometricFactorWeights() const
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, int order_,
-                                         FaceType face_type_)
-   : QuadratureSpaceBase(mesh_, order_),
-     face_type(face_type_),
-     num_faces(mesh.GetNFbyType(face_type))
+                                         FaceType face_type_,
+                                         QSpaceStorage storage)
+   : QuadratureSpaceBase(mesh_, order_, storage), face_type(face_type_),
+     face_indices(mesh.GetFaceIndices(face_type_)),
+     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
 {
-   Construct();
+   ne = face_indices.Size();
+   ConstructIntRules(mesh.Dimension() - 1);
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                                         FaceType face_type_)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir),
-     face_type(face_type_),
-     num_faces(mesh.GetNFbyType(face_type))
+                                         FaceType face_type_,
+                                         QSpaceStorage storage)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir, storage),
+     face_type(face_type_), face_indices(mesh.GetFaceIndices(face_type_)),
+     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension() - 1) <= 1,
                "Constructor not valid for mixed meshes");
-   ConstructOffsets();
+   ne = face_indices.Size();
 }
 
 void FaceQuadratureSpace::ConstructOffsets()
 {
-   face_indices.SetSize(num_faces);
-   offsets.SetSize(num_faces + 1);
-   ne = num_faces;
-   int offset = 0;
-   int f_idx = 0;
-   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
+   if (storage == QSpaceStorage::COMPRESSED &&
+       mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
    {
-      const Mesh::FaceInformation face = mesh.GetFaceInformation(i);
-      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(face_type))
-      {
-         continue;
-      }
-      face_indices[f_idx] = i;
-      face_indices_inv[i] = f_idx;
-      offsets[f_idx] = offset;
-      Geometry::Type geom = mesh.GetFaceGeometry(i);
-      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
-      offset += int_rule[geom]->GetNPoints();
-
-      f_idx++;
+      Array<Geometry::Type> geoms;
+      mesh.GetGeometries(mesh.Dimension() - 1, geoms);
+      offsets.SetSize(1);
+      offsets.HostWrite();
+      offsets[0] = int_rule[geoms[0]]->GetNPoints();
+      size = ne * offsets[0];
+   }
+   else
+   {
+      offsets.SetSize(face_indices.Size() + 1);
+      int offset = 0;
+      for (int i = 0; i < mesh.GetNFbyType(face_type); ++i)
+      {
+         offsets[i] = offset;
+         Geometry::Type geom = mesh.GetFaceGeometry(face_indices[i]);
+         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
+         offset += int_rule[geom]->GetNPoints();
+      }
+      offsets[face_indices.Size()] = size = offset;
    }
-   offsets[num_faces] = size = offset;
-}
-
-void FaceQuadratureSpace::Construct()
-{
-   ConstructIntRules(mesh.Dimension() - 1);
-   ConstructOffsets();
 }
 
 int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const

fem/qspace.hpp

@@ -19,39 +19,49 @@
 namespace mfem
 {
 
+enum class QSpaceStorage
+{
+   FULL,
+   COMPRESSED
+};
+
 /// Abstract base class for QuadratureSpace and FaceQuadratureSpace.
 /** This class represents the storage layout for QuadratureFunction%s, that may
     be defined either on mesh elements or mesh faces. */
 class QuadratureSpaceBase
 {
 protected:
-   friend class QuadratureFunction; // Uses the offsets.
-
    Mesh &mesh; ///< The underlying mesh.
    int order; ///< The order of integration rule.
-   int size; ///< Total number of quadrature points.
+   int size = -1; ///< Total number of quadrature points. -1 indicates
+   ///< offsets/size not computed yet.
    int ne; ///< Actual number of entities
    mutable Vector weights; ///< Integration weights.
    mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.
 
+   QSpaceStorage storage;
+
    /// @brief Entity quadrature point offset array.
    ///
    /// Supports a constant compression scheme for meshes which have a single
    /// geometry type. When compressed, will have a single value. The true offset
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
-   /// has size num_entities + 1.
+   /// has size num_entities + 1. Lazily constructed.
    ///
    Array<int> offsets;
    /// The quadrature rules used for each geometry type.
    const IntegrationRule *int_rule[Geometry::NumGeom];
 
    /// Protected constructor. Used by derived classes.
-   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0)
-      : mesh(mesh_), order(order_) { }
+   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED)
+      : mesh(mesh_), order(order_), storage(storage)
+   {}
 
    /// Protected constructor. Used by derived classes.
    QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                       const IntegrationRule &ir);
+                       const IntegrationRule &ir,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Fill the @ref int_rule array for each geometry type using @ref order.
    void ConstructIntRules(int dim);
@@ -62,13 +72,21 @@ protected:
    /// Compute the integration weights.
    void ConstructWeights() const;
 
+   virtual void ConstructOffsets() = 0;
+
 public:
+   QSpaceStorage StorageType() const { return storage; }
+
    /// @brief Gets the offset for a given entity @a idx.
    ///
    /// The quadrature point values for entity i are stored in the indices
    /// between Offset(i) and Offset(i+1)
    int Offset(int idx) const
    {
+      if (size < 0)
+      {
+         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
+      }
       return (offsets.Size() == 1) ? (idx * offsets[0]) : offsets[idx];
    }
 
@@ -79,10 +97,24 @@ public:
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
    /// has size num_entities + 1.
    ///
-   const Array<int> &Offsets() const { return offsets; }
+   const Array<int> &Offsets() const
+   {
+      if (size < 0)
+      {
+         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
+      }
+      return offsets;
+   }
 
    /// Return the total number of quadrature points.
-   int GetSize() const { return size; }
+   int GetSize() const
+   {
+      if (size < 0)
+      {
+         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
+      }
+      return size;
+   }
 
    /// Return the order of the quadrature rule(s) used by all elements.
    int GetOrder() const { return order; }
@@ -142,19 +174,20 @@ class QuadratureSpace : public QuadratureSpaceBase
 {
 protected:
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets();
-   void Construct();
+   void ConstructOffsets() override;
 public:
    /// Create a QuadratureSpace based on the global rules from #IntRules.
-   QuadratureSpace(Mesh *mesh_, int order_)
-      : QuadratureSpaceBase(*mesh_, order_) { Construct(); }
+   QuadratureSpace(Mesh *mesh_, int order_,
+                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// @brief Create a QuadratureSpace with an IntegrationRule, valid only when
    /// the mesh has one element type.
-   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir);
+   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
+                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Read a QuadratureSpace from the stream @a in.
-   QuadratureSpace(Mesh *mesh_, std::istream &in);
+   QuadratureSpace(Mesh *mesh_, std::istream &in,
+                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Returns number of elements in the mesh.
    inline int GetNE() const { return mesh.GetNE(); }
@@ -191,29 +224,29 @@ public:
 class FaceQuadratureSpace : public QuadratureSpaceBase
 {
    FaceType face_type; ///< Is the space defined on interior or boundary faces?
-   const int num_faces; ///< Number of faces.
 
    /// Map from boundary or interior face indices to mesh face indices.
-   Array<int> face_indices;
+   const Array<int> &face_indices;
 
    /// Inverse of the map @a face_indices.
-   std::unordered_map<int,int> face_indices_inv;
+   const std::unordered_map<int,int> &face_indices_inv;
 
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets();
-   void Construct();
+   void ConstructOffsets() override;
 
 public:
    /// Create a FaceQuadratureSpace based on the global rules from #IntRules.
-   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_);
+   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// @brief Create a FaceQuadratureSpace with an IntegrationRule, valid only
    /// when the mesh has one type of face geometry.
    FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                       FaceType face_type_);
+                       FaceType face_type_,
+                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
 
    /// Returns number of faces in the mesh.
-   inline int GetNumFaces() const { return num_faces; }
+   inline int GetNumFaces() const { return face_indices.Size(); }
 
    /// Returns the face type (boundary or interior).
    FaceType GetFaceType() const { return face_type; }

fem/quadinterpolator.cpp

@@ -503,7 +503,6 @@ 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/restriction.cpp

@@ -25,7 +25,7 @@ namespace mfem
 
 ElementRestriction::ElementRestriction(const FiniteElementSpace &f,
                                        ElementDofOrdering e_ordering)
-   : fes((MFEM_PERF_BEGIN(_MFEM_FUNC_NAME), f)),
+   : fes(f),
      ne(fes.GetNE()),
      vdim(fes.GetVDim()),
      byvdim(fes.GetOrdering() == Ordering::byVDIM),
@@ -104,13 +104,10 @@ ElementRestriction::ElementRestriction(const FiniteElementSpace &f,
       offsets[i] = offsets[i - 1];
    }
    offsets[0] = 0;
-   MFEM_PERF_END(_MFEM_FUNC_NAME);
 }
 
 void ElementRestriction::Mult(const Vector& x, Vector& y) const
 {
-   MFEM_PERF_FUNCTION;
-
    // Assumes all elements have the same number of dofs
    const int nd = dof;
    const int vd = vdim;
@@ -155,8 +152,6 @@ void ElementRestriction::AbsMult(const Vector& x, Vector& y) const
 template <bool ADD>
 void ElementRestriction::TAddMultTranspose(const Vector& x, Vector& y) const
 {
-   MFEM_PERF_FUNCTION;
-
    // Assumes all elements have the same number of dofs
    const int nd = dof;
    const int vd = vdim;

fem/transfer.cpp

@@ -13,7 +13,6 @@
 #include "bilinearform.hpp"
 #include "pbilinearform.hpp"
 #include "../general/forall.hpp"
-#include "kernels.hpp"
 
 namespace mfem
 {
@@ -2323,76 +2322,6 @@ 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)
-{
-   MFEM_PERF_FUNCTION;
-
-   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)
@@ -2474,9 +2403,9 @@ void Prolongation3D(const int NE, const int D1D, const int Q1D,
    });
 }
 
-void ProlongationTranspose2D(const int NE, const int D1D, const int Q1D,
-                             const Vector& localH, Vector& localL,
-                             const Array<real_t>& Bt, const Vector& mask)
+void Restriction2D(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);
@@ -2519,80 +2448,9 @@ void ProlongationTranspose2D(const int NE, const int D1D, const int Q1D,
       }
    });
 }
-
-template <int DLO, int DHI>
-static void SmemProlongationTranspose3D(
-   const int NE, const Vector& localH, Vector& localL,
-   const Array<real_t>& bt, const Vector& mask)
-{
-   MFEM_PERF_FUNCTION;
-
-   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)
+void Restriction3D(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);
@@ -2660,15 +2518,11 @@ void ProlongationTranspose3D(const int NE, const int D1D, const int Q1D,
       }
    });
 }
-
 } // namespace TransferKernels
 
 void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
                                                     Vector& y) const
 {
-   MFEM_PERF_FUNCTION;
-   using namespace TransferKernels;
-
    if (lFESpace.GetMesh()->GetNE() == 0)
    {
       return;
@@ -2677,25 +2531,11 @@ void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
    elem_restrict_lex_l->Mult(x, localL);
    if (dim == 2)
    {
-      Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
+      TransferKernels::Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
    }
    else if (dim == 3)
    {
-      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;
-      }
+      TransferKernels::Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask);
    }
    else
    {
@@ -2709,9 +2549,6 @@ void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
 void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
                                                              Vector& y) const
 {
-   MFEM_PERF_FUNCTION;
-   using namespace TransferKernels;
-
    if (lFESpace.GetMesh()->GetNE() == 0)
    {
       return;
@@ -2720,25 +2557,11 @@ void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
    elem_restrict_lex_h->Mult(x, localH);
    if (dim == 2)
    {
-      ProlongationTranspose2D(NE, D1D, Q1D, localH, localL, Bt, mask);
+      TransferKernels::Restriction2D(NE, D1D, Q1D, localH, localL, Bt, mask);
    }
    else if (dim == 3)
    {
-      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;
-      }
+      TransferKernels::Restriction3D(NE, D1D, Q1D, localH, localL, Bt, mask);
    }
    else
    {
@@ -2760,20 +2583,20 @@ TrueTransferOperator::TrueTransferOperator(const FiniteElementSpace& lFESpace_,
 
    P = lFESpace.GetProlongationMatrix();
    R = hFESpace.IsVariableOrder() ? hFESpace.GetHpRestrictionMatrix() :
-       hFESpace.GetRestrictionOperator();
+       hFESpace.GetRestrictionMatrix();
 
    // P and R can be both null
    // P can be null and R not null
    // If P is not null it is assumed that R is not null as well
    if (P) { MFEM_VERIFY(R, "Both P and R have to be not NULL") }
 
-   if (!IsIdentityProlongation(P))
+   if (P)
    {
       tmpL.SetSize(lFESpace_.GetVSize());
       tmpH.SetSize(hFESpace_.GetVSize());
    }
    // P can be null and R not null
-   else if (!IsIdentityProlongation(R))
+   else if (R)
    {
       tmpH.SetSize(hFESpace_.GetVSize());
    }
@@ -2786,14 +2609,13 @@ TrueTransferOperator::~TrueTransferOperator()
 
 void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
 {
-   MFEM_PERF_FUNCTION;
-   if (!IsIdentityProlongation(P))
+   if (P)
    {
       P->Mult(x, tmpL);
       localTransferOperator->Mult(tmpL, tmpH);
       R->Mult(tmpH, y);
    }
-   else if (!IsIdentityProlongation(R))
+   else if (R)
    {
       localTransferOperator->Mult(x, tmpH);
       R->Mult(tmpH, y);
@@ -2806,14 +2628,13 @@ void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
 
 void TrueTransferOperator::MultTranspose(const Vector& x, Vector& y) const
 {
-   MFEM_PERF_FUNCTION;
-   if (!IsIdentityProlongation(P))
+   if (P)
    {
       R->MultTranspose(x, tmpH);
       localTransferOperator->MultTranspose(tmpH, tmpL);
       P->MultTranspose(tmpL, y);
    }
-   else if (!IsIdentityProlongation(R))
+   else if (R)
    {
       R->MultTranspose(x, tmpH);
       localTransferOperator->MultTranspose(tmpH, y);

fem/transfer.hpp

@@ -621,6 +621,9 @@ 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;
@@ -639,7 +642,7 @@ private:
    const FiniteElementSpace& lFESpace;
    const FiniteElementSpace& hFESpace;
    const Operator * P = nullptr;
-   const Operator * R = nullptr;
+   const SparseMatrix * R = nullptr;
    TransferOperator* localTransferOperator;
    mutable Vector tmpL;
    mutable Vector tmpH;

general/annotation.hpp

@@ -14,98 +14,24 @@
 
 #include "../config/config.hpp"
 
-#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 <caliper/cali.h>
 #include <caliper/cali-manager.h>
-#endif
-
-namespace mfem
-{
-
-namespace internal
-{
-
-extern int annotation_sync_stream; // defined in globals.cpp
-extern int annotation_sync_mpi; // defined in globals.cpp
-
-#ifdef MFEM_USE_CALIPER
-
-inline void AnnotationSync()
-{
-   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
-{
-   ::cali::Function cali_func;
-
-   FunctionAnnotation(const char *fname)
-      : cali_func((AnnotationSync(), fname)) { }
-
-   ~FunctionAnnotation() { AnnotationSync(); }
-};
-
-struct ScopeAnnotation
-{
-   ::cali::ScopeAnnotation cali_scope;
-
-   ScopeAnnotation(const char *name)
-      : cali_scope((AnnotationSync(), 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_FUNCTION CALI_CXX_MARK_FUNCTION
+#define MFEM_PERF_BEGIN(s) CALI_MARK_BEGIN(s)
+#define MFEM_PERF_END(s) CALI_MARK_END(s)
 #define MFEM_PERF_SCOPE(name) \
-   mfem::internal::ScopeAnnotation \
-      MFEM_CONCAT(mfem_scope_annotation_,__LINE__)(name)
+  cali::Annotation::Guard cali_autogenerated_guard_name(cali::Annotation("function").begin(std::string(name).c_str()))
 
-#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))
 
-#else // #ifdef MFEM_USE_CALIPER
+#else
 
 #define MFEM_PERF_FUNCTION
-#define MFEM_PERF_BEGIN(s) ((void)(0))
+#define MFEM_PERF_BEGIN(s)
 #define MFEM_PERF_END(s)
 #define MFEM_PERF_SCOPE(name)
 
-#define MFEM_PERF_SYNC_STREAM(b)
-#define MFEM_PERF_SYNC_MPI(b)
-#define MFEM_PERF_SYNC(b)
+#endif
 
-#endif // #ifdef MFEM_USE_CALIPER
-
-#endif // MFEM_ANNOTATION_HPP
+#endif

general/array.hpp

@@ -211,6 +211,9 @@ public:
    /// Delete the first entry with value == 'el'.
    inline void DeleteFirst(const T &el);
 
+   /// Delete entries at @a indices, and resize.
+   inline void DeleteAt(const Array<int> &indices);
+
    /// Delete the whole array.
    inline void DeleteAll();
 
@@ -935,6 +938,30 @@ inline void Array<T>::DeleteFirst(const T &el)
    }
 }
 
+template <class T>
+inline void Array<T>::DeleteAt(const Array<int> &indices)
+{
+   // Make a copy of the indices, sorted.
+   Array<int> sorted_indices(indices);
+   sorted_indices.Sort();
+
+   int rm_count = 0;
+   for (int i = 0; i < size; i++)
+   {
+      if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
+      {
+         rm_count++;
+      }
+      else
+      {
+         data[i-rm_count] = data[i]; // shift data rm_count
+      }
+   }
+
+   // Resize to remove tail
+   size -= rm_count;
+}
+
 template <class T>
 inline void Array<T>::DeleteAll()
 {

general/communication.hpp

@@ -23,9 +23,13 @@
 #include <mpi.h>
 #include <cstdint>
 
-// 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
+// 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
 
 namespace mfem
 {

general/device.cpp

@@ -151,22 +151,6 @@ 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 // MFEM_ERROR_HPP
+#endif

general/globals.cpp

@@ -31,9 +31,6 @@ 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
 }
 
 void OutStream::Init()

general/mem_manager.hpp

@@ -657,8 +657,7 @@ 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,
@@ -742,7 +741,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,
@@ -982,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_((void*)h_ptr);
+      MemoryType h_ptr_mt = MemoryManager::GetHostMemoryType_(h_ptr);
       MFEM_VERIFY(h_mt == h_ptr_mt,
                   "h_mt = " << (int)h_mt << ", h_ptr_mt = " << (int)h_ptr_mt);
    }
@@ -990,8 +989,7 @@ 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_((void*)ptr, (void*)ptr, bytes, h_mt, own, false,
-                               flags);
+      MemoryManager::Register_(ptr, ptr, bytes, h_mt, own, false, flags);
    }
 }
 

general/scan.hpp

@@ -0,0 +1,176 @@
+// 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 MFEM_SCAN_HPP
+#define MFEM_SCAN_HPP
+
+#ifdef MFEM_USE_CUDA
+#include <cub/device/device_scan.cuh>
+#define MFEM_CUB_NAMESPACE cub
+#elif MFEM_USE_HIP
+#include <hipcub/device/device_scan.hpp>
+#define MFEM_CUB_NAMESPACE hipcub
+#endif
+
+#include <functional>
+#include <numeric>
+
+namespace mfem
+{
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
+/// std::plus<>{})
+template <class InputIt, class OutputIt>
+void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   Array<char> &workspace)
+{
+   // forward to InclusiveSum for potentially faster kernels
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                       workspace.Write(), bytes, d_in, d_out, num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                        nullptr, bytes, d_in, d_out, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
+                        workspace.Write(), bytes, d_in, d_out, num_items));
+      return;
+   }
+#endif
+   std::inclusive_scan(d_in, d_in + num_items, d_out);
+}
+
+/// Performs an inclusive scan of [d_in, d_in+num_items) -> [d_out,
+/// d_out+num_items). This call is potentially asynchronous on the device.
+/// @a d_in input start.
+/// @a d_out output start. Can perform in-place scans with d_out = d_in
+/// @a workspace temporary workspace used for device scans. TODO: replace with
+/// internal temporary workspace once that's added to the memory manager.
+/// @a scan_op binary scan functor. Must be associative. If only weakly
+/// associative (i.e. floating point addition) results are not deterministic. On
+/// device this must also be commutative.
+template <class InputIt, class OutputIt, class ScanOp>
+void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   Array<char> &workspace, ScanOp scan_op)
+{
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                        nullptr, bytes, d_in, d_out, scan_op, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
+      return;
+   }
+#endif
+   std::inclusive_scan(d_in, d_in + num_items, d_out, scan_op);
+}
+
+/// Performs an exclusive scan of [d_in, d_in+num_items) -> [d_out,
+/// d_out+num_items). This call is potentially asynchronous on the device.
+/// @a d_in input start.
+/// @a d_out output start. Can perform in-place scans with d_out = d_in
+/// @a workspace temporary workspace used for device scans. TODO: replace with
+/// internal temporary workspace once that's added to the memory manager.
+/// @a scan_op binary scan functor. Must be associative. If only weakly
+/// associative (i.e. floating point addition) results are not deterministic. On
+/// device this must also be commutative.
+template <class InputIt, class OutputIt, class T, class ScanOp>
+void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   T init_value, Array<char> &workspace, ScanOp scan_op)
+{
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      size_t bytes = workspace.Size();
+      if (bytes)
+      {
+         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
+#if defined(MFEM_USE_CUDA)
+         if (err == cudaSuccess)
+         {
+            return;
+         }
+#elif defined(MFEM_USE_HIP)
+         if (err == hipSuccess)
+         {
+            return;
+         }
+#endif
+      }
+      // try allocating a larger buffer
+      bytes = 0;
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                        nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
+      workspace.SetSize(bytes);
+      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
+      return;
+   }
+#endif
+   std::exclusive_scan(d_in, d_in + num_items, d_out, init_value, scan_op);
+}
+
+/// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
+/// workspace, std::plus<>{})
+template <class InputIt, class OutputIt, class T>
+void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
+                   T init_value, Array<char> &workspace)
+{
+   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, workspace,
+                 std::plus<> {});
+}
+
+} // namespace mfem
+
+#undef MFEM_CUB_NAMESPACE
+
+#endif

linalg/batched/gpu_blas.cpp

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

linalg/hypre.cpp

@@ -1868,8 +1868,6 @@ HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
 
 void HypreParMatrix::Mult(real_t a, const Vector &x, real_t b, Vector &y) const
 {
-   MFEM_PERF_FUNCTION;
-
    MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
                << ", expected size = " << Width());
    MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
@@ -1928,8 +1926,6 @@ void HypreParMatrix::Mult(real_t a, const Vector &x, real_t b, Vector &y) const
 void HypreParMatrix::MultTranspose(real_t a, const Vector &x,
                                    real_t b, Vector &y) const
 {
-   MFEM_PERF_FUNCTION;
-
    MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
                << ", expected size = " << Height());
    MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
@@ -4095,8 +4091,6 @@ 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);
@@ -4122,8 +4116,6 @@ 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,19 +50,17 @@ void Operator::InitTVectors(const Operator *Po, const Operator *Ri,
 
 void Operator::AddMult(const Vector &x, Vector &y, const real_t a) const
 {
-   z_am.SetSize(y.Size());
-   z_am.UseDevice(true);
-   Mult(x, z_am);
-   y.Add(a, z_am);
+   mfem::Vector z(y.Size());
+   Mult(x, z);
+   y.Add(a, z);
 }
 
 void Operator::AddMultTranspose(const Vector &x, Vector &y,
                                 const real_t a) const
 {
-   z_am.SetSize(y.Size());
-   z_am.UseDevice(true);
-   MultTranspose(x, z_am);
-   y.Add(a, z_am);
+   mfem::Vector z(y.Size());
+   MultTranspose(x, z);
+   y.Add(a, z);
 }
 
 void Operator::ArrayMult(const Array<const Vector *> &X,
@@ -588,8 +586,6 @@ void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
 void ConstrainedOperator::ConstrainedMult(const Vector &x, Vector &y,
                                           const bool transpose) const
 {
-   MFEM_PERF_FUNCTION;
-
    const int csz = constraint_list.Size();
    if (csz == 0)
    {
@@ -789,8 +785,6 @@ void RectangularConstrainedOperator::EliminateRHS(const Vector &x,
 
 void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
 {
-   MFEM_PERF_FUNCTION;
-
    const int trial_csz = trial_constraints.Size();
    const int test_csz = test_constraints.Size();
    if (trial_csz == 0)
@@ -826,8 +820,6 @@ void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
 void RectangularConstrainedOperator::MultTranspose(const Vector &x,
                                                    Vector &y) const
 {
-   MFEM_PERF_FUNCTION;
-
    const int trial_csz = trial_constraints.Size();
    const int test_csz = test_constraints.Size();
    if (test_csz == 0)

linalg/operator.hpp

@@ -13,7 +13,6 @@
 #define MFEM_OPERATOR
 
 #include "vector.hpp"
-#include "../general/annotation.hpp"
 
 namespace mfem
 {
@@ -24,13 +23,6 @@ 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.
@@ -826,12 +818,10 @@ public:
    explicit IdentityOperator(int n) : Operator(n) { }
 
    /// Operator application
-   void Mult(const Vector &x, Vector &y) const override
-   { MFEM_PERF_FUNCTION; y = x; }
+   void Mult(const Vector &x, Vector &y) const override { y = x; }
 
    /// Application of the transpose
-   void MultTranspose(const Vector &x, Vector &y) const override
-   { MFEM_PERF_FUNCTION; y = x; }
+   void MultTranspose(const Vector &x, Vector &y) const override { y = x; }
 };
 
 /// Returns true if P is the identity prolongation, i.e. if it is either NULL or

linalg/solvers.cpp

@@ -55,8 +55,6 @@ IterativeSolver::IterativeSolver(MPI_Comm comm_)
 
 real_t IterativeSolver::Dot(const Vector &x, const Vector &y) const
 {
-   MFEM_PERF_FUNCTION;
-
 #ifndef MFEM_USE_MPI
    return (x * y);
 #else
@@ -316,29 +314,25 @@ 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
-      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]);
-      });
+      subtract(x, residual, residual); // r = x - A y
    }
    else
    {
-      auto Y = y.Write();
-      // y = D^{-1} x
-      mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
-      {
-         Y[i] = DI[i] * X[i];
-      });
+      residual = x;
+      y.UseDevice(true);
+      y = 0.0;
    }
+   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_,
@@ -354,8 +348,7 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
    diag(d),
    coeffs(order),
    ess_tdof_list(ess_tdofs),
-   residual(order > 1 ? N : 0),
-   z(order > 1 ? N : 0),
+   residual(N),
    oper(&oper_) { Setup(); }
 
 #ifdef MFEM_USE_MPI
@@ -375,15 +368,14 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
                                                      real_t power_tolerance,
                                                      int power_seed)
 #endif
-   : Solver((MFEM_PERF_BEGIN(_MFEM_FUNC_NAME), d.Size())),
+   : Solver(d.Size()),
      order(order_),
      N(d.Size()),
      dinv(N),
      diag(d),
      coeffs(order),
      ess_tdof_list(ess_tdofs),
-     residual(order > 1 ? N : 0),
-     z(order > 1 ? N : 0),
+     residual(N),
      oper(&oper_)
 {
    OperatorJacobiSmoother invDiagOperator(diag, ess_tdofs, 1.0);
@@ -402,7 +394,6 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
                                                             power_seed);
 
    Setup();
-   MFEM_PERF_END(_MFEM_FUNC_NAME);
 }
 
 OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator* oper_,
@@ -431,7 +422,7 @@ void OperatorChebyshevSmoother::Setup()
 {
    // Invert diagonal
    residual.UseDevice(true);
-   z.UseDevice(true);
+   helperVector.UseDevice(true);
    auto D = diag.Read();
    auto X = dinv.Write();
    mfem::forall(N, [=] MFEM_HOST_DEVICE (int i) { X[i] = 1.0 / D[i]; });
@@ -441,20 +432,6 @@ 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.
@@ -524,8 +501,6 @@ void OperatorChebyshevSmoother::SetOrder(int new_order)
 
 void OperatorChebyshevSmoother::Mult(const Vector& x, Vector &y) const
 {
-   MFEM_PERF_FUNCTION;
-
    if (iterative_mode)
    {
       MFEM_ABORT("Chebyshev smoother not implemented for iterative mode");
@@ -536,55 +511,32 @@ void OperatorChebyshevSmoother::Mult(const Vector& x, Vector &y) const
       MFEM_ABORT("Chebyshev smoother requires operator");
    }
 
-   // 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)
-   {
-      mfem::forall(N, [=] MFEM_HOST_DEVICE (int i)
-      {
-         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]);
-      });
-   }
+   residual = x;
+   helperVector.SetSize(x.Size());
+   helperVector.UseDevice(true);
 
-   for (int k = 1; k < order; ++k)
-   {
-      // Apply: z = A r
-      oper->Mult(residual, z);
+   y.UseDevice(true);
+   y = 0.0;
 
-      // 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();
+   for (int k = 0; k < order; ++k)
+   {
+      // Apply
+      if (k > 0)
+      {
+         oper->Mult(residual, helperVector);
+         residual = helperVector;
+      }
+
+      // 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]; });
+
+      // Add weighted contribution to y
       auto Y = y.ReadWrite();
-      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];
-         });
-      }
+      auto C = coeffs.Read();
+      mfem::forall(n, [=] MFEM_HOST_DEVICE (int i) { Y[i] += C[k] * R[i]; });
    }
 }
 
@@ -3261,7 +3213,7 @@ void ResidualBCMonitor::MonitorResidual(
    MPI_Comm comm = iter_solver->GetComm();
    if (comm != MPI_COMM_NULL)
    {
-      real_t glob_bc_norm_squared = 0.0;
+      double 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 override;
+   void Mult(const Vector &x, Vector &y) const;
 
    /** @brief Approach the solution of the transposed linear system by applying
        Jacobi smoothing. */
-   void MultTranspose(const Vector &x, Vector &y) const override { Mult(x, y); }
+   void MultTranspose(const Vector &x, Vector &y) const { 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) override;
+   void SetOperator(const Operator &op);
 
 private:
    Vector dinv;
@@ -481,22 +481,21 @@ public:
 
    /** @brief Approach the solution of the linear system by applying Chebyshev
        smoothing. */
-   void Mult(const Vector &x, Vector &y) const override;
+   void Mult(const Vector &x, Vector &y) const;
 
    /** @brief Approach the solution of the transposed linear system by applying
        Chebyshev smoothing. */
-   void MultTranspose(const Vector &x, Vector &y) const override { Mult(x, y); }
+   void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
 
-   void SetOperator(const Operator &op_) override
+   void SetOperator(const Operator &op_)
    {
       oper = &op_;
    }
 
    void Setup();
-   void SetOrder(int new_order);
 
 private:
-   int order;
+   const int order;
    real_t max_eig_estimate;
    const int N;
    Vector dinv;
@@ -504,7 +503,7 @@ private:
    Array<real_t> coeffs;
    const Array<int>& ess_tdof_list;
    mutable Vector residual;
-   mutable Vector z;
+   mutable Vector helperVector;
    const Operator* oper;
 };
 

linalg/sparsemat.cpp

@@ -764,8 +764,6 @@ void SparseMatrix::Mult(const Vector &x, Vector &y) const
 
 void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
 {
-   MFEM_PERF_FUNCTION;
-
    MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
                << ") must match matrix width (" << width << ")");
    MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
@@ -966,8 +964,6 @@ void SparseMatrix::MultTranspose(const Vector &x, Vector &y) const
 void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
                                     const real_t a) const
 {
-   MFEM_PERF_FUNCTION;
-
    MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
                << ") must match matrix height (" << height << ")");
    MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()

linalg/vector.cpp

@@ -14,6 +14,7 @@
 #include "../general/forall.hpp"
 #include "../general/reducers.hpp"
 #include "../general/hash.hpp"
+#include "../general/scan.hpp"
 #include "vector.hpp"
 
 #ifdef MFEM_USE_OPENMP
@@ -205,16 +206,14 @@ Vector &Vector::operator=(const Vector &v)
    data.CopyFrom(v.data, v.Size());
    UseDevice(v.UseDevice());
 #else
+   SetSize(v.Size());
    const bool vuse = v.UseDevice();
    const bool use_dev = UseDevice() || vuse;
-   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
-   SetSize(v.Size());
    v.UseDevice(use_dev);
    // keep 'data' where it is, unless 'use_dev' is true
    if (use_dev) { Write(); }
    data.CopyFrom(v.data, v.Size());
    v.UseDevice(vuse);
-   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
 #endif
    return *this;
 }
@@ -229,11 +228,9 @@ Vector &Vector::operator=(Vector &&v)
 Vector &Vector::operator=(real_t value)
 {
    const bool use_dev = UseDevice();
-   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
    const int N = size;
    auto y = Write(use_dev);
    mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = value; });
-   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
    return *this;
 }
 
@@ -294,12 +291,10 @@ Vector &Vector::operator-=(const Vector &v)
    MFEM_ASSERT(size == v.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
    const int N = size;
    const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
    mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= x[i]; });
-   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
    return *this;
 }
 
@@ -317,12 +312,10 @@ Vector &Vector::operator+=(const Vector &v)
    MFEM_ASSERT(size == v.size, "incompatible Vectors!");
 
    const bool use_dev = UseDevice() || v.UseDevice();
-   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
    const int N = size;
    const auto x = v.Read(use_dev);
    auto y = ReadWrite(use_dev);
    mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += x[i]; });
-   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
    return *this;
 }
 
@@ -334,11 +327,9 @@ Vector &Vector::Add(const real_t a, const Vector &Va)
    {
       const int N = size;
       const bool use_dev = UseDevice() || Va.UseDevice();
-      if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
       const auto x = Va.Read(use_dev);
       auto y = ReadWrite(use_dev);
       mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += a * x[i]; });
-      if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
    }
    return *this;
 }
@@ -455,7 +446,6 @@ void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v)
    {
 #if !defined(MFEM_USE_LEGACY_OPENMP)
       const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
-      if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
       const int N = v.size;
       // Note: get read access first, in case v is the same as v1/v2.
       const auto d_x = v1.Read(use_dev);
@@ -465,7 +455,6 @@ void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v)
       {
          d_z[i] = d_x[i] + alpha * d_y[i];
       });
-      if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
 #else
       const real_t *v1p = v1.data, *v2p = v2.data;
       real_t *vp = v.data;
@@ -581,7 +570,6 @@ void subtract(const Vector &x, const Vector &y, Vector &z)
 
 #if !defined(MFEM_USE_LEGACY_OPENMP)
    const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
-   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
    const int N = x.size;
    // Note: get read access first, in case z is the same as x/y.
    const auto xd = x.Read(use_dev);
@@ -591,7 +579,6 @@ void subtract(const Vector &x, const Vector &y, Vector &z)
    {
       zd[i] = xd[i] - yd[i];
    });
-   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
 #else
    const real_t *xp = x.data;
    const real_t *yp = y.data;
@@ -1266,4 +1253,44 @@ real_t Vector::Sum() const
    return res;
 }
 
+void Vector::DeleteAt(const Array<int> &indices)
+{
+   const bool use_dev = UseDevice();
+
+   Array<int> flag(size);
+   const auto d_flag = flag.Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      d_flag[i] = true;
+   });
+   const auto d_indices = indices.Read(use_dev);
+   mfem::forall_switch(use_dev, indices.Size(), [=] MFEM_HOST_DEVICE (int i)
+   {
+      d_flag[d_indices[i]] = false;
+   });
+
+   Array<int> out_idx(size);
+   auto d_out_idx = out_idx.Write(use_dev);
+   Array<char> workspace;
+   // Perform inclusive scan so that the last entry is the new size.
+   InclusiveScan(use_dev, d_flag, d_out_idx, size, workspace);
+
+   Vector copy(*this);
+   auto d_in = copy.Read(use_dev);
+   auto d_out = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      if (d_flag[i])
+      {
+         // Transform inclusive scan to exclusive by shifting.
+         const int j = (i > 0) ? d_out_idx[i - 1] : 0;
+         d_out[j] = d_in[i];
+      }
+   });
+
+   // Get the new size of the vector. Copy only the last entry.
+   Memory<int> submem(out_idx.GetMemory(), out_idx.Size() - 1, 1);
+   size = submem.Read(MemoryClass::HOST, 1)[0];
+}
+
 } // namespace mfem

linalg/vector.hpp

@@ -171,6 +171,12 @@ public:
    /// Resize the vector to size @a s using the MemoryType of @a v.
    void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
 
+   /// Update \ref Capacity() to @a res (if less than current), keeping existing entries.
+   void Reserve(int res);
+
+   /// Delete entries at @a indices and resize vector accordingly.
+   void DeleteAt(const Array<int> &indices);
+
    /// Set the Vector data.
    /// @warning This method should be called only when OwnsData() is false.
    void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
@@ -621,6 +627,18 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
+inline void Vector::Reserve(int res)
+{
+   if (res > Capacity())
+   {
+      Memory<real_t> p(res, data.GetMemoryType());
+      p.CopyFrom(data, size);
+      p.UseDevice(data.UseDevice());
+      data.Delete();
+      data = p;
+   }
+}
+
 inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
                                      bool own_mem)
 {

makefile

@@ -125,8 +125,7 @@ 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 \
- benchmarks/ceed-solver-bps
+ hooke multidomain dpg hdiv-linear-solver spde diag-smoothers
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \

mesh/mesh.cpp

@@ -979,9 +979,48 @@ const Array<int>& Mesh::GetElementAttributes() const
    return elem_attrs_cache;
 }
 
+void Mesh::ComputeFaceInfo(FaceType ftype) const
+{
+   auto &fidcs = face_indices[static_cast<int>(ftype)];
+   auto &ifidcs = inv_face_indices[static_cast<int>(ftype)];
+   fidcs.SetSize(GetNFbyType(ftype));
+   fidcs.HostWrite();
+   ifidcs.reserve(fidcs.Size());
+   int f_idx = 0;
+   for (int i = 0; i < GetNumFacesWithGhost(); ++i)
+   {
+      const FaceInformation face = GetFaceInformation(i);
+      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(ftype))
+      {
+         continue;
+      }
+      fidcs[f_idx] = i;
+      ifidcs[i] = f_idx;
+      ++f_idx;
+   }
+}
+
+const Array<int> &Mesh::GetFaceIndices(FaceType ftype) const
+{
+   if (face_indices[static_cast<int>(ftype)].Size() == 0)
+   {
+      ComputeFaceInfo(ftype);
+   }
+   return face_indices[static_cast<int>(ftype)];
+}
+
+const std::unordered_map<int, int> &
+Mesh::GetInvFaceIndices(FaceType ftype) const
+{
+   if (inv_face_indices[static_cast<int>(ftype)].empty())
+   {
+      ComputeFaceInfo(ftype);
+   }
+   return inv_face_indices[static_cast<int>(ftype)];
+}
+
 void Mesh::DeleteGeometricFactors()
 {
-   MFEM_PERF_FUNCTION;
    for (int i = 0; i < geom_factors.Size(); i++)
    {
       delete geom_factors[i];
@@ -1867,6 +1906,11 @@ void Mesh::Destroy()
    bdr_face_attrs_cache.DeleteAll();
    attributes.DeleteAll();
    bdr_attributes.DeleteAll();
+
+   face_indices[0].DeleteAll();
+   face_indices[1].DeleteAll();
+   inv_face_indices[0] = std::unordered_map<int, int>();
+   inv_face_indices[1] = std::unordered_map<int, int>();
 }
 
 void Mesh::ResetLazyData()
@@ -6580,7 +6624,6 @@ void XYZ_VectorFunction(const Vector &p, Vector &v)
 
 void Mesh::GetNodes(GridFunction &nodes) const
 {
-   MFEM_PERF_FUNCTION;
    if (Nodes == NULL || Nodes->FESpace() != nodes.FESpace())
    {
       const int newSpaceDim = nodes.FESpace()->GetVDim();
@@ -6601,7 +6644,6 @@ void Mesh::SetNodalFESpace(FiniteElementSpace *nfes)
 
 void Mesh::EnsureNodes()
 {
-   MFEM_PERF_FUNCTION;
    if (Nodes)
    {
       const FiniteElementCollection *fec = GetNodalFESpace()->FEColl();
@@ -6654,7 +6696,6 @@ const FiniteElementSpace *Mesh::GetNodalFESpace() const
 
 void Mesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
 {
-   MFEM_PERF_FUNCTION;
    if (order <= 0)
    {
       delete Nodes;
@@ -8141,6 +8182,12 @@ void Mesh::GenerateFaces()
       FreeElement(f);
    }
 
+   // delete caches
+   face_indices[0].SetSize(0);
+   face_indices[1].SetSize(0);
+   inv_face_indices[0].clear();
+   inv_face_indices[1].clear();
+
    // (re)generate the interior faces and the info for them
    faces.SetSize(nfaces);
    faces_info.SetSize(nfaces);
@@ -10940,6 +10987,11 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
    // copy attribute caches
    mfem::Swap(elem_attrs_cache, other.elem_attrs_cache);
    mfem::Swap(bdr_face_attrs_cache, other.bdr_face_attrs_cache);
+
+   mfem::Swap(face_indices[0], other.face_indices[0]);
+   mfem::Swap(face_indices[1], other.face_indices[1]);
+   inv_face_indices[0].swap(other.inv_face_indices[0]);
+   inv_face_indices[1].swap(other.inv_face_indices[1]);
 }
 
 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,
@@ -14739,7 +14791,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();
    const int dim  = fe->GetDim();

mesh/mesh.hpp

@@ -278,6 +278,13 @@ protected:
 
    // used during NC mesh initialization only
    Array<Triple<int, int, int> > tmp_vertex_parents;
+   /// cache for FaceIndices(ftype)
+   mutable Array<int> face_indices[2];
+   /// cache for FaceIndices(ftype)
+   mutable std::unordered_map<int, int> inv_face_indices[2];
+
+   /// compute face_indices[ftype] and inv_face_indices[type]
+   void ComputeFaceInfo(FaceType ftype) const;
 
 public:
    typedef Geometry::Constants<Geometry::SEGMENT>     seg_t;
@@ -312,6 +319,11 @@ public:
    // (true) is set in mesh_readers.cpp.
    static bool remove_unused_vertices;
 
+   /// Map from boundary or interior face indices to mesh face indices.
+   const Array<int>& GetFaceIndices(FaceType ftype) const;
+   /// Inverse of the map FaceIndices(ftype)
+   const std::unordered_map<int, int>& GetInvFaceIndices(FaceType ftype) const;
+
 protected:
    Operation last_operation;
 

mesh/pmesh.cpp

@@ -2016,7 +2016,6 @@ void ParMesh::DeleteFaceNbrData()
 
 void ParMesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
 {
-   MFEM_PERF_FUNCTION;
    DeleteFaceNbrData();
    space_dim = (space_dim == -1) ? spaceDim : space_dim;
    FiniteElementCollection* nfec;

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

@@ -1,156 +0,0 @@
-// 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

@@ -1,77 +0,0 @@
-# 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

@@ -1,129 +0,0 @@
-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

@@ -1,241 +0,0 @@
-// 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

@@ -1,100 +0,0 @@
-// 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

@@ -1,364 +0,0 @@
-// 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
-
-using BPSFunctionType = real_t(*)(int n, const real_t *xyz);
-
-template <BPSFunctionType F>
-void ProjectBPSFunction(int n, QuadratureFunction &qf)
-{
-   MFEM_PERF_FUNCTION;
-   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);
-   });
-}
-
-MFEM_HOST_DEVICE inline
-real_t sol_1d(const int n, const real_t *xyz)
-{
-   return w(n, xyz[0]);
-}
-
-MFEM_HOST_DEVICE inline
-real_t sol_2d(const int n, const real_t *xyz)
-{
-   return w(n, xyz[0])*w(n, xyz[1]);
-}
-
-MFEM_HOST_DEVICE inline
-real_t sol_3d(const int n, const real_t *xyz)
-{
-   return w(n, xyz[0])*w(n, xyz[1])*w(n, xyz[2]);
-}
-
-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]);
-      }
-   }
-
-   void Project(QuadratureFunction &qf) override
-   {
-      switch (dim)
-      {
-         case 1: ProjectBPSFunction<sol_1d>(n, qf); break;
-         case 2: ProjectBPSFunction<sol_2d>(n, qf); break;
-         case 3: ProjectBPSFunction<sol_3d>(n, qf); break;
-         default: MFEM_ABORT("Unsupported dimension.");
-      }
-   }
-};
-
-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]);
-}
-
-void ProjectRHS(int n, QuadratureFunction &qf)
-{
-   const int dim = qf.GetSpace()->GetMesh()->Dimension();
-   switch (dim)
-   {
-      case 1: ProjectBPSFunction<rhs_1d>(n, qf); break;
-      case 2: ProjectBPSFunction<rhs_2d>(n, qf); break;
-      case 3: ProjectBPSFunction<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

@@ -1,125 +0,0 @@
-bsep="============================================================"
-ssep="----------------------------------------"
-# Enable GPU-aware MPI:
-# gpu_aware_mpi_env_cmd="env MPICH_GPU_SUPPORT_ENABLED=1"
-# gpu_aware_mpi="-g"
-# number of MPI ranks, number of ranks per node, number of nodes:
-np=1
-nrnode=4
-((nnodes = (np+nrnode-1)/nrnode))
-# dev="-d gpu ${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 -x -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 " *** np = ${np}, 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} -nrn ${nrnode}"
-         printf " -ey ${eps} -mg \"${mg}\" -cs 1 ${p_mg_opts}"
-         printf " -nx ${nx} -rp ${rp}\n"
-         echo "${ssep}"
-         $mpirun_np "${np}" ./solver-bp ${dev} -nrn ${nrnode} \
-            -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} -nrn ${nrnode}"
-         printf " -ey ${eps} -mg \"${p}\" -cs 2 ${lor_opts}"
-         printf " -nx ${nx} -rp ${rp}\n"
-         echo "${ssep}"
-         $mpirun_np "${np}" ./solver-bp ${dev} -nrn ${nrnode} \
-            -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} -nrn ${nrnode}"
-         printf " -ey ${eps} -mg \"${p}\" -cs 0 -nx ${nx} -rp ${rp}\n"
-         echo "${ssep}"
-         $mpirun_np "${np}" ./solver-bp ${dev} -nrn ${nrnode} \
-            -ey ${eps} -mg "${p}" -cs 0 -nx "${nx}" -rp "${rp}"
-      fi
-   done
-done

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

@@ -1,845 +0,0 @@
-// 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_PERF_FUNCTION;
-
-      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>();
-   DiffusionIntegrator::AddSpecialization<3,7,7>();
-
-   Mpi::Init(argc, argv);
-   Hypre::Init();
-
-   const char *device_config = "cpu";
-   int nrnode = 4; // number of ranks per node, used for CSV output only
-   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(&nrnode, "-nrn", "--num-ranks-per-node",
-                  "Number of ranks per compute node. Used to compute the number"
-                  " of nodes written in CSV output.");
-   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
-   MFEM_PERF_BEGIN("CreateKershawMesh");
-   ParMesh mesh_coarse = CreateKershawMesh(nx, ny, nz, epsy, epsz);
-   MFEM_PERF_END("CreateKershawMesh");
-   const int dim = mesh_coarse.Dimension();
-   for (int i=0; i<ref_par; ++i)
-   {
-      MFEM_PERF_SCOPE("Mesh UniformRefinement");
-      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
-
-   MFEM_PERF_BEGIN("Setup [hierarchy]");
-   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());
-      }
-   }
-   MFEM_PERF_END("Setup [hierarchy]");
-
-   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();
-   MFEM_PERF_BEGIN("ParMesh PrintInfo");
-   mesh.PrintInfo(cout);
-   MFEM_PERF_END("ParMesh PrintInfo");
-   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; }
-   MFEM_PERF_BEGIN("Setup [RHS]");
-   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();
-   MFEM_PERF_END("Setup [RHS]");
-   if (Mpi::Root()) { cout << "Assembling right-hand side... Done." << endl; }
-
-   // Free device memory: the geometric facros computed so far are:
-   // * the coordinates, for the rhs coefficient evaluation, and
-   // * the detJ, for the DomainLFIntegrator.
-   // These are no-longer needed (?), so we can free the memory.
-   mesh.DeleteGeometricFactors();
-
-   // 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);
-   MFEM_PERF_BEGIN("Setup [DiffusionMultigrid]");
-   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 double t_setup = tic_toc.RealTime();
-   MFEM_PERF_END("Setup [DiffusionMultigrid]");
-
-   ParGridFunction x(&fes);
-   x = 0.0;
-
-   OperatorPtr A;
-   Vector X, B;
-   MFEM_PERF_BEGIN("Setup [MG.FormFineLinearSystem]");
-   MG.FormFineLinearSystem(x, b, A, X, B);
-   MFEM_PERF_END("Setup [MG.FormFineLinearSystem]");
-
-   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(3);
-   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"; }
-   MFEM_PERF_BEGIN("Warm-up");
-   cg.SetMaxIter(2);
-   {
-      Vector X_save(X);
-      cg.Mult(B, X);
-      X = X_save;
-   }
-   MFEM_PERF_END("Warm-up");
-   if (coarse_solver.inner_sli &&
-       ((coarse_solver.type == SolverConfig::FA_HYPRE /* && order > 1 */) ||
-        coarse_solver.type == SolverConfig::LOR_HYPRE))
-   {
-      MFEM_PERF_SCOPE("Auto-tuning");
-      // 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++)
-         {
-            MFEM_PERF_SCOPE(("Timing [" + to_string(sli_it) + "," +
-                             to_string(cheby_order) + "]").c_str());
-            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);
-   MFEM_PERF_BEGIN("Final CG Solve");
-   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();
-   MFEM_PERF_END("Final CG Solve");
-
-   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);
-
-   MFEM_PERF_BEGIN("Compute L2 Error");
-   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);
-   MFEM_PERF_END("Compute L2 Error");
-   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() + (nrnode-1))/nrnode; // 4
-      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)
-{
-#if defined(HYPRE_WITH_GPU_AWARE_MPI) || defined(HYPRE_USING_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: "
-           << (HypreUsingGPU() ? "enabled" : "disabled");
-      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/shifted/distance.cpp

@@ -384,10 +384,8 @@ int main(int argc, char *argv[])
    dacol.Save();
 
    ConstantCoefficient zero(0.0);
-   Vector zero_vec(dim); zero_vec = 0_r;
-   VectorConstantCoefficient vzero(zero_vec);
    const real_t s_norm  = distance_s.ComputeL2Error(zero),
-                v_norm  = distance_v.ComputeL2Error(vzero);
+                v_norm  = distance_v.ComputeL2Error(zero);
    if (myid == 0)
    {
       cout << fixed << setprecision(10) << "Norms: "

miniapps/solvers/block-solvers.cpp

@@ -365,7 +365,6 @@ 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,7 +32,6 @@ 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,7 +57,6 @@ 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
-   mutable MINRESSolver solver_;
+   MINRESSolver solver_;
    Array<int> ess_zero_dofs_;
 public:
    BDPMinresSolver(const HypreParMatrix& M,

miniapps/solvers/div_free_solver.cpp

@@ -84,6 +84,7 @@ 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_);
@@ -173,9 +174,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);
-   HypreParMatrix *res =
-      data_.C[level_+1].As<HypreParMatrix>()->EliminateCols(ess_hcurl_tdof);
-   delete res;
+   data_.Ae[level_+1].reset(
+      data_.C[level_+1].As<HypreParMatrix>()
+      ->EliminateCols(ess_hcurl_tdof));
 
    ++level_;
 

miniapps/solvers/div_free_solver.hpp

@@ -36,6 +36,7 @@ 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
@@ -45,6 +46,7 @@ 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;
 };
 

tests/unit/CMakeLists.txt

@@ -22,6 +22,7 @@ set(UNIT_TESTS_SRCS
   dfem/test_mass.cpp
   general/test_array.cpp
   general/test_reduction.cpp
+  general/test_scan.cpp
   general/test_arrays_by_name.cpp
   general/test_error.cpp
   general/test_mem.cpp

tests/unit/general/test_array.cpp

@@ -124,3 +124,19 @@ TEST_CASE("Array stl-interactions", "[Array]")
       CHECK(x[i] == y[i]);
    }
 }
+
+TEST_CASE("Array delete at indices", "[Array]")
+{
+   Array<int>       test({0,1,2,3,4,5,6,7,8});
+   Array<int> rm_indices({0,    3,4,  6,  8});
+   Array<int>     result({  1,2,    5,  7  });
+
+   test.DeleteAt(rm_indices);
+
+   REQUIRE(test.Size() == result.Size());
+
+   for (int i = 0; i < test.Size(); i++)
+   {
+      CHECK(test[i] == result[i]);
+   }
+}

tests/unit/general/test_scan.cpp

@@ -0,0 +1,102 @@
+// 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 <algorithm>
+#include <limits>
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+// must be included after mfem.hpp
+#include "general/scan.hpp"
+
+using namespace mfem;
+
+TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
+{
+   Array<char> workspace;
+   Array<int> a(10);
+
+   for (int use_dev = 0; use_dev < 2; ++use_dev)
+   {
+      CAPTURE(use_dev);
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i;
+      }
+      auto dptr = a.ReadWrite(use_dev);
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace);
+      a.HostRead();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         int expected = (i + 1) * i / 2;
+         CAPTURE(i);
+         REQUIRE(a[i] == expected);
+      }
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i + 1;
+      }
+      a.ReadWrite(use_dev);
+      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace, std::multiplies<> {});
+      a.HostRead();
+      int expected = 1;
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         expected *= i + 1;
+         CAPTURE(i);
+         REQUIRE(a[i] == expected);
+      }
+   }
+}
+
+TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
+{
+   Array<char> workspace;
+   Array<int> a(10);
+
+   for (int use_dev = 0; use_dev < 2; ++use_dev)
+   {
+      CAPTURE(use_dev);
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i;
+      }
+      auto dptr = a.ReadWrite(use_dev);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace);
+      a.HostRead();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         int expected = (i + 1) * i / 2 - i + 5;
+         CAPTURE(i);
+         REQUIRE(a[i] == expected);
+      }
+      a.HostReadWrite();
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         a[i] = i + 1;
+      }
+      a.ReadWrite(use_dev);
+      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace,
+                    std::multiplies<> {});
+      a.HostRead();
+      int expected = 5;
+      for (int i = 0; i < a.Size(); ++i)
+      {
+         CAPTURE(i);
+         REQUIRE(a[i] == expected);
+         expected *= i + 1;
+      }
+   }
+}

tests/unit/linalg/test_chebyshev.cpp

@@ -14,47 +14,53 @@
 
 using namespace mfem;
 
-TEST_CASE("Chebyshev symmetry", "[OperatorChebyshevSmoother]")
+TEST_CASE("OperatorChebyshevSmoother", "[Chebyshev symmetry]")
 {
-   const int order = GENERATE(2, 3, 4);
-   const int cheb_order = GENERATE(2, 3);
+   for (int order = 2; order < 5; ++order)
+   {
+      const int cheb_order = 2;
 
-   Mesh mesh = Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
-   H1_FECollection fec(order, 3);
-   FiniteElementSpace fespace(&mesh, &fec);
+      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);
 
-   Array<int> ess_tdof_list;
-   fespace.GetBoundaryTrueDofs(ess_tdof_list);
+      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);
 
-   BilinearForm aform(&fespace);
-   aform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   aform.AddDomainIntegrator(new DiffusionIntegrator);
-   aform.Assemble();
+      Solver* smoother = new OperatorChebyshevSmoother(*opr, diag, ess_tdof_list,
+                                                       cheb_order);
 
-   OperatorPtr opr;
-   opr.SetType(Operator::ANY_TYPE);
-   aform.FormSystemMatrix(ess_tdof_list, opr);
+      int n = smoother->Width();
+      Vector left(n);
+      Vector right(n);
+      int seed = (int) time(0);
+      left.Randomize(seed);
+      right.Randomize(seed + 2);
 
-   Vector diag(fespace.GetTrueVSize());
-   aform.AssembleDiagonal(diag);
+      // 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;
 
-   OperatorChebyshevSmoother smoother(*opr, diag, ess_tdof_list, cheb_order);
+      double error = fabs(forward_val - transpose_val) / fabs(forward_val);
+      CAPTURE(order, error);
+      REQUIRE(error < 1.e-13);
 
-   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));
+      delete smoother;
+      delete fec;
+   }
 }

tests/unit/linalg/test_vector.cpp

@@ -247,3 +247,21 @@ TEST_CASE("Vector Sum", "[Vector],[GPU]")
 
    REQUIRE(sum_1 == MFEM_Approx(sum_2));
 }
+
+TEST_CASE("Vector delete at indices", "[Vector][GPU]")
+{
+   Vector           test({0,1,2,3,4,5,6,7,8});
+   Array<int> rm_indices({0,    3,4,  6,  8});
+   Vector         result({  1,2,    5,  7  });
+
+   test.UseDevice(true);
+   test.DeleteAt(rm_indices);
+
+   REQUIRE(test.Size() == result.Size());
+
+   test.HostReadWrite();
+   for (int i = 0; i < test.Size(); i++)
+   {
+      CHECK(test[i] == result[i]);
+   }
+}