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

@@ -48,6 +48,9 @@
 #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

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

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

@@ -141,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 =
@@ -378,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
 

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

@@ -25,6 +25,8 @@ int main(int argc, char *argv[])
    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",
@@ -32,6 +34,10 @@ int main(int argc, char *argv[])
    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",
@@ -63,8 +69,7 @@ int main(int argc, char *argv[])
    //    'ref_levels' to be the largest number that gives a final mesh with no
    //    more than 10,000 elements.
    {
-      int ref_levels = 2;
-      for (int l = 0; l < ref_levels; l++)
+      for (int l = 0; l < sr; l++)
       {
          mesh.UniformRefinement();
       }
@@ -76,8 +81,7 @@ int main(int argc, char *argv[])
    ParMesh pmesh(MPI_COMM_WORLD, mesh);
    mesh.Clear();
    {
-      int par_ref_levels = 2;
-      for (int l = 0; l < par_ref_levels; l++)
+      for (int l = 0; l < pr; l++)
       {
          pmesh.UniformRefinement();
       }
@@ -136,36 +140,62 @@ int main(int argc, char *argv[])
    // // 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(1);
+   cg.SetPrintLevel(0);
    if (prec) { cg.SetPreconditioner(*prec); }
    cg.SetOperator(A);
    cg.Mult(B, X);
    delete prec;
 
-   StopWatch chrono;
+   if (myid == 0)
+   {
+      cout << "PCG-AMG time: " << chrono.RealTime() << endl;
+   }
+
    chrono.Clear();
    chrono.Start();
-
    {
       MUMPSSolver MA;
-      MA.SetMatrixSymType(1);
+      MA.SetMatrixSymType(0);
       MA.SetOperator(A);
       Vector Y(X.Size());
       MA.Mult(B,Y);
       Y-=X;
       cout << "Mumps Diff norm = " << Y.Norml2() << endl;
    }
-   cout << "mumps time: " << chrono.RealTime() << endl;
+   if (myid == 0)
+   {
+      cout << "mumps time: " << chrono.RealTime() << endl;
+   }
 
    chrono.Clear();
    chrono.Start();
 
    {
-      // Transform to monolithic HypreParMatrix
+      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);
@@ -178,8 +208,10 @@ int main(int argc, char *argv[])
       cout << "Superlu Diff norm = " << Y.Norml2() << endl;
    }
 
-   cout << "superlu time: " << chrono.RealTime() << endl;
-
+   if (myid == 0)
+   {
+      cout << "superlu time: " << chrono.RealTime() << endl;
+   }
 
    a.RecoverFEMSolution(X, b, x);
 

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

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

@@ -65,7 +65,9 @@
 #include "strumpack.hpp"
 #endif
 
-
+#ifdef MFEM_USE_MKL_CPARDISO
+#include "cpardiso.hpp"
+#endif
 
 #endif // MFEM_USE_MPI
 

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),))
@@ -261,7 +265,7 @@ endif
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 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))
 
@@ -325,7 +329,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  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\
@@ -648,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))