Merge pull request #1510 from norihiro-w/eigen-pardiso

support Pardiso in Eigen

CMakeLists.txt

@@ -84,6 +84,9 @@ option(OGS_USE_EIGEN "Use Eigen linear solver" ON)
 option(OGS_EIGEN_DYNAMIC_SHAPE_MATRICES "Use dynamically allocated shape matrices" ON)
 option(EIGEN_NO_DEBUG "Disables Eigen's assertions" OFF)
 
+# MKL
+option(OGS_USE_MKL "Use Intel MKL" OFF)
+
 # Logging
 set(OGS_LOG_LEVEL "LOGOG_LEVEL_DEBUG" CACHE STRING "Set logging level included in compilation.")
 set_property(CACHE OGS_LOG_LEVEL PROPERTY STRINGS
@@ -143,6 +146,11 @@ if(OGS_USE_LIS)
     set(OGS_USE_EIGEN ON)
 endif()
 
+if(OGS_USE_MKL)
+    add_definitions(-DUSE_MKL)
+    include_directories(SYSTEM ${MKL_INCLUDE_DIR})
+endif()
+
 if(OGS_USE_PETSC)
     add_definitions(-DUSE_PETSC)
     set(OGS_USE_MPI ON CACHE BOOL "Use MPI" FORCE)

MathLib/CMakeLists.txt

@@ -43,6 +43,10 @@ if (OGS_USE_LIS)
     target_link_libraries(MathLib ${LIS_LIBRARIES})
 endif()
 
+if (OGS_USE_MKL)
+    target_link_libraries(MathLib ${MKL_LIBRARIES})
+endif()
+
 if (OGS_USE_PETSC)
     target_link_libraries(MathLib ${PETSC_LIBRARIES})
 endif()

MathLib/LinAlg/Eigen/EigenLinearSolver.cpp

@@ -11,6 +11,10 @@
 
 #include <logog/include/logog.hpp>
 
+#ifdef USE_MKL
+#include <Eigen/PardisoSupport>
+#endif
+
 #include "BaseLib/ConfigTree.h"
 #include "EigenVector.h"
 #include "EigenMatrix.h"
@@ -174,6 +178,17 @@ EigenLinearSolver::EigenLinearSolver(
             _solver = details::createIterativeSolver(_option.solver_type,
                                                      _option.precon_type);
             return;
+        case EigenOption::SolverType::PardisoLU: {
+#ifdef USE_MKL
+            using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
+            _solver.reset(new details::EigenDirectLinearSolver<SolverType>);
+            return;
+#else
+            OGS_FATAL(
+                "The code is not compiled with Intel MKL. Linear solver type "
+                "PardisoLU is not available.");
+#endif
+        }
     }
 
     OGS_FATAL("Invalid Eigen linear solver type. Aborting.");

MathLib/LinAlg/Eigen/EigenOption.cpp

@@ -29,6 +29,8 @@ EigenOption::SolverType EigenOption::getSolverType(const std::string &solver_nam
         return SolverType::BiCGSTAB;
     if (solver_name == "SparseLU")
         return SolverType::SparseLU;
+    if (solver_name == "PardisoLU")
+        return SolverType::PardisoLU;
 
     OGS_FATAL("Unknown Eigen solver type `%s'", solver_name.c_str());
 }

MathLib/LinAlg/Eigen/EigenOption.h

@@ -23,7 +23,8 @@ struct EigenOption final
     {
         CG,
         BiCGSTAB,
-        SparseLU
+        SparseLU,
+        PardisoLU
     };
 
     /// Preconditioner type

scripts/cmake/Find.cmake

@@ -115,6 +115,10 @@ if(OGS_USE_LIS)
     find_package( LIS REQUIRED )
 endif()
 
+if(OGS_USE_MKL)
+    find_package( MKL REQUIRED )
+endif()
+
 if(OGS_USE_PETSC)
     message(STATUS "Configuring for PETSc")
 

scripts/cmake/cmake/FindMKL.cmake

+# Find Intel Math Karnel Library (MKL)
+#
+# Options
+# - MKL_DIR      MKL root directory
+# - MKL_OPENMP   use OpenMP threading
+#
+# Results
+# - MKL_INCLUDE_DIR
+# - MKL_LIBRARIES
+#
+# Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+# Distributed under a Modified BSD License.
+# See accompanying file LICENSE.txt or
+# http://www.opengeosys.org/project/license
+
+# Lookg for MKL root dir
+if (NOT MKL_DIR)
+    find_path(MKL_DIR
+        include/mkl.h
+        PATHS
+        $ENV{MKL_DIR}
+        /opt/intel/mkl/
+        )
+endif()
+message("MKL_DIR : ${MKL_DIR}")
+
+# Find MKL include dir
+find_path(MKL_INCLUDE_DIR NAMES mkl.h
+    PATHS
+    ${MKL_DIR}/include
+)
+
+# Set the directory path storing MKL libraries
+if (NOT MKL_LIB_DIR)
+    if(APPLE)
+        set(MKL_LIB_DIR ${MKL_DIR}/lib)
+    else()
+        if (${CMAKE_HOST_SYSTEM_PROCESSOR} STREQUAL "x86_64")
+            set(MKL_LIB_DIR ${MKL_DIR}/lib/intel64)
+        else()
+            set(MKL_LIB_DIR ${MKL_DIR}/lib/ia32)
+        endif()
+    endif()
+endif()
+
+# Find MKL libs
+find_library(MKL_LIB_CORE mkl_core PATHS ${MKL_LIB_DIR})
+
+if (${CMAKE_HOST_SYSTEM_PROCESSOR} STREQUAL "x86_64")
+    set(MKL_INTEL_LIB_NAME mkl_intel_lp64)
+else()
+    set(MKL_INTEL_LIB_NAME mkl_intel)
+endif()
+
+find_library(MKL_LIB_INTEL ${MKL_INTEL_LIB_NAME} PATHS ${MKL_LIB_DIR})
+
+if(OPENMP_FOUND)
+    set(MKL_THREAD_LIB_NAME "mkl_gnu_thread")
+else()
+    set(MKL_THREAD_LIB_NAME "mkl_sequential")
+endif()
+find_library(MKL_LIB_THREAD ${MKL_THREAD_LIB_NAME} PATHS ${MKL_LIB_DIR})
+
+
+set(MKL_LIBRARIES "${MKL_LIB_INTEL}" "${MKL_LIB_THREAD}" "${MKL_LIB_CORE}")
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(MKL DEFAULT_MSG MKL_INCLUDE_DIR MKL_LIBRARIES)
+mark_as_advanced(MKL_INCLUDE_DIR MKL_LIBRARIES)