style

Add toggle for SuperLU_DIST device offload

fem/dgmassinv_kernels.hpp

@@ -387,7 +387,7 @@ void DGMassInverse::DGMassCGIteration(const Vector &b_, Vector &u_) const
 
    static constexpr int NB = Q1D ? Q1D : 1; // block size
 
-   mfem::forall_2D(NE, NB, NB, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall_2D<NB*NB>(NE, NB, NB, [=] MFEM_HOST_DEVICE (int e)
    {
       // Perform change of basis if needed
       if (CHANGE_BASIS)

fem/integ/bilininteg_diffusion_kernels.hpp

@@ -1004,13 +1004,16 @@ inline void SmemPADiffusionApply3D(const int NE,
    const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
    MFEM_VERIFY(D1D <= max_d1d, "");
    MFEM_VERIFY(Q1D <= max_q1d, "");
-   auto b = Reshape(b_.Read(), Q1D, D1D);
-   auto g = Reshape(g_.Read(), Q1D, D1D);
-   auto d = Reshape(d_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
-   auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
+   const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto g = Reshape(g_.Read(), Q1D, D1D);
+   const auto d = Reshape(d_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
+   const auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
    auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
    MFEM_VERIFY(D1D <= Q1D, "THREAD_DIRECT requires D1D <= Q1D");
-   mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+
+   mfem::forall_3D<T_Q1D*T_Q1D*T_Q1D>(NE,
+                                      Q1D, Q1D, Q1D,
+                                      [=] MFEM_HOST_DEVICE (int e)
    {
       const int D1D = T_D1D ? T_D1D : d1d;
       const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/integ/bilininteg_mass_kernels.hpp

@@ -1133,11 +1133,11 @@ inline void SmemPAMassApply3D(const int NE,
    const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
    MFEM_VERIFY(D1D <= max_d1d, "");
    MFEM_VERIFY(Q1D <= max_q1d, "");
-   auto b = b_.Read();
-   auto d = d_.Read();
-   auto x = x_.Read();
+   const auto b = b_.Read();
+   const auto d = d_.Read();
+   const auto x = x_.Read();
    auto y = y_.ReadWrite();
-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
    {
       internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
    });
@@ -1156,8 +1156,8 @@ inline void EAMassAssemble1D(const int NE,
    const int Q1D = T_Q1D ? T_Q1D : q1d;
    MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
    MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
-   auto B = Reshape(basis.Read(), Q1D, D1D);
-   auto D = Reshape(padata.Read(), Q1D, NE);
+   const auto B = Reshape(basis.Read(), Q1D, D1D);
+   const auto D = Reshape(padata.Read(), Q1D, NE);
    auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, NE);
    mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
    {

fem/integ/bilininteg_vecdiffusion_pa.hpp

@@ -54,7 +54,7 @@ void SmemPAVectorDiffusionApply2D(const int NE,
    const auto XE = Reshape(x.Read(), D1D, D1D, SDIM, NE);
    auto YE = Reshape(y.ReadWrite(), D1D, D1D, SDIM, NE);
 
-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
    {
       constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
       constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
@@ -120,7 +120,7 @@ void SmemPAVectorDiffusionApply3D(const int NE,
    const auto XE = Reshape(x.Read(), D1D, D1D, D1D, SDIM, NE);
    auto YE = Reshape(y.ReadWrite(), D1D, D1D, D1D, SDIM, NE);
 
-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
    {
       constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
       constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;

fem/integ/bilininteg_vecmass_pa.hpp

@@ -51,7 +51,7 @@ void SmemPAVectorMassApply2D(const int NE,
    const auto X = Reshape(x.Read(), D1D, D1D, VDIM, NE);
    auto Y = Reshape(y.ReadWrite(), D1D, D1D, VDIM, NE);
 
-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
    {
       constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
       constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
@@ -119,7 +119,7 @@ void SmemPAVectorMassApply3D(const int NE,
    const auto X = Reshape(x.Read(), D1D, D1D, D1D, VDIM, NE);
    auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, VDIM, NE);
 
-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
    {
       constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
       constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;

fem/transfer.cpp

@@ -1030,52 +1030,12 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAProlongateTranspose(
    BatchedLinAlg::MultTranspose(P_dt, x, y);
 }
 
-namespace
-{
-
-class H1ConsistentMassOperator : public Operator
-{
-private:
-   const Operator &M_LH;
-   const Solver &M_L_solver;
-
-public:
-   H1ConsistentMassOperator(const Operator &M_LH_, const Solver &M_L_solver_)
-      : Operator(M_LH_.Height(), M_LH_.Width()),
-        M_LH(M_LH_),
-        M_L_solver(M_L_solver_)
-   {
-      MFEM_VERIFY(M_LH.Height() == M_L_solver.Height() &&
-                  M_LH.Height() == M_L_solver.Width(),
-                  "incompatible consistent mass operator dimensions");
-   }
-
-   void Mult(const Vector &x, Vector &y) const override
-   {
-      Vector tmp(M_LH.Height());
-      M_LH.Mult(x, tmp);
-      M_L_solver.Mult(tmp, y);
-   }
-
-   void MultTranspose(const Vector &x, Vector &y) const override
-   {
-      Vector tmp(M_LH.Height());
-      M_L_solver.Mult(x, tmp);
-      M_LH.MultTranspose(tmp, y);
-   }
-};
-
-}
-
 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_,
-   const bool use_ea_, const bool use_consistent_mass_, MemoryType d_mt_)
+   const bool use_ea_, MemoryType d_mt_)
    : L2Projection(fes_ho_, fes_lor_, d_mt_),
-     use_ea(use_ea_),
-     use_consistent_mass(use_consistent_mass_)
+     use_ea(use_ea_)
 {
-   MFEM_VERIFY(!(use_ea && use_consistent_mass),
-               "consistent mass is not supported with element assembly");
 
    // need scalar to keep dimensions matching (operators are built to apply
    // individually on each vdim)
@@ -1093,7 +1053,7 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
 
    std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
 
-   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH(!use_consistent_mass);
+   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH();
 
    const SparseMatrix *P_ho = fes_ho_scalar->GetConformingProlongation();
    const SparseMatrix *P_lor = fes_lor_scalar->GetConformingProlongation();
@@ -1102,68 +1062,40 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    {
       if (P_ho && P_lor)
       {
-         if (R_mat) { R_mat.reset(RAP(*P_lor, *R_mat, *P_ho)); }
+         R_mat.reset(RAP(*P_lor, *R_mat, *P_ho));
          M_LH_mat.reset(RAP(*P_lor, *M_LH_mat, *P_ho));
       }
       else if (P_ho)
       {
-         if (R_mat) { R_mat.reset(mfem::Mult(*R_mat, *P_ho)); }
+         R_mat.reset(mfem::Mult(*R_mat, *P_ho));
          M_LH_mat.reset(mfem::Mult(*M_LH_mat, *P_ho));
       }
       else // P_lor != nullptr
       {
-         if (R_mat) { R_mat.reset(mfem::Mult(*P_lor, *R_mat)); }
+         R_mat.reset(mfem::Mult(*P_lor, *R_mat));
          M_LH_mat.reset(mfem::Mult(*P_lor, *M_LH_mat));
       }
    }
 
-   if (use_consistent_mass)
-   {
-      BilinearForm M_lor(fes_lor_scalar.get());
-      M_lor.AddDomainIntegrator(new MassIntegrator);
-      M_lor.Assemble();
-      M_lor.Finalize(0);
-      SparseMatrix *M_L_mat = new SparseMatrix(M_lor.SpMat());
+   SparseMatrix *RTxM_LH_mat = TransposeMult(*R_mat, *M_LH_mat);
+   precon.reset(new DSmoother(*RTxM_LH_mat));
 
-      ML_precon.reset(new DSmoother(*M_L_mat));
-      ML_pcg.SetPrintLevel(0);
-      ML_pcg.SetMaxIter(1000);
-      ML_pcg.SetRelTol(1e-13);
-      ML_pcg.SetAbsTol(1e-13);
-      ML_pcg.SetPreconditioner(*ML_precon);
-      ML_pcg.SetOperator(*M_L_mat);
+   // Set ownership
+   RTxM_LH.reset(RTxM_LH_mat);
+   R = std::move(R_mat);
+   M_LH = std::move(M_LH_mat);
 
-      M_L.reset(M_L_mat);
-      M_LH = std::move(M_LH_mat);
-      R.reset(new H1ConsistentMassOperator(*M_LH, ML_pcg));
-   }
-   else
-   {
-      SparseMatrix *RTxM_LH_mat = TransposeMult(*R_mat, *M_LH_mat);
-      precon.reset(new DSmoother(*RTxM_LH_mat));
-
-      // Set ownership
-      RTxM_LH.reset(RTxM_LH_mat);
-      R = std::move(R_mat);
-      M_LH = std::move(M_LH_mat);
-
-      SetupPCG();
-   }
+   SetupPCG();
 }
 
 #ifdef MFEM_USE_MPI
 
 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor,
-   const bool use_ea_, const bool use_consistent_mass_, MemoryType d_mt_)
+   const bool use_ea_, MemoryType d_mt_)
    : L2Projection(pfes_ho, pfes_lor, d_mt_),
-     use_ea(use_ea_),
-     use_consistent_mass(use_consistent_mass_),
-     pcg(pfes_ho.GetComm()),
-     ML_pcg(pfes_ho.GetComm())
+     use_ea(use_ea_), pcg(pfes_ho.GetComm())
 {
-   MFEM_VERIFY(!(use_ea && use_consistent_mass),
-               "consistent mass is not supported with element assembly");
 
    // need scalar to keep dimensions matching (operators are built to apply
    // individually on each vdim)
@@ -1179,58 +1111,26 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
       return;
    }
 
-   std::unique_ptr<SparseMatrix> R_local_sp, M_LH_local_sp;
-   std::tie(R_local_sp, M_LH_local_sp) =
-      ComputeSparseRAndM_LH(!use_consistent_mass);
+   std::tie(R, M_LH) = ComputeSparseRAndM_LH();
 
-   M_LH_local_sp->Finalize(0);
-   if (R_local_sp) { R_local_sp->Finalize(0); }
-
-   HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho.GetComm(),
-                                              pfes_lor_scalar->GlobalVSize(),
-                                              pfes_ho_scalar->GlobalVSize(),
-                                              pfes_lor_scalar->GetDofOffsets(),
-                                              pfes_ho_scalar->GetDofOffsets(),
-                                              M_LH_local_sp.get());
-   HypreParMatrix *M_LH_mat = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(),
-                                  &M_LH_local, pfes_ho_scalar->Dof_TrueDof_Matrix());
-
-   if (use_consistent_mass)
-   {
-      ParBilinearForm M_lor(pfes_lor_scalar.get());
-      M_lor.AddDomainIntegrator(new MassIntegrator);
-      M_lor.Assemble();
-      M_lor.Finalize(0);
-      HypreParMatrix *M_L_mat = M_lor.ParallelAssemble();
-
-      ML_pcg.SetPrintLevel(0);
-      ML_pcg.SetMaxIter(1000);
-      ML_pcg.SetRelTol(1e-13);
-      ML_pcg.SetAbsTol(1e-13);
-      ML_pcg.SetOperator(*M_L_mat);
-
-      M_L.reset(M_L_mat);
-      M_LH.reset(M_LH_mat);
-      HyprePCG *ML_hypre_pcg = new HyprePCG(*M_L_mat);
-      ML_hypre_pcg->SetPrintLevel(0);
-      ML_hypre_pcg->SetMaxIter(1000);
-      ML_hypre_pcg->SetTol(1e-13);
-      ML_hypre_pcg->SetAbsTol(1e-13);
-      ML_hypre_pcg->SetZeroInitialIterate();
-      ML_solver.reset(ML_hypre_pcg);
-      R.reset(new H1ConsistentMassOperator(*M_LH, *ML_solver));
-      return;
-   }
 
    HypreParMatrix R_local = HypreParMatrix(pfes_ho.GetComm(),
                                            pfes_lor_scalar->GlobalVSize(),
                                            pfes_ho_scalar->GlobalVSize(),
                                            pfes_lor_scalar->GetDofOffsets(),
                                            pfes_ho_scalar->GetDofOffsets(),
-                                           R_local_sp.get());
+                                           static_cast<SparseMatrix*>(R.get()));
+   HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho.GetComm(),
+                                              pfes_lor_scalar->GlobalVSize(),
+                                              pfes_ho_scalar->GlobalVSize(),
+                                              pfes_lor_scalar->GetDofOffsets(),
+                                              pfes_ho_scalar->GetDofOffsets(),
+                                              static_cast<SparseMatrix*>(M_LH.get()));
 
    HypreParMatrix *R_mat = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(),
                                &R_local, pfes_ho_scalar->Dof_TrueDof_Matrix());
+   HypreParMatrix *M_LH_mat = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(),
+                                  &M_LH_local, pfes_ho_scalar->Dof_TrueDof_Matrix());
 
    std::unique_ptr<HypreParMatrix> R_T(R_mat->Transpose());
    HypreParMatrix *RTxM_LH_mat = ParMult(R_T.get(), M_LH_mat, true);
@@ -1538,8 +1438,6 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
 void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
    const Vector& x, Vector& y) const
 {
-   MFEM_VERIFY(!use_consistent_mass,
-               "BackwardOperator is not supported with consistent mass");
 
    Vector X(fes_lor.GetTrueVSize());
    Vector X_dim(M_LH->Height());
@@ -1571,9 +1469,6 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
 void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
    const Vector& x, Vector& y) const
 {
-   MFEM_VERIFY(!use_consistent_mass,
-               "BackwardOperator is not supported with consistent mass");
-
    Vector X(fes_ho.GetTrueVSize());
    Vector X_dim(pcg.Width());
    Vector Xbar(pcg.Height());
@@ -1604,34 +1499,17 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(real_t p_rtol_)
 {
    pcg.SetRelTol(p_rtol_);
-   ML_pcg.SetRelTol(p_rtol_);
-#ifdef MFEM_USE_MPI
-   if (ML_solver)
-   {
-      HyprePCG *hypre_pcg = dynamic_cast<HyprePCG*>(ML_solver.get());
-      if (hypre_pcg) { hypre_pcg->SetTol(p_rtol_); }
-   }
-#endif
 }
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(real_t p_atol_)
 {
    pcg.SetAbsTol(p_atol_);
-   ML_pcg.SetAbsTol(p_atol_);
-#ifdef MFEM_USE_MPI
-   if (ML_solver)
-   {
-      HyprePCG *hypre_pcg = dynamic_cast<HyprePCG*>(ML_solver.get());
-      if (hypre_pcg) { hypre_pcg->SetAbsTol(p_atol_); }
-   }
-#endif
 }
 
 std::pair<
 std::unique_ptr<SparseMatrix>,
 std::unique_ptr<SparseMatrix>>
-L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH(
-   bool build_R)
+                            L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH()
 {
    std::pair<std::unique_ptr<SparseMatrix>,
        std::unique_ptr<SparseMatrix>> r_and_mlh;
@@ -1664,76 +1542,69 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH(
 
    BuildHo2Lor(nel_ho, nel_lor, cf_tr);
 
+   // ML_inv contains the inverse lumped (row sum) mass matrix. Note that the
+   // method will also work with a full (consistent) mass matrix, though this is
+   // not implemented here. L refers to the low-order refined mesh
    Vector ML_inv(ndof_lor);
-   if (build_R)
+   ML_inv = 0.0;
+
+   // Compute ML_inv
+   for (int iho = 0; iho < nel_ho; ++iho)
    {
-      // ML_inv contains the inverse lumped (row sum) mass matrix. L refers to
-      // the low-order refined mesh.
-      ML_inv = 0.0;
+      Array<int> lor_els;
+      ho2lor.GetRow(iho, lor_els);
+      int nref = ho2lor.RowSize(iho);
 
-      // Compute ML_inv
-      for (int iho = 0; iho < nel_ho; ++iho)
+      Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
+      const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
+      int nedof_lor = fe_lor.GetDof();
+
+      // Instead of using a MassIntegrator, manually loop over integration
+      // points so we can row sum and store the diagonal as a Vector.
+      Vector ML_el(nedof_lor);
+      Vector shape_lor(nedof_lor);
+      Array<int> dofs_lor(nedof_lor);
+
+      for (int iref = 0; iref < nref; ++iref)
       {
-         Array<int> lor_els;
-         ho2lor.GetRow(iho, lor_els);
-         int nref = ho2lor.RowSize(iho);
+         int ilor = lor_els[iref];
+         ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
 
-         Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
-         const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
-         int nedof_lor = fe_lor.GetDof();
-
-         // Instead of using a MassIntegrator, manually loop over integration
-         // points so we can row sum and store the diagonal as a Vector.
-         Vector ML_el(nedof_lor);
-         Vector shape_lor(nedof_lor);
-         Array<int> dofs_lor(nedof_lor);
-
-         for (int iref = 0; iref < nref; ++iref)
+         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
+         const IntegrationRule* ir = &IntRules.Get(geom, order);
+         ML_el = 0.0;
+         for (int i = 0; i < ir->GetNPoints(); ++i)
          {
-            int ilor = lor_els[iref];
-            ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
-
-            int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
-            const IntegrationRule* ir = &IntRules.Get(geom, order);
-            ML_el = 0.0;
-            for (int i = 0; i < ir->GetNPoints(); ++i)
-            {
-               const IntegrationPoint& ip_lor = ir->IntPoint(i);
-               fe_lor.CalcShape(ip_lor, shape_lor);
-               el_tr->SetIntPoint(&ip_lor);
-               ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
-            }
-            fes_lor.GetElementDofs(ilor, dofs_lor);
-            ML_inv.AddElementVector(dofs_lor, ML_el);
+            const IntegrationPoint& ip_lor = ir->IntPoint(i);
+            fe_lor.CalcShape(ip_lor, shape_lor);
+            el_tr->SetIntPoint(&ip_lor);
+            ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
          }
+         fes_lor.GetElementDofs(ilor, dofs_lor);
+         ML_inv.AddElementVector(dofs_lor, ML_el);
       }
-      // DOF by DOF inverse of non-zero entries
-      LumpedMassInverse(ML_inv);
    }
+   // DOF by DOF inverse of non-zero entries
+   LumpedMassInverse(ML_inv);
 
    // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
-   std::unique_ptr<SparseMatrix> pattern = AllocR();
-   if (build_R)
-   {
-      r_and_mlh.first = std::move(pattern);
-   }
+   r_and_mlh.first = AllocR();
    // Allocate M_LH (same sparsity pattern as R)
    // L refers to the low-order refined mesh (DOFs correspond to rows)
    // H refers to the higher-order mesh (DOFs correspond to columns)
-   SparseMatrix &pattern_mat = build_R ? *r_and_mlh.first : *pattern;
-   Memory<int> I(pattern_mat.Height() + 1);
-   for (int icol = 0; icol < pattern_mat.Height() + 1; ++icol)
+   Memory<int> I(r_and_mlh.first->Height() + 1);
+   for (int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol)
    {
-      I[icol] = pattern_mat.GetI()[icol];
+      I[icol] = r_and_mlh.first->GetI()[icol];
    }
-   Memory<int> J(pattern_mat.NumNonZeroElems());
-   for (int jcol = 0; jcol < pattern_mat.NumNonZeroElems(); ++jcol)
+   Memory<int> J(r_and_mlh.first->NumNonZeroElems());
+   for (int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol)
    {
-      J[jcol] = pattern_mat.GetJ()[jcol];
+      J[jcol] = r_and_mlh.first->GetJ()[jcol];
    }
    r_and_mlh.second = std::unique_ptr<SparseMatrix>(
-                         new SparseMatrix(I, J, NULL, pattern_mat.Height(),
-                                          pattern_mat.Width(), true, true, true));
+                         new SparseMatrix(I, J, NULL, r_and_mlh.first->Height(),
+                                          r_and_mlh.first->Width(), true, true, true));
 
    IntegrationPointTransformation ip_tr;
    IsoparametricTransformation& emb_tr = ip_tr.Transf;
@@ -1776,21 +1647,15 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH(
          Array<int> dofs_lor(nedof_lor);
          fes_lor.GetElementDofs(ilor, dofs_lor);
          Vector R_row;
-         if (build_R)
+         for (int i = 0; i < nedof_lor; ++i)
          {
-            for (int i = 0; i < nedof_lor; ++i)
-            {
-               M_LH_el.GetRow(i, R_row);
-               R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
-            }
+            M_LH_el.GetRow(i, R_row);
+            R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
          }
          Array<int> dofs_ho(nedof_ho);
          fes_ho.GetElementDofs(iho, dofs_ho);
          r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
-         if (build_R)
-         {
-            r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
-         }
+         r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
 
       }
    }
@@ -2144,10 +2009,6 @@ const Operator &L2ProjectionGridTransfer::ForwardOperator()
 
 const Operator &L2ProjectionGridTransfer::BackwardOperator()
 {
-   MFEM_VERIFY(!(use_consistent_mass && !force_l2_space &&
-                 dom_fes.FEColl()->GetContType() ==
-                 FiniteElementCollection::CONTINUOUS),
-               "BackwardOperator is not supported with consistent mass");
    if (!B)
    {
       if (!F) { BuildF(); }
@@ -2156,24 +2017,15 @@ const Operator &L2ProjectionGridTransfer::BackwardOperator()
    return *B;
 }
 
-void L2ProjectionGridTransfer::UseConsistentMass(bool use_consistent_mass_)
-{
-   MFEM_VERIFY(!F && !B,
-               "UseConsistentMass must be called before constructing operators");
-   use_consistent_mass = use_consistent_mass_;
-}
-
 void L2ProjectionGridTransfer::BuildF()
 {
-   MFEM_VERIFY(!(use_ea && use_consistent_mass),
-               "consistent mass is not supported with element assembly");
    if (!force_l2_space &&
        dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS)
    {
       if (!Parallel())
       {
          F = new L2ProjectionH1Space(dom_fes, ran_fes,
-                                     use_ea, use_consistent_mass, d_mt);
+                                     use_ea, d_mt);
       }
       else
       {
@@ -2183,7 +2035,7 @@ void L2ProjectionGridTransfer::BuildF()
          const mfem::ParFiniteElementSpace& ran_pfes =
             static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
          F = new L2ProjectionH1Space(dom_pfes, ran_pfes,
-                                     use_ea, use_consistent_mass, d_mt);
+                                     use_ea, d_mt);
 #endif
       }
    }
@@ -2196,11 +2048,6 @@ void L2ProjectionGridTransfer::BuildF()
 
 bool L2ProjectionGridTransfer::SupportsBackwardsOperator() const
 {
-   if (use_consistent_mass && !force_l2_space &&
-       dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS)
-   {
-      return false;
-   }
    return ran_fes.GetTrueVSize() >= dom_fes.GetTrueVSize();
 }
 

fem/transfer.hpp

@@ -169,10 +169,8 @@ public:
     is the forward transfer matrix, and M_f is the mass matrix on the coarse
     element. For L2 spaces, M_f is the mass matrix on the union of all fine
     elements comprising the coarse element. For H1 spaces, M_f is a diagonal
-    (lumped) mass matrix computed through row-summation, unless
-    UseConsistentMass() is enabled for the forward H1 operator. Note that the
-    backward transfer operator, B, is a left inverse of the forward transfer
-    operator, F,
+    (lumped) mass matrix computed through row-summation. Note that the backward
+    transfer operator, B, is a left inverse of the forward transfer operator, F,
     i.e. B F = I. Both F and B are defined in physical space and, generally for
     L2 spaces, vary between different mesh elements.
 
@@ -354,19 +352,16 @@ public:
    class L2ProjectionH1Space : public L2Projection
    {
       const bool use_ea;
-      const bool use_consistent_mass;
 
    public:
       L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
                           const FiniteElementSpace &fes_lor_,
                           const bool use_ea_,
-                          const bool use_consistent_mass_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
 #ifdef MFEM_USE_MPI
       L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
                           const ParFiniteElementSpace &pfes_lor_,
                           const bool use_ea_,
-                          const bool use_consistent_mass_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
 #endif
       /// Same as above but assembles action of R through 4 parts:
@@ -428,8 +423,7 @@ public:
       /// @brief Computes on-rank R and M_LH matrices. If true, computes mixed mass and/or
       /// inverse lumped mass matrix error when compared to device implementation.
       std::pair<std::unique_ptr<SparseMatrix>,
-          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH(
-             bool build_R = true);
+          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH();
 
       /// @brief Recovers vector of tdofs given a vector of dofs and a finite
       /// element space
@@ -460,11 +454,7 @@ public:
       std::unique_ptr<SparseMatrix> AllocR();
 
       CGSolver pcg;
-      CGSolver ML_pcg;
       std::unique_ptr<Solver> precon;
-      std::unique_ptr<Solver> ML_precon;
-      std::unique_ptr<Solver> ML_solver;
-      std::unique_ptr<Operator> M_L;
       // The restriction operator is represented as an Operator R. The
       // prolongation operator is a dense matrix computed as the inverse of (R^T
       // M_L R), and hence, is not stored.
@@ -488,6 +478,7 @@ public:
       Vector M_LH_ea;
       // Element Assembled lumped M_L inverse built via EA. Stores diagonal as a Ldof vector.
       Vector ML_inv_ea;
+
 #ifdef MFEM_USE_MPI
       std::unique_ptr<ParFiniteElementSpace> pfes_ho_scalar;
       std::unique_ptr<ParFiniteElementSpace> pfes_lor_scalar;
@@ -520,7 +511,6 @@ public:
    L2Projection   *F; ///< Forward, coarse-to-fine, operator
    L2Prolongation *B; ///< Backward, fine-to-coarse, operator
    bool force_l2_space;
-   bool use_consistent_mass;
 
 public:
    L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
@@ -528,15 +518,10 @@ public:
                             bool force_l2_space_ = false,
                             MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
       : GridTransfer(coarse_fes_, fine_fes_),
-        F(NULL), B(NULL), force_l2_space(force_l2_space_),
-        use_consistent_mass(false)
+        F(NULL), B(NULL), force_l2_space(force_l2_space_)
    { }
    virtual ~L2ProjectionGridTransfer();
 
-   /** Use the consistent low-order mass matrix in H1 non-EA Mult() and
-       MultTranspose(). This option does not support BackwardOperator(). */
-   void UseConsistentMass(bool use_consistent_mass_ = true);
-
    const Operator &ForwardOperator() override;
 
    const Operator &BackwardOperator() override;

general/backends.hpp

@@ -44,6 +44,7 @@
 #endif
 
 #if !defined(MFEM_USE_CUDA_OR_HIP)
+constexpr bool mfem_use_gpu = false;
 #define MFEM_DEVICE
 #define MFEM_HOST
 #define MFEM_LAMBDA
@@ -52,6 +53,7 @@
 #define MFEM_DEVICE_SYNC
 // MFEM_STREAM_SYNC is used for UVM and MPI GPU-Aware kernels
 #define MFEM_STREAM_SYNC
+#define MFEM_LAUNCH_BOUNDS(...)
 #endif
 
 #if !((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \

general/cuda.hpp

@@ -20,9 +20,11 @@
 
 #if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_USE_CUDA_OR_HIP
+constexpr bool mfem_use_gpu = true;
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
 #define MFEM_LAMBDA __host__
+#define MFEM_LAUNCH_BOUNDS __launch_bounds__
 // #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
 #define MFEM_DEVICE_SYNC MFEM_GPU_CHECK(cudaDeviceSynchronize())
 #define MFEM_STREAM_SYNC MFEM_GPU_CHECK(cudaStreamSynchronize(0))

general/forall.hpp

@@ -295,11 +295,12 @@ using hip_threads_z =
 #endif
 
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_CUDA) && defined(__CUDACC__)
-template <const int BLOCKS = MFEM_CUDA_BLOCKS, typename DBODY>
+template <typename DBODY>
 void RajaCuWrap1D(const int N, DBODY &&d_body)
 {
    //true denotes asynchronous kernel
-   RAJA::forall<RAJA::cuda_exec<BLOCKS,true>>(RAJA::RangeSegment(0,N),d_body);
+   RAJA::forall<RAJA::cuda_exec<MFEM_CUDA_BLOCKS,true>>(RAJA::RangeSegment(0,N),
+                                                        d_body);
 }
 
 template <typename DBODY>
@@ -362,18 +363,18 @@ struct RajaCuWrap;
 template <>
 struct RajaCuWrap<1>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
-      RajaCuWrap1D<BLCK>(N, d_body);
+      RajaCuWrap1D(N, d_body);
    }
 };
 
 template <>
 struct RajaCuWrap<2>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -384,7 +385,7 @@ struct RajaCuWrap<2>
 template <>
 struct RajaCuWrap<3>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -395,11 +396,12 @@ struct RajaCuWrap<3>
 #endif
 
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_HIP) && defined(__HIP__)
-template <const int BLOCKS = MFEM_HIP_BLOCKS, typename DBODY>
+template <typename DBODY>
 void RajaHipWrap1D(const int N, DBODY &&d_body)
 {
    //true denotes asynchronous kernel
-   RAJA::forall<RAJA::hip_exec<BLOCKS,true>>(RAJA::RangeSegment(0,N),d_body);
+   RAJA::forall<RAJA::hip_exec<MFEM_HIP_BLOCKS,true>>(RAJA::RangeSegment(0,N),
+                                                      d_body);
 }
 
 template <typename DBODY>
@@ -462,18 +464,18 @@ struct RajaHipWrap;
 template <>
 struct RajaHipWrap<1>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
-      RajaHipWrap1D<BLCK>(N, d_body);
+      RajaHipWrap1D(N, d_body);
    }
 };
 
 template <>
 struct RajaHipWrap<2>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -484,7 +486,7 @@ struct RajaHipWrap<2>
 template <>
 struct RajaHipWrap<3>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -584,12 +586,31 @@ void CuKernel2D(const int N, BODY body)
    body(k);
 }
 
+// __launch_bounds__ second argument is omitted to get the default behavior
+template <int MAX_THREADS_PER_BLOCK, typename BODY>
+__global__
+MFEM_LAUNCH_BOUNDS(MAX_THREADS_PER_BLOCK)
+static void CuKernel2DLaunchBounds(const int N, BODY body)
+{
+   const int k = blockIdx.x*blockDim.z + threadIdx.z;
+   if (k >= N) { return; }
+   body(k);
+}
+
 template <typename BODY> __global__ static
 void CuKernel3D(const int N, BODY body)
 {
    for (int k = blockIdx.x; k < N; k += gridDim.x) { body(k); }
 }
 
+template <int MAX_THREADS_PER_BLOCK, typename BODY>
+__global__
+MFEM_LAUNCH_BOUNDS(MAX_THREADS_PER_BLOCK)
+static void CuKernel3DLaunchBounds(const int N, BODY body)
+{
+   for (int k = blockIdx.x; k < N; k += gridDim.x) { body(k); }
+}
+
 template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
 void CuWrap1D(const int N, DBODY &&d_body)
 {
@@ -604,6 +625,8 @@ void CuWrap2D(const int N, DBODY &&d_body,
               const int X, const int Y, const int BZ)
 {
    if (N==0) { return; }
+   // required for optimized GCC/NVCC builds to prevent runtime
+   // ODR/linkage violations of inlined templated kernel helpers
    MFEM_VERIFY(BZ>0, "");
    const int GRID = (N+BZ-1)/BZ;
    const dim3 BLCK(X,Y,BZ);
@@ -611,6 +634,19 @@ void CuWrap2D(const int N, DBODY &&d_body,
    MFEM_GPU_CHECK(cudaGetLastError());
 }
 
+template <int MAX_THREADS_PER_BLOCK, typename DBODY>
+void CuWrap2DLaunchBounds(const int N, DBODY &&d_body,
+                          const int X, const int Y, const int BZ)
+{
+   if (N==0) { return; }
+   MFEM_VERIFY(BZ>0, "");
+   const int GRID = (N+BZ-1)/BZ;
+   const dim3 BLCK(X,Y,BZ);
+   static_assert(MAX_THREADS_PER_BLOCK > 0);
+   CuKernel2DLaunchBounds<MAX_THREADS_PER_BLOCK><<<GRID,BLCK>>>(N, d_body);
+   MFEM_GPU_CHECK(cudaGetLastError());
+}
+
 template <typename DBODY>
 void CuWrap3D(const int N, DBODY &&d_body,
               const int X, const int Y, const int Z, const int G)
@@ -622,24 +658,35 @@ void CuWrap3D(const int N, DBODY &&d_body,
    MFEM_GPU_CHECK(cudaGetLastError());
 }
 
-template <int Dim>
-struct CuWrap;
-
-template <>
-struct CuWrap<1>
+template <int MAX_THREADS_PER_BLOCK, typename DBODY>
+void CuWrap3DLaunchBounds(const int N, DBODY &&d_body,
+                          const int X, const int Y, const int Z, const int G)
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   if (N==0) { return; }
+   const int GRID = G == 0 ? N : G;
+   const dim3 BLCK(X,Y,Z);
+   static_assert(MAX_THREADS_PER_BLOCK > 0);
+   CuKernel3DLaunchBounds<MAX_THREADS_PER_BLOCK><<<GRID, BLCK>>>(N, d_body);
+   MFEM_GPU_CHECK(cudaGetLastError());
+}
+
+template <int Dim, int MAX_THREADS_PER_BLOCK> struct CuWrap;
+
+template <int MAX_THREADS_PER_BLOCK>
+struct CuWrap<1, MAX_THREADS_PER_BLOCK>
+{
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
-      CuWrap1D<BLCK>(N, d_body);
+      CuWrap1D<MFEM_CUDA_BLOCKS>(N, d_body);
    }
 };
 
 template <>
-struct CuWrap<2>
+struct CuWrap<2, 0>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -647,10 +694,22 @@ struct CuWrap<2>
    }
 };
 
-template <>
-struct CuWrap<3>
+template <int MAX_THREADS_PER_BLOCK>
+struct CuWrap<2, MAX_THREADS_PER_BLOCK>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
+   static void run(const int N, DBODY &&d_body,
+                   const int X, const int Y, const int Z, const int G)
+   {
+      static_assert(MAX_THREADS_PER_BLOCK > 0);
+      CuWrap2DLaunchBounds<MAX_THREADS_PER_BLOCK>(N, d_body, X, Y, Z);
+   }
+};
+
+template <>
+struct CuWrap<3, 0>
+{
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -658,6 +717,17 @@ struct CuWrap<3>
    }
 };
 
+template <int MAX_THREADS_PER_BLOCK>
+struct CuWrap<3, MAX_THREADS_PER_BLOCK>
+{
+   template <typename DBODY>
+   static void run(const int N, DBODY &&d_body,
+                   const int X, const int Y, const int Z, const int G)
+   {
+      CuWrap3DLaunchBounds<MAX_THREADS_PER_BLOCK>(N, d_body, X, Y, Z, G);
+   }
+};
+
 #endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 
 
@@ -680,13 +750,31 @@ void HipKernel2D(const int N, BODY body)
    body(k);
 }
 
+template <int MAX_THREADS_PER_BLOCK, typename BODY>
+__global__
+MFEM_LAUNCH_BOUNDS(MAX_THREADS_PER_BLOCK)
+static void HipKernel2DLaunchBounds(const int N, BODY body)
+{
+   const int k = hipBlockIdx_x*hipBlockDim_z + hipThreadIdx_z;
+   if (k >= N) { return; }
+   body(k);
+}
+
 template <typename BODY> __global__ static
 void HipKernel3D(const int N, BODY body)
 {
    for (int k = hipBlockIdx_x; k < N; k += hipGridDim_x) { body(k); }
 }
 
-template <const int BLCK = MFEM_HIP_BLOCKS, typename DBODY>
+template <int MAX_THREADS_PER_BLOCK, typename BODY>
+__global__
+MFEM_LAUNCH_BOUNDS(MAX_THREADS_PER_BLOCK)
+static void HipKernel3DLaunchBounds(const int N, BODY body)
+{
+   for (int k = hipBlockIdx_x; k < N; k += hipGridDim_x) { body(k); }
+}
+
+template <int BLCK = MFEM_HIP_BLOCKS, typename DBODY>
 void HipWrap1D(const int N, DBODY &&d_body)
 {
    if (N==0) { return; }
@@ -700,12 +788,27 @@ void HipWrap2D(const int N, DBODY &&d_body,
                const int X, const int Y, const int BZ)
 {
    if (N==0) { return; }
+   MFEM_VERIFY(BZ>0, "");
    const int GRID = (N+BZ-1)/BZ;
    const dim3 BLCK(X,Y,BZ);
    hipLaunchKernelGGL(HipKernel2D,GRID,BLCK,0,nullptr,N,d_body);
    MFEM_GPU_CHECK(hipGetLastError());
 }
 
+template <int MAX_THREADS_PER_BLOCK, typename DBODY>
+void HipWrap2DLaunchBounds(const int N, DBODY &&d_body,
+                           const int X, const int Y, const int BZ)
+{
+   if (N==0) { return; }
+   MFEM_VERIFY(BZ>0, "");
+   const int GRID = (N+BZ-1)/BZ;
+   const dim3 BLCK(X,Y,BZ);
+   static_assert(MAX_THREADS_PER_BLOCK > 0);
+   HipKernel2DLaunchBounds<MAX_THREADS_PER_BLOCK><<<dim3(GRID), dim3(BLCK), 0, 0>>>
+   (N, d_body);
+   MFEM_GPU_CHECK(hipGetLastError());
+}
+
 template <typename DBODY>
 void HipWrap3D(const int N, DBODY &&d_body,
                const int X, const int Y, const int Z, const int G)
@@ -717,24 +820,36 @@ void HipWrap3D(const int N, DBODY &&d_body,
    MFEM_GPU_CHECK(hipGetLastError());
 }
 
-template <int Dim>
-struct HipWrap;
-
-template <>
-struct HipWrap<1>
+template <int MAX_THREADS_PER_BLOCK, typename DBODY>
+void HipWrap3DLaunchBounds(const int N, DBODY &&d_body,
+                           const int X, const int Y, const int Z, const int G)
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   if (N==0) { return; }
+   const int GRID = G == 0 ? N : G;
+   const dim3 BLCK(X,Y,Z);
+   static_assert(MAX_THREADS_PER_BLOCK > 0);
+   HipKernel3DLaunchBounds<MAX_THREADS_PER_BLOCK><<<dim3(GRID), dim3(BLCK), 0, 0>>>
+   (N, d_body);
+   MFEM_GPU_CHECK(hipGetLastError());
+}
+
+template <int Dim, int MAX_THREADS_PER_BLOCK> struct HipWrap;
+
+template <int MAX_THREADS_PER_BLOCK>
+struct HipWrap<1, MAX_THREADS_PER_BLOCK>
+{
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
-      HipWrap1D<BLCK>(N, d_body);
+      HipWrap1D<MFEM_HIP_BLOCKS>(N, d_body);
    }
 };
 
 template <>
-struct HipWrap<2>
+struct HipWrap<2, 0>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -742,10 +857,21 @@ struct HipWrap<2>
    }
 };
 
-template <>
-struct HipWrap<3>
+template <int MAX_THREADS_PER_BLOCK>
+struct HipWrap<2, MAX_THREADS_PER_BLOCK>
 {
-   template <const int BLCK = MFEM_CUDA_BLOCKS, typename DBODY>
+   template <typename DBODY>
+   static void run(const int N, DBODY &&d_body,
+                   const int X, const int Y, const int Z, const int G)
+   {
+      HipWrap2DLaunchBounds<MAX_THREADS_PER_BLOCK>(N, d_body, X, Y, Z);
+   }
+};
+
+template <>
+struct HipWrap<3, 0>
+{
+   template <typename DBODY>
    static void run(const int N, DBODY &&d_body,
                    const int X, const int Y, const int Z, const int G)
    {
@@ -753,11 +879,24 @@ struct HipWrap<3>
    }
 };
 
+template <int MAX_THREADS_PER_BLOCK>
+struct HipWrap<3, MAX_THREADS_PER_BLOCK>
+{
+   template <typename DBODY>
+   static void run(const int N, DBODY &&d_body,
+                   const int X, const int Y, const int Z, const int G)
+   {
+      HipWrap3DLaunchBounds<MAX_THREADS_PER_BLOCK>(N, d_body, X, Y, Z, G);
+   }
+};
+
 #endif // defined(MFEM_USE_HIP) && defined(__HIP__)
 
 
-/// The forall kernel body wrapper
-template <const int DIM, typename d_lambda, typename h_lambda>
+///////////////////////////////////////////////////////////////////////////////
+/// Forall host & device kernel dispatch
+template <int DIM, int MAX_THREADS_PER_BLOCK = 0,
+          typename d_lambda, typename h_lambda>
 inline void ForallWrap(const bool use_dev, const int N,
                        d_lambda &&d_body, h_lambda &&h_body,
                        const int X=0, const int Y=0, const int Z=0,
@@ -790,7 +929,7 @@ inline void ForallWrap(const bool use_dev, const int N,
    // If Backend::CUDA is allowed, use it
    if (Device::Allows(Backend::CUDA))
    {
-      return CuWrap<DIM>::run(N, d_body, X, Y, Z, G);
+      return CuWrap<DIM, MAX_THREADS_PER_BLOCK>::run(N, d_body, X, Y, Z, G);
    }
 #endif
 
@@ -798,7 +937,7 @@ inline void ForallWrap(const bool use_dev, const int N,
    // If Backend::HIP is allowed, use it
    if (Device::Allows(Backend::HIP))
    {
-      return HipWrap<DIM>::run(N, d_body, X, Y, Z, G);
+      return HipWrap<DIM, MAX_THREADS_PER_BLOCK>::run(N, d_body, X, Y, Z, G);
    }
 #endif
 
@@ -827,7 +966,9 @@ backend_cpu:
    for (int k = 0; k < N; k++) { h_body(k); }
 }
 
-template <const int DIM, typename lambda>
+///////////////////////////////////////////////////////////////////////////////
+/// Forall host & device kernel wrappers
+template <int DIM, typename lambda>
 inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
                        const int X=0, const int Y=0, const int Z=0,
                        const int G=0)
@@ -835,6 +976,16 @@ inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
    ForallWrap<DIM>(use_dev, N, body, body, X, Y, Z, G);
 }
 
+template <int DIM, int MAX_THREADS_PER_BLOCK, typename lambda>
+inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
+                       const int X=0, const int Y=0, const int Z=0,
+                       const int G=0)
+{
+   ForallWrap<DIM, MAX_THREADS_PER_BLOCK>(use_dev, N, body, body, X, Y, Z, G);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// forall interfaces
 template<typename lambda>
 inline void forall(int N, lambda &&body) { ForallWrap<1>(true, N, body); }
 
@@ -843,7 +994,7 @@ inline void forall(int Nx, int Ny, lambda &&body)
 {
    if (Device::Allows(Backend::DEVICE_MASK))
    {
-      forall(Nx * Ny, [=] MFEM_HOST_DEVICE(int idx)
+      mfem::forall(Nx * Ny, [=] MFEM_HOST_DEVICE(int idx)
       {
          int j = idx / Nx;
          int i = idx % Nx;
@@ -879,7 +1030,7 @@ inline void forall(int Nx, int Ny, int Nz, lambda &&body)
 {
    if (Device::Allows(Backend::DEVICE_MASK))
    {
-      forall(Nx * Ny * Nz, [=] MFEM_HOST_DEVICE(int idx)
+      mfem::forall(Nx * Ny * Nz, [=] MFEM_HOST_DEVICE(int idx)
       {
          int i = idx % Nx;
          int j = idx / Nx;
@@ -927,6 +1078,12 @@ inline void forall_2D(int N, int X, int Y, lambda &&body)
    ForallWrap<2>(true, N, body, X, Y, 1);
 }
 
+template<int MAX_THREADS_PER_BLOCK, typename lambda>
+inline void forall_2D(int N, int X, int Y, lambda &&body)
+{
+   ForallWrap<2, MAX_THREADS_PER_BLOCK>(true, N, body, X, Y, 1);
+}
+
 template<typename lambda>
 inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
 {
@@ -939,6 +1096,12 @@ inline void forall_3D(int N, int X, int Y, int Z, lambda &&body)
    ForallWrap<3>(true, N, body, X, Y, Z, 0);
 }
 
+template<int MAX_THREADS_PER_BLOCK, typename lambda>
+inline void forall_3D(int N, int X, int Y, int Z, lambda &&body)
+{
+   ForallWrap<3, MAX_THREADS_PER_BLOCK>(true, N, body, X, Y, Z, 0);
+}
+
 template<typename lambda>
 inline void forall_3D_grid(int N, int X, int Y, int Z, int G, lambda &&body)
 {

general/hip.hpp

@@ -20,9 +20,11 @@
 
 #if defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_USE_CUDA_OR_HIP
+constexpr bool mfem_use_gpu = true;
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
 #define MFEM_LAMBDA __host__ __device__
+#define MFEM_LAUNCH_BOUNDS __launch_bounds__
 // #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
 #define MFEM_DEVICE_SYNC MFEM_GPU_CHECK(hipDeviceSynchronize())
 #define MFEM_STREAM_SYNC MFEM_GPU_CHECK(hipStreamSynchronize(0))

linalg/CMakeLists.txt

@@ -55,6 +55,7 @@ list(APPEND HDRS
   dinvariants.hpp
   dtensor.hpp
   dual.hpp
+  eigensolver.hpp
   filteredsolver.hpp
   handle.hpp
   invariants.hpp

linalg/eigensolver.hpp

@@ -0,0 +1,203 @@
+// 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.
+
+/**
+ * @file eigensolver.hpp
+ *
+ * @brief This file contains a common interface for all eigensolver classes
+ */
+
+#ifndef MFEM_EIGENSOLVER
+#define MFEM_EIGENSOLVER
+
+#ifdef MFEM_HYPRE
+#include "hypre.hpp"
+#endif
+
+#ifdef MFEM_SLEPC
+#include "slepc.hpp"
+#endif
+
+namespace mfem
+{
+
+enum class EigenSolverType
+{
+   HYPRE,
+   SLEPC,
+   INVALID_TYPE
+};
+
+/// Provides base class for MFEM Eigensolvers
+class EigenSolverBase
+{
+public:
+   EigenSolverBase() {}
+
+   /// Destructor
+   virtual ~EigenSolverBase() = default;
+
+   /// Solves the eigenvalue problem
+   virtual void Solve() = 0;
+
+   /// Set the required number of modes
+   virtual void SetNumModes(int num_Modes)
+   {
+      numModes=num_Modes;
+   }
+
+   /// @brief  Set the operator to the eigenvalue problem
+   /// @param A - operator
+   virtual void SetOperator(Operator& A) = 0;
+
+   /// @brief Sets operators for the generalized eigenvalue problem
+   /// @param A - operator
+   /// @param M - mass matrix
+   virtual void SetOperator(Operator& A, Operator& M)
+   {
+      MFEM_ABORT("Generalized eigensolver is not supported!");
+   }
+
+   /// Optional method - sets preconditioner for the
+   /// eigenvalue solver.
+   virtual void SetPreconditioner(Solver& precond)
+   {
+      MFEM_ABORT("Preconditioner is not supported!");
+   }
+
+   /// Returns the converged eigenvalues
+   virtual void GetEigenvalues(Array<real_t>& eigen_vals) = 0;
+
+   /// Returns the vec_index eigenvector.
+   virtual void GetEigenvector(int vec_index, Vector& vector) = 0;
+
+   /// Returns the eigensolver type.
+   EigenSolverType GetSolverType() { return eigSolverType; }
+
+protected:
+   int numModes = 0;
+   EigenSolverType eigSolverType = EigenSolverType::INVALID_TYPE;
+};
+
+#ifdef MFEM_HYPRE
+class EigenSolverHypreLOBPCG : public EigenSolverBase
+{
+public:
+   EigenSolverHypreLOBPCG(MPI_Comm comm)
+   {
+      eigenSolver = std::make_unique<HypreLOBPCG>(comm);
+      eigSolverType = EigenSolverType::HYPRE;
+   }
+
+   ~EigenSolverHypreLOBPCG() {}
+
+   void Solve() override { eigenSolver->Solve(); }
+   void SetNumModes(int num_Modes) override
+   {
+      eigenSolver->SetNumModes(num_Modes);
+      numModes = num_Modes;
+   }
+
+   void SetOperator(Operator& A) override { eigenSolver->SetOperator(A); }
+
+   void SetOperator(Operator& A, Operator& M) override
+   {
+      eigenSolver->SetOperator(A);
+      eigenSolver->SetMassMatrix(M);
+   }
+
+   void SetPreconditioner(Solver& precond) override { eigenSolver->SetPreconditioner(precond); }
+   void GetEigenvalues(Array<real_t>& eigen_vals) override { eigenSolver->GetEigenvalues(eigen_vals); }
+   void GetEigenvector(int vec_index, Vector& vector) override
+   {
+      const HypreParVector& eigenvec = eigenSolver->GetEigenvector(vec_index);
+      vector = eigenvec;
+   }
+
+   void SetTol(real_t tol) { eigenSolver->SetTol(tol); }
+   void SetRelTol(real_t rel_tol) { eigenSolver->SetRelTol(rel_tol); }
+   void SetMaxIter(int max_iter) { eigenSolver->SetMaxIter(max_iter); }
+   void SetPrintLevel(int logging) { eigenSolver->SetPrintLevel(logging); }
+   void SetRandomSeed(int seed) { eigenSolver->SetRandomSeed(seed); }
+   void SetPrecondUsageMode(int usage_mode) { eigenSolver->SetPrecondUsageMode(usage_mode); }
+
+private:
+   std::unique_ptr<HypreLOBPCG> eigenSolver = nullptr;
+};
+#endif
+
+#ifdef MFEM_SLEPC
+class EigenSolverSlepc : public EigenSolverBase
+{
+public:
+   EigenSolverSlepc(MPI_Comm comm)
+   {
+      eigSolverType = EigenSolverType::SLEPC;
+      eigenSolver = std::make_unique<SlepcEigenSolver>(comm);
+
+      eigenSolver->SetWhichEigenpairs(SlepcEigenSolver::TARGET_REAL);
+      eigenSolver->SetTarget(0.0);
+      eigenSolver->SetSpectralTransformation(SlepcEigenSolver::SHIFT_INVERT);
+   }
+
+   ~EigenSolverSlepc() {}
+
+   void Solve() override { eigenSolver->Solve(); }
+   void SetNumModes(int num_Modes) override
+   {
+      eigenSolver->SetNumModes(num_Modes);
+      numModes = num_Modes;
+   }
+   /// @brief  Set the operator to the slepc eigenvalue problem. This method deep copies data to create a PetscParMatrix
+   /// @param A - operator, must be of type HypreParMatrix.
+   void SetOperator(Operator& A) override
+   {
+      petscMatA = std::make_unique<PetscParMatrix>
+                  (dynamic_cast<HypreParMatrix*>(&A));
+      eigenSolver->SetOperator(*petscMatA);
+   }
+   /// @brief  Set the operators to the slepc eigenvalue problem. This method deep copies data to create a PetscParMatrix
+   /// @param A - operator, must be of type HypreParMatrix.
+   /// @param M - operator, must be of type HypreParMatrix.
+   void SetOperator(Operator& A, Operator& M) override
+   {
+      petscMatA = std::make_unique<PetscParMatrix>
+                  (dynamic_cast<const HypreParMatrix*>(&A));
+      petscMatM = std::make_unique<PetscParMatrix>
+                  (dynamic_cast<const HypreParMatrix*>(&M));
+
+      eigenSolver->SetOperators(*petscMatA, *petscMatM);
+   }
+   void SetPreconditioner([[maybe_unused]] Solver& precond) override {}
+   void GetEigenvalues(Array<real_t>& eigen_vals) override
+   {
+      eigen_vals.SetSize(numModes);
+      for (int ik = 0; ik < numModes; ik++)
+      {
+         eigenSolver->GetEigenvalue(static_cast<unsigned int>(ik), eigen_vals[ik]);
+      }
+   }
+   void GetEigenvector( int vec_index, Vector& vector) override
+   { eigenSolver->GetEigenvector(vec_index, vector); }
+
+   void SetTol(real_t tol) { eigenSolver->SetTol(tol); }
+   void SetMaxIter(int max_iter) { eigenSolver->SetMaxIter(max_iter); }
+
+private:
+   std::unique_ptr<SlepcEigenSolver> eigenSolver = nullptr;
+   std::unique_ptr<PetscParMatrix> petscMatA = nullptr;
+   std::unique_ptr<PetscParMatrix> petscMatM = nullptr;
+};
+#endif
+
+}  // namespace mfem
+
+#endif

linalg/superlu.cpp

@@ -363,14 +363,19 @@ void SuperLUSolver::Init(MPI_Comm comm)
    // Set default options:
    //    options.Fact = DOFACT;
    //    options.Equil = YES;
+   //    options.ParSymbFact = NO;
    //    options.ColPerm = METIS_AT_PLUS_A;
    //    options.RowPerm = LargeDiag_MC64;
    //    options.ReplaceTinyPivot = NO;
-   //    options.Trans = NOTRANS;
    //    options.IterRefine = SLU_DOUBLE;
+   //    options.Trans = NOTRANS;
    //    options.SolveInitialized = NO;
    //    options.RefineInitialized = NO;
    //    options.PrintStat = YES;
+   //    options.lookahead_etree = NO;
+   //    options.num_lookaheads = 10;
+   //    options.superlu_acc_offload = 1;
+   //    options.SymPattern = NO;
    superlu_dist_options_t *options = (superlu_dist_options_t *)optionsPtr_;
    set_default_options_dist(options);
 #if SUPERLU_DIST_MAJOR_VERSION > 7 || \
@@ -472,6 +477,12 @@ void SuperLUSolver::SetFact(superlu::Fact fact)
    options->Fact = opt;
 }
 
+void SuperLUSolver::SetDeviceOffload(bool offload)
+{
+   superlu_dist_options_t *options = (superlu_dist_options_t *)optionsPtr_;
+   options->superlu_acc_offload = offload;
+}
+
 void SuperLUSolver::SetOperator(const Operator &op)
 {
    // Verify that we have a compatible operator

linalg/superlu.hpp

@@ -250,7 +250,8 @@ public:
        work (default false) */
    void SetSymmetricPattern(bool sym);
 
-   /** @brief Specify whether to perform parallel symbolic factorization.
+   /** @brief Specify whether to perform parallel symbolic factorization
+       (default false)
        @note If true SuperLU will use superlu::PARMETIS for the Column
        Permutation regardless of the setting */
    void SetParSymbFact(bool par);
@@ -263,6 +264,10 @@ public:
        superlu::FACTORED*/
    void SetFact(superlu::Fact fact);
 
+   /** @brief Specify whether to offload numerical factorization onto the device
+       (default true if SuperLU_DIST has been compiled with GPU support) */
+   void SetDeviceOffload(bool offload);
+
    // Processor grid for SuperLU_DIST.
    const int nprow_, npcol_, npdep_;
 

tests/benchmarks/CMakeLists.txt

@@ -31,11 +31,6 @@ function(add_benchmark name)
        set_property(SOURCE ${${NAME}_BENCH_SRCS} PROPERTY LANGUAGE CUDA)
     endif(MFEM_USE_CUDA)
 
-    if (MFEM_USE_HIP)
-        set_property(SOURCE ${${NAME}_BENCH_SRCS} PROPERTY LANGUAGE 
-                     HIP_SOURCE_PROPERTY_FORMAT TRUE)
-    endif(MFEM_USE_HIP)
-
     add_executable(bench_${name} ${${NAME}_BENCH_SRCS})
     target_link_libraries(bench_${name} mfem pthread)
     add_dependencies(${MFEM_ALL_BENCHMARKS_TARGET_NAME} bench_${name})

tests/benchmarks/bench_ceed.cpp

@@ -8,23 +8,89 @@
 // 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.
+//
+//
+//  This benchmark contains the implementation of the CEED's bake-off problems:
+//  high-order kernels/benchmarks designed to test and compare the performance
+//  of high-order codes.
+//
+//  See: https://ceed.exascaleproject.org/bps
 
-#include "bench.hpp"
+#include "bench.hpp" // IWYU pragma: keep
 
 #ifdef MFEM_USE_BENCHMARK
 
-/*
-  This benchmark contains the implementation of the CEED's bake-off problems:
-  high-order kernels/benchmarks designed to test and compare the performance
-  of high-order codes.
+#include <cassert>
+#include <string>
 
-  See: ceed.exascaleproject.org/bps and github.com/CEED/benchmarks
-*/
-template <int VDIM, bool GLL>
+#include "fem/qinterp/det.hpp" // IWYU pragma: keep
+#include "fem/qinterp/grad.hpp" // IWYU pragma: keep
+#include "fem/integ/lininteg_domain_kernels.hpp" // IWYU pragma: keep
+#include "fem/integ/bilininteg_vecdiffusion_pa.hpp" // IWYU pragma: keep
+
+// Custom benchmark arguments generator
+static void CustomArguments(bmi::Benchmark *b) noexcept
+{
+   constexpr int MAX_NDOFS = 16 * 1024 * (mfem_use_gpu ? 1024 : 8);
+
+   const auto orders = { 7, 6, 5, 4, 3, 2, 1 };
+
+   constexpr auto ndofs = [](int n) constexpr noexcept -> int
+   {
+      return (n + 1) * (n + 1) * (n + 1);
+   };
+
+   constexpr auto inc = [](int n) constexpr noexcept -> int
+   {
+      return n < 160 ?  4 : n < 240 ?  8 : n < 320 ? 16 : 32;
+   };
+
+   for (auto p : orders)
+   {
+      for (int n = 16; ndofs(n) <= MAX_NDOFS; n += inc(n))
+      {
+         b->Args({p, n});
+      }
+   }
+}
+
+// Register kernel specializations used in the benchmarks
+static void AddKernelSpecializations()
+{
+   using DET = QuadratureInterpolator::DetKernels;
+   DET::Specialization<3, 3, 2, 2>::Add();
+   DET::Specialization<3, 3, 2, 3>::Add();
+   DET::Specialization<3, 3, 2, 5>::Add();
+   DET::Specialization<3, 3, 2, 6>::Add();
+   DET::Specialization<3, 3, 5, 5>::Add();
+   // Others might exceed memory limits
+
+   using GRAD = QuadratureInterpolator::GradKernels;
+   GRAD::Specialization<3, QVectorLayout::byNODES, false, 3, 2, 2>::Add();
+   GRAD::Specialization<3, QVectorLayout::byNODES, false, 3, 2, 7>::Add();
+   GRAD::Specialization<3, QVectorLayout::byNODES, false, 3, 2, 8>::Add();
+   GRAD::Specialization<3, QVectorLayout::byNODES, false, 3, 2, 9>::Add();
+
+   using LIN = DomainLFIntegrator::AssembleKernels;
+   LIN::Specialization<3, 7, 7>::Add();
+   LIN::Specialization<3, 6, 6>::Add();
+   LIN::Specialization<3, 8, 8>::Add();
+
+   using VDIFF = VectorDiffusionIntegrator::ApplyPAKernels;
+   VDIFF::Specialization<3, 3, 3, 3>::Add();
+   VDIFF::Specialization<3, 3, 4, 4>::Add();
+   VDIFF::Specialization<3, 3, 5, 5>::Add();
+   VDIFF::Specialization<3, 3, 6, 6>::Add();
+   VDIFF::Specialization<3, 3, 7, 7>::Add();
+   VDIFF::Specialization<3, 3, 8, 8>::Add();
+}
+
+// Bake-off base class
+template <int BFI, int VDIM, bool GLL>
 struct BakeOff
 {
-   static constexpr int DIM = 3;
-   const int N, p, q;
+   inline static constexpr int DIM = 3;
+   const int p, c, q, n, nx, ny, nz;
    Mesh mesh;
    H1_FECollection fec;
    FiniteElementSpace fes;
@@ -38,12 +104,15 @@ struct BakeOff
    GridFunction x, y;
    BilinearForm a;
    double mdofs{};
+   BilinearFormIntegrator *bfi;
 
-   BakeOff(int p):
-      N(Device::IsEnabled() ? 32 : 4),
-      p(p),
-      q(2 * p + (GLL ? -1 : 3)),
-      mesh(Mesh::MakeCartesian3D(N, N, N, Element::HEXAHEDRON)),
+   BakeOff(int p, int side):
+      p(p), c(side), q(2 * p + (GLL ? -1 : 3)),
+      n((assert(c >= p), c / p)),
+      nx(n + (p * (n + 1) * p * n * p * n < c * c * c ? 1 : 0)),
+      ny(n + (p * (n + 1) * p * (n + 1) * p * n < c * c * c ? 1 : 0)),
+      nz(n),
+      mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON)),
       fec(p, DIM, BasisType::GaussLobatto),
       fes(&mesh, &fec, VDIM, VDIM == 3 ? Ordering::byVDIM : Ordering::byNODES),
       geom_type(mesh.GetTypicalElementGeometry()),
@@ -58,22 +127,41 @@ struct BakeOff
       a(&fes)
    {
       x = 0.0;
+      if constexpr (BFI == 1)
+      {
+         bfi = new MassIntegrator(one, ir);
+      }
+      else if constexpr (BFI == 2)
+      {
+         bfi = new VectorMassIntegrator(one, ir);
+      }
+      else if constexpr (BFI == 3 || BFI == 5)
+      {
+         bfi = new DiffusionIntegrator(one, ir);
+      }
+      else if constexpr (BFI == 4 || BFI == 6)
+      {
+         bfi = new VectorDiffusionIntegrator(one, ir);
+      }
+      else
+      {
+         static_assert(BFI >= 1 && BFI <= 6, "Invalid BilinearFormIntegrator");
+      }
+      a.AddDomainIntegrator(bfi);
    }
 
    virtual void benchmark() = 0;
 
-   double SumMdofs() const { return mdofs; }
+   [[nodiscard]] double SumMdofs() const noexcept { return mdofs; }
 
-   double MDofs() const { return 1e-6 * dofs; }
+   [[nodiscard]] double MDofs() const noexcept { return 1e-6 * dofs; }
 };
 
-/// Bake-off Problems (BPs)
-template <typename BFI, int VDIM, bool GLL>
-struct Problem : public BakeOff<VDIM, GLL>
+// Bake-off Problems (BPs)
+template <int BFI, int VDIM, bool GLL>
+struct BP : public BakeOff<BFI, VDIM, GLL>
 {
-   const double rtol = 1e-12;
-   const int max_it = 32;
-   const int print_lvl = -1;
+   const int max_it = 32, print_lvl = -1;
 
    Array<int> ess_tdof_list;
    Array<int> ess_bdr;
@@ -82,44 +170,56 @@ struct Problem : public BakeOff<VDIM, GLL>
    Vector B, X;
    CGSolver cg;
 
-   using BakeOff<VDIM, GLL>::a;
-   using BakeOff<VDIM, GLL>::ir;
-   using BakeOff<VDIM, GLL>::one;
-   using BakeOff<VDIM, GLL>::mesh;
-   using BakeOff<VDIM, GLL>::fes;
-   using BakeOff<VDIM, GLL>::x;
-   using BakeOff<VDIM, GLL>::y;
-   using BakeOff<VDIM, GLL>::mdofs;
+   using base = BakeOff<BFI, VDIM, GLL>;
+   using base::a;
+   using base::ir;
+   using base::one;
+   using base::mesh;
+   using base::fes;
+   using base::x;
+   using base::y;
+   using base::mdofs;
+   using base::unit_vec;
+   using base::bfi;
 
-   Problem(int order):
-      BakeOff<VDIM, GLL>(order),
+   BP(int p, int side) noexcept: base(p, side),
       ess_bdr(mesh.bdr_attributes.Max()),
       b(&fes)
    {
       ess_bdr = 1;
       fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-      if (VDIM == 1)
+
+      if constexpr (VDIM == 1)
       {
-         b.AddDomainIntegrator(new DomainLFIntegrator(this->one));
+         b.AddDomainIntegrator(new DomainLFIntegrator(one));
       }
       else
       {
-         b.AddDomainIntegrator(new VectorDomainLFIntegrator(this->unit_vec));
+         b.AddDomainIntegrator(new VectorDomainLFIntegrator(unit_vec));
       }
       b.UseFastAssembly(true);
       b.Assemble();
 
       a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-      a.AddDomainIntegrator(new BFI(one, ir));
       a.Assemble();
       a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
 
-      cg.SetRelTol(rtol);
       cg.SetOperator(*A);
+      cg.SetAbsTol(0.0);
+      cg.iterative_mode = false;
+      {
+         cg.SetPrintLevel(-1);
+         cg.SetMaxIter(1000);
+         cg.SetRelTol(1e-8);
+         cg.Mult(B, X);
+         MFEM_VERIFY(cg.GetConverged(), "CG solver did not converge!");
+      }
+      cg.SetRelTol(0.0);
       cg.SetMaxIter(max_it);
       cg.SetPrintLevel(print_lvl);
-      cg.iterative_mode = false;
-      MFEM_DEVICE_SYNC;
+
+      benchmark();
+      mdofs = 0.0;
    }
 
    void benchmark() override
@@ -130,104 +230,115 @@ struct Problem : public BakeOff<VDIM, GLL>
    }
 };
 
-/// Bake-off Problems (BPs)
-#define BakeOff_Problem(i, Kernel, VDIM, p_eq_q)                     \
-   static void BP##i(bm::State &state)                               \
-   {                                                                 \
-      Problem<Kernel##Integrator, VDIM, p_eq_q> ker(state.range(0)); \
-      while (state.KeepRunning()) { ker.benchmark(); }               \
-      state.counters["MDof/s"] =                                     \
-         bm::Counter(ker.SumMdofs(), bm::Counter::kIsRate);          \
-   }                                                                 \
-   BENCHMARK(BP##i)->DenseRange(1, 6)->Unit(bm::kMillisecond);
-
-/// BP1: scalar PCG with mass matrix, q=p+2
-BakeOff_Problem(1, Mass, 1, false)
-
-/// BP2: vector PCG with mass matrix, q=p+2
-BakeOff_Problem(2, VectorMass, 3, false)
-
-/// BP3: scalar PCG with stiffness matrix, q=p+2
-BakeOff_Problem(3, Diffusion, 1, false)
-
-/// BP4: vector PCG with stiffness matrix, q=p+2
-BakeOff_Problem(4, VectorDiffusion, 3, false)
-
-/// BP5: scalar PCG with stiffness matrix, q=p+1
-BakeOff_Problem(5, Diffusion, 1, true)
-
-/// BP6: vector PCG with stiffness matrix, q=p+1
-BakeOff_Problem(6, VectorDiffusion, 3, true)
-
-/// Bake-off Kernels (BKs)
-template <typename BFI, int VDIM, bool GLL>
-struct Kernel : public BakeOff<VDIM, GLL>
+// Bake-off Kernels (BKs)
+template <int BFI, int VDIM, bool GLL>
+struct BK : public BakeOff<BFI, VDIM, GLL>
 {
-   using BakeOff<VDIM, GLL>::a;
-   using BakeOff<VDIM, GLL>::ir;
-   using BakeOff<VDIM, GLL>::one;
-   using BakeOff<VDIM, GLL>::fes;
-   using BakeOff<VDIM, GLL>::x;
-   using BakeOff<VDIM, GLL>::y;
-   using BakeOff<VDIM, GLL>::mdofs;
+   Vector xe, ye;
 
-   Kernel(int order): BakeOff<VDIM, GLL>(order)
+   using base = BakeOff<BFI, VDIM, GLL>;
+   using base::ir;
+   using base::one;
+   using base::bfi;
+   using base::fes;
+   using base::mdofs;
+
+   BK(int order, int side) noexcept: base(order, side)
    {
-      x.Randomize(1);
-      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-      a.AddDomainIntegrator(new BFI(one, ir));
-      a.Assemble();
-      a.Mult(x, y);
-      MFEM_DEVICE_SYNC;
+      bfi->AssemblePA(fes);
+
+      const Table &el2dof = fes.GetElementToDofTable();
+      const int e_size = el2dof.Size_of_connections()*fes.GetVDim();
+      const auto R = fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      MFEM_VERIFY(e_size == R->Height(), "Input/Output E-vector size mismatch!");
+
+      xe.SetSize(R->Height());
+      ye.SetSize(R->Height());
+      xe.UseDevice(true);
+      ye.UseDevice(true);
+
+      xe.Randomize(1);
+      xe.Read();
+      ye = 0.0;
+
+      benchmark();
+      mdofs = 0.0;
    }
 
    void benchmark() override
    {
-      a.Mult(x, y);
+      bfi->AddMultPA(xe, ye);
       MFEM_DEVICE_SYNC;
       mdofs += this->MDofs();
    }
 };
 
-/// Generic CEED BKi
-#define BakeOff_Kernel(i, KER, VDIM, GLL)                      \
-   static void BK##i(bm::State &state)                         \
-   {                                                           \
-      Kernel<KER##Integrator, VDIM, GLL> ker(state.range(0));  \
-      while (state.KeepRunning()) { ker.benchmark(); }         \
-      state.counters["MDof/s"] =                               \
-         bm::Counter(ker.SumMdofs(), bm::Counter::kIsRate);    \
-   }                                                           \
-   BENCHMARK(BK##i)->DenseRange(1, 6)->Unit(bm::kMillisecond);
+// Benchmarks
+template <typename T>
+static void Benchmark(bm::State& state) noexcept
+{
+   T run(state.range(0), state.range(1));
+   while (state.KeepRunning()) { run.benchmark(); }
+   state.counters["Dofs"] = bm::Counter(run.dofs);
+   state.counters["MDof/s"] = bm::Counter(run.SumMdofs(), bm::Counter::kIsRate);
+   state.counters["Order"] = bm::Counter(state.range(0));
+}
 
-/// BK1: scalar E-vector-to-E-vector evaluation of mass matrix, q=p+2
-BakeOff_Kernel(1, Mass, 1, false)
+#define REGISTER(PK, BFI, VDIM, GLL) \
+   BENCHMARK_TEMPLATE(Benchmark, PK<BFI, VDIM, GLL>) \
+   ->Name(#PK #BFI)->Apply(CustomArguments)->Unit(bm::kMillisecond)
 
-/// BK2: vector E-vector-to-E-vector evaluation of mass matrix, q=p+2
-BakeOff_Kernel(2, VectorMass, 3, false)
+// BP1: scalar PCG with mass matrix, q=p+2
+REGISTER(BP, 1, 1, false);
 
-/// BK3: scalar E-vector-to-E-vector evaluation of stiffness matrix, q=p+2
-BakeOff_Kernel(3, Diffusion, 1, false)
+// BP2: vector PCG with mass matrix, q=p+2
+REGISTER(BP, 2, 3, false);
 
-/// BK4: vector E-vector-to-E-vector evaluation of stiffness matrix, q=p+2
-BakeOff_Kernel(4, VectorDiffusion, 3, false)
+// BP3: scalar PCG with stiffness matrix, q=p+2
+REGISTER(BP, 3, 1, false);
 
-/// BK5: scalar E-vector-to-E-vector evaluation of stiffness matrix, q=p+1
-BakeOff_Kernel(5, Diffusion, 1, true)
+// BP4: vector PCG with stiffness matrix, q=p+2
+REGISTER(BP, 4, 3, false);
 
-/// BK6: vector E-vector-to-E-vector evaluation of stiffness matrix, q=p+1
-BakeOff_Kernel(6, VectorDiffusion, 3, true)
+// BP5: scalar PCG with stiffness matrix, q=p+1
+REGISTER(BP, 5, 1, true);
+
+// BP6: vector PCG with stiffness matrix, q=p+1
+REGISTER(BP, 6, 3, true);
+
+// BK1: scalar E-vector-to-E-vector evaluation of mass matrix, q=p+2
+REGISTER(BK, 1, 1, false);
+
+// BK2: vector E-vector-to-E-vector evaluation of mass matrix, q=p+2
+REGISTER(BK, 2, 3, false);
+
+// BK3: scalar E-vector-to-E-vector evaluation of stiffness matrix, q=p+2
+REGISTER(BK, 3, 1, false);
+
+// BK4: vector E-vector-to-E-vector evaluation of stiffness matrix, q=p+2
+REGISTER(BK, 4, 3, false);
+
+// BK5: scalar E-vector-to-E-vector evaluation of stiffness matrix, q=p+1
+REGISTER(BK, 5, 1, true);
+
+// BK6: vector E-vector-to-E-vector evaluation of stiffness matrix, q=p+1
+REGISTER(BK, 6, 3, true);
 
 /**
- * @brief main entry point
- * --benchmark_filter=BK1/6
- * --benchmark_context=device=cpu
+ * @brief CEED Bake-off Problems main entry point
+ * Command line options:
+ *    --benchmark_context=device=gpu
+ *    --benchmark_filter=BP1
+ *    --benchmark_out_format=csv
+ *    --benchmark_out=bp1.csv
  */
 int main(int argc, char *argv[])
 {
    bm::ConsoleReporter CR;
    bm::Initialize(&argc, argv);
 
+   AddKernelSpecializations();
+
    // Device setup, cpu by default
    std::string device_config = "cpu";
    auto global_context = bmi::GetGlobalContext();
@@ -240,12 +351,16 @@ int main(int argc, char *argv[])
          device_config = device->second;
       }
    }
+
    Device device(device_config.c_str());
    device.Print();
 
-   if (bm::ReportUnrecognizedArguments(argc, argv)) { return 1; }
+   if (bm::ReportUnrecognizedArguments(argc, argv)) { return EXIT_FAILURE; }
+
    bm::RunSpecifiedBenchmarks(&CR);
-   return 0;
+   bm::Shutdown();
+
+   return EXIT_SUCCESS;
 }
 
 #endif // MFEM_USE_BENCHMARK

tests/unit/fem/test_transfer.cpp

@@ -425,95 +425,6 @@ TEST_CASE("Variable Order True Transfer", "[Transfer][VariableOrder]")
    delete c_fec;
 }
 
-TEST_CASE("H1 L2 transfer with consistent mass", "[Transfer]")
-{
-   auto vectorspace = GENERATE(VecSpace::H1, VecSpace::VectorH1);
-   dimension = GENERATE(2, 3);
-
-   const int order = 2;
-   const int ne = 2;
-   const int vdim = (vectorspace == VecSpace::VectorH1) ? dimension : 1;
-
-   CAPTURE(VecSpaceName(vectorspace), dimension, order);
-
-   Mesh mesh;
-   if (dimension == 2)
-   {
-      mesh = Mesh::MakeCartesian2D(ne, ne, Element::QUADRILATERAL,
-                                   1, 1.0, 1.0);
-   }
-   else
-   {
-      mesh = Mesh::MakeCartesian3D(ne, ne, ne, Element::HEXAHEDRON,
-                                   1.0, 1.0, 1.0);
-   }
-
-   Mesh fineMesh(mesh);
-   fineMesh.UniformRefinement();
-
-   H1_FECollection fec(order, dimension);
-   FiniteElementSpace c_fespace(&mesh, &fec, vdim);
-   FiniteElementSpace f_fespace(&fineMesh, &fec, vdim);
-
-   L2ProjectionGridTransfer transfer(c_fespace, f_fespace);
-   transfer.UseConsistentMass();
-   const Operator &R = transfer.ForwardOperator();
-
-   GridFunction X(&c_fespace);
-   GridFunction Y(&f_fespace);
-   GridFunction Y_ref(&f_fespace);
-   coeff_order = 1;
-
-   LinearForm rhs(&f_fespace);
-   BilinearForm mass(&f_fespace);
-   FunctionCoefficient funcCoeff(&coeff);
-   VectorFunctionCoefficient vecCoeff(dimension, &vectorcoeff);
-   if (vectorspace == VecSpace::H1)
-   {
-      X.ProjectCoefficient(funcCoeff);
-      rhs.AddDomainIntegrator(new DomainLFIntegrator(funcCoeff));
-      mass.AddDomainIntegrator(new MassIntegrator);
-   }
-   else
-   {
-      X.ProjectCoefficient(vecCoeff);
-      rhs.AddDomainIntegrator(new VectorDomainLFIntegrator(vecCoeff));
-      mass.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-
-   rhs.Assemble();
-   mass.Assemble();
-   SparseMatrix M;
-   Array<int> empty;
-   mass.FormSystemMatrix(empty, M);
-
-   GSSmoother M_prec(M);
-   Y_ref = 0.0;
-   PCG(M, M_prec, rhs, Y_ref, 0, 500, 1e-24, 0.0);
-
-   Y = 0.0;
-   R.Mult(X, Y);
-   Y -= Y_ref;
-   REQUIRE(Y.Norml2() < 1e-11 * Y_ref.Norml2());
-
-   Vector x(c_fespace.GetVSize());
-   Vector y(f_fespace.GetVSize());
-   Vector Ry(f_fespace.GetVSize());
-   Vector Rtx(c_fespace.GetVSize());
-   x.Randomize(1);
-   y.Randomize(2);
-
-   R.Mult(x, Ry);
-   R.MultTranspose(y, Rtx);
-
-   const real_t ip1 = InnerProduct(Ry, y);
-   const real_t ip2 = InnerProduct(x, Rtx);
-   REQUIRE(std::abs(ip1 - ip2) <
-           1e-10 * std::max(std::abs(ip1), std::abs(ip2)));
-
-   REQUIRE_FALSE(transfer.SupportsBackwardsOperator());
-}
-
 TEST_CASE("Restriction Transpose Operator")
 {
    int order = GENERATE(1, 2);
@@ -673,57 +584,4 @@ TEST_CASE("Parallel Transfer", "[Transfer][Parallel]")
    delete pmesh;
 }
 
-TEST_CASE("Parallel H1 L2 transfer with consistent mass",
-          "[Transfer][Parallel]")
-{
-   dimension = GENERATE(2, 3);
-
-   const int order = 2;
-   const int ne = 2;
-   const int vdim = 1;
-
-   CAPTURE(dimension, order);
-
-   Mesh mesh;
-   if (dimension == 2)
-   {
-      mesh = Mesh::MakeCartesian2D(ne, ne, Element::QUADRILATERAL,
-                                   1, 1.0, 1.0);
-   }
-   else
-   {
-      mesh = Mesh::MakeCartesian3D(ne, ne, ne, Element::HEXAHEDRON,
-                                   1.0, 1.0, 1.0);
-   }
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   ParMesh pfineMesh(MPI_COMM_WORLD, mesh);
-   pfineMesh.UniformRefinement();
-
-   H1_FECollection fec(order, dimension);
-   ParFiniteElementSpace c_fespace(&pmesh, &fec, vdim);
-   ParFiniteElementSpace f_fespace(&pfineMesh, &fec, vdim);
-
-   L2ProjectionGridTransfer transfer(c_fespace, f_fespace);
-   transfer.UseConsistentMass();
-   const Operator &R = transfer.TrueForwardOperator();
-
-   Vector x(c_fespace.GetTrueVSize());
-   Vector y(f_fespace.GetTrueVSize());
-   Vector Rx(f_fespace.GetTrueVSize());
-   Vector Rty(c_fespace.GetTrueVSize());
-   x.Randomize(1);
-   y.Randomize(2);
-
-   R.Mult(x, Rx);
-   R.MultTranspose(y, Rty);
-
-   const real_t ip1 = InnerProduct(MPI_COMM_WORLD, Rx, y);
-   const real_t ip2 = InnerProduct(MPI_COMM_WORLD, x, Rty);
-   REQUIRE(std::abs(ip1 - ip2) <
-           1e-10 * std::max(std::abs(ip1), std::abs(ip2)));
-
-   REQUIRE_FALSE(transfer.SupportsBackwardsOperator());
-}
-
 #endif