minor

CHANGELOG

@@ -137,6 +137,9 @@ Linear and nonlinear solvers
 - Added partially assembled convergent diagonal preconditioner for adaptively
   refined meshes (i.e. non-conforming finite element spaces), see Example 6/6p.
 
+- Added an interface to the Intel MKL Parallel Direct Sparse Solver for
+  Clusters. An example usage of the interface is shown in Example 11p.
+
 New and updated examples and miniapps
 -------------------------------------
 - Added a new example, Example 25/25p, to demonstrate the use of a Perfectly

CMakeLists.txt

@@ -346,6 +346,12 @@ if (MFEM_USE_ADIOS2)
    find_package(ADIOS2 REQUIRED)
 endif()
 
+if (MFEM_USE_MKL_CPARDISO)
+  if (MFEM_USE_MPI)
+   find_package(MKL_CPARDISO REQUIRED MKL_SEQUENTIAL MKL_LP64 MKL_MPI_WRAPPER)
+  endif()
+endif()
+
 # MFEM_TIMER_TYPE
 if (NOT DEFINED MFEM_TIMER_TYPE)
   if (APPLE)
@@ -372,7 +378,7 @@ endif()
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     SLEPC MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
     MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2
-    CUSPARSE)
+    CUSPARSE MKL_CPARDISO)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")

INSTALL

@@ -486,6 +486,13 @@ MFEM_USE_CEED = YES/NO
    library for performant high-order operator evaluation developed by the Center
    for Efficient Exascale Discretizations in the Exascale Computing Project.
 
+MFEM_USE_MKL_CPARDISO = YES/NO
+   Enables the interface to the Intel MKL Parallel Direct Sparse Solver for
+   Clusters. Make sure to set the correct values for MKL_MPI_WRAPPER and
+   MKL_LIBRARY_SUBDIR as shown in defaults.mk. If you configure MFEM with
+   MFEM_USE_LAPACK=YES, verify that the MKL LAPACK libraries are used. The
+   OpenMP capabilities are disabled at link time.
+
 MFEM_BUILD_TAG = (any value)
    An optional tag to characterize the build. Exported to config/config.mk.
    Can be used to identify the MFEM build from other makefiles.

config/cmake/config.hpp.in

@@ -156,4 +156,7 @@
 // library.
 #cmakedefine MFEM_USE_SIMMETRIX
 
+// Enable interface to the MKL CPardiso library.
+#cmakedefine MFEM_USE_MKL_CPARDISO
+
 #endif // MFEM_CONFIG_HEADER

config/cmake/modules/FindMKL_CPARDISO.cmake

@@ -0,0 +1,106 @@
+# 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.
+
+# Defines the following variables:
+#   - MKL_CPARDISO_FOUND
+#   - MKL_CPARDISO_LIBRARIES
+#   - MKL_CPARDISO_INCLUDE_DIRS
+
+if(NOT MKL_MPI_WRAPPER_LIB)
+  message(FATAL_ERROR "MKL CPardiso enabled but no MKL MPI Wrapper lib specified")
+endif()
+
+if(NOT MKL_LIBRARY_DIR)
+  message(WARNING "Using default MKL library path. Double check the variable MKL_LIBRARY_DIR")
+  set(MKL_LIBRARY_DIR "lib")
+endif()
+
+include(MfemCmakeUtilities)
+mfem_find_package(MKL_CPARDISO MKL_CPARDISO
+    MKL_CPARDISO_DIR "include" mkl_cluster_sparse_solver.h ${MKL_LIBRARY_DIR} mkl_core
+  "Paths to headers required by MKL CPardiso." "Libraries required by MKL CPARDISO."
+  ADD_COMPONENT MKL_LP64 "include" "" ${MKL_LIBRARY_DIR} mkl_intel_lp64
+  ADD_COMPONENT MKL_SEQUENTIAL "include" "" ${MKL_LIBRARY_DIR} mkl_sequential
+  ADD_COMPONENT MKL_MPI_WRAPPER "include" "" ${MKL_LIBRARY_DIR} ${MKL_MPI_WRAPPER_LIB}
+  CHECK_BUILD MKL_CPARDISO_VERSION_OK TRUE
+  "
+#include <mpi.h>
+#include <mkl.h>
+#include <mkl_cluster_sparse_solver.h>
+int main (void)
+{
+    MKL_INT n = 5;
+    MKL_INT ia[6]  = { 1, 4, 6, 9, 12, 14};
+    MKL_INT ja[13] = {  1, 2, 4,        /* index of non-zeros in 1 row*/
+                        1, 2,           /* index of non-zeros in 2 row*/
+                        3, 4, 5,        /* index of non-zeros in 3 row*/
+                        1, 3, 4,        /* index of non-zeros in 4 row*/
+                        2, 5            /* index of non-zeros in 5 row*/
+    };
+    double a[13] = {
+                         1.0, -1.0, /*0*/ -3.0, /*0*/
+                        -2.0,  5.0, /*0*/ /*0*/ /*0*/
+                        /*0*/  4.0,  6.0,  4.0, /*0*/
+                        -4.0, /*0*/  2.0,  7.0, /*0*/
+                        /*0*/  8.0, /*0*/ /*0*/ -5.0
+                    };
+
+    MKL_INT mtype = 11; /* set matrix type to \"real unsymmetric matrix\" */
+    MKL_INT nrhs  = 1;  /* Number of right hand sides. */
+    double b[5], x[5], bs[5], res, res0; /* RHS and solution vectors. */
+
+    /* Internal solver memory pointer pt
+     *       32-bit:      int pt[64] or void *pt[64];
+     *       64-bit: long int pt[64] or void *pt[64]; */
+    void *pt[64] = { 0 };
+
+    /* Cluster Sparse Solver control parameters. */
+    MKL_INT iparm[64] = { 0 };
+    MKL_INT maxfct, mnum, phase, msglvl, error;
+
+    /* Auxiliary variables. */
+    double  ddum; /* Double dummy   */
+    MKL_INT idum; /* Integer dummy. */
+    MKL_INT i, j;
+    int     mpi_stat = 0;
+    int     argc = 0;
+    int     comm, rank;
+    char*   uplo;
+    char**  argv;
+
+    mpi_stat = MPI_Init( &argc, &argv );
+    mpi_stat = MPI_Comm_rank( MPI_COMM_WORLD, &rank );
+    comm =  MPI_Comm_c2f( MPI_COMM_WORLD );
+
+    iparm[ 0] =  1; /* Solver default parameters overriden with provided by iparm */
+    iparm[ 1] =  2; /* Use METIS for fill-in reordering */
+    iparm[ 5] =  0; /* Write solution into x */
+    iparm[ 7] =  2; /* Max number of iterative refinement steps */
+    iparm[ 9] = 13; /* Perturb the pivot elements with 1E-13 */
+    iparm[10] =  1; /* Use nonsymmetric permutation and scaling MPS */
+    iparm[12] =  1; /* Switch on Maximum Weighted Matching algorithm (default for non-symmetric) */
+    iparm[17] = -1; /* Output: Number of nonzeros in the factor LU */
+    iparm[18] = -1; /* Output: Mflops for LU factorization */
+    iparm[26] =  1; /* Check input data for correctness */
+    iparm[39] =  0; /* Input: matrix/rhs/solution stored on master */
+    maxfct = 1; /* Maximum number of numerical factorizations. */
+    mnum   = 1; /* Which factorization to use. */
+    msglvl = 1; /* Print statistical information in file */
+    error  = 0; /* Initialize error flag */
+
+    phase = 11;
+    cluster_sparse_solver ( pt, &maxfct, &mnum, &mtype, &phase,
+                &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &comm, &error );
+
+    mpi_stat = MPI_Finalize();
+    return error;
+}
+")

config/config.hpp

@@ -42,9 +42,15 @@
 #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
+#ifdef MFEM_USE_MKL_CPARDISO
+#error Building with MKL CPARDISO (MFEM_USE_MKL_CPARDISO=YES) requires MPI (MFEM_USE_MPI=YES)
+#endif
 #ifdef MFEM_USE_PETSC
 #error Building with PETSc (MFEM_USE_PETSC=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
 
@@ -163,4 +167,7 @@
 // library.
 // #define MFEM_USE_SIMMETRIX
 
+// Enable interface to the MKL CPardiso library.
+// #define MFEM_USE_MKL_CPARDISO
+
 #endif // MFEM_CONFIG_HEADER

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@
@@ -52,6 +53,7 @@ MFEM_USE_CEED          = @MFEM_USE_CEED@
 MFEM_USE_UMPIRE        = @MFEM_USE_UMPIRE@
 MFEM_USE_SIMD          = @MFEM_USE_SIMD@
 MFEM_USE_ADIOS2        = @MFEM_USE_ADIOS2@
+MFEM_USE_MKL_CPARDISO  = @MFEM_USE_MKL_CPARDISO@
 
 # Compiler, compile options, and link options
 MFEM_CXX       = @MFEM_CXX@

config/defaults.cmake

@@ -52,6 +52,7 @@ option(MFEM_USE_CEED "Enable CEED" OFF)
 option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
 option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" OFF)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
+option(MFEM_USE_MKL_CPARDISO "Enable MKL CPardiso" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -180,6 +181,10 @@ set(HIOP_DIR "${MFEM_DIR}/../hiop/install" CACHE STRING
 set(HIOP_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
     "Packages that HiOp depends on.")
 
+set(MKL_CPARDISO_DIR "" CACHE STRING "MKL installation path.")
+set(MKL_MPI_WRAPPER_LIB "mkl_blacs_mpich_lp64" CACHE STRING "MKL MPI wrapper library")
+set(MKL_LIBRARY_DIR "" CACHE STRING "Custom library subdirectory")
+
 set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
 set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
 set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")

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
@@ -140,6 +141,7 @@ MFEM_USE_CEED          = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
+MFEM_USE_MKL_CPARDISO  = NO
 
 # Compile and link options for zlib.
 ZLIB_DIR =
@@ -156,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 =
@@ -232,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:
@@ -240,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)
@@ -372,6 +379,15 @@ UMPIRE_DIR = @MFEM_DIR@/../umpire
 UMPIRE_OPT = -I$(UMPIRE_DIR)/include
 UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -lumpire
 
+# MKL CPardiso library configuration
+MKL_CPARDISO_DIR ?=
+MKL_MPI_WRAPPER ?= mkl_blacs_mpich_lp64
+MKL_LIBRARY_SUBDIR ?= lib
+MKL_CPARDISO_OPT = -I$(MKL_CPARDISO_DIR)/include
+MKL_CPARDISO_LIB = -Wl,-rpath,$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR)\
+                   -L$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR) -l$(MKL_MPI_WRAPPER)\
+                   -lmkl_intel_lp64 -lmkl_sequential -lmkl_core
+
 # If YES, enable some informational messages
 VERBOSE = NO
 

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/ex11p.cpp

@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
    int seed = 75;
    bool slu_solver  = false;
    bool sp_solver = false;
+   bool pardiso_solver = false;
    bool visualization = 1;
 
    OptionsParser args(argc, argv);
@@ -95,6 +96,14 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_STRUMPACK
    args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
                   "--no-strumpack", "Use the STRUMPACK Solver.");
+#endif
+#ifdef MFEM_USE_MKL_CPARDISO
+   args.AddOption(&pardiso_solver,
+                  "-pardiso",
+                  "--pardiso",
+                  "-no-pardiso",
+                  "--no-pardiso",
+                  "Use the MKL Cluster Pardiso Solver.");
 #endif
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
@@ -236,7 +245,7 @@ int main(int argc, char *argv[])
    //    preconditioner for A to be used within the solver. Set the matrices
    //    which define the generalized eigenproblem A x = lambda M x.
    Solver * precond = NULL;
-   if (!slu_solver && !sp_solver)
+   if (!slu_solver && !sp_solver && !pardiso_solver)
    {
       HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
       amg->SetPrintLevel(0);
@@ -268,10 +277,19 @@ int main(int argc, char *argv[])
          strumpack->SetFromCommandLine();
          precond = strumpack;
       }
+#endif
+#ifdef MFEM_USE_MKL_CPARDISO
+      if (pardiso_solver)
+      {
+         auto pardiso = new CPardisoSolver(A->GetComm());
+         pardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
+         pardiso->SetPrintLevel(1);
+         pardiso->SetOperator(*A);
+         precond = pardiso;
+      }
 #endif
    }
 
-
    HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
    lobpcg->SetNumModes(nev);
    lobpcg->SetRandomSeed(seed);

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

examples/makefile

@@ -119,6 +119,11 @@ ex11p-test-superlu: ex11p
 	@$(call mfem-test,$<, $(RUN_MPI), SuperLU_DIST example,--superlu)
 test-par-YES: ex11p-test-superlu
 endif
+ifeq ($(MFEM_USE_MKL_CPARDISO),YES)
+ex11p-test-superlu: ex11p
+	@$(call mfem-test,$<, $(RUN_MPI), MKL_CPARDISO example,--pardiso)
+test-par-YES: ex11p-test-pardiso
+endif
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

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)
@@ -98,6 +105,11 @@ if (MFEM_USE_HIOP)
   list(APPEND HDRS hiop.hpp)
 endif()
 
+if (MFEM_USE_MKL_CPARDISO)
+  list(APPEND SRCS cpardiso.cpp)
+  list(APPEND HDRS cpardiso.hpp)
+endif()
+
 convert_filenames_to_full_paths(SRCS)
 convert_filenames_to_full_paths(HDRS)
 

linalg/cpardiso.cpp

@@ -0,0 +1,233 @@
+
+#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>
+
+namespace mfem
+{
+CPardisoSolver::CPardisoSolver(MPI_Comm comm) : comm_(comm)
+{
+   // Solver default parameters overridden with provided by iparm
+   iparm[0] = 1;
+   // Use METIS for fill-in reordering
+   iparm[1] = 2;
+   // Write solution into x
+   iparm[5] = 0;
+   // Max number of iterative refinement steps
+   iparm[7] = 2;
+   // Perturb the pivot elements with 1E-13
+   iparm[9] = 13;
+   // Use non-symmetric permutation and scaling MPS
+   iparm[10] = 1;
+   // Switch on Maximum Weighted Matching algorithm (default for non-symmetric)
+   iparm[12] = 1;
+   // Output: Number of non-zeros in the factor LU
+   iparm[17] = -1;
+   // Output: Mflops for LU factorization
+   iparm[18] = -1;
+   // Check input data for correctness
+   iparm[26] = 1;
+   // 0-based indexing
+   iparm[34] = 1;
+   // All inputs are distributed between MPI processes
+   iparm[39] = 2;
+   // Maximum number of numerical factorizations
+   maxfct = 1;
+   // Which factorization to use
+   mnum = 1;
+   // Print statistical information in file
+   msglvl = 0;
+   // Initialize error flag
+   error = 0;
+   // Real unsymmetric matrix
+   mtype = MatType::REAL_UNSYMMETRIC;
+   // Number of right hand sides
+   nrhs = 1;
+};
+
+void CPardisoSolver::SetOperator(const Operator &op)
+{
+   auto hypreParMat = dynamic_cast<const HypreParMatrix &>(op);
+
+   MFEM_ASSERT(hypreParMat, "Must pass HypreParMatrix as Operator");
+
+   auto parcsr_op = static_cast<hypre_ParCSRMatrix *>(
+                       const_cast<HypreParMatrix &>(hypreParMat));
+
+   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+#if MFEM_HYPRE_VERSION >= 21600
+   hypre_CSRMatrixBigJtoJ(csr_op);
+#endif
+
+   m = parcsr_op->global_num_rows;
+   first_row = parcsr_op->first_row_index;
+   nnz_loc = csr_op->num_nonzeros;
+   m_loc = csr_op->num_rows;
+
+   height = m_loc;
+   width = m_loc;
+
+   double *csr_nzval = csr_op->data;
+   int *csr_colind = csr_op->j;
+
+   delete[] csr_rowptr;
+   delete[] reordered_csr_colind;
+   delete[] reordered_csr_nzval;
+   csr_rowptr = new int[m_loc + 1];
+   reordered_csr_colind = new int[nnz_loc];
+   reordered_csr_nzval = new double[nnz_loc];
+
+   for (int i = 0; i <= m_loc; i++)
+   {
+      csr_rowptr[i] = (csr_op->i)[i];
+   }
+
+   // CPardiso expects the column indices to be sorted for each row
+   std::vector<int> permutation_idx(nnz_loc);
+   std::iota(permutation_idx.begin(), permutation_idx.end(), 0);
+   for (int i = 0; i < m_loc; i++)
+   {
+      std::sort(permutation_idx.begin() + csr_rowptr[i],
+                permutation_idx.begin() + csr_rowptr[i + 1],
+                [csr_colind](int i1, int i2)
+      {
+         return csr_colind[i1] < csr_colind[i2];
+      });
+   }
+
+   for (int i = 0; i < nnz_loc; i++)
+   {
+      reordered_csr_colind[i] = csr_colind[permutation_idx[i]];
+      reordered_csr_nzval[i] = csr_nzval[permutation_idx[i]];
+   }
+
+   hypre_CSRMatrixDestroy(csr_op);
+
+   // The number of row in global matrix, rhs element and solution vector that
+   // begins the input domain belonging to this MPI process
+   iparm[40] = first_row;
+
+   // The number of row in global matrix, rhs element and solution vector that
+   // ends the input domain belonging to this MPI process
+   iparm[41] = first_row + m_loc - 1;
+
+   // Analyze inputs
+   phase = 11;
+   cluster_sparse_solver(pt,
+                         &maxfct,
+                         &mnum,
+                         &mtype,
+                         &phase,
+                         &m,
+                         reordered_csr_nzval,
+                         csr_rowptr,
+                         reordered_csr_colind,
+                         &idum,
+                         &nrhs,
+                         iparm,
+                         &msglvl,
+                         &ddum,
+                         &ddum,
+                         &comm_,
+                         &error);
+
+   MFEM_ASSERT(error == 0, "Pardiso analyze input error");
+
+   // Numerical factorization
+   phase = 22;
+   cluster_sparse_solver(pt,
+                         &maxfct,
+                         &mnum,
+                         &mtype,
+                         &phase,
+                         &m,
+                         reordered_csr_nzval,
+                         csr_rowptr,
+                         reordered_csr_colind,
+                         &idum,
+                         &nrhs,
+                         iparm,
+                         &msglvl,
+                         &ddum,
+                         &ddum,
+                         &comm_,
+                         &error);
+
+   MFEM_ASSERT(error == 0, "Pardiso factorization input error");
+}
+
+void CPardisoSolver::Mult(const Vector &b, Vector &x) const
+{
+   // Solve
+   phase = 33;
+   cluster_sparse_solver(pt,
+                         &maxfct,
+                         &mnum,
+                         &mtype,
+                         &phase,
+                         &m,
+                         reordered_csr_nzval,
+                         csr_rowptr,
+                         reordered_csr_colind,
+                         &idum,
+                         &nrhs,
+                         iparm,
+                         &msglvl,
+                         b.GetData(),
+                         x.GetData(),
+                         &comm_,
+                         &error);
+
+   MFEM_ASSERT(error == 0, "Pardiso solve error");
+}
+
+void CPardisoSolver::SetPrintLevel(int print_level)
+{
+   msglvl = print_level;
+}
+
+void CPardisoSolver::SetMatrixType(MatType mat_type)
+{
+   mtype = mat_type;
+}
+
+CPardisoSolver::~CPardisoSolver()
+{
+   // Release all internal memory
+   phase = -1;
+   cluster_sparse_solver(pt,
+                         &maxfct,
+                         &mnum,
+                         &mtype,
+                         &phase,
+                         &m,
+                         reordered_csr_nzval,
+                         csr_rowptr,
+                         reordered_csr_colind,
+                         &idum,
+                         &nrhs,
+                         iparm,
+                         &msglvl,
+                         &ddum,
+                         &ddum,
+                         &comm_,
+                         &error);
+
+   MFEM_ASSERT(error == 0, "Pardiso free error");
+
+   delete[] csr_rowptr;
+   delete[] reordered_csr_colind;
+   delete[] reordered_csr_nzval;
+}
+
+} // namespace mfem
+
+#endif // MFEM_USE_MKL_CPARDISO
+#endif // MFEM_USE_MPI

linalg/cpardiso.hpp

@@ -0,0 +1,125 @@
+// 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_CPARDISO
+#define MFEM_CPARDISO
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MKL_CPARDISO
+#ifdef MFEM_USE_MPI
+#include "operator.hpp"
+#include <mpi.h>
+#include "mkl_cluster_sparse_solver.h"
+
+namespace mfem
+{
+/**
+ * @brief MKL Parallel Direct Sparse Solver for Clusters
+ *
+ * Interface to the MPI enabled MKL version of Pardiso
+ */
+class CPardisoSolver : public Solver
+{
+public:
+   enum MatType
+   {
+      REAL_STRUCTURE_SYMMETRIC = 1,
+      REAL_UNSYMMETRIC = 11
+   };
+
+   /**
+    * @brief Construct a new CPardisoSolver object
+    *
+    * @param comm MPI Communicator
+    */
+   CPardisoSolver(MPI_Comm comm);
+
+   /**
+    * @brief Set the Operator object and perform factorization
+    *
+    * @a op needs to be of type HypreParMatrix. The contents are copied and
+    * reordered in an internal CSR structure.
+    *
+    * @param op Operator to use in factorization and solve
+    */
+   void SetOperator(const Operator &op) override;
+
+   /**
+    * @brief Solve
+    *
+    * @param b RHS vector
+    * @param x Solution vector
+    */
+   void Mult(const Vector &b, Vector &x) const override;
+
+   /**
+    * @brief Set the print level for Pardiso
+    *
+    * Prints statistics after the factorization and after each solve.
+    *
+    * @param print_lvl Print level
+    */
+   void SetPrintLevel(int print_lvl);
+
+   /**
+    * @brief Set the matrix type
+    *
+    * The matrix type supported is either real and symmetric or real and
+    * non-symmetric.
+    *
+    * @param mat_type Matrix type
+    */
+   void SetMatrixType(MatType mat_type);
+
+   ~CPardisoSolver();
+
+private:
+   MPI_Comm comm_;
+
+   // Global number of rows
+   int m;
+
+   // First row index of the global matrix on the local MPI rank
+   int first_row;
+
+   // Local number of nonzero entries
+   int nnz_loc;
+
+   // Local number of rows, obtained from a ParCSR matrix
+   int m_loc;
+
+   // CSR data structure for the copy data of the local CSR matrix
+   int *csr_rowptr = nullptr;
+   double *reordered_csr_nzval = nullptr;
+   int *reordered_csr_colind = nullptr;
+
+   // Internal solver memory pointer pt,
+   // 32-bit: int pt[64]
+   // 64-bit: long int pt[64] or void *pt[64] should be OK on both architectures
+   mutable void *pt[64] = {0};
+
+   // Solver control parameters, detailed description can be found in the
+   // constructor.
+   mutable int iparm[64] = {0};
+   mutable int maxfct, mnum, msglvl, phase, error;
+   int mtype;
+   int nrhs;
+
+   // Dummy variables
+   mutable int idum;
+   mutable double ddum;
+};
+} // namespace mfem
+
+#endif
+#endif // MFEM_USE_MKL_CPARDISO
+#endif // MFEM_USE_MPI

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
@@ -61,6 +65,10 @@
 #include "strumpack.hpp"
 #endif
 
+#ifdef MFEM_USE_MKL_CPARDISO
+#include "cpardiso.hpp"
+#endif
+
 #endif // MFEM_USE_MPI
 
 #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)
-   PKGS_NEED_MPI = SUPERLU STRUMPACK PETSC PUMI SLEPC
+<<<<<<< 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,9 +263,9 @@ 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
+ OCCA CEED RAJA UMPIRE MKL_CPARDISO
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
 
@@ -321,11 +325,11 @@ 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\
- MFEM_USE_ADIOS2 MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_SOURCE_DIR MFEM_INSTALL_DIR
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -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))
@@ -647,6 +652,7 @@ status info:
 	$(info MFEM_USE_UMPIRE        = $(MFEM_USE_UMPIRE))
 	$(info MFEM_USE_SIMD          = $(MFEM_USE_SIMD))
 	$(info MFEM_USE_ADIOS2        = $(MFEM_USE_ADIOS2))
+	$(info MFEM_USE_MKL_CPARDISO  = $(MFEM_USE_MKL_CPARDISO))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_HOST_CXX          = $(value MFEM_HOST_CXX))
 	$(info MFEM_CPPFLAGS          = $(value MFEM_CPPFLAGS))