Working with Lemma 0.2

examples/EMSchur3D.cpp

 int main( int argc, char** argv ) {
 
     // BiCGSTAB Diagonal preconditioner
-    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::NewSP();
+    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, Eigen::IncompleteLUT<Complex> > >::NewSP();
 
     if (argc < 3) {
         std::cout << "EMSchur3D  <rgrid>  <1dmod>  <aemsurvey> " << std::endl;
         exit( EXIT_SUCCESS );
     }
 
+    #ifdef LEMMAUSEOMP
+    std::cout << "Eigen is using OpenMP!!!\n";
+    Eigen::initParallel();
+    #endif
+
     ///////////////////////////////////////////////////
     // Read in Grid
     auto GridRead = RectilinearGridReader::NewSP();

include/EMSchur3D.h

          */
         static std::shared_ptr< EMSchur3D > NewSP() {
             return std::make_shared< EMSchur3D<Solver> >( ctor_key() );
-            //return std::make_shared< EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > > >( ctor_key() ) ;
+            //return std::make_shared< EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > > >( ctor_key() ) ;
         }
 
         /** Default protected constructor, use New */
     template < class Solver >
     void EMSchur3D<Solver>::SolveSource ( std::shared_ptr<DipoleSource> Source, const int& isource ) {
 
+        std::cout << "In vanilla SolveSource" << std::endl;
+
         //  figure out which omega we are working with
         int iw = -1;
         for (int iiw=0; iiw<Omegas.size(); ++iiw) {
         FillPoints(dpoint);
         PrimaryField(Source, dpoint);
 
+        std::cout << "Done with primary field" << std::endl;
+
         // Allocate a ton of memory
         VectorXcr Phi    = VectorXcr::Zero(uns);
         VectorXcr ms(unx+uny+unz);  // mu sigma
         VectorXcr E0     = VectorXcr::Zero(unx+uny+unz);
 
         // Lets get cracking
+        std::cout << "Filling source terms" << std::endl;
         FillSourceTerms(ms, Se, E0, dpoint, Omegas[iw]);
+        std::cout << "Done source terms" << std::endl;
 
         /////////////////////////////////////////////////
         // LOG File
         logfile += to_string(isource) + std::string(".log");
         ofstream logio(logfile.c_str());
 
-        logio << *Source << std::endl;
+        std::cout << "just logging" << std::endl;
+//        logio << *Source << std::endl;
         logio << *Grid << std::endl;
         logio << *LayModel << std::endl;
+        std::cout << "dun logging" << std::endl;
 
         // solve for RHS
         int max_it(nx*ny*nz), iter_done(0);
 
     #ifdef HAVE_SUPERLUMT
     template<>
-    void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::BuildCDirectSolver() {
+    void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > >::BuildCDirectSolver() {
 
-        CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
+        CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [Omegas.size()];
 
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
     #endif
 
     template<>
-    void EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > >::BuildCDirectSolver() {
-        CSolver = new Eigen::SparseLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > [Omegas.size()];
+    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()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             jsw_timer timer;
             timer.begin();
     }
 
 //     template<>
-//     void EMSchur3D< Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
+//     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()];
 //         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::ColMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::NaturalOrdering<int> > [Omegas.size()];
+//     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()];
 //         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::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CSymSimplicialLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
+//     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()];
 //         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::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::CSymSimplicialLDLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::AMDOrdering<int> > [Omegas.size()];
+//     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()];
 //         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::ColMajor>, Eigen::IncompleteLUT<Complex> > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
+    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()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            CSolver[iw].preconditioner().setDroptol(1e-3);
+            //CSolver[iw].preconditioner().setFillfactor(1e3);
             jsw_timer timer;
             timer.begin();
             /*  Complex system */
     }
 
     template<>
-    void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
+    void EMSchur3D< Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::BiCGSTAB< Eigen::SparseMatrix<Complex, Eigen::RowMajor> > [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::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
-        for (int iw=0; iw<Omegas.size(); ++iw) {
-            //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
-            jsw_timer timer;
-            timer.begin();
-            /*  Complex system */
-            std::cout << "ConjugateGradient pattern analyzing C_" << iw << ",";
-            std::cout.flush();
-            CSolver[iw].analyzePattern( Cvec[iw] );
-            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
-            // factorize
-            timer.begin();
-            std::cout << "ConjugateGradient factorising C_" << iw << ", ";
-            std::cout.flush();
-            CSolver[iw].factorize( Cvec[iw] );
-            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
-        }
-    }
+//     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()];
+//         for (int iw=0; iw<Omegas.size(); ++iw) {
+//             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+//             jsw_timer timer;
+//             timer.begin();
+//             /*  Complex system */
+//             std::cout << "ConjugateGradient pattern analyzing C_" << iw << ",";
+//             std::cout.flush();
+//             CSolver[iw].analyzePattern( Cvec[iw] );
+//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+//             // factorize
+//             timer.begin();
+//             std::cout << "ConjugateGradient factorising C_" << iw << ", ";
+//             std::cout.flush();
+//             CSolver[iw].factorize( Cvec[iw] );
+//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+//         }
+//     }
 
 //     template<>
-//     void EMSchur3D<   Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
+//     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()];
 //         //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::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
+//     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()];
 //         //MPI_Init(NULL, NULL);
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
 //     }
 //
 //     template<>
-//     void EMSchur3D<   Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > > ::BuildCDirectSolver() {
-//         CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > [Omegas.size()];
+//     void EMSchur3D<   Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::RowMajor>, true > [Omegas.size()];
 //         //MPI_Init(NULL, NULL);
 //         for (int iw=0; iw<Omegas.size(); ++iw) {
 //             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();

include/EMSchur3DBase.h

 
     /** Matrix objects representing discrete operators
      */
-    Eigen::SparseMatrix<Complex, Eigen::ColMajor>*                  Cvec;
-    Eigen::SparseMatrix<Complex, Eigen::ColMajor>                   D;
+    Eigen::SparseMatrix<Complex, Eigen::RowMajor>*                  Cvec;
+    Eigen::SparseMatrix<Complex, Eigen::RowMajor>                   D;
 
     /** Squeezed matrices for reduced phi grid
      */
-    Eigen::SparseMatrix<Complex, Eigen::ColMajor>*                  Cvec_s;
-    Eigen::SparseMatrix<Complex, Eigen::ColMajor>                   D_s;
+    Eigen::SparseMatrix<Complex, Eigen::RowMajor>*                  Cvec_s;
+    Eigen::SparseMatrix<Complex, Eigen::RowMajor>                   D_s;
 
     /** number of cells in x, set when RectilinearGrid is attached */
     int nx;

src/EMSchur3DBase.cpp

     // ====================  FRIEND METHODS  =====================
 
     std::ostream &operator<<(std::ostream &stream, const EMSchur3DBase &ob) {
-        stream << ob.Serialize()  << "\n---\n"; // End of doc ---
+        stream << ob.Serialize()  << "\n";
         return stream;
     }
 
         // First set up grid stuff
         if (Cvec) delete [] Cvec;
         //if (CSolver) delete [] CSolver;
-        Cvec    = new Eigen::SparseMatrix<Complex> [Omegas.size()];
+        Cvec    = new Eigen::SparseMatrix<Complex, Eigen::RowMajor> [Omegas.size()];
 
         //#ifdef LEMMAUSEOMP
         //#pragma omp parallel for schedule(static, 1)
         Setup();
 
         std::cout << "Entering parallel solve" << std::endl;
-        #ifdef LEMMAUSEOMP
-        #pragma omp parallel for schedule(static,1)
-        #endif
+        //#ifdef LEMMAUSEOMP
+        //#pragma omp parallel for schedule(static,1)
+        //#endif
         for (int isource=0; isource<Survey->GetNumberOfSources(); ++isource) {
             SolveSource(Survey->GetSource(isource), isource);
         }