minor

config/config.hpp

@@ -42,6 +42,9 @@
 #ifdef MFEM_USE_SUPERLU
 #error Building with SuperLU_DIST (MFEM_USE_SUPERLU=YES) requires MPI (MFEM_USE_MPI=YES)
 #endif
+#ifdef MFEM_USE_MUMPS
+#error Building with MUMPS (MFEM_USE_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

@@ -94,6 +94,10 @@
 // Enable MFEM functionality based on the SuperLU library.
 // #define MFEM_USE_SUPERLU
 
+// Enable MFEM functionality based on the MUMPS library.
+// #define MFEM_USE_MUMPS
+// #define MFEM_MUMPS_VERSION @MFEM_MUMPS_VERSION@
+
 // Enable MFEM functionality based on the STRUMPACK library.
 // #define MFEM_USE_STRUMPACK
 

config/config.mk.in

@@ -32,6 +32,7 @@ MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
 MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
+MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@
 MFEM_USE_STRUMPACK     = @MFEM_USE_STRUMPACK@
 MFEM_USE_GINKGO        = @MFEM_USE_GINKGO@
 MFEM_USE_GNUTLS        = @MFEM_USE_GNUTLS@

config/defaults.mk

@@ -120,6 +120,7 @@ MFEM_USE_SUNDIALS      = NO
 MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
+MFEM_USE_MUMPS         = NO
 MFEM_USE_STRUMPACK     = NO
 MFEM_USE_GINKGO        = NO
 MFEM_USE_GNUTLS        = NO
@@ -157,7 +158,7 @@ HYPRE_OPT = -I$(HYPRE_DIR)/include
 HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE
 
 # METIS library configuration
-ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK),NONO)
+ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
    ifeq ($(MFEM_USE_METIS_5),NO)
      METIS_DIR = @MFEM_DIR@/../metis-4.0
      METIS_OPT =
@@ -233,7 +234,7 @@ SCALAPACK_DIR = @MFEM_DIR@/../scalapack-2.0.2
 SCALAPACK_OPT = -I$(SCALAPACK_DIR)/SRC
 SCALAPACK_LIB = -L$(SCALAPACK_DIR)/lib -lscalapack $(LAPACK_LIB)
 
-# MPI Fortran library, needed e.g. by STRUMPACK
+# MPI Fortran library, needed e.g. by STRUMPACK or MUMPS
 # MPICH:
 MPI_FORTRAN_LIB = -lmpifort
 # OpenMPI:
@@ -241,6 +242,11 @@ MPI_FORTRAN_LIB = -lmpifort
 # Additional Fortan library:
 # MPI_FORTRAN_LIB += -lgfortran
 
+# MUMPS library configuration
+MUMPS_DIR =
+MUMPS_OPT = -I$(MUMPS_DIR)/include
+MUMPS_LIB = -Wl,-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
+
 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
 STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)

config/get_mumps_version.cpp

@@ -0,0 +1,33 @@
+// Copyright (c) 2010-2020, 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 "dmumps_c.h"
+#include <string>
+#include <iostream>
+#include <algorithm>
+
+// Macros to expand a macro as a string
+#define STR_EXPAND(s) #s
+#define STR(s) STR_EXPAND(s)
+
+int main()
+{
+#ifdef MUMPS_VERSION
+   const char *ptr = STR(MUMPS_VERSION);
+   std::string s(ptr);
+   s.erase(std::remove(s.begin(), s.end(), '"'), s.end());
+   s.erase(std::remove(s.begin(), s.end(), '.'), s.end());
+   std::cout << s << "\n";
+   return 0;
+#else
+   return -1;
+#endif
+}

config/makefile

@@ -42,6 +42,10 @@ GHV_FLAGS = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(HYPRE_OPT))
 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=)
+GMV_CXX ?= $(MFEM_CXX)
+GMV = get_mumps_version
+GMV_FLAGS = $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(MUMPS_OPT))
 
 $(GHV): $(SRC)$(GHV).cpp
 	$(call mfem-info, Determining HYPRE version ...)
@@ -50,6 +54,13 @@ $(GHV).out: $(GHV)
 	./$(GHV) > $(GHV).out
 .INTERMEDIATE: $(GHV) $(GHV).out
 
+$(GMV): $(SRC)$(GMV).cpp
+	$(call mfem-info, Determining MUMPS version ...)
+	$(GMV_CXX) ${GMV_FLAGS} $(SRC)$(GMV).cpp -o $(GMV)
+$(GMV).out: $(GMV)
+	./$(GMV) > $(GMV).out
+.INTERMEDIATE: $(GMV) $(GMV).out
+
 get-hypre-version: $(GHV).out
 	$(eval MFEM_HYPRE_VERSION:=$(shell cat $(GHV).out))
 	$(if $(MFEM_HYPRE_VERSION),$(eval export MFEM_HYPRE_VERSION)\
@@ -62,10 +73,16 @@ check-smx:
 	$(call mfem-info, MFEM_USE_SIMMETRIX = $(MFEM_USE_SIMMETRIX))
 	$(eval export MFEM_USE_SIMMETRIX)
 
-header: $(if $(MPI),get-hypre-version,) $(if $(SMX),check-smx)
+get-mumps-version: $(GMV).out
+	$(eval MFEM_MUMPS_VERSION:=$(shell cat $(GMV).out))
+	$(if $(MFEM_MUMPS_VERSION),$(eval export MFEM_MUMPS_VERSION)\
+	$(info MUMPS version: $(MFEM_MUMPS_VERSION)),\
+	$(error Unable to determine MUMPS version))
+
+header: $(if $(MPI),get-hypre-version,) $(if $(SMX),check-smx,) $(if $(MUMPS),get-mumps-version,)
 	$(call mfem-info, Writing $(CONFIG_HPP) ...)
 	@set -- && \
-	for def in $${MFEM_DEFINES} $(if $(MPI),MFEM_HYPRE_VERSION) $(SMX); do \
+	for def in $${MFEM_DEFINES} $(if $(MPI),MFEM_HYPRE_VERSION) $(SMX) $(if $(MUMPS),MFEM_MUMPS_VERSION); do \
 	   eval var=\$$$$def && \
 	   if [ "NO" != "$${var}" ]; then \
 	      set -- "$$@" -e "s|// \(#define $${def} \)|\1|" && \

examples/ex1p_mumps.cpp

@@ -0,0 +1,234 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool visualization = true;
+   int sr = 1;
+   int pr = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&sr, "-sr", "--serial_ref",
+                  "Number of serial refinements");        
+   args.AddOption(&pr, "-pr", "--parallel_ref",
+                  "Number of parallel refinements");                               
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      for (int l = 0; l < sr; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      for (int l = 0; l < pr; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, fec);
+   HYPRE_Int size = fespace.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //     domain integrator.
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   if (static_cond) { a.EnableStaticCondensation(); }
+   a.Assemble();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // // 13. Solve the linear system A X = B.
+   // //     * With full assembly, use the BoomerAMG preconditioner from hypre.
+   // //     * With partial assembly, use Jacobi smoothing, for now.
+   StopWatch chrono;
+   chrono.Clear();
+   chrono.Start();
+   HypreBoomerAMG *prec = new HypreBoomerAMG;
+   prec->SetPrintLevel(0);
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-13);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(0);
+   if (prec) { cg.SetPreconditioner(*prec); }
+   cg.SetOperator(A);
+   cg.Mult(B, X);
+   delete prec;
+
+   if (myid == 0)
+   {
+      cout << "PCG-AMG time: " << chrono.RealTime() << endl;
+   }
+
+   chrono.Clear();
+   chrono.Start();
+   {
+      MUMPSSolver MA;
+      MA.SetMatrixSymType(0);
+      MA.SetOperator(A);
+      Vector Y(X.Size());
+      MA.Mult(B,Y);
+      Y-=X;
+      cout << "Mumps Diff norm = " << Y.Norml2() << endl;
+   }
+   if (myid == 0)
+   {
+      cout << "mumps time: " << chrono.RealTime() << endl;
+   }
+
+   chrono.Clear();
+   chrono.Start();
+
+   {
+      CPardisoSolver pardiso(A.GetComm());
+      // pardiso.SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
+      pardiso.SetMatrixType(CPardisoSolver::MatType::REAL_UNSYMMETRIC);
+      pardiso.SetPrintLevel(0);
+      pardiso.SetOperator(A);
+      Vector Y(X.Size());
+      pardiso.Mult(B, Y);
+      Y-=X;
+      cout << "Pardiso Diff norm = " << Y.Norml2() << endl;
+   }
+
+   if (myid == 0)
+   {
+      cout << "pardiso time: " << chrono.RealTime() << endl;
+   }
+
+   {
+      SuperLURowLocMatrix SA(A);
+      SuperLUSolver superlu(MPI_COMM_WORLD);
+      superlu.SetPrintStatistics(false);
+      superlu.SetSymmetricPattern(false);
+      superlu.SetColumnPermutation(superlu::PARMETIS);
+      superlu.SetOperator(SA);
+      Vector Y(X.Size());
+      superlu.Mult(B, Y);
+      Y-=X;
+      cout << "Superlu Diff norm = " << Y.Norml2() << endl;
+   }
+
+   if (myid == 0)
+   {
+      cout << "superlu time: " << chrono.RealTime() << endl;
+   }
+
+   a.RecoverFEMSolution(X, b, x);
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   delete fec;
+   MPI_Finalize();
+
+   return 0;
+}

linalg/CMakeLists.txt

@@ -76,6 +76,13 @@ if (MFEM_USE_GINKGO)
   list(APPEND HDRS ginkgo.hpp)
 endif()
 
+if (MFEM_USE_MUMPS)
+  list(APPEND SRCS mumps.cpp)
+  # If this list (HDRS -> HEADERS) is used for install, we probably want the
+  # header added all the time.
+  list(APPEND HDRS mumps.hpp)
+endif()
+
 if (MFEM_USE_SUNDIALS)
   list(APPEND SRCS sundials.cpp)
   list(APPEND HDRS sundials.hpp)

linalg/cpardiso.cpp

@@ -1,12 +1,15 @@
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MKL_CPARDISO
+#ifdef MFEM_USE_MPI
+
 #include "cpardiso.hpp"
 #include "hypre.hpp"
 #include <algorithm>
 #include <vector>
 #include <numeric>
 
-#ifdef MFEM_USE_MPI
-#ifdef MFEM_USE_MKL_CPARDISO
-
 namespace mfem
 {
 CPardisoSolver::CPardisoSolver(MPI_Comm comm) : comm_(comm)

linalg/cpardiso.hpp

@@ -12,12 +12,13 @@
 #ifndef MFEM_CPARDISO
 #define MFEM_CPARDISO
 
-#ifdef MFEM_USE_MPI
-#ifdef MFEM_USE_MKL_CPARDISO
-
 #include "../config/config.hpp"
-#include "mkl_cluster_sparse_solver.h"
+
+#ifdef MFEM_USE_MKL_CPARDISO
+#ifdef MFEM_USE_MPI
 #include "operator.hpp"
+#include <mpi.h>
+#include "mkl_cluster_sparse_solver.h"
 
 namespace mfem
 {

linalg/linalg.hpp

@@ -45,6 +45,10 @@
 #include "hypre_parcsr.hpp"
 #include "hypre.hpp"
 
+#ifdef MFEM_USE_MUMPS
+#include "mumps.hpp"
+#endif
+
 #ifdef MFEM_USE_PETSC
 #include "petsc.hpp"
 #endif

linalg/mumps.cpp

@@ -0,0 +1,387 @@
+// Copyright (c) 2010-2020, 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 "../config/config.hpp"
+
+#ifdef MFEM_USE_MUMPS
+#ifdef MFEM_USE_MPI
+
+#include "mumps.hpp"
+
+namespace mfem
+{
+MUMPSSolver::~MUMPSSolver()
+{
+   if (id)
+   {
+      id->job = -2;
+      dmumps_c(id);
+      delete[] J;
+      delete[] I;
+      delete [] data;
+   }
+}
+
+void MUMPSSolver::SetParameters()
+{
+   // output messages
+   id->ICNTL(1) = -1;
+
+   // Diagnosting printing
+   id->ICNTL(2) = -1;
+
+   // Global info on host
+   id->ICNTL(3) = -1;
+
+   // Level of error printing
+   id->ICNTL(4) = 0;
+
+   //input matrix format (assembled)
+   id->ICNTL(5) = 0;
+
+   // Use A or A^T
+   id->ICNTL(9) = 1;
+
+   // Iterative refinement (disabled)
+   id->ICNTL(10) = 0;
+
+   // Error analysis-statistics (disabled)
+   id->ICNTL(11) = 0;
+
+   // Use of ScaLAPACK (Parallel factorization on root)
+   id->ICNTL(13) = 0;
+
+   // Percentage increase of estimated workspace (default = 20%)
+   id->ICNTL(14) = 20;
+
+   // Number of OpenMP threads (default)
+   id->ICNTL(16) = 0;
+
+   // Matrix input format (distributed)
+   id->ICNTL(18) = 3;
+
+   // Schur complement (no Schur complement matrix returned)
+   id->ICNTL(19) = 0;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Distributed RHS and Sol
+   id->ICNTL(20) = 10;
+
+   id->ICNTL(21) = 1;
+#else
+   // Centralized RHS and Sol
+   id->ICNTL(20) = 0;
+
+   id->ICNTL(21) = 0;
+#endif
+   // Out of core factorization and solve (disabled)
+   id->ICNTL(22) = 0;
+
+   // Max size of working memory (default = based on estimates)
+   id->ICNTL(23) = 0;
+}
+
+void MUMPSSolver::SetOperator(const Operator &op)
+{
+   // Verify that the operator is a HypreParMatrix
+   auto APtr = dynamic_cast<const HypreParMatrix *>(&op);
+   MFEM_VERIFY(APtr, "Not compatible matrix type");
+   height = op.Height();
+   width = op.Width();
+
+   comm = APtr->GetComm();
+   MPI_Comm_size(comm, &numProcs);
+   MPI_Comm_rank(comm, &myid);
+
+   hypre_ParCSRMatrix *parcsr_op
+      = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>(*APtr);
+   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+#if MFEM_HYPRE_VERSION >= 21600
+   hypre_CSRMatrixBigJtoJ(csr_op);
+#endif
+
+   int *Iptr = csr_op->i;
+   int *Jptr = csr_op->j;
+   int n_loc = csr_op->num_rows;
+
+   row_start = parcsr_op->first_row_index;
+
+   int nnz;
+   if (sym)
+   {
+      // count nnz;
+      nnz = 0;
+      int k = 0;
+      for (int i = 0; i < n_loc; i++)
+      {
+         for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+         {
+            int ii = row_start + i + 1;
+            int jj = Jptr[k] + 1;
+            k++;
+            if (ii>=jj) { nnz++; }
+         }
+      }
+   }
+   else
+   {
+      nnz = csr_op->num_nonzeros;
+   }
+
+   I = new int[nnz];
+   J = new int[nnz];
+
+   int k = 0;
+   if (sym)
+   {
+      int l = 0;
+      data = new double[nnz];
+      for (int i = 0; i < n_loc; i++)
+      {
+         for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+         {
+            // Global I and J indices in 1-based index (for fortran)
+            int ii = row_start + i + 1;
+            int jj = Jptr[k] + 1;
+            if (ii>=jj)
+            {
+               I[l] = ii;
+               J[l] = jj;
+               data[l++] = csr_op->data[k];
+            }
+            k++;
+         }
+      }
+   }
+   else
+   {
+      for (int i = 0; i < n_loc; i++)
+      {
+         for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+         {
+            // Global I and J indices in 1-based index (for fortran)
+            I[k] = row_start + i + 1;
+            J[k] = Jptr[k] + 1;
+            k++;
+         }
+      }
+      data = csr_op->data;
+   }
+
+   // new MUMPS object
+   id = new DMUMPS_STRUC_C;
+
+   // C to Fortran communicator
+   id->comm_fortran = (MUMPS_INT) MPI_Comm_c2f(comm);
+
+   // Host is involved in computation
+   id->par = 1;
+
+   // Unsymmetric matrix
+   id->sym = sym;
+
+   // Mumps init
+   id->job = -1;
+   dmumps_c(id);
+
+   // Set MUMPS default parameters
+   SetParameters();
+
+   // Global number of rows
+   id->n = parcsr_op->global_num_rows;
+
+   // Number of non zeros on the processor
+   id->nnz_loc = nnz;
+
+   // Distributed row array
+   id->irn_loc = I;
+
+   // Distributed column array
+   id->jcn_loc = J;
+
+   // Distributed data array
+   id->a_loc = data;
+
+   // MUMPS Analysis
+   id->job = 1;
+   dmumps_c(id);
+
+   // MUMPS Factorization
+   id->job = 2;
+   dmumps_c(id);
+
+   // matrix can be destroyed now
+   hypre_CSRMatrixDestroy(csr_op);
+   if (!sym) { data = nullptr; }
+
+#if MFEM_MUMPS_VERSION >= 530
+   irhs_loc.SetSize(n_loc);
+   for (int i = 0; i < n_loc; i++)
+   {
+      irhs_loc[i] = row_start + i + 1;
+   }
+   row_starts.SetSize(numProcs);
+   MPI_Allgather(&row_start, 1, MPI_INT, row_starts, 1, MPI_INT, comm);
+   sol_loc.SetSize(id->INFO(23));
+   isol_loc.SetSize(id->INFO(23));
+#else
+   if (myid == 0)
+   {
+      rhs_glob.SetSize(parcsr_op->global_num_rows);
+      recv_counts.SetSize(numProcs);
+   }
+   MPI_Gather(&n_loc, 1, MPI_INT, recv_counts, 1, MPI_INT, 0, comm);
+   if (myid == 0)
+   {
+      displs.SetSize(numProcs); displs[0] = 0;
+      int s = 0;
+      for (int k = 0; k < numProcs-1; k++)
+      {
+         s += recv_counts[k];
+         displs[k+1] = s;
+      }
+   }
+#endif
+}
+
+void MUMPSSolver::Mult(const Vector &x, Vector &y) const
+{
+#if MFEM_MUMPS_VERSION >= 530
+   id->nloc_rhs = x.Size();
+   id->lrhs_loc = x.Size();
+   id->rhs_loc = x.GetData();
+   id->irhs_loc = const_cast<int *>(irhs_loc.GetData());
+   id->sol_loc = sol_loc.GetData();
+   id->lsol_loc = id->INFO(23);
+   id->isol_loc = const_cast<int *>(isol_loc.GetData());
+   id->job = 3;
+   dmumps_c(id);
+   RedistributeSol(isol_loc, sol_loc, y);
+#else
+   MPI_Gatherv(x.GetData(), x.Size(), MPI_DOUBLE,
+               rhs_glob.GetData(), recv_counts,
+               displs, MPI_DOUBLE, 0, comm);
+   if (myid == 0)
+   {
+      id->rhs = rhs_glob.GetData();
+   }
+   id->job = 3;
+   dmumps_c(id);
+   MPI_Scatterv(rhs_glob.GetData(), recv_counts, displs,
+                MPI_DOUBLE, y.GetData(), y.Size(),
+                MPI_DOUBLE, 0, comm);
+#endif
+}
+
+void MUMPSSolver::MultTranspose(const Vector &x, Vector &y) const
+{
+   id->ICNTL(9) = 0;
+   Mult(x,y);
+}
+
+#if MFEM_MUMPS_VERSION >= 530
+int MUMPSSolver::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 std::distance(row_starts_.begin(), up) - 1;
+}
+
+void MUMPSSolver::RedistributeSol(const Array<int> &row_map,
+                                  const Vector &x,
+                                  Vector &y) const
+{
+   MFEM_VERIFY(row_map.Size() == x.Size(), "Inconcistent sizes");
+   int size = x.Size();
+
+   // compute send_count
+   Array<int> send_count(numProcs);
+   send_count = 0;
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1; //fix to 0-based indexing
+      int row_rank = GetRowRank(j, row_starts);
+      send_count[row_rank]++; // both for val and global index
+   }
+
+   // compute recv_count
+   Array<int> recv_count(numProcs);
+   MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
+
+   // compute offsets
+   Array<int> send_displ(numProcs);
+   send_displ[0] = 0;
+   Array<int> recv_displ(numProcs);
+   recv_displ[0] = 0;
+   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];
+   }
+   int sbuff_size = send_count.Sum();
+   int rbuff_size = recv_count.Sum();
+
+   Array<int> sendbuf_index(sbuff_size);
+   sendbuf_index = 0;
+   Array<double> sendbuf_value(sbuff_size);
+   sendbuf_value = 0;
+   Array<int> soffs(numProcs);
+   soffs = 0;
+
+   // Fill in send buffers
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1; //fix to 0-based indexing
+      int row_rank = GetRowRank(j, row_starts);
+      int k = send_displ[row_rank] + soffs[row_rank];
+      sendbuf_index[k] = j;
+      sendbuf_value[k] = x(i);
+      soffs[row_rank]++;
+   }
+
+   // communicate
+   Array<int> recvbuf_index(rbuff_size);
+   Array<double> recvbuf_value(rbuff_size);
+   MPI_Alltoallv(sendbuf_index,
+                 send_count,
+                 send_displ,
+                 MPI_INT,
+                 recvbuf_index,
+                 recv_count,
+                 recv_displ,
+                 MPI_INT,
+                 comm);
+   MPI_Alltoallv(sendbuf_value,
+                 send_count,
+                 send_displ,
+                 MPI_DOUBLE,
+                 recvbuf_value,
+                 recv_count,
+                 recv_displ,
+                 MPI_DOUBLE,
+                 comm);
+
+   // Unpack recv buffer
+   for (int i = 0; i < rbuff_size; i++)
+   {
+      int local_index = recvbuf_index[i] - row_start;
+      y(local_index) = recvbuf_value[i];
+   }
+}
+#endif
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MUMPS

linalg/mumps.hpp

@@ -0,0 +1,104 @@
+// Copyright (c) 2010-2020, 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_MUMPS
+#define MFEM_MUMPS
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MUMPS
+#ifdef MFEM_USE_MPI
+#include "operator.hpp"
+#include "hypre.hpp"
+
+#include <mpi.h>
+#include "dmumps_c.h"
+#include <vector>
+
+namespace mfem
+{
+class MUMPSSolver : public mfem::Solver
+{
+public:
+   // Default Constructor.
+   MUMPSSolver() {}
+
+   void SetMatrixSymType(int sym_) { sym = (sym_>2) ? 0 : sym_ ; }
+
+   // Factor and solve the linear system y = Op^{-1} x.
+   void Mult(const Vector &x, Vector &y) const;
+
+   void MultTranspose(const Vector &x, Vector &y) const;
+
+   // Set the operator.
+   void SetOperator(const Operator &op);
+
+   // Default destructor.
+   ~MUMPSSolver();
+
+private:
+
+   MPI_Comm comm;
+
+   int numProcs;
+
+   int myid;
+
+   int sym=0;
+
+   int row_start;
+
+   int *I;
+
+   int *J;
+
+   double * data;
+
+   // MUMPS workspace
+   // macro s.t. indices match MUMPS documentation
+#define ICNTL(I) icntl[(I) -1]
+#define INFO(I) info[(I) -1]
+
+   DMUMPS_STRUC_C *id=nullptr;
+
+   void SetParameters();
+
+#if MFEM_MUMPS_VERSION >= 530
+
+   Array<int> row_starts;
+
+   Array<int> irhs_loc;
+
+   Array<int> isol_loc;
+
+   Vector sol_loc;
+
+   int GetRowRank(int i, const Array<int> &row_starts_) const;
+
+   void RedistributeSol(const Array<int> &row_map,
+                        const Vector &x,
+                        Vector &y) const;
+#else
+   Array<int> recv_counts;
+
+   Array<int> displs;
+
+   Vector rhs_glob;
+
+#endif
+
+}; // mfem::MUMPSSolver class
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+#endif // MFEM_USE_MUMPS
+#endif // MFEM_MUMPS

makefile

@@ -204,7 +204,11 @@ CXXFLAGS ?= $(OPTIM_FLAGS)
 # MPI configuration
 ifneq ($(MFEM_USE_MPI),YES)
    MFEM_HOST_CXX = $(CXX)
+<<<<<<< HEAD
+   PKGS_NEED_MPI = SUPERLU MUMPS STRUMPACK PETSC PUMI SLEPC
+=======
    PKGS_NEED_MPI = SUPERLU STRUMPACK PETSC PUMI SLEPC MKL_CPARDISO
+>>>>>>> mkl-pardiso
    $(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),))
@@ -259,7 +263,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 METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
  OCCA CEED RAJA UMPIRE MKL_CPARDISO
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -321,7 +325,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_ZLIB MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE\
  MFEM_USE_OPENMP MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE\
  MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE MFEM_USE_GINKGO\
- MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GNUTLS\
+ MFEM_USE_SUPERLU MFEM_USE_MUMPS MFEM_USE_STRUMPACK MFEM_USE_GNUTLS\
  MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD\
@@ -359,7 +363,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
+   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT
 endif
 
 # If we have 'install' target, export variables used by config/makefile
@@ -627,6 +631,7 @@ status info:
 	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
+	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))
 	$(info MFEM_USE_GINKGO        = $(MFEM_USE_GINKGO))
 	$(info MFEM_USE_GNUTLS        = $(MFEM_USE_GNUTLS))