Merge branch 'master' into delete-indices-dev

Use `HipMemcpyDtoD` instead of `HipMemcpyDtoDAsync`

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

@@ -289,10 +289,10 @@ jobs:
         run: |
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
-          brew install llvm@19 enzyme
-          echo "LLVM_PREFIX=$(brew --prefix llvm@19)" >> $GITHUB_ENV
-          echo "OMPI_CC=$(brew --prefix llvm@19)/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$(brew --prefix llvm@19)/bin/clang++" >> $GITHUB_ENV
+          brew install llvm@20 enzyme
+          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
+          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
 
       # MFEM build and test
       - name: build

.gitignore

@@ -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

.gitlab/configs/ruby-config.yml

@@ -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

CHANGELOG

@@ -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
 -------------------------------------

CMakeLists.txt

@@ -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()

INSTALL

@@ -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>
@@ -842,6 +847,7 @@ The specific libraries and their options are:
 - HIP (optional), used when MFEM_USE_HIP = YES.
   URL: https://rocmdocs.amd.com
   Options: HIP_CXX, HIP_ARCH, HIP_OPT, HIP_LIB.
+  Versions: ROCm >= 5.6.1.
 
 - OCCA (optional), used when MFEM_USE_OCCA = YES.
   URL: https://libocca.org
@@ -1074,6 +1080,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 +1144,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
 ==================================

config/cmake/modules/FindHYPRE.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."

config/cmake/modules/FindMETIS.cmake

@@ -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"

config/config.hpp

@@ -27,7 +27,7 @@ namespace mfem
 {
 
 #if (defined(MFEM_USE_CUDA) && defined(__CUDACC__)) || \
-    (defined(MFEM_USE_HIP) && defined(__HIPCC__))
+    (defined(MFEM_USE_HIP) && defined(__HIP__))
 #define MFEM_HOST_DEVICE __host__ __device__
 #else
 #define MFEM_HOST_DEVICE

config/defaults.cmake

@@ -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++)

examples/ex22.cpp

@@ -62,9 +62,14 @@ static real_t epsilon_ = 1.0;
 static real_t sigma_ = 20.0;
 static real_t omega_ = 10.0;
 
-complex<real_t> u0_exact(const Vector &x);
-void u1_exact(const Vector &, ComplexVector &);
-void u2_exact(const Vector &, ComplexVector &);
+real_t u0_real_exact(const Vector &);
+real_t u0_imag_exact(const Vector &);
+
+void u1_real_exact(const Vector &, Vector &);
+void u1_imag_exact(const Vector &, Vector &);
+
+void u2_real_exact(const Vector &, Vector &);
+void u2_imag_exact(const Vector &, Vector &);
 
 bool check_for_inline_mesh(const char * mesh_file);
 
@@ -210,48 +215,54 @@ int main(int argc, char *argv[])
    ComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
 
-   ComplexFunctionCoefficient u0(u0_exact);
-   ComplexVectorFunctionCoefficient u1(dim, u1_exact);
-   ComplexVectorFunctionCoefficient u2(dim, u2_exact);
+   FunctionCoefficient u0_r(u0_real_exact);
+   FunctionCoefficient u0_i(u0_imag_exact);
+   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
+   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
+   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
+   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
 
-   ComplexConstantCoefficient oneCoef(1.0);
+   ConstantCoefficient zeroCoef(0.0);
+   ConstantCoefficient oneCoef(1.0);
 
+   Vector zeroVec(dim); zeroVec = 0.0;
    Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
-   ComplexVectorConstantCoefficient oneVecCoef(oneVec);
+   VectorConstantCoefficient zeroVecCoef(zeroVec);
+   VectorConstantCoefficient oneVecCoef(oneVec);
 
    switch (prob)
    {
       case 0:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficient(u0, ess_bdr);
-            u_exact->ProjectCoefficient(u0);
+            u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
+            u_exact->ProjectCoefficient(u0_r, u0_i);
          }
          else
          {
-            u.ProjectBdrCoefficient(oneCoef, ess_bdr);
+            u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
          }
          break;
       case 1:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientTangent(u1, ess_bdr);
-            u_exact->ProjectCoefficient(u1);
+            u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
+            u_exact->ProjectCoefficient(u1_r, u1_i);
          }
          else
          {
-            u.ProjectBdrCoefficientTangent(oneVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
          }
          break;
       case 2:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientNormal(u2, ess_bdr);
-            u_exact->ProjectCoefficient(u2);
+            u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
+            u_exact->ProjectCoefficient(u2_r, u2_i);
          }
          else
          {
-            u.ProjectBdrCoefficientNormal(oneVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
          }
          break;
       default: break; // This should be unreachable
@@ -289,24 +300,27 @@ int main(int argc, char *argv[])
    ConstantCoefficient lossCoef(omega_ * sigma_);
    ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
 
-   ComplexConstantCoefficient complexMassCoef(-omega_ * omega_ * epsilon_,
-                                              omega_ * sigma_);
-
    SesquilinearForm *a = new SesquilinearForm(fespace, conv);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    switch (prob)
    {
       case 0:
-         a->AddDomainIntegrator<DiffusionIntegrator>(stiffnessCoef);
-         a->AddDomainIntegrator<MassIntegrator>(complexMassCoef);
+         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new MassIntegrator(massCoef),
+                                new MassIntegrator(lossCoef));
          break;
       case 1:
-         a->AddDomainIntegrator<CurlCurlIntegrator>(stiffnessCoef);
-         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
+         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
          break;
       case 2:
-         a->AddDomainIntegrator<DivDivIntegrator>(stiffnessCoef);
-         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
+         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
          break;
       default: break; // This should be unreachable
    }
@@ -422,24 +436,29 @@ int main(int argc, char *argv[])
 
    if (exact_sol)
    {
-      real_t err_u = -1.0;
+      real_t err_r = -1.0;
+      real_t err_i = -1.0;
 
       switch (prob)
       {
          case 0:
-            err_u = u.ComputeL2Error(u0);
+            err_r = u.real().ComputeL2Error(u0_r);
+            err_i = u.imag().ComputeL2Error(u0_i);
             break;
          case 1:
-            err_u = u.ComputeL2Error(u1);
+            err_r = u.real().ComputeL2Error(u1_r);
+            err_i = u.imag().ComputeL2Error(u1_i);
             break;
          case 2:
-            err_u = u.ComputeL2Error(u2);
+            err_r = u.real().ComputeL2Error(u2_r);
+            err_i = u.imag().ComputeL2Error(u2_i);
             break;
          default: break; // This should be unreachable
       }
 
       cout << endl;
-      cout << "|| u_h - u ||_{L^2} = " << err_u << endl;
+      cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
+      cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
       cout << endl;
    }
 
@@ -545,14 +564,36 @@ complex<real_t> u0_exact(const Vector &x)
    return std::exp(-i * kappa * x[dim - 1]);
 }
 
-void u1_exact(const Vector &x, ComplexVector &v)
+real_t u0_real_exact(const Vector &x)
 {
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_exact(x);
+   return u0_exact(x).real();
 }
 
-void u2_exact(const Vector &x, ComplexVector &v)
+real_t u0_imag_exact(const Vector &x)
+{
+   return u0_exact(x).imag();
+}
+
+void u1_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_exact(x);
+   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+}
+
+void u1_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
+}
+
+void u2_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
+}
+
+void u2_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
 }

examples/ex22p.cpp

@@ -62,10 +62,6 @@ static real_t epsilon_ = 1.0;
 static real_t sigma_ = 20.0;
 static real_t omega_ = 10.0;
 
-complex<real_t> u0_exact(const Vector &x);
-void u1_exact(const Vector &, ComplexVector &);
-void u2_exact(const Vector &, ComplexVector &);
-
 real_t u0_real_exact(const Vector &);
 real_t u0_imag_exact(const Vector &);
 
@@ -248,22 +244,13 @@ int main(int argc, char *argv[])
    ParComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
 
-   ComplexFunctionCoefficient u0(u0_exact);
-   ComplexVectorFunctionCoefficient u1(dim, u1_exact);
-   ComplexVectorFunctionCoefficient u2(dim, u2_exact);
-
-   ComplexConstantCoefficient oneCoef(1.0);
-
-   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
-   ComplexVectorConstantCoefficient oneVecCoef(oneVec);
-
    FunctionCoefficient u0_r(u0_real_exact);
    FunctionCoefficient u0_i(u0_imag_exact);
    VectorFunctionCoefficient u1_r(dim, u1_real_exact);
    VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
    VectorFunctionCoefficient u2_r(dim, u2_real_exact);
    VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
-   /*
+
    ConstantCoefficient zeroCoef(0.0);
    ConstantCoefficient oneCoef(1.0);
 
@@ -271,40 +258,40 @@ int main(int argc, char *argv[])
    Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
    VectorConstantCoefficient zeroVecCoef(zeroVec);
    VectorConstantCoefficient oneVecCoef(oneVec);
-   */
+
    switch (prob)
    {
       case 0:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficient(u0, ess_bdr);
-            u_exact->ProjectCoefficient(u0);
+            u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
+            u_exact->ProjectCoefficient(u0_r, u0_i);
          }
          else
          {
-            u.ProjectBdrCoefficient(oneCoef, ess_bdr);
+            u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
          }
          break;
       case 1:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientTangent(u1, ess_bdr);
-            u_exact->ProjectCoefficient(u1);
+            u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
+            u_exact->ProjectCoefficient(u1_r, u1_i);
          }
          else
          {
-            u.ProjectBdrCoefficientTangent(oneVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
          }
          break;
       case 2:
          if (exact_sol)
          {
-            u.ProjectBdrCoefficientNormal(u2, ess_bdr);
-            u_exact->ProjectCoefficient(u2);
+            u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
+            u_exact->ProjectCoefficient(u2_r, u2_i);
          }
          else
          {
-            u.ProjectBdrCoefficientNormal(oneVecCoef, ess_bdr);
+            u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
          }
          break;
       default: break; // This should be unreachable
@@ -344,24 +331,27 @@ int main(int argc, char *argv[])
    ConstantCoefficient lossCoef(omega_ * sigma_);
    ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
 
-   ComplexConstantCoefficient complexMassCoef(-omega_ * omega_ * epsilon_,
-                                              omega_ * sigma_);
-
    ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    switch (prob)
    {
       case 0:
-         a->AddDomainIntegrator<DiffusionIntegrator>(stiffnessCoef);
-         a->AddDomainIntegrator<MassIntegrator>(complexMassCoef);
+         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new MassIntegrator(massCoef),
+                                new MassIntegrator(lossCoef));
          break;
       case 1:
-         a->AddDomainIntegrator<CurlCurlIntegrator>(stiffnessCoef);
-         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
+         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
          break;
       case 2:
-         a->AddDomainIntegrator<DivDivIntegrator>(stiffnessCoef);
-         a->AddDomainIntegrator<VectorFEMassIntegrator>(complexMassCoef);
+         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
          break;
       default: break; // This should be unreachable
    }
@@ -485,18 +475,22 @@ int main(int argc, char *argv[])
 
    if (exact_sol)
    {
-      real_t err_u = -1.0;
+      real_t err_r = -1.0;
+      real_t err_i = -1.0;
 
       switch (prob)
       {
          case 0:
-            err_u = u.ComputeL2Error(u0);
+            err_r = u.real().ComputeL2Error(u0_r);
+            err_i = u.imag().ComputeL2Error(u0_i);
             break;
          case 1:
-            err_u = u.ComputeL2Error(u1);
+            err_r = u.real().ComputeL2Error(u1_r);
+            err_i = u.imag().ComputeL2Error(u1_i);
             break;
          case 2:
-            err_u = u.ComputeL2Error(u2);
+            err_r = u.real().ComputeL2Error(u2_r);
+            err_i = u.imag().ComputeL2Error(u2_i);
             break;
          default: break; // This should be unreachable
       }
@@ -504,7 +498,8 @@ int main(int argc, char *argv[])
       if ( myid == 0 )
       {
          cout << endl;
-         cout << "|| u_h - u ||_{L^2} = " << err_u << endl;
+         cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
+         cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
          cout << endl;
       }
    }
@@ -632,12 +627,6 @@ real_t u0_imag_exact(const Vector &x)
    return u0_exact(x).imag();
 }
 
-void u1_exact(const Vector &x, ComplexVector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[0] = u0_exact(x);
-}
-
 void u1_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();
@@ -650,12 +639,6 @@ void u1_imag_exact(const Vector &x, Vector &v)
    v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
 }
 
-void u2_exact(const Vector &x, ComplexVector &v)
-{
-   int dim = x.Size();
-   v.SetSize(dim); v = 0.0; v[dim-1] = u0_exact(x);
-}
-
 void u2_real_exact(const Vector &x, Vector &v)
 {
    int dim = x.Size();

fem/CMakeLists.txt

@@ -59,7 +59,6 @@ set(SRCS
   integ/nonlininteg_vecconvection_pa.cpp
   integ/nonlininteg_vecconvection_mf.cpp
   coefficient.cpp
-  complex_coefficient.cpp
   complex_fem.cpp
   convergence.cpp
   datacollection.cpp
@@ -83,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
@@ -177,7 +178,6 @@ set(HDRS
   integ/bilininteg_hcurlhdiv_kernels.hpp
   integ/bilininteg_mass_kernels.hpp
   coefficient.hpp
-  complex_coefficient.hpp
   complex_fem.hpp
   convergence.hpp
   datacollection.hpp

fem/bilininteg.hpp

@@ -812,7 +812,7 @@ protected:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -884,7 +884,7 @@ protected:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV  &&
+      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -919,7 +919,7 @@ protected:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV  &&
+      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -1600,7 +1600,7 @@ public:
    {
       return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -1635,7 +1635,7 @@ public:
    {
       return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 
@@ -1669,7 +1669,7 @@ public:
    {
       return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -1760,7 +1760,7 @@ public:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
    }
 
@@ -1793,7 +1793,7 @@ public:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::DIV   );
    }
@@ -1832,7 +1832,7 @@ public:
       const FiniteElement & test_fe) const
    {
       return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::DIV    &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::DIV   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::GRAD
              );
@@ -1973,7 +1973,7 @@ protected:
       const FiniteElement & test_fe) const override
    {
       return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
 

fem/complex_coefficient.cpp

@@ -1,217 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "complex_fem.hpp"
-#include "../general/forall.hpp"
-
-using namespace std;
-
-namespace mfem
-{
-
-real_t
-RealPartCoefficient::Eval(ElementTransformation &T,
-                          const IntegrationPoint &ip)
-{
-   complex_t val = complex_coef_.Eval(T, ip);
-   return val.real();
-}
-
-real_t
-ImagPartCoefficient::Eval(ElementTransformation &T,
-                          const IntegrationPoint &ip)
-{
-   complex_t val = complex_coef_.Eval(T, ip);
-   return val.imag();
-}
-
-RealPartVectorCoefficient::RealPartVectorCoefficient(ComplexVectorCoefficient &
-                                                     complex_vcoef)
-   : VectorCoefficient(complex_vcoef.GetVDim()),
-     complex_vcoef_(complex_vcoef),
-     val_(vdim)
-{}
-
-void
-RealPartVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                const IntegrationPoint &ip)
-{
-   complex_vcoef_.Eval(val_, T, ip);
-   V = val_.real();
-}
-
-ImagPartVectorCoefficient::ImagPartVectorCoefficient(ComplexVectorCoefficient &
-                                                     complex_vcoef)
-   : VectorCoefficient(complex_vcoef.GetVDim()),
-     complex_vcoef_(complex_vcoef),
-     val_(vdim)
-{}
-
-void
-ImagPartVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                const IntegrationPoint &ip)
-{
-   complex_vcoef_.Eval(val_, T, ip);
-   V = val_.imag();
-}
-
-RealPartMatrixCoefficient::RealPartMatrixCoefficient(ComplexMatrixCoefficient &
-                                                     complex_mcoef)
-   : MatrixCoefficient(complex_mcoef.GetHeight(), complex_mcoef.GetWidth()),
-     complex_mcoef_(complex_mcoef),
-     val_(height, width)
-{}
-
-void
-RealPartMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                const IntegrationPoint &ip)
-{
-   complex_mcoef_.Eval(val_, T, ip);
-   M = val_.real();
-}
-
-ImagPartMatrixCoefficient::ImagPartMatrixCoefficient(ComplexMatrixCoefficient &
-                                                     complex_mcoef)
-   : MatrixCoefficient(complex_mcoef.GetHeight(), complex_mcoef.GetWidth()),
-     complex_mcoef_(complex_mcoef),
-     val_(height, width)
-{}
-
-void
-ImagPartMatrixCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
-                                const IntegrationPoint &ip)
-{
-   complex_mcoef_.Eval(val_, T, ip);
-   M = val_.imag();
-}
-
-ComplexCoefficient::ComplexCoefficient()
-   : time(0.),
-     re_part_coef_(*this), im_part_coef_(*this),
-     real_coef_(re_part_coef_), imag_coef_(im_part_coef_)
-{ }
-
-ComplexCoefficient::ComplexCoefficient(Coefficient &c_r,
-                                       Coefficient &c_i)
-   : time(c_r.GetTime()),
-     re_part_coef_(*this), im_part_coef_(*this),
-     real_coef_(c_r), imag_coef_(c_i)
-{
-   c_i.SetTime(time);
-}
-
-complex_t
-ComplexCoefficient::Eval(ElementTransformation &T,
-                         const IntegrationPoint &ip)
-{
-   // Avoid circular dependency
-   MFEM_VERIFY(std::addressof(real_coef_) != std::addressof(re_part_coef_) &&
-               std::addressof(imag_coef_) != std::addressof(im_part_coef_),
-               "Classes dervied from ComplexCoefficient must either "
-               "implement an Eval method or supply Coefficients "
-               "for both the real and imaginary parts of the field.");
-
-   return complex_t(real_coef_.Eval(T, ip), imag_coef_.Eval(T, ip));
-}
-
-ComplexVectorCoefficient::ComplexVectorCoefficient(VectorCoefficient &v_r,
-                                                   VectorCoefficient &v_i)
-   : vdim(v_r.GetVDim()), time(v_r.GetTime()),
-     re_part_vcoef_(*this), im_part_vcoef_(*this),
-     real_vcoef_(v_r), imag_vcoef_(v_i)
-{
-   MFEM_ASSERT(v_r.GetVDim() == v_i.GetVDim(), "ComplexVectorCoefficient"
-               " - incompatible vector dimensions of real and imaginary parts.");
-
-   v_i.SetTime(time);
-}
-
-void ComplexVectorCoefficient::Eval(ComplexVector &V, ElementTransformation &T,
-                                    const IntegrationPoint &ip)
-{
-   // Avoid circular dependency
-   MFEM_VERIFY(std::addressof(real_vcoef_) != std::addressof(re_part_vcoef_) &&
-               std::addressof(imag_vcoef_) != std::addressof(im_part_vcoef_),
-               "Classes dervied from ComplexVectorCoefficient must either "
-               "implement an Eval method or supply VectorCoefficients "
-               "for both the real and imaginary parts of the field.");
-
-   V_r_.SetSize(vdim);
-   V_i_.SetSize(vdim);
-
-   real_vcoef_.Eval(V_r_, T, ip);
-   imag_vcoef_.Eval(V_i_, T, ip);
-
-   V.Set(V_r_, V_i_);
-}
-
-ComplexConstantCoefficient::ComplexConstantCoefficient(
-   const complex_t z)
-   : val(z), real_coef(z.real()), imag_coef(z.imag())
-{
-   real_coef_ = real_coef;
-   imag_coef_ = imag_coef;
-}
-
-ComplexConstantCoefficient::ComplexConstantCoefficient(
-   real_t z_r, real_t z_i)
-   : real_coef(z_r), imag_coef(z_i)
-{
-   val = complex_t(z_r, z_i);
-
-   real_coef_ = real_coef;
-   imag_coef_ = imag_coef;
-}
-
-complex_t ComplexFunctionCoefficient::Eval(ElementTransformation & T,
-                                           const IntegrationPoint & ip)
-{
-   real_t x[3];
-   Vector transip(x, 3);
-
-   T.Transform(ip, transip);
-
-   if (Function)
-   {
-      return Function(transip);
-   }
-   else
-   {
-      return TDFunction(transip, GetTime());
-   }
-}
-
-void ComplexVectorFunctionCoefficient::Eval(ComplexVector &V,
-                                            ElementTransformation &T,
-                                            const IntegrationPoint &ip)
-{
-   real_t x[3];
-   Vector transip(x, 3);
-
-   T.Transform(ip, transip);
-
-   V.SetSize(vdim);
-   if (Function)
-   {
-      Function(transip, V);
-   }
-   else
-   {
-      TDFunction(transip, GetTime(), V);
-   }
-   if (Q)
-   {
-      V *= Q->Eval(T, ip, GetTime());
-   }
-}
-
-} // end namespace mfem
-

fem/complex_coefficient.hpp

@@ -1,523 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_COMPLEX_COEFFICIENT
-#define MFEM_COMPLEX_COEFFICIENT
-
-#include "../config/config.hpp"
-#include "../linalg/linalg.hpp"
-#include "coefficient.hpp"
-#include "intrules.hpp"
-#include "eltrans.hpp"
-
-namespace mfem
-{
-
-class ComplexCoefficient;
-class ComplexVectorCoefficient;
-class ComplexMatrixCoefficient;
-
-/// Standard Coefficient which returns the real part of a ComplexCoefficient
-class RealPartCoefficient : public Coefficient
-{
-private:
-   ComplexCoefficient &complex_coef_;
-
-public:
-   RealPartCoefficient(ComplexCoefficient & complex_coef)
-      : complex_coef_(complex_coef) {}
-
-   real_t Eval(ElementTransformation &T,
-               const IntegrationPoint &ip);
-};
-
-/// Standard Coefficient which returns the imaginary part of a
-/// ComplexCoefficient
-class ImagPartCoefficient : public Coefficient
-{
-private:
-   ComplexCoefficient &complex_coef_;
-
-public:
-   ImagPartCoefficient(ComplexCoefficient & complex_coef)
-      : complex_coef_(complex_coef) {}
-
-   real_t Eval(ElementTransformation &T,
-               const IntegrationPoint &ip);
-};
-
-typedef ImagPartCoefficient ImaginaryPartCoefficient;
-
-class RealPartVectorCoefficient : public VectorCoefficient
-{
-private:
-   ComplexVectorCoefficient &complex_vcoef_;
-   mutable ComplexVector val_;
-
-public:
-   RealPartVectorCoefficient(ComplexVectorCoefficient & complex_vcoef);
-
-   void Eval(Vector &V, ElementTransformation &T,
-             const IntegrationPoint &ip);
-};
-
-class ImagPartVectorCoefficient : public VectorCoefficient
-{
-private:
-   ComplexVectorCoefficient &complex_vcoef_;
-   mutable ComplexVector val_;
-
-public:
-   ImagPartVectorCoefficient(ComplexVectorCoefficient & complex_vcoef);
-
-   void Eval(Vector &V, ElementTransformation &T,
-             const IntegrationPoint &ip);
-};
-
-typedef ImagPartVectorCoefficient ImaginaryPartVectorCoefficient;
-
-class RealPartMatrixCoefficient : public MatrixCoefficient
-{
-private:
-   ComplexMatrixCoefficient &complex_mcoef_;
-   mutable ComplexTypeDenseMatrix val_;
-
-public:
-   RealPartMatrixCoefficient(ComplexMatrixCoefficient & complex_mcoef);
-
-   void Eval(DenseMatrix &M, ElementTransformation &T,
-             const IntegrationPoint &ip);
-};
-
-class ImagPartMatrixCoefficient : public MatrixCoefficient
-{
-private:
-   ComplexMatrixCoefficient &complex_mcoef_;
-   mutable ComplexTypeDenseMatrix val_;
-
-public:
-   ImagPartMatrixCoefficient(ComplexMatrixCoefficient & complex_mcoef);
-
-   void Eval(DenseMatrix &V, ElementTransformation &T,
-             const IntegrationPoint &ip);
-};
-
-typedef ImagPartMatrixCoefficient ImaginaryPartMatrixCoefficient;
-
-/** @brief Base class ComplexCoefficients that optionally depend on space and
-    time. These are used by the SesquilinearForm, ComplexLinearForm, and
-    ComplexGridFunction classes to represent the physical coefficients in
-    the PDEs that are being discretized. This class can also be used in a more
-    general way to represent functions that don't necessarily belong to a FE
-    space, e.g., to project onto ComplexGridFunctions to use as initial
-    conditions, exact solutions, etc. See, e.g., ex22 for these uses. */
-class ComplexCoefficient
-{
-protected:
-   real_t time;
-
-private:
-   RealPartCoefficient re_part_coef_;
-   ImagPartCoefficient im_part_coef_;
-
-protected:
-   Coefficient &real_coef_;
-   Coefficient &imag_coef_;
-
-public:
-
-   ComplexCoefficient();
-   ComplexCoefficient(Coefficient &c_r, Coefficient &c_i);
-
-   /// Set the time for time dependent coefficients
-   virtual void SetTime(real_t t)
-   { time = t; real_coef_.SetTime(t); imag_coef_.SetTime(t); }
-
-   /// Get the time for time dependent coefficients
-   real_t GetTime() { return time; }
-
-   /** @brief Evaluate the coefficient in the element described by @a T at the
-       point @a ip. */
-   /** @note When this method is called, the caller must make sure that the
-       IntegrationPoint associated with @a T is the same as @a ip. This can be
-       achieved by calling T.SetIntPoint(&ip). */
-   virtual complex_t Eval(ElementTransformation &T,
-                          const IntegrationPoint &ip);
-
-   /** @brief Evaluate the coefficient in the element described by @a T at the
-       point @a ip at time @a t. */
-   /** @note When this method is called, the caller must make sure that the
-       IntegrationPoint associated with @a T is the same as @a ip. This can be
-       achieved by calling T.SetIntPoint(&ip). */
-   complex_t Eval(ElementTransformation &T,
-                  const IntegrationPoint &ip, real_t t)
-   {
-      SetTime(t);
-      return Eval(T, ip);
-   }
-
-   /** @brief Access a standard Coefficient object reproducing the real part of
-       the complex-valued field */
-   /** @note By default this method returns an internal object which
-       computes the complex value using the above Eval method and
-       returns its real part. Custom implementations may choose to
-       override this method with a more efficient real-valued
-       coefficient. */
-   virtual Coefficient & real() { return real_coef_; }
-
-   /** @brief Access a standard Coefficient object reproducing the imaginary
-       part of the complex-valued field */
-   /** @note By default this method returns an internal object which
-       computes the complex value using the above Eval method and
-       returns its imaginary part. Custom implementations may choose to
-       override this method with a more efficient real-valued
-       coefficient. */
-   virtual Coefficient & imag() { return imag_coef_; }
-
-   virtual ~ComplexCoefficient() { }
-};
-
-/** @brief Base class ComplexVectorCoefficients that optionally depend
-    on space and time. These are used by the SesquilinearForm,
-    ComplexLinearForm, and ComplexGridFunction classes to represent
-    the physical vector-valued coefficients in the PDEs that are being
-    discretized. This class can also be used in a more general way to
-    represent functions that don't necessarily belong to a FE space,
-    e.g., to project onto ComplexGridFunctions to use as initial
-    conditions, exact solutions, etc. See, e.g., ex22 for these
-    uses. */
-class ComplexVectorCoefficient
-{
-protected:
-   int vdim;
-   real_t time;
-
-private:
-   RealPartVectorCoefficient re_part_vcoef_;
-   ImagPartVectorCoefficient im_part_vcoef_;
-
-protected:
-   VectorCoefficient &real_vcoef_;
-   VectorCoefficient &imag_vcoef_;
-
-   mutable Vector V_r_;
-   mutable Vector V_i_;
-
-public:
-   ComplexVectorCoefficient(int vd)
-      : vdim(vd), time(0.),
-        re_part_vcoef_(*this), im_part_vcoef_(*this),
-        real_vcoef_(re_part_vcoef_), imag_vcoef_(im_part_vcoef_)
-   { }
-
-   ComplexVectorCoefficient(VectorCoefficient &v_r, VectorCoefficient &v_i);
-
-
-   /// Set the time for time dependent coefficients
-   virtual void SetTime(real_t t)
-   { time = t; real_vcoef_.SetTime(t); imag_vcoef_.SetTime(t); }
-
-   /// Get the time for time dependent coefficients
-   real_t GetTime() { return time; }
-
-   /// Returns dimension of the vector.
-   int GetVDim() { return vdim; }
-
-   /** @brief Evaluate the vector coefficient in the element described by @a T
-       at the point @a ip, storing the result in @a V. */
-   /** @note When this method is called, the caller must make sure that the
-       IntegrationPoint associated with @a T is the same as @a ip. This can be
-       achieved by calling T.SetIntPoint(&ip). */
-   virtual void Eval(ComplexVector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip);
-
-   /** @brief Evaluate the vector coefficient in the element described by @a T
-       at the point @a ip at time @a t, storing the result in @a V. */
-   /** @note When this method is called, the caller must make sure that the
-       IntegrationPoint associated with @a T is the same as @a ip. This can be
-       achieved by calling T.SetIntPoint(&ip). */
-   void Eval(ComplexVector &V, ElementTransformation &T,
-             const IntegrationPoint &ip, real_t t)
-   {
-      SetTime(t);
-      Eval(V, T, ip);
-   }
-
-   /** @brief Access a standard Coefficient object reproducing the real part of
-       the complex-valued field */
-   /** @note By default this method returns an internal object which
-       computes the complex value using the above Eval method and
-       returns its real part. Custom implementations may choose to
-       override this method with a more efficient real-valued
-       coefficient. */
-   virtual VectorCoefficient & real() { return real_vcoef_; }
-
-   /** @brief Access a standard Coefficient object reproducing the imaginary
-       part of the complex-valued field */
-   /** @note By default this method returns an internal object which
-       computes the complex value using the above Eval method and
-       returns its imaginary part. Custom implementations may choose to
-       override this method with a more efficient real-valued
-       coefficient. */
-   virtual VectorCoefficient & imag() { return imag_vcoef_; }
-
-   virtual ~ComplexVectorCoefficient() { }
-};
-
-/** @brief Base class ComplexMatrixCoefficients that optionally depend
-    on space and time. These are used by the SesquilinearForm,
-    ComplexLinearForm, and ComplexGridFunction classes to represent
-    the physical matrix-valued coefficients in the PDEs that are being
-    discretized. This class can also be used in a more general way to
-    represent functions that don't necessarily belong to a FE space.
-    See, e.g., ex22 for these uses. */
-class ComplexMatrixCoefficient
-{
-protected:
-   int height, width;
-   real_t time;
-
-private:
-   RealPartMatrixCoefficient re_part_mcoef_;
-   ImagPartMatrixCoefficient im_part_mcoef_;
-
-protected:
-   MatrixCoefficient &real_mcoef_;
-   MatrixCoefficient &imag_mcoef_;
-
-   mutable DenseMatrix M_r_;
-   mutable DenseMatrix M_i_;
-
-public:
-   /// Construct a dim x dim matrix coefficient.
-   explicit ComplexMatrixCoefficient(int dim)
-      : height(dim), width(dim), time(0.),
-        re_part_mcoef_(*this), im_part_mcoef_(*this),
-        real_mcoef_(re_part_mcoef_), imag_mcoef_(im_part_mcoef_)
-   { }
-
-   /// Construct a h x w matrix coefficient.
-   ComplexMatrixCoefficient(int h, int w) :
-      height(h), width(w), time(0.),
-      re_part_mcoef_(*this), im_part_mcoef_(*this),
-      real_mcoef_(re_part_mcoef_), imag_mcoef_(im_part_mcoef_)
-   { }
-
-   /// Set the time for time dependent coefficients
-   virtual void SetTime(real_t t) { time = t; }
-
-   /// Get the time for time dependent coefficients
-   real_t GetTime() { return time; }
-
-   /// Get the height of the matrix.
-   int GetHeight() const { return height; }
-
-   /// Get the width of the matrix.
-   int GetWidth() const { return width; }
-
-   /// For backward compatibility get the width of the matrix.
-   int GetVDim() const { return width; }
-
-   /** @brief Evaluate the matrix coefficient in the element described by @a T
-       at the point @a ip, storing the result in @a K. */
-   /** @note When this method is called, the caller must make sure that the
-       IntegrationPoint associated with @a T is the same as @a ip. This can be
-       achieved by calling T.SetIntPoint(&ip). */
-   virtual void Eval(ComplexTypeDenseMatrix &K, ElementTransformation &T,
-                     const IntegrationPoint &ip) = 0;
-
-   /** @brief Access a standard Coefficient object reproducing the real part of
-       the complex-valued field */
-   /** @note By default this method returns an internal object which
-       computes the complex value using the above Eval method and
-       returns its real part. Custom implementations may choose to
-       override this method with a more efficient real-valued
-       coefficient. */
-   virtual MatrixCoefficient & real() { return real_mcoef_; }
-
-   /** @brief Access a standard Coefficient object reproducing the imaginary
-       part of the complex-valued field */
-   /** @note By default this method returns an internal object which
-       computes the complex value using the above Eval method and
-       returns its imaginary part. Custom implementations may choose to
-       override this method with a more efficient real-valued
-       coefficient. */
-   virtual MatrixCoefficient & imag() { return imag_mcoef_; }
-
-   virtual ~ComplexMatrixCoefficient() { }
-};
-
-/// A complex-valued coefficient that is constant across space and time
-class ComplexConstantCoefficient : public ComplexCoefficient
-{
-private:
-   complex_t val;
-
-   ConstantCoefficient real_coef;
-   ConstantCoefficient imag_coef;
-
-public:
-   ComplexConstantCoefficient(const complex_t z);
-
-   ComplexConstantCoefficient(real_t z_r, real_t z_i = 0.);
-
-   complex_t Eval(ElementTransformation &T,
-                  const IntegrationPoint &ip) { return val; }
-};
-
-/// Complex-valued vector coefficient that is constant in space and time.
-class ComplexVectorConstantCoefficient : public ComplexVectorCoefficient
-{
-private:
-   ComplexVector vec;
-
-public:
-   /// Construct the coefficient with constant vector @a v.
-   ComplexVectorConstantCoefficient(const ComplexVector &v)
-      : ComplexVectorCoefficient(v.Size()), vec(v) { }
-
-   /// Construct the coefficient with constant vector @a v.
-   ComplexVectorConstantCoefficient(const Vector &v)
-      : ComplexVectorCoefficient(v.Size()), vec(v) { }
-
-   using ComplexVectorCoefficient::Eval;
-
-   ///  Evaluate the vector coefficient at @a ip.
-   void Eval(ComplexVector &V, ElementTransformation &T,
-             const IntegrationPoint &ip) override { V = vec; }
-
-   /// Return a reference to the constant vector in this class.
-   const ComplexVector& GetVec() const { return vec; }
-};
-
-/// Complex-valued vector coefficient that is constant in space and time.
-class ComplexMatrixConstantCoefficient : public ComplexMatrixCoefficient
-{
-private:
-   ComplexTypeDenseMatrix mat;
-
-public:
-   /// Construct the coefficient with constant vector @a v.
-   ComplexMatrixConstantCoefficient(const ComplexTypeDenseMatrix &m)
-      : ComplexMatrixCoefficient(m.Height(), m.Width()), mat(m) { }
-
-   /// Construct the coefficient with constant vector @a v.
-   ComplexMatrixConstantCoefficient(const DenseMatrix &m)
-      : ComplexMatrixCoefficient(m.Height(), m.Width()), mat(m) { }
-
-   using ComplexMatrixCoefficient::Eval;
-
-   ///  Evaluate the matrix coefficient at @a ip.
-   void Eval(ComplexTypeDenseMatrix &M, ElementTransformation &T,
-             const IntegrationPoint &ip) override { M = mat; }
-
-   /// Return a reference to the constant matrix in this class.
-   const ComplexTypeDenseMatrix& GetMat() const { return mat; }
-};
-
-/// A general complex-valued function coefficient
-class ComplexFunctionCoefficient : public ComplexCoefficient
-{
-protected:
-   std::function<complex_t(const Vector &)> Function;
-   std::function<complex_t(const Vector &, real_t)> TDFunction;
-
-public:
-   /// Define a time-independent coefficient from a std function
-   /** \param F time-independent std::function */
-   ComplexFunctionCoefficient(std::function<complex_t
-                              (const Vector &)> F)
-      : Function(std::move(F))
-   { }
-
-   /// Define a time-dependent coefficient from a std function
-   /** \param TDF time-dependent function */
-   ComplexFunctionCoefficient(std::function<complex_t
-                              (const Vector &, real_t)> TDF)
-      : TDFunction(std::move(TDF))
-   { }
-
-   /// (DEPRECATED) Define a time-independent coefficient from a C-function
-   /** @deprecated Use the method where the C-function, @a f, uses a const
-       Vector argument instead of Vector. */
-   MFEM_DEPRECATED ComplexFunctionCoefficient(complex_t
-                                              (*f)(Vector &))
-   {
-      // Cast first to (void*) to suppress a warning from newer version of
-      // Clang when using -Wextra.
-      Function = reinterpret_cast<complex_t(*)
-                 (const Vector&)>((void*)f);
-      TDFunction = NULL;
-   }
-
-   /// (DEPRECATED) Define a time-dependent coefficient from a C-function
-   /** @deprecated Use the method where the C-function, @a tdf, uses a const
-       Vector argument instead of Vector. */
-   MFEM_DEPRECATED ComplexFunctionCoefficient(complex_t
-                                              (*tdf)(Vector &, real_t))
-   {
-      Function = NULL;
-      // Cast first to (void*) to suppress a warning from newer version of
-      // Clang when using -Wextra.
-      TDFunction =
-         reinterpret_cast<complex_t(*)(const Vector&,
-                                       real_t)>((void*)tdf);
-   }
-
-   /// Evaluate the coefficient at @a ip.
-   complex_t Eval(ElementTransformation &T,
-                  const IntegrationPoint &ip) override;
-};
-
-/// A general vector function coefficient
-class ComplexVectorFunctionCoefficient : public ComplexVectorCoefficient
-{
-private:
-   std::function<void(const Vector &, ComplexVector &)> Function;
-   std::function<void(const Vector &, real_t, ComplexVector &)> TDFunction;
-   ComplexCoefficient *Q;
-
-public:
-   /// Define a time-independent complex-valued vector coefficient
-   /// from a std function
-   /** \param dim - the size of the vector
-       \param F - time-independent function
-       \param q - optional scalar Coefficient to scale the vector coefficient */
-   ComplexVectorFunctionCoefficient(int dim,
-                                    std::function<void(const Vector &,
-                                                       ComplexVector &)> F,
-                                    ComplexCoefficient *q = nullptr)
-      : ComplexVectorCoefficient(dim), Function(std::move(F)), Q(q)
-   { }
-
-   /// Define a time-dependent complex-valued vector coefficient from
-   /// a std function
-   /** \param dim - the size of the vector
-       \param TDF - time-dependent function
-       \param q - optional scalar ComplexCoefficient to scale the vector coefficient */
-   ComplexVectorFunctionCoefficient(int dim,
-                                    std::function<void(const Vector &, real_t,
-                                                       ComplexVector &)> TDF,
-                                    ComplexCoefficient *q = nullptr)
-      : ComplexVectorCoefficient(dim), TDFunction(std::move(TDF)), Q(q)
-   { }
-
-   using ComplexVectorCoefficient::Eval;
-   /// Evaluate the vector coefficient at @a ip.
-   void Eval(ComplexVector &V, ElementTransformation &T,
-             const IntegrationPoint &ip) override;
-
-   virtual ~ComplexVectorFunctionCoefficient() { }
-};
-
-} // end namespace mfem
-
-#endif

fem/complex_fem.cpp

@@ -96,23 +96,6 @@ ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
    gfi->SyncAliasMemory(*this);
 }
 
-void
-ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff)
-{
-   gfr->SyncMemory(*this);
-   gfi->SyncMemory(*this);
-   gfr->ProjectCoefficient(real_coeff);
-   *gfi = 0.0;
-   gfr->SyncAliasMemory(*this);
-   gfi->SyncAliasMemory(*this);
-}
-
-void
-ComplexGridFunction::ProjectCoefficient(ComplexCoefficient &coeff)
-{
-   this->ProjectCoefficient(coeff.real(), coeff.imag());
-}
-
 void
 ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                         VectorCoefficient &imag_vcoeff)
@@ -125,23 +108,6 @@ ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
    gfi->SyncAliasMemory(*this);
 }
 
-void
-ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff)
-{
-   gfr->SyncMemory(*this);
-   gfi->SyncMemory(*this);
-   gfr->ProjectCoefficient(real_vcoeff);
-   *gfi = 0.0;
-   gfr->SyncAliasMemory(*this);
-   gfi->SyncAliasMemory(*this);
-}
-
-void
-ComplexGridFunction::ProjectCoefficient(ComplexVectorCoefficient &vcoeff)
-{
-   this->ProjectCoefficient(vcoeff.real(), vcoeff.imag());
-}
-
 void
 ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                            Coefficient &imag_coeff,
@@ -155,26 +121,6 @@ ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
    gfi->SyncAliasMemory(*this);
 }
 
-void
-ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
-                                           Array<int> &attr)
-{
-   ConstantCoefficient zero_coeff(0.0);
-   gfr->SyncMemory(*this);
-   gfi->SyncMemory(*this);
-   gfr->ProjectBdrCoefficient(real_coeff, attr);
-   gfi->ProjectBdrCoefficient(zero_coeff, attr);
-   gfr->SyncAliasMemory(*this);
-   gfi->SyncAliasMemory(*this);
-}
-
-void
-ComplexGridFunction::ProjectBdrCoefficient(ComplexCoefficient &coeff,
-                                           Array<int> &attr)
-{
-   this->ProjectBdrCoefficient(coeff.real(), coeff.imag(), attr);
-}
-
 void
 ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
                                                  VectorCoefficient &imag_vcoeff,
@@ -188,28 +134,6 @@ ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
    gfi->SyncAliasMemory(*this);
 }
 
-void
-ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
-                                                 Array<int> &attr)
-{
-   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
-   VectorConstantCoefficient zero_vcoeff(zero_vec);
-   gfr->SyncMemory(*this);
-   gfi->SyncMemory(*this);
-   gfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
-   gfi->ProjectBdrCoefficientNormal(zero_vcoeff, attr);
-   gfr->SyncAliasMemory(*this);
-   gfi->SyncAliasMemory(*this);
-}
-
-void
-ComplexGridFunction::ProjectBdrCoefficientNormal(
-   ComplexVectorCoefficient &vcoeff,
-   Array<int> &attr)
-{
-   this->ProjectBdrCoefficientNormal(vcoeff.real(), vcoeff.imag(), attr);
-}
-
 void
 ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                   &real_vcoeff,
@@ -225,80 +149,6 @@ ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
    gfi->SyncAliasMemory(*this);
 }
 
-void
-ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
-                                                  &real_vcoeff,
-                                                  Array<int> &attr)
-{
-   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
-   VectorConstantCoefficient zero_vcoeff(zero_vec);
-   gfr->SyncMemory(*this);
-   gfi->SyncMemory(*this);
-   gfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
-   gfi->ProjectBdrCoefficientTangent(zero_vcoeff, attr);
-   gfr->SyncAliasMemory(*this);
-   gfi->SyncAliasMemory(*this);
-}
-
-void
-ComplexGridFunction::ProjectBdrCoefficientTangent(
-   ComplexVectorCoefficient &vcoeff,
-   Array<int> &attr)
-{
-   this->ProjectBdrCoefficientTangent(vcoeff.real(), vcoeff.imag(), attr);
-}
-
-real_t
-ComplexGridFunction::ComputeL2Error(Coefficient &re_exsol,
-                                    Coefficient &im_exsol,
-                                    const IntegrationRule *irs[],
-                                    const Array<int> *elems) const
-{
-   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
-   real_t err_i = gfi->ComputeL2Error(im_exsol, irs, elems);
-
-   return sqrt(err_r * err_r + err_i * err_i);
-}
-
-real_t
-ComplexGridFunction::ComputeL2Error(Coefficient &re_exsol,
-                                    const IntegrationRule *irs[],
-                                    const Array<int> *elems) const
-{
-   ConstantCoefficient zero_coef(0.0);
-
-   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
-   real_t err_i = gfi->ComputeL2Error(zero_coef, irs, elems);
-
-   return sqrt(err_r * err_r + err_i * err_i);
-}
-
-real_t
-ComplexGridFunction::ComputeL2Error(VectorCoefficient &re_exsol,
-                                    VectorCoefficient &im_exsol,
-                                    const IntegrationRule *irs[],
-                                    const Array<int> *elems) const
-{
-   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
-   real_t err_i = gfi->ComputeL2Error(im_exsol, irs, elems);
-
-   return sqrt(err_r * err_r + err_i * err_i);
-}
-
-real_t
-ComplexGridFunction::ComputeL2Error(VectorCoefficient &re_exsol,
-                                    const IntegrationRule *irs[],
-                                    const Array<int> *elems) const
-{
-   Vector zero_vec(re_exsol.GetVDim()); zero_vec = 0.0;
-   VectorConstantCoefficient zero_coef(zero_vec);
-
-   real_t err_r = gfr->ComputeL2Error(re_exsol, irs, elems);
-   real_t err_i = gfi->ComputeL2Error(zero_coef, irs, elems);
-
-   return sqrt(err_r * err_r + err_i * err_i);
-}
-
 
 ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
                                      ComplexOperator::Convention convention)
@@ -881,17 +731,6 @@ ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
    pgfi->SyncAliasMemory(*this);
 }
 
-void
-ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff)
-{
-   pgfr->SyncMemory(*this);
-   pgfi->SyncMemory(*this);
-   pgfr->ProjectCoefficient(real_coeff);
-   *pgfi = 0.0;
-   pgfr->SyncAliasMemory(*this);
-   pgfi->SyncAliasMemory(*this);
-}
-
 void
 ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                            VectorCoefficient &imag_vcoeff)
@@ -904,17 +743,6 @@ ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
    pgfi->SyncAliasMemory(*this);
 }
 
-void
-ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff)
-{
-   pgfr->SyncMemory(*this);
-   pgfi->SyncMemory(*this);
-   pgfr->ProjectCoefficient(real_vcoeff);
-   *pgfi = 0.0;
-   pgfr->SyncAliasMemory(*this);
-   pgfi->SyncAliasMemory(*this);
-}
-
 void
 ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                               Coefficient &imag_coeff,
@@ -928,19 +756,6 @@ ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
    pgfi->SyncAliasMemory(*this);
 }
 
-void
-ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
-                                              Array<int> &attr)
-{
-   ConstantCoefficient zero_coeff(0.0);
-   pgfr->SyncMemory(*this);
-   pgfi->SyncMemory(*this);
-   pgfr->ProjectBdrCoefficient(real_coeff, attr);
-   pgfi->ProjectBdrCoefficient(zero_coeff, attr);
-   pgfr->SyncAliasMemory(*this);
-   pgfi->SyncAliasMemory(*this);
-}
-
 void
 ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
                                                     &real_vcoeff,
@@ -956,21 +771,6 @@ ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
    pgfi->SyncAliasMemory(*this);
 }
 
-void
-ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
-                                                    &real_vcoeff,
-                                                    Array<int> &attr)
-{
-   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
-   VectorConstantCoefficient zero_vcoeff(zero_vec);
-   pgfr->SyncMemory(*this);
-   pgfi->SyncMemory(*this);
-   pgfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
-   pgfi->ProjectBdrCoefficientNormal(zero_vcoeff, attr);
-   pgfr->SyncAliasMemory(*this);
-   pgfi->SyncAliasMemory(*this);
-}
-
 void
 ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                      &real_vcoeff,
@@ -986,21 +786,6 @@ ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
    pgfi->SyncAliasMemory(*this);
 }
 
-void
-ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
-                                                     &real_vcoeff,
-                                                     Array<int> &attr)
-{
-   Vector zero_vec(real_vcoeff.GetVDim()); zero_vec = 0.;
-   VectorConstantCoefficient zero_vcoeff(zero_vec);
-   pgfr->SyncMemory(*this);
-   pgfi->SyncMemory(*this);
-   pgfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
-   pgfi->ProjectBdrCoefficientTangent(zero_vcoeff, attr);
-   pgfr->SyncAliasMemory(*this);
-   pgfi->SyncAliasMemory(*this);
-}
-
 void
 ParComplexGridFunction::Distribute(const Vector *tv)
 {
@@ -1040,31 +825,6 @@ ParComplexGridFunction::ParallelProject(Vector &tv) const
    tvi.SyncAliasMemory(tv);
 }
 
-real_t
-ParComplexGridFunction::ComputeL2Error(Coefficient &exsolr,
-                                       const IntegrationRule *irs[],
-                                       Array<int> *elems) const
-{
-   ConstantCoefficient zeroCoef(0.0);
-
-   real_t err_r = pgfr->ComputeL2Error(exsolr, irs, elems);
-   real_t err_i = pgfi->ComputeL2Error(zeroCoef, irs, elems);
-   return sqrt(err_r * err_r + err_i * err_i);
-}
-
-real_t
-ParComplexGridFunction::ComputeL2Error(VectorCoefficient &exsolr,
-                                       const IntegrationRule *irs[],
-                                       Array<int> *elems) const
-{
-   Vector zeroVec(exsolr.GetVDim()); zeroVec = 0.0;
-   VectorConstantCoefficient zeroCoef(zeroVec);
-
-   real_t err_r = pgfr->ComputeL2Error(exsolr, irs, elems);
-   real_t err_i = pgfi->ComputeL2Error(zeroCoef, irs, elems);
-   return sqrt(err_r * err_r + err_i * err_i);
-}
-
 
 ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
                                            ComplexOperator::Convention

fem/conduitdatacollection.cpp

@@ -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;
       }
    }
 

fem/derefmat_op.cpp

@@ -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

fem/derefmat_op.hpp

@@ -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

fem/dfem/util.hpp

@@ -568,7 +568,7 @@ struct ThreadBlocks
    int z = 1;
 };
 
-#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA_OR_HIP)
 template <typename func_t>
 __global__ void forall_kernel_shmem(func_t f, int n)
 {
@@ -591,7 +591,7 @@ void forall(func_t f,
    if (Device::Allows(Backend::CUDA_MASK) ||
        Device::Allows(Backend::HIP_MASK))
    {
-#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA_OR_HIP)
       // int gridsize = (N + Z - 1) / Z;
       int num_bytes = num_shmem * sizeof(decltype(shmem));
       dim3 block_size(blocks.x, blocks.y, blocks.z);

fem/fes_kernels.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

fem/fespace.cpp

@@ -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)

fem/fespace.hpp

@@ -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).

fem/intrules_cut.hpp

@@ -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));
 

fem/pderefmat_op.cpp

@@ -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

fem/pderefmat_op.hpp

@@ -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

fem/pfespace.cpp

@@ -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);

fem/pfespace.hpp

@@ -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;

general/CMakeLists.txt

@@ -39,7 +39,6 @@ list(APPEND HDRS
   arrays_by_name.hpp
   backends.hpp
   binaryio.hpp
-  complex_type.hpp
   cuda.hpp
   device.hpp
   error.hpp

general/array.hpp

@@ -190,6 +190,12 @@ public:
    /// Prepend an 'el' to the array, resize if necessary.
    inline int Prepend(const T &el);
 
+   /// Insert @a els into the array at index @a i
+   inline int Insert(int i, const Array<T> &els);
+
+   /// Insert @a el into the array at index @a i
+   inline int Insert(int i, const T &el) { return Insert(i, Array<T>({el})); }
+
    /// Return the last element in the array.
    inline T &Last();
 
@@ -211,6 +217,9 @@ public:
    /// Delete the first entry with value == 'el'.
    inline void DeleteFirst(const T &el);
 
+   /// Delete entries at @a indices, and resize
+   inline void DeleteAt(const Array<int> &indices);
+
    /// Delete the whole array.
    inline void DeleteAll();
 
@@ -249,6 +258,9 @@ public:
    /// Copy sub array starting from @a offset out to the provided @a sa.
    inline void GetSubArray(int offset, int sa_size, Array<T> &sa) const;
 
+   /// Set from sub array @sa at @a offset
+   inline void SetSubArray(int offset, const Array<T> &sa);
+
    /// Prints array to stream with width elements per row.
    void Print(std::ostream &out = mfem::out, int width = 4) const;
 
@@ -873,6 +885,22 @@ inline int Array<T>::Prepend(const T &el)
    return size;
 }
 
+template<class T>
+inline int Array<T>::Insert(int i, const Array<T> &els)
+{
+   MFEM_ASSERT(i < size, "Insert index is out-of-bounds.");
+
+   const int old_size = size;
+
+   SetSize(size + els.Size());
+   for (int j = old_size-1; j >= i; j--)
+   {
+      data[j+els.Size()] = data[j];
+   }
+   SetSubArray(i, els);
+   return size;
+}
+
 template <class T>
 inline T &Array<T>::Last()
 {
@@ -935,6 +963,30 @@ inline void Array<T>::DeleteFirst(const T &el)
    }
 }
 
+template <class T>
+inline void Array<T>::DeleteAt(const Array<int> &indices)
+{
+   // Make a copy of the indices, sorted.
+   Array<int> sorted_indices(indices);
+   sorted_indices.Sort();
+
+   int rm_count = 0;
+   for (int i = 0; i < size; i++)
+   {
+      if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
+      {
+         rm_count++;
+      }
+      else
+      {
+         data[i-rm_count] = data[i]; // shift data rm_count
+      }
+   }
+
+   // Resize to remove tail
+   SetSize(size - rm_count);
+}
+
 template <class T>
 inline void Array<T>::DeleteAll()
 {
@@ -987,6 +1039,18 @@ inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
    }
 }
 
+template<class T>
+inline void Array<T>::SetSubArray(int offset, const Array<T> &sa)
+{
+   MFEM_ASSERT(offset + sa.Size() < size,
+               "Sub-array with size " << sa.Size() << " is too large to set at offset " <<
+               offset << ", given array size " << size);
+   for (int i = 0; i < sa.Size(); i++)
+   {
+      data[offset + i] = sa[i];
+   }
+}
+
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {

general/backends.hpp

@@ -14,7 +14,7 @@
 
 #include "../config/config.hpp"
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #include <cusparse.h>
 #include <library_types.h>
 #include <cuda_runtime.h>
@@ -22,7 +22,7 @@
 #endif
 #include "cuda.hpp"
 
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 #include <hip/hip_runtime.h>
 #endif
 #include "hip.hpp"
@@ -43,7 +43,7 @@
 #endif
 #endif
 
-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if !defined(MFEM_USE_CUDA_OR_HIP)
 #define MFEM_DEVICE
 #define MFEM_HOST
 #define MFEM_LAMBDA
@@ -55,7 +55,7 @@
 #endif
 
 #if !((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
-      (defined(MFEM_USE_HIP)  && defined(__HIP_DEVICE_COMPILE__)))
+      (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
 #define MFEM_SHARED
 #define MFEM_SYNC_THREAD
 #define MFEM_BLOCK_ID(k) 0
@@ -66,7 +66,7 @@
 #endif
 
 // 'double' and 'float' atomicAdd implementation for previous versions of CUDA
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && (__CUDA_ARCH__ < 600)
 MFEM_DEVICE inline mfem::real_t atomicAdd(mfem::real_t *add, mfem::real_t val)
 {
    unsigned long long int *ptr = (unsigned long long int *) add;
@@ -94,7 +94,7 @@ template <typename T>
 MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
 {
 #if ((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
-     (defined(MFEM_USE_HIP)  && defined(__HIP_DEVICE_COMPILE__)))
+     (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
    return atomicAdd(&add,val);
 #else
    T old = add;

general/communication.hpp

@@ -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.

general/complex_type.hpp

@@ -1,125 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_COMPLEX_TYPE
-#define MFEM_COMPLEX_TYPE
-
-#include "../config/config.hpp"
-
-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
-#include <complex>
-#include <utility>
-#endif
-
-#if defined(MFEM_USE_CUDA)
-#include <cuComplex.h>
-#endif
-
-#if defined(MFEM_USE_HIP)
-#include <hip/hip_complex.h>
-#endif
-
-
-namespace mfem
-{
-
-/// @brief Complex number type for device.
-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
-
-#define zAbs std::abs
-#define zExp std::exp
-#define zNorm std::norm
-using complex_t = std::complex<real_t>;
-
-#else // CUDA or HIP
-
-#if defined(MFEM_USE_CUDA)
-using DoubleComplex_t = cuDoubleComplex;
-#endif
-
-#if defined(MFEM_USE_HIP)
-using DoubleComplex_t = hipDoubleComplex;
-#endif
-
-struct Complex : public DoubleComplex_t
-{
-   MFEM_HOST_DEVICE Complex() = default;
-   MFEM_HOST_DEVICE Complex(real_t r) { x = r, y = 0.0; }
-   MFEM_HOST_DEVICE Complex(real_t r, real_t i) { x = r, y = i; }
-   MFEM_HOST_DEVICE real_t real() const { return x; }
-   MFEM_HOST_DEVICE void real(real_t r) { x = r; }
-   MFEM_HOST_DEVICE real_t imag() const { return y; }
-   MFEM_HOST_DEVICE void imag(real_t i) { y = i; }
-
-   template <typename U>
-   MFEM_HOST_DEVICE inline Complex &operator*=(const U &z)
-   {
-      return *this = *this * z, *this;
-   }
-
-   template <typename U>
-   MFEM_HOST_DEVICE inline Complex &operator/=(const U &z)
-   {
-      return *this = *this / z, *this;
-   }
-};
-
-MFEM_HOST_DEVICE inline Complex operator*(const Complex &x, const real_t &y)
-{
-   return Complex(x.real() * y, x.imag() * y);
-}
-
-MFEM_HOST_DEVICE inline Complex operator+(const Complex &a, const Complex &b)
-{
-   return Complex(a.real() + b.real(), a.imag() + b.imag());
-}
-
-MFEM_HOST_DEVICE inline Complex operator*(const real_t d, const Complex &z)
-{
-   return Complex(z.real() * d, z.imag() * d);
-}
-
-MFEM_HOST_DEVICE inline Complex operator*(const Complex &a, const Complex &b)
-{
-   return Complex(a.real() * b.real() - a.imag() * b.imag(),
-                  a.real() * b.imag() + a.imag() * b.real());
-}
-
-MFEM_HOST_DEVICE inline Complex operator/(const Complex &z, const real_t &d)
-{
-   return Complex(z.real() / d, z.imag() / d);
-}
-
-MFEM_HOST_DEVICE inline real_t zAbs(const Complex &z)
-{
-   return std::hypot(z.real(), z.imag());
-}
-
-MFEM_HOST_DEVICE inline Complex zExp(const Complex &q)
-{
-   Complex z;
-   real_t s, c, e = std::exp(q.real());
-   sincos(q.imag(), &s, &c);
-   z.real(c * e), z.imag(s * e);
-   return z;
-}
-
-MFEM_HOST_DEVICE inline real_t zNorm(const Complex &z)
-{
-   return z.real() * z.real() + z.imag() * z.imag();
-}
-
-using complex_t = Complex;
-#endif // MFEM_USE_CUDA || MFEM_USE_HIP
-
-} // namespace mfem
-
-#endif // MFEM_COMPLEX_TYPE

general/cuda.cpp

@@ -24,7 +24,7 @@ void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                      const char *file, int line)
 {
    mfem::err << "\n\nCUDA error: (" << expr << ") failed with error:\n --> "
-             << cudaGetErrorString(err)
+             << cudaGetErrorString(err) << " [code: " << (int)err << ']'
              << "\n ... in function: " << func
              << "\n ... in file: " << file << ':' << line << '\n';
    mfem_error();

general/cuda.hpp

@@ -18,7 +18,7 @@
 // CUDA block size used by MFEM.
 #define MFEM_CUDA_BLOCKS 256
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,10 +37,9 @@
                               __FILE__, __LINE__);                             \
     }                                                                          \
   } while (0)
-#endif // MFEM_USE_CUDA
 
 // Define the MFEM inner threading macros
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
+#if defined(__CUDA_ARCH__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) blockIdx.k
@@ -48,12 +47,13 @@
 #define MFEM_THREAD_SIZE(k) blockDim.k
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=threadIdx.k; i<N; i+=blockDim.k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) if(const int i=threadIdx.k; i<N)
-#endif
+#endif // defined(__CUDA_ARCH__)
+#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 
 namespace mfem
 {
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 // Function used by the macro MFEM_GPU_CHECK.
 void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                      const char *file, int line);

general/device.cpp

@@ -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;

general/device.hpp

@@ -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. */

general/error.hpp

@@ -176,7 +176,7 @@ __device__ void abort_msg(T & msg)
       printf(__VA_ARGS__);     \
       asm("trap;");            \
    }
-#elif defined(MFEM_USE_HIP)
+#elif defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_ABORT_KERNEL(...) \
    {                           \
       printf(__VA_ARGS__);     \

general/forall.hpp

@@ -158,8 +158,8 @@ private:
 #define MFEM_PRAGMA(X) _Pragma(#X)
 
 // MFEM_UNROLL pragma macro that can be used inside MFEM_FORALL macros.
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
-#ifdef __NVCC__
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) // Clang cuda or nvcc
+#ifdef __NVCC__ // nvcc specifically
 #define MFEM_UNROLL(N) MFEM_PRAGMA(unroll(N))
 #else // Assuming Clang CUDA
 #define MFEM_UNROLL(N) MFEM_PRAGMA(unroll N)
@@ -169,12 +169,12 @@ private:
 #endif
 
 // MFEM_GPU_FORALL: "parallel for" executed with CUDA or HIP based on the MFEM
-// build-time configuration (MFEM_USE_CUDA or MFEM_USE_HIP). If neither CUDA nor
-// HIP is enabled, this macro is a no-op.
-#if defined(MFEM_USE_CUDA)
+// build-time configuration (MFEM_USE_CUDA or MFEM_USE_HIP), and if compiling
+// with CUDA/HIP language. Otherwise, this macro is a no-op.
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_GPU_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
                                        (int i) {__VA_ARGS__})
-#elif defined(MFEM_USE_HIP)
+#elif defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_GPU_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
                                         (int i) {__VA_ARGS__})
 #else
@@ -481,7 +481,7 @@ void RajaSeqWrap(const int N, HBODY &&h_body)
 
 
 /// CUDA backend
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 
 template <typename BODY> __global__ static
 void CuKernel1D(const int N, BODY body)
@@ -573,11 +573,11 @@ struct CuWrap<3>
    }
 };
 
-#endif // MFEM_USE_CUDA
+#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 
 
 /// HIP backend
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 
 template <typename BODY> __global__ static
 void HipKernel1D(const int N, BODY body)
@@ -668,7 +668,7 @@ struct HipWrap<3>
    }
 };
 
-#endif // MFEM_USE_HIP
+#endif // defined(MFEM_USE_HIP) && defined(__HIP__)
 
 
 /// The forall kernel body wrapper
@@ -701,7 +701,7 @@ inline void ForallWrap(const bool use_dev, const int N,
    }
 #endif
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
    // If Backend::CUDA is allowed, use it
    if (Device::Allows(Backend::CUDA))
    {
@@ -709,7 +709,7 @@ inline void ForallWrap(const bool use_dev, const int N,
    }
 #endif
 
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
    // If Backend::HIP is allowed, use it
    if (Device::Allows(Backend::HIP))
    {

general/hip.cpp

@@ -24,7 +24,7 @@ void mfem_hip_error(hipError_t err, const char *expr, const char *func,
                     const char *file, int line)
 {
    mfem::err << "\n\nHIP error: (" << expr << ") failed with error:\n --> "
-             << hipGetErrorString(err)
+             << hipGetErrorString(err) << " [code: " << (int)err << ']'
              << "\n ... in function: " << func
              << "\n ... in file: " << file << ':' << line << '\n';
    mfem_error();

general/hip.hpp

@@ -18,7 +18,7 @@
 // HIP block size used by MFEM.
 #define MFEM_HIP_BLOCKS 256
 
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,10 +37,9 @@
                              __FILE__, __LINE__);                              \
     }                                                                          \
   } while (0)
-#endif // MFEM_USE_HIP
 
 // Define the MFEM inner threading macros
-#if defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)
+#if defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) hipBlockIdx_ ##k
@@ -50,7 +49,8 @@
    for(int i=hipThreadIdx_ ##k; i<N; i+=hipBlockDim_ ##k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) \
    if(const int i=hipThreadIdx_ ##k; i<N)
-#endif
+#endif // defined(__HIP_DEVICE_COMPILE__)
+#endif // defined(MFEM_USE_HIP) && defined(__HIP__)
 
 namespace mfem
 {

general/mem_manager.cpp

@@ -513,10 +513,7 @@ public:
    void *HtoD(void *dst, const void *src, size_t bytes) override
    { return HipMemcpyHtoD(dst, src, bytes); }
    void *DtoD(void* dst, const void* src, size_t bytes) override
-   // Unlike cudaMemcpy(DtoD), hipMemcpy(DtoD) causes a host-side synchronization so
-   // instead we use hipMemcpyAsync to get similar behavior.
-   // for more info see: https://github.com/mfem/mfem/pull/2780
-   { return HipMemcpyDtoDAsync(dst, src, bytes); }
+   { return HipMemcpyDtoD(dst, src, bytes); }
    void *DtoH(void *dst, const void *src, size_t bytes) override
    { return HipMemcpyDtoH(dst, src, bytes); }
 };
@@ -658,10 +655,7 @@ public:
       return CuMemcpyDtoD(dst, src, bytes);
 #endif
 #ifdef MFEM_USE_HIP
-      // Unlike cudaMemcpy(DtoD), hipMemcpy(DtoD) causes a host-side synchronization so
-      // instead we use hipMemcpyAsync to get similar behavior.
-      // for more info see: https://github.com/mfem/mfem/pull/2780
-      return HipMemcpyDtoDAsync(dst, src, bytes);
+      return HipMemcpyDtoD(dst, src, bytes);
 #endif
       // rm.copy(dst, const_cast<void*>(src), bytes); return dst;
    }

general/mem_manager.hpp

@@ -896,6 +896,7 @@ inline HYPRE_MemoryLocation GetHypreMemoryLocation()
 #elif MFEM_HYPRE_VERSION < 23100
    return HYPRE_MEMORY_DEVICE;
 #else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
+   if (!HYPRE_Initialized()) { return HYPRE_MEMORY_HOST; }
    HYPRE_MemoryLocation loc;
    HYPRE_GetMemoryLocation(&loc);
    return loc;
@@ -1057,7 +1058,8 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
          // register the 'base' if the MemoryManager::Exists():
          MemoryManager::Exists()
 #else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
-         MemoryManager::Exists() && HypreUsingGPU()
+         IsDeviceMemory(MemoryManager::GetDeviceMemoryType()) ||
+         (MemoryManager::Exists() && HypreUsingGPU())
 #endif
       )
       {

general/reducers.hpp

@@ -537,7 +537,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
       return;
    }
 
-#if defined(MFEM_USE_HIP) || defined(MFEM_USE_CUDA)
+#if defined(MFEM_USE_CUDA_OR_HIP)
    if (use_dev &&
        mfem::Device::Allows(Backend::CUDA | Backend::HIP | Backend::RAJA_CUDA |
                             Backend::RAJA_HIP))

linalg/CMakeLists.txt

@@ -21,7 +21,6 @@ list(APPEND SRCS
   blockvector.cpp
   complex_densemat.cpp
   complex_operator.cpp
-  complex_vector.cpp
   constraints.cpp
   densemat.cpp
   symmat.cpp
@@ -48,7 +47,6 @@ list(APPEND HDRS
   blockvector.hpp
   complex_densemat.hpp
   complex_operator.hpp
-  complex_vector.hpp
   constraints.hpp
   densemat.hpp
   dinvariants.hpp

linalg/complex_densemat.cpp

@@ -9,7 +9,6 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "../general/forall.hpp"
 #include "complex_densemat.hpp"
 #include "lapack.hpp"
 #include <complex>
@@ -17,8 +16,6 @@
 namespace mfem
 {
 
-using namespace std;
-
 DenseMatrix & ComplexDenseMatrix::real()
 {
    MFEM_ASSERT(Op_Real_, "ComplexDenseMatrix has no real part!");
@@ -1020,303 +1017,4 @@ void ComplexCholeskyFactors::GetInverseMatrix(int m, real_t * X_r,
    delete [] X;
 }
 
-ComplexTypeDenseMatrix::ComplexTypeDenseMatrix()
-   : height(0), width(0)
-{}
-
-ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(const ComplexTypeDenseMatrix &m)
-   : height(m.Height()), width(m.Width())
-{
-   const int hw = height * width;
-   if (hw > 0)
-   {
-      MFEM_ASSERT(m.data, "invalid source matrix");
-      data.New(hw);
-      std::memcpy(data, m.data, sizeof(complex_t)*hw);
-   }
-}
-
-ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(const DenseMatrix &m)
-   : height(m.Height()), width(m.Width())
-{
-   const int hw = height * width;
-   if (hw > 0)
-   {
-      MFEM_ASSERT(m.data, "invalid source matrix");
-      data.New(hw);
-
-      for (int i = 0; i < hw; i++)
-      {
-         data[i] = m.data[i];
-      }
-   }
-}
-
-ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(int s)
-   : height(s), width(s)
-{
-   MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
-   if (s > 0)
-   {
-      data.New(s*s);
-      *this = 0.0; // init with zeroes
-   }
-}
-
-ComplexTypeDenseMatrix::ComplexTypeDenseMatrix(int m, int n)
-   : height(m), width(n)
-{
-   MFEM_ASSERT(m >= 0 && n >= 0,
-               "invalid DenseMatrix size: " << m << " x " << n);
-   const int capacity = m*n;
-   if (capacity > 0)
-   {
-      data.New(capacity);
-      *this = 0.0; // init with zeroes
-   }
-}
-
-void ComplexTypeDenseMatrix::SetSize(int h, int w)
-{
-   MFEM_ASSERT(h >= 0 && w >= 0,
-               "invalid ComplexTypeDenseMatrix size: " << h << " x " << w);
-   if (Height() == h && Width() == w)
-   {
-      return;
-   }
-   height = h;
-   width = w;
-   const int hw = h*w;
-   if (hw > data.Capacity())
-   {
-      data.Delete();
-      data.New(hw);
-      *this = 0.0; // init with zeroes
-   }
-}
-
-/// Returns reference to a_{ij}.
-complex_t &ComplexTypeDenseMatrix::Elem(int i, int j)
-{
-   return (*this)(i,j);
-}
-
-/// Returns constant reference to a_{ij}.
-const complex_t &ComplexTypeDenseMatrix::Elem(int i, int j) const
-{
-   return (*this)(i,j);
-}
-
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(real_t c)
-{
-   const int s = Height()*Width();
-   for (int i = 0; i < s; i++)
-   {
-      data[i] = c;
-   }
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(complex_t c)
-{
-   const int s = Height()*Width();
-   for (int i = 0; i < s; i++)
-   {
-      data[i] = c;
-   }
-   return *this;
-}
-
-/// Copy the matrix entries from the given array
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(const real_t *d)
-{
-   const int s = Height()*Width();
-   for (int i = 0; i < s; i++)
-   {
-      data[i] = d[i];
-   }
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=
-(const complex_t *d)
-{
-   const int s = Height()*Width();
-   for (int i = 0; i < s; i++)
-   {
-      data[i] = d[i];
-   }
-   return *this;
-}
-
-/// Sets the matrix size and elements equal to those of m
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=(const DenseMatrix &m)
-{
-   SetSize(m.height, m.width);
-
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] = m.data[i];
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator=
-(const ComplexTypeDenseMatrix &m)
-{
-   SetSize(m.height, m.width);
-
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] = m.data[i];
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=(const real_t *m)
-{
-   const int s = Height()*Width();
-   for (int i = 0; i < s; i++)
-   {
-      data[i] += m[i];
-   }
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=
-(const complex_t *m)
-{
-   const int s = Height()*Width();
-   for (int i = 0; i < s; i++)
-   {
-      data[i] += m[i];
-   }
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=(const DenseMatrix &m)
-{
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] += m.data[i];
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator+=
-(const ComplexTypeDenseMatrix &m)
-{
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] += m.data[i];
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator-=(const DenseMatrix &m)
-{
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] -= m.data[i];
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator-=
-(const ComplexTypeDenseMatrix &m)
-{
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] -= m.data[i];
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator*=(real_t c)
-{
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] *= c;
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::operator*=(complex_t c)
-{
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] *= c;
-   }
-
-   return *this;
-}
-
-ComplexTypeDenseMatrix &ComplexTypeDenseMatrix::Set(const DenseMatrix &Mr,
-                                                    const DenseMatrix &Mi)
-{
-   MFEM_ASSERT(height == Mr.Height() && height == Mi.Height() &&
-               width == Mr.Width() && width == Mi.Width(),
-               "incompatible Matrices!");
-
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      data[i] = complex_t(Mr.data[i], Mi.data[i]);
-   }
-
-   return *this;
-}
-
-void ComplexTypeDenseMatrix::Swap(ComplexTypeDenseMatrix &other)
-{
-   mfem::Swap(width, other.width);
-   mfem::Swap(height, other.height);
-   mfem::Swap(data, other.data);
-}
-
-ComplexTypeDenseMatrix::~ComplexTypeDenseMatrix()
-{
-   data.Delete();
-}
-
-const DenseMatrix &ComplexTypeDenseMatrix::real() const
-{
-   re_part.SetSize(height, width);
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      re_part.data[i] = data[i].real();
-   }
-
-   return re_part;
-}
-
-const DenseMatrix &ComplexTypeDenseMatrix::imag() const
-{
-   im_part.SetSize(height, width);
-   const int hw = height * width;
-   for (int i = 0; i < hw; i++)
-   {
-      im_part.data[i] = data[i].imag();
-   }
-
-   return im_part;
-}
-
 } // mfem namespace

linalg/complex_densemat.hpp

@@ -13,7 +13,6 @@
 #define MFEM_COMPLEX_DENSEMAT
 
 #include "complex_operator.hpp"
-#include "../general/complex_type.hpp"
 #include <complex>
 
 namespace mfem
@@ -242,220 +241,6 @@ public:
 
 };
 
-class ComplexTypeDenseMatrix
-{
-protected:
-   int height; ///< Dimension of the output / number of rows in the matrix.
-   int width;  ///< Dimension of the input / number of columns in the matrix.
-
-private:
-   Memory<complex_t > data;
-
-   mutable DenseMatrix re_part;
-   mutable DenseMatrix im_part;
-
-public:
-   /** Default constructor for DenseMatrix.
-       Sets data = NULL and height = width = 0. */
-   ComplexTypeDenseMatrix();
-
-   /// Copy constructor
-   ComplexTypeDenseMatrix(const ComplexTypeDenseMatrix &);
-   ComplexTypeDenseMatrix(const DenseMatrix &);
-
-   /// Creates square matrix of size s.
-   explicit ComplexTypeDenseMatrix(int s);
-
-   /// Creates rectangular matrix of size m x n.
-   ComplexTypeDenseMatrix(int m, int n);
-
-   /// Construct a ComplexTypeDenseMatrix using an existing data array.
-   /** The ComplexTypeDenseMatrix does not assume ownership of the data array,
-       i.e. it will not delete the array. */
-   ComplexTypeDenseMatrix(complex_t *d, int h, int w)
-      : height(h), width(w) { UseExternalData(d, h, w); }
-
-   /// Create a dense matrix using a braced initializer list
-   /// The inner lists correspond to rows of the matrix
-   template <int M, int N, typename T = real_t>
-   explicit ComplexTypeDenseMatrix(const T (&values)[M][N]) :
-      ComplexTypeDenseMatrix(
-         M, N)
-   {
-      // DenseMatrix is column-major so copies have to be element-wise
-      for (int i = 0; i < M; i++)
-      {
-         for (int j = 0; j < N; j++)
-         {
-            (*this)(i,j) = values[i][j];
-         }
-      }
-   }
-
-   /// Change the data array and the size of the DenseMatrix.
-   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
-       not delete the data array @a d. This method should not be used with
-       DenseMatrix that owns its current data array. */
-   void UseExternalData(complex_t *d, int h, int w)
-   {
-      data.Wrap(d, h*w, false);
-      height = h; width = w;
-   }
-
-   /// Change the data array and the size of the DenseMatrix.
-   /** The DenseMatrix does not assume ownership of the data array, i.e. it will
-       not delete the new array @a d. This method will delete the current data
-       array, if owned. */
-   void Reset(complex_t *d, int h, int w)
-   { if (OwnsData()) { data.Delete(); } UseExternalData(d, h, w); }
-
-   /** Clear the data array and the dimensions of the DenseMatrix. This method
-       should not be used with DenseMatrix that owns its current data array. */
-   void ClearExternalData() { data.Reset(); height = width = 0; }
-
-   /// Delete the matrix data array (if owned) and reset the matrix state.
-   void Clear()
-   { if (OwnsData()) { data.Delete(); } ClearExternalData(); }
-
-   /// Get the height (size of output) of the Operator. Synonym with NumRows().
-   inline int Height() const { return height; }
-   /** @brief Get the number of rows (size of output) of the Operator. Synonym
-       with Height(). */
-   inline int NumRows() const { return height; }
-
-   /// Get the width (size of input) of the Operator. Synonym with NumCols().
-   inline int Width() const { return width; }
-   /** @brief Get the number of columns (size of input) of the Operator. Synonym
-       with Width(). */
-   inline int NumCols() const { return width; }
-
-   /// For backward compatibility define Size to be synonym of Width()
-   int Size() const { return Width(); }
-
-   // Total size = width*height
-   int TotalSize() const { return width*height; }
-
-   /// Change the size of the DenseMatrix to s x s.
-   void SetSize(int s) { SetSize(s, s); }
-
-   /// Change the size of the DenseMatrix to h x w.
-   void SetSize(int h, int w);
-
-   /// Returns the matrix data array.
-   inline complex_t *Data() const
-   {
-      return const_cast<complex_t*>
-             ((const complex_t*)data);
-   }
-
-   /// Returns the matrix data array.
-   inline complex_t *GetData() const { return Data(); }
-
-   Memory<complex_t > &GetMemory() { return data; }
-   const Memory<complex_t > &GetMemory() const { return data; }
-
-   /// Return the DenseMatrix data (host pointer) ownership flag.
-   inline bool OwnsData() const { return data.OwnsHostPtr(); }
-
-   /// Returns reference to a_{ij}.
-   inline complex_t &operator()(int i, int j);
-
-   /// Returns constant reference to a_{ij}.
-   inline const complex_t &operator()(int i, int j) const;
-
-   /// Returns reference to a_{ij}.
-   complex_t &Elem(int i, int j);
-
-   /// Returns constant reference to a_{ij}.
-   const complex_t &Elem(int i, int j) const;
-
-   /// Sets the matrix elements equal to constant c
-   ComplexTypeDenseMatrix &operator=(real_t c);
-   ComplexTypeDenseMatrix &operator=(complex_t c);
-
-   /// Copy the matrix entries from the given array
-   ComplexTypeDenseMatrix &operator=(const real_t *d);
-   ComplexTypeDenseMatrix &operator=(const complex_t *d);
-
-   /// Sets the matrix size and elements equal to those of m
-   ComplexTypeDenseMatrix &operator=(const DenseMatrix &m);
-   ComplexTypeDenseMatrix &operator=(const ComplexTypeDenseMatrix &m);
-
-   ComplexTypeDenseMatrix &operator+=(const real_t *m);
-   ComplexTypeDenseMatrix &operator+=(const complex_t *m);
-   ComplexTypeDenseMatrix &operator+=(const DenseMatrix &m);
-   ComplexTypeDenseMatrix &operator+=(const ComplexTypeDenseMatrix &m);
-
-   ComplexTypeDenseMatrix &operator-=(const DenseMatrix &m);
-   ComplexTypeDenseMatrix &operator-=(const ComplexTypeDenseMatrix &m);
-
-   ComplexTypeDenseMatrix &operator*=(real_t c);
-   ComplexTypeDenseMatrix &operator*=(complex_t c);
-
-   /// (*this) = x + i * y
-   ComplexTypeDenseMatrix &Set(const DenseMatrix &x, const DenseMatrix &y);
-
-   std::size_t MemoryUsage() const
-   { return data.Capacity() * sizeof(complex_t); }
-
-   /// Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
-   const complex_t *Read(bool on_dev = true) const
-   { return mfem::Read(data, Height()*Width(), on_dev); }
-
-   /// Shortcut for mfem::Read(GetMemory(), TotalSize(), false).
-   const complex_t *HostRead() const
-   { return mfem::Read(data, Height()*Width(), false); }
-
-   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
-   complex_t *Write(bool on_dev = true)
-   { return mfem::Write(data, Height()*Width(), on_dev); }
-
-   /// Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
-   complex_t *HostWrite()
-   { return mfem::Write(data, Height()*Width(), false); }
-
-   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), on_dev).
-   complex_t *ReadWrite(bool on_dev = true)
-   { return mfem::ReadWrite(data, Height()*Width(), on_dev); }
-
-   /// Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
-   complex_t *HostReadWrite()
-   { return mfem::ReadWrite(data, Height()*Width(), false); }
-
-   void Swap(ComplexTypeDenseMatrix &other);
-
-   /// Return a reference to the real part of this matrix
-   const DenseMatrix &real() const;
-
-   /// Return a reference to the imaginary part of this matrix
-   const DenseMatrix &imag() const;
-
-   /// Destroys dense matrix.
-   virtual ~ComplexTypeDenseMatrix();
-};
-
-/// Specialization of the template function Swap<> for class ComplexTypeDenseMatrix
-template<> inline void Swap<ComplexTypeDenseMatrix>(ComplexTypeDenseMatrix &a,
-                                                    ComplexTypeDenseMatrix &b)
-{
-   a.Swap(b);
-}
-
-// Inline methods
-
-inline complex_t &ComplexTypeDenseMatrix::operator()(int i, int j)
-{
-   MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
-   return data[i+j*height];
-}
-
-inline const complex_t &ComplexTypeDenseMatrix::operator()
-(int i, int j) const
-{
-   MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
-   return data[i+j*height];
-}
-
 } // namespace mfem
 
 #endif // MFEM_COMPLEX_DENSEMAT

linalg/complex_vector.cpp

@@ -1,424 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "../general/forall.hpp"
-#include "../general/reducers.hpp"
-#include "complex_vector.hpp"
-
-using namespace std;
-
-namespace mfem
-{
-
-ComplexVector::ComplexVector(const ComplexVector &v)
-{
-   const int s = v.Size();
-   size = s;
-   if (s > 0)
-   {
-      MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
-      data.New(s, v.data.GetMemoryType());
-      data.CopyFrom(v.data, s);
-   }
-   UseDevice(v.UseDevice());
-}
-
-ComplexVector::ComplexVector(const Vector &v)
-{
-   const int s = v.Size();
-   size = s;
-   if (s > 0)
-   {
-      MFEM_ASSERT(!v.data.Empty(), "invalid source vector");
-      data.New(s, v.data.GetMemoryType());
-      MFEM_FORALL(i, size, data[i] = v.data[i]; );
-   }
-   UseDevice(v.UseDevice());
-}
-
-ComplexVector::ComplexVector(ComplexVector &&v)
-{
-   *this = std::move(v);
-}
-
-complex_t &ComplexVector::Elem(int i)
-{
-   return operator()(i);
-}
-
-const complex_t &ComplexVector::Elem(int i) const
-{
-   return operator()(i);
-}
-
-complex_t ComplexVector::operator*(const complex_t *v) const
-{
-   HostRead();
-   complex_t dot = 0.0;
-#ifdef MFEM_USE_LEGACY_OPENMP
-   #pragma omp parallel for reduction(+:dot)
-#endif
-   for (int i = 0; i < size; i++)
-   {
-      dot += data[i] * v[i];
-   }
-   return dot;
-}
-
-complex_t ComplexVector::operator*(const real_t *v) const
-{
-   HostRead();
-   complex_t dot = 0.0;
-#ifdef MFEM_USE_LEGACY_OPENMP
-   #pragma omp parallel for reduction(+:dot)
-#endif
-   for (int i = 0; i < size; i++)
-   {
-      dot += data[i] * v[i];
-   }
-   return dot;
-}
-
-complex_t ComplexVector::operator*(const ComplexVector &v) const
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   if (size == 0) { return 0.0; }
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const auto m_data = Read(use_dev), v_data = v.Read(use_dev);
-
-   // The standard way of computing the dot product is non-deterministic
-   complex_t prod = 0.0;
-   for (int i = 0; i < size; i++)
-   {
-      prod += m_data[i] * v_data[i];
-   }
-   return prod;
-}
-
-complex_t ComplexVector::operator*(const Vector &v) const
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   if (size == 0) { return 0.0; }
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const auto m_data = Read(use_dev);
-   const auto v_data = v.Read(use_dev);
-
-   // The standard way of computing the dot product is non-deterministic
-   complex_t prod = 0.0;
-   for (int i = 0; i < size; i++)
-   {
-      prod += m_data[i] * v_data[i];
-   }
-   return prod;
-}
-
-ComplexVector &ComplexVector::operator=(const complex_t *v)
-{
-   HostRead();
-   MFEM_FORALL(i, size, data[i] = v[i]; );
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator=(const real_t *v)
-{
-   HostRead();
-   MFEM_FORALL(i, size, data[i] = v[i]; );
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator=(const ComplexVector &v)
-{
-#if 0
-   SetSize(v.Size(), v.data.GetMemoryType());
-   data.CopyFrom(v.data, v.Size());
-   UseDevice(v.UseDevice());
-#else
-   SetSize(v.Size());
-   const bool vuse = v.UseDevice();
-   const bool use_dev = UseDevice() || vuse;
-   v.UseDevice(use_dev);
-   // keep 'data' where it is, unless 'use_dev' is true
-   if (use_dev) { Write(); }
-   data.CopyFrom(v.data, v.Size());
-   v.UseDevice(vuse);
-#endif
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator=(const Vector &v)
-{
-   SetSize(v.Size());
-   const bool vuse = v.UseDevice();
-   const bool use_dev = UseDevice() || vuse;
-   v.UseDevice(use_dev);
-   // keep 'data' where it is, unless 'use_dev' is true
-   if (use_dev) { Write(); }
-   MFEM_FORALL(i, size, data[i] = v[i]; );
-   v.UseDevice(vuse);
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator=(ComplexVector &&v)
-{
-   v.Swap(*this);
-   if (this != &v) { v.Destroy(); }
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator=(complex_t value)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] = value; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator=(real_t value)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] = value; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator*=(complex_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] *= c; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator*=(real_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] *= c; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator*=(const ComplexVector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] *= x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator*=(const Vector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] *= x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator/=(complex_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   const complex_t m = conj(c) / norm(c);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] *= m; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator/=(real_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   const real_t m = 1.0/c;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] *= m; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator/=(const ComplexVector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] /= x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator/=(const Vector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] /= x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator-=(complex_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] -= c; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator-=(real_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] -= c; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator-=(const ComplexVector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] -= x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator-=(const Vector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] -= x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator+=(complex_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] += c; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator+=(real_t c)
-{
-   const bool use_dev = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] += c; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator+=(const ComplexVector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] += x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::operator+=(const Vector &v)
-{
-   MFEM_ASSERT(size == v.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || v.UseDevice();
-   const int N = size;
-   const auto x = v.Read(use_dev);
-   auto y = ReadWrite(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] += x[i]; });
-   return *this;
-}
-
-ComplexVector &ComplexVector::Set(const Vector &Vr, const Vector &Vi)
-{
-   MFEM_ASSERT(size == Vr.size && size == Vi.size, "incompatible Vectors!");
-
-   const bool use_dev = UseDevice() || Vr.UseDevice() || Vi.UseDevice();
-   const int N = size;
-   const auto x = Vr.Read(use_dev);
-   const auto y = Vi.Read(use_dev);
-   auto z = Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { z[i] = complex_t(x[i], y[i]); });
-   return *this;
-}
-
-const Vector &ComplexVector::real() const
-{
-   re_part.SetSize(size);
-   const bool use_dev = UseDevice();
-   const int N = size;
-   const auto z = Read(use_dev);
-   auto x = re_part.Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { x[i] = z[i].real(); });
-   return re_part;
-}
-
-const Vector &ComplexVector::imag() const
-{
-   im_part.SetSize(size);
-   const bool use_dev = UseDevice();
-   const int N = size;
-   const auto z = Read(use_dev);
-   auto y = im_part.Write(use_dev);
-   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i)
-   { y[i] = z[i].imag(); });
-   return im_part;
-}
-
-}

linalg/complex_vector.hpp

@@ -1,479 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_COMPLEX_VECTOR
-#define MFEM_COMPLEX_VECTOR
-
-#include "vector.hpp"
-#include "../general/complex_type.hpp"
-
-namespace mfem
-{
-
-class ComplexVector
-{
-private:
-
-   Memory<complex_t > data;
-   int size;
-
-   mutable Vector re_part;
-   mutable Vector im_part;
-
-public:
-
-   /// Default constructor for ComplexVector. Sets size = 0
-   ComplexVector() : size(0) { }
-
-   /// Copy constructor. Allocates a new data array and copies the data.
-   ComplexVector(const ComplexVector &);
-
-   /// Copy constructor. Allocates a new data array and copies the
-   /// data into real part of this vector.
-   ComplexVector(const Vector &);
-
-   /// Move constructor. "Steals" data from its argument.
-   ComplexVector(ComplexVector&& v);
-
-   /// @brief Creates vector of size s.
-   /// @warning Entries are not initialized to zero!
-   explicit ComplexVector(int s);
-
-   /// Creates a vector referencing an array of complex<doubles>,
-   /// owned by someone else.
-   /// The pointer @a data_ can be NULL. The data array can be replaced later
-   /// with SetData().
-   ComplexVector(complex_t *data_, int size_)
-   { data.Wrap(data_, size_, false); size = size_; }
-
-   /// @brief Create a ComplexVector referencing a sub-vector of the
-   //  ComplexVector @a base starting at the given offset, @a
-   //  base_offset, and size @a size_.
-   ComplexVector(ComplexVector &base, int base_offset, int size_)
-      : data(base.data, base_offset, size_), size(size_) { }
-
-   /// Create a ComplexVector of size @a size_ using MemoryType @a mt.
-   ComplexVector(int size_, MemoryType mt)
-      : data(size_, mt), size(size_) { }
-
-   /// @brief Create a ComplexVector of size @a size_ using host
-   /// MemoryType @a h_mt and device MemoryType @a d_mt.
-   ComplexVector(int size_, MemoryType h_mt, MemoryType d_mt)
-      : data(size_, h_mt, d_mt), size(size_) { }
-
-   /// Create a vector from a statically sized C-style array of convertible type
-   template <typename CT, int N>
-   explicit ComplexVector(const CT (&values)[N]) : ComplexVector(N)
-   { std::copy(values, values + N, begin()); }
-
-   /// Create a vector using a braced initializer list
-   template <typename CT, typename std::enable_if<
-                std::is_convertible<CT,complex_t >::value,bool>::type = true>
-   explicit ComplexVector(std::initializer_list<CT> values) : ComplexVector(
-         values.size())
-   { std::copy(values.begin(), values.end(), begin()); }
-
-   /// Enable execution of Vector operations using the mfem::Device.
-   /// The default is to use Backend::CPU (serial execution on each MPI rank),
-   /// regardless of the mfem::Device configuration.
-   ///
-   /// When appropriate, MFEM functions and class methods will enable the use
-   /// of the mfem::Device for their Vector parameters.
-   ///
-   /// Some derived classes, e.g. GridFunction, enable the use of the
-   /// mfem::Device by default.
-   virtual void UseDevice(bool use_dev) const { data.UseDevice(use_dev); }
-
-   /// Return the device flag of the Memory object used by the Vector
-   virtual bool UseDevice() const { return data.UseDevice(); }
-
-   /// @brief Resize the vector to size @a s.
-   /// If the new size is less than or equal to Capacity() then the internal
-   /// data array remains the same. Otherwise, the old array is deleted, if
-   /// owned, and a new array of size @a s is allocated without copying the
-   /// previous content of the ComplexVector.
-   /// @warning In the second case above (new size greater than current one),
-   /// the vector will allocate new data array, even if it did not own the
-   /// original data! Also, new entries are not initialized!
-   void SetSize(int s);
-
-   /// Resize the vector to size @a s using MemoryType @a mt.
-   void SetSize(int s, MemoryType mt);
-
-   /// Resize the vector to size @a s using the MemoryType of @a v.
-   void SetSize(int s, const ComplexVector &v)
-   { SetSize(s, v.GetMemory().GetMemoryType()); }
-
-   /// Resize the vector to size @a s using the MemoryType of @a v.
-   void SetSize(int s, const Vector &v)
-   { SetSize(s, v.GetMemory().GetMemoryType()); }
-
-   /// Set the Vector data.
-   /// @warning This method should be called only when OwnsData() is false.
-   void SetData(complex_t *d)
-   { data.Wrap(d, data.Capacity(), false); }
-
-   /// Set the Vector data and size.
-   /// The Vector does not assume ownership of the new data. The new size is
-   /// also used as the new Capacity().
-   /// @warning This method should be called only when OwnsData() is false.
-   /// @sa NewDataAndSize().
-   void SetDataAndSize(complex_t *d, int s)
-   { data.Wrap(d, s, false); size = s; }
-
-   /// Set the Vector data and size, deleting the old data, if owned.
-   /// The Vector does not assume ownership of the new data. The new size is
-   /// also used as the new Capacity().
-   /// @sa SetDataAndSize().
-   void NewDataAndSize(complex_t *d, int s)
-   {
-      data.Delete();
-      SetDataAndSize(d, s);
-   }
-
-   /// Reset the Vector to use the given external Memory @a mem and size @a s.
-   /// If @a own_mem is false, the Vector will not own any of the pointers of
-   /// @a mem.
-   ///
-   /// Note that when @a own_mem is true, the @a mem object can be destroyed
-   /// immediately by the caller but `mem.Delete()` should NOT be called since
-   /// the Vector object takes ownership of all pointers owned by @a mem.
-   ///
-   /// @sa NewDataAndSize().
-   inline void NewMemoryAndSize(const Memory<complex_t > &mem,
-                                int s, bool own_mem);
-
-   /// Reset the Vector to be a reference to a sub-vector of @a base.
-   inline void MakeRef(ComplexVector &base, int offset, int size);
-
-   /// @brief Reset the Vector to be a reference to a sub-vector of @a base
-   /// without changing its current size.
-   inline void MakeRef(ComplexVector &base, int offset);
-
-   /// Set the Vector data (host pointer) ownership flag.
-   void MakeDataOwner() const { data.SetHostPtrOwner(true); }
-
-   /// Destroy a vector
-   void Destroy();
-
-   /// @brief Delete the device pointer, if owned. If @a copy_to_host is true
-   /// and the data is valid only on device, move it to host before deleting.
-   /// Invalidates the device memory.
-   void DeleteDevice(bool copy_to_host = true)
-   { data.DeleteDevice(copy_to_host); }
-
-   /// Returns the size of the vector.
-   inline int Size() const { return size; }
-
-   /// Return the size of the currently allocated data array.
-   /// It is always true that Capacity() >= Size().
-   inline int Capacity() const { return data.Capacity(); }
-
-   /// Return a pointer to the beginning of the ComplexVector data.
-   /// @warning This method should be used with caution as it gives write access
-   /// to the data of const-qualified ComplexVector%s.
-   inline complex_t *GetData() const
-   { return const_cast<complex_t*>((const complex_t*)data); }
-
-   /// STL-like begin.
-   inline complex_t *begin() { return data; }
-
-   /// STL-like end.
-   inline complex_t *end() { return data + size; }
-
-   /// STL-like begin (const version).
-   inline const complex_t *begin() const { return data; }
-
-   /// STL-like end (const version).
-   inline const complex_t *end() const { return data + size; }
-
-   /// Return a reference to the Memory object used by the Vector.
-   Memory<complex_t > &GetMemory() { return data; }
-
-   /// @brief Return a reference to the Memory object used by the
-   /// ComplexVector, const version.
-   const Memory<complex_t > &GetMemory() const { return data; }
-
-   /// Update the memory location of the vector to match @a v.
-   void SyncMemory(const ComplexVector &v) const
-   { GetMemory().Sync(v.GetMemory()); }
-
-   /// Update the alias memory location of the vector to match @a v.
-   void SyncAliasMemory(const ComplexVector &v) const
-   { GetMemory().SyncAlias(v.GetMemory(),Size()); }
-
-   /// Read the Vector data (host pointer) ownership flag.
-   inline bool OwnsData() const { return data.OwnsHostPtr(); }
-
-   /// Changes the ownership of the data; after the call the Vector is empty
-   inline void StealData(complex_t **p)
-   { *p = data; data.Reset(); size = 0; }
-
-   /// Changes the ownership of the data; after the call the Vector is empty
-   inline complex_t *StealData()
-   { complex_t *p; StealData(&p); return p; }
-
-   /// Access Vector entries. Index i = 0 .. size-1.
-   complex_t &Elem(int i);
-
-   /// Read only access to Vector entries. Index i = 0 .. size-1.
-   const complex_t &Elem(int i) const;
-
-   /// Access Vector entries using () for 0-based indexing.
-   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
-   inline complex_t &operator()(int i);
-
-   /// Read only access to Vector entries using () for 0-based indexing.
-   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
-   inline const complex_t &operator()(int i) const;
-
-   /// Access Vector entries using [] for 0-based indexing.
-   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
-   inline complex_t &operator[](int i) { return (*this)(i); }
-
-   /// Read only access to Vector entries using [] for 0-based indexing.
-   /// @note If MFEM_DEBUG is enabled, bounds checking is performed.
-   inline const complex_t &operator[](int i) const
-   { return (*this)(i); }
-
-   /// Dot product with a `complex<double> *` array.
-   /// @note No complex conjugate is performed
-   complex_t operator*(const complex_t *v) const;
-   complex_t operator*(const real_t *v) const;
-
-   /// Return the inner-product.
-   /// @note No complex conjugate is performed
-   complex_t operator*(const ComplexVector &v) const;
-   complex_t operator*(const Vector &v) const;
-
-   /// Copy Size() entries from @a v.
-   ComplexVector &operator=(const complex_t *v);
-   ComplexVector &operator=(const real_t *v);
-
-   /// Copy assignment.
-   /// @note Defining this method overwrites the implicitly defined copy
-   /// assignment operator.
-   ComplexVector &operator=(const ComplexVector &v);
-   ComplexVector &operator=(const Vector &v);
-
-   /// Move assignment
-   ComplexVector &operator=(ComplexVector&& v);
-
-   /// Redefine '=' for vector = constant.
-   ComplexVector &operator=(complex_t value);
-   ComplexVector &operator=(real_t value);
-
-   /// Scale vector by a constant
-   ComplexVector &operator*=(complex_t c);
-   ComplexVector &operator*=(real_t c);
-
-   /// Component-wise scaling: (*this)(i) *= v(i)
-   ComplexVector &operator*=(const ComplexVector &v);
-   ComplexVector &operator*=(const Vector &v);
-
-   /// Divide vector by a consant
-   ComplexVector &operator/=(complex_t c);
-   ComplexVector &operator/=(real_t c);
-
-   /// Component-wise division: (*this)(i) /= v(i)
-   ComplexVector &operator/=(const ComplexVector &v);
-   ComplexVector &operator/=(const Vector &v);
-
-   /// Subtract a constant from this vector
-   ComplexVector &operator-=(complex_t c);
-   ComplexVector &operator-=(real_t c);
-
-   /// Subtract a vector from this vector
-   ComplexVector &operator-=(const ComplexVector &v);
-   ComplexVector &operator-=(const Vector &v);
-
-   /// Add a constant to this vector
-   ComplexVector &operator+=(complex_t c);
-   ComplexVector &operator+=(real_t c);
-
-   /// Add a vector to this vector
-   ComplexVector &operator+=(const ComplexVector &v);
-   ComplexVector &operator+=(const Vector &v);
-
-   /// (*this) = x + i * y
-   ComplexVector &Set(const Vector &x, const Vector &y);
-
-   /// Swap the contents of two Vectors
-   inline void Swap(ComplexVector &other);
-
-   /// Return a reference to the real part of this vector
-   const Vector &real() const;
-
-   /// Return a reference to the imaginary part of this vector
-   const Vector &imag() const;
-
-   /// Destroys vector.
-   virtual ~ComplexVector();
-
-   /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
-   virtual const complex_t *Read(bool on_dev = true) const
-   { return mfem::Read(data, size, on_dev); }
-
-   /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
-   virtual const complex_t *HostRead() const
-   { return mfem::Read(data, size, false); }
-
-   /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
-   virtual complex_t *Write(bool on_dev = true)
-   { return mfem::Write(data, size, on_dev); }
-
-   /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
-   virtual complex_t *HostWrite()
-   { return mfem::Write(data, size, false); }
-
-   /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
-   virtual complex_t *ReadWrite(bool on_dev = true)
-   { return mfem::ReadWrite(data, size, on_dev); }
-
-   /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
-   virtual complex_t *HostReadWrite()
-   { return mfem::ReadWrite(data, size, false); }
-};
-
-inline ComplexVector::ComplexVector(int s)
-{
-   MFEM_ASSERT(s>=0,"Unexpected negative size.");
-   size = s;
-   if (s > 0)
-   {
-      data.New(s);
-   }
-}
-
-inline void ComplexVector::SetSize(int s)
-{
-   if (s == size)
-   {
-      return;
-   }
-   if (s <= data.Capacity())
-   {
-      size = s;
-      return;
-   }
-   // preserve a valid MemoryType and device flag
-   const MemoryType mt = data.GetMemoryType();
-   const bool use_dev = data.UseDevice();
-   data.Delete();
-   size = s;
-   data.New(s, mt);
-   data.UseDevice(use_dev);
-}
-
-inline void ComplexVector::SetSize(int s, MemoryType mt)
-{
-   if (mt == data.GetMemoryType())
-   {
-      if (s == size)
-      {
-         return;
-      }
-      if (s <= data.Capacity())
-      {
-         size = s;
-         return;
-      }
-   }
-   const bool use_dev = data.UseDevice();
-   data.Delete();
-   if (s > 0)
-   {
-      data.New(s, mt);
-      size = s;
-   }
-   else
-   {
-      data.Reset();
-      size = 0;
-   }
-   data.UseDevice(use_dev);
-}
-
-inline void ComplexVector::NewMemoryAndSize(
-   const Memory<complex_t > &mem,
-   int s,
-   bool own_mem)
-{
-   data.Delete();
-   size = s;
-   if (own_mem)
-   {
-      data = mem;
-   }
-   else
-   {
-      data.MakeAlias(mem, 0, s);
-   }
-}
-
-inline void ComplexVector::MakeRef(ComplexVector &base, int offset, int s)
-{
-   data.Delete();
-   size = s;
-   data.MakeAlias(base.GetMemory(), offset, s);
-}
-
-inline void ComplexVector::MakeRef(ComplexVector &base, int offset)
-{
-   data.Delete();
-   data.MakeAlias(base.GetMemory(), offset, size);
-}
-
-inline void ComplexVector::Destroy()
-{
-   const bool use_dev = data.UseDevice();
-   data.Delete();
-   size = 0;
-   data.Reset();
-   data.UseDevice(use_dev);
-}
-
-inline complex_t &ComplexVector::operator()(int i)
-{
-   MFEM_ASSERT(data && i >= 0 && i < size,
-               "index [" << i << "] is out of range [0," << size << ")");
-
-   return data[i];
-}
-
-inline const complex_t &ComplexVector::operator()(int i) const
-{
-   MFEM_ASSERT(data && i >= 0 && i < size,
-               "index [" << i << "] is out of range [0," << size << ")");
-
-   return data[i];
-}
-
-inline void ComplexVector::Swap(ComplexVector &other)
-{
-   mfem::Swap(data, other.data);
-   mfem::Swap(size, other.size);
-}
-
-/// Specialization of the template function Swap<> for class ComplexVector
-template<> inline void Swap<ComplexVector>(ComplexVector &a, ComplexVector &b)
-{
-   a.Swap(b);
-}
-
-inline ComplexVector::~ComplexVector()
-{
-   data.Delete();
-}
-
-} // namespace mfem
-
-#endif

linalg/densemat.hpp

@@ -24,7 +24,6 @@ class DenseMatrix : public Matrix
 {
    friend class DenseTensor;
    friend class DenseMatrixInverse;
-   friend class ComplexTypeDenseMatrix;
 
 private:
    Memory<real_t> data;

linalg/vector.cpp

@@ -668,6 +668,28 @@ void Vector::median(const Vector &lo, const Vector &hi)
    });
 }
 
+void Vector::Insert(int offset, const Vector &sv)
+{
+   const int old_size = size;
+
+   if (sv.Size() + old_size > Capacity())
+   {
+      Vector copy = *this;
+      SetSize(size + sv.Size());
+      SetVector(copy, 0);
+   }
+   else
+   {
+      SetSize(size + sv.Size());
+   }
+
+   for (int j = old_size-1; j >= offset; j--)
+   {
+      data[j+sv.Size()] = data[j];
+   }
+   SetVector(sv, offset);
+}
+
 void Vector::GetSubVector(const Array<int> &dofs, Vector &elemvect) const
 {
    const int n = dofs.Size();

linalg/vector.hpp

@@ -80,8 +80,6 @@ inline real_t rand_real()
 /// Vector data type.
 class Vector
 {
-   friend class ComplexVector;
-
 protected:
 
    Memory<real_t> data;
@@ -173,6 +171,9 @@ public:
    /// Resize the vector to size @a s using the MemoryType of @a v.
    void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
 
+   /// Delete elements at @a indices and resize vector accordingly
+   void DeleteAt(const Array<int> &indices);
+
    /// Set the Vector data.
    /// @warning This method should be called only when OwnsData() is false.
    void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
@@ -398,6 +399,12 @@ public:
    /// v = median(v,lo,hi) entrywise.  Implementation assumes lo <= hi.
    void median(const Vector &lo, const Vector &hi);
 
+   /// Insert sub Vector @a sv at @a offset and resize
+   void Insert(int offset, const Vector &sv);
+
+   /// Insert @a value at @a offset and resize
+   void Insert(int offset, const real_t value) { Insert(offset, Vector({value})); }
+
    /// Extract entries listed in @a dofs to the output Vector @a elemvect.
    /** Negative dof values cause the -dof-1 position in @a elemvect to receive
        the -val in from this Vector. */
@@ -623,6 +630,29 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
+inline void Vector::DeleteAt(const Array<int> &indices)
+{
+   // Make copy of the indices, sorted.
+   Array<int> sorted_indices(indices);
+   sorted_indices.Sort();
+
+   int rm_count = 0;
+   for (int i = 0; i < size; i++)
+   {
+      if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
+      {
+         rm_count++;
+      }
+      else
+      {
+         data[i-rm_count] = data[i]; // shift data rm_count
+      }
+   }
+
+   // Resize to remove tail
+   SetSize(size - rm_count);
+}
+
 inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
                                      bool own_mem)
 {

mesh/mesh.cpp

@@ -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

mesh/mesh.hpp

@@ -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

miniapps/performance/makefile

@@ -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

tests/unit/CMakeLists.txt

@@ -32,7 +32,6 @@ set(UNIT_TESTS_SRCS
   linalg/test_chebyshev.cpp
   linalg/test_complex_dense_matrix.cpp
   linalg/test_complex_operator.cpp
-  linalg/test_complex_vector.cpp
   linalg/test_constrainedsolver.cpp
   linalg/test_direct_solvers.cpp
   linalg/test_hypre_ilu.cpp

tests/unit/fem/test_derefine.cpp

@@ -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)
    {

tests/unit/general/test_array.cpp

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

tests/unit/linalg/test_complex_vector.cpp

@@ -1,145 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-#include <numeric>
-
-using namespace mfem;
-
-TEST_CASE("Complex Vector init-list and C-style array constructors",
-          "[ComplexVector]")
-{
-   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
-                                         std::complex<real_t>(5.0,2.0),
-                                         std::complex<real_t>(4.0,3.0),
-                                         std::complex<real_t>(3.0,4.0),
-                                         std::complex<real_t>(2.0,5.0),
-                                         std::complex<real_t>(1.0,6.0)
-                                        };
-   // Point and size constructor
-   ComplexVector a(ContigData, 6);
-   // Braced-list constructor
-   ComplexVector b({std::complex<real_t>(6.0,1.0),
-                    std::complex<real_t>(5.0,2.0),
-                    std::complex<real_t>(4.0,3.0),
-                    std::complex<real_t>(3.0,4.0),
-                    std::complex<real_t>(2.0,5.0),
-                    std::complex<real_t>(1.0,6.0)});
-   // Statically sized C-style array constructor
-   ComplexVector c(ContigData);
-
-   for (int i = 0; i < a.Size(); i++)
-   {
-      REQUIRE(a[i] == b[i]);
-      REQUIRE(a[i] == c[i]);
-   }
-}
-
-TEST_CASE("Complex Vector Move Constructor", "[ComplexVector]")
-{
-   constexpr int N = 6;
-   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
-                                         std::complex<real_t>(5.0,2.0),
-                                         std::complex<real_t>(4.0,3.0),
-                                         std::complex<real_t>(3.0,4.0),
-                                         std::complex<real_t>(2.0,5.0),
-                                         std::complex<real_t>(1.0,6.0)
-                                        };
-   ComplexVector a(ContigData, N);
-   ComplexVector b(N);
-   for (int i = 0; i < N; i++)
-   {
-      b(i) = std::complex<real_t>(N - i, i + 1);
-   }
-
-   std::complex<real_t>* a_data = a.GetData();
-   std::complex<real_t>* b_data = b.GetData();
-
-   ComplexVector move_non_owning(std::move(a));
-   ComplexVector move_owning(std::move(b));
-
-   REQUIRE(a.Size() == 0);
-   REQUIRE(a.GetData() == nullptr);
-   REQUIRE(b.Size() == 0);
-   REQUIRE(b.GetData() == nullptr);
-
-   // Should both be no-ops
-   a.Destroy();
-   b.Destroy();
-
-   REQUIRE(move_non_owning.OwnsData() == false);
-   REQUIRE(move_owning.OwnsData() == true);
-
-   REQUIRE(move_non_owning.Size() == N);
-   REQUIRE(move_owning.Size() == N);
-
-   // Make sure that the pointers were reused
-   REQUIRE(move_non_owning.GetData() == a_data);
-   REQUIRE(move_owning.GetData() == b_data);
-
-   for (int i = 0; i < N; i++)
-   {
-      REQUIRE(move_non_owning(i) == std::complex<real_t>(N - i, i + 1));
-      REQUIRE(move_owning(i) == std::complex<real_t>(N - i, i + 1));
-   }
-}
-
-TEST_CASE("Complex Vector Move Assignment", "[ComplexVector]")
-{
-   constexpr int N = 6;
-   std::complex<real_t> ContigData[6] = {std::complex<real_t>(6.0,1.0),
-                                         std::complex<real_t>(5.0,2.0),
-                                         std::complex<real_t>(4.0,3.0),
-                                         std::complex<real_t>(3.0,4.0),
-                                         std::complex<real_t>(2.0,5.0),
-                                         std::complex<real_t>(1.0,6.0)
-                                        };
-   ComplexVector a(ContigData, N);
-   ComplexVector b(N);
-   for (int i = 0; i < N; i++)
-   {
-      b(i) = std::complex<real_t>(N - i, i + 1);
-   }
-
-   std::complex<real_t>* a_data = a.GetData();
-   std::complex<real_t>* b_data = b.GetData();
-
-   ComplexVector move_non_owning;
-   move_non_owning = std::move(a);
-   ComplexVector move_owning;
-   move_owning = std::move(b);
-
-   REQUIRE(a.Size() == 0);
-   REQUIRE(a.GetData() == nullptr);
-   REQUIRE(b.Size() == 0);
-   REQUIRE(b.GetData() == nullptr);
-
-   // Should both be no-ops
-   a.Destroy();
-   b.Destroy();
-
-   REQUIRE(move_non_owning.OwnsData() == false);
-   REQUIRE(move_owning.OwnsData() == true);
-
-   REQUIRE(move_non_owning.Size() == N);
-   REQUIRE(move_owning.Size() == N);
-
-   // Make sure that the pointers were reused
-   REQUIRE(move_non_owning.GetData() == a_data);
-   REQUIRE(move_owning.GetData() == b_data);
-
-   for (int i = 0; i < N; i++)
-   {
-      REQUIRE(move_non_owning(i) == std::complex<real_t>(N - i, i + 1));
-      REQUIRE(move_owning(i) == std::complex<real_t>(N - i, i + 1));
-   }
-}

tests/unit/linalg/test_vector.cpp

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