MACRO to switch ColMajor vs RowMajor

CMakeLists.txt

 	INCLUDE_DIRECTORIES(${SUPERLU_INCLUDES})
 endif()
 
-find_package( PaSTiX )
+#find_package( PaSTiX )
+find_package( PaSTiX COMPONENTS METIS SEQ )
+
+find_path(PASTIX_pastix_nompi.h_INCLUDE_DIRS
+  #NAMES pastix/nompi.h
+  NAMES pastix_nompi.h
+  HINTS ${PASTIX_INCLUDE_DIRS}
+)
+if (NOT PASTIX_pastix_nompi.h_INCLUDE_DIRS)
+  message(STATUS "A version of Pastix has been found but pastix_nompi.h does not exist in the include directory."
+                 " Because Eigen tests require a version without MPI, we disable the Pastix backend.")
+endif()
+if(PASTIX_FOUND AND PASTIX_pastix_nompi.h_INCLUDE_DIRS)
+  add_definitions("-DEIGEN_PASTIX_SUPPORT")
+  include_directories(${PASTIX_INCLUDE_DIRS_DEP})
+  if(SCOTCH_FOUND)
+    include_directories(${SCOTCH_INCLUDE_DIRS})
+    set(PASTIX_LIBRARIES ${PASTIX_LIBRARIES} ${SCOTCH_LIBRARIES})
+  elseif(METIS_FOUND)
+    include_directories(${METIS_INCLUDE_DIRS})
+    set(PASTIX_LIBRARIES ${PASTIX_LIBRARIES} ${METIS_LIBRARIES})
+  else(SCOTCH_FOUND)
+#    ei_add_property(EIGEN_MISSING_BACKENDS  "PaStiX, ")
+  endif(SCOTCH_FOUND)
+  set(SPARSE_LIBS ${SPARSE_LIBS} ${PASTIX_LIBRARIES_DEP} ${ORDERING_LIBRARIES})
+  set(PASTIX_ALL_LIBS ${PASTIX_LIBRARIES_DEP})
+#  ei_add_property(EIGEN_TESTED_BACKENDS  "PaStiX, ")
+else()
+#  ei_add_property(EIGEN_MISSING_BACKENDS  "PaStiX, ")
+endif()
+
 if ( PASTIX_FOUND )
 	message( STATUS "PaSTiX was found" )
 	add_compile_options(-DHAVE_PASTIX)
 	INCLUDE_DIRECTORIES(${PASTIX_INCLUDES})
 endif()
 
+
+
 find_package( Umfpack )
 if ( UMFPACK_FOUND )
 	message( STATUS "Umfpack was found" )

examples/EMSchur3D-vtk.cpp

     // And solve
 
     // 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, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::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();
+    //auto EM3D = EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Symmetric > >::NewSP();
 
     // LS CG
-    //auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     // CG...dangerous
-    //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > >::NewSP();
-    //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > >::NewSP();
 
     // LU direct
-    //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::COLAMDOrdering<int> > >::NewSP();
 
     EM3D->SetRectilinearGrid( gimp->GetGrid() );
 

examples/EMSchur3D.cpp

 int main( int argc, char** argv ) {
 
     // BiCGSTAB Diagonal preconditioner
-    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > >::NewSP();
 
-    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     // SUPERLU
-    //auto EM3D = EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
-    //auto EM3D = EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
-    //auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Symmetric&Eigen::Lower > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Symmetric&Eigen::Lower > >::NewSP();
 
     // Eigen built-in
-    //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     // UmfPack
-    //auto EM3D = EMSchur3D< Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::NewSP();
+    auto EM3D = EMSchur3D< Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     // PaStiX
-    //auto EM3D = EMSchur3D< Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::PastixLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, true > >::NewSP();
 
 
     if (argc < 3) {

include/EMSchur3D.h

 
 #ifdef HAVE_SUPERLU
     template<>
-    void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::BuildCDirectSolver() {
-        CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
+    void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::BuildCDirectSolver() {
+        CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > [Omegas.size()];
 
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
 #endif
 
     template<>
-    void EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::COLAMDOrdering<int> > >::BuildCDirectSolver() {
-        CSolver = new Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::COLAMDOrdering<int> > [Omegas.size()];
+    void EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::COLAMDOrdering<int> > >::BuildCDirectSolver() {
+        CSolver = new Eigen::SparseLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::COLAMDOrdering<int> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
             timer.begin();
 
 #ifdef HAVE_PARDISO
     template<>
-    void EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::BuildCDirectSolver() {
-        CSolver = new Eigen::PardisoLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
+    void EMSchur3D< Eigen::PardisoLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::BuildCDirectSolver() {
+        CSolver = new Eigen::PardisoLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
             timer.begin();
     }
 
     template<>
-    void EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Symmetric > >::BuildCDirectSolver() {
-        CSolver = new Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Symmetric > [Omegas.size()];
+    void EMSchur3D< Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Symmetric > >::BuildCDirectSolver() {
+        CSolver = new Eigen::PardisoLDLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Symmetric > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
             timer.begin();
 #ifdef HAVE_UMFPACK
     // Umfpack only seems to work when LOWER and UPPER are set to 1. Workarounds this have not been found.
     template<>
-    void EMSchur3D< Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::BuildCDirectSolver() {
-        CSolver = new Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
+    void EMSchur3D< Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::BuildCDirectSolver() {
+        CSolver = new Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
             timer.begin();
 #endif
 
 //     template<>
-//     void EMSchur3D< Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
+//     void EMSchur3D< Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > [Omegas.size()];
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 //             jsw_timer timer;
 //     }
 
 //     template<>
-//     void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > [Omegas.size()];
+//     void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::NaturalOrdering<int> > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::NaturalOrdering<int> > [Omegas.size()];
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 //             jsw_timer timer;
 //     }
 //
 //     template<>
-//     void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
+//     void EMSchur3D< Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 //             jsw_timer timer;
 //     }
 //
 //     template<>
-//     void EMSchur3D< Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
+//     void EMSchur3D< Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 //             jsw_timer timer;
 
 
     template<>
-    void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
+    void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
             CSolver[iw].preconditioner().setDroptol(1e-6);       //1e-5);      // 1e-12
     }
 
     template<>
-    void EMSchur3D< Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
+    void EMSchur3D< Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, SPARSEMAJOR> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, SPARSEMAJOR> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
             jsw_timer timer;
     }
 
     template<>
-    void EMSchur3D< Eigen::LeastSquaresConjugateGradient< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::LeastSquaresConjugateGradient< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
+    void EMSchur3D< Eigen::LeastSquaresConjugateGradient< Eigen::SparseMatrix<Complex, SPARSEMAJOR> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::LeastSquaresConjugateGradient< Eigen::SparseMatrix<Complex, SPARSEMAJOR> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
             jsw_timer timer;
     }
 
     template<>
-    void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
+    void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
             jsw_timer timer;
     }
 
     template<>
-    void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > [Omegas.size()];
+    void EMSchur3D<   Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
             jsw_timer timer;
     }
 
 //     template<>
-//     void EMSchur3D<   Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
+//     void EMSchur3D<   Eigen::PastixLLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > [Omegas.size()];
 //         //MPI_Init(NULL, NULL);
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 //     }
 //
 //     template<>
-//     void EMSchur3D<   Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::Lower > [Omegas.size()];
+//     void EMSchur3D<   Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::Lower > [Omegas.size()];
 //         //MPI_Init(NULL, NULL);
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 
 #ifdef HAVE_PASTIX
     template<>
-    void EMSchur3D<   Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > [Omegas.size()];
+    void EMSchur3D<   Eigen::PastixLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, true > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, true > [Omegas.size()];
 
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;

include/EMSchur3DBase.h

 
 #include "timer.h"
 #include <Eigen/Sparse>
+
+//#define  SPARSEMAJOR Eigen::ColMajor
+#define  SPARSEMAJOR Eigen::RowMajor
+
+// SPARSEMAJOR;
+//typedef  Eigen::RowMajor  SPARSEMAJOR;
+
 //#include "bicgstab.h"
 
 // // Solvers
 
     /** Matrix objects representing discrete operators
      */
-    Eigen::SparseMatrix<Complex, Eigen::RowMajor>*                  Cvec;
-    Eigen::SparseMatrix<Complex, Eigen::RowMajor>                   D;
+    Eigen::SparseMatrix<Complex, SPARSEMAJOR>*                  Cvec;
+    Eigen::SparseMatrix<Complex, SPARSEMAJOR>                   D;
 
     /** Squeezed matrices for reduced phi grid
      */
-    Eigen::SparseMatrix<Complex, Eigen::RowMajor>*                  Cvec_s;
-    Eigen::SparseMatrix<Complex, Eigen::RowMajor>                   D_s;
+    Eigen::SparseMatrix<Complex, SPARSEMAJOR>*                  Cvec_s;
+    Eigen::SparseMatrix<Complex, SPARSEMAJOR>                   D_s;
 
     /** number of cells in x, set when RectilinearGrid is attached */
     int nx;
     std::vector<int> idx;
 
     /** Dirichlet on low X */
-    bool DirichletXLO = true;
-    //bool DirichletXLO = false;
+    //bool DirichletXLO = true;
+    bool DirichletXLO = false;
 
-    bool DirichletXHI = true;
-    //bool DirichletXHI = false;
+    //bool DirichletXHI = true;
+    bool DirichletXHI = false;
 
-    bool DirichletYLO = true;
-    //bool DirichletYLO = false;
+    //bool DirichletYLO = true;
+    bool DirichletYLO = false;
 
-    bool DirichletYHI = true;
-    //bool DirichletYHI = false;
+    //bool DirichletYHI = true;
+    bool DirichletYHI = false;
 
 //    bool DirichletZLO = true;
     bool DirichletZLO = false;

src/EMSchur3DBase.cpp

         // First set up grid stuff
         if (Cvec) delete [] Cvec;
         //if (CSolver) delete [] CSolver;
-        Cvec    = new Eigen::SparseMatrix<Complex, Eigen::RowMajor> [Omegas.size()];
+        Cvec    = new Eigen::SparseMatrix<Complex, SPARSEMAJOR> [Omegas.size()];
 
         //#ifdef LEMMAUSEOMP
         //#pragma omp parallel for schedule(static, 1)