Pardiso on kingspeak

CMakeLists.txt

@@ -48,6 +48,12 @@ if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT OR LEMMA_VTK8_SUPPORT AND LEMMA_MO
     	CXX_STANDARD_REQUIRED_ON
 	)
 
+        # No OpenMP
+	#set (MKLLINK  "-L$MKLROOT/lib/intel64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -Wl,-rpath=$MKLROOT/lib/intel64")
+        # OpenMP 
+        set (MKLLINK  "-L$MKLROOT/lib/intel64 -Wl,-rpath=$MKLROOT/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread")
+	target_link_libraries(emschur3d "lemmacore" ${MKLLINK})
+
 	# Linking
 	target_link_libraries(emschur3d "lemmacore" "fdem1d" )
 

examples/EMSchur3D-vtk.cpp

@@ -115,7 +115,10 @@ int main( int argc, char** argv ) {
 
     // Use BiCGSTAB Diagonal preconditioner
     //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
-    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    
+    auto EM3D = EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > >::NewSP();
 
     // LS CG
     //auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();

examples/EMSchur3D.cpp

@@ -33,7 +33,10 @@ int main( int argc, char** argv ) {
     //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
 
     // SUPERLU
-    auto EM3D = EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    
+    //auto EM3D = EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > >::NewSP();
 
     // Eigen built-in
     //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();

include/EMSchur3D.h

@@ -37,6 +37,8 @@
 #include <Eigen/PaStiXSupport>
 #endif
 
+#include <Eigen/PardisoSupport>
+
 //#include "CSymSimplicialCholesky.h"
 
 namespace Lemma {
@@ -437,6 +439,57 @@ namespace Lemma {
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
         }
     }
+    
+// Pardiso 
+    template<>
+    void EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::BuildCDirectSolver() {
+        CSolver = new Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            jsw_timer timer;
+            timer.begin();
+        
+            //CSolver[iw].pardisoParameterArray()[2] = 28; // OMP_NUM_THREADS?
+            
+            /*  Complex system */
+            std::cout << "PardisoLU pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+
+            // factorize
+            timer.begin();
+            std::cout << "PardisoLU factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
+
+    template<>
+    void EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > >::BuildCDirectSolver() {
+        CSolver = new Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            jsw_timer timer;
+            timer.begin();
+        
+            //CSolver[iw].pardisoParameterArray()[2] = 28; // OMP_NUM_THREADS?
+            
+            /*  Complex system */
+            std::cout << "PardisoLDLT pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            //CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            CSolver[iw].analyzePattern( Cvec[iw] );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+
+            // factorize
+            timer.begin();
+            std::cout << "PardisoLDLT factorising C_" << iw << ", ";
+            std::cout.flush();
+            //CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            CSolver[iw].factorize( Cvec[iw] );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
 
 #ifdef HAVE_UMFPACK
     // Umfpack only seems to work when LOWER and UPPER are set to 1. Workarounds this have not been found.