EMSchur3D/include/EMSchur3D.h

/* 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 "bicgstab.h"

#ifdef HAVE_SUPERLU
#include "Eigen/SuperLUSupport"
#endif

#ifdef HAVE_UMFPACK
#include <Eigen/UmfPackSupport>
#endif

#ifdef HAVE_PASTIX
#include <Eigen/PaStiXSupport>
#endif

//#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) {
            stream << ob.Serialize()  << "\n"; // End of doc
            return stream;
        }

        //friend std::ostream &operator<<(std::ostream &stream,
        //        const EMSchur3D &ob);

        public:

        // ====================  LIFECYCLE     =======================

        /**
         * @copybrief LemmaObject::New()
         * @copydetails LemmaObject::New()
         */
        static std::shared_ptr< EMSchur3D > NewSP() {
            return std::make_shared< EMSchur3D<Solver> >( ctor_key() );
        }

        /** Default protected constructor, use New */
        explicit EMSchur3D ( const ctor_key& key ) : EMSchur3DBase( key ), CSolver( nullptr ) {

        }

        /** Locked DeDerializing constructor, use factory DeSerialize  method*/
        EMSchur3D (const YAML::Node& node, const ctor_key& key): EMSchur3DBase(node, key), CSolver( nullptr ) {
        }

        /** Default protected destructor, use Delete */
        virtual ~EMSchur3D () {
            // TODO delete arrays
        }

        /**
         *  Uses YAML to serialize this object.
         *  @return a YAML::Node
         */
        YAML::Node Serialize() const {
            YAML::Node node = EMSchur3DBase::Serialize();
            //node["NumberOfLayers"] = NumberOfLayers;
            node.SetTag( this->GetName() );
            return node;
        }

        /**
         *   Constructs an object from a YAML::Node.
         */
        static EMSchur3D* DeSerialize(const YAML::Node& node);


        // ====================  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( std::shared_ptr<DipoleSource> Source , const int& isource);

        /** Builds the solver for the C matrix */
        void BuildCDirectSolver(  );

        // ====================  ACCESS        =======================

        virtual std::string GetName() const {
            return this->CName;
        }

        // ====================  INQUIRY       =======================

        protected:

        // ====================  LIFECYCLE     =======================

        private:

        /** Copy constructor */
        EMSchur3D( const EMSchur3D& ) = delete;

        // ====================  DATA MEMBERS  =========================

        /** The templated solver for C */
        Solver*     CSolver;

        Eigen::SparseMatrix<Complex>  Csym;

        static constexpr auto CName = "EMSchur3D";


    }; // -----  end of class  EMSchur3D  -----


    ////////////////////////////////////////////////////////////////////////////////////////
    //                    Implimentation and Specialisations                              //
    ////////////////////////////////////////////////////////////////////////////////////////

    //--------------------------------------------------------------------------------------
    //       Class:  EMSchur3D
    //      Method:  SolveSource
    //--------------------------------------------------------------------------------------
    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) {
           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();
        std::shared_ptr< FieldPoints > dpoint = FieldPoints::NewSP();

        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

        // 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
        std::cout << "Filling source terms" << std::endl;
        FillSourceTerms(ms, Se, E0, dpoint, Omegas[iw]);
        std::cout << "Done source terms" << std::endl;

        /////////////////////////////////////////////////
        // LOG File
        std::string logfile (ResFile);
        logfile += to_string(isource) + std::string(".log");
        ofstream logio(logfile.c_str());

        std::cout << "just logging, TODO fix source" << 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);
        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
        //CSolver[iw].setMaxIterations(750);
        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());
        VectorXcr Error = ((Cvec[iw].selfadjointView<Eigen::Lower>()*A).array() - Se.array());
        logio << "|| Div(A) || = " << ADiv.norm()
              << "\tInital solution error="<<   Error.norm()  // Iteritive info
//              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
              << "\ttime " << timer.end() / 60. << " [m]   "
//              <<  CSolver[iw].iterations() << "  iterations"
              << 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-20;
        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()/60. << " [m]" << 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());
        Error = ((Cvec[iw].selfadjointView<Eigen::Lower>()*A).array() + E.array() - Se.array());

        //      << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
        //      << "\ttime " << timer.end() / 60. << " [m]   "
        //      <<  CSolver[iw].iterations() << "  iterations"


        logio << "|| Div(A) || = " << ADiv2.norm()
              << "\tSolution error="<<   Error.norm()  // Iteritive info
//              << "\tSolver reported error="<<   CSolver[iw].error()  // Iteritive info
              << "\ttime " << timer.end() / 60. << " [m]   "
//              <<  CSolver[iw].iterations() << "  iterations"
              << std::endl << std::endl;

        logio << std::fixed << std::setprecision(2) << "\ttime " << timer.end()/60. << "\ttotal time " << timer2.end()/60. << 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_SUPERLU
    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 */
            // Recommended values for symmetric mode
            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 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 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); // OK for symmetric complex systems with real and imaginary positive definite parts.
            //                                  // but our imaginary part is negative definite
            //http://www.ams.org/journals/mcom/1998-67-224/S0025-5718-98-00978-8/S0025-5718-98-00978-8.pdf

            /*  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;
        }
    }

#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::ColMajor> > >::BuildCDirectSolver() {
        CSolver = new Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor> > [Omegas.size()];
        for (int iw=0; iw<Omegas.size(); ++iw) {
            jsw_timer timer;
            timer.begin();
            /*  Complex system */
            std::cout << "UmfPackLU pattern analyzing C_" << iw << ",";
            std::cout.flush();

            // Doesn't work, seg faults in solve
            //Eigen::SparseMatrix<Complex> Csym = Cvec[iw].selfadjointView< Eigen::Lower >();

            // Compiler errors get thrown with the view setting, good performance if LOWER and UPPER are set
            // CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() ); // Compiler error
            CSolver[iw].analyzePattern( Cvec[iw] ); // requires LOWER=1 UPPER=1, double memory
            // CSolver[iw].analyzePattern( Csym );     // seg faults
            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;

            // factorize
            timer.begin();
            std::cout << "UmfPackLU factorising C_" << iw << ", ";
            std::cout.flush();
            // CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
            CSolver[iw].factorize( Cvec[iw] );
            // CSolver[iw].factorize( Csym );
            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
        }
    }
#endif

//     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>();
            CSolver[iw].preconditioner().setDroptol(1e-6);       //1e-5);      // 1e-12
            CSolver[iw].preconditioner().setFillfactor(5e1);     // 1e2
            jsw_timer timer;
            timer.begin();
            /*  Complex system */
            std::cout << "BiCGSTAB(ILU) pattern analyzing C_" << iw << ",";
            std::cout.flush();
            CSolver[iw].analyzePattern( Csym );
            //CSolver[iw].analyzePattern( Cvec[iw]);
            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 );
            //CSolver[iw].factorize( Cvec[iw] );
            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::LeastSquaresConjugateGradient< Eigen::SparseMatrix<Complex, Eigen::ColMajor> > > ::BuildCDirectSolver() {
        CSolver = new Eigen::LeastSquaresConjugateGradient< 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 << "LeastSquaresConjugateGradient 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 << "LeastSquaresConjugateGradient 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::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::IncompleteCholesky<Complex> > > ::BuildCDirectSolver() {
        CSolver = new Eigen::ConjugateGradient<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, Eigen::Lower, Eigen::IncompleteCholesky<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 << "ConjugateGradient<IncompleteCholesky> 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<IncompleteCholesky> 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;
//         }
//     }
//

#ifdef HAVE_PASTIX
    template<>
    void EMSchur3D<   Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > > ::BuildCDirectSolver() {
        CSolver = new Eigen::PastixLU<Eigen::SparseMatrix<Complex, Eigen::ColMajor>, true > [Omegas.size()];

        for (int iw=0; iw<Omegas.size(); ++iw) {
            jsw_timer timer;
            timer.begin();
            /*  Complex system */
            std::cout << "PastixLU pattern analyzing C_" << iw << ",";
            std::cout.flush();
            CSolver[iw].analyzePattern( Cvec[iw].selfadjointView< Eigen::Lower>() );
            //CSolver[iw].analyzePattern( Cvec[iw] );
            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;

            // factorize
            timer.begin();
            std::cout << "PastixLU factorising C_" << iw << ", ";
            std::cout.flush();
            CSolver[iw].factorize( Cvec[iw].selfadjointView< Eigen::Lower>() );
            //CSolver[iw].factorize( Cvec[iw] );
            std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
        }
    }
#endif


}		// -----  end of Lemma  name  -----

#endif   // ----- #ifndef EMSCHUR3D_INC  -----