Adding files from SVN and porting over to CMake and latest Lemma

CMakeLists.txt

@@ -0,0 +1,35 @@
+add_subdirectory("src")
+add_library( emschur3d ${EMSCHUR3DSOURCE} )  
+target_include_directories( emschur3d PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include" )
+
+set_target_properties(emschur3d PROPERTIES 
+	VERSION  "${LEMMA_VERSION_NOQUOTES}"
+	SOVERSION "${LEMMA_VERSION_MAJOR}.${LEMMA_VERSION_MINOR}"
+	PROJECT_LABEL "FDEM1D ${LABEL_SUFFIX}"
+)
+
+# Linking
+target_link_libraries(emschur3d "lemmacore" "fdem1d" )
+
+# Linking
+if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT ) 
+	target_link_libraries(emschur3d ${VTK_LIBRARIES})
+endif()
+
+# Testing
+if (LEMMA_ENABLE_TESTING)
+	add_subdirectory(testing)
+endif()
+
+# Install
+install ( TARGETS emschur3d DESTINATION ${CMAKE_INSTALL_PREFIX}/lib )
+install ( FILES include/EMSchur3D  DESTINATION ${CMAKE_INSTALL_PREFIX}/include/Lemma ) 
+install ( DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_PREFIX}/include/Lemma  FILES_MATCHING PATTERN "*.h")
+
+#install ( DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_PREFIX}/include/Lemma/  FILES_MATCHING PATTERN "FDEM1D")
+#install ( DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_PREFIX}/include/Lemma/FDEM1D  FILES_MATCHING PATTERN "*.h")
+
+# Examples
+if (LEMMA_BUILD_EXAMPLES)
+	add_subdirectory(examples)
+endif()

include/EMSchur3D.h

@@ -0,0 +1,616 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      02/19/2015 01:10:39 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2015, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2015, Trevor Irons
+ * @copyright Copyright (c) 2011, Trevor Irons
+ * @copyright Copyright (c) 2011, Colorado School of Mines
+ */
+
+#ifndef  EMSCHUR3D_INC
+#define  EMSCHUR3D_INC
+
+#include "EMSchur3DBase.h"
+#include "CSymSimplicialCholesky.h"
+
+namespace Lemma {
+
+
+    /**
+      \brief   Templated concrete classes of EMSChur3DBase.
+      \details
+     */
+    template < class Solver >
+    class EMSchur3D : public EMSchur3DBase {
+
+        //friend std::ostream &operator<<(std::ostream &stream,
+        //        const EMSchur3D &ob);
+
+        public:
+
+        // ====================  LIFECYCLE     =======================
+
+        /**
+         * @copybrief LemmaObject::New()
+         * @copydetails LemmaObject::New()
+         */
+        static EMSchur3D* New() {
+            EMSchur3D<Solver>*  Obj = new EMSchur3D<Solver>("EMSchur3D");
+            Obj->AttachTo(Obj);
+            return Obj;
+        }
+
+        /**
+         *  @copybrief   LemmaObject::Delete()
+         *  @copydetails LemmaObject::Delete()
+         */
+        void Delete() {
+            this->DetachFrom(this);
+        }
+
+        // ====================  OPERATORS     =======================
+
+        // ====================  OPERATIONS    =======================
+
+        /** Solves a single source problem. This method is thread safe.
+         *  @param[in] Source is the source term for generating primary fields
+         *  @param[in] isource is the source index
+         */
+        void SolveSource( DipoleSource* Source , const int& isource);
+
+        /** Builds the solver for the C matrix */
+        void BuildCDirectSolver(  );
+
+        // ====================  ACCESS        =======================
+
+        // ====================  INQUIRY       =======================
+
+#ifdef HAVE_YAMLCPP
+//         /**
+//          *  Uses YAML to serialize this object.
+//          *  @return a YAML::Node
+//          */
+//         YAML::Node Serialize() const;
+//
+//         /**
+//          *   Constructs an object from a YAML::Node.
+//          */
+//         static EMSchur3D* DeSerialize(const YAML::Node& node);
+#endif
+
+        protected:
+
+        // ====================  LIFECYCLE     =======================
+
+        /** Default protected constructor, use New */
+        EMSchur3D (const std::string& name) : EMSchur3DBase(name), CSolver(NULL) {
+        }
+
+// #ifdef HAVE_YAMLCPP
+//         /** Protected DeDerializing constructor, use factory DeSerialize  method*/
+//         EMSchur3D (const YAML::Node& node): EMSchur3DBase(node), CSolver(NULL) {
+//         }
+// #endif
+
+        /** Default protected destructor, use Delete */
+        ~EMSchur3D () {
+            // TODO delete arrays
+        }
+
+        /**
+         *  @copybrief   LemmaObject::Release()
+         *  @copydetails LemmaObject::Release()
+         */
+        void Release() {
+            delete this;
+        }
+
+        private:
+
+        // ====================  DATA MEMBERS  =========================
+
+        /** The templated solver for C */
+        Solver*     CSolver;
+
+        Eigen::SparseMatrix<Complex>  Csym;
+
+    }; // -----  end of class  EMSchur3D  -----
+
+
+    ////////////////////////////////////////////////////////////////////////////////////////
+    //                    Implimentation and Specialisations                              //
+    ////////////////////////////////////////////////////////////////////////////////////////
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  EMSchur3D
+    //      Method:  SolveSource
+    //--------------------------------------------------------------------------------------
+    template < class Solver >
+    void EMSchur3D<Solver>::SolveSource ( DipoleSource* Source, const int& isource ) {
+
+        //  figure out which omega we are working with
+        int iw = -1;
+        for (int iiw=0; iiw<Omegas.size(); ++iiw) {
+           if (Omegas[iiw] - Source->GetAngularFrequency(0) < 1e-3 ) {
+               iw = iiw;
+           }
+        }
+        if (iw == -1) {
+            std::cerr << "FREQUENCY DOOM IN EMSchur3D::SolveSource \n";
+            exit(EXIT_FAILURE);
+        }
+
+        ///////////////////////////////////
+        // Set up primary fields
+        // TODO, this is a little stupid as they all share the same points. We need to extend
+        //       EmEARTH to be able to input a grid so that points are not explicitly needed like
+        //       this. This requires some care as calcs are made on faces.
+        //       Alternatively, the bins function of ReceiverPoints could be extended quite easily.
+        //       This may be the way to do this.
+
+        Lemma::ReceiverPoints* dpoint = Lemma::ReceiverPoints::New();
+
+        FillPoints(dpoint);
+        PrimaryField(Source, dpoint);
+
+        // Allocate a ton of memory
+        VectorXcr Phi    = VectorXcr::Zero(uns);
+        VectorXcr ms(unx+uny+unz);  // mu sigma
+
+        // Vector potential (A) Vector and phi
+        VectorXcr Se     = VectorXcr::Zero(unx+uny+unz);
+        //VectorXcr A      = VectorXcr::Zero(unx+uny+unz);
+        VectorXcr E      = VectorXcr::Zero(unx+uny+unz);
+        VectorXcr E0     = VectorXcr::Zero(unx+uny+unz);
+
+        // Lets get cracking
+        FillSourceTerms(ms, Se, E0, dpoint, Omegas[iw]);
+
+        /////////////////////////////////////////////////
+        // LOG File
+        std::string logfile (ResFile);
+        logfile += to_string(isource) + std::string(".log");
+        ofstream logio(logfile.c_str());
+
+        logio << *Source << std::endl;
+        logio << *Grid << std::endl;
+        logio << *LayModel << std::endl;
+
+        // solve for RHS
+        int max_it(nx*ny*nz), iter_done(0);
+        Real tol(3e-16), errorn(0);
+        logio << "solving RHS for source " << isource << std::endl;
+
+        // TODO, this is stupid, try and get rid of this copy!
+        Eigen::SparseMatrix<Complex>  Cc  = Cvec[iw];
+
+        jsw_timer timer;
+        jsw_timer timer2;
+
+        timer.begin();
+        timer2.begin();
+
+        /////////////////////////////////////////
+        // Solve for RHS
+        VectorXcr A = CSolver[iw].solve(Se);
+
+//         // Solve Real system instead
+           // The Real system is quasi-definite, though an LDLT decomposition exists, CHOLMOD doesn't find it.
+           // An LU can be done on this, but compute performance is very similiar to the complex system, and diagonal pivoting
+           // cannot be assumed to be best, hurting solve time.
+//         /* EXPERIMENTAL */
+//         VectorXr b2 = VectorXr::Zero(2*(unx+uny+unz));
+//         b2.head(unx+uny+unz) = Se.real();
+//         b2.tail(unx+uny+unz) = Se.imag();
+//         VectorXr A2 = CReSolver[iw].solve(b2);
+//         A.real() =   A2.head( unx+uny+unz );
+//         A.imag() =  -A2.tail( unx+uny+unz ); // Due to decomp. negative!
+//         /* END EXPERIMENTAL */
+
+        VectorXcr ADiv = D*A;  // ADiv == RHS == D C^I Se
+        VectorXcr Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() - Se.array());
+        logio << "|| Div(A) || = " << ADiv.norm()
+                // << " in " << iter_done << " iterations"
+              //<<  " with error " << errorn << "\t"
+              << "\tInital solution error "<<   Error.norm()  // Iteritive info
+              << "\ttime " << timer.end() << std::endl;
+
+        //VectorXcr ADivMAC = ADiv.array() * MAC.array().cast<Complex>();
+        //logio << "|| Div(A) || on MAC grid " << ADivMAC.norm() << std::endl;
+
+        /////////////////////
+        // Solve for Phi
+        logio << "Solving for Phi " << std::flush;
+        timer.begin();
+        tol = 1e-18;
+        int success(2);
+
+        success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
+        //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
+
+        /* Restart if necessary */
+        int nrestart(1);
+        // TODO send MAC to implicitbicgstab?
+        while (success == 2 && nrestart < 18 && iter_done > 1) {
+            success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
+            //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
+            nrestart += 1;
+        }
+
+        logio << "Implicit BiCGStab solution in " << iter_done << " iterations."
+                << " with error " << std::setprecision(8) << std::scientific << errorn << std::endl;
+        logio << "time "<< timer.end() << " [s]" << std::endl;
+
+
+        E = ms.array()*(D.transpose()*Phi).array(); // Temp, field due to charge
+
+        /////////////////////////////////////
+        // Compute A
+        /////////////////////////////////////
+        logio << "Solving for A using phi" << std::endl;
+        std::cout << "Solving for A" << std::endl;
+        max_it = nx*ny*nz;
+        tol = 5e-16;
+        errorn = 0;
+        iter_done = 0;
+
+        timer.begin();
+
+        A = CSolver[iw].solve( (Se-E).eval() ); // UmfPack requires eval?
+
+        VectorXcr ADiv2 = D*A;
+        logio << "|| Div(A) || = " << ADiv2.norm() ;
+              //" in " << iter_done << " iterations"
+              //<<  " with error " << errorn << "\t";
+
+        // Report error of solutions
+        Error = ((Cc.selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());
+        logio << "\tsolution error " << Error.norm()
+              << std::fixed << std::setprecision(2) << "\ttime " << timer.end() << "\ttotal time " << timer2.end() << std::endl;
+        logio.close();
+
+        //////////////////////////////////////
+        // Update Fields and report
+        E.array() = Complex(0,-Omegas[iw])*A.array() - (D.transpose()*Phi).array();   // Secondary Field Only
+        VectorXcr B = StaggeredGridCurl(A);
+
+        WriteVTKResults( ResFile+ to_string(isource), A, Se, E0, E , Phi, ADiv, ADiv2, B);
+
+        dpoint->Delete();
+        return ;
+
+    }		// -----  end of method EMSchur3D::SolveSource  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  EMSchur3DBase
+    //      Method:  BuildCDirectSolver
+    //--------------------------------------------------------------------------------------
+    template < class Solver >
+    void EMSchur3D<Solver>::BuildCDirectSolver (  ) {
+
+        CSolver = new Solver[Omegas.size()];
+
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+
+            jsw_timer timer;
+            timer.begin();
+
+            /*  Complex system */
+            /*
+            std::cout << "Generic solver 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 << "Generic solver factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            */
+
+            std::cerr << "No solver Specified!" << iw << ",";
+            exit(EXIT_FAILURE);
+            //CSolver[iw].compute( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+
+        }
+    }
+
+    #ifdef HAVE_SUPERLUMT
+    template<>
+    void EMSchur3D< Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > >::BuildCDirectSolver() {
+
+        CSolver = new Eigen::SuperLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
+
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            jsw_timer timer;
+            timer.begin();
+
+            /* SuperLU */
+            //CSolver[iw].options().DiagPivotThresh = 0.01;
+            //CSolver[iw].options().SymmetricMode = YES;
+            //CSolver[iw].options().ColPerm = MMD_AT_PLUS_A;
+            //CSolver[iw].options().Trans = NOTRANS;
+            //CSolver[iw].options().ConditionNumber = NO;
+            //std::cout << "SuperLU options:\n";
+            //std::cout << "\tPivot Threshold: " << CSolver[iw].options().DiagPivotThresh << std::endl;
+            //std::cout << "\tSymmetric mode: " << CSolver[iw].options().SymmetricMode << std::endl;
+            //std::cout << "\tEquilibrate: " << CSolver[iw].options().Equil << std::endl;
+            //std::cout << "\tCol Permutation: " << CSolver[iw].options().ColPerm << std::endl;
+            //std::cout << "\tTrans: " << CSolver[iw].options().Trans << std::endl;
+            //std::cout << "\tCondition Number: " << CSolver[iw].options().ConditionNumber << std::endl;
+
+            /*  Complex system */
+            std::cout << "SuperLU_MT 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 << "SuperLU_MT factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+
+        }
+    }
+    #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()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            jsw_timer timer;
+            timer.begin();
+
+            CSolver[iw].isSymmetric(true);
+            CSolver[iw].setPivotThreshold(0.0);
+
+            /*  Complex system */
+            std::cout << "SparseLU 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 << "SparseLU 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::CholmodSupernodalLLT< Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::CholmodSupernodalLLT< 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 << "CholmodSupernodalLLT pattern analyzing C_" << iw << ",";
+//             std::cout.flush();
+//             CSolver[iw].analyzePattern( Csym );
+//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+//             /* factorize */
+//             timer.begin();
+//             std::cout << "CholmodSupernodalLLT factorising C_" << iw << ", ";
+//             std::cout.flush();
+//             CSolver[iw].factorize( Csym );
+//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+//         }
+//     }
+
+    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()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            jsw_timer timer;
+            timer.begin();
+            /*  Complex system */
+            std::cout << "CSymSimplicialLLT<NaturalOrdering> pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            CSolver[iw].analyzePattern( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+            /* factorize */
+            timer.begin();
+            std::cout << "CSymSimplicialLLT<NaturalOrdering> factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
+
+    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()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            jsw_timer timer;
+            timer.begin();
+            /*  Complex system */
+            std::cout << "CSymSimplicialLLT<AMDOrdering> 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 << "CSymSimplicialLLT<AMDOrdering> 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::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()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            jsw_timer timer;
+            timer.begin();
+            /*  Complex system */
+            std::cout << "CSymSimplicialLDLT<AMDOrdering> pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            CSolver[iw].analyzePattern( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+            /* factorize */
+            timer.begin();
+            std::cout << "CSymSimplicialLDLT<AMDOrdering> factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
+
+    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()];
+        for (int iw=0; iw<Omegas.size(); ++iw) {
+            Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            jsw_timer timer;
+            timer.begin();
+            /*  Complex system */
+            std::cout << "BiCGSTAB(ILU) pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            CSolver[iw].analyzePattern( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+            /* factorize */
+            timer.begin();
+            std::cout << "BiCGSTAB(ILU) factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
+
+    template<>
+    void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > > ::BuildCDirectSolver() {
+        CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [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 << "BiCGSTAB pattern analyzing C_" << iw << ",";
+            std::cout.flush();
+            CSolver[iw].analyzePattern( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+            // factorize
+            timer.begin();
+            std::cout << "BiCGSTAB factorising C_" << iw << ", ";
+            std::cout.flush();
+            CSolver[iw].factorize( Csym );
+            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+        }
+    }
+
+    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::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::PastixLLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
+//         //MPI_Init(NULL, NULL);
+//         for (int iw=0; iw<Omegas.size(); ++iw) {
+//             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+//             jsw_timer timer;
+//             timer.begin();
+//             /*  Complex system */
+//             std::cout << "PaStiX LLT 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 << "PaStiX LLT 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::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::PastixLDLT<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower > [Omegas.size()];
+//         //MPI_Init(NULL, NULL);
+//         for (int iw=0; iw<Omegas.size(); ++iw) {
+//             //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+//             jsw_timer timer;
+//             timer.begin();
+//             /*  Complex system */
+//             std::cout << "PaStiX LDLT 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 << "PaStiX LDLT factorising C_" << iw << ", ";
+//             std::cout.flush();
+//             CSolver[iw].factorize( Cvec[iw] );
+//             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+//             std::cout << "INFO " << CSolver[iw].info(  ) << std::endl;
+//         }
+//     }
+//
+//     template<>
+//     void EMSchur3D<   Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > > ::BuildCDirectSolver() {
+//         CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > [Omegas.size()];
+//         //MPI_Init(NULL, NULL);
+//         for (int iw=0; iw<Omegas.size(); ++iw) {
+//             Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+//             jsw_timer timer;
+//             timer.begin();
+//             /*  Complex system */
+//             std::cout << "PaStiX LU pattern analyzing C_" << iw << ",";
+//             std::cout.flush();
+//             CSolver[iw].compute( Csym );
+//             std::cout << "PaStiX LU Done C_" << iw << std::endl;;
+// //             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+// //             // factorize
+// //             timer.begin();
+// //             std::cout << "PaStiX LU factorising C_" << iw << ", ";
+// //             std::cout.flush();
+// //             CSolver[iw].factorize( Csym );
+// //             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
+//         }
+//     }
+
+}		// -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef EMSCHUR3D_INC  -----
+

include/EMSchur3DBase.h

@@ -0,0 +1,438 @@
+// ===========================================================================
+//
+//       Filename:  EMSchur3DBase.h
+//
+//        Created:  09/20/2013 04:35:57 PM
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  University of Utah,
+//                  Colorado School of Mines
+//                  US Geological Survey
+//
+//          Email:  tirons@egi.utah.edu
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     09/20/2013
+  @version  $Id$
+ **/
+
+
+#ifndef  EMSCHUR3DBASE_INC
+#define  EMSCHUR3DBASE_INC
+
+#include <LemmaCore>
+#include <FDEM1D>
+
+//#include "LemmaObject.h"
+//#include "rectilineargrid.h"
+//#include "RectilinearGridVTKExporter.h"
+//#include "ASCIIParser.h"
+//#include "AEMSurvey.h"
+//#include "receiverpoints.h"
+//#include "layeredearthem.h"
+//#include "emearth1d.h"
+
+#include "timer.h"
+#include <Eigen/Sparse>
+#include "bicgstab.h"
+
+// Solvers
+#ifdef HAVE_PASTIX
+#include <Eigen/PaStiXSupport>
+#endif
+
+#ifdef HAVE_METIS
+#include <Eigen/MetisSupport>
+#endif
+
+#ifdef HAVE_SUPERLU
+#include <Eigen/SuperLUSupport>
+#endif
+
+#ifdef HAVE_SUPERLUMT
+#include <Eigen/SuperLUMTSupport>
+#endif
+
+#ifdef HAVE_SPQR
+#include <Eigen/SPQRSupport>
+#endif
+
+// Cholmod Support won't compile typedef issue
+// #ifdef HAVE_CHOLMOD
+// #include <Eigen/CholmodSupport>
+// #endif
+//
+// // Cholmod Support won't compile
+// #ifdef HAVE_UMFPACK
+// #include <Eigen/UmfPackSupport>
+// #endif
+
+namespace Lemma {
+
+/**
+ \defgroup EMSchur3DBase EMSchur3DBase
+  Provides 3D solution to Maxwell's equations.
+ */
+
+enum SOLVER{ SPARSELU, SimplicialLLT, SimplicialLDLT, BiCGStab, SparseQR };
+
+
+/**
+  @class EMSchur3DBase
+  \ingroup EMSchur3DBase
+  \brief Provides a 3D solution to Maxwell's equations.
+  \details 3D finite difference solution to maxwells equations
+            using a SCHUR decomposition on a staggered grid.
+  Performs a Schur decomposition on the vector scalar formulation of
+  Maxwell's equations.
+   \f[
+  -\nabla^2 (\mathbf{A}) - \jmath \omega \mu \sigma  \mathbf{A} - \mu \sigma \nabla (\phi) = -  \mu  \mathbf{J}_s
+   \f]
+
+Which can be written in the functional form
+\f[ \begin{pmatrix}
+      -\nabla^2 + \jmath \omega \mu \sigma   & \mu \sigma \nabla  \\
+       \nabla \cdot  & 0
+    \end{pmatrix}
+    \begin{pmatrix}   \mathbf{A} \\  \phi   \end{pmatrix}
+     = \begin{pmatrix}      \mathbf{s}_E \\   0      \end{pmatrix}
+\f]
+Using the notation
+\f[ \begin{pmatrix}
+        \mathbf{C} & \mathbf{B} \\
+        \mathbf{D} & \mathbf{0}
+    \end{pmatrix}  \begin{pmatrix} \mathbf{A} \\ \phi \end{pmatrix} =
+    \begin{pmatrix} \mathbf{s}_E \\ 0 \end{pmatrix}
+\f]
+Which is decomposed for seperate solutions to \f$ \mathbf{A}, \phi \f$ using a Schur decomposition
+  \f[   \begin{matrix}
+        \mathbf{D}\mathbf{C}^{-1}\mathbf{B} \phi & = \mathbf{D} \mathbf{C}^{-1} \mathbf{s}_E \\
+        \mathbf{C}\mathbf{A}                     & = \mathbf{s}_E - \mathbf{G} \phi
+        \end{matrix}
+  \f]
+
+Where \f$ \mathbf{B} = \mu \sigma \mathbf{D}^T \f$. Additional algorithmic details may be found at
+@verbatim
+@inproceedings{doi:10.1190/segam2012-0896.1,
+  author = {Trevor Irons and Yaoguo Li and Jason R. McKenna},
+  title = {3D frequency-domain electromagnetics modeling using decoupled scalar and vector potentials},
+  booktitle = {SEG Technical Program Expanded Abstracts 2012},
+  chapter = {112},
+  year = {2012},
+  pages = {1-6},
+  doi = {10.1190/segam2012-0896.1},
+  URL = {http://library.seg.org/doi/abs/10.1190/segam2012-0896.1},
+  eprint = {http://library.seg.org/doi/pdf/10.1190/segam2012-0896.1}
+}
+@endverbatim
+*/
+
+//template< class Solver >
+class EMSchur3DBase : public LemmaObject {
+
+    friend std::ostream &operator<<(std::ostream &stream,
+            const EMSchur3DBase &ob);
+
+    protected:
+    struct ctor_key {};
+
+    public:
+
+    // ====================  LIFECYCLE     =======================
+
+    /** Default protected constructor, use New */
+    explicit EMSchur3DBase ( const ctor_key& );
+
+    /** Default protected constructor, use New */
+    explicit EMSchur3DBase ( const YAML::Node& node, const ctor_key& );
+
+    /** Default protected destructor, use Delete */
+    virtual ~EMSchur3DBase ();
+
+    /**
+     * Initialises antenna to contain no points, with no current
+     * and no frequency. NumberOfTurns set to 1
+     */
+    static std::shared_ptr<EMSchur3DBase> NewSP();
+
+    /**
+     * Provides deep copy
+     */
+    virtual std::shared_ptr<EMSchur3DBase> Clone() const ;
+
+    /**
+     *  Uses YAML to serialize this object.
+     *  @return a YAML::Node
+     */
+    YAML::Node Serialize() const;
+
+    /**
+     *   Constructs an object from a YAML::Node.
+     */
+     static std::shared_ptr<EMSchur3DBase> DeSerialize( const YAML::Node& node );
+
+    // ====================  OPERATORS     =======================
+
+    // ====================  OPERATIONS    =======================
+
+    /* Performs the solution
+     * This is thread safe. TODO, but where should the results go? Just to file?
+     * Who does the parsing? Actually I think this method is the wrong place to talk
+     * about that. This is just a big red button. The details should be worked out in private
+     * methods, that this could well call. Still though, where should the damn results go. What if
+     * someone wants to use them *right now*, and not go through file IO. This is a good cause for
+     * fixing the model class. So the result will be a final RectGrid (or Model) where the results live.
+     * THEN users can call a seperate WriteToVTK or whatever based on that.
+     */
+    void Solve( );
+
+    // ====================  ACCESS        =======================
+
+    /** Sets the RectilinearGrid to use
+     *  @param[in] Grid is a pointer to the Grid to be used.
+     */
+    void SetRectilinearGrid( std::shared_ptr<RectilinearGrid> Grid);
+
+    /** Sets the RectilinearGrid to use
+     *  @param[in] Grid is a pointer to the Grid to be used.
+     */
+    //void SetAEMSurvey(AEMSurvey* Survey);
+
+    /** Sets the prefix for results files (.log and .vtr) the source fiducial is added as well
+     */
+    void SetResFileName(const std::string& filename);
+
+    /** Sets the solver to use to invert the C matrix. This is done many times. If you are reusing the same matrix
+        for numerous simulations, it may be benefitial to use the direct (SPARSELU) solver. For one off calculations the BiCGSTAB
+        is a good choice. Default is SPARSELU.
+     */
+    void SetCSolver(const SOLVER& CSOLVER);
+
+    /**
+     *  Sets the LayredEarthEM model that will be used for the primary field calculation as well as deterimining the
+     *  bacground conductivity file.
+     @  @param[in] LayModel is a pointer to the LayeredEarthEM model to use.
+     */
+    void SetLayeredEarthEM( std::shared_ptr<LayeredEarthEM> LayModel );
+
+    /**
+     * Loads a MeshTools conductivity model.
+     * @param[in] fname is the file to load.
+     */
+    void LoadMeshToolsConductivityModel( const std::string& fname );
+
+    /**
+     *  Writes out results to into a vtkRectilinearGrid file
+     *  @param[in] fname is the file name that is created, the .vtr suffix is added.
+     */
+    void WriteVTKResults( const std::string& fname, Eigen::Ref<VectorXcr> A,
+          Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0, Eigen::Ref<VectorXcr> E,
+          Eigen::Ref<VectorXcr> Phi, Eigen::Ref<VectorXcr> ADiv, Eigen::Ref<VectorXcr> ADiv2,
+          Eigen::Ref<VectorXcr> B
+        );
+
+    // ====================  INQUIRY       =======================
+
+    /** Returns the name of the underlying class, similiar to Python's type */
+    virtual std::string GetName() const;
+
+    protected:
+
+    // ====================  LIFECYCLE     =======================
+
+    //private:
+
+    /** Builds the C matrix */
+    void BuildC( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw);
+
+    /* Builds the C matrix as a block real system. Benefits of this are the availability of an
+     *  LDL^T decomposition. Also, as complex number in C++ are templates and will necessarily have
+     *  real and imaginary parts, this formulation will have a reduced cost, due to less computations
+     *  with the zero valued imaginary parts (off-diagonals)
+     *  The \f$C \in I^3\f$ matrix is instead written as
+     *   [  C_r   C_i ] [ x_r ] = [ b_r ]
+     *   [  C_i  -C_r ] [ x_i ]   [ b_i ]
+     */
+    //void BuildCReal( Real*** sigmax, Real*** sigmay, Real*** sigmaz, const int& iw);
+
+    /** Builds the C matrix */
+    void BuildCPreconditioner( const int& iw);
+
+    /** Builds the C matrix */
+    virtual void BuildCDirectSolver( )=0;
+
+    /** Fills the actual points on the grid that 1D primary field calculations need to be made.
+        @todo  This is a little stupid as all threads share the same points. Stupid in that right now
+               this is done for every calculation. Not a huge amount of time is spent here, I suppose some extra memory
+               though. We need to extend
+               EmEARTH to be able to input a grid so that points are not explicitly needed like
+               this. This requires some care as calcs are made on faces.
+               Alternatively, the bins function of FieldPoints could be extended quite easily.
+               This may be the way to do this.
+     */
+    void FillPoints( std::shared_ptr<FieldPoints> Points );
+
+    /** Builds D/G
+     */
+    void BuildD( );
+
+    /** Used to manage tradititional C 3D array */
+    template <typename T>
+    void Allocate3DScalar(T ***&Array, T init) {
+        Array  = new T**[nx];
+        for (int ix=0; ix<nx; ix++){
+            Array[ix] = new T*[ny];
+            for (int iy=0; iy<ny; iy++){
+                Array[ix][iy]  = new T[nz];
+                for (int iz=0; iz<nz; iz++) Array[ix][iy][iz] = init;
+            }
+        }
+        return;
+    }
+
+    /** Used to manage tradititional C 3D array */
+    template <typename T>
+    void Delete3DScalar(T ***&Array) {
+        for (int ix=0; ix<nx; ix++){
+            for (int iy=0; iy<ny; iy++){
+                delete [] Array[ix][iy];
+            }
+            delete [] Array[ix];
+        }
+        delete [] Array;
+        Array = NULL;
+        return;
+    }
+
+    /**
+     * This is called just before solve and gets all shared objects ready to go
+     */
+    void Setup( );
+
+    /** Solves a single source problem. This method is thread safe.
+     *  @param[in] Source is the source term for generating primary fields
+     *  @param[in] isource is the source index
+     */
+    virtual void SolveSource( DipoleSource* Source , const int& isource)=0;
+
+    /** Computes the primary field
+     */
+    void PrimaryField( std::shared_ptr<DipoleSource> Source, std::shared_ptr<FieldPoints> dpoint);
+
+    /**
+     *  Fills the vectors that are called in
+     */
+    void FillSourceTerms(  Eigen::Ref<VectorXcr>  ms,
+                                        Eigen::Ref<VectorXcr> Se, Eigen::Ref<VectorXcr> E0,
+                                        std::shared_ptr<FieldPoints> dpoint, const Real& omega );
+
+    /** Computes the curl of A on the staggered grid
+     */
+    VectorXcr StaggeredGridCurl(Eigen::Ref<VectorXcr> A);
+
+    // ====================  DATA MEMBERS  =========================
+
+    /** Grid over which operators are active */
+    std::shared_ptr<RectilinearGrid>                    Grid;
+
+    /* Used to help dump results */
+    //std::shared_ptr<RectilinearGridVTKExporter>         VTKGridExporter;
+
+    /** Class containing information about the AEM survey */
+    //AEMSurvey*                      Survey;
+
+    std::shared_ptr<LayeredEarthEM>                     LayModel;
+
+    /** Matrix objects representing discrete operators
+     */
+    Eigen::SparseMatrix<Complex, Eigen::ColMajor>*                  Cvec;
+    Eigen::SparseMatrix<Complex, Eigen::ColMajor>                   D;
+
+    /** Squeezed matrices for reduced phi grid
+     */
+    Eigen::SparseMatrix<Complex, Eigen::ColMajor>*                  Cvec_s;
+    Eigen::SparseMatrix<Complex, Eigen::ColMajor>                   D_s;
+
+    /** number of cells in x, set when RectilinearGrid is attached */
+    int nx;
+
+    /** number of cells in y, set when RectilinearGrid is attached */
+    int ny;
+
+    /** number of cells in z, set when RectilinearGrid is attached */
+    int nz;
+
+    /** number of fields/faces in x, unwrapped x */
+    int unx;
+
+    /** number of fields/faces in y, unwrapped y */
+    int uny;
+
+    /** number of fields/faces in z, unwrapped z */
+    int unz;
+
+    /** number of cell centres, unwrapped scalars */
+    int uns;
+
+    /** name for log files and VTK files */
+    std::string ResFile;
+
+    /** frequency of source */
+    VectorXr Omegas;
+
+    /** Conductivity model */
+    //Complex ***sigma_jwe;
+    Real ***sigma;
+
+    /** Conductivity model minus reference model */
+    //Complex ***sigmap;
+    Real ***sigmap;
+
+    /** rectilinear grid spacing in x */
+    VectorXr hx;
+
+    /** rectilinear grid spacing in y */
+    VectorXr hy;
+
+    /** rectilinear grid spacing in z */
+    VectorXr hz;
+
+    /** inverse of hx */
+    VectorXr ihx;
+
+    /** inverse of hx squared */
+    VectorXr ihx2;
+
+    /** inverse of hy */
+    VectorXr ihy;
+
+    /** inverse of hy squared */
+    VectorXr ihy2;
+
+    /** inverse of hz */
+    VectorXr ihz;
+
+    /** inverse of hz squared */
+    VectorXr ihz2;
+
+    /** Marker for air cells */
+    VectorXi MAC;
+
+    /** Marker for air cells */
+    std::vector<int> idx;
+
+    private:
+
+        static constexpr auto CName = "EMSchur3DBase";
+
+}; // -----  end of class  EMSchur3DBase  -----
+
+}
+
+#endif   // ----- #ifndef EMSCHUR3BASE_INC  -----

include/bicgstab.h

@@ -0,0 +1,936 @@
+// ===========================================================================
+//
+//       Filename:  bicgstab.h
+//
+//    Description:
+//
+//        Version:  0.0
+//        Created:  10/27/2009 03:15:06 PM
+//       Revision:  none
+//       Compiler:  g++ (c++)
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+#ifdef CHOLMODPRECONDITION
+#include <Eigen/CholmodSupport>
+#endif // CHOLMODPRECONDITION
+
+//#include "unsupported/Eigen/IterativeSolvers"
+//#include <unsupported/Eigen/SuperLUSupport>
+
+#include <iostream>
+#include <string>
+#include <complex>
+#include <fstream>
+#include "lemma.h"
+#include "timer.h"
+
+using namespace Eigen;
+using namespace Lemma;
+
+//typedef Eigen::VectorXcd VectorXcr;
+typedef Eigen::SparseMatrix<Complex> SparseMat;
+
+
+// On Input
+// A = Matrix
+// B = Right hand side
+// X = initial guess, and solution
+// maxit = maximum Number of iterations
+// tol = error tolerance
+// On Output
+// X real solution vector
+// errorn = Real error norm
+int bicgstab(const SparseMat &A, const SparseMat &M, const VectorXcr &b, VectorXcr &x,
+                int &max_it, Real &tol, Real &errorn, int &iter_done,
+                const bool& banner = true) {
+
+    Complex omega, rho, rho_1, alpha, beta;
+    Real bnrm2, eps, errmin;
+    int n, iter; //, istat;
+
+    // Determine size of system and init vectors
+    n = x.size();
+    VectorXcr r(n);
+    VectorXcr r_tld(n);
+    VectorXcr p(n);
+    VectorXcr v(n);
+    VectorXcr p_hat(n);
+    VectorXcr s(n);
+    VectorXcr s_hat(n);
+    VectorXcr t(n);
+    VectorXcr xmin(n);
+
+    if (banner) {
+        std::cout << "Start BiCGStab, memory needed: "
+                  <<  (sizeof(Complex)*(9+2)*n/(1024.*1024*1024)) << " [Gb]\n";
+    }
+
+    // Initialise
+    iter_done = 0;
+    v.setConstant(0.); // not necessary I don't think
+    t.setConstant(0.);
+    eps = 1e-100;
+
+    bnrm2 = b.norm();
+    if (bnrm2 == 0) {
+        x.setConstant(0.0);
+        errorn = 0;
+        std::cerr << "Trivial case of Ax = b, where b is 0\n";
+        return (0);
+    }
+
+    // If there is an initial guess
+    if ( x.norm() ) {
+        r = b - A.selfadjointView<Eigen::Upper>()*x;
+        //r = b - A*x;
+    } else {
+        r = b;
+    }
+
+    errorn = r.norm() / bnrm2;
+    omega = 1.;
+    r_tld = r;
+    errmin = 1e30;
+
+    // Get down to business
+    for (iter=0; iter<max_it; ++iter) {
+
+        rho = r_tld.dot(r);
+        if ( abs(rho) < eps) return (0);
+
+        if (iter > 0) {
+            beta = (rho/rho_1) * (alpha/omega);
+            p = r.array() + beta*(p.array()-omega*v.array()).array();
+        } else {
+            p = r;
+        }
+
+        // Use pseudo inverse to get approximate answer
+        //#pragma omp sections
+        p_hat = M*p;
+        //v = A*p_hat; // TODO double check
+        v = A.selfadjointView<Eigen::Upper>()*p_hat; // TODO double check
+
+        alpha = rho / r_tld.dot(v);
+        s = r.array() - alpha*v.array();
+        errorn = s.norm()/bnrm2;
+
+        if (errorn < tol && iter > 1) {
+            x.array() += alpha*p_hat.array();
+            return (0);
+        }
+
+        s_hat = M*s;
+        t = A.selfadjointView<Eigen::Upper>()*s_hat;
+        //t = A*s_hat;
+
+        omega = t.dot(s)  / t.dot(t);
+        x.array() += alpha*p_hat.array() + omega*s_hat.array();
+        r = s.array() - omega*t.array();
+        errorn = r.norm() / bnrm2;
+        iter_done = iter;
+
+        if (errorn < errmin) {
+            // remember the model with the smallest norm
+            errmin = errorn;
+            xmin = x;
+        }
+
+        if ( errorn <= tol ) return (0);
+        if ( abs(omega) < eps ) return (0);
+        rho_1 = rho;
+
+    }
+    return (0);
+}
+
+template <typename Preconditioner>
+bool preconditionedBiCGStab(const SparseMat &A, const Preconditioner &M,
+        const Ref< VectorXcr const > b,
+        Ref <VectorXcr > x,
+        const int &max_it, const Real &tol,
+        Real &errorn, int &iter_done) {
+
+    Complex omega, rho, rho_1, alpha, beta;
+    Real bnrm2, eps;
+    int n, iter;
+    Real tol2 = tol*tol;
+
+    // Determine size of system and init vectors
+    n = x.size();
+
+    VectorXcd r(n);
+    VectorXcd r_tld(n);
+    VectorXcd p(n);
+    VectorXcd s(n);
+    VectorXcd s_hat(n);
+    VectorXcd p_hat(n);
+    VectorXcd v = VectorXcr::Zero(n);
+    VectorXcd t = VectorXcr::Zero(n);
+
+    //std::cout << "Start BiCGStab, memory needed: "
+    //          <<  (sizeof(Complex)*(8+2)*n/(1024.*1024)) / (1024.) << " [Gb]\n";
+
+    // Initialise
+    iter_done = 0;
+    eps = 1e-100;
+
+    bnrm2 = b.squaredNorm();
+    if (bnrm2 == 0) {
+        x.setConstant(0.0);
+        errorn = 0;
+        std::cerr << "Trivial case of Ax = b, where b is 0\n";
+        return (false);
+    }
+
+    // If there is an initial guess
+    if ( x.squaredNorm() ) {
+        r = b - A.selfadjointView<Eigen::Upper>()*x;
+    } else {
+        r = b;
+    }
+
+    errorn = r.squaredNorm() / bnrm2;
+    omega = 1.;
+    r_tld = r;
+
+    // Get down to business
+    for (iter=0; iter<max_it; ++iter) {
+
+        rho = r_tld.dot(r);
+        if (abs(rho) < eps) {
+            std::cerr << "arbitrary orthogonality issue in bicgstab\n";
+            std::cerr << "consider eigen restarting\n";
+            return (false);
+        }
+
+        if (iter > 0) {
+            beta = (rho/rho_1) * (alpha/omega);
+            p = r + beta*(p-omega*v);
+        } else {
+            p = r;
+        }
+
+        p_hat = M.solve(p);
+        v.noalias() = A.selfadjointView<Eigen::Upper>()*p_hat;
+
+        alpha = rho / r_tld.dot(v);
+        s = r - alpha*v;
+        errorn = s.squaredNorm()/bnrm2;
+
+        if (errorn < tol2 && iter > 1) {
+            x = x + alpha*p_hat;
+            errorn = std::sqrt(errorn);
+            return (true);
+        }
+
+        s_hat = M.solve(s);
+        t.noalias() = A.selfadjointView<Eigen::Upper>()*s_hat;
+
+        omega = t.dot(s)  / t.dot(t);
+        x += alpha*p_hat + omega*s_hat;
+        r = s - omega*t;
+        errorn = r.squaredNorm() / bnrm2;
+        iter_done = iter;
+
+        if ( errorn <= tol2 || abs(omega) < eps) {
+            errorn = std::sqrt(errorn);
+            return (true);
+        }
+
+        rho_1 = rho;
+    }
+    return (false);
+}
+
+template <typename Preconditioner>
+bool preconditionedSCBiCG(const SparseMat &A, const Preconditioner &M,
+        const Ref< VectorXcr const > b,
+        Ref <VectorXcr > x,
+        const int &max_iter, const Real &tol,
+        Real &errorn, int &iter_done) {
+
+    Real resid;
+    VectorXcr p, z, q;
+    Complex alpha, beta, rho, rho_1;
+
+    Real normb = b.norm( );
+    VectorXcr r = b - A*x;
+
+    if (normb == 0.0) normb = 1;
+
+    if ((resid = r.norm( ) / normb) <= tol) {
+        errorn = resid;
+        iter_done = 0;
+        return 0;
+    }
+
+    for (int i = 1; i <= max_iter; i++) {
+        z = M.solve(r);
+        rho = r.dot(z);
+
+        if (i == 1)  p = z;
+        else {
+            beta = rho / rho_1;
+            p = z + beta * p;
+        }
+
+        q = A*p;
+        alpha = rho / p.dot(q);
+
+        x += alpha * p;
+        r -= alpha * q;
+        std::cout << "resid\t" << resid << std::endl;
+        if ((resid = r.norm( ) / normb) <= tol) {
+            errorn = resid;
+            iter_done = i;
+            return 0;
+        }
+
+        rho_1 = rho;
+    }
+
+    errorn = resid;
+
+    return (false);
+}
+
+
+/** \internal Low-level conjugate gradient algorithm
+  * \param mat The matrix A
+  * \param rhs The right hand side vector b
+  * \param x On input and initial solution, on output the computed solution.
+  * \param precond A preconditioner being able to efficiently solve for an
+  *                approximation of Ax=b (regardless of b)
+  * \param iters On input the max number of iteration, on output the number of performed iterations.
+  * \param tol_error On input the tolerance error, on output an estimation of the relative error.
+  */
+template<typename Rhs, typename Dest, typename Preconditioner>
+EIGEN_DONT_INLINE
+void conjugateGradient(const SparseMat& mat, const Rhs& rhs, Dest& x,
+                        const Preconditioner& precond, int& iters,
+                        typename Dest::RealScalar& tol_error)
+{
+  using std::sqrt;
+  using std::abs;
+  typedef typename Dest::RealScalar RealScalar;
+  typedef typename Dest::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1> VectorType;
+
+  RealScalar tol = tol_error;
+  int maxIters = iters;
+
+  int n = mat.cols();
+
+  VectorType residual = rhs - mat.selfadjointView<Eigen::Upper>() * x; //initial residual
+
+  RealScalar rhsNorm2 = rhs.squaredNorm();
+  if(rhsNorm2 == 0)
+  {
+    x.setZero();
+    iters = 0;
+    tol_error = 0;
+    return;
+  }
+  RealScalar threshold = tol*tol*rhsNorm2;
+  RealScalar residualNorm2 = residual.squaredNorm();
+  if (residualNorm2 < threshold)
+  {
+    iters = 0;
+    tol_error = sqrt(residualNorm2 / rhsNorm2);
+    return;
+  }
+
+  VectorType p(n);
+  p = precond.solve(residual);      //initial search direction
+
+  VectorType z(n), tmp(n);
+  RealScalar absNew = numext::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
+  int i = 0;
+  while(i < maxIters)
+  {
+    tmp.noalias() = mat.selfadjointView<Eigen::Upper>() * p;              // the bottleneck of the algorithm
+
+    Scalar alpha = absNew / p.dot(tmp);   // the amount we travel on dir
+    x += alpha * p;                       // update solution
+    residual -= alpha * tmp;              // update residue
+
+    residualNorm2 = residual.squaredNorm();
+    if(residualNorm2 < threshold)
+      break;
+
+    z = precond.solve(residual);          // approximately solve for "A z = residual"
+
+    RealScalar absOld = absNew;
+    absNew = numext::real(residual.dot(z));     // update the absolute value of r
+    RealScalar beta = absNew / absOld;            // calculate the Gram-Schmidt value used to create the new search direction
+    p = z + beta * p;                             // update search direction
+    i++;
+  }
+  tol_error = sqrt(residualNorm2 / rhsNorm2);
+  iters = i;
+}
+
+// // Computes implicit
+// VectorXcr implicitDCInvBPhi (const SparseMat& D, const SparseMat& C,
+//                         const SparseMat& B, const SparseMat& MC,
+//                         const VectorXcr& Phi, Real& tol,
+//                         int& max_it) {
+//     int iter_done(0);
+//     Real errorn(0);
+//     VectorXcr b = B*Phi;
+//     VectorXcr y = VectorXcr::Zero(C.rows()) ; // = C^1*b;
+//     bicgstab(C, MC, b, y, max_it, tol, errorn, iter_done, false);
+//     //std::cout << "Temp " << errorn << std::endl;
+//     return  D*y;
+// }
+
+// Computes implicit
+VectorXcr implicitDCInvBPhi (const SparseMat& D, const SparseMat& C,
+                        const VectorXcr& ioms, const SparseMat& MC,
+                        const VectorXcr& Phi, Real& tol,
+                        int& max_it) {
+    int iter_done(0);
+    Real errorn(0);
+    VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
+    VectorXcr y = VectorXcr::Zero(C.rows()) ; // = C^1*b;
+    bicgstab(C, MC, b, y, max_it, tol, errorn, iter_done, false);
+    //std::cout << "Temp " << errorn << std::endl;
+    max_it = iter_done;
+    return  D*y;
+}
+
+// Computes implicit
+template <typename Preconditioner>
+VectorXcr implicitDCInvBPhi2 (const SparseMat& D, const SparseMat& C,
+                        const Ref<VectorXcr const> ioms, const Preconditioner& solver,
+                        const Ref<VectorXcr const> Phi, Real& tol,
+                        int& max_it) {
+
+    VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
+    VectorXcr y = VectorXcr::Zero(C.rows()) ; // = C^1*b;
+
+    // Home Made
+    //int iter_done(0);
+    //Real errorn(0);
+    //preconditionedBiCGStab(C, solver, b, y, max_it, tol, errorn, iter_done); //, false); // Jacobi M
+    //max_it = iter_done;
+
+    // Eigen BiCGStab
+    Eigen::BiCGSTAB<SparseMatrix<Complex> > BiCG;
+    BiCG.compute( C ); // TODO move this out of this loop!
+    y = BiCG.solve(b);
+    max_it = BiCG.iterations();
+    tol = BiCG.error();
+
+    // Direct
+/*
+    std::cout << "Computing LLT" << std::endl;
+    Eigen::SimplicialLLT<SparseMatrix<Complex>, Eigen::Upper, Eigen::AMDOrdering<int> >  LLT;
+    LLT.compute(C);
+    max_it = 1;
+    std::cout << "Computed LLT" << std::endl;
+    y = LLT.solve(b);
+*/
+
+    return  D*y;
+}
+
+// Computes implicit
+//template <typename Solver>
+template < typename Solver >
+inline VectorXcr implicitDCInvBPhi3 (const SparseMat& D, const Solver& solver,
+                        const Ref<VectorXcr const> ioms,
+                        const Ref<VectorXcr const> Phi, Real& tol,
+                        int& max_it) {
+    VectorXcr b = (ioms).asDiagonal() * (D.transpose()*Phi);
+    VectorXcr y = solver.solve(b);
+    max_it = 0;
+    //max_it = solver.iterations();
+    //errorn = solver.error();
+    return  D*y;
+}
+
+
+// // Simple extraction of indices in idx into reduceed array x1
+// void vmap( const Ref<VectorXcr const> x0, Ref<VectorXcr> x1, const std::vector<int>& idx ) {
+//     for (unsigned int ii=0; ii<idx.size(); ++ii) {
+//         x1(ii) = x0(idx[ii]);
+//     }
+// }
+
+// Simple extraction of indices in idx into reduceed array x1
+VectorXcr vmap( const Ref<VectorXcr const> x0, const std::vector<int>& idx ) {
+    VectorXcr x1 = VectorXcr::Zero( idx.size() );
+    for (unsigned int ii=0; ii<idx.size(); ++ii) {
+        x1(ii) = x0(idx[ii]);
+    }
+    return x1;
+}
+
+// reverse of above
+void ivmap( Ref<VectorXcr > x0, const Ref<VectorXcr const> x1, const std::vector<int>& idx ) {
+    for (unsigned int ii=0; ii<idx.size(); ++ii) {
+        x0(idx[ii]) = x1(ii);
+    }
+}
+
+
+// On Input
+// A = Matrix
+// B = Right hand side
+// X = initial guess, and solution
+// maxit = maximum Number of iterations
+// tol = error tolerance
+// On Output
+// X real solution vector
+// errorn = Real error norm
+template < typename CSolver >
+int implicitbicgstab(//const SparseMat& D,
+                     //const SparseMat& C,
+                     const Ref< Eigen::SparseMatrix<Complex> const > D,
+                     const std::vector<int>& idx,
+                     const Ref< VectorXcr const > ioms,
+                     const Ref< VectorXcr const > rhs,
+                     Ref <VectorXcr> phi,
+                     CSolver& solver,
+                     int &max_it, const Real &tol, Real &errorn, int &iter_done, ofstream& logio) {
+
+    logio << "using the preconditioned implicit solver" << std::endl;
+
+    Complex omega, rho, rho_1, alpha, beta;
+    Real    tol2;
+    int     iter, max_it2, max_it1;
+
+    // Look at reduced problem
+    VectorXcr rhs2 = vmap(rhs, idx);
+    VectorXcr phi2 = vmap(phi, idx);
+
+    // Determine size of system and init vectors
+    int n = idx.size();        // was phi.size();
+    VectorXcr r(n);
+    VectorXcr r_tld(n);
+    VectorXcr p(n);
+    VectorXcr s(n);
+    VectorXcr v = VectorXcr::Zero(n);
+    VectorXcr t = VectorXcr::Zero(n);
+
+//     TODO, refigure for implicit large system
+//     std::cout << "Start BiCGStab, memory needed: "
+//               <<  (sizeof(Complex)*(9+2)*n/(1024.*1024*1024)) << " [Gb]\n";
+
+    // Initialise
+    iter_done = 0;
+    Real eps = 1e-100;
+
+    Real bnrm2 = rhs.norm();
+    if (bnrm2 == 0) {
+        phi.setConstant(0.0);
+        errorn = 0;
+        std::cerr << "Trivial case of Ax = b, where b is 0\n";
+        return (0);
+    }
+
+    // If there is an initial guess
+    if ( phi.norm() ) {
+        tol2 = tol;
+        max_it2 = 50000;
+        //r = rhs - implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2);
+        //r = rhs - implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2);
+        r = rhs2 - vmap( implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2), idx );
+    } else {
+        r = vmap(rhs, idx);
+    }
+
+    jsw_timer timer;
+
+    errorn = r.norm() / bnrm2;
+    omega = 1.;
+    r_tld = r;
+    Real errornold = 1e14;
+    // Get down to business
+    for (iter=0; iter<max_it; ++iter) {
+
+        timer.begin();
+
+        rho = r_tld.dot(r);
+        if (abs(rho) < eps) {
+            ivmap( phi, phi2, idx );
+            return (0);
+        }
+
+        if (iter > 0) {
+            beta = (rho/rho_1) * (alpha/omega);
+            p = r.array() + beta*(p.array()-omega*v.array()).array();
+        } else {
+            p = r;
+        }
+
+        tol2 = tol;
+
+        max_it2 = 500000;
+        //v = implicitDCInvBPhi2(D, C, ioms, solver, p, tol2, max_it2);
+        ivmap(phi, p, idx);
+        v = vmap(implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2), idx);
+
+        alpha = rho / r_tld.dot(v);
+        s = r.array() - alpha*v.array();
+        errorn = s.norm()/bnrm2;
+
+        if (errorn < tol && iter > 1) {
+            phi2.array() += alpha*p.array();
+            ivmap( phi, phi2, idx );
+            return (0);
+        }
+
+        tol2 = tol;
+
+        max_it1 = 500000;
+        //t = implicitDCInvBPhi2(D, C, ioms, solver, s, tol2, max_it1);
+        //t = implicitDCInvBPhi3(D, solver, ioms, s, tol2, max_it1);
+        ivmap(phi, s, idx);
+        t = vmap(implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it1), idx);
+        omega = t.dot(s)  / t.dot(t);
+
+        r = s.array() - omega*t.array();
+        errorn = r.norm() / bnrm2;
+        iter_done = iter;
+
+        if (errorn <= tol) {
+            ivmap( phi, phi2, idx );
+            return (0);
+        }
+
+        if (abs(omega) < eps) {
+            ivmap( phi, phi2, idx );
+            return (0);
+        }
+
+        rho_1 = rho;
+
+        logio << "iteration " << std::setw(3) << iter
+              << "  errorn " << std::setw(6) << std::setprecision(4) << std::scientific << errorn
+              //<< "\timplicit iterations " << std::setw(5) << max_it1+max_it2
+              << "  time " << std::setw(6) << std::fixed << std::setprecision(2) << timer.end() << std::endl;
+
+        // Check to see how progress is going
+
+        if (errornold - errorn < 0) {
+            logio << "Irregular non-monotonic (negative) convergence. Recommend restart. \n";
+            ivmap( phi, phi2, idx );
+            return (2);
+        }
+
+        /*
+        if (errornold - errorn < 1e-14) {
+            logio << "not making any progress. Giving up\n";
+            return (1);
+        }
+        */
+
+        //std::cout << "|| p-s ||" << (alpha*p - omega*s).norm() << std::endl;
+
+        // only update phi if good things are happening
+        phi2.array() += alpha*p.array() + omega*s.array();
+        errornold = errorn;
+
+    }
+    ivmap( phi, phi2, idx );
+    return (0);
+}
+
+// On Input
+// A = Matrix
+// B = Right hand side
+// X = initial guess, and solution
+// maxit = maximum Number of iterations
+// tol = error tolerance
+// On Output
+// X real solution vector
+// errorn = Real error norm
+template < typename Solver >
+int implicitbicgstab_ei(const SparseMat&  D,
+                        const Ref< VectorXcr const > ioms,
+                        const Ref< VectorXcr const > rhs,
+                        Ref <VectorXcr> phi,
+                        Solver& solver,
+                        int &max_it, const Real &tol, Real &errorn, int &iter_done, ofstream& logio) {
+
+    logio << "using the preconditioned Eigen implicit solver" << std::endl;
+
+    Complex omega, rho, rho_1, alpha, beta;
+    Real tol2;
+    int  iter, max_it2,max_it1;
+
+    // Determine size of system and init vectors
+    int n = phi.size();
+    VectorXcr r(n);
+    VectorXcr r_tld(n);
+    VectorXcr p(n);
+    VectorXcr v(n);
+    VectorXcr s(n);
+    VectorXcr t(n);
+
+    // Initialise
+    iter_done = 0;
+    Real eps = 1e-100;
+
+    Real bnrm2 = rhs.norm();
+    if (bnrm2 == 0) {
+        phi.setConstant(0.0);
+        errorn = 0;
+        std::cerr << "Trivial case of Ax = b, where b is 0\n";
+        return (0);
+    }
+
+    // If there is an initial guess
+    if ( phi.norm() ) {
+        tol2 = tol;
+        max_it2 = 50000;
+        r = rhs - implicitDCInvBPhi3(D, solver, ioms, phi, tol2, max_it2);
+    } else {
+        r = rhs;
+    }
+
+    jsw_timer timer;
+
+    errorn = r.norm() / bnrm2;
+    omega = 1.;
+    r_tld = r;
+    Real errornold = 1e14;
+
+    // Get down to business
+    for (iter=0; iter<max_it; ++iter) {
+
+        timer.begin();
+
+        rho = r_tld.dot(r);
+        if (abs(rho) < eps) return (0);
+
+        if (iter > 0) {
+            beta = (rho/rho_1) * (alpha/omega);
+            p = r.array() + beta*(p.array()-omega*v.array()).array();
+        } else {
+            p = r;
+        }
+
+        tol2 = tol;
+        max_it2 = 500000;
+        v = implicitDCInvBPhi3(D, solver, ioms, p, tol2, max_it2);
+        max_it2 = 0; // solver.iterations();
+
+        alpha = rho / r_tld.dot(v);
+        s = r.array() - alpha*v.array();
+        errorn = s.norm()/bnrm2;
+
+        if (errorn < tol && iter > 1) {
+            phi.array() += alpha*p.array();
+            return (0);
+        }
+
+        tol2 = tol;
+        max_it1 = 500000;
+        t = implicitDCInvBPhi3(D, solver, ioms, s, tol2, max_it1);
+        max_it1 = 0; //solver.iterations();
+        omega = t.dot(s)  / t.dot(t);
+
+        r = s.array() - omega*t.array();
+        errorn = r.norm() / bnrm2;
+        iter_done = iter;
+
+        if (errorn <= tol ) return (0);
+        if (abs(omega) < eps) return (0);
+        rho_1 = rho;
+
+        logio << "iteration " << std::setw(4) << iter
+              << "\terrorn " << std::setw(6) << std::setprecision(4) << std::scientific << errorn
+              << "\timplicit iterations " << std::setw(5) << max_it1+max_it2
+              << "\ttime " << std::setw(10) << std::fixed << std::setprecision(2) << timer.end() << std::endl;
+
+        // Check to see how progress is going
+        if (errornold - errorn < 0) {
+            logio << "irregular (negative) convergence. Try again? \n";
+            return (2);
+        }
+
+        // only update phi if good things are happening
+        phi.array() += alpha*p.array() + omega*s.array();
+        errornold = errorn;
+
+    }
+    return (0);
+}
+
+
+// On Input
+// A = Matrix
+// B = Right hand side
+// X = initial guess, and solution
+// maxit = maximum Number of iterations
+// tol = error tolerance
+// On Output
+// X real solution vector
+// errorn = Real error norm
+int implicitbicgstab(const SparseMat& D,
+                     const SparseMat& C,
+                     const VectorXcr& ioms,
+                     const SparseMat& MC,
+                     Eigen::Ref< VectorXcr > rhs,
+                     VectorXcr& phi,
+                     int &max_it, Real &tol, Real &errorn, int &iter_done) {
+
+    Complex omega, rho, rho_1, alpha, beta;
+    Real errmin, tol2;
+    int  iter, max_it2;
+
+//     // Cholesky decomp
+//     SparseLLT<SparseMatrix<Complex>, Cholmod>
+//         CholC(SparseMatrix<Complex> (C.real()) );
+//     if(!CholC.succeeded()) {
+//         std::cerr << "decomposiiton failed\n";
+//         return EXIT_FAILURE;
+//     }
+
+    // Determine size of system and init vectors
+    int n = phi.size();
+    VectorXcr r(n);
+    VectorXcr r_tld(n);
+    VectorXcr p(n);
+    VectorXcr v(n);
+    //VectorXcr p_hat(n);
+    VectorXcr s(n);
+    //VectorXcr s_hat(n);
+    VectorXcr t(n);
+    VectorXcr xmin(n);
+
+//     TODO, refigure for implicit large system
+//     std::cout << "Start BiCGStab, memory needed: "
+//               <<  (sizeof(Complex)*(9+2)*n/(1024.*1024*1024)) << " [Gb]\n";
+
+    // Initialise
+    iter_done = 0;
+    v.setConstant(0.); // not necessary I don't think
+    t.setConstant(0.);
+    Real eps = 1e-100;
+
+    Real bnrm2 = rhs.norm();
+    if (bnrm2 == 0) {
+        phi.setConstant(0.0);
+        errorn = 0;
+        std::cerr << "Trivial case of Ax = b, where b is 0\n";
+        return (0);
+    }
+
+    // If there is an initial guess
+    if ( phi.norm() ) {
+        //r = rhs - A*phi;
+        tol2 = tol;
+        max_it2 = 50000;
+        std::cout << "Initial guess " << std::endl;
+        r = rhs - implicitDCInvBPhi(D, C, ioms, MC, phi, tol2, max_it2);
+        //r = rhs - implicitDCInvBPhi (D, C, B, CholC, phi, tol2, max_it2);
+    } else {
+        r = rhs;
+    }
+
+
+    errorn = r.norm() / bnrm2;
+    //std::cout << "Initial |r|  " << r.norm() << "\t" << errorn<< std::endl;
+    omega = 1.;
+    r_tld = r;
+    errmin = 1e30;
+    Real errornold = 1e6;
+    // Get down to business
+    for (iter=0; iter<max_it; ++iter) {
+
+        rho = r_tld.dot(r);
+        if (abs(rho) < eps) return (0);
+
+        if (iter > 0) {
+            beta = (rho/rho_1) * (alpha/omega);
+            p = r.array() + beta*(p.array()-omega*v.array()).array();
+        } else {
+            p = r;
+        }
+
+        // Use pseudo inverse to get approximate answer
+        //p_hat = p;
+        tol2  = std::max(1e-4*errorn, tol);
+        tol2 = tol;
+        max_it2 = 500000;
+        //v = A*p_hat;
+        v = implicitDCInvBPhi(D, C, ioms, MC, p, tol2, max_it2);
+        //v = implicitDCInvBPhi(D, C, B, CholC, p, tol2, max_it2);
+
+        alpha = rho / r_tld.dot(v);
+        s = r.array() - alpha*v.array();
+        errorn = s.norm()/bnrm2;
+
+        if (errorn < tol && iter > 1) {
+            phi.array() += alpha*p.array();
+            return (0);
+        }
+
+        // s_hat = M*s;
+        //tol2 = tol;
+        tol2  = std::max(1e-4*errorn, tol);
+        tol2 = tol;
+        max_it2 = 50000;
+        // t = A*s_hat;
+        t = implicitDCInvBPhi(D, C, ioms, MC, s, tol2, max_it2);
+        //t = implicitDCInvBPhi(D, C, B, CholC, s, tol2, max_it2);
+        omega = t.dot(s)  / t.dot(t);
+        phi.array() += alpha*p.array() + omega*s.array();
+        r = s.array() - omega*t.array();
+        errorn = r.norm() / bnrm2;
+        iter_done = iter;
+        if (errorn < errmin) {
+            // remember the model with the smallest norm
+            errmin = errorn;
+            xmin = phi;
+        }
+
+        if (errorn <= tol ) return (0);
+        if (abs(omega) < eps) return (0);
+        rho_1 = rho;
+
+        std::cout << "iteration " << std::setw(4) << iter << "\terrorn "  << std::setw(6) << std::scientific << errorn
+                  << "\timplicit iterations " << std::setw(5) << max_it2 << std::endl;
+        if (errornold - errorn < 1e-14) {
+            std::cout << "not making any progress. Giving up\n";
+            return (2);
+        }
+        errornold = errorn;
+
+    }
+    return (0);
+}
+
+
+//int bicgstab(const SparseMat &A, Eigen::SparseLU< Eigen::SparseMatrix<Complex, RowMajor> ,
+
+

src/CMakeLists.txt

@@ -0,0 +1,4 @@
+set (EMSCHUR3DSOURCE
+	${CMAKE_CURRENT_SOURCE_DIR}/EMSchur3DBase.cpp
+	PARENT_SCOPE
+)