resolving conflicts with master

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>

CHANGELOG

@@ -68,6 +68,10 @@ Linear and nonlinear solvers
   PRefinement multigrid methods for problems posed on trace spaces (see e.g. the
   DPG miniapps).
 
+- Added interface to MUMPS direct solver for complex-valued problems.
+  Its usage is demonstrated in ex25p. See http://mumps.enseeiht.fr/ for more details.
+  Supported versions >= 5.1.1.
+
 GPU computing
 -------------
 - Added device assembly support for 3D H(curl) VectorFEDomainLFIntegrator.

CMakeLists.txt

@@ -76,7 +76,9 @@ set(XSDK_ENABLE_Fortran OFF)
 # Check if we need to enable C or Fortran.
 if (MFEM_USE_CONDUIT OR
     MFEM_USE_SIDRE OR
-    MFEM_USE_PETSC)
+    MFEM_USE_PETSC OR
+    MFEM_USE_MUMPS OR
+    MFEM_USE_COMPLEX_MUMPS)
    # This seems to be needed by:
    #  * find_package(BLAS REQUIRED) and
    #  * find_package(HDF5 REQUIRED) needed, in turn, by:
@@ -84,7 +86,7 @@ if (MFEM_USE_CONDUIT OR
    #  * find_package(PETSc REQUIRED)
    set(XSDK_ENABLE_C ON)
 endif()
-if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS)
+if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS OR MFEM_USE_COMPLEX_MUMPS)
   # Just needed to find the MPI_Fortran libraries to link with
   set(XSDK_ENABLE_Fortran ON)
 endif()
@@ -342,7 +344,7 @@ if (MFEM_USE_MPI)
     endif()
   endif()
 else()
-  set(PKGS_NEED_MPI SUPERLU MUMPS PETSC SLEPC STRUMPACK PUMI)
+  set(PKGS_NEED_MPI SUPERLU MUMPS COMPLEX_MUMPS PETSC SLEPC STRUMPACK PUMI)
   foreach(PKG IN LISTS PKGS_NEED_MPI)
     if (MFEM_USE_${PKG})
       message(STATUS "Disabling package ${PKG} - requires MPI")
@@ -415,7 +417,7 @@ if (MFEM_USE_SUPERLU)
 endif()
 
 # MUMPS can only be enabled in parallel
-if (MFEM_USE_MUMPS)
+if (MFEM_USE_MUMPS OR MFEM_USE_COMPLEX_MUMPS)
   if (MFEM_USE_MPI)
     find_package(MUMPS REQUIRED mumps_common pord)
     set(MFEM_MUMPS_VERSION ${MUMPS_VERSION})

config/cmake/MFEMConfig.cmake.in

@@ -34,6 +34,7 @@ set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
 set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
+set(MFEM_USE_COMPLEX_MUMPS @MFEM_USE_COMPLEX_MUMPS@)
 set(MFEM_USE_STRUMPACK @MFEM_USE_STRUMPACK@)
 set(MFEM_USE_CUDSS @MFEM_USE_CUDSS@)
 set(MFEM_USE_GINKGO @MFEM_USE_GINKGO@)

config/cmake/config.hpp.in

@@ -103,6 +103,7 @@
 
 // Enable MFEM functionality based on the MUMPS library.
 #cmakedefine MFEM_USE_MUMPS
+#cmakedefine MFEM_USE_COMPLEX_MUMPS
 #cmakedefine MFEM_MUMPS_VERSION @MFEM_MUMPS_VERSION@
 
 // Enable MFEM functionality based on the STRUMPACK library.

config/cmake/modules/FindMUMPS.cmake

@@ -17,13 +17,31 @@
 
 include(MfemCmakeUtilities)
 
-# Toggle which precision of MUMPS to use depending on the precision of MFEM.
+# Decide headers/libs by MFEM precision
 if (MFEM_USE_DOUBLE)
-  set(_mumps_header dmumps_c.h)
-  set(_mumps_lib dmumps)
-elseif(MFEM_USE_SINGLE)
-  set(_mumps_header smumps_c.h)
-  set(_mumps_lib smumps)
+  set(_rmumps_header dmumps_c.h)
+  set(_rmumps_lib    dmumps)
+  set(_cmumps_header zmumps_c.h)
+  set(_cmumps_lib    zmumps)
+elseif (MFEM_USE_SINGLE)
+  set(_rmumps_header smumps_c.h)
+  set(_rmumps_lib    smumps)
+  set(_cmumps_header cmumps_c.h)
+  set(_cmumps_lib    cmumps)
+endif()
+
+# Choose which header/lib mfem_find_package should use as the "primary" one.
+# If both enabled, prefer the real one as primary (either is fine).
+if (MFEM_USE_MUMPS)
+  set(_mumps_header ${_rmumps_header})
+  set(_mumps_lib    ${_rmumps_lib})
+elseif (MFEM_USE_COMPLEX_MUMPS)
+  set(_mumps_header ${_cmumps_header})
+  set(_mumps_lib    ${_cmumps_lib})
+else()
+  # Should not happen in practice because FindMUMPS is only called when enabled,
+  set(_mumps_header ${_rmumps_header})
+  set(_mumps_lib    ${_rmumps_lib})
 endif()
 
 mfem_find_package(MUMPS MUMPS MUMPS_DIR
@@ -31,8 +49,35 @@ mfem_find_package(MUMPS MUMPS MUMPS_DIR
   "Paths to headers required by MUMPS."
   "Libraries required by MUMPS."
   ADD_COMPONENT mumps_common "include" ${_mumps_header} "lib" mumps_common
-  ADD_COMPONENT pord "include" ${_mumps_header} "lib" pord)
+  ADD_COMPONENT pord        "include" ${_mumps_header} "lib" pord)
 
+# If BOTH real and complex are enabled, ensure BOTH solver libs are linked.
+if (MUMPS_FOUND AND MFEM_USE_MUMPS AND MFEM_USE_COMPLEX_MUMPS)
+  # Find the "other" solver library and append it.
+  find_library(_mfem_other_mumps_solver
+    NAMES ${_cmumps_lib}   
+    HINTS ${MUMPS_DIR}
+    PATH_SUFFIXES lib lib64
+    NO_DEFAULT_PATH)
+
+  if (NOT _mfem_other_mumps_solver)
+    # Fall back to system search
+    find_library(_mfem_other_mumps_solver NAMES ${_cmumps_lib})
+  endif()
+
+  if (NOT _mfem_other_mumps_solver)
+    message(FATAL_ERROR
+      "MFEM_USE_MUMPS=ON and MFEM_USE_COMPLEX_MUMPS=ON, but could not find "
+      "the complex solver library '${_cmumps_lib}' in MUMPS_DIR='${MUMPS_DIR}'.")
+  endif()
+
+  # Put solver libs first (important for static link order)
+  # MUMPS_LIBRARIES contains the primary solver already + common + pord.
+  # We prepend the other solver.
+  list(INSERT MUMPS_LIBRARIES 0 ${_mfem_other_mumps_solver})
+endif()
+
+# Version detection 
 if (MUMPS_FOUND AND (NOT MUMPS_VERSION))
   try_run(MUMPS_VERSION_RUN_RESULT MUMPS_VERSION_COMPILE_RESULT
           ${CMAKE_CURRENT_BINARY_DIR}/config

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -876,12 +876,12 @@ function(mfem_export_mk_files)
       MFEM_USE_ZLIB MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE
       MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS
       MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
-      MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_MAGMA
-      MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC
-      MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
-      MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
-      MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
-      MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
+      MFEM_USE_COMPLEX_MUMPS MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX
+      MFEM_USE_MAGMA MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC
+      MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT
+      MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP
+      MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE
+      MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
       MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG
       MFEM_USE_TRIBOL MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME
       MFEM_USE_HDF5)

config/config.hpp

@@ -140,6 +140,9 @@ constexpr real_t operator""_r(unsigned long long v)
 #ifdef MFEM_USE_MUMPS
 #error Building with MUMPS (MFEM_USE_MUMPS=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif
+#ifdef MFEM_USE_COMPLEX_MUMPS
+#error Building with COMPLEX_MUMPS (MFEM_USE_COMPLEX_MUMPS=YES) requires MPI (MFEM_USE_MPI=YES)
+#endif
 #ifdef MFEM_USE_STRUMPACK
 #error Building with STRUMPACK (MFEM_USE_STRUMPACK=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif

config/config.hpp.in

@@ -103,6 +103,7 @@
 
 // Enable MFEM functionality based on the MUMPS library.
 // #define MFEM_USE_MUMPS
+// #define MFEM_USE_COMPLEX_MUMPS
 // #define MFEM_MUMPS_VERSION @MFEM_MUMPS_VERSION@
 
 // Enable MFEM functionality based on the STRUMPACK library.

config/defaults.cmake

@@ -37,6 +37,7 @@ option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
 option(MFEM_USE_MUMPS "Enable MUMPS usage" OFF)
+option(MFEM_USE_COMPLEX_MUMPS "Enable COMPLEX_MUMPS usage" OFF)
 option(MFEM_USE_STRUMPACK "Enable STRUMPACK usage" OFF)
 option(MFEM_USE_CUDSS "Enable cuDSS usage" OFF)
 option(MFEM_USE_GINKGO "Enable Ginkgo usage" OFF)
@@ -153,7 +154,7 @@ set(SuperLUDist_REQUIRED_PACKAGES "MPI" "ParMETIS" "METIS"
 set(MUMPS_DIR "${MFEM_DIR}/../MUMPS_5.5.0" CACHE PATH
     "Path to the MUMPS library.")
 # MUMPS may also depend on "OpenMP", depending on how it was compiled.
-set(MUMPS_REQUIRED_PACKAGES "MPI" "MPI_Fortran" "ParMETIS" "METIS"
+set(MUMPS_REQUIRED_PACKAGES "MPI" "MPI_Fortran" "METIS"
     "ScaLAPACK" "LAPACK" "BLAS" CACHE STRING
     "Additional packages required by MUMPS.")
 # If the MPI package does not find all required Fortran libraries:

config/defaults.mk

@@ -152,6 +152,7 @@ MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
 MFEM_USE_MUMPS         = NO
+MFEM_USE_COMPLEX_MUMPS = NO
 MFEM_USE_STRUMPACK     = NO
 MFEM_USE_CUDSS         = NO
 MFEM_USE_GINKGO        = NO
@@ -249,7 +250,7 @@ ifeq (YES,$(MFEM_USE_HIP))
 endif
 
 # METIS library configuration
-ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
+ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS)$(MFEM_USE_COMPLEX_MUMPS),NONONONO)
    ifeq ($(MFEM_USE_METIS_5),NO)
      METIS_DIR = @MFEM_DIR@/../metis-4.0
      METIS_OPT =
@@ -353,13 +354,19 @@ MPI_FORTRAN_LIB = -lmpifort
 # MUMPS library configuration
 MUMPS_DIR = @MFEM_DIR@/../MUMPS_5.5.0
 MUMPS_OPT = -I$(MUMPS_DIR)/include
-MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib
+COMPLEX_MUMPS_OPT = $(MUMPS_OPT)
+
+MUMPS_COMMON_PATH = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib
+MUMPS_COMMON_LIB = -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
 ifeq ($(MFEM_USE_SINGLE),YES)
-   MUMPS_LIB += -lsmumps
+  MUMPS_SOLVER_LIB         = -lsmumps
+  COMPLEX_MUMPS_SOLVER_LIB = -lcmumps
 else
-   MUMPS_LIB += -ldmumps
+  MUMPS_SOLVER_LIB         = -ldmumps
+  COMPLEX_MUMPS_SOLVER_LIB = -lzmumps
 endif
-MUMPS_LIB += -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
+MUMPS_LIB         = $(MUMPS_COMMON_PATH) $(MUMPS_SOLVER_LIB) $(MUMPS_COMMON_LIB)
+COMPLEX_MUMPS_LIB = $(MUMPS_COMMON_PATH) $(COMPLEX_MUMPS_SOLVER_LIB) $(MUMPS_COMMON_LIB)
 
 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build

config/get_mumps_version.cpp

@@ -10,10 +10,18 @@
 // CONTRIBUTING.md for details.
 
 #ifdef MFEM_USE_SINGLE
+#ifdef MFEM_USE_COMPLEX_MUMPS
+#include "cmumps_c.h"
+#else
 #include "smumps_c.h"
+#endif
+#else
+#ifdef MFEM_USE_COMPLEX_MUMPS
+#include "zmumps_c.h"
 #else
 #include "dmumps_c.h"
 #endif
+#endif
 #include <string>
 #include <iostream>
 #include <algorithm>

config/makefile

@@ -42,7 +42,7 @@ GHV_FLAGS = $(MFEM_CXXFLAGS) $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..)
 SMX = $(if $(MFEM_USE_PUMI:NO=),MFEM_USE_SIMMETRIX)
 SMX_PATH = $(PUMI_DIR)/include/gmi_sim.h
 SMX_FILE = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(SMX_PATH))
-MUMPS = $(MFEM_USE_MUMPS:NO=)
+MUMPS = $(if $(filter YES,$(MFEM_USE_MUMPS) $(MFEM_USE_COMPLEX_MUMPS)),YES,)
 GMV_CXX ?= $(MFEM_CXX)
 GMV = get_mumps_version
 GMV_FLAGS = $(MFEM_CXXFLAGS) $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(MUMPS_OPT))

examples/ex25p.cpp

@@ -200,7 +200,7 @@ int main(int argc, char *argv[])
    args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
                   "--no-superlu", "Use the SuperLU Solver.");
 #endif
-#ifdef MFEM_USE_MUMPS
+#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
    args.AddOption(&mumps_solver, "-mumps", "--mumps-solver", "-no-mumps",
                   "--no-mumps-solver", "Use the MUMPS Solver.");
 #endif
@@ -502,16 +502,24 @@ int main(int argc, char *argv[])
       delete A;
    }
 #endif
-#ifdef MFEM_USE_MUMPS
+#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
    if (!pa && mumps_solver)
    {
-      HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+      ComplexHypreParMatrix *Ahc = Ah.As<ComplexHypreParMatrix>();
+#ifdef MFEM_USE_COMPLEX_MUMPS
+      ComplexMUMPSSolver mumps(MPI_COMM_WORLD);
+      mumps.SetPrintLevel(0);
+      mumps.SetOperator(*Ahc);
+      mumps.Mult(B, X);
+#else
+      HypreParMatrix *A = Ahc->GetSystemMatrix();
       MUMPSSolver mumps(A->GetComm());
       mumps.SetPrintLevel(0);
       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
       mumps.SetOperator(*A);
       mumps.Mult(B, X);
       delete A;
+#endif
    }
 #endif
    // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner

linalg/linalg.hpp

@@ -65,7 +65,7 @@
 #include "hypre_parcsr.hpp"
 #include "hypre.hpp"
 
-#ifdef MFEM_USE_MUMPS
+#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
 #include "mumps.hpp"
 #endif
 

linalg/mumps.cpp

@@ -12,11 +12,9 @@
 #include "../config/config.hpp"
 #include "../general/communication.hpp"
 
-#ifdef MFEM_USE_MUMPS
-#ifdef MFEM_USE_MPI
-
+#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
 #include "mumps.hpp"
-
+#include <unordered_map>
 #include <algorithm>
 
 #if MFEM_MUMPS_VERSION >= 530
@@ -35,9 +33,14 @@
 #define MUMPS_INFO(I) info[(I) -1]
 #define MUMPS_INFOG(I) infog[(I) -1]
 
+#endif // MFEM_USE_MUMPS || MFEM_USE_COMPLEX_MUMPS
+
+
 namespace mfem
 {
 
+#ifdef MFEM_USE_MUMPS
+
 MUMPSSolver::MUMPSSolver(MPI_Comm comm_)
 {
    Init(comm_);
@@ -667,7 +670,760 @@ void MUMPSSolver::RedistributeSol(const int *rmap, const real_t *x,
 }
 #endif
 
-} // namespace mfem
-
-#endif // MFEM_USE_MPI
 #endif // MFEM_USE_MUMPS
+
+
+#ifdef MFEM_USE_COMPLEX_MUMPS
+
+ComplexMUMPSSolver::ComplexMUMPSSolver(MPI_Comm comm_)
+{
+   Init(comm_);
+}
+
+ComplexMUMPSSolver::ComplexMUMPSSolver(const Operator &op)
+{
+   auto APtr = dynamic_cast<const ComplexHypreParMatrix *>(&op);
+   MFEM_VERIFY(APtr,
+               "ComplexMUMPSSolver requires a ComplexHypreParMatrix operator");
+   SetOperator(op);
+}
+
+void ComplexMUMPSSolver::Init(MPI_Comm comm_)
+{
+   comm = comm_;
+   MPI_Comm_size(comm, &numProcs);
+   MPI_Comm_rank(comm, &myid);
+
+   print_level = 2;
+   row_start = 0;
+
+   id = nullptr;
+
+#if MFEM_MUMPS_VERSION >= 530
+   irhs_loc = nullptr;
+   isol_loc = nullptr;
+   rhs_loc  = nullptr;
+   sol_loc  = nullptr;
+#else
+   global_num_rows = 0;
+   recv_counts = nullptr;
+   displs = nullptr;
+   rhs_glob = nullptr;
+   rhs_glob_r = nullptr;
+   rhs_glob_i = nullptr;
+#endif
+}
+
+ComplexMUMPSSolver::~ComplexMUMPSSolver()
+{
+#if MFEM_MUMPS_VERSION >= 530
+   delete [] irhs_loc;
+   delete [] isol_loc;
+   delete [] rhs_loc;
+   delete [] sol_loc;
+#else
+   delete [] recv_counts;
+   delete [] displs;
+   delete [] rhs_glob;
+   delete [] rhs_glob_r;
+   delete [] rhs_glob_i;
+#endif
+
+   if (id)
+   {
+      id->job = -2;
+      mumps_call();
+      delete id;
+      id = nullptr;
+   }
+}
+
+void ComplexMUMPSSolver::SetOperator(const Operator &op)
+{
+   auto APtr = dynamic_cast<const ComplexHypreParMatrix *>(&op);
+   MFEM_VERIFY(APtr,
+               "ComplexMUMPSSolver requires a ComplexHypreParMatrix operator");
+
+   height = op.Height();
+   width  = op.Width();
+
+   const HypreParMatrix *Ar = (APtr->hasRealPart()) ? &APtr->real() : nullptr;
+   const HypreParMatrix *Ai = (APtr->hasImagPart()) ? &APtr->imag() : nullptr;
+
+   MFEM_VERIFY(Ar || Ai, "ComplexMUMPSSolver: both real and imag parts are null.");
+
+   // Pick communicator from the non-null part
+   MPI_Comm op_comm = (Ar ? Ar->GetComm() : Ai->GetComm());
+
+   // Comm setup/check
+   if (comm == MPI_COMM_NULL) { Init(op_comm); }
+   else
+   {
+      int cmp = MPI_UNEQUAL;
+      MPI_Comm_compare(comm, op_comm, &cmp);
+      MFEM_VERIFY(cmp != MPI_UNEQUAL, "MPI Comm mismatch");
+   }
+
+   // HostRead only if non-null
+   if (Ar) { Ar->HostRead(); }
+   if (Ai) { Ai->HostRead(); }
+
+   // hypre parcsr pointers
+   hypre_ParCSRMatrix *parcsr_op_r = nullptr;
+   hypre_ParCSRMatrix *parcsr_op_i = nullptr;
+
+   if (Ar) { parcsr_op_r = (hypre_ParCSRMatrix*) const_cast<HypreParMatrix&>(*Ar); }
+   if (Ai) { parcsr_op_i = (hypre_ParCSRMatrix*) const_cast<HypreParMatrix&>(*Ai); }
+
+   // Merge diag+offd for whichever exists
+   hypre_CSRMatrix *csr_op_r = nullptr;
+   hypre_CSRMatrix *csr_op_i = nullptr;
+
+   if (parcsr_op_r) { csr_op_r = hypre_MergeDiagAndOffd(parcsr_op_r); }
+   if (parcsr_op_i) { csr_op_i = hypre_MergeDiagAndOffd(parcsr_op_i); }
+
+#if MFEM_HYPRE_VERSION >= 21600
+   if (csr_op_r) { hypre_CSRMatrixBigJtoJ(csr_op_r); }
+   if (csr_op_i) { hypre_CSRMatrixBigJtoJ(csr_op_i); }
+#endif
+
+   // Determine local/global sizes and row_start from an existing part
+   const int n_loc = internal::to_int((csr_op_r ? csr_op_r->num_rows :
+                                       csr_op_i->num_rows));
+   row_start = internal::to_int((parcsr_op_r ? parcsr_op_r->first_row_index
+                                 : parcsr_op_i->first_row_index));
+   const int global_n = internal::to_int((parcsr_op_r ?
+                                          parcsr_op_r->global_num_rows
+                                          : parcsr_op_i->global_num_rows));
+
+   // Use nullptr checks
+   const int *Ir = csr_op_r ? csr_op_r->i : nullptr;
+   const int *Jr = csr_op_r ? csr_op_r->j : nullptr;
+   const real_t *Vr = csr_op_r ? (const real_t*)csr_op_r->data : nullptr;
+
+   const int *Ii = csr_op_i ? csr_op_i->i : nullptr;
+   const int *Ji = csr_op_i ? csr_op_i->j : nullptr;
+   const real_t *Vi = csr_op_i ? (const real_t*)csr_op_i->data : nullptr;
+
+   // Build union COO
+   std::vector<int> Icoo, Jcoo;
+   std::vector<mumps_complex_t> Zcoo;
+
+   size_t nnz_r = csr_op_r ? (size_t)csr_op_r->num_nonzeros : 0;
+   size_t nnz_i = csr_op_i ? (size_t)csr_op_i->num_nonzeros : 0;
+   Icoo.reserve(nnz_r + nnz_i);
+   Jcoo.reserve(nnz_r + nnz_i);
+   Zcoo.reserve(nnz_r + nnz_i);
+
+   BuildUnionCOO(n_loc, row_start, Ir, Jr, Vr, Ii, Ji, Vi, Icoo, Jcoo, Zcoo);
+
+   const int nnz = (int)Icoo.size();
+   int *I = new int[nnz];
+   int *J = new int[nnz];
+   mumps_complex_t *A = new mumps_complex_t[nnz];
+
+   std::copy(Icoo.begin(), Icoo.end(), I);
+   std::copy(Jcoo.begin(), Jcoo.end(), J);
+   std::copy(Zcoo.begin(), Zcoo.end(), A);
+
+   // New ComplexMUMPS object or reuse an existing one
+   if (!id || !reorder_reuse)
+   {
+      if (id)
+      {
+         id->job = -2;
+         mumps_call();
+         delete id;
+         id = nullptr;
+      }
+
+#ifdef MFEM_USE_SINGLE
+      id = new CMUMPS_STRUC_C();
+#else
+      id = new ZMUMPS_STRUC_C();
+#endif
+
+      id->sym = 0; // general complex
+      id->par = 1;
+      id->comm_fortran = (MUMPS_INT)MPI_Comm_c2f(comm);
+
+      // Init
+      id->job = -1;
+      mumps_call();
+
+      // Set parameters
+      SetParameters();
+
+      // Attach matrix
+      id->n       = global_n;
+      id->nnz_loc = nnz;
+      id->irn_loc = I;
+      id->jcn_loc = J;
+      id->a_loc   = A;
+
+      // Analysis (ordering + symbolic)
+      id->job = 1;
+      mumps_call();
+   }
+   else
+   {
+      // Reuse symbolic factorization / ordering
+      MFEM_VERIFY(id->n == global_n,
+                  "ReorderingReuse requires same global size (id->n mismatch)");
+
+      // Update matrix pointers (pattern is assumed compatible)
+      id->nnz_loc = nnz;
+      id->irn_loc = I;
+      id->jcn_loc = J;
+      id->a_loc   = A;
+   }
+
+   // Factorization
+   id->job = 2;
+   {
+      const int mem_relax_lim = 200;
+      while (true)
+      {
+         mumps_call();
+         if (id->MUMPS_INFOG(1) < 0)
+         {
+            if (id->MUMPS_INFOG(1) == -8 || id->MUMPS_INFOG(1) == -9)
+            {
+               id->MUMPS_ICNTL(14) += 20;
+               MFEM_VERIFY(id->MUMPS_ICNTL(14) <= mem_relax_lim,
+                           "Memory relaxation limit reached for ComplexMUMPSSolver factorization");
+               if (myid == 0 && print_level > 0)
+               {
+                  out << "Re-running ComplexMUMPSSolver factorization with memory relaxation "
+                      << id->MUMPS_ICNTL(14) << '\n';
+               }
+            }
+            else
+            {
+               MFEM_ABORT("Error during ComplexMUMPSSolver numerical factorization");
+            }
+         }
+         else { break; }
+      }
+   }
+
+   // Done with input storage
+   if (csr_op_r) { hypre_CSRMatrixDestroy(csr_op_r);}
+   if (csr_op_i) { hypre_CSRMatrixDestroy(csr_op_i);}
+   delete [] I;
+   delete [] J;
+   delete [] A;
+
+   // Post-factorization RHS/SOL setup
+   id->nrhs = -1;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Distributed RHS/SOL sizes
+   id->nloc_rhs = n_loc;
+   id->lrhs_loc = n_loc;
+   id->lsol_loc = id->MUMPS_INFO(23);
+
+   delete [] irhs_loc;
+   irhs_loc = new int[id->lrhs_loc];
+   for (int i = 0; i < n_loc; i++)
+   {
+      irhs_loc[i] = row_start + i + 1;
+   }
+   id->irhs_loc = irhs_loc;
+
+   delete [] isol_loc;
+   isol_loc = new int[id->lsol_loc];
+   id->isol_loc = isol_loc;
+
+   row_starts.SetSize(numProcs);
+   MPI_Allgather(&row_start, 1, MPI_INT, row_starts, 1, MPI_INT, comm);
+
+   // Reset cached buffers
+   delete [] rhs_loc; rhs_loc = nullptr;
+   delete [] sol_loc; sol_loc = nullptr;
+   rhs1_buf.clear();
+
+#else
+   // Centralized RHS/SOL on root
+   id->lrhs = id->n;
+
+   global_num_rows = id->n;
+
+   if (myid == 0)
+   {
+      delete [] recv_counts;
+      delete [] displs;
+      recv_counts = new int[numProcs];
+      displs = new int[numProcs];
+
+      delete [] rhs_glob;   rhs_glob = nullptr;
+      delete [] rhs_glob_r; rhs_glob_r = nullptr;
+      delete [] rhs_glob_i; rhs_glob_i = nullptr;
+   }
+
+   MPI_Gather(&n_loc, 1, MPI_INT, recv_counts, 1, MPI_INT, 0, comm);
+
+   if (myid == 0)
+   {
+      displs[0] = 0;
+      int s = 0;
+      for (int k = 0; k < numProcs - 1; k++)
+      {
+         s += recv_counts[k];
+         displs[k+1] = s;
+      }
+   }
+#endif
+}
+
+void ComplexMUMPSSolver::InitRhsSol(int nrhs) const
+{
+#if MFEM_MUMPS_VERSION >= 530
+
+   MFEM_VERIFY(id, "InitRhsSol called before SetOperator");
+
+   if (id->nrhs != nrhs)
+   {
+      delete [] rhs_loc;
+      delete [] sol_loc;
+
+      rhs_loc = new mumps_complex_t[(size_t)nrhs * (size_t)id->lrhs_loc];
+      sol_loc = new mumps_complex_t[(size_t)nrhs * (size_t)id->lsol_loc];
+
+      id->rhs_loc = rhs_loc;
+      id->sol_loc = sol_loc;
+   }
+   id->nrhs = nrhs;
+
+#else
+   MFEM_VERIFY(id, "InitRhsSol called before SetOperator");
+
+   id->nrhs = nrhs;
+   id->lrhs = id->n;
+
+   if (myid == 0)
+   {
+      const size_t N = (size_t)nrhs * (size_t)global_num_rows;
+
+      delete [] rhs_glob;
+      delete [] rhs_glob_r;
+      delete [] rhs_glob_i;
+
+      rhs_glob   = new mumps_complex_t[N];
+      rhs_glob_r = new real_t[N];
+      rhs_glob_i = new real_t[N];
+
+      id->rhs = rhs_glob;
+   }
+#endif
+}
+
+void ComplexMUMPSSolver::Mult(const Vector &x, Vector &y) const
+{
+   Array<const Vector *> X(1);
+   Array<Vector *> Y(1);
+   X[0] = &x;
+   Y[0] = &y;
+   ArrayMult(X, Y);
+}
+
+void ComplexMUMPSSolver::ArrayMult(const Array<const Vector *> &X,
+                                   Array<Vector *> &Y) const
+{
+   MFEM_ASSERT(X.Size() == Y.Size(),
+               "Number of columns mismatch in ComplexMUMPSSolver::Mult!");
+   MFEM_VERIFY(id, "ComplexMUMPSSolver::ArrayMult called before SetOperator");
+
+   InitRhsSol(X.Size());
+
+#if MFEM_MUMPS_VERSION >= 530
+   MFEM_VERIFY(irhs_loc && isol_loc, "RHS/SOL maps not initialized");
+   MFEM_VERIFY(rhs_loc && sol_loc, "RHS/SOL buffers not initialized");
+   const int n_loc = id->lrhs_loc;
+   const int nrhs  = id->nrhs;
+
+   // Pack all RHS
+   int xisign = (conv == ComplexOperator::BLOCK_SYMMETRIC) ? -1 : 1;
+   for (int i = 0; i < nrhs; i++)
+   {
+      MFEM_ASSERT(X[i], "Missing Vector in Mult!");
+      X[i]->HostRead();
+      MFEM_VERIFY(X[i]->Size() == 2*n_loc, "RHS size mismatch");
+
+      const real_t *xdata = X[i]->GetData();
+      const real_t *xr = xdata;
+      const real_t *xi = xdata + n_loc;
+
+      mumps_complex_t *dst = rhs_loc + i * n_loc;
+      for (int j = 0; j < n_loc; j++)
+      {
+         dst[j].r = xr[j];
+         dst[j].i = xisign * xi[j];
+      }
+   }
+
+   id->rhs_loc  = rhs_loc;
+   id->sol_loc  = sol_loc;
+   id->irhs_loc = irhs_loc;
+   id->isol_loc = isol_loc;
+
+   // MUMPS solve
+   id->job = 3;
+   mumps_call();
+
+   const int lsol = id->lsol_loc;
+
+   // Redistribute each solution column into Y
+   for (int i = 0; i < nrhs; i++)
+   {
+      MFEM_ASSERT(Y[i], "Missing output Vector in Mult!");
+      Y[i]->HostWrite();
+      MFEM_VERIFY(Y[i]->Size() == 2*n_loc, "Output size mismatch");
+
+      const mumps_complex_t *xcol = sol_loc + i * lsol;
+      RedistributeSol(isol_loc, xcol, Y[i]->GetData(), n_loc, lsol);
+   }
+
+#else // MFEM_MUMPS_VERSION < 530
+
+   const int nrhs = id->nrhs;
+
+   MFEM_VERIFY(X.Size() > 0 && X[0], "Missing RHS");
+   const int n_loc = X[0]->Size()/2;
+
+   for (int i = 0; i < nrhs; i++)
+   {
+      MFEM_ASSERT(X[i], "Missing Vector in Mult!");
+      X[i]->HostRead();
+      MFEM_VERIFY(X[i]->Size() == 2*n_loc, "RHS size mismatch");
+   }
+
+   // Gather each RHS column (real+imag separately) into root staging
+   for (int i = 0; i < nrhs; i++)
+   {
+      const real_t *xdata = X[i]->GetData();
+
+      MPI_Gatherv(xdata, n_loc, MPITypeMap<real_t>::mpi_type,
+                  rhs_glob_r + i * global_num_rows,
+                  recv_counts, displs, MPITypeMap<real_t>::mpi_type,
+                  0, comm);
+
+      MPI_Gatherv(xdata + n_loc, n_loc, MPITypeMap<real_t>::mpi_type,
+                  rhs_glob_i + i * global_num_rows,
+                  recv_counts, displs, MPITypeMap<real_t>::mpi_type,
+                  0, comm);
+   }
+
+   // Pack into MUMPS complex RHS on root: id->rhs is in-place
+   if (myid == 0)
+   {
+      for (int i = 0; i < nrhs; i++)
+      {
+         mumps_complex_t *dst = rhs_glob + i * global_num_rows;
+         const real_t *rr = rhs_glob_r + i * global_num_rows;
+         const real_t *ri = rhs_glob_i + i * global_num_rows;
+
+         for (int j = 0; j < global_num_rows; j++)
+         {
+            dst[j].r = rr[j];
+            dst[j].i = ri[j];
+         }
+      }
+      id->rhs = rhs_glob;
+   }
+
+   // Solve
+   id->job = 3;
+   mumps_call();
+
+   // Unpack to real/imag
+   if (myid == 0)
+   {
+      for (int i = 0; i < nrhs; i++)
+      {
+         const mumps_complex_t *src = rhs_glob + i * global_num_rows;
+         real_t *rr = rhs_glob_r + i * global_num_rows;
+         real_t *ri = rhs_glob_i + i * global_num_rows;
+
+         for (int j = 0; j < global_num_rows; j++)
+         {
+            rr[j] = src[j].r;
+            ri[j] = src[j].i;
+         }
+      }
+   }
+
+   // Scatter each RHS solution
+   for (int i = 0; i < nrhs; i++)
+   {
+      MFEM_ASSERT(Y[i], "Missing Vector in Mult!");
+      Y[i]->HostWrite();
+      MFEM_VERIFY(Y[i]->Size() == 2*n_loc, "Output size mismatch");
+
+      real_t *ydata = Y[i]->GetData();
+
+      MPI_Scatterv(rhs_glob_r + i * global_num_rows,
+                   recv_counts, displs, MPITypeMap<real_t>::mpi_type,
+                   ydata, n_loc, MPITypeMap<real_t>::mpi_type,
+                   0, comm);
+
+      MPI_Scatterv(rhs_glob_i + i * global_num_rows,
+                   recv_counts, displs, MPITypeMap<real_t>::mpi_type,
+                   ydata + n_loc, n_loc, MPITypeMap<real_t>::mpi_type,
+                   0, comm);
+   }
+
+#endif
+}
+
+void ComplexMUMPSSolver::MultTranspose(const Vector &x, Vector &y) const
+{
+   MFEM_VERIFY(id, "MultTranspose called before SetOperator");
+
+   // Transpose solve
+   id->MUMPS_ICNTL(9) = 0;
+   Mult(x, y);
+   id->MUMPS_ICNTL(9) = 1;
+}
+
+void ComplexMUMPSSolver::ArrayMultTranspose(const Array<const Vector *> &X,
+                                            Array<Vector *> &Y) const
+{
+   MFEM_VERIFY(id, "ArrayMultTranspose called before SetOperator");
+
+   // Transpose solve
+   id->MUMPS_ICNTL(9) = 0;
+   ArrayMult(X, Y);
+   id->MUMPS_ICNTL(9) = 1;
+}
+
+void ComplexMUMPSSolver::SetParameters()
+{
+   // Output stream for error messages
+   id->MUMPS_ICNTL(1) = 6;
+   // Output stream for diagnostic printing local to each proc
+   id->MUMPS_ICNTL(2) = 0;
+   // Output stream for global info
+   id->MUMPS_ICNTL(3) = 6;
+   // Level of error printing
+   id->MUMPS_ICNTL(4) = print_level;
+
+   // Input matrix format (assembled)
+   id->MUMPS_ICNTL(5) = 0;
+   // Use A or A^T
+   id->MUMPS_ICNTL(9) = 1;
+   // Iterative refinement (disabled)
+   id->MUMPS_ICNTL(10) = 0;
+   // Error analysis-statistics (disabled)
+   id->MUMPS_ICNTL(11) = 0;
+   // Use of ScaLAPACK (disabled)
+   id->MUMPS_ICNTL(13) = 0;
+   // Workspace relaxation (% increase)
+   id->MUMPS_ICNTL(14) = 20;
+   // OpenMP threads (default)
+   id->MUMPS_ICNTL(16) = 0;
+   // Matrix input format (distributed)
+   id->MUMPS_ICNTL(18) = 3;
+   // Schur complement (none)
+   id->MUMPS_ICNTL(19) = 0;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Distributed RHS
+   id->MUMPS_ICNTL(20) = 10;
+   // Distributed Sol
+   id->MUMPS_ICNTL(21) = 1;
+#else
+   // Centralized RHS
+   id->MUMPS_ICNTL(20) = 0;
+   // Centralized Sol
+   id->MUMPS_ICNTL(21) = 0;
+#endif
+
+   // Out-of-core (disabled)
+   id->MUMPS_ICNTL(22) = 0;
+   // Max size of working memory (default)
+   id->MUMPS_ICNTL(23) = 0;
+
+   switch (reorder_method)
+   {
+      case ReorderingStrategy::AUTOMATIC:
+         id->MUMPS_ICNTL(28) = 0;
+         id->MUMPS_ICNTL(7) = 7;
+         id->MUMPS_ICNTL(29) = 0;
+         break;
+      case ReorderingStrategy::AMD:
+         id->MUMPS_ICNTL(28) = 1;
+         id->MUMPS_ICNTL(7) = 0;
+         break;
+      case ReorderingStrategy::AMF:
+         id->MUMPS_ICNTL(28) = 1;
+         id->MUMPS_ICNTL(7) = 2;
+         break;
+      case ReorderingStrategy::PORD:
+         id->MUMPS_ICNTL(28) = 1;
+         id->MUMPS_ICNTL(7) = 4;
+         break;
+      case ReorderingStrategy::METIS:
+         id->MUMPS_ICNTL(28) = 1;
+         id->MUMPS_ICNTL(7) = 5;
+         break;
+      case ReorderingStrategy::PARMETIS:
+         id->MUMPS_ICNTL(28) = 2;
+         id->MUMPS_ICNTL(29) = 2;
+         break;
+      case ReorderingStrategy::SCOTCH:
+         id->MUMPS_ICNTL(28) = 1;
+         id->MUMPS_ICNTL(7) = 3;
+         break;
+      case ReorderingStrategy::PTSCOTCH:
+         id->MUMPS_ICNTL(28) = 2;
+         id->MUMPS_ICNTL(29) = 1;
+         break;
+      default:
+         break; // This should be unreachable
+   }
+
+}
+
+void ComplexMUMPSSolver::BuildUnionCOO(const int n_loc,
+                                       const int row_start_,
+                                       const int *Ir, const int *Jr, const real_t *Vr,
+                                       const int *Ii, const int *Ji, const real_t *Vi,
+                                       std::vector<int> &Icoo,
+                                       std::vector<int> &Jcoo,
+                                       std::vector<mumps_complex_t> &Zcoo) const
+{
+   for (int r = 0; r < n_loc; ++r)
+   {
+      std::unordered_map<int, std::pair<real_t, real_t>> row;
+
+      const int rr0 = Ir ? Ir[r]   : 0;
+      const int rr1 = Ir ? Ir[r+1] : 0;
+      const int ii0 = Ii ? Ii[r]   : 0;
+      const int ii1 = Ii ? Ii[r+1] : 0;
+
+      row.reserve((rr1 - rr0) + (ii1 - ii0));
+
+      if (Ir)
+      {
+         for (int p = rr0; p < rr1; ++p) { row[Jr[p]].first += Vr[p]; }
+      }
+      if (Ii)
+      {
+         for (int p = ii0; p < ii1; ++p) { row[Ji[p]].second += Vi[p]; }
+      }
+
+      for (const auto &kv : row)
+      {
+         Icoo.push_back(row_start_ + r + 1);
+         Jcoo.push_back(kv.first + 1);
+         Zcoo.push_back(mumps_complex_t{kv.second.first, kv.second.second});
+      }
+   }
+}
+
+#if MFEM_MUMPS_VERSION >= 530
+int ComplexMUMPSSolver::GetRowRank(int i, const Array<int> &row_starts_) const
+{
+   if (row_starts_.Size() == 1) { return 0; }
+   auto up = std::upper_bound(row_starts_.begin(), row_starts_.end(), i);
+   return (int)std::distance(row_starts_.begin(), up) - 1;
+}
+
+void ComplexMUMPSSolver::RedistributeSol(const int *row_map,
+                                         const mumps_complex_t *x,
+                                         real_t *y_ri,
+                                         int n_loc,
+                                         int lsol_loc) const
+{
+   int *send_count = new int[numProcs]();
+   for (int i = 0; i < lsol_loc; i++)
+   {
+      const int j = row_map[i] - 1;
+      const int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank) { continue; }
+      send_count[row_rank]++;
+   }
+
+   int *recv_count = new int[numProcs];
+   MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
+
+   int *send_displ = new int[numProcs]; send_displ[0] = 0;
+   int *recv_displ = new int[numProcs]; recv_displ[0] = 0;
+
+   int sbuff_size = send_count[numProcs-1];
+   int rbuff_size = recv_count[numProcs-1];
+   for (int k = 0; k < numProcs - 1; k++)
+   {
+      send_displ[k+1] = send_displ[k] + send_count[k];
+      recv_displ[k+1] = recv_displ[k] + recv_count[k];
+      sbuff_size += send_count[k];
+      rbuff_size += recv_count[k];
+   }
+
+   int   *sendbuf_index = new int[sbuff_size];
+   real_t *sendbuf_r    = new real_t[sbuff_size];
+   real_t *sendbuf_i    = new real_t[sbuff_size];
+   int   *soffs         = new int[numProcs]();
+
+   for (int i = 0; i < lsol_loc; i++)
+   {
+      const int j = row_map[i] - 1;
+      const int row_rank = GetRowRank(j, row_starts);
+
+      const real_t xr = (real_t)x[i].r;
+      const real_t xi = (real_t)x[i].i;
+
+      if (myid == row_rank)
+      {
+         const int local_index = j - row_start;
+         y_ri[local_index]        = xr;
+         y_ri[local_index+n_loc]  = xi;
+      }
+      else
+      {
+         const int k = send_displ[row_rank] + soffs[row_rank];
+         sendbuf_index[k] = j;
+         sendbuf_r[k] = xr;
+         sendbuf_i[k] = xi;
+         soffs[row_rank]++;
+      }
+   }
+
+   int    *recvbuf_index = new int[rbuff_size];
+   real_t *recvbuf_r     = new real_t[rbuff_size];
+   real_t *recvbuf_i     = new real_t[rbuff_size];
+
+   MPI_Alltoallv(sendbuf_index, send_count, send_displ, MPI_INT,
+                 recvbuf_index, recv_count, recv_displ, MPI_INT, comm);
+
+   MPI_Alltoallv(sendbuf_r, send_count, send_displ, MPITypeMap<real_t>::mpi_type,
+                 recvbuf_r, recv_count, recv_displ, MPITypeMap<real_t>::mpi_type, comm);
+
+   MPI_Alltoallv(sendbuf_i, send_count, send_displ, MPITypeMap<real_t>::mpi_type,
+                 recvbuf_i, recv_count, recv_displ, MPITypeMap<real_t>::mpi_type, comm);
+
+   for (int i = 0; i < rbuff_size; i++)
+   {
+      const int local_index = recvbuf_index[i] - row_start;
+      y_ri[local_index]       = recvbuf_r[i];
+      y_ri[local_index+n_loc] = recvbuf_i[i];
+   }
+
+   delete [] recvbuf_i;
+   delete [] recvbuf_r;
+   delete [] recvbuf_index;
+   delete [] soffs;
+   delete [] sendbuf_i;
+   delete [] sendbuf_r;
+   delete [] sendbuf_index;
+   delete [] recv_displ;
+   delete [] send_displ;
+   delete [] recv_count;
+   delete [] send_count;
+}
+#endif // MFEM_MUMPS_VERSION >= 530
+
+#endif // MFEM_USE_COMPLEX_MUMPS
+
+} // namespace mfem

linalg/mumps.hpp

@@ -14,22 +14,34 @@
 
 #include "../config/config.hpp"
 
-#ifdef MFEM_USE_MUMPS
-#ifdef MFEM_USE_MPI
-
+#if defined(MFEM_USE_MUMPS) || defined(MFEM_USE_COMPLEX_MUMPS)
 #include "operator.hpp"
 #include "hypre.hpp"
 #include <mpi.h>
+#endif
 
+#ifdef MFEM_USE_MUMPS
 #ifdef MFEM_USE_SINGLE
 #include "smumps_c.h"
 #else
 #include "dmumps_c.h"
 #endif
+#endif
+
+#ifdef MFEM_USE_COMPLEX_MUMPS
+#include "complex_operator.hpp"
+#ifdef MFEM_USE_SINGLE
+#include "cmumps_c.h"
+#else
+#include "zmumps_c.h"
+#endif
+#include <vector>
+#endif
 
 namespace mfem
 {
 
+#ifdef MFEM_USE_MUMPS
 /**
  * @brief MUMPS: A Parallel Sparse Direct Solver
  *
@@ -52,7 +64,7 @@ public:
    /// Specify the reordering strategy for the MUMPS solver
    enum ReorderingStrategy
    {
-      /// Let MUMPS automatically decide the reording strategy
+      /// Let MUMPS automatically decide the reordering strategy
       AUTOMATIC = 0,
       /// Approximate Minimum Degree with auto quasi-dense row detection is used
       AMD,
@@ -98,7 +110,7 @@ public:
    void Mult(const Vector &x, Vector &y) const;
 
    /**
-    * @brief Solve $ Y_i = Op^{-T} X_i $
+    * @brief Solve $ Y_i = Op^{-1} X_i $
     *
     * @param X Array of RHS vectors
     * @param Y Array of Solution vectors
@@ -129,8 +141,8 @@ public:
     * - 0:  No output printed
     * - 1:  Only errors printed
     * - 2:  Errors, warnings, and main stats printed
-    * - 3:  Errors, warning, main stats, and terse diagnostics printed
-    * - 4:  Errors, warning, main stats, diagnostics, and input/output printed
+    * - 3:  Errors, warnings, main stats, and terse diagnostics printed
+    * - 4:  Errors, warnings, main stats, diagnostics, and input/output printed
     *
     * @param print_lvl Print level, default is 2
     *
@@ -258,8 +270,247 @@ private:
 #endif
 }; // mfem::MUMPSSolver class
 
+#endif // MFEM_USE_MUMPS
+
+#ifdef MFEM_USE_COMPLEX_MUMPS
+/**
+ * @brief Complex MUMPS: Parallel sparse direct solver for ComplexHypreParMatrix
+ *
+ * Notes:
+ *  - Expects Operator to be a ComplexHypreParMatrix.
+ *  - Complex vectors are assumed packed as [Re; Im] in a real Vector.
+ *  - SetOperator(): analysis + factorization
+ *  - Mult()       : solve
+ */
+class ComplexMUMPSSolver : public Solver
+{
+public:
+   /// Specify the reordering strategy
+   enum ReorderingStrategy
+   {
+      /// Let MUMPS automatically decide the reordering strategy
+      AUTOMATIC = 0,
+      /// Approximate Minimum Degree with auto quasi-dense row detection is used
+      AMD,
+      /// Approximate Minimum Fill method will be used
+      AMF,
+      /// The PORD library will be used
+      PORD,
+      /// The METIS library will be used
+      METIS,
+      /// The ParMETIS library will be used
+      PARMETIS,
+      /// The Scotch library will be used
+      SCOTCH,
+      /// The PTScotch library will be used
+      PTSCOTCH
+   };
+
+   /**
+    * @brief Constructor with MPI_Comm parameter.
+    */
+   ComplexMUMPSSolver(MPI_Comm comm_);
+   /**
+    * @brief Constructor with a ComplexHypreParMatrix Operator.
+    */
+   ComplexMUMPSSolver(const Operator &op);
+
+   /**
+    * @brief Set the Operator and perform factorization
+    *
+    * @a op needs to be of type ComplexHypreParMatrix.
+    *
+    * @param op Operator used in factorization and solve
+    */
+   void SetOperator(const Operator &op);
+
+   /**
+    * @brief Solve $ y = Op^{-1} x $
+    *
+    * @param x RHS vector
+    * @param y Solution vector
+    */
+   void Mult(const Vector &x, Vector &y) const;
+   /**
+    * @brief Solve $ Y_i = Op^{-1} X_i $
+    *
+    * @param X Array of RHS vectors
+    * @param Y Array of Solution vectors
+    */
+   void ArrayMult(const Array<const Vector *> &X, Array<Vector *> &Y) const;
+   /**
+    * @brief Transpose Solve $ y = Op^{-T} x $
+    * @note This is not a Hermitian/conjugate-transpose solve.
+    *
+    * @param x RHS vector
+    * @param y Solution vector
+    */
+   void MultTranspose(const Vector &x, Vector &y) const;
+
+   /**
+    * @brief Transpose Solve $ Y_i = Op^{-T} X_i $
+    * @note This is not a Hermitian/conjugate-transpose solve.
+    *
+    * @param X Array of RHS vectors
+    * @param Y Array of Solution vectors
+    */
+   void ArrayMultTranspose(const Array<const Vector *> &X,
+                           Array<Vector *> &Y) const;
+
+   /**
+    * @brief Set the error print level for MUMPS
+    *
+    * Supported values are:
+    * - 0:  No output printed
+    * - 1:  Only errors printed
+    * - 2:  Errors, warnings, and main stats printed
+    * - 3:  Errors, warnings, main stats, and terse diagnostics printed
+    * - 4:  Errors, warnings, main stats, diagnostics, and input/output printed
+    *
+    * @param print_lvl Print level, default is 2
+    *
+    * @note This method has to be called before SetOperator
+    */
+   void SetPrintLevel(int print_lvl) { print_level = print_lvl;}
+
+   /**
+    * @brief Set the reordering strategy
+    *
+    * Supported reorderings are: ComplexMUMPSSolver::AUTOMATIC,
+    * ComplexMUMPSSolver::AMD, ComplexMUMPSSolver::AMF,
+    * ComplexMUMPSSolver::PORD, ComplexMUMPSSolver::METIS,
+    * ComplexMUMPSSolver::PARMETIS, ComplexMUMPSSolver::SCOTCH,
+    * and ComplexMUMPSSolver::PTSCOTCH
+    *
+    * @param method Reordering method
+    *
+    * @note This method has to be called before SetOperator
+    */
+   void SetReorderingStrategy(ReorderingStrategy method) { reorder_method = method; }
+
+   /**
+    * @brief Set the flag controlling reuse of the symbolic factorization
+    * for multiple operators
+    *
+    * @param reuse Flag to reuse symbolic factorization
+    *
+    * @note This method has to be called before repeated calls to SetOperator
+    */
+   void SetReorderingReuse(bool reuse) { reorder_reuse = reuse; }
+
+   ~ComplexMUMPSSolver();
+
+private:
+   // MPI communicator
+   MPI_Comm comm = MPI_COMM_NULL;
+
+   // Number of procs
+   int numProcs;
+
+   // MPI rank
+   int myid;
+
+   // Parameter controlling the printing level
+   int print_level = 0;
+
+   // Parameter controlling the reordering strategy
+   ReorderingStrategy reorder_method = ReorderingStrategy::AUTOMATIC;
+
+   // Parameter controlling whether or not to reuse the symbolic factorization
+   // for multiple calls to SetOperator
+   bool reorder_reuse = false;
+
+   // Local row offsets
+   int row_start;
+
+   // ComplexMUMPS object
+#ifdef MFEM_USE_SINGLE
+   CMUMPS_STRUC_C *id = nullptr;
+   using mumps_complex_t = mumps_complex;
+#else
+   ZMUMPS_STRUC_C *id = nullptr;
+   using mumps_complex_t = mumps_double_complex;
+#endif
+
+   /// Method for initialization
+   void Init(MPI_Comm comm_);
+
+   /// Method for setting ComplexMUMPS internal parameters
+   void SetParameters();
+
+   /// Method for configuring storage for distributed/centralized
+   /// RHS and solution
+   void InitRhsSol(int nrhs) const;
+
+   /// Method for calling the single/double ComplexMUMPS solver
+   inline void mumps_call() const
+   {
+#ifdef MFEM_USE_SINGLE
+      cmumps_c(id);
+#else
+      zmumps_c(id);
+#endif
+   }
+
+   /// Method for building the COO format of the combined complex operator
+   /// from the real and imaginary parts. This is particularly useful when
+   /// real and imaginary parts have different sparsity patterns.
+   void BuildUnionCOO(const int n_loc,
+                      const int row_start,
+                      const int *Ir, const int *Jr, const real_t *Vr,
+                      const int *Ii, const int *Ji, const real_t *Vi,
+                      std::vector<int> &Icoo,
+                      std::vector<int> &Jcoo,
+                      std::vector<mumps_complex_t> &Zcoo) const;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Row offsets on all procs
+   Array<int> row_starts;
+
+   // Local RHS row indices
+   int *irhs_loc = nullptr;
+
+   // Local solution row map returned by MUMPS
+   int *isol_loc = nullptr;
+
+   // Cached buffers
+   mutable mumps_complex_t *rhs_loc = nullptr;
+   mutable mumps_complex_t *sol_loc = nullptr;
+
+   // RHS buffers
+   mutable std::vector<mumps_complex_t> rhs1_buf;
+
+   // These two methods are needed to distribute the local solution
+   // vectors returned by MUMPS to the original MFEM parallel partition
+   int GetRowRank(int i, const Array<int> &row_starts_) const;
+
+   void RedistributeSol(const int *row_map,
+                        const mumps_complex_t *x,
+                        real_t *y_ri,
+                        int n_loc,
+                        int lsol_loc) const;
+
+   ComplexOperator::Convention conv = ComplexOperator::Convention::HERMITIAN;
+
+#else
+   // Root-gather path
+   int global_num_rows;
+
+   // Arrays needed for MPI_Gatherv and MPI_Scatterv
+   int *recv_counts = nullptr;
+   int *displs = nullptr;
+
+   // Complex RHS/solution on root
+   mutable mumps_complex_t *rhs_glob = nullptr;
+
+   // Cached real/imag staging on root
+   mutable real_t *rhs_glob_r = nullptr;
+   mutable real_t *rhs_glob_i = nullptr;
+#endif
+};
+
+#endif // MFEM_USE_COMPLEX_MUMPS
+
 } // namespace mfem
 
-#endif // MFEM_USE_MPI
-#endif // MFEM_USE_MUMPS
 #endif // MFEM_MUMPS

makefile

@@ -229,7 +229,7 @@ CXXFLAGS ?= $(OPTIM_FLAGS)
 # MPI configuration
 ifneq ($(MFEM_USE_MPI),YES)
    HOST_CXX = $(CXX)
-   PKGS_NEED_MPI = SUPERLU MUMPS STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
+   PKGS_NEED_MPI = SUPERLU MUMPS COMPLEX_MUMPS STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
    $(foreach mpidep,$(PKGS_NEED_MPI),$(if $(MFEM_USE_$(mpidep):NO=),\
      $(warning *** [MPI is OFF] setting MFEM_USE_$(mpidep) = NO)\
      $(eval override MFEM_USE_$(mpidep)=NO),))
@@ -304,7 +304,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS COMPLEX_MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
  SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
  TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS
@@ -374,7 +374,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
- MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
+ MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_COMPLEX_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
  MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
  MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
  MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB
@@ -412,7 +412,7 @@ MFEM_INSTALL_DIR = $(abspath $(MFEM_PREFIX))
 # If we have 'config' target, export variables used by config/makefile
 ifneq (,$(filter config,$(MAKECMDGOALS)))
    export $(MFEM_DEFINES) MFEM_DEFINES $(MFEM_CONFIG_VARS) MFEM_CONFIG_VARS
-   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT GSLIB_OPT CUDSS_OPT
+   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT COMPLEX_MUMPS_OPT GSLIB_OPT CUDSS_OPT
 endif
 
 # If we have 'install' target, export variables used by config/makefile
@@ -742,6 +742,7 @@ status info:
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
+	$(info MFEM_USE_COMPLEX_MUMPS = $(MFEM_USE_COMPLEX_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))
 	$(info MFEM_USE_CUDSS         = $(MFEM_USE_CUDSS))
 	$(info MFEM_USE_GINKGO        = $(MFEM_USE_GINKGO))

tests/unit/linalg/test_direct_solvers.cpp

@@ -188,8 +188,6 @@ TEST_CASE("Serial Direct Solvers", "[GPU]")
 
 TEST_CASE("Parallel Direct Solvers", "[Parallel], [GPU]")
 {
-   int rank;
-   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    const int ne = 4;
    for (int dim = 1; dim < 4; ++dim)
    {
@@ -357,3 +355,113 @@ TEST_CASE("Parallel Direct Solvers", "[Parallel], [GPU]")
 }
 
 #endif
+
+
+#ifdef MFEM_USE_COMPLEX_MUMPS
+
+TEST_CASE("ComplexMUMPS Solver", "[Parallel], [GPU]")
+{
+   const int ne = 4;
+   const int order = 3;
+   auto dim = GENERATE(1, 2, 3);
+   auto conv = GENERATE(ComplexOperator::HERMITIAN,
+                        ComplexOperator::BLOCK_SYMMETRIC);
+   CAPTURE(dim, conv);
+
+   Mesh mesh;
+   if (dim == 1)
+   {
+      mesh = Mesh::MakeCartesian1D(ne, 1.0);
+   }
+   else if (dim == 2)
+   {
+      mesh = Mesh::MakeCartesian2D(
+                ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
+   }
+   else
+   {
+      mesh = Mesh::MakeCartesian3D(
+                ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+   Array<int> ess_tdof_list, ess_bdr;
+   if (pmesh.bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient negone(-1.0);
+   ConstantCoefficient two(2.0);
+
+   ComplexLinearForm b(&fespace, conv);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one), new DomainLFIntegrator(two));
+   b.Assemble();
+
+   ParSesquilinearForm a_r(&fespace, conv);
+   a_r.AddDomainIntegrator(new DiffusionIntegrator(one), nullptr);
+   a_r.AddDomainIntegrator(new MassIntegrator(one), nullptr);
+   a_r.Assemble();
+
+   ParSesquilinearForm a_i(&fespace, conv);
+   a_i.AddDomainIntegrator(nullptr, new DiffusionIntegrator(one));
+   a_i.AddDomainIntegrator(nullptr, new MassIntegrator(one));
+   a_i.Assemble();
+
+   ParSesquilinearForm a_c(&fespace, conv);
+   a_c.AddDomainIntegrator(new DiffusionIntegrator(one), new MassIntegrator(two));
+   a_c.AddDomainIntegrator(new MassIntegrator(negone),nullptr);
+   a_c.Assemble();
+
+   ParComplexGridFunction x_c(&fespace);
+   ParComplexGridFunction x_r(&fespace);
+   ParComplexGridFunction x_i(&fespace);
+   x_c = 0.0; x_r = 0.0; x_i = 0.0;
+
+   OperatorPtr Ac, Ar, Ai;
+   Vector Bc, Br, Bi, Xr, Xc, Xi;
+   a_c.FormLinearSystem(ess_tdof_list, x_c, b, Ac, Xc, Bc);
+   a_r.FormLinearSystem(ess_tdof_list, x_r, b, Ar, Xr, Br);
+   a_i.FormLinearSystem(ess_tdof_list, x_i, b, Ai, Xi, Bi);
+
+   ComplexHypreParMatrix *Ahc = Ac.As<ComplexHypreParMatrix>();
+   ComplexHypreParMatrix *Ahr = Ar.As<ComplexHypreParMatrix>();
+   ComplexHypreParMatrix *Ahi = Ai.As<ComplexHypreParMatrix>();
+
+   ComplexMUMPSSolver cmumps(MPI_COMM_WORLD);
+   cmumps.SetPrintLevel(0);
+
+   cmumps.SetOperator(*Ahc);
+   cmumps.Mult(Bc, Xc);
+
+   cmumps.SetOperator(*Ahr);
+   cmumps.Mult(Br, Xr);
+
+   cmumps.SetOperator(*Ahi);
+   cmumps.Mult(Bi, Xi);
+
+   Vector Yc(Xc.Size()), Yr(Xr.Size()), Yi(Xi.Size());
+   Ahc->Mult(Xc, Yc); Ahr->Mult(Xr, Yr); Ahi->Mult(Xi, Yi);
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      int n = Bc.Size()/2;
+      real_t * dataBc = Bc.GetData(); Vector Bcimag(&dataBc[n],n);
+      real_t * dataBr = Br.GetData(); Vector Brimag(&dataBr[n],n);
+      real_t * dataBi = Bi.GetData(); Vector Biimag(&dataBi[n],n);
+      Bcimag *= -1; Brimag *= -1; Biimag *= -1;
+   }
+
+   Yc -= Bc; REQUIRE(Yc.Norml2() < 1.e-12);
+   Yr -= Br; REQUIRE(Yr.Norml2() < 1.e-12);
+   Yi -= Bi; REQUIRE(Yi.Norml2() < 1.e-12);
+
+} // Test case "ComplexMUMPS Solver"
+
+#endif // MFEM_USE_COMPLEX_MUMPS