Merge branch '71-load-mkl-pardiso-at-run-time-if-available' into 'mas…

…ter' Resolve "Load MKL-Pardiso at run-time if available" Closes #71 See merge request jschoeberl/ngsolve!262
NGSolve · Nov 9, 2017 · 45aff082b07605cb9ad479c02d282e54b6fa29bd · 45aff08
2 parents fa40e7c + f940cd2
commit 45aff082b07605cb9ad479c02d282e54b6fa29bd
diff --git a/fem/code_generation.cpp b/fem/code_generation.cpp
@@ -24,20 +24,9 @@ namespace ngfem
        return name;
    }

-    void Library::Load( string alib_name )
-    {
-      lib_name = alib_name;
-#ifdef WIN32
-      lib = LoadLibrary(lib_name.c_str());
-      if (!lib) throw std::runtime_error(string("Could not load library ") + lib_name);
-#else // WIN32
-      lib = dlopen(lib_name.c_str(), RTLD_NOW);
-      if(lib == nullptr) throw std::runtime_error(dlerror());
-#endif // WIN32
-    }
-
-    void Library::Compile(const std::vector<string> &codes, const std::vector<string> &link_flags )
+    unique_ptr<SharedLibrary> CompileCode(const std::vector<string> &codes, const std::vector<string> &link_flags )
    {
+      static int counter = 0;
      static ngstd::Timer tcompile("CompiledCF::Compile");
      static ngstd::Timer tlink("CompiledCF::Link");
      string object_files;
@@ -75,41 +64,13 @@ namespace ngfem
      if (err) throw Exception ("problem calling linker");      
      tlink.Stop();
      cout << IM(3) << "done" << endl;
+      auto library = make_unique<SharedLibrary>();
 #ifdef WIN32
-      Load(prefix+".dll");
+      library->Load(prefix+".dll");
 #else
-      Load("./"+prefix+".so");
+      library->Load("./"+prefix+".so");
 #endif
+      return library;
    }

-    void* Library::GetRawFunction( string func_name )
-    {
-#ifdef WIN32
-      void* func = GetProcAddress(lib, func_name.c_str());
-      if(func == nullptr)
-        throw std::runtime_error(string("Could not find function ") + func_name + " in library " + lib_name);
-
-#else // WIN32
-      void* func = dlsym(lib, func_name.c_str());
-      if(func == nullptr)
-          throw std::runtime_error(dlerror());
-
-#endif // WIN32
-
-      return func;
-    }
-
-    Library::~Library()
-    {
-      if(lib)
-      {
-
-#ifdef WIN32
-        FreeLibrary(lib);
-#else // WIN32
-        int rc = dlclose(lib);
-        if(rc != 0) cerr << "Failed to close library " << lib_name << endl;
-#endif // WIN32
-      }
-    }
 }
diff --git a/fem/code_generation.hpp b/fem/code_generation.hpp
@@ -7,12 +7,6 @@
 /* Date:   13. Apr. 2016                                             */
 /*********************************************************************/

-#ifdef WIN32
-#include <windows.h>
-#else // WIN32
-#include <dlfcn.h>
-#endif //WIN32
-
 #include <map>

 namespace ngfem
@@ -177,34 +171,7 @@ namespace ngfem
    }
  }

-  class Library
-  {
-    static int counter; // defined in coefficient.cpp
-    string lib_name;
-
-#ifdef WIN32
-    HINSTANCE lib;
-#else // WIN32
-    void *lib;
-#endif // WIN32
-
-    void *GetRawFunction( string func_name );
-
-  public:
-    Library() : lib(nullptr) {}
-    // Compile a given string and load the library
-    void Compile(const std::vector<string> &codes, const std::vector<string> &libraries );
-
-    void Load( string alib_name );
+  unique_ptr<SharedLibrary> CompileCode(const std::vector<string> &codes, const std::vector<string> &libraries );

-    template <typename TFunc>
-    TFunc GetFunction( string func_name )
-    {
-      return reinterpret_cast<TFunc>(GetRawFunction(func_name));
-    }
-
-    ~Library();
-
-  };
 }
 #endif // FILE_NGS_CODE_GENERATION___
diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp
@@ -4686,7 +4686,7 @@ shared_ptr<CoefficientFunction> MakeCoordinateCoefficientFunction (int comp)
    int totdim;
    Array<bool> is_complex;
    // Array<Timer*> timers;
-    Library library;
+    unique_ptr<SharedLibrary> library;
    lib_function compiled_function = nullptr;
    lib_function_simd compiled_function_simd = nullptr;
    lib_function_deriv compiled_function_deriv = nullptr;
@@ -4843,25 +4843,25 @@ shared_ptr<CoefficientFunction> MakeCoordinateCoefficientFunction (int comp)
        }

        auto compile_func = [this, codes, link_flags, maxderiv] () {
-              library.Compile( codes, link_flags );
+              library = CompileCode( codes, link_flags );
              if(cf->IsComplex())
              {
-                  compiled_function_simd_complex = library.GetFunction<lib_function_simd_complex>("CompiledEvaluateSIMD");
-                  compiled_function_complex = library.GetFunction<lib_function_complex>("CompiledEvaluate");
+                  compiled_function_simd_complex = library->GetFunction<lib_function_simd_complex>("CompiledEvaluateSIMD");
+                  compiled_function_complex = library->GetFunction<lib_function_complex>("CompiledEvaluate");
              }
              else
              {
-                  compiled_function_simd = library.GetFunction<lib_function_simd>("CompiledEvaluateSIMD");
-                  compiled_function = library.GetFunction<lib_function>("CompiledEvaluate");
+                  compiled_function_simd = library->GetFunction<lib_function_simd>("CompiledEvaluateSIMD");
+                  compiled_function = library->GetFunction<lib_function>("CompiledEvaluate");
                  if(maxderiv>0)
                  {
-                      compiled_function_simd_deriv = library.GetFunction<lib_function_simd_deriv>("CompiledEvaluateDerivSIMD");
-                      compiled_function_deriv = library.GetFunction<lib_function_deriv>("CompiledEvaluateDeriv");
+                      compiled_function_simd_deriv = library->GetFunction<lib_function_simd_deriv>("CompiledEvaluateDerivSIMD");
+                      compiled_function_deriv = library->GetFunction<lib_function_deriv>("CompiledEvaluateDeriv");
                  }
                  if(maxderiv>1)
                  {
-                      compiled_function_simd_dderiv = library.GetFunction<lib_function_simd_dderiv>("CompiledEvaluateDDerivSIMD");
-                      compiled_function_dderiv = library.GetFunction<lib_function_dderiv>("CompiledEvaluateDDeriv");
+                      compiled_function_simd_dderiv = library->GetFunction<lib_function_simd_dderiv>("CompiledEvaluateDDerivSIMD");
+                      compiled_function_dderiv = library->GetFunction<lib_function_dderiv>("CompiledEvaluateDDeriv");
                  }
              }
              cout << IM(7) << "Compilation done" << endl;
@@ -5289,8 +5289,5 @@ shared_ptr<CoefficientFunction> MakeCoordinateCoefficientFunction (int comp)
    return make_shared<CompiledCoefficientFunction> (c, realcompile, maxderiv, wait);
  }
  
-
-  int Library::counter = 0;
- 
 }

diff --git a/linalg/la.hpp b/linalg/la.hpp
@@ -35,8 +35,8 @@ namespace ngla
 #include "sparsematrix.hpp"
 #include "order.hpp"
 #include "sparsecholesky.hpp"
+#include "pardisoinverse.hpp"
 // include these only from c++-files
-// #include "pardisoinverse.hpp"
 // #include "umfpackinverse.hpp"
 // #include "superluinverse.hpp"
 // #include "mumpsinverse.hpp"
diff --git a/linalg/pardisoinverse.cpp b/linalg/pardisoinverse.cpp
@@ -7,7 +7,6 @@
 #include <la.hpp>
 #include "pardisoinverse.hpp"

-#ifdef USE_PARDISO

 // #include "/opt/intel/mkl/include/mkl_service.h"

@@ -19,29 +18,45 @@ extern "C"
 {
  /* PARDISO prototype. */

+#ifdef USE_PARDISO
+
 #ifdef USE_MKL
+#define MKL_PARDISO

  void mkl_free_buffers (void);
  void F77_FUNC(pardiso)
    (void * pt, integer * maxfct, integer * mnum, integer * mtype, integer * phase, integer * n, 
     double * a, integer * ia, integer * ja, integer * perm, integer * nrhs, integer * iparam, 
     integer * msglvl, double * b, double * x, integer * error);

-#else
+#else // USE_MKL

 #ifdef USE_PARDISO400
 extern  integer F77_FUNC(pardisoinit)
    (void *, integer *, integer *, integer *, double *, integer *);
-#else
+#else // USE_PARDISO400
 extern  integer F77_FUNC(pardisoinit)
    (void *, integer *, integer *);
-#endif
+#endif // USE_PARDISO400
  integer F77_FUNC(pardiso)
    (void * pt, integer * maxfct, integer * mnum, integer * mtype, integer * phase, integer * n, 
     double * a, integer * ia, integer * ja, integer * perm, integer * nrhs, integer * iparam, 
     integer * msglvl, double * b, double * x, integer * error);

 #endif
+
+#else // USE_PARDISO
+  // Neither MKL nor PARDISO linked at compile-time
+  // check for MKL at run-time and set function pointers if available
+#define MKL_PARDISO
+
+  void (*mkl_free_buffers) (void) = nullptr;
+  void (*F77_FUNC(pardiso))
+    (void * pt, integer * maxfct, integer * mnum, integer * mtype, integer * phase, integer * n,
+     double * a, integer * ia, integer * ja, integer * perm, integer * nrhs, integer * iparam,
+     integer * msglvl, double * b, double * x, integer * error) = nullptr;
+
+#endif // USE_PARDISO
 }


@@ -51,6 +66,31 @@ extern  integer F77_FUNC(pardisoinit)

 namespace ngla
 {
+#ifdef USE_PARDISO
+  bool is_pardiso_available = true;
+#else
+  static SharedLibrary libmkl;
+  static bool LoadMKL()
+  {
+      try
+      {
+#ifdef WIN32
+          libmkl.Load("mkl_rt.dll");
+#else
+          libmkl.Load("libmkl_rt.so");
+#endif
+          mkl_free_buffers = libmkl.GetFunction<decltype(mkl_free_buffers)>("mkl_free_buffers");
+          F77_FUNC(pardiso) = libmkl.GetFunction<decltype(pardiso_)>("pardiso_");
+          return mkl_free_buffers && F77_FUNC(pardiso);
+      }
+      catch(const std::runtime_error &)
+      {
+          return false;
+      }
+  };
+  bool is_pardiso_available = LoadMKL();
+#endif
+
  int pardiso_msg = 0;


@@ -162,7 +202,7 @@ namespace ngla

    for (int i = 0; i < 128; i++) pt[i] = 0;

-#ifdef USE_MKL
+#ifdef MKL_PARDISO
    //    no init in MKL PARDISO
 #else

@@ -555,9 +595,9 @@ namespace ngla
    F77_FUNC(pardiso) ( pt, &maxfct, &mnum, &matrixtype, &phase, &compressed_height, NULL,
 			&rowstart[0], &indices[0], NULL, &nrhs, params, &msglevel,
 			NULL, NULL, &error );
-#ifdef USE_MKL
+#ifdef MKL_PARDISO
    mkl_free_buffers();
-#endif // USE_MKL
+#endif // MKL_PARDISO
    if (task_manager) task_manager -> StartWorkers();
    if (error != 0)
      cout << "Clean Up: PARDISO returned error " << error << "!" << endl;
@@ -721,6 +761,3 @@ namespace ngla

 }

-
-
-#endif
diff --git a/linalg/pardisoinverse.hpp b/linalg/pardisoinverse.hpp
@@ -19,6 +19,7 @@
 namespace ngla
 {
  using ngbla::integer;
+  extern bool is_pardiso_available;


  /*
@@ -188,4 +189,4 @@ namespace ngla

 }

-#endif
+#endif
diff --git a/linalg/python_linalg.cpp b/linalg/python_linalg.cpp
@@ -417,8 +417,21 @@ void NGS_DLL_HEADER ExportNgla(py::module &m) {
                                       if (inverse != "") m.SetInverseType(inverse);
                                       return m.InverseMatrix(freedofs);
                                     }
-         ,"Inverse", py::arg("freedofs")=nullptr, py::arg("inverse")=py::str("")
-         )
+         ,"Inverse", py::arg("freedofs")=nullptr, py::arg("inverse")=py::str(""), 
+         docu_string(R"raw_string(Calculate inverse of sparse matrix"
+Parameters
+
+freedofs : BitArray
+  If set, invert only the rows/columns the matrix defined by the bit array, otherwise invert the whole matrix
+
+inverse : string
+  Solver to use, allowed values are:
+    sparsecholesky - internal solver of NGSolve for symmetric matrices
+    umfpack        - solver by Suitesparse/UMFPACK (if NGSolve was configured with USE_UMFPACK=ON)
+    pardiso        - PARDISO, either provided by libpardiso (USE_PARDISO=ON) or Intel MKL (USE_MKL=ON).
+                     If neither Pardiso nor Intel MKL was linked at compile-time, NGSolve will look
+                     for libmkl_rt in LD_LIBRARY_PATH (Unix) or PATH (Windows) at run-time.
+)raw_string"))
    // .def("Inverse", [](BM &m)  { return m.InverseMatrix(); })
    .def("Update", [](BM &m) { m.Update(); })
    ;