minor

AI policy

.gitlab/jobs/dane.yml

@@ -85,3 +85,8 @@ opt_par_gcc_10_pumi:
   extends: .mfem_job_on_dane
   variables:
     SPEC: "%gcc@10.3.1 +pumi"
+
+opt_par_gcc_10_gslib:
+  extends: .mfem_job_on_dane
+  variables:
+    SPEC: "%gcc@10.3.1 +gslib"

.gitlab/jobs/matrix.yml

@@ -63,3 +63,8 @@ opt_mpi_cuda_hypre_cuda_gcc:
   extends: .mfem_job_on_matrix
   variables:
     SPEC: "%gcc@10.3.1 +mpi +cuda cuda_arch=90 ^hypre+cuda"
+
+opt_mpi_cuda_gcc_gslib:
+  extends: .mfem_job_on_matrix
+  variables:
+    SPEC: "%gcc@10.3.1 +mpi +cuda +gslib cuda_arch=90 ^hypre+cuda"

.gitlab/scripts/baseline

@@ -32,9 +32,9 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
 
 # run
 if [[ "${MACHINE_NAME}" == "dane" ]]; then
-  salloc --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  srun --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
-  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  srun --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
   echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
   exit 1

CHANGELOG

@@ -8,8 +8,10 @@
                                https://mfem.org
 
 
-Version 4.10 (development)
-==========================
+Version 4.9.1 (development)
+===========================
+
+- Policy for AI-assisted contribution added to CONTRIBUTING.md
 
 Discretization improvements
 ---------------------------
@@ -23,12 +25,6 @@ Discretization improvements
   Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
   2022.
 
-
-Version 4.9.1 (development)
-===========================
-
-Discretization improvements
----------------------------
 - Improved the gridfunction projection routines. Projections work for Scalar,
   Vector and VectorFE, also NURBS versions. Optionally different types of
   projections can be selected, default behaviour has not changed.

CONTRIBUTING.md

@@ -24,6 +24,14 @@ must be made under this license.
 Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
 in the MFEM community, you agree to abide by its rules.
 
+## AI Policy
+- Use of AI code generation in MFEM is allowed but must be disclosed, e.g. by
+  selecting the `AI-assisted` label on the PR.
+- By submitting a PR, the author acknowledges that they have reviewed and
+  understand the changes they are proposing.
+- PR authors are still responsible for correctness, licensing, and attribution
+  of all changes.
+
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull

examples/CMakeLists.txt

@@ -215,7 +215,7 @@ if (MFEM_ENABLE_TESTING)
     add_test(NAME ex1p_ceed_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:ex1p> "-no-vis" "-d ceed-cpu" "-pa" "-a"
+      $<TARGET_FILE:ex1p> "-no-vis" "-d" "ceed-cpu" "-pa" "-a"
       ${MPIEXEC_POSTFLAGS})
   endif()
 endif()

examples/amgx/makefile

@@ -64,7 +64,7 @@ PARALLEL_NAME := Parallel AMGX example
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not build)
 
-clean: clean-build
+clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)

examples/caliper/makefile

@@ -64,12 +64,12 @@ ex1p-test-par: ex1p
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-clean: clean-build clean-exec $(SUBDIRS_CLEAN)
+clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
+	@rm -f refined.mesh mesh.*
+	@rm -f sol.*

examples/ginkgo/makefile

@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh sol.gf
+	@rm -f refined.mesh sol.gf mesh.* sol.*

examples/makefile

@@ -71,6 +71,7 @@ endif
 
 SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
 SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
+SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(SUBDIRS))
 SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
 SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
@@ -87,8 +88,9 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
 all: $(EXAMPLES) $(SUBDIRS_ALL)
 
-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
+   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -107,6 +109,7 @@ endif
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)
+test-noclean: $(SUBDIRS_TEST_NOCLEAN)
 test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Parallel vs. serial runs

fem/bounds.cpp

@@ -41,9 +41,14 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
    tol = tol_i;
    lbound.SetSize(ncp, nb);
    ubound.SetSize(ncp, nb);
+   lbound_t.SetSize(nb, ncp);
+   ubound_t.SetSize(nb, ncp);
    nodes.SetSize(nb);
    weights.SetSize(nb);
    control_points.SetSize(ncp);
+   xhat.SetSize(nb);
+   what.SetSize(nb);
+   cphat.SetSize(ncp);
 
    auto scalenodes = [](const Vector &in, const real_t a, const real_t b) -> Vector
    {
@@ -90,6 +95,10 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
       MFEM_ABORT("Unsupported interval points. Use [0,1].\n");
    }
    control_points = scalenodes(control_points, 0.0, 1.0); // rescale to [0,1]
+   for (int i = 0; i < ncp; i++)
+   {
+      cphat(i) = 2.0*control_points(i) - 1.0;
+   }
 
    Poly_1D::Basis &basis1d(poly1d.GetBasis(nb-1, b_type));
 
@@ -145,6 +154,8 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
                lbound(j,i) = std::max(lbound(j,i),0_r);
             }
          }
+         lbound_t(i,j) = lbound(j,i);
+         ubound_t(i,j) = ubound(j,i);
       }
    }
 
@@ -176,6 +187,11 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
          nodes(i) = irule.IntPoint(i).x;
       }
    }
+   for (int i = 0; i < nb; i++)
+   {
+      xhat(i) = 2.0*nodes(i) - 1.0;
+      what(i) = 2.0*weights(i);
+   }
 
    if (b_type == 2)
    {
@@ -755,4 +771,4 @@ void PLBound::Print(std::ostream &outp) const
    ubound.Print(outp);
 }
 
-}
+}

fem/bounds.hpp

@@ -13,6 +13,7 @@
 #define MFEM_BOUNDS
 
 #include "../config/config.hpp"
+#include "../general/forall.hpp"
 #include "fespace.hpp"
 
 namespace mfem
@@ -60,7 +61,9 @@ private:
    bool proj = true; // Use linear projection to compute bounds.
    real_t tol = 0.0; // offset bounds to avoid round-off errors
    Vector nodes, weights, control_points;
+   Vector xhat, what, cphat;
    DenseMatrix lbound, ubound; // ncp x nb matrices with bounds of all bases
+   DenseMatrix lbound_t, ubound_t; // nb x ncp transposes for device kernel
    // Some auxillary storage for computing the bounds with Bernstein
    DenseMatrix basisMatNodes; // Bernstein bases at equispaced nodes
    DenseMatrix basisMatInt;   // Bernstein bases at GLL nodes
@@ -113,7 +116,10 @@ public:
     *  @details This projection increases the computational cost but results in
     *  tighter bounds.
     */
-   void SetProjectionFlagForBounding(bool proj_) { proj = proj_; }
+   void SetProjectionFlagForBounding(bool proj_)
+   {
+      proj = proj_;
+   }
 
    /** @brief Compute piecewise linear bounds for the lexicographically-ordered
     *  nodal coefficients in @a coeff in 1D/2D/3D.
@@ -137,9 +143,23 @@ public:
    /// Get number of control points used to compute the bounds.
    int GetNControlPoints() const { return ncp; }
 
+   /// Get the underlying 1D basis type.
+   int GetBasisType() const { return b_type; }
+
    /// Get 1D control point locations (lexicographic order) in [0,1].
    const Vector &GetControlPoints() const { return control_points; }
 
+   /** @brief Compute element-wise bounds from a lexicographic E-vector.
+    *
+    *  @details The expected layout of @a e_vec is `ND x VDIM x NE`, where
+    *  `ND = nb^rdim`, `VDIM = fes_vdim`, and `NE` is the number of elements.
+    *  The output layout matches GridFunction::GetElementBounds:
+    *  `NE x active_vdim`, with the element index varying fastest.
+    */
+   void GetElementBoundsKernel(const int rdim, const int fes_vdim,
+                               const Vector &e_vec, Vector &lower,
+                               Vector &upper, const int vdim = 0) const;
+
    /** @brief Get lower and upper bounding matrix (ncp^dim x nb^dim)
     *
     *  @details The matrices can be used to compute the bounds at control points
@@ -183,6 +203,600 @@ private:
               const int cp_type_i, const real_t tol_i);
 };
 
+namespace internal
+{
+
+struct PLBoundDeviceData
+{
+   int nb;
+   int ncp;
+   const real_t *xhat;
+   const real_t *what;
+   const real_t *cphat;
+   const real_t *lbound;
+   const real_t *ubound;
+};
+
+template<int T_NB = 0, bool T_PROJ = true>
+inline void GetElementBoundsKernel1D(const PLBoundDeviceData &data,
+                                     const int fes_vdim,
+                                     const int ne,
+                                     const Vector &e_vec,
+                                     Vector &lower,
+                                     Vector &upper,
+                                     const int comp0,
+                                     const int ncomp)
+{
+   constexpr int GENERIC_MAX_ND = 32;
+   constexpr int MAX_ND = T_NB ? T_NB : GENERIC_MAX_ND;
+   constexpr int BLOCK_X = 2*MAX_ND;
+
+   const int nd = T_NB ? T_NB : data.nb;
+   MFEM_VERIFY(nd <= MAX_ND,
+               "Device element bounds kernel supports up to 32 "
+               "1D degrees of freedom.");
+
+   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
+   auto L = Reshape(lower.Write(), ne, ncomp);
+   auto U = Reshape(upper.Write(), ne, ncomp);
+
+   mfem::forall_2D<BLOCK_X>(ne*ncomp, BLOCK_X, 1,
+                            [=] MFEM_HOST_DEVICE (int ec)
+   {
+      const int e = ec % ne;
+      const int c = ec / ne;
+      const int vc = comp0 + c;
+      const real_t *coeff = &E(0, vc, e);
+      const int tid = MFEM_THREAD_ID(x);
+
+      MFEM_SHARED real_t sproj[MAX_ND];
+      MFEM_SHARED real_t ssum0[MAX_ND];
+      MFEM_SHARED real_t ssum1[MAX_ND];
+      MFEM_SHARED real_t smin[BLOCK_X];
+      MFEM_SHARED real_t smax[BLOCK_X];
+      MFEM_SHARED real_t sa0;
+      MFEM_SHARED real_t sa1;
+
+      MFEM_FOREACH_THREAD(i, x, nd)
+      {
+         if constexpr (T_PROJ)
+         {
+            const real_t x = data.xhat[i];
+            const real_t w = data.what[i];
+            ssum0[i] = 0.5*coeff[i]*w;
+            ssum1[i] = 1.5*coeff[i]*w*x;
+         }
+         else
+         {
+            ssum0[i] = 0.0;
+            ssum1[i] = 0.0;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(ii, x, 1)
+      {
+         sa0 = 0.0;
+         sa1 = 0.0;
+         for (int i = 0; i < nd; i++)
+         {
+            sa0 += ssum0[i];
+            sa1 += ssum1[i];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(i, x, nd)
+      {
+         if constexpr (T_PROJ)
+         {
+            const real_t x = data.xhat[i];
+            sproj[i] = coeff[i] - sa0 - sa1*x;
+         }
+         else
+         {
+            sproj[i] = coeff[i];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      real_t lower_local = HUGE_VAL;
+      real_t upper_local = -HUGE_VAL;
+      MFEM_FOREACH_THREAD(j, x, data.ncp)
+      {
+         real_t lo = 0.0;
+         real_t hi = 0.0;
+         if constexpr (T_PROJ)
+         {
+            const real_t xcp = data.cphat[j];
+            lo = sa0 + sa1*xcp;
+            hi = lo;
+         }
+
+         for (int i = 0; i < nd; i++)
+         {
+            const real_t val = sproj[i];
+            const real_t lv = data.lbound[j + i*data.ncp]*val;
+            const real_t uv = data.ubound[j + i*data.ncp]*val;
+            lo += lv < uv ? lv : uv;
+            hi += lv > uv ? lv : uv;
+         }
+         lower_local = lower_local < lo ? lower_local : lo;
+         upper_local = upper_local > hi ? upper_local : hi;
+      }
+
+      smin[tid] = lower_local;
+      smax[tid] = upper_local;
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(ii, x, 1)
+      {
+         real_t lower_ec = smin[0];
+         real_t upper_ec = smax[0];
+         const int nthreads = MFEM_THREAD_SIZE(x);
+         const int nactive = data.ncp < nthreads ? data.ncp : nthreads;
+         for (int t = 1; t < nactive; t++)
+         {
+            lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
+            upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
+         }
+         L(e, c) = lower_ec;
+         U(e, c) = upper_ec;
+      }
+   });
+}
+
+template<int T_NB = 0, int T_NCP = 0, bool T_PROJ = true>
+inline void GetElementBoundsKernel2D(const PLBoundDeviceData &data,
+                                     const int fes_vdim,
+                                     const int ne,
+                                     const Vector &e_vec,
+                                     Vector &lower,
+                                     Vector &upper,
+                                     const int comp0,
+                                     const int ncomp)
+{
+   constexpr int DEFAULT_MAX_NB = 8;
+   constexpr int DEFAULT_MAX_CP = 3*DEFAULT_MAX_NB;
+   constexpr int MAX_NB = T_NB ? T_NB : DEFAULT_MAX_NB;
+   constexpr int MAX_CP = T_NCP ? T_NCP : DEFAULT_MAX_CP;
+   constexpr int MAX_THREADS = MAX_CP*MAX_CP;
+
+   const int nb = data.nb;
+   const int ncp = data.ncp;
+   const int nd = nb*nb;
+   MFEM_VERIFY(nb <= MAX_NB,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "1D degree bound.");
+   MFEM_VERIFY(ncp <= MAX_CP,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "control-point bound.");
+   MFEM_VERIFY(ncp*ncp <= MAX_THREADS,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "thread-block bound.");
+
+   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
+   auto L = Reshape(lower.Write(), ne, ncomp);
+   auto U = Reshape(upper.Write(), ne, ncomp);
+
+   mfem::forall_2D<MAX_THREADS>(ne*ncomp, ncp, ncp,
+                                [=] MFEM_HOST_DEVICE (int ec)
+   {
+      const int e = ec % ne;
+      const int c = ec / ne;
+      const int vc = comp0 + c;
+      const real_t *coeff = &E(0, vc, e);
+      const int tx = MFEM_THREAD_ID(x);
+      const int ty = MFEM_THREAD_ID(y);
+
+      MFEM_SHARED real_t sproj[MAX_NB*MAX_NB];
+      MFEM_SHARED real_t srow_min[MAX_NB*MAX_CP];
+      MFEM_SHARED real_t srow_max[MAX_NB*MAX_CP];
+      MFEM_SHARED real_t srow_a0[MAX_NB];
+      MFEM_SHARED real_t srow_a1[MAX_NB];
+      MFEM_SHARED real_t sa0[MAX_CP];
+      MFEM_SHARED real_t sa1[MAX_CP];
+      MFEM_SHARED real_t smin[MAX_THREADS];
+      MFEM_SHARED real_t smax[MAX_THREADS];
+
+      // Stage 1a: for each nodal row, form the per-node contributions to the
+      // row-wise linear fit used by the first 1D bounding solve.
+      MFEM_FOREACH_THREAD(jrow, y, nb)
+      {
+         const real_t *row_coeff = coeff + jrow*nb;
+         const int row_ncp_off = jrow*MAX_CP;
+         MFEM_FOREACH_THREAD(i, x, nb)
+         {
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[i];
+               const real_t w = data.what[i];
+               srow_min[row_ncp_off + i] = 0.5*row_coeff[i]*w;
+               srow_max[row_ncp_off + i] = 1.5*row_coeff[i]*w*x;
+            }
+            else
+            {
+               srow_min[row_ncp_off + i] = 0.0;
+               srow_max[row_ncp_off + i] = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 1b: reduce the row-wise projection coefficients a0/a1.
+      if constexpr (T_PROJ)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_ncp_off = jrow*MAX_CP;
+            real_t a0 = 0.0;
+            real_t a1 = 0.0;
+            MFEM_FOREACH_THREAD(ii, x, 1)
+            {
+               for (int i = 0; i < nb; i++)
+               {
+                  a0 += srow_min[row_ncp_off + i];
+                  a1 += srow_max[row_ncp_off + i];
+               }
+               srow_a0[jrow] = a0;
+               srow_a1[jrow] = a1;
+            }
+         }
+         MFEM_SYNC_THREAD;
+      }
+
+      // Stage 1c: subtract the row-wise linear fit once and cache the
+      // projected row coefficients for reuse across all x-control points.
+      MFEM_FOREACH_THREAD(jrow, y, nb)
+      {
+         const real_t *row_coeff = coeff + jrow*nb;
+         MFEM_FOREACH_THREAD(i, x, nb)
+         {
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[i];
+               sproj[jrow*MAX_NB + i] = row_coeff[i]
+                                        - srow_a0[jrow] - srow_a1[jrow]*x;
+            }
+            else
+            {
+               sproj[jrow*MAX_NB + i] = row_coeff[i];
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 1d: solve the first 1D bounding problem along each nodal row and
+      // store bounds at every x-direction control point.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_cp_off = jrow*ncp;
+            real_t lo = 0.0;
+            real_t hi = 0.0;
+            if constexpr (T_PROJ)
+            {
+               const real_t xcp = data.cphat[icp];
+               lo = srow_a0[jrow] + srow_a1[jrow]*xcp;
+               hi = lo;
+            }
+            for (int i = 0; i < nb; i++)
+            {
+               const real_t val = sproj[jrow*MAX_NB + i];
+               const real_t lv = data.lbound[icp + i*data.ncp]*val;
+               const real_t uv = data.ubound[icp + i*data.ncp]*val;
+               lo += lv < uv ? lv : uv;
+               hi += lv > uv ? lv : uv;
+            }
+            srow_min[row_cp_off + icp] = lo;
+            srow_max[row_cp_off + icp] = hi;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2a: from the row bounds, form the per-row contributions to the
+      // second 1D projection solve in the y-direction.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_cp_off = jrow*ncp;
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[jrow];
+               const real_t w = data.what[jrow];
+               const real_t t = 0.5*(srow_min[row_cp_off + icp] +
+                                     srow_max[row_cp_off + icp]);
+               smin[row_cp_off + icp] = 0.5*t*w;
+               smax[row_cp_off + icp] = 1.5*t*w*x;
+            }
+            else
+            {
+               smin[row_cp_off + icp] = 0.0;
+               smax[row_cp_off + icp] = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2b: reduce the y-direction projection coefficients for each
+      // x-control-point column.
+      MFEM_FOREACH_THREAD(jj, y, 1)
+      {
+         MFEM_FOREACH_THREAD(icp, x, ncp)
+         {
+            real_t a0 = 0.0;
+            real_t a1 = 0.0;
+            for (int jrow = 0; jrow < nb; jrow++)
+            {
+               a0 += smin[jrow*ncp + icp];
+               a1 += smax[jrow*ncp + icp];
+            }
+            sa0[icp] = a0;
+            sa1[icp] = a1;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2c: subtract the y-direction linear fit from the intermediate
+      // row bounds so the final tensor-product bound uses the perturbation.
+      if constexpr (T_PROJ)
+      {
+         MFEM_FOREACH_THREAD(icp, x, ncp)
+         {
+            MFEM_FOREACH_THREAD(jrow, y, nb)
+            {
+               const int row_cp_off = jrow*ncp;
+               const real_t x = data.xhat[jrow];
+               const real_t t = sa0[icp] + sa1[icp]*x;
+               srow_min[row_cp_off + icp] -= t;
+               srow_max[row_cp_off + icp] -= t;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 3: each thread now owns one 2D control point (icp, kcp) and
+      // accumulates its final lower/upper bound from the row-bound data.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(kcp, y, ncp)
+         {
+            real_t lo = 0.0;
+            real_t hi = 0.0;
+            if constexpr (T_PROJ)
+            {
+               const real_t xcp = data.cphat[kcp];
+               lo = sa0[icp] + sa1[icp]*xcp;
+               hi = lo;
+            }
+            for (int jrow = 0; jrow < nb; jrow++)
+            {
+               const real_t w0 = srow_min[jrow*ncp + icp];
+               const real_t w1 = srow_max[jrow*ncp + icp];
+               const real_t lb = data.lbound[kcp + jrow*data.ncp];
+               const real_t ub = data.ubound[kcp + jrow*data.ncp];
+               const real_t v0 = lb*w0;
+               const real_t v1 = ub*w0;
+               const real_t v2 = lb*w1;
+               const real_t v3 = ub*w1;
+               real_t vlo = v0 < v1 ? v0 : v1;
+               real_t vhi = v0 > v1 ? v0 : v1;
+               vlo = vlo < v2 ? vlo : v2;
+               vlo = vlo < v3 ? vlo : v3;
+               vhi = vhi > v2 ? vhi : v2;
+               vhi = vhi > v3 ? vhi : v3;
+               lo += vlo;
+               hi += vhi;
+            }
+            const int slot = kcp*ncp + icp;
+            smin[slot] = lo;
+            smax[slot] = hi;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      const int lane = ty*ncp + tx;
+      const int nactive = ncp*ncp;
+      const int nthreads = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
+
+      // Reduce all 2D control-point bounds to one lower/upper pair per
+      // (element, component).
+      if (nthreads == 1)
+      {
+         if (tx == 0 && ty == 0)
+         {
+            real_t lower_ec = smin[0];
+            real_t upper_ec = smax[0];
+            for (int t = 1; t < nactive; t++)
+            {
+               lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
+               upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
+            }
+            L(e, c) = lower_ec;
+            U(e, c) = upper_ec;
+         }
+      }
+      else
+      {
+         for (int stride = (nactive + 1)/2; stride > 0;
+              stride = (stride + 1)/2)
+         {
+            if (lane < stride && lane + stride < nactive)
+            {
+               smin[lane] = smin[lane] < smin[lane + stride] ?
+                            smin[lane] : smin[lane + stride];
+               smax[lane] = smax[lane] > smax[lane + stride] ?
+                            smax[lane] : smax[lane + stride];
+            }
+            MFEM_SYNC_THREAD;
+            if (stride == 1) { break; }
+         }
+
+         if (lane == 0)
+         {
+            L(e, c) = smin[0];
+            U(e, c) = smax[0];
+         }
+      }
+   });
+}
+
+} // namespace internal
+
+inline void PLBound::GetElementBoundsKernel(const int rdim, const int fes_vdim,
+                                            const Vector &e_vec,
+                                            Vector &lower, Vector &upper,
+                                            const int vdim) const
+{
+   MFEM_VERIFY(b_type != BasisType::Positive,
+               "Bernstein device bounds are not implemented.");
+   if (rdim == 3)
+   {
+      MFEM_ABORT("Device element bounds kernel currently only supports 1D/2D.");
+   }
+   MFEM_VERIFY(rdim == 1 || rdim == 2, "Invalid element dimension.");
+   MFEM_VERIFY(vdim >= -1 && vdim <= fes_vdim, "Invalid vector component.");
+   const int nd = static_cast<int>(std::pow(nb, rdim));
+   const int ne = e_vec.Size()/(nd*fes_vdim);
+   const int ncomp = (vdim > 0) ? 1 : fes_vdim;
+
+   lower.SetSize(ne*ncomp, e_vec);
+   upper.SetSize(ne*ncomp, e_vec);
+   lower.UseDevice(true);
+   upper.UseDevice(true);
+
+   if (!proj)
+   {
+      MFEM_ABORT("Device element bounds kernel currently requires projection "
+                 "enabled.");
+   }
+
+   const real_t *dxhat = xhat.Read();
+   const real_t *dwhat = what.Read();
+   const real_t *dcphat = cphat.Read();
+   const real_t *dlbound = lbound.Read();
+   const real_t *dubound = ubound.Read();
+
+   internal::PLBoundDeviceData data
+   {
+      nb,
+      ncp,
+      dxhat,
+      dwhat,
+      dcphat,
+      dlbound,
+      dubound
+   };
+
+   const int comp0 = (vdim > 0) ? (vdim - 1) : 0;
+
+   if (rdim == 1)
+   {
+      switch (nb)
+      {
+         case 2: return internal::GetElementBoundsKernel1D<2, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 3: return internal::GetElementBoundsKernel1D<3, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 4: return internal::GetElementBoundsKernel1D<4, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 5: return internal::GetElementBoundsKernel1D<5, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 6: return internal::GetElementBoundsKernel1D<6, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 7: return internal::GetElementBoundsKernel1D<7, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 8: return internal::GetElementBoundsKernel1D<8, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 9: return internal::GetElementBoundsKernel1D<9, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 10: return internal::GetElementBoundsKernel1D<10, true>(data, fes_vdim, ne,
+                                                                       e_vec, lower, upper,
+                                                                       comp0, ncomp);
+         default: return internal::GetElementBoundsKernel1D<0, true>(data, fes_vdim, ne,
+                                                                      e_vec, lower, upper,
+                                                                      comp0, ncomp);
+      }
+   }
+#define MFEM_PLBOUND_2D_DISPATCH(NB, NCP) \
+   return internal::GetElementBoundsKernel2D<NB, NCP, true>(data, fes_vdim, ne, \
+                                                            e_vec, lower, upper, \
+                                                            comp0, ncomp)
+   switch (nb)
+   {
+      case 2:
+         switch (ncp)
+         {
+            case 4: MFEM_PLBOUND_2D_DISPATCH(2, 4);
+            case 6: MFEM_PLBOUND_2D_DISPATCH(2, 6);
+            case 8: MFEM_PLBOUND_2D_DISPATCH(2, 8);
+         }
+         break;
+      case 3:
+         switch (ncp)
+         {
+            case 6: MFEM_PLBOUND_2D_DISPATCH(3, 6);
+            case 9: MFEM_PLBOUND_2D_DISPATCH(3, 9);
+            case 12: MFEM_PLBOUND_2D_DISPATCH(3, 12);
+         }
+         break;
+      case 4:
+         switch (ncp)
+         {
+            case 8: MFEM_PLBOUND_2D_DISPATCH(4, 8);
+            case 12: MFEM_PLBOUND_2D_DISPATCH(4, 12);
+            case 16: MFEM_PLBOUND_2D_DISPATCH(4, 16);
+         }
+         break;
+      case 5:
+         switch (ncp)
+         {
+            case 10: MFEM_PLBOUND_2D_DISPATCH(5, 10);
+            case 15: MFEM_PLBOUND_2D_DISPATCH(5, 15);
+            case 20: MFEM_PLBOUND_2D_DISPATCH(5, 20);
+         }
+         break;
+      case 6:
+         switch (ncp)
+         {
+            case 12: MFEM_PLBOUND_2D_DISPATCH(6, 12);
+            case 18: MFEM_PLBOUND_2D_DISPATCH(6, 18);
+            case 24: MFEM_PLBOUND_2D_DISPATCH(6, 24);
+         }
+         break;
+      case 7:
+         switch (ncp)
+         {
+            case 14: MFEM_PLBOUND_2D_DISPATCH(7, 14);
+            case 21: MFEM_PLBOUND_2D_DISPATCH(7, 21);
+            case 28: MFEM_PLBOUND_2D_DISPATCH(7, 28);
+         }
+         break;
+      case 8:
+         switch (ncp)
+         {
+            case 16: MFEM_PLBOUND_2D_DISPATCH(8, 16);
+            case 24: MFEM_PLBOUND_2D_DISPATCH(8, 24);
+            case 32: MFEM_PLBOUND_2D_DISPATCH(8, 32);
+         }
+         break;
+   }
+#undef MFEM_PLBOUND_2D_DISPATCH
+   return internal::GetElementBoundsKernel2D<0, 0, true>(data, fes_vdim, ne,
+                                                         e_vec, lower, upper,
+                                                         comp0, ncomp);
+}
+
 } // namespace mfem
 
 #endif // MFEM_BOUNDS

fem/dgmassinv_kernels.hpp

@@ -57,7 +57,7 @@ void DGMassApply(const int e,
       }
       else if (DIM == 3)
       {
-         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
+         SmemPAMassApply3D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
       }
       else
       {

fem/gridfunc.cpp

@@ -5252,6 +5252,30 @@ void GridFunction::GetElementBounds(const PLBound &plb,
                                     Vector &lower, Vector &upper,
                                     const int vdim) const
 {
+   if (UseDevice() && Device::Allows(Backend::DEVICE_MASK) &&
+       plb.GetBasisType() != BasisType::Positive &&
+       UsesTensorBasis(*fes))
+   {
+      const FiniteElement &fe = *fes->GetTypicalFE();
+      const int rdim = fe.GetDim();
+      const int fes_dim = fes->GetVDim();
+      const int nel = fes->GetNE();
+      const int nd = fe.GetDof();
+
+      Vector e_vec(nd*fes_dim*nel, Device::GetDeviceMemoryType());
+      e_vec.UseDevice(true);
+      const ElementRestrictionOperator *elem_restr =
+         fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      MFEM_VERIFY(elem_restr != nullptr,
+                  "Element restriction is required for device bounds.");
+      elem_restr->Mult(*this, e_vec);
+
+      plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower, upper, vdim);
+      lower.HostRead();
+      upper.HostRead();
+      return;
+   }
+
    int nel  = fes->GetNE();
    int fes_dim  = fes->GetVDim();
    lower.SetSize(nel*(vdim > 0 ? 1 :fes_dim));

fem/gslib.cpp

@@ -106,9 +106,7 @@ FindPointsGSLIB::FindPointsGSLIB()
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized = 0;
-   MPI_Initialized(&initialized);
-   if (!initialized) { MPI_Init(NULL, NULL); }
+   if (!Mpi::IsInitialized()) { Mpi::Init(); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else
@@ -2624,9 +2622,7 @@ GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized;
-   MPI_Initialized(&initialized);
-   if (!initialized) { MPI_Init(NULL, NULL); }
+   if (!Mpi::IsInitialized()) { Mpi::Init(); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else

fem/integ/bilininteg_mass_kernels.hpp

@@ -181,6 +181,12 @@ constexpr int NBZ(int D1D)
 {
    return ipow(2, D(D1D) >= 0 ? D(D1D) : 0);
 }
+constexpr int NBZ3D(int MDQ)
+{
+   return MDQ > 0 ? std::min<int>(
+             (128 + MDQ * MDQ * MDQ - 1) / (MDQ * MDQ * MDQ), 64)
+          : 1;
+}
 }
 
 // Shared memory PA Mass Diagonal 2D kernel
@@ -804,19 +810,23 @@ void PAMassApply3D_Element(const int e,
    }
 }
 
-template<int T_D1D, int T_Q1D, bool ACCUMULATE = true>
-MFEM_HOST_DEVICE inline
-void SmemPAMassApply3D_Element(const int e,
-                               const int NE,
-                               const real_t *b_,
-                               const real_t *d_,
-                               const real_t *x_,
-                               real_t *y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+template <int T_D1D, int T_Q1D, int TBATCH, bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline void
+SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
+                          const real_t *d_, const real_t *x_, real_t *y_,
+                          int d1d = 0, int q1d = 0)
 {
-   constexpr int D1D = T_D1D ? T_D1D : d1d;
-   constexpr int Q1D = T_Q1D ? T_Q1D : q1d;
+   static_assert(TBATCH > 0, "TBATCH must be positive");
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+   constexpr int tbatch = TBATCH;
+   const int tidz = MFEM_THREAD_ID(z);
+#else
+   // host always batch size 1
+   constexpr int tbatch = 1;
+   constexpr int tidz = 0;
+#endif
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
    constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
    constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
    constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
@@ -829,33 +839,37 @@ void SmemPAMassApply3D_Element(const int e,
    MFEM_SHARED real_t sDQ[MQ1*MD1];
    real_t (*B)[MD1] = (real_t (*)[MD1]) sDQ;
    real_t (*Bt)[MQ1] = (real_t (*)[MQ1]) sDQ;
-   MFEM_SHARED real_t sm0[MDQ*MDQ*MDQ];
-   MFEM_SHARED real_t sm1[MDQ*MDQ*MDQ];
-   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) sm0;
-   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) sm1;
-   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm0;
-   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm1;
-   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) sm0;
-   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) sm1;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_SHARED real_t sm0[tbatch][MDQ*MDQ*MDQ];
+   MFEM_SHARED real_t sm1[tbatch][MDQ*MDQ*MDQ];
+   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) (sm0+tidz);
+   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) (sm1+tidz);
+   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm0+tidz);
+   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm1+tidz);
+   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) (sm0+tidz);
+   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) (sm1+tidz);
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          MFEM_UNROLL(MD1)
          for (int dz = 0; dz < D1D; ++dz)
          {
-            X[dz][dy][dx] = x(dx,dy,dz,e);
+            X[dz][dy][dx] = x(dx, dy, dz, e);
          }
       }
-      MFEM_FOREACH_THREAD(dx,x,Q1D)
+      MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
+   }
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(dy, y, D1D)
       {
-         B[dx][dy] = b(dx,dy);
+         MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -880,9 +894,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -907,9 +921,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -929,22 +943,22 @@ void SmemPAMassApply3D_Element(const int e,
          MFEM_UNROLL(MQ1)
          for (int qz = 0; qz < Q1D; qz++)
          {
-            QQQ[qz][qy][qx] = u[qz] * d(qx,qy,qz,e);
+            QQQ[qz][qy][qx] = u[qz] * d(qx, qy, qz, e);
          }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(di,y,D1D)
+   if (tidz == 0)
    {
-      MFEM_FOREACH_THREAD(q,x,Q1D)
+      MFEM_FOREACH_THREAD(di, y, D1D)
       {
-         Bt[di][q] = b(q,di);
+         MFEM_FOREACH_THREAD(q, x, Q1D) { Bt[di][q] = b(q, di); }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -969,9 +983,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -996,9 +1010,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -1020,11 +1034,11 @@ void SmemPAMassApply3D_Element(const int e,
          {
             if (ACCUMULATE)
             {
-               y(dx,dy,dz,e) += u[dz];
+               y(dx, dy, dz, e) += u[dz];
             }
             else
             {
-               y(dx,dy,dz,e) = u[dz];
+               y(dx, dy, dz, e) = u[dz];
             }
          }
       }
@@ -1115,8 +1129,8 @@ inline void PAMassApply3D(const int NE,
    });
 }
 
-// Shared memory PA Mass Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
+// Shared memory PA Mass Apply 3D kernel
+template<int T_D1D = 0, int T_Q1D = 0, int TBATCH=1>
 inline void SmemPAMassApply3D(const int NE,
                               const Array<real_t> &b_,
                               const Array<real_t> &bt_,
@@ -1126,6 +1140,9 @@ inline void SmemPAMassApply3D(const int NE,
                               const int d1d = 0,
                               const int q1d = 0)
 {
+   static_assert(T_D1D > 0, "T_D1D must be positive");
+   static_assert(T_Q1D > 0, "T_Q1D must be positive");
+   static_assert(TBATCH > 0, "TBATCH must be positive");
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1137,9 +1154,11 @@ inline void SmemPAMassApply3D(const int NE,
    const auto d = d_.Read();
    const auto x = x_.Read();
    auto y = y_.ReadWrite();
-   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall_2D_batch<T_Q1D * T_Q1D * TBATCH>(NE, Q1D, Q1D, TBATCH,
+                                                 [=] MFEM_HOST_DEVICE(int e)
    {
-      internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
+      internal::SmemPAMassApply3D_Element<T_D1D, T_Q1D, TBATCH>(e, NE, b, d, x,
+                                                                y, d1d, q1d);
    });
 }
 
@@ -1394,7 +1413,16 @@ ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
 {
    if constexpr (DIM == 1) { return internal::PAMassApply1D; }
    else if constexpr (DIM == 2) { return internal::SmemPAMassApply2D<T_D1D,T_Q1D>; }
-   else if constexpr (DIM == 3) { return internal::SmemPAMassApply3D<T_D1D, T_Q1D>; }
+   else if constexpr (DIM == 3)
+   {
+      constexpr int MDQ = T_D1D >= T_Q1D ? T_D1D : T_Q1D;
+      // max 64 threads in z limit in cuda and hip
+      if constexpr (MDQ > 0)
+      {
+         return internal::SmemPAMassApply3D<T_D1D, T_Q1D,
+                internal::mass::NBZ3D(MDQ)>;
+      }
+   }
    MFEM_ABORT("");
 }
 

general/forall.hpp

@@ -1090,6 +1090,12 @@ inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
    ForallWrap<2>(true, N, body, X, Y, BZ);
 }
 
+template<int MAX_THREADS_PER_BLOCK, typename lambda>
+inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
+{
+   ForallWrap<2, MAX_THREADS_PER_BLOCK>(true, N, body, X, Y, BZ);
+}
+
 template<typename lambda>
 inline void forall_3D(int N, int X, int Y, int Z, lambda &&body)
 {

makefile

@@ -123,15 +123,20 @@ EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
-MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
+MINIAPP_ALL_SUBDIRS = common electromagnetics meshing performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
  hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
  fluids/navier fluids/schrodinger-flow plasma plasma/pic
+MINIAPP_RECURSIVE_SUBDIRS = plasma/pic
+MINIAPP_SUBDIRS := $(filter-out \
+   $(MINIAPP_RECURSIVE_SUBDIRS),$(MINIAPP_ALL_SUBDIRS))
+MINIAPP_ALL_DIRS := $(addprefix miniapps/,$(MINIAPP_ALL_SUBDIRS))
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
+ toys gslib shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
 
+EM_ALL_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_ALL_DIRS)
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 
 TEST_SUBDIRS = unit
@@ -146,7 +151,7 @@ MFEM_BUILD_DIR ?= .
 BUILD_DIR := $(MFEM_BUILD_DIR)
 BUILD_REAL_DIR := $(abspath $(BUILD_DIR))
 ifneq ($(BUILD_REAL_DIR),$(MFEM_REAL_DIR))
-   BUILD_SUBDIRS = $(DIRS) config $(EM_DIRS) doc $(TEST_DIRS)
+   BUILD_SUBDIRS = $(DIRS) config $(EM_ALL_DIRS) doc $(TEST_DIRS)
    CONFIG_FILE_DEF = -DMFEM_CONFIG_FILE='"$(BUILD_REAL_DIR)/config/_config.hpp"'
    BLD := $(if $(BUILD_REAL_DIR:$(CURDIR)=),$(BUILD_DIR)/,)
    $(if $(word 2,$(BLD)),$(error Spaces in BLD = "$(BLD)" are not supported))
@@ -483,10 +488,10 @@ $(OBJECT_FILES): $(BLD)%.o: $(SRC)%.cpp $(CONFIG_MK)
 
 all: examples miniapps $(TEST_DIRS)
 
-.PHONY: miniapps $(EM_DIRS) $(TEST_DIRS)
+.PHONY: miniapps $(EM_ALL_DIRS) $(TEST_DIRS)
 miniapps: $(MINIAPP_DIRS)
 $(MINIAPP_USE_COMMON): miniapps/common
-$(EM_DIRS) $(TEST_DIRS): lib
+$(EM_ALL_DIRS) $(TEST_DIRS): lib
 	$(MAKE) -C $(BLD)$(@)
 
 .PHONY: doc
@@ -694,7 +699,7 @@ local-config:
 .PHONY: build-config
 build-config:
 	for d in $(BUILD_SUBDIRS); do mkdir -p $(BLD)$${d}; done
-	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS)); do \
+	for dir in "" $(addsuffix /,config $(EM_ALL_DIRS) doc $(TEST_DIRS)); do\
 	   printf "# Auto-generated file.\n%s\n%s\n" \
 	      "MFEM_DIR = $(MFEM_REAL_DIR)" \
 	      "include \$$(MFEM_DIR)/$${dir}makefile" \
@@ -796,13 +801,15 @@ status info:
 
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
-FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
+FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_ALL_DIRS) config,$(dir)/*.?pp)
 TESTS_SUBDIRS = unit benchmarks convergence mem_manager par-mesh-format
-UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme
-MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners hooke/materials hooke/kernels
+UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme dfem
+MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
+   hooke/materials hooke/kernels
 FORMAT_FILES += $(foreach dir,$(TESTS_SUBDIRS),tests/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
+FORMAT_FILES += config/cmake/config.hpp.in config/config.hpp.in mfem*.hpp
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))
 
@@ -833,14 +840,29 @@ mfem_check_command = \
 # Verify the C++ code styling in MFEM and check that std::cout and std::cerr are
 # not used in the library (use mfem::out and mfem::err instead).
 style:
-	@echo "Applying C++ code style..."
 	@astyle_version="$$($(ASTYLE_BIN) --version)";\
 	 if [ "$$astyle_version" != $(ASTYLE_VER) ]; then\
 	    printf "%s\n" "Invalid astyle version: '$$astyle_version'"\
 	           "Please use: '"$(ASTYLE_VER)"'";\
 	    exit 1;\
 	 fi
-	@err_code=0;\
+	@err_code=0; \
+	if command -v git 2>&1 > /dev/null && [ -d $(MFEM_DIR)/.git ]; then \
+	   echo "Checking if all git files are selected for formatting ..."; \
+	   ls -1 $(FORMAT_FILES) | sort > format-files-make.txt; \
+	   git -C $(MFEM_DIR) ls-files '*.[ch]pp*' | sort \
+	      > format-files-git.txt; \
+	   cat format-files-make.txt format-files-git.txt | sort | uniq \
+	      > format-files-make-plus-git.txt; \
+	   rm -f format-files-git.txt; \
+	   $(call mfem_check_command,\
+	      diff format-files-make.txt format-files-make-plus-git.txt | \
+	      grep "^> ",\
+	      "All git files are selected for formatting",\
+	      "The above git files are NOT selected for formatting"); \
+	   rm -f format-files-make.txt format-files-make-plus-git.txt; \
+	fi; \
+	echo "Applying C++ code style...";\
 	$(call mfem_check_command,\
 	    $(ASTYLE) $(FORMAT_LIST) | grep Formatted,\
 	    "No source files were changed",\

miniapps/gslib/CMakeLists.txt

@@ -46,8 +46,16 @@ if (MFEM_USE_GSLIB)
 
   if (MFEM_ENABLE_TESTING)
     foreach (test "schwarz_ex1" "field-diff" "findpts" "field-interp")
-      add_test(NAME ${test}
-        COMMAND $<TARGET_FILE:${test}> -no-vis)
+      if (MFEM_USE_MPI)
+        add_test(NAME ${test}
+          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1
+          ${MPIEXEC_PREFLAGS}
+          $<TARGET_FILE:${test}> -no-vis
+          ${MPIEXEC_POSTFLAGS})
+      else()
+        add_test(NAME ${test}
+          COMMAND $<TARGET_FILE:${test}> -no-vis)
+      endif()
     endforeach()
   endif()
 

miniapps/gslib/makefile

@@ -83,11 +83,16 @@ include $(MFEM_TEST_MK)
 
 # Testing: Parallel vs. serial runs
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+ifeq ($(MFEM_USE_MPI),YES)
+   RUN_MPI_SER = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) 1
+else
+   RUN_MPI_SER =
+endif
 TEST_NAME := GSLIB miniapp
 %-test-par: %
 	@$(call mfem-test,$<, $(RUN_MPI), $(TEST_NAME))
 %-test-seq: %
-	@$(call mfem-test,$<,, $(TEST_NAME))
+	@$(call mfem-test,$<, $(RUN_MPI_SER), $(TEST_NAME))
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

miniapps/plasma/makefile

@@ -28,6 +28,12 @@ endif
 
 PLASMA_SUBDIRS = pic
 
+SUBDIRS_ALL = $(addsuffix /all,$(PLASMA_SUBDIRS))
+SUBDIRS_TEST = $(addsuffix /test,$(PLASMA_SUBDIRS))
+SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(PLASMA_SUBDIRS))
+SUBDIRS_CLEAN = $(addsuffix /clean,$(PLASMA_SUBDIRS))
+SUBDIRS_TPRINT = $(addsuffix /test-print,$(PLASMA_SUBDIRS))
+
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all lib-common clean clean-build clean-exec
@@ -37,31 +43,24 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
-
-COMMON_O=
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 # Remove built-in rules
 %: %.cpp
 %.o: %.cpp
 
-all: $(MINIAPPS) subdirs
+all: $(MINIAPPS) $(SUBDIRS_ALL) 
 
-.PHONY: subdirs $(PLASMA_SUBDIRS)
-subdirs: $(PLASMA_SUBDIRS)
-$(PLASMA_SUBDIRS): lib-common
-	$(MAKE) -C $(BLD)$(@)
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
+   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
+	$(MAKE) -C $(@D) $(@F)
+$(SUBDIRS_TPRINT):
+	@$(MAKE) -C $(@D) $(@F)
 
 # Rules for building the miniapps
-%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_O) $(COMMON_LIB) \
-        $(MFEM_LIBS)
-
-# Rules for compiling miniapp dependencies
-$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
-%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_LIB) $(MFEM_LIBS)
 
 # Rule for building lib-common
 lib-common:
@@ -69,6 +68,9 @@ lib-common:
 
 MFEM_TESTS = MINIAPPS
 include $(MFEM_TEST_MK)
+test: $(SUBDIRS_TEST)
+test-noclean: $(SUBDIRS_TEST_NOCLEAN)
+test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Specific execution options
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
@@ -79,14 +81,9 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-ALL_CLEAN_SUBDIRS = $(addsuffix /clean,$(PLASMA_SUBDIRS))
-.PHONY: $(ALL_CLEAN_SUBDIRS)
-$(ALL_CLEAN_SUBDIRS):
-	$(MAKE) -C $(BLD)$(@D) $(@F)
+clean: clean-build clean-exec $(SUBDIRS_CLEAN)
 
-clean: clean-build clean-exec
-
-clean-build: $(addsuffix /clean,$(PLASMA_SUBDIRS))
+clean-build:
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints
 

miniapps/plasma/pic/makefile

@@ -20,6 +20,7 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
+SEQ_MINIAPPS =
 PAR_MINIAPPS =
 
 ifeq ($(MFEM_USE_GSLIB),YES)
@@ -27,9 +28,9 @@ ifeq ($(MFEM_USE_GSLIB),YES)
 endif
 
 ifeq ($(MFEM_USE_MPI),NO)
-   MINIAPPS =
+   MINIAPPS = $(SEQ_MINIAPPS)
 else
-   MINIAPPS = $(PAR_MINIAPPS)
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
 endif
 
 .SUFFIXES:
@@ -50,7 +51,7 @@ COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
 all: $(MINIAPPS)
 
 # Rules for building the miniapps
-electrostatic-pic: electrostatic-pic.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
 

miniapps/tools/CMakeLists.txt

@@ -61,6 +61,10 @@ if (MFEM_USE_MPI)
       LIBRARIES mfem-common)
       add_dependencies(gridfunction-bounds copy_miniapps_tools_data)
 
+    add_mfem_miniapp(random-gridfunction-bounds
+      MAIN random-gridfunction-bounds.cpp
+      LIBRARIES mfem)
+
     add_mfem_miniapp(plor-transfer
       MAIN plor-transfer.cpp LIBRARIES mfem)
 

miniapps/tools/makefile

@@ -23,7 +23,8 @@ MFEM_LIB_FILE = mfem_is_not_built
 SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
 	tmop-check-metric tmop-metric-magnitude compare-dc
 
-PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds
+PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds \
+	random-gridfunction-bounds
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -79,7 +80,7 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 # Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer, plor-transfer
 NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer \
 	plor-transfer tmop-check-metric tmop-metric-magnitude gridfunction-bounds \
-	compare-dc
+	random-gridfunction-bounds compare-dc
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

miniapps/tools/random-gridfunction-bounds.cpp

@@ -0,0 +1,263 @@
+// 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.
+//
+//    ---------------------------------------------------------------------
+//     Compute bounds of a random grid function on a generated tensor mesh
+//    ---------------------------------------------------------------------
+//
+// This miniapp generates a 1D segment mesh or 2D quad mesh, builds a random
+// discontinuous grid function, computes element-wise piecewise linear bounds,
+// and visualizes the input field together with the lower and upper bounds.
+//
+// Compile with: make random-gridfunction-bounds
+//
+// Sample runs:
+//   mpirun -np 4 random-gridfunction-bounds
+//   mpirun -np 4 random-gridfunction-bounds -nx 64 -o 6 -ref 3 -d hip
+
+#include "mfem.hpp"
+
+#include <algorithm>
+#include <type_traits>
+
+using namespace mfem;
+using namespace std;
+
+void VisualizeField(ParMesh &pmesh, ParGridFunction &input,
+                    char *title, int pos_x, int pos_y);
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   int dim = 2;
+   int nx = 16;
+   int order = 4;
+   int num_comp = 2;
+   int ref = 2;
+   int niter = 1000;
+   int seed = 12345;
+   bool kernel_only = true;
+   bool visualization = false;
+   const char *device_config = "cpu";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&dim, "-dim", "--dimension",
+                  "Dimension of the generated tensor-product mesh (1 or 2).");
+   args.AddOption(&nx, "-nx", "--num-elements",
+                  "Number of elements in each mesh direction.");
+   args.AddOption(&order, "-o", "--order",
+                  "Polynomial degree of the random discontinuous field.");
+   args.AddOption(&num_comp, "-nc", "--num-components",
+                  "Number of vector components in the ParFiniteElementSpace.");
+   args.AddOption(&ref, "-ref", "--piecewise-linear-ref-factor",
+                  "Scaling factor for the resolution of the piecewise linear "
+                  "bounds. If less than 2, the resolution is picked "
+                  "automatically.");
+   args.AddOption(&niter, "-ni", "--num-iters",
+                  "Number of times to evaluate the bounds.");
+   args.AddOption(&seed, "-rs", "--random-seed",
+                  "Random seed used to initialize the field.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&kernel_only, "-ko", "--kernel-only",
+                  "-no-ko", "--no-kernel-only",
+                  "Run only PLBound::GetElementBoundsKernel on a prebuilt "
+                  "element E-vector.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   MFEM_VERIFY(dim == 1 || dim == 2, "dim must be 1 or 2.");
+   MFEM_VERIFY(nx > 0, "nx must be positive.");
+   MFEM_VERIFY(order >= 0, "order must be non-negative.");
+   MFEM_VERIFY(num_comp > 0, "num_comp must be positive.");
+   MFEM_VERIFY(niter > 0, "niter must be positive.");
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   Mesh mesh = (dim == 1) ?
+               Mesh::MakeCartesian1D(nx, 1.0) :
+               Mesh::MakeCartesian2D(nx, nx, Element::QUADRILATERAL, true,
+                                     1.0, 1.0);
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+   const int mesh_dim = pmesh.Dimension();
+   L2_FECollection fec(order, mesh_dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&pmesh, &fec, num_comp, Ordering::byNODES);
+   ParGridFunction input(&fes);
+   input.Randomize(seed + Mpi::WorldRank());
+   input.UseDevice(true);
+
+   L2_FECollection fec_pc(0, mesh_dim);
+   ParFiniteElementSpace fes_pc(&pmesh, &fec_pc, num_comp, Ordering::byNODES);
+   ParGridFunction lowerb(&fes_pc), upperb(&fes_pc);
+   Vector lower_vec, upper_vec;
+
+   PLBound plb(&fes, ref*(fes.GetMaxElementOrder() + 1));
+   if (kernel_only)
+   {
+      const FiniteElement &fe = *fes.GetTypicalFE();
+      const int rdim = fe.GetDim();
+      const int nd = fe.GetDof();
+      const int fes_dim = fes.GetVDim();
+      Vector e_vec(nd*fes_dim*fes.GetNE(), Device::GetDeviceMemoryType());
+      e_vec.UseDevice(true);
+
+      const ElementRestrictionOperator *elem_restr =
+         fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      MFEM_VERIFY(elem_restr != nullptr,
+                  "Element restriction is required for kernel-only mode.");
+      elem_restr->Mult(input, e_vec);
+
+      for (int i = 0; i < niter; i++)
+      {
+         plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower_vec, upper_vec);
+      }
+   }
+   else
+   {
+      for (int i = 0; i < niter; i++)
+      {
+         input.GetElementBounds(plb, lower_vec, upper_vec);
+      }
+   }
+
+   const real_t *lower_data = lower_vec.HostRead();
+   const real_t *upper_data = upper_vec.HostRead();
+
+   // Build a host reference from the lexicographic E-vector and the scalar
+   // PLBound::GetNDBounds path to avoid re-entering the device dispatch.
+   const bool use_dev = input.UseDevice();
+   PLBound plb_host(&fes, ref*(fes.GetMaxElementOrder() + 1));
+   Vector lower_ref, upper_ref;
+   const FiniteElement &fe = *fes.GetTypicalFE();
+   const int rdim = fe.GetDim();
+   const int nd = fe.GetDof();
+   const int nel = fes.GetNE();
+   const int fes_dim = fes.GetVDim();
+   Vector e_vec_ref(nd*fes_dim*nel);
+   lower_ref.SetSize(nel*fes_dim);
+   upper_ref.SetSize(nel*fes_dim);
+   const ElementRestrictionOperator *elem_restr =
+      fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+   MFEM_VERIFY(elem_restr != nullptr,
+               "Element restriction is required for host reference.");
+   input.UseDevice(false);
+   input.HostRead();
+   elem_restr->Mult(input, e_vec_ref);
+   input.UseDevice(use_dev);
+   const real_t *e_ref_data = e_vec_ref.HostRead();
+
+   for (int d = 0; d < fes_dim; d++)
+   {
+      for (int e = 0; e < nel; e++)
+      {
+         Vector coeff(nd);
+         for (int i = 0; i < nd; i++)
+         {
+            coeff(i) = e_ref_data[i + nd*(d + fes_dim*e)];
+         }
+         Vector lower_c, upper_c;
+         plb_host.GetNDBounds(rdim, coeff, lower_c, upper_c);
+         lower_ref(e + d*nel) = lower_c.Min();
+         upper_ref(e + d*nel) = upper_c.Max();
+      }
+   }
+   const real_t *lower_ref_data = lower_ref.HostRead();
+   const real_t *upper_ref_data = upper_ref.HostRead();
+
+   MFEM_VERIFY(lower_vec.Size() == lower_ref.Size() &&
+               upper_vec.Size() == upper_ref.Size(),
+               "Reference element-bound vectors have inconsistent sizes.");
+
+   real_t lower_diff = 0.0;
+   real_t upper_diff = 0.0;
+   for (int i = 0; i < lower_vec.Size(); i++)
+   {
+      lower_diff = std::max(lower_diff,
+                            std::abs(lower_data[i] - lower_ref_data[i]));
+   }
+   for (int i = 0; i < upper_vec.Size(); i++)
+   {
+      upper_diff = std::max(upper_diff,
+                            std::abs(upper_data[i] - upper_ref_data[i]));
+   }
+   MPI_Allreduce(MPI_IN_PLACE, &lower_diff, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MAX, pmesh.GetComm());
+   MPI_Allreduce(MPI_IN_PLACE, &upper_diff, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MAX, pmesh.GetComm());
+
+   const real_t verify_tol = std::is_same<real_t, float>::value ?
+                             real_t(1.0e-5) : real_t(1.0e-12);
+   MFEM_VERIFY(lower_diff <= verify_tol && upper_diff <= verify_tol,
+               "Device element bounds do not match host reference.");
+
+   lowerb = lower_vec;
+   upperb = upper_vec;
+
+   real_t lower_min = lowerb.Min();
+   real_t upper_max = upperb.Max();
+   MPI_Allreduce(MPI_IN_PLACE, &lower_min, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MIN, pmesh.GetComm());
+   MPI_Allreduce(MPI_IN_PLACE, &upper_max, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MAX, pmesh.GetComm());
+
+   if (Mpi::Root())
+   {
+      cout << "dim: " << mesh_dim << '\n'
+           << "nx: " << nx << '\n'
+           << "order: " << order << '\n'
+           << "num components: " << num_comp << '\n'
+           << "PL bound control-point factor: " << ref << '\n'
+           << "iterations: " << niter << '\n'
+           << "kernel-only mode: " << (kernel_only ? "yes" : "no") << '\n'
+           << "host/device lower max diff: " << lower_diff << '\n'
+           << "host/device upper max diff: " << upper_diff << '\n'
+           << "global lower bound minimum: " << lower_min << '\n'
+           << "global upper bound maximum: " << upper_max << endl;
+   }
+
+   if (visualization)
+   {
+      char title1[] = "Random input gridfunction";
+      char title2[] = "Element-wise lower bound";
+      char title3[] = "Element-wise upper bound";
+      VisualizeField(pmesh, input, title1, 0, 0);
+      VisualizeField(pmesh, lowerb, title2, 450, 0);
+      VisualizeField(pmesh, upperb, title3, 900, 0);
+   }
+
+   return 0;
+}
+
+void VisualizeField(ParMesh &pmesh, ParGridFunction &input,
+                    char *title, int pos_x, int pos_y)
+{
+   socketstream sock;
+   if (pmesh.GetMyRank() == 0)
+   {
+      sock.open("localhost", 19916);
+      sock << "solution\n";
+   }
+   pmesh.PrintAsOne(sock);
+   input.SaveAsOne(sock);
+   if (pmesh.GetMyRank() == 0)
+   {
+      sock << "window_title '" << title << "'\n"
+           << "window_geometry "
+           << pos_x << " " << pos_y << " " << 400 << " " << 400 << "\n"
+           << "keys jRmclApppppppppppp//]]]]]]]]" << endl;
+   }
+}