fix year

merge master and resolve conflicts
Merge pull request #4960 from mfem/fix-warnings-algoim
2025-08-07 11:18:41 -07:00 · 2025-08-07 10:59:40 -07:00 · 2025-08-05 12:05:40 -07:00 · 2025-08-04 07:44:47 -07:00 · 2025-08-02 09:27:45 -07:00 · 2025-07-29 12:49:00 -07:00
30 changed files with 1986 additions and 105 deletions
@@ -211,7 +211,7 @@ miniapps/electromagnetics/joule
 miniapps/electromagnetics/Volta-AMR*
 miniapps/electromagnetics/Tesla-AMR*
 miniapps/electromagnetics/Maxwell-Parallel*
-miniapps/electromagnetics/Joule_*
+miniapps/electromagnetics/Joule_[0-9]*

 miniapps/gslib/field-diff
 miniapps/gslib/field-interp
@@ -267,9 +267,9 @@ miniapps/meshing/bounding-box*
 miniapps/meshing/jacobian-determinant*

 miniapps/mtop/parheat
-miniapps/mtop/ParHeat*
+miniapps/mtop/ParHeat/*
 miniapps/mtop/seqheat
-miniapps/mtop/SeqHeat*
+miniapps/mtop/SeqHeat/*

 miniapps/autodiff/paradiff
 miniapps/autodiff/seqadiff
@@ -277,7 +277,7 @@ miniapps/autodiff/seqtest
 miniapps/autodiff/par_example
 miniapps/autodiff/seq_example
 miniapps/autodiff/seq_test
-miniapps/autodiff/Exampl*
+miniapps/autodiff/Example/*

 miniapps/navier/navier_mms
 miniapps/navier/navier_kovasznay
@@ -52,4 +52,4 @@ variables:
    - echo ${JOBID}
    - echo ${MFEM_DATA_DIR}
    - echo ${SPEC}
-    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) --reservation=ci -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) --reservation=ci -t 60 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
@@ -46,6 +46,11 @@ GPU computing
  where host execution is always needed (e.g. when the DOFs array is small).
 - Introduced MFEM_FOREACH_THREAD_DIRECT, which directly maps loop tasks to GPU
  threads, assigning one task per thread.
+- Implemented a GPU-accelerated matrix-free AMR derefinement `GridFunction`
+  update operator. This supports mixed geometry meshes and variable order
+  spaces, and is the default derefinement operator constructed by
+  `FiniteElementSpace::Update` and `ParFiniteElementSpace::Update`.
+  The operator requires `FiniteElementSpace::Nonconforming() == true`.

 New and updated examples and miniapps
 -------------------------------------
@@ -598,14 +598,20 @@ set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
    ALGOIM ENZYME)

-# Add all *_FOUND libraries in the variable TPL_LIBRARIES.
+# Add all created targets and *_FOUND libraries in the variables TPL_TARGETS and
+# TPL_LIBRARIES, respectively.
+set(TPL_TARGETS)
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")
 foreach(TPL IN LISTS MFEM_TPLS)
-  if (${TPL}_FOUND)
+  if (${TPL}_FOUND OR TARGET ${TPL})
    message(STATUS "MFEM: using package ${TPL}")
-    list(APPEND TPL_LIBRARIES ${${TPL}_LIBRARIES})
-    list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
+    if (TARGET ${TPL})
+      list(APPEND TPL_TARGETS ${TPL})
+    else()
+      list(APPEND TPL_LIBRARIES ${${TPL}_LIBRARIES})
+      list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
+    endif()
  endif()
 endforeach(TPL)
 list(REVERSE TPL_LIBRARIES)
@@ -680,7 +686,10 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})
 # Declaring the library
 mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
-target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES})
+target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
+if (TPL_TARGETS)
+  add_dependencies(mfem ${TPL_TARGETS})
+endif()
 if (MINGW)
  target_link_libraries(mfem PRIVATE ws2_32)
 endif()
@@ -121,6 +121,11 @@ Parallel build:
   make -j 4
 (For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

+Parallel build with fetching of hypre and METIS:
+   mkdir <mfem-buil-dir> ; cd <mfem-build-dir>
+   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES -DFETCH_TPLS=YES
+   make -j 4
+
 CUDA build:
   (this build requires CMake 3.17 or newer)
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
@@ -1074,6 +1079,9 @@ The following options are CMake specific:
 MFEM_ENABLE_TESTING  - Enable the ctest framework for testing.
 MFEM_ENABLE_EXAMPLES - Build all of the examples by default.
 MFEM_ENABLE_MINIAPPS - Build all of the miniapps by default.
+FETCH_TPLS           - Enable fetching of all supported third-party libraries.
+HYPRE_FETCH          - Enable fetching of hypre.
+METIS_FETCH          - Enable fetching of metis.

 External libraries (CMake):
 ---------------------------
@@ -1135,6 +1143,12 @@ The following built-in CMake packages are also used:
      set the <LIBNAME>_LIBRARIES option directly; the configuration option
      <LIBNAME>_DIR is not supported.

+The MFEM CMake build system also provides fetching (automated building) for the
+packages/libraries listed below. Note that when fetching is enabled, any related
+auto-detection functionality is disabled.
+
+ - HYPRE
+ - METIS

 Building without GNU make or CMake
 ==================================
@@ -9,15 +9,18 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

-# Defines the following variables:
+# Defines the following variables if fetching of TPLs is disabled (default):
 #   - HYPRE_FOUND
 #   - HYPRE_LIBRARIES
 #   - HYPRE_INCLUDE_DIRS
 #   - HYPRE_VERSION
 #   - HYPRE_USING_CUDA (internal)
 #   - HYPRE_USING_HIP (internal)
+# otherwise, the following are defined:
+#   - HYPRE (imported library target)
+#   - HYPRE_VERSION (cache variable)

-if (HYPRE_FOUND)
+if (HYPRE_FOUND OR TARGET HYPRE)
  if (HYPRE_USING_CUDA)
    find_package(CUDAToolkit REQUIRED)
  endif()
@@ -33,6 +36,55 @@ if (HYPRE_FOUND)
  endif()
 endif()

+if (HYPRE_FETCH OR FETCH_TPLS)
+  set(HYPRE_FETCH_VERSION 2.33.0)
+  add_library(HYPRE STATIC IMPORTED)
+  # set options and associated dependencies
+  set(CMAKE_OPTIONS)
+  list(APPEND CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
+  if (MFEM_USE_CUDA)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_CUDA:BOOL=ON)
+    find_package(CUDAToolkit REQUIRED)
+    target_link_libraries(HYPRE INTERFACE CUDA::cusparse CUDA::curand CUDA::cublas)
+  elseif (MFEM_USE_HIP)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_HIP:BOOL=ON)
+    find_package(rocsparse REQUIRED)
+    find_package(rocrand REQUIRED)
+    target_link_libraries(HYPRE INTERFACE rocsparse rocrand)
+  endif()
+  if (MFEM_USE_SINGLE)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_ENABLE_SINGLE:BOOL=ON)
+  endif()
+  # define external project and create future include directory so it is present
+  # to pass CMake checks at end of MFEM configuration step
+  message(STATUS "Will fetch HYPRE ${HYPRE_FETCH_VERSION} to be built with ${CMAKE_OPTIONS}")
+  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/hypre)
+  include(ExternalProject)
+  ExternalProject_Add(hypre
+    GIT_REPOSITORY https://github.com/hypre-space/hypre.git
+    GIT_TAG v${HYPRE_FETCH_VERSION}
+    GIT_SHALLOW TRUE
+    UPDATE_DISCONNECTED TRUE
+    SOURCE_SUBDIR src
+    PREFIX ${PREFIX}
+    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${PREFIX} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${CMAKE_OPTIONS})
+  file(MAKE_DIRECTORY ${PREFIX}/include)
+  # set imported library target properties
+  add_dependencies(HYPRE hypre)
+  set_target_properties(HYPRE PROPERTIES
+    IMPORTED_LOCATION ${PREFIX}/lib/libHYPRE.a
+    INTERFACE_INCLUDE_DIRECTORIES ${PREFIX}/include)
+  # convert HYPRE version to integer
+  string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION ${HYPRE_FETCH_VERSION})
+  list(GET HYPRE_SPLIT_VERSION 0 HYPRE_MAJOR_VERSION)
+  list(GET HYPRE_SPLIT_VERSION 1 HYPRE_MINOR_VERSION)
+  list(GET HYPRE_SPLIT_VERSION 2 HYPRE_PATCH_VERSION)
+  math(EXPR HYPRE_VERSION "10000*${HYPRE_MAJOR_VERSION} + 100*${HYPRE_MINOR_VERSION} + ${HYPRE_PATCH_VERSION}")
+  # set cache variables that would otherwise be set after mfem_find_package call
+  set(HYPRE_VERSION ${HYPRE_VERSION} CACHE STRING "HYPRE version." FORCE)
+  return()
+endif()
+
 include(MfemCmakeUtilities)
 mfem_find_package(HYPRE HYPRE HYPRE_DIR "include" "HYPRE.h" "lib" "HYPRE"
  "Paths to headers required by HYPRE." "Libraries required by HYPRE."
@@ -9,10 +9,38 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

-# Defines the following variables:
+# Defines the following variables if fetching of TPLs is disabled (default):
 #   - METIS_FOUND
 #   - METIS_LIBRARIES
 #   - METIS_INCLUDE_DIRS
+#   - METIS_VERSION_5
+# otherwise, the following are defined:
+#   - METIS (imported library target)
+#   - METIS_VERSION_5 (cache variable)
+
+if (METIS_FETCH OR FETCH_TPLS)
+  set(METIS_FETCH_VERSION 4.0.3)
+  add_library(METIS STATIC IMPORTED)
+  # define external project
+  message(STATUS "Will fetch METIS ${METIS_FETCH_VERSION} to be built with default options")
+  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/metis)
+  include(ExternalProject)
+  ExternalProject_Add(metis
+    GIT_REPOSITORY https://github.com/mfem/tpls
+    GIT_TAG b60352fbe9675d374b00828055e55be4584c7995 # tag from 1/16/25
+    GIT_SHALLOW TRUE
+    UPDATE_DISCONNECTED TRUE
+    PREFIX ${PREFIX}
+    CONFIGURE_COMMAND tar -xzf ../metis/metis-${METIS_FETCH_VERSION}-mac.tgz --strip=1
+    INSTALL_COMMAND mkdir -p ${PREFIX}/lib && cp libmetis.a ${PREFIX}/lib/)
+  # set imported library target properties
+  add_dependencies(METIS metis)
+  set_target_properties(METIS PROPERTIES
+    IMPORTED_LOCATION ${PREFIX}/lib/libmetis.a)
+  # set cache variables that would otherwise be set after mfem_find_package call
+  set(METIS_VERSION_5 FALSE CACHE BOOL "Is METIS version 5?")
+  return()
+endif()

 include(MfemCmakeUtilities)
 mfem_find_package(METIS METIS METIS_DIR "include;Lib" "metis.h"
@@ -89,6 +89,12 @@ option(MFEM_ENABLE_EXAMPLES "Build all of the examples" OFF)
 option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
 option(MFEM_ENABLE_BENCHMARKS "Build all of the benchmarks" OFF)

+# Allow a user to specify fetching of certain third-party libraries instead of
+# searching for existing installations.
+option(FETCH_TPLS "Enable fetching of all supported third-party libraries" OFF)
+option(HYPRE_FETCH "Enable fetching of hypre" OFF)
+option(METIS_FETCH "Enable fetching of METIS" OFF)
+
 # Setting CXX/MPICXX on the command line or in user.cmake will overwrite the
 # autodetected C++ compiler.
 # set(CXX g++)
@@ -82,6 +82,8 @@ set(SRCS
  fe/fe_ser.cpp
  fe_coll.cpp
  fespace.cpp
+  derefmat_op.cpp
+  pderefmat_op.cpp
  geom.cpp
  gridfunc.cpp
  hybridization.cpp
@@ -163,6 +165,7 @@ set(SRCS
  hyperbolic.cpp
  integrator.cpp
  bounds.cpp
+  blockform.cpp
  )

 set(HDRS
@@ -275,6 +278,7 @@ set(HDRS
  hyperbolic.hpp
  integrator.hpp
  bounds.hpp
+  blockform.hpp
  )

 if (MFEM_USE_SIDRE)
@@ -0,0 +1,353 @@
+// 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.
+
+// Implementation of class BlockForm
+
+#include "../config/config.hpp"
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+void BlockForm::BuildProlongation()
+{
+   P = new BlockMatrix(dof_offsets, tdof_offsets);
+   R = new BlockMatrix(tdof_offsets, dof_offsets);
+   P->owns_blocks = 0;
+   R->owns_blocks = 0;
+
+   for (int i = 0; i<nblocks; i++)
+   {
+      const SparseMatrix * P_ = fes[i]->GetConformingProlongation();
+      P->SetBlock(i,i,const_cast<SparseMatrix*>(P_));
+      const SparseMatrix * R_ = fes[i]->GetRestrictionMatrix();
+      R->SetBlock(i,i,const_cast<SparseMatrix*>(R_));
+   }
+}
+
+void BlockForm::Finalize(int skip_zeros)
+{
+   if (mat) { mat->Finalize(skip_zeros); }
+   if (mat_e) { mat_e->Finalize(skip_zeros); }
+}
+
+void BlockForm::ConformingAssemble()
+{
+   Finalize(0);
+   if (!P) { BuildProlongation(); }
+
+   BlockMatrix * Pt = Transpose(*P);
+   BlockMatrix * PtA = mfem::Mult(*Pt, *mat);
+   // mat->owns_blocks = 0;
+   for (int i = 0; i<nblocks; i++)
+   {
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (mat->IsZeroBlock(i,j)) { continue; }
+         if (Pt->IsZeroBlock(i,i))
+         {
+            PtA->SetBlock(i,j,&mat->GetBlock(i,j));
+         }
+      }
+   }
+   delete mat;
+   if (mat_e)
+   {
+      BlockMatrix *PtAe = mfem::Mult(*Pt, *mat_e);
+      mat_e->owns_blocks = 0;
+      for (int i = 0; i<nblocks; i++)
+      {
+         for (int j = 0; j<nblocks; j++)
+         {
+            if (mat_e->IsZeroBlock(i,j)) { continue; }
+            SparseMatrix * tmp = &mat_e->GetBlock(i,j);
+            if (Pt->IsZeroBlock(i,i))
+            {
+               PtAe->SetBlock(i,j,tmp);
+            }
+            else
+            {
+               delete tmp;
+            }
+         }
+      }
+      delete mat_e;
+      mat_e = PtAe;
+      mat_e->owns_blocks = 1;
+   }
+   delete Pt;
+
+   mat = mfem::Mult(*PtA, *P);
+   PtA->owns_blocks = 0;
+   for (int i = 0; i<nblocks; i++)
+   {
+      for (int j = 0; j<nblocks; j++)
+      {
+         if (PtA->IsZeroBlock(j,i)) { continue; }
+         SparseMatrix * tmp = &PtA->GetBlock(j,i);
+         if (P->IsZeroBlock(i,i))
+         {
+            mat->SetBlock(j,i,tmp);
+         }
+         else
+         {
+            delete tmp;
+         }
+      }
+   }
+   delete PtA;
+
+   if (mat_e)
+   {
+      BlockMatrix *PtAeP = mfem::Mult(*mat_e, *P);
+      mat_e->owns_blocks = 0;
+      for (int i = 0; i<nblocks; i++)
+      {
+         for (int j = 0; j<nblocks; j++)
+         {
+            if (mat_e->IsZeroBlock(j,i)) { continue; }
+            SparseMatrix * tmp = &mat_e->GetBlock(j,i);
+            if (P->IsZeroBlock(i,i))
+            {
+               PtAeP->SetBlock(j,i,tmp);
+            }
+            else
+            {
+               delete tmp;
+            }
+         }
+      }
+
+      delete mat_e;
+      mat_e = PtAeP;
+   }
+   height = mat->Height();
+   width = mat->Width();
+}
+
+
+BlockForm::BlockForm(const Array<FiniteElementSpace*> fes_ ): fes(
+      fes_)
+{
+   nblocks = fes.Size();
+   bforms.SetSize(nblocks,nblocks);
+   mforms.SetSize(nblocks,nblocks);
+   dof_offsets.Append(0);
+   tdof_offsets.Append(0);
+   for (int i = 0; i<nblocks; i++)
+   {
+      dof_offsets.Append(fes[i]->GetVSize());
+      tdof_offsets.Append(fes[i]->GetTrueVSize());
+      for (int j = 0; j<nblocks; j++)
+      {
+         bforms(i,j) = nullptr;
+         mforms(i,j) = nullptr;
+      }
+   }
+   dof_offsets.PartialSum();
+   tdof_offsets.PartialSum();
+   diag_policy = mfem::Operator::DIAG_ONE;
+
+}
+
+void BlockForm::SetBlock(BilinearForm * bform, int row_idx, int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   bforms(row_idx,col_idx) = bform;
+}
+void BlockForm::SetBlock(MixedBilinearForm * mform, int row_idx,
+                         int col_idx)
+{
+   MFEM_VERIFY((row_idx >=0 && row_idx < nblocks), "row index out of bounds");
+   MFEM_VERIFY((col_idx >=0 && col_idx < nblocks), "col index out of bounds");
+   MFEM_VERIFY(!mforms(row_idx,col_idx), "Entry has already been set");
+   MFEM_VERIFY(!bforms(row_idx,col_idx), "Entry has already been set");
+   mforms(row_idx,col_idx) = mform;
+}
+
+/// Assemble the local matrix
+void BlockForm::Assemble(int skip_zeros)
+{
+   mat = new BlockMatrix(dof_offsets);
+   for (int i = 0; i<nblocks; i++)
+   {
+      int h = dof_offsets[i+1]-dof_offsets[i];
+      for (int j = 0; j<nblocks; j++)
+      {
+         int w = dof_offsets[j+1]-dof_offsets[j];
+         if (bforms(i,j))
+         {
+            bforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = bforms(i,j)->Height(), "inconsistent height of bilinear form");
+            MFEM_VERIFY(w = bforms(i,j)->Width(), "inconsistent width of bilinear form");
+            mat->SetBlock(i,j,&bforms(i,j)->SpMat());
+         }
+         else if (mforms(i,j))
+         {
+            mforms(i,j)->Assemble(skip_zeros);
+            MFEM_VERIFY(h = mforms(i,j)->Height(),
+                        "inconsistent height of MixedBilinear form");
+            MFEM_VERIFY(w = mforms(i,j)->Width(),
+                        "inconsistent width of Mixedbilinear form");
+            mat->SetBlock(i,j,&mforms(i,j)->SpMat());
+         }
+         else
+         {
+            mat->SetBlock(i,j,nullptr);
+         }
+      }
+   }
+}
+
+void BlockForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
+                                 Vector & b,
+                                 OperatorHandle &A, Vector &X,
+                                 Vector &B, int copy_interior)
+{
+   FormSystemMatrix(ess_tdof_list, A);
+
+   if (!P)
+   {
+      EliminateVDofsInRHS(ess_tdof_list, x, b);
+      X.MakeRef(x, 0, x.Size());
+      B.MakeRef(b, 0, b.Size());
+      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   }
+   else // non conforming space
+   {
+      B.SetSize(P->Width());
+
+      P->MultTranspose(b, B);
+      real_t *data = b.GetData();
+      Vector tmp;
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (P->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            B.SetVector(tmp,offset);
+         }
+      }
+
+      X.SetSize(R->Height());
+
+      R->Mult(x, X);
+      data = x.GetData();
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (R->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            X.SetVector(tmp,offset);
+         }
+      }
+
+      EliminateVDofsInRHS(ess_tdof_list, X, B);
+      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   }
+}
+
+void BlockForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
+                                 OperatorHandle &A)
+{
+   if (!mat_e)
+   {
+      bool conforming = true;
+      for (int i = 0; i<nblocks; i++)
+      {
+         const SparseMatrix *P_ = fes[i]->GetConformingProlongation();
+         if (P_)
+         {
+            conforming = false;
+            break;
+         }
+      }
+      if (!conforming) { ConformingAssemble(); }
+      const int remove_zeros = 0;
+      EliminateVDofs(ess_tdof_list, diag_policy);
+      Finalize(remove_zeros);
+   }
+   A.Reset(mat, false);
+}
+
+void BlockForm::RecoverFEMSolution(const Vector &X, Vector &x)
+{
+   if (!P)
+   {
+      x.SyncMemory(X);
+   }
+   else
+   {
+      x.SetSize(P->Height());
+      P->Mult(X, x);
+      real_t *data = X.GetData();
+      Vector tmp;
+      for (int i = 0; i<nblocks; i++)
+      {
+         if (P->IsZeroBlock(i,i))
+         {
+            int offset = tdof_offsets[i];
+            tmp.SetDataAndSize(&data[offset],tdof_offsets[i+1]-tdof_offsets[i]);
+            x.SetVector(tmp,offset);
+         }
+      }
+   }
+}
+
+void BlockForm::EliminateVDofs(const Array<int> &vdofs,
+                               Operator::DiagonalPolicy dpolicy)
+{
+   if (mat_e == NULL)
+   {
+      Array<int> offsets;
+
+      offsets.MakeRef( (P) ? tdof_offsets : dof_offsets);
+
+      mat_e = new BlockMatrix(offsets);
+      mat_e->owns_blocks = 1;
+      for (int i = 0; i<mat_e->NumRowBlocks(); i++)
+      {
+         int h = offsets[i+1] - offsets[i];
+         for (int j = 0; j<mat_e->NumColBlocks(); j++)
+         {
+            int w = offsets[j+1] - offsets[j];
+            mat_e->SetBlock(i,j,new SparseMatrix(h, w));
+         }
+      }
+   }
+
+   mat->EliminateRowCols(vdofs,mat_e,diag_policy);
+}
+
+void BlockForm::EliminateVDofsInRHS(
+   const Array<int> &vdofs, const Vector &x, Vector &b)
+{
+   mat_e->AddMult(x,b,-1.);
+   mat->PartMult(vdofs,x,b);
+}
+
+BlockForm::~BlockForm()
+{
+   delete mat_e;
+   mat_e = nullptr;
+   delete mat;
+   mat = nullptr;
+   delete P;
+   delete R;
+}
+
+};
@@ -0,0 +1,88 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_BLOCKFORM
+#define MFEM_BLOCKFORM
+
+#include "../config/config.hpp"
+
+#include "fespace.hpp"
+#include "gridfunc.hpp"
+#include "bilinearform.hpp"
+
+namespace mfem
+{
+// square block forms
+class BlockForm
+{
+
+private:
+   int height, width;
+   int nblocks;
+   Array<int> dof_offsets;
+   Array<int> tdof_offsets;
+   // BilinearForms
+   Array2D<BilinearForm * > bforms;
+   Array2D<MixedBilinearForm * > mforms;
+   Array<FiniteElementSpace *> fes;
+
+   // Block Prolongation
+   BlockMatrix * P = nullptr;
+   // Block Restriction
+   BlockMatrix * R = nullptr;
+
+   BlockMatrix * mat = nullptr;
+   BlockMatrix * mat_e = nullptr;
+
+   void BuildProlongation();
+   void ConformingAssemble();
+
+   mfem::Operator::DiagonalPolicy diag_policy;
+
+public:
+
+   BlockForm(const Array<FiniteElementSpace*> pfes_ );
+
+   void SetBlock(BilinearForm * bform, int row_idx, int col_idx);
+   void SetBlock(MixedBilinearForm * mform, int row_idx, int col_idx);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector & b,
+                         OperatorHandle &A, Vector &X,
+                         Vector &B, int copy_interior = 0);
+
+   void FormSystemMatrix(const Array<int> &ess_tdof_list,
+                         OperatorHandle &A);
+
+   void RecoverFEMSolution(const Vector &X, Vector &x);
+
+   ///  Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   void EliminateVDofs(const Array<int> &vdofs,
+                       Operator::DiagonalPolicy dpolicy = Operator::DIAG_ONE);
+
+   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
+
+   void SetDiagonalPolicy(Operator::DiagonalPolicy policy)
+   {
+      diag_policy = policy;
+   }
+
+   /// Destroys bilinear form.
+   ~BlockForm();
+};
+
+}
+
+#endif
@@ -912,7 +912,7 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,

   if (vdim == 1) // scalar case
   {
-      n_field["values"].set_external(gf->GetData(),
+      n_field["values"].set_external(const_cast<real_t *>(gf->HostRead()),
                                     ndofs);
   }
   else // vector case
@@ -925,18 +925,18 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,
      int vdim_stride  = (ordering == Ordering::byNODES ? ndofs : 1);

      index_t offset = 0;
-      index_t stride = sizeof(double) * entry_stride;
+      index_t stride = sizeof(real_t) * entry_stride;

      for (int d = 0;  d < vdim; d++)
      {
         std::ostringstream oss;
         oss << "v" << d;
         std::string comp_name = oss.str();
-         n_field["values"][comp_name].set_external(gf->GetData(),
+         n_field["values"][comp_name].set_external(const_cast<real_t *>(gf->HostRead()),
                                                   ndofs,
                                                   offset,
                                                   stride);
-         offset +=  sizeof(double) * vdim_stride;
+         offset +=  sizeof(real_t) * vdim_stride;
      }
   }

@@ -0,0 +1,266 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "derefmat_op.hpp"
+#include "fes_kernels.hpp"
+
+/// \cond DO_NOT_DOCUMENT
+namespace mfem
+{
+
+namespace internal
+{
+template <Ordering::Type Order, bool Atomic>
+static void DerefMultKernelImpl(const DerefineMatrixOp &op, const Vector &x,
+                                Vector &y)
+{
+   DerefineMatrixOpMultFunctor<Order, Atomic> func;
+   func.xptr = x.Read();
+   y.UseDevice();
+   y = 0.;
+   func.yptr = y.ReadWrite();
+   func.bsptr = op.block_storage.Read();
+   func.boptr = op.block_offsets.Read();
+   func.brptr = op.block_row_idcs_offsets.Read();
+   func.bcptr = op.block_col_idcs_offsets.Read();
+   func.rptr = op.row_idcs.Read();
+   func.cptr = op.col_idcs.Read();
+   func.vdims = op.fespace->GetVDim();
+   func.nblocks = op.block_offsets.Size();
+   func.width = op.Width() / func.vdims;
+   func.height = op.Height() / func.vdims;
+   func.Run(op.max_rows);
+}
+
+} // namespace internal
+
+DerefineMatrixOp::DerefineMatrixOp(FiniteElementSpace &fespace_, int old_ndofs,
+                                   const Table *old_elem_dof,
+                                   const Table *old_elem_fos)
+   : Operator(fespace_.GetVSize(), old_ndofs * fespace_.GetVDim()),
+     fespace(&fespace_)
+{
+   static Kernels kernels;
+   constexpr int max_team_size = 256;
+   /// TODO: Implement DofTransformation support
+
+   MFEM_VERIFY(fespace->Nonconforming(),
+               "Not implemented for conforming meshes.");
+   MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
+   MFEM_VERIFY(fespace->GetNDofs() <= old_ndofs,
+               "Previous space is not finer.");
+
+   const CoarseFineTransformations &dtrans =
+      fespace->GetMesh()->ncmesh->GetDerefinementTransforms();
+
+   MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
+
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   DenseMatrix localRVO; // for variable-order only
+
+   DenseTensor localR[Geometry::NumGeom];
+   int total_rows = 0;
+   int total_cols = 0;
+   block_offsets.SetSize(dtrans.embeddings.Size());
+   block_offsets.HostWrite();
+   if (fespace->IsVariableOrder())
+   {
+      // TODO: any potential for some compression here?
+      // determine storage size and offsets
+      block_offsets[0] = 0;
+      int total_size = 0;
+      for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+      {
+         const Embedding &emb = dtrans.embeddings[k];
+         const FiniteElement *fe = fespace->GetFE(emb.parent);
+         const int ldof = fe->GetDof();
+         if (k + 1 < dtrans.embeddings.Size())
+         {
+            block_offsets[k + 1] = block_offsets[k] + ldof * ldof;
+         }
+         total_rows += ldof;
+         total_cols += ldof;
+         total_size += ldof * ldof;
+      }
+      block_storage.SetSize(total_size);
+   }
+   else
+   {
+      // compression scheme:
+      // block_offsets is the start of each block, potentially repeated
+      // only need to store localR for used shapes
+      Mesh::GeometryList elem_geoms(*fespace->GetMesh());
+
+      int geom_offsets[Geometry::NumGeom];
+      {
+         int size = 0;
+         for (int i = 0; i < elem_geoms.Size(); ++i)
+         {
+            fespace->GetLocalDerefinementMatrices(elem_geoms[i],
+                                                  localR[elem_geoms[i]]);
+            geom_offsets[elem_geoms[i]] = size;
+            size += localR[elem_geoms[i]].TotalSize();
+         }
+         block_storage.SetSize(size);
+         // copy blocks into block_storage
+         auto bs_ptr = block_storage.HostWrite();
+         for (int i = 0; i < elem_geoms.Size(); ++i)
+         {
+            std::copy(localR[elem_geoms[i]].Data(),
+                      localR[elem_geoms[i]].Data()
+                      + localR[elem_geoms[i]].TotalSize(),
+                      bs_ptr);
+            bs_ptr += localR[elem_geoms[i]].TotalSize();
+         }
+      }
+      for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+      {
+         const Embedding &emb = dtrans.embeddings[k];
+         Geometry::Type geom =
+            fespace->GetMesh()->GetElementBaseGeometry(emb.parent);
+
+         auto size = localR[geom].SizeI() * localR[geom].SizeJ();
+         total_rows += localR[geom].SizeI();
+         total_cols += localR[geom].SizeJ();
+         // set block offsets and sizes
+         block_offsets[k] = geom_offsets[geom] + size * emb.matrix;
+      }
+   }
+   row_idcs.SetSize(total_rows);
+   row_idcs.HostWrite();
+   col_idcs.SetSize(total_cols);
+   col_idcs.HostWrite();
+   block_row_idcs_offsets.SetSize(dtrans.embeddings.Size() + 1);
+   block_row_idcs_offsets.HostWrite();
+   block_col_idcs_offsets.SetSize(dtrans.embeddings.Size() + 1);
+   block_col_idcs_offsets.HostWrite();
+   block_row_idcs_offsets[0] = 0;
+   block_col_idcs_offsets[0] = 0;
+
+   // compute index information
+   Array<int> dofs, old_dofs;
+   max_rows = 1;
+
+   {
+      Array<int> mark(fespace->GetNDofs());
+      mark = 0;
+      auto bs_ptr = block_storage.HostWrite();
+      int ridx = 0;
+      int cidx = 0;
+      int num_marked = 0;
+      for (int k = 0; k < dtrans.embeddings.Size(); k++)
+      {
+         const Embedding &emb = dtrans.embeddings[k];
+         Geometry::Type geom =
+            fespace->GetMesh()->GetElementBaseGeometry(emb.parent);
+
+         if (fespace->IsVariableOrder())
+         {
+            const FiniteElement *fe = fespace->GetFE(emb.parent);
+            const DenseTensor &pmats = dtrans.point_matrices[geom];
+            const int ldof = fe->GetDof();
+
+            IsoparametricTransformation isotr;
+            isotr.SetIdentityTransformation(geom);
+
+            localRVO.SetSize(ldof, ldof);
+            isotr.SetPointMat(pmats(emb.matrix));
+            // Local restriction is size ldofxldof assuming that the parent
+            // and child are of same polynomial order.
+            fe->GetLocalRestriction(isotr, localRVO);
+            // copy block
+            auto size = localRVO.Height() * localRVO.Width();
+            std::copy(localRVO.Data(), localRVO.Data() + size, bs_ptr);
+            bs_ptr += size;
+         }
+         DenseMatrix &lR =
+            fespace->IsVariableOrder() ? localRVO : localR[geom](emb.matrix);
+         block_row_idcs_offsets[k + 1] =
+            block_row_idcs_offsets[k] + lR.Height();
+         block_col_idcs_offsets[k + 1] = block_col_idcs_offsets[k] + lR.Width();
+         max_rows = std::max(lR.Height(), max_rows);
+         // index information
+         fespace->elem_dof->GetRow(emb.parent, dofs);
+         old_elem_dof->GetRow(k, old_dofs);
+         MFEM_VERIFY(old_dofs.Size() == dofs.Size(),
+                     "Parent and child must have same #dofs.");
+         for (int i = 0; i < lR.Height(); ++i, ++ridx)
+         {
+            if (!std::isfinite(lR(i, 0)))
+            {
+               row_idcs[ridx] = INT_MAX;
+               continue;
+            }
+            int r = dofs[i];
+            int m = (r >= 0) ? r : (-1 - r);
+            if (is_dg || !mark[m])
+            {
+               row_idcs[ridx] = r;
+               mark[m] = 1;
+               ++num_marked;
+            }
+            else
+            {
+               row_idcs[ridx] = INT_MAX;
+            }
+         }
+         for (int i = 0; i < lR.Width(); ++i, ++cidx)
+         {
+            col_idcs[cidx] = old_dofs[i];
+         }
+      }
+      if (!is_dg && !fespace->IsVariableOrder())
+      {
+         MFEM_VERIFY(num_marked * fespace->GetVDim() == Height(),
+                     "internal error: not all rows were set.");
+      }
+   }
+   // if not using GPU, set max_rows/max_cols to zero
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      max_rows = std::min(max_rows, max_team_size);
+   }
+   else
+   {
+      max_rows = 1;
+   }
+}
+
+void DerefineMatrixOp::Mult(const Vector &x, Vector &y) const
+{
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   // DG needs atomic summation
+   MultKernel::Run(fespace->GetOrdering(), is_dg, *this, x, y);
+}
+
+DerefineMatrixOp::Kernels::Kernels()
+{
+   MultKernel::Specialization<Ordering::byNODES, false>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, false>::Add();
+   MultKernel::Specialization<Ordering::byNODES, true>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, true>::Add();
+}
+
+template <Ordering::Type Order, bool Atomic>
+DerefineMatrixOp::MultKernelType DerefineMatrixOp::MultKernel::Kernel()
+{
+   return internal::DerefMultKernelImpl<Order, Atomic>;
+}
+
+DerefineMatrixOp::MultKernelType
+DerefineMatrixOp::MultKernel::Fallback(Ordering::Type, bool)
+{
+   MFEM_ABORT("invalid MultKernel parameters");
+}
+} // namespace mfem
+/// \endcond DO_NOT_DOCUMENT
@@ -0,0 +1,65 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_DEREFMAT_OP
+#define MFEM_DEREFMAT_OP
+
+#include "fespace.hpp"
+
+#include "kernel_dispatch.hpp"
+
+namespace mfem
+{
+
+/// \cond DO_NOT_DOCUMENT
+
+struct DerefineMatrixOp : public Operator
+{
+   FiniteElementSpace *fespace;
+   /// offsets into block_storage
+   Array<int> block_offsets;
+   /// offsets into row_idcs
+   Array<int> block_row_idcs_offsets;
+   /// offsets into col_idcs
+   Array<int> block_col_idcs_offsets;
+   /// mapping for row dofs, INT_MAX indicates the block row should be ignored.
+   /// negative means the row data should be negated.
+   Array<int> row_idcs;
+   /// mapping for col dofs, negative means the col data should be negated.
+   Array<int> col_idcs;
+   /// dense block matrices which can be reused to construct the full matrix
+   /// operation. These are stored contiguously and blocks have no restrictions
+   /// on shape (can be rectangle and differ from block to block).
+   Vector block_storage;
+   /// maximum height of any block in block_storage for GPU
+   /// parallelization, or 1 for CPU runs.
+   int max_rows;
+
+   using MultKernelType = void (*)(const DerefineMatrixOp &, const Vector &,
+                                   Vector &);
+   /// template args: ordering, atomic
+   MFEM_REGISTER_KERNELS(MultKernel, MultKernelType, (Ordering::Type, bool));
+
+   struct Kernels
+   {
+      Kernels();
+   };
+
+   void Mult(const Vector &x, Vector &y) const;
+
+   DerefineMatrixOp(FiniteElementSpace &fespace_, int old_ndofs,
+                    const Table *old_elem_dof, const Table *old_elem_fos);
+};
+
+/// \endcond DO_NOT_DOCUMENT
+
+} // namespace mfem
+#endif
@@ -50,8 +50,8 @@
 #include "dgmassinv.hpp"
 #include "hyperbolic.hpp"
 #include "bounds.hpp"
-
 #include "dfem/doperator.hpp"
+#include "blockform.hpp"

 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -0,0 +1,249 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_FES_KERNELS_HPP
+#define MFEM_FES_KERNELS_HPP
+
+#include "../general/forall.hpp"
+
+#include <climits>
+
+namespace mfem
+{
+/// \cond DO_NOT_DOCUMENT
+namespace internal
+{
+
+///
+/// Implements matrix-vector multiply $y = A x$ for a sparse matrix composed of
+/// a sum of smaller dense blocks. There is additional permutation/sign
+/// information associated with each block. The base class only implements
+/// helper routines such as computing block widths, index into x, index into y,
+/// and column in A given sub-block information.
+/// @sa DerefineMatrixOpMultFunctor
+///
+/// @tparam Order vdim ordering for x and y. Note that for Diag = false this is
+/// ignored for x as x has a special interleaved order.
+/// @tparam Base used for the curious recurring template pattern (CRTP) so the
+/// base class can access child class fields without virtual functions
+/// @tparam Diag true if this corresponds to the diagonal block (coarse element
+/// and fine element are on our rank), false otherwise (coarse element is on our
+/// rank, fine element is on a different rank).
+///
+template <Ordering::Type Order, class Base, bool Diag = true>
+struct DerefineMatrixOpFunctorBase;
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byNODES, Base, true>
+{
+   /// block column indices offsets
+   const int *bcptr;
+   /// column indices
+   const int *cptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const
+   {
+      return bcptr[k + 1] - bcptr[k];
+   }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = cptr[bcptr[k] + j];
+      if (col < 0)
+      {
+         col = -1 - col;
+         sign = -sign;
+      }
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int) const
+   {
+      return col + vdim * static_cast<const Base *>(this)->width;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return row + vdim * static_cast<const Base *>(this)->height;
+   }
+};
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byVDIM, Base, true>
+{
+   /// block column indices offsets
+   const int *bcptr;
+   /// column indices
+   const int *cptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const
+   {
+      return bcptr[k + 1] - bcptr[k];
+   }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = cptr[bcptr[k] + j];
+      if (col < 0)
+      {
+         col = -1 - col;
+         sign = -sign;
+      }
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int) const
+   {
+      return vdim + col * static_cast<const Base *>(this)->vdims;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return vdim + row * static_cast<const Base *>(this)->vdims;
+   }
+};
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byNODES, Base, false>
+{
+   /// receive segment offsets
+   const int *segptr;
+   /// receive segment index
+   const int *rsptr;
+   /// off-diagonal block column offsets
+   const int *coptr;
+   /// off-diagonal block widths
+   const int *bwptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const { return bwptr[k]; }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = coptr[k] + j;
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int k) const
+   {
+      int tmp = rsptr[k];
+      int segwidth = segptr[tmp + 1] - segptr[tmp];
+      return segptr[tmp] * static_cast<const Base *>(this)->vdims + col +
+             vdim * segwidth;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return row + vdim * static_cast<const Base *>(this)->height;
+   }
+};
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byVDIM, Base, false>
+{
+   /// receive segment offsets
+   const int *segptr;
+   /// receive segment index
+   const int *rsptr;
+   /// off-diagonal block column offsets
+   const int *coptr;
+   /// off-diagonal block widths
+   const int *bwptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const { return bwptr[k]; }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = coptr[k] + j;
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int k) const
+   {
+      int tmp = rsptr[k];
+      int segwidth = segptr[tmp + 1] - segptr[tmp];
+      return segptr[tmp] * static_cast<const Base *>(this)->vdims + col +
+             vdim * segwidth;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return vdim + row * static_cast<const Base *>(this)->vdims;
+   }
+};
+
+/// internally used to implement the derefinement operator Mult diagonal
+/// block
+template <Ordering::Type Order, bool Atomic, bool Diag = true>
+struct DerefineMatrixOpMultFunctor
+   : public DerefineMatrixOpFunctorBase<
+     Order, DerefineMatrixOpMultFunctor<Order, Atomic, Diag>, Diag>
+{
+   const real_t *xptr;
+   real_t *yptr;
+   /// block storage
+   const real_t *bsptr;
+   /// block offsets
+   const int *boptr;
+   /// block row index offsets
+   const int *brptr;
+   /// row indices
+   const int *rptr;
+
+   // number of blocks
+   int nblocks;
+   // number of components
+   int vdims;
+   /// overall operator height (for vdim = 1)
+   int height;
+   /// overall operator width (for vdim = 1)
+   int width;
+   void MFEM_HOST_DEVICE operator()(int kidx) const
+   {
+      int k = kidx % nblocks;
+      int vdim = kidx / nblocks;
+
+      int block_height = brptr[k + 1] - brptr[k];
+      int block_width = this->BlockWidth(k);
+      MFEM_FOREACH_THREAD(i, x, block_height)
+      {
+         int row = rptr[brptr[k] + i];
+         int rsign = 1;
+         if (row < 0)
+         {
+            row = -1 - row;
+            rsign = -1;
+         }
+         if (row < INT_MAX)
+         {
+            // row not marked as unused
+            real_t sum = 0;
+            for (int j = 0; j < block_width; ++j)
+            {
+               int col, sign = rsign;
+               this->Col(j, k, col, sign);
+               sum += sign * bsptr[boptr[k] + i + j * block_height] *
+                      xptr[this->IndexX(col, vdim, k)];
+            }
+#if defined(__CUDA_ARCH__) or defined(__HIP_DEVICE_COMPILE__)
+            if (Atomic)
+            {
+               atomicAdd(yptr + this->IndexY(row, vdim), sum);
+            }
+            else
+#endif
+            {
+               yptr[this->IndexY(row, vdim)] += sum;
+            }
+         }
+      }
+   }
+
+   /// N is the max block row size (doesn't have to be a power of 2)
+   void Run(int N) const { forall_2D(nblocks * vdims, N, 1, *this); }
+};
+
+} // namespace internal
+/// \endcond DO_NOT_DOCUMENT
+} // namespace mfem
+
+#endif
@@ -17,6 +17,9 @@
 #include "fem.hpp"
 #include "ceed/interface/util.hpp"

+#include "derefmat_op.hpp"
+
+#include <algorithm>
 #include <cmath>
 #include <cstdarg>

@@ -24,9 +27,9 @@ using namespace std;

 namespace mfem
 {
-
-template <> void Ordering::
-DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
+template <>
+void Ordering::DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim,
+                                              Array<int> &dofs)
 {
   // static method
   int size = dofs.Size();
@@ -40,8 +43,9 @@ DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
   }
 }

-template <> void Ordering::
-DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
+template <>
+void Ordering::DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim,
+                                             Array<int> &dofs)
 {
   // static method
   int size = dofs.Size();
@@ -55,7 +59,6 @@ DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
   }
 }

-
 FiniteElementSpace::FiniteElementSpace()
   : mesh(NULL), fec(NULL), vdim(0), ordering(Ordering::byNODES),
     ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0),
@@ -4244,7 +4247,11 @@ void FiniteElementSpace::Update(bool want_transform)
         case Mesh::DEREFINE:
         {
            BuildConformingInterpolation();
+#if 0
            Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof, old_elem_fos));
+#else
+            Th.Reset(new DerefineMatrixOp(*this, old_ndofs, old_elem_dof, old_elem_fos));
+#endif
            if (IsVariableOrder())
            {
               if (cP && cR_hp)
@@ -113,7 +113,7 @@ class QuadratureSpace;
 class QuadratureInterpolator;
 class FaceQuadratureInterpolator;
 class PRefinementTransferOperator;
-
+struct DerefineMatrixOp;

 /** @brief Class FiniteElementSpace - responsible for providing FEM view of the
    mesh, mainly managing the set of degrees of freedom.
@@ -246,6 +246,7 @@ class FiniteElementSpace
   friend class PRefinementTransferOperator;
   friend void Mesh::Swap(Mesh &, bool);
   friend class LORBase;
+   friend struct DerefineMatrixOp;

 protected:
   /// The mesh that FE space lives on (not owned).
@@ -346,13 +346,13 @@ private:
      template<typename T>
      T operator() (const blitz::TinyVector<T,3>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
-         TmplPoly_1D::CalcBernstein(el_order, x[2], u3);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[2], u3.data());

         const Array<int>& dof_map=el->GetDofMap();

@@ -370,17 +370,17 @@ private:
      template<typename T>
      blitz::TinyVector<T,3> grad(const blitz::TinyVector<T,3>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
-         T d3[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);
+         std::vector<T> d3(el_order+1);

-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
-         TmplPoly_1D::CalcBernstein(el_order,x[2], u3, d3);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[2], u3.data(), d3.data());

         blitz::TinyVector<T,3> res(T(0.0),T(0.0),T(0.0));

@@ -415,11 +415,11 @@ private:
      template<typename T>
      T operator() (const blitz::TinyVector<T,2>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());

         const Array<int>& dof_map=el->GetDofMap();

@@ -437,14 +437,14 @@ private:
      template<typename T>
      blitz::TinyVector<T,2> grad(const blitz::TinyVector<T,2>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);

-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());

         blitz::TinyVector<T,2> res(T(0.0),T(0.0));

@@ -0,0 +1,591 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "pderefmat_op.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "fes_kernels.hpp"
+/// \cond DO_NOT_DOCUMENT
+namespace mfem
+{
+namespace internal
+{
+template <Ordering::Type Order, bool Atomic>
+static void ParDerefMultKernelImpl(const ParDerefineMatrixOp &op,
+                                   const Vector &x, Vector &y)
+{
+   // pack sends
+   if (op.xghost_send.Size())
+   {
+      auto src = x.Read();
+      auto idcs = op.send_permutations.Read();
+      auto dst = Device::GetGPUAwareMPI() ? op.xghost_send.Write()
+                 : op.xghost_send.HostWrite();
+      auto vdims = op.fespace->GetVDim();
+      auto sptr = op.send_segment_idcs.Read();
+      auto lptr = op.send_segments.Read();
+      auto old_ndofs = x.Size() / vdims;
+
+      forall(op.send_permutations.Size(), [=] MFEM_HOST_DEVICE(int i)
+      {
+         int seg = sptr[i];
+         int width = lptr[seg + 1] - lptr[seg];
+         auto tdst = dst + i + lptr[seg] * vdims;
+         int sign = 1;
+         int col = idcs[i];
+         if (col < 0)
+         {
+            sign = -1;
+            col = -1 - col;
+         }
+         for (int vdim = 0; vdim < vdims; ++vdim)
+         {
+            tdst[vdim * width] =
+               sign
+               * src[Order == Ordering::byNODES ? (col + vdim * old_ndofs)
+                           : (col * vdims + vdim)];
+         }
+      });
+      // TODO: is this needed so we can send the packed data correctly?
+      // unclear for GPU-aware MPI, definitely required otherwise
+      MFEM_DEVICE_SYNC;
+   }
+   // initialize off-diagonal receive and send
+   op.requests.clear();
+   if (op.xghost_recv.Size())
+   {
+      auto vdims = op.fespace->GetVDim();
+      auto rcv = Device::GetGPUAwareMPI() ? op.xghost_recv.Write()
+                 : op.xghost_recv.HostWrite();
+      for (int i = 0; i < op.recv_ranks.Size(); ++i)
+      {
+         op.requests.emplace_back();
+         MPI_Irecv(rcv + op.recv_segments[i] * vdims,
+                   (op.recv_segments[i + 1] - op.recv_segments[i]) * vdims,
+                   MPITypeMap<real_t>::mpi_type, op.recv_ranks[i],
+                   MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA,
+                   op.fespace->GetComm(), &op.requests.back());
+      }
+   }
+   if (op.xghost_send.Size())
+   {
+      auto vdims = op.fespace->GetVDim();
+      // only is a GPU mem ptr if GPU-aware MPI is enabled
+      auto dst = Device::GetGPUAwareMPI() ? op.xghost_send.Write()
+                 : op.xghost_send.HostWrite();
+      for (int i = 0; i < op.send_ranks.Size(); ++i)
+      {
+         op.requests.emplace_back();
+         MPI_Isend(dst + op.send_segments[i] * vdims,
+                   (op.send_segments[i + 1] - op.send_segments[i]) * vdims,
+                   MPITypeMap<real_t>::mpi_type, op.send_ranks[i],
+                   MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA,
+                   op.fespace->GetComm(), &op.requests.back());
+      }
+   }
+   {
+      // diagonal
+      DerefineMatrixOpMultFunctor<Order, Atomic, true> func;
+      func.xptr = x.Read();
+      y.UseDevice();
+      y = 0.;
+      func.yptr = y.ReadWrite();
+      func.bsptr = op.block_storage.Read();
+      func.boptr = op.block_offsets.Read();
+      func.brptr = op.block_row_idcs_offsets.Read();
+      func.bcptr = op.block_col_idcs_offsets.Read();
+      func.rptr = op.row_idcs.Read();
+      func.cptr = op.col_idcs.Read();
+      func.vdims = op.fespace->GetVDim();
+      func.nblocks = op.block_offsets.Size();
+      func.width = op.Width() / func.vdims;
+      func.height = op.Height() / func.vdims;
+      func.Run(op.max_rows);
+   }
+   // wait for comm to finish, if any
+   if (op.requests.size())
+   {
+      MPI_Waitall(op.requests.size(), op.requests.data(), MPI_STATUSES_IGNORE);
+      if (op.xghost_recv.Size())
+      {
+         // off-diagonal kernel
+         DerefineMatrixOpMultFunctor<Order, Atomic, false> func;
+         // directly read from host-pinned memory if not using GPU-aware MPI
+         func.xptr = Device::GetGPUAwareMPI() ? op.xghost_recv.Read()
+                     : op.xghost_recv.HostRead();
+         func.yptr = y.ReadWrite();
+         func.bsptr = op.block_storage.Read();
+         func.boptr = op.off_diag_block_offsets.Read();
+         func.brptr = op.block_off_diag_row_idcs_offsets.Read();
+         func.rsptr = op.recv_segment_idcs.Read();
+         func.segptr = op.recv_segments.Read();
+         func.coptr = op.block_off_diag_col_offsets.Read();
+         func.bwptr = op.block_off_diag_widths.Read();
+         func.rptr = op.row_off_diag_idcs.Read();
+         func.vdims = op.fespace->GetVDim();
+         func.nblocks = op.off_diag_block_offsets.Size();
+         func.width = op.xghost_recv.Size() / func.vdims;
+         func.height = op.Height() / func.vdims;
+         func.Run(op.max_rows);
+      }
+   }
+}
+} // namespace internal
+
+template <Ordering::Type Order, bool Atomic>
+ParDerefineMatrixOp::MultKernelType ParDerefineMatrixOp::MultKernel::Kernel()
+{
+   return internal::ParDerefMultKernelImpl<Order, Atomic>;
+}
+
+ParDerefineMatrixOp::MultKernelType
+ParDerefineMatrixOp::MultKernel::Fallback(Ordering::Type, bool)
+{
+   MFEM_ABORT("invalid MultKernel parameters");
+}
+
+ParDerefineMatrixOp::Kernels::Kernels()
+{
+   MultKernel::Specialization<Ordering::byNODES, false>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, false>::Add();
+   MultKernel::Specialization<Ordering::byNODES, true>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, true>::Add();
+}
+
+void ParDerefineMatrixOp::Mult(const Vector &x, Vector &y) const
+{
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   // DG needs atomic summation
+   MultKernel::Run(fespace->GetOrdering(), is_dg, *this, x, y);
+   // use this to prevent xghost* from being re-purposed for subsequent Mult
+   // calls
+   MFEM_DEVICE_SYNC;
+}
+
+ParDerefineMatrixOp::ParDerefineMatrixOp(ParFiniteElementSpace &fespace_,
+                                         int old_ndofs,
+                                         const Table *old_elem_dof,
+                                         const Table *old_elem_fos)
+   : Operator(fespace_.GetVSize(), old_ndofs * fespace_.GetVDim()),
+     fespace(&fespace_)
+{
+   static Kernels kernels;
+   constexpr int max_team_size = 256;
+
+   const int NRanks = fespace->GetNRanks();
+
+   const int nrk = HYPRE_AssumedPartitionCheck() ? 2 : NRanks;
+
+   MFEM_VERIFY(fespace->Nonconforming(),
+               "Not implemented for conforming meshes.");
+   MFEM_VERIFY(fespace->old_dof_offsets[nrk],
+               "Missing previous (finer) space.");
+
+   const int MyRank = fespace->GetMyRank();
+   ParNCMesh *old_pncmesh = fespace->GetParMesh()->pncmesh;
+   const CoarseFineTransformations &dtrans =
+      old_pncmesh->GetDerefinementTransforms();
+   const Array<int> &old_ranks = old_pncmesh->GetDerefineOldRanks();
+
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   DenseMatrix localRVO; // for variable-order only
+
+   DenseTensor localR[Geometry::NumGeom];
+   int diag_rows = 0;
+   int off_diag_rows = 0;
+   int diag_cols = 0;
+
+   auto get_ldofs = [&](int k) -> int
+   {
+      const Embedding &emb = dtrans.embeddings[k];
+      if (fespace->IsVariableOrder())
+      {
+         const FiniteElement *fe = fespace->GetFE(emb.parent);
+         return fe->GetDof();
+      }
+      else
+      {
+         Geometry::Type geom =
+         fespace->GetParMesh()->GetElementBaseGeometry(emb.parent);
+         return fespace->FEColl()->FiniteElementForGeometry(geom)->GetDof();
+      }
+   };
+   Array<int> dofs, old_dofs;
+   max_rows = 1;
+   // first pass:
+   // - determine memory block lengths
+   // - identify dofs in x we need to send/receive
+   // don't need to send the indices, fine rank will re-arrange and sign
+   // change x before transmitting the ghost data
+
+   // key: coarse rank to send to
+   // value: old dofs to send (with sign)
+   std::map<int, std::vector<int>> to_send;
+   // key: fine rank
+   // value: indices into dtrans.embeddings
+   std::map<int, std::vector<int>> od_ks;
+   // key: fine rank
+   // value: recv segment length
+   std::map<int, int> od_seg_lens;
+   int send_len = 0;
+   int recv_len = 0;
+   // size of block_storage, if fespace->IsVariableOrder()
+   // otherwise unused
+   int total_size = 0;
+   int num_diagonal_blocks = 0;
+   int num_offdiagonal_blocks = 0;
+   for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+   {
+      const Embedding &emb = dtrans.embeddings[k];
+      int fine_rank = old_ranks[k];
+      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
+                        : old_pncmesh->ElementRank(emb.parent);
+      if (coarse_rank != MyRank && fine_rank == MyRank)
+      {
+         // this rank needs to send data in x to course_rank
+         old_elem_dof->GetRow(k, old_dofs);
+         auto &tmp = to_send[coarse_rank];
+         send_len += old_dofs.Size();
+         for (int i = 0; i < old_dofs.Size(); ++i)
+         {
+            tmp.emplace_back(old_dofs[i]);
+         }
+      }
+      else if (coarse_rank == MyRank && fine_rank != MyRank)
+      {
+         // this rank needs to receive data in x from fine_rank
+         MFEM_ASSERT(emb.parent >= 0, "");
+         auto ldofs = get_ldofs(k);
+         off_diag_rows += ldofs;
+         recv_len += ldofs;
+         od_ks[fine_rank].emplace_back(k);
+         od_seg_lens[fine_rank] += ldofs;
+         ++num_offdiagonal_blocks;
+         if (fespace->IsVariableOrder())
+         {
+            total_size += ldofs * ldofs;
+         }
+      }
+      else if (coarse_rank == MyRank && fine_rank == MyRank)
+      {
+         MFEM_ASSERT(emb.parent >= 0, "");
+         // diagonal
+         ++num_diagonal_blocks;
+         auto ldofs = get_ldofs(k);
+         diag_rows += ldofs;
+         diag_cols += ldofs;
+         if (fespace->IsVariableOrder())
+         {
+            total_size += ldofs * ldofs;
+         }
+      }
+   }
+   send_segments.SetSize(to_send.size() + 1);
+   send_segments.HostWrite();
+   send_ranks.SetSize(to_send.size());
+   send_ranks.HostWrite();
+   {
+      int idx = 0;
+      send_segments[0] = 0;
+      for (auto &tmp : to_send)
+      {
+         send_ranks[idx] = tmp.first;
+         send_segments[idx + 1] = send_segments[idx] + tmp.second.size();
+         ++idx;
+      }
+   }
+   recv_segment_idcs.SetSize(off_diag_rows);
+   recv_segment_idcs.HostWrite();
+   recv_segments.SetSize(od_ks.size() + 1);
+   recv_segments.HostWrite();
+   recv_ranks.SetSize(od_ks.size());
+   recv_ranks.HostWrite();
+
+   // set sizes
+   row_idcs.SetSize(diag_rows);
+   row_idcs.HostWrite();
+   row_off_diag_idcs.SetSize(off_diag_rows);
+   row_off_diag_idcs.HostWrite();
+   col_idcs.SetSize(diag_cols);
+   col_idcs.HostWrite();
+   block_row_idcs_offsets.SetSize(num_diagonal_blocks + 1);
+   block_row_idcs_offsets.HostWrite();
+   block_col_idcs_offsets.SetSize(num_diagonal_blocks + 1);
+   block_col_idcs_offsets.HostWrite();
+   block_off_diag_row_idcs_offsets.SetSize(num_offdiagonal_blocks + 1);
+   block_off_diag_row_idcs_offsets.HostWrite();
+   block_off_diag_col_offsets.SetSize(num_offdiagonal_blocks);
+   block_off_diag_col_offsets.HostWrite();
+   block_off_diag_widths.SetSize(num_offdiagonal_blocks);
+   block_off_diag_widths.HostWrite();
+   pack_col_idcs.SetSize(send_len);
+   // memory manager doesn't appear to have a graceful fallback for
+   // HOST_PINNED if not built with CUDA or HIP
+#if defined(MFEM_USE_CUDA) or defined(MFEM_USE_HIP)
+   xghost_send.SetSize(send_len * fespace->GetVDim(),
+                       Device::GetGPUAwareMPI() ? MemoryType::DEFAULT
+                       : MemoryType::HOST_PINNED);
+   xghost_recv.SetSize(recv_len * fespace->GetVDim(),
+                       Device::GetGPUAwareMPI() ? MemoryType::DEFAULT
+                       : MemoryType::HOST_PINNED);
+#else
+   xghost_send.SetSize(send_len * fespace->GetVDim());
+   xghost_recv.SetSize(recv_len * fespace->GetVDim());
+#endif
+   send_permutations.SetSize(send_len);
+   send_segment_idcs.SetSize(send_len);
+   block_offsets.SetSize(num_diagonal_blocks);
+   block_offsets.HostWrite();
+   off_diag_block_offsets.SetSize(num_offdiagonal_blocks);
+   off_diag_block_offsets.HostWrite();
+   int geom_offsets[Geometry::NumGeom];
+   real_t *bs_ptr;
+
+   if (fespace->IsVariableOrder())
+   {
+      block_storage.SetSize(total_size);
+      bs_ptr = block_storage.HostWrite();
+      // compute block data later
+   }
+   else
+   {
+      // compression scheme:
+      // block_offsets is the start of each block, potentially repeated
+      // only need to store localR for used shapes
+      Mesh::GeometryList elem_geoms(*fespace->GetMesh());
+
+      int size = 0;
+      for (int i = 0; i < elem_geoms.Size(); ++i)
+      {
+         fespace->GetLocalDerefinementMatrices(elem_geoms[i],
+                                               localR[elem_geoms[i]]);
+         geom_offsets[elem_geoms[i]] = size;
+         size += localR[elem_geoms[i]].TotalSize();
+      }
+      block_storage.SetSize(size);
+      bs_ptr = block_storage.HostWrite();
+      // copy blocks into block_storage
+      for (int i = 0; i < elem_geoms.Size(); ++i)
+      {
+         std::copy(localR[elem_geoms[i]].Data(),
+                   localR[elem_geoms[i]].Data()
+                   + localR[elem_geoms[i]].TotalSize(),
+                   bs_ptr);
+         bs_ptr += localR[elem_geoms[i]].TotalSize();
+      }
+   }
+
+   // second pass:
+   // - initialize buffers
+
+   {
+      auto ptr = send_permutations.HostWrite();
+      auto ptr2 = send_segment_idcs.HostWrite();
+      int i = 0;
+      for (auto &v : to_send)
+      {
+         ptr = std::copy(v.second.begin(), v.second.end(), ptr);
+         for (size_t idx = 0; idx < v.second.size(); ++idx)
+         {
+            *ptr2 = i;
+            ++ptr2;
+         }
+         ++i;
+      }
+   }
+
+   block_row_idcs_offsets[0] = 0;
+   block_col_idcs_offsets[0] = 0;
+   block_off_diag_row_idcs_offsets[0] = 0;
+   Array<int> mark(fespace->GetNDofs());
+   mark = 0;
+   {
+      int idx = 0;
+      recv_segments[0] = 0;
+      for (auto &v : od_seg_lens)
+      {
+         recv_ranks[idx] = v.first;
+         recv_segments[idx + 1] = recv_segments[idx] + v.second;
+         ++idx;
+      }
+   }
+   // key: index into dtrans.embeddings
+   // value: off-diagonal block offset, od_ridx, seg id
+   std::unordered_map<int, std::array<int, 3>> ks_map;
+   {
+      int od_ridx = 0;
+      int seg_id = 0;
+      for (auto &v1 : od_ks)
+      {
+         for (auto k : v1.second)
+         {
+            auto &tmp = ks_map[k];
+            tmp[0] = ks_map.size() - 1;
+            tmp[1] = od_ridx;
+            tmp[2] = seg_id;
+            od_ridx += get_ldofs(k);
+         }
+         ++seg_id;
+      }
+   }
+   int diag_idx = 0;
+   int var_offset = 0;
+   int ridx = 0;
+   int cidx = 0;
+   // can't break this up into separate diagonals/off-diagonals loops because
+   // of mark
+   for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+   {
+      const Embedding &emb = dtrans.embeddings[k];
+      if (emb.parent < 0)
+      {
+         continue;
+      }
+      int fine_rank = old_ranks[k];
+      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
+                        : old_pncmesh->ElementRank(emb.parent);
+      if (coarse_rank == MyRank)
+      {
+         // either diagonal or off-diagonal
+         Geometry::Type geom =
+            fespace->GetMesh()->GetElementBaseGeometry(emb.parent);
+         if (fespace->IsVariableOrder())
+         {
+            const FiniteElement *fe = fespace->GetFE(emb.parent);
+            const DenseTensor &pmats = dtrans.point_matrices[geom];
+            const int ldof = fe->GetDof();
+
+            IsoparametricTransformation isotr;
+            isotr.SetIdentityTransformation(geom);
+
+            localRVO.SetSize(ldof, ldof);
+            isotr.SetPointMat(pmats(emb.matrix));
+            // Local restriction is size ldofxldof assuming that the parent
+            // and child are of same polynomial order.
+            fe->GetLocalRestriction(isotr, localRVO);
+            // copy block
+            auto s = localRVO.Height() * localRVO.Width();
+            std::copy(localRVO.Data(), localRVO.Data() + s, bs_ptr);
+            bs_ptr += s;
+         }
+         DenseMatrix &lR =
+            fespace->IsVariableOrder() ? localRVO : localR[geom](emb.matrix);
+         max_rows = std::max(lR.Height(), max_rows);
+         auto size = lR.Height() * lR.Width();
+         fespace->elem_dof->GetRow(emb.parent, dofs);
+         if (fine_rank == MyRank)
+         {
+            // diagonal
+            old_elem_dof->GetRow(k, old_dofs);
+            MFEM_VERIFY(old_dofs.Size() == dofs.Size(),
+                        "Parent and child must have same #dofs.");
+            block_row_idcs_offsets[diag_idx + 1] =
+               block_row_idcs_offsets[diag_idx] + lR.Height();
+            block_col_idcs_offsets[diag_idx + 1] =
+               block_col_idcs_offsets[diag_idx] + lR.Width();
+
+            if (fespace->IsVariableOrder())
+            {
+               block_offsets[diag_idx] = var_offset;
+               var_offset += size;
+            }
+            else
+            {
+               block_offsets[diag_idx] = geom_offsets[geom] + size * emb.matrix;
+            }
+            for (int i = 0; i < lR.Height(); ++i, ++ridx)
+            {
+               if (!std::isfinite(lR(i, 0)))
+               {
+                  row_idcs[ridx] = INT_MAX;
+                  continue;
+               }
+               int r = dofs[i];
+               int m = (r >= 0) ? r : (-1 - r);
+               if (is_dg || !mark[m])
+               {
+                  row_idcs[ridx] = r;
+                  mark[m] = 1;
+               }
+               else
+               {
+                  row_idcs[ridx] = INT_MAX;
+               }
+            }
+            for (int i = 0; i < lR.Width(); ++i, ++cidx)
+            {
+               col_idcs[cidx] = old_dofs[i];
+            }
+            ++diag_idx;
+         }
+         else
+         {
+            // off-diagonal
+            auto &tmp = ks_map.at(k);
+            auto od_idx = tmp[0];
+            auto od_ridx = tmp[1];
+            block_off_diag_row_idcs_offsets[od_idx + 1] =
+               block_off_diag_row_idcs_offsets[od_idx] + lR.Height();
+            block_off_diag_col_offsets[od_idx] = od_ridx;
+            block_off_diag_widths[od_idx] = lR.Width();
+            recv_segment_idcs[od_idx] = tmp[2];
+
+            if (fespace->IsVariableOrder())
+            {
+               off_diag_block_offsets[od_idx] = var_offset;
+               var_offset += size;
+            }
+            else
+            {
+               off_diag_block_offsets[od_idx] =
+                  geom_offsets[geom] + size * emb.matrix;
+            }
+            for (int i = 0; i < lR.Height(); ++i, ++od_ridx)
+            {
+               if (!std::isfinite(lR(i, 0)))
+               {
+                  row_off_diag_idcs[od_ridx] = INT_MAX;
+                  continue;
+               }
+               int r = dofs[i];
+               int m = (r >= 0) ? r : (-1 - r);
+               if (is_dg || !mark[m])
+               {
+                  row_off_diag_idcs[od_ridx] = r;
+                  mark[m] = 1;
+               }
+               else
+               {
+                  row_off_diag_idcs[od_ridx] = INT_MAX;
+               }
+            }
+            ++od_idx;
+         }
+      }
+   }
+   // if not using GPU, set max_rows/max_cols to zero
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      max_rows = std::min(max_rows, max_team_size);
+   }
+   else
+   {
+      max_rows = 1;
+   }
+   requests.reserve(recv_ranks.Size() + send_ranks.Size());
+}
+} // namespace mfem
+/// \endcond DO_NOT_DOCUMENT
+
+#endif
@@ -0,0 +1,111 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PDEREFMAT_OP
+#define MFEM_PDEREFMAT_OP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "pfespace.hpp"
+
+#include "kernel_dispatch.hpp"
+
+#include <vector>
+
+namespace mfem
+{
+
+/// \cond DO_NOT_DOCUMENT
+struct ParDerefineMatrixOp : public Operator
+{
+   ParFiniteElementSpace *fespace;
+   /// offsets into block_storage for diagonal
+   Array<int> block_offsets;
+   /// offsets into row_idcs for diagonal
+   Array<int> block_row_idcs_offsets;
+   /// offsets into col_idcs for diagonal
+   Array<int> block_col_idcs_offsets;
+
+   /// offsets into block_storage for off-diagonal
+   Array<int> off_diag_block_offsets;
+   /// offsets into row_idcs for off-diagonal
+   Array<int> block_off_diag_row_idcs_offsets;
+   Array<int> block_off_diag_col_offsets;
+   Array<int> block_off_diag_widths;
+   /// mapping for row dofs, INT_MAX indicates the block row should be ignored.
+   /// negative means the row data should be negated.
+   /// only for diagonal blocks
+   Array<int> row_idcs;
+   /// mapping for col dofs, negative means the col data should be negated.
+   /// only for diagonal blocks
+   Array<int> col_idcs;
+
+   Array<int> pack_col_idcs;
+
+   /// mapping for row dofs, INT_MAX indicates the block row should be ignored.
+   /// negative means the row data should be negated.
+   /// only for off-diagonal blocks
+   Array<int> row_off_diag_idcs;
+   /// dense block matrices which can be reused to construct the full matrix
+   /// operation. These are stored contiguously and blocks have no restrictions
+   /// on shape (can be rectangle and differ from block to block).
+   /// This is only for the diagonal block.
+   Vector block_storage;
+   /// maximum height of any block in block_storage for GPU
+   /// parallelization, or 1 for CPU runs.
+   int max_rows;
+
+   /// quasi Ordering::byNODES, broken into sections by ranks we need to send
+   /// the data to
+   mutable Vector xghost_send;
+   /// quasi Ordering::byNODES, broken into sections by ranks we received
+   /// the data from
+   mutable Vector xghost_recv;
+   /// maps off-diagonal k to segment
+   Array<int> recv_segment_idcs;
+   /// cumulative count of dofs which will be received from other ranks
+   Array<int> recv_segments;
+   /// Source rank of each recv segment
+   Array<int> recv_ranks;
+   /// What send segment each entry in send_permutations corresponds to
+   Array<int> send_segment_idcs;
+   /// cumulative count of dofs which will be sent to other ranks
+   Array<int> send_segments;
+   /// Destination rank of each send segment
+   Array<int> send_ranks;
+   /// how to permute/sign change values from our local x to send to other ranks
+   Array<int> send_permutations;
+   /// internal buffer for MPI requests
+   mutable std::vector<MPI_Request> requests;
+
+   using MultKernelType = void (*)(const ParDerefineMatrixOp &, const Vector &,
+                                   Vector &);
+   /// template args: ordering, atomic
+   MFEM_REGISTER_KERNELS(MultKernel, MultKernelType, (Ordering::Type, bool));
+
+   struct Kernels
+   {
+      Kernels();
+   };
+
+   void Mult(const Vector &x, Vector &y) const;
+
+   ParDerefineMatrixOp(ParFiniteElementSpace &fespace_, int old_ndofs,
+                       const Table *old_elem_dof, const Table *old_elem_fos);
+};
+/// \endcond DO_NOT_DOCUMENT
+} // namespace mfem
+
+#endif
+
+#endif
@@ -22,12 +22,13 @@
 #include "../mesh/mesh_headers.hpp"
 #include "../general/binaryio.hpp"

+#include "pderefmat_op.hpp"
+
 #include <limits>
 #include <list>

 namespace mfem
 {
-
 ParFiniteElementSpace::ParFiniteElementSpace(
   const ParFiniteElementSpace &orig, ParMesh *pmesh,
   const FiniteElementCollection *fec)
@@ -4487,13 +4488,6 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
   return M;
 }

-
-struct DerefDofMessage
-{
-   std::vector<HYPRE_BigInt> dofs;
-   MPI_Request request;
-};
-
 HypreParMatrix*
 ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
                                                  const Table* old_elem_dof,
@@ -4536,7 +4530,13 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
      old_pncmesh->GetDerefinementTransforms();
   const Array<int> &old_ranks = old_pncmesh->GetDerefineOldRanks();

-   std::map<int, DerefDofMessage> messages;
+   // key: other rank
+   // value: send or recieve buffer
+   std::map<int, std::vector<HYPRE_BigInt>> to_send;
+   std::map<int, std::vector<HYPRE_BigInt>> to_recv;
+   // key: index into dtrans.embeddings
+   // value: [start, stop]
+   std::unordered_map<int, std::array<size_t, 2>> recv_messages;

   HYPRE_BigInt old_offset = HYPRE_AssumedPartitionCheck()
                             ? old_dof_offsets[0] : old_dof_offsets[MyRank];
@@ -4556,30 +4556,46 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
         old_elem_dof->GetRow(k, dofs);
         DofsToVDofs(dofs, old_ndofs);

-         DerefDofMessage &msg = messages[k];
-         msg.dofs.resize(dofs.Size());
+         std::vector<HYPRE_BigInt>& send_buf = to_send[coarse_rank];
+         auto pos = send_buf.size();
+         send_buf.resize(pos + dofs.Size());
         for (int i = 0; i < dofs.Size(); i++)
         {
-            msg.dofs[i] = old_offset + dofs[i];
+            send_buf[pos + i] = old_offset + dofs[i];
         }
-
-         MPI_Isend(&msg.dofs[0], static_cast<int>(msg.dofs.size()), HYPRE_MPI_BIG_INT,
-                   coarse_rank, 291, MyComm, &msg.request);
      }
      else if (coarse_rank == MyRank && fine_rank != MyRank)
      {
         MFEM_ASSERT(emb.parent >= 0, "");
         Geometry::Type geom = mesh->GetElementBaseGeometry(emb.parent);

-         DerefDofMessage &msg = messages[k];
-         msg.dofs.resize(ldof[geom]*vdim);
-
-         MPI_Irecv(&msg.dofs[0], ldof[geom]*vdim, HYPRE_MPI_BIG_INT,
-                   fine_rank, 291, MyComm, &msg.request);
+         std::vector<HYPRE_BigInt>& recv_buf = to_recv[fine_rank];
+         auto& msg = recv_messages[k];
+         msg[0] = recv_buf.size();
+         recv_buf.resize(recv_buf.size() + ldof[geom] * vdim);
+         msg[1] = recv_buf.size();
      }
-      // TODO: coalesce Isends/Irecvs to the same rank. Typically, on uniform
-      // derefinement, there should be just one send to MyRank-1 and one recv
-      // from MyRank+1
+   }
+
+   // assume embedding orders are consistent (i.e. what we expect to receive
+   // first from a given rank is sent first, etc.)
+   std::vector<MPI_Request> requests;
+   requests.reserve(to_send.size() + to_recv.size());
+   // enqueue recvs
+   for (auto &v : to_recv)
+   {
+      requests.emplace_back();
+      MPI_Irecv(v.second.data(), v.second.size(), HYPRE_MPI_BIG_INT, v.first,
+                MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA, MyComm,
+                &requests.back());
+   }
+   // enqueue sends
+   for (auto &v : to_send)
+   {
+      requests.emplace_back();
+      MPI_Isend(v.second.data(), v.second.size(), HYPRE_MPI_BIG_INT, v.first,
+                MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA, MyComm,
+                &requests.back());
   }

   DenseTensor localR[Geometry::NumGeom];
@@ -4637,10 +4653,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
   diag->Finalize();

   // wait for all sends/receives to complete
-   for (auto it = messages.begin(); it != messages.end(); ++it)
-   {
-      MPI_Wait(&it->second.request, MPI_STATUS_IGNORE);
-   }
+   MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);

   // create the off-diagonal part of the derefinement matrix
   SparseMatrix *offd = new SparseMatrix(ndofs*vdim, 1);
@@ -4661,13 +4674,14 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,

         elem_dof->GetRow(emb.parent, dofs);

-         DerefDofMessage &msg = messages[k];
-         MFEM_ASSERT(msg.dofs.size(), "");
+         auto& odofs = to_recv.at(fine_rank);
+         auto &msg = recv_messages[k];
+         MFEM_ASSERT(msg[1] > msg[0], "");

         for (int vd = 0; vd < vdim; vd++)
         {
            MFEM_ASSERT(ldof[geom], "");
-            HYPRE_BigInt* remote_dofs = &msg.dofs[vd*ldof[geom]];
+            HYPRE_BigInt *remote_dofs = odofs.data() + msg[0] + vd * ldof[geom];

            for (int i = 0; i < lR.Height(); i++)
            {
@@ -4694,7 +4708,6 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
      }
   }

-   messages.clear();
   offd->Finalize(0);
   offd->SetWidth(static_cast<int>(col_map.size()));

@@ -4946,8 +4959,13 @@ void ParFiniteElementSpace::Update(bool want_transform)

         case Mesh::DEREFINE:
         {
+#if 0
            Th.Reset(ParallelDerefinementMatrix(old_ndofs, old_elem_dof,
                                                old_elem_fos));
+#else
+            Th.Reset(new ParDerefineMatrixOp(*this, old_ndofs, old_elem_dof,
+                                             old_elem_fos));
+#endif
            if (Nonconforming())
            {
               Th.SetOperatorOwner(false);
@@ -24,9 +24,12 @@
 namespace mfem
 {

+struct ParDerefineMatrixOp;
+
 /// Abstract parallel finite element space.
 class ParFiniteElementSpace : public FiniteElementSpace
 {
+   friend struct ParDerefineMatrixOp;
 private:
   /// MPI data.
   MPI_Comm MyComm;
@@ -425,16 +425,24 @@ public:
   ~GroupCommunicator();
 };

+/// General MPI message tags used by MFEM
+enum MessageTag
+{
+   DEREFINEMENT_MATRIX_CONSTRUCTION_DATA =
+      291, /// ParFiniteElementSpace ParallelDerefinementMatrix and
+   /// ParDerefineMatrixOp
+};
+
 enum VarMessageTag
 {
-   NEIGHBOR_ELEMENT_RANK_VM,  ///< NeighborElementRankMessage
-   NEIGHBOR_ORDER_VM,         ///< NeighborOrderMessage
-   NEIGHBOR_DEREFINEMENT_VM,  ///< NeighborDerefinementMessage
-   NEIGHBOR_REFINEMENT_VM,    ///< NeighborRefinementMessage
-   NEIGHBOR_PREFINEMENT_VM,   ///< NeighborPRefinementMessage
-   NEIGHBOR_ROW_VM,           ///< NeighborRowMessage
-   REBALANCE_VM,              ///< RebalanceMessage
-   REBALANCE_DOF_VM           ///< RebalanceDofMessage
+   NEIGHBOR_ELEMENT_RANK_VM, ///< NeighborElementRankMessage
+   NEIGHBOR_ORDER_VM,        ///< NeighborOrderMessage
+   NEIGHBOR_DEREFINEMENT_VM, ///< NeighborDerefinementMessage
+   NEIGHBOR_REFINEMENT_VM,   ///< NeighborRefinementMessage
+   NEIGHBOR_PREFINEMENT_VM,  ///< NeighborPRefinementMessage
+   NEIGHBOR_ROW_VM,          ///< NeighborRowMessage
+   REBALANCE_VM,             ///< RebalanceMessage
+   REBALANCE_DOF_VM,         ///< RebalanceDofMessage
 };

 /// \brief Variable-length MPI message containing unspecific binary data.
@@ -615,7 +615,7 @@ void Device::Setup(const std::string &device_option, const int device_id)
   if (Allows(Backend::DEBUG_DEVICE)) { ngpu = 1; }
 }

-MemoryType Device::QueryMemoryType(void *ptr)
+MemoryType Device::QueryMemoryType(const void* ptr)
 {
   // from HYPRE's hypre_GetPointerLocation
   MemoryType res = MemoryType::HOST;
@@ -297,9 +297,9 @@ public:
   /// Get the status of GPU-aware MPI flag.
   static bool GetGPUAwareMPI() { return Get().mpi_gpu_aware; }

-   /** @brief Query the device driver for what memory type a given @a ptr is
-       allocated with. */
-   static MemoryType QueryMemoryType(void *ptr);
+   /** Query the device driver for what memory type a given @a ptr is allocated
+    * with. */
+   static MemoryType QueryMemoryType(const void* ptr);

   /** @brief The number of hardware compute units/streaming multiprocessors
       available on a given compute device @a device_id. */
@@ -13257,7 +13257,7 @@ void Mesh::ScaleElements(real_t sf)
   delete [] vn;
 }

-void Mesh::Transform(void (*f)(const Vector&, Vector&))
+void Mesh::Transform(std::function<void(const Vector &, Vector&)> f)
 {
   // TODO: support for different new spaceDim.
   if (Nodes == NULL)
@@ -13270,7 +13270,7 @@ void Mesh::Transform(void (*f)(const Vector&, Vector&))
            vold(j) = vertices[i](j);
         }
         vnew.SetData(vertices[i]());
-         (*f)(vold, vnew);
+         f(vold, vnew);
      }
   }
   else
@@ -2254,7 +2254,7 @@ public:
   void ScaleSubdomains (real_t sf);
   void ScaleElements (real_t sf);

-   void Transform(void (*f)(const Vector&, Vector&));
+   void Transform(std::function<void(const Vector &, Vector&)> f);
   void Transform(VectorCoefficient &deformation);

   /** @brief This function should be called after the mesh node coordinates
@@ -105,8 +105,11 @@ MFEM_PERF_CXXFLAGS_xlc = -mcpu=native
 # - Clang extra options:
 ifeq ($(MFEM_MACHINE),riscv64)
  MFEM_PERF_CXXFLAGS_clang += -march=rv64gc
-else ifneq ($(MFEM_MACHINE),arm64)
-  # -march=native is unavailable on clang/ARM64 as of 05/2021: support could be added later.
+else ifneq (,$(findstring ppc,$(MFEM_MACHINE)))
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else ifeq ($(MFEM_MACHINE),arm64)
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else
  MFEM_PERF_CXXFLAGS_clang += -march=native
 endif
 MFEM_PERF_CXXFLAGS_clang += $(PEDANTIC_FLAG) -Wall
@@ -235,7 +235,7 @@ void test_derefine_L2_element(int order, Element::Type el_type, int basis_type)
   }
 }

-TEST_CASE("AMR Coarsen L2 Element","[AMR][Coarsen]")
+TEST_CASE("AMR Coarsen L2 Element","[AMR][Coarsen][CUDA]")
 {
   std::vector<int> orders_1d{0,1,2,3};
   std::vector<int> orders_2d{0,1,2,3};
@@ -420,7 +420,7 @@ void stress_parallel_coarsen(int order, Element::Type el_type, int basis_type)
 #endif
 }

-TEST_CASE("Parallel AMR Coarsen Stress Test", "[AMR][Coarsen][Parallel]")
+TEST_CASE("Parallel AMR Coarsen Stress Test", "[AMR][Coarsen][Parallel][CUDA]")
 {
   std::vector<int> orders_1d{0,1,2,3};
   std::vector<int> orders_2d{0,1,2,3};
@@ -482,7 +482,7 @@ TEST_CASE("Parallel AMR Coarsen Stress Test", "[AMR][Coarsen][Parallel]")

 }

-TEST_CASE("ParDerefine", "[Parallel]")
+TEST_CASE("ParDerefine", "[Parallel][CUDA]")
 {
   for (dimension = 2; dimension <= 3; ++dimension)
   {
Author	SHA1	Message	Date
Socratis Petrides	2f3e1b9c61	fix year	2025-08-07 11:18:41 -07:00
Socratis Petrides	442e793540	merge master and resolve conflicts	2025-08-07 10:59:40 -07:00
Veselin Dobrev	bd7281bdfa	Merge pull request #4960 from mfem/fix-warnings-algoim Algoim compilation warnings	2025-08-05 12:05:40 -07:00
Tzanio Kolev	6dfbe25c07	Merge pull request #4955 from mfem/clang-perf-flags In `miniapps/performance/makefile`, adjust the flags for Clang	2025-08-04 07:44:47 -07:00
Tzanio Kolev	02edb95a31	Merge pull request #4957 from farscape-project/conduit Guarantee data on host when writing conduit files	2025-08-02 09:27:45 -07:00
Will Pazner	597cba8d37	Merge pull request #4954 from mfem/ruby-60 Extend time allocation for testing on ruby from 45 to 60 minutes	2025-07-29 12:49:00 -07:00
Veselin Dobrev	075a160187	Merge pull request #4827 from mfem/feature/guthrey1/transform_lambda Allow capturing lambdas for mesh::Transform	2025-07-29 09:34:22 -07:00
Vladimir Z Tomov	8d9579057f	Simplified to std:vector.	2025-07-28 15:55:21 -07:00
Vladimir Z Tomov	76d65f5866	Fixed compilation warnings when Algoim is on.	2025-07-28 15:39:50 -07:00
Tzanio Kolev	720fb72372	Merge pull request #4884 from adam-sim-dev/gitignore Fix gitignore for Windows	2025-07-26 10:45:44 -07:00
Nuno Nobre	bbde976378	Guarantee data on host when writing conduit files	2025-07-26 01:31:22 +01:00
Tzanio Kolev	58ad1c79cd	Merge pull request #4673 from mfem/cmake-replace-library-files-with-targets Support fetching of third party libraries using CMake's ExternalProject	2025-07-24 12:23:03 -07:00
Veselin Dobrev	447cd0f500	In miniapps/performance/makefile, adjust the flags for Clang	2025-07-22 22:14:35 -07:00
adam-sim-dev	2d68b877f2	Merge branch 'master' into gitignore	2025-07-23 10:01:26 +08:00
Tzanio Kolev	cce301f92b	Extend time allocation for testing on ruby from 45 to 60 minutes	2025-07-22 12:45:09 -07:00
Veselin Dobrev	46ccaee493	Merge branch 'master' into feature/guthrey1/transform_lambda	2025-07-21 19:36:55 -07:00
Tzanio Kolev	7aa327174a	Merge pull request #4811 from mfem/dev-derefine AMR Matrix-free derefinement operator	2025-07-15 12:22:19 -07:00
adam-sim-dev	8c861b4602	Merge branch 'master' into gitignore	2025-07-15 09:20:33 +08:00
Andrew Ho	2d5f00635f	fixed derefinement for ordering byVDIM	2025-07-10 17:00:49 -07:00
Arturo Vargas	e629502dca	Merge branch 'master' into dev-derefine	2025-07-07 09:55:58 -07:00
Arturo Vargas	48dd78b154	Merge branch 'master' into dev-derefine	2025-07-03 09:23:21 -07:00
Andrew Ho	dc08fc65ce	updated changelog	2025-07-01 16:01:38 -07:00
Andrew Ho	3276cb955f	Merge branch 'master' into dev-derefine	2025-07-01 15:33:44 -07:00
Andrew Ho	bcdbda1597	Merge branch 'master' into dev-derefine	2025-07-01 11:24:42 -07:00
adam-sim-dev	0239b86304	Merge branch 'master' into gitignore	2025-06-30 12:15:37 +08:00
Andrew Ho	aca8d51f95	Merge branch 'master' into dev-derefine	2025-06-27 09:25:25 -07:00
Andrew Ho	368724d456	fix missing climit header	2025-06-26 09:30:45 -07:00
Tzanio Kolev	a99efcdb63	Merge branch 'master' into dev-derefine	2025-06-26 08:25:48 -07:00
adam-sim-dev	b5a7cf6470	Merge branch 'master' into gitignore	2025-06-24 12:11:35 +08:00
Andrew Ho	e6fc7893be	Merge branch 'master' into dev-derefine	2025-06-16 12:20:13 -07:00
adam-sim-devandJan Nikl	78df8dffb2	Update Co-authored-by: Jan Nikl <nikl1@llnl.gov>	2025-06-10 05:04:44 +08:00
adam-sim-dev	d0f5cb73e7	Fix for Windows	2025-06-09 12:17:58 +08:00
Andrew Ho	db55be84ad	Merge branch 'master' into dev-derefine	2025-06-02 11:37:29 -07:00
Andrew Ho	091c6272fc	Merge branch 'master' into dev-derefine	2025-05-27 11:57:28 -07:00
Andrew Ho	101af72047	Merge branch 'master' into dev-derefine	2025-05-19 12:17:39 -07:00
Andrew Ho	01990d4605	don't need _VM	2025-05-13 10:32:01 -07:00
Andrew Ho	4a3cd44695	move tag from VarMessageTag to MessageTag	2025-05-13 04:46:00 -07:00
Andrew Ho	3b059c3995	struct instead of class	2025-05-12 17:21:05 -07:00
Andrew Ho	89dc5ce5c0	parallel version needs to compile to nothing if not building with MPI	2025-05-12 17:19:16 -07:00
Andrew Ho	eb205caaa2	fixed collision in ifdef guard	2025-05-12 17:08:39 -07:00
Andrew Ho	f50601f1aa	fixed using struct vs. class	2025-05-12 16:52:32 -07:00
Andrew Ho	529d6e3bc2	Split deref matrix ops into separate files	2025-05-12 16:34:40 -07:00
Andrew Ho	ef29b1dd30	Merge remote-tracking branch 'base/dev-derefine' into dev-derefine	2025-05-12 11:58:50 -07:00
Andrew Ho	33b3ef6045	Enable GPU testing of fespace/pfespace derefinement	2025-05-12 11:40:26 -07:00
Andrew Ho	657f6ef5fe	Merge branch 'master' into dev-derefine	2025-05-12 11:19:22 -07:00
Andrew Ho	96cfef783c	Review suggestions - const on some local variables - use VarMessageTag instead of hard-coded integers for MPI messages - comment updates	2025-05-12 10:59:23 -07:00
Andrew Ho	603f43e775	updated fes_kernels member variable documentation	2025-05-12 10:47:13 -07:00
Andrew Ho	a2e33109eb	Merge branch 'master' into dev-derefine	2025-05-02 13:04:12 -06:00
Andrew Ho	b277a0cb2d	Merge branch 'master' into dev-derefine	2025-05-02 10:29:36 -06:00
Tzanio Kolev	ac26026ff1	Merge branch 'master' into cmake-replace-library-files-with-targets	2025-04-26 12:19:30 -07:00
Pierson Guthrey	d439781751	resolve ambiguity, thanks @v-dobrev	2025-04-23 18:08:06 -07:00
Pierson Guthrey	16c9123f98	allow lambdas for mesh Transformation	2025-04-23 17:21:16 -07:00
Andrew Ho	e73ea1fb77	documentation for parallel version	2025-04-23 12:11:54 -07:00
Andrew Ho	ea6c104475	added const and documentation	2025-04-23 12:05:21 -07:00
Andrew Ho	e7762d126e	Merge branch 'master' into dev-derefine	2025-04-23 09:02:36 -07:00
Andrew Ho	c61d0e07d3	switch to using INT_MAX for rows to skip	2025-04-22 16:51:04 -07:00
Andrew Ho	720051d641	unused testing var	2025-04-22 08:55:41 -07:00
Andrew Ho	5c25fa49a7	graceful fallback for HOST_PINNED	2025-04-22 08:50:23 -07:00
Andrew Ho	d9ef117ee9	Merge branch 'master' into dev-derefine	2025-04-21 23:21:29 -07:00
Andrew Ho	b7a8cf6290	CI fixes	2025-04-21 23:20:43 -07:00
Andrew Ho	bed1823d72	unused variables	2025-04-21 23:12:54 -07:00
Andrew Ho	93a3648712	fixed a few performance bugs	2025-04-21 20:12:35 -07:00
Christopher vogl	01eeccee7f	updated hypre fetch version; removed unnecessary dependencies on install step; leverage METIS_FETCH_VERSION variable	2025-04-21 15:31:52 -07:00
Christopher vogl	39b6112e69	reverted use of MPI target instead of files	2025-04-21 13:49:18 -07:00
Andrew Ho	c70d2b7309	off-diagonals should be implemented now	2025-04-21 00:34:35 -07:00
Andrew Ho	c98e8bfe07	Think I have ghost comm data being transferred correctly now	2025-04-20 23:03:10 -07:00
Andrew Ho	ef9137bd1f	const correctness fix	2025-04-20 22:59:08 -07:00
Andrew Ho	8e0c42a132	fixed a few bugs still need to finish implementing off-diagonals	2025-04-19 23:04:11 -07:00
Andrew Ho	d755037904	Re-arranged DerefineMatrixOpFunctor Now it's a common implementation with the indexing order separated. Still need to make it work for off-diagonals.	2025-04-19 20:42:17 -07:00
Andrew Ho	b577d2322e	diagonals for pfespace derefop working still need to work out the off-diagonals	2025-04-19 17:46:10 -07:00
Andrew Ho	6f99204d37	code cleanup	2025-04-19 01:32:52 -07:00
Andrew Ho	c6e93f5474	Coalesce ParallelDerefineMatrix MPI messages Not strictly required for correctness, but fewer MPI messages should always be good	2025-04-18 22:17:24 -07:00
Andrew Ho	21b6ade630	removed derefine multT I don't think this is ever used or has any meaning	2025-04-17 14:31:44 -07:00
Andrew Ho	2e8685d80a	started working on ParDerefineMatrixOp	2025-04-17 11:59:01 -07:00
Andrew Ho	6d8c32db80	only need atomic for mult with DG	2025-04-17 09:01:07 -07:00
Andrew Ho	4adb51df50	fixed documentation and p derefinement oob read	2025-04-16 14:37:47 -07:00
Andrew Ho	6cec49e220	fix compiler warnings	2025-04-16 14:23:35 -07:00
Andrew Ho	9006171327	Merge remote-tracking branch 'base/master' into dev-derefine	2025-04-16 13:18:46 -07:00
Andrew Ho	9e10cd87b1	fixed verify for vdim > 1	2025-04-15 23:15:23 -07:00
Andrew Ho	3f7c551d66	implemented MultTranspose	2025-04-15 22:45:35 -07:00
Andrew Ho	a7a5ce5cb1	implemented Mult for VDim ordering	2025-04-15 21:48:53 -07:00
Andrew Ho	4218c11b41	have at least one kernel working	2025-04-15 21:31:55 -07:00
Andrew Ho	3c4a860f3a	re-arranged the run operation so functor is responsible need to change parallelization so vdim can be allocated into shared memory as well	2025-04-15 10:55:39 -07:00
Andrew Ho	e64300f84b	Re-arranged DerefineMatrixOp to use kernel specialization infrastructure all that's left is to implement Mult and MultT for specializations	2025-04-10 13:31:26 -07:00
Andrew Ho	3854184620	have setup for prototype derefinematrixop implemented still not sure if DerefinementOperator is acceptable/modifiable or not for this	2025-04-09 16:22:13 -07:00
Christopher vogl	bfc6b0708d	corrected some unncecessary whitespace added	2025-02-20 16:48:54 -08:00
Christopher vogl	3d5c6d58e8	added options to defaults.cmake and updated INSTALL document	2025-02-20 16:45:40 -08:00
Christopher vogl	6efbd4cc27	applied change suggested by nmnobre to use patched metis to support ARM macs	2025-02-20 16:18:35 -08:00
Christopher vogl	369b9365ad	added check for hypre target to top of FindHYPRE	2025-02-20 12:35:19 -08:00
Christopher vogl	b35ff2270a	applying nmnobre's suggestion to clean up fetched metis directory	2025-02-20 10:24:42 -08:00
Chris VoglandNuno Nobre	7a54a543b0	Apply update_disconnect fix Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2025-02-20 10:06:24 -08:00
Christopher vogl	15420f7049	restored support for FETCH_TPLS	2025-02-20 09:56:32 -08:00
Chris VoglandNuno Nobre	dd20722733	Applying suggestions for individual fetching and dynamic hypre versioning Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2025-02-20 09:42:50 -08:00
Christopher vogl	5577edf446	Added status messages suggested by nmnobre	2025-02-19 09:08:40 -08:00
Chris VoglandNuno Nobre	230a345ef3	Apply suggestions from nmnobre 1. use git tag instead of commit 2. explicitly specify hypre library install directory Co-authored-by: Nuno Nobre <nuno.nobre@stfc.ac.uk>	2025-02-19 08:42:45 -08:00
Christopher vogl	89f1dce705	added support for fetching hypre with CUDA, HIP, and floats	2025-01-23 10:26:20 -08:00
Christopher vogl	7f17c1f049	see prior message (clearly I need more coffee)	2025-01-22 17:48:30 -08:00
Christopher vogl	7e6edcddec	changes that should have been part of prior commit (more commenting)	2025-01-22 17:47:38 -08:00
Christopher vogl	d89b01e24c	updated commenting of new code	2025-01-22 17:46:04 -08:00
Christopher vogl	33f1150c34	updated FindMETIS to also use ExternalProject	2025-01-22 16:14:09 -08:00
Christopher vogl	c6a05d371b	fixed capitalization and made fetch directory consistent	2025-01-22 16:13:41 -08:00
Christopher vogl	f97a29f6e8	refactored FindHYPRE to use ExternalProject instead of FetchContent	2025-01-22 15:00:06 -08:00
Christopher vogl	85ae72fada	changed lib to Lib to be consistent with METIS capitalization	2025-01-21 16:16:30 -08:00
Christopher vogl	27a5bd71ee	whitespace	2025-01-21 15:27:14 -08:00
Christopher vogl	643db22bef	minor tweaks, commenting, and cleanup	2025-01-21 15:25:18 -08:00
Christopher vogl	97fb6be029	Merge remote-tracking branch 'origin/master' into cmake-replace-library-files-with-targets	2025-01-21 13:41:24 -08:00
Christopher vogl	9587ef212d	added fetching of HYPRE and METIS with FETCH_TPLS flag	2025-01-21 13:41:08 -08:00
Christopher vogl	35462b8645	refactored CMakeLists to support targets for TPLs	2025-01-15 13:13:18 -08:00
Socratis Petrides	461208444a	fix double free	2024-04-12 18:42:26 -07:00
Socratis Petrides	7c78c34af6	adding convinince class for blockforms	2024-04-12 15:38:55 -07:00