Pardiso support on Kingspeak

examples/EMSchur3D-vtk.cpp

@@ -117,8 +117,8 @@ int main( int argc, char** argv ) {
     //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::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
-    //auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Symmetric > >::NewSP();
 
     // LS CG
     //auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();

examples/EMSchur3D.cpp

@@ -36,7 +36,7 @@ int main( int argc, char** argv ) {
     //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();
+    auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Symmetric&Eigen::Lower > >::NewSP();
 
     // Eigen built-in
     //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();

include/EMSchur3D.h

@@ -12,10 +12,9 @@
  * @date      02/19/2015 01:10:39 PM
  * @version   $Id$
  * @author    Trevor Irons (ti)
- * @email     Trevor.Irons@xri-geo.com
+ * @email     Trevor.Irons@Lemmasoftware.org
  * @copyright Copyright (c) 2015, XRI Geophysics, LLC
- * @copyright Copyright (c) 2015, Trevor Irons
- * @copyright Copyright (c) 2011, Trevor Irons
+ * @copyright Copyright (c) 2011, 2015, 2017, 2018  Trevor Irons
  * @copyright Copyright (c) 2011, Colorado School of Mines
  */
 
@@ -247,7 +246,13 @@ namespace Lemma {
 
         VectorXcr ADiv = D*A;  // ADiv == RHS == D C^I Se
         //VectorXcr Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() - Se.array());
+#if LOWER==1
+        std::cout << "Using Eigen::Lower to calculate Error" << std::endl;
         VectorXcr Error = ((Cvec[iw].selfadjointView<Eigen::Lower>()*A).array() - Se.array());
+#elif UPPER==1
+        std::cout << "Using Eigen::Upper to calculate Error" << std::endl;
+        VectorXcr Error = ((Cvec[iw].selfadjointView<Eigen::Upper>()*A).array() - Se.array());
+#endif
         logio << "|| Div(A) || = " << ADiv.norm()
               << "\tInital solution error="<<   Error.norm()  // Iteritive info
 //              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
@@ -306,8 +311,13 @@ namespace Lemma {
 
         // Report error of solutions
         //Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
+#if LOWER==1
+        std::cout << "Using Eigen::Lower to calculate Error" << std::endl;
         Error = ((Cvec[iw].selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
-
+#elif UPPER==1
+        std::cout << "Using Eigen::Upper to calculate Error" << std::endl;
+        Error = ((Cvec[iw].selfadjointView<Eigen::Upper>()*A).array() + E.array() - Se.array());
+#endif
         //      << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
         //      << "\ttime " << timer.end() / 60. << " [m]   "
         //      <<  CSolver[iw].iterations() << "  iterations"
@@ -466,13 +476,14 @@ namespace Lemma {
     }
 
     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()];
+    void EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Symmetric > >::BuildCDirectSolver() {
+        CSolver = new Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Symmetric > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
             timer.begin();
         
             //CSolver[iw].pardisoParameterArray()[2] = 28; // OMP_NUM_THREADS?
+            //Eigen::SparseMatrix<Complex> Csym = Cvec[iw].selfadjointView< Eigen::Lower >();
             
             /*  Complex system */
             std::cout << "PardisoLDLT pattern analyzing C_" << iw << ",";

include/EMSchur3DBase.h

@@ -77,6 +77,9 @@
 
 namespace Lemma {
 
+#define UPPER 0  // LOWER WAS 0
+#define LOWER 1  // 1=true, 0=false
+
 enum SOLVER{ SPARSELU, SimplicialLLT, SimplicialLDLT, BiCGStab, SparseQR };
 
 /**

src/EMSchur3DBase.cpp

@@ -26,8 +26,9 @@
 typedef Eigen::Triplet<Lemma::Complex> Tc;
 typedef Eigen::Triplet<Lemma::Real> Tr;
 
-#define UPPER 0  // LOWER WAS 0
-#define LOWER 1  // 1=true, 0=false
+// Moved to header 
+//#define UPPER 0  // LOWER WAS 0
+//#define LOWER 1  // 1=true, 0=false
 
 namespace Lemma {