mirror of
				https://github.com/paboyle/Grid.git
				synced 2025-11-04 05:54:32 +00:00 
			
		
		
		
	added DistilPar-module and cleaned up some code
This commit is contained in:
		@@ -70,6 +70,7 @@
 | 
			
		||||
#include <Hadrons/Modules/MDistil/Noises.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MDistil/LapEvec.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MDistil/Perambulator.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MDistil/DistilPar.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MAction/ScaledDWF.hpp>
 | 
			
		||||
 
 | 
			
		||||
@@ -63,42 +63,14 @@ BEGIN_MODULE_NAMESPACE(MDistil)
 | 
			
		||||
 | 
			
		||||
struct DistilParameters: Serializable {
 | 
			
		||||
    GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParameters,
 | 
			
		||||
                                    int, nvec,
 | 
			
		||||
                                    int, nnoise,
 | 
			
		||||
                                    int, tsrc,
 | 
			
		||||
                                    std::string, TI,
 | 
			
		||||
                                    std::string, LI,
 | 
			
		||||
                                    std::string, SI )
 | 
			
		||||
    DistilParameters() = default;
 | 
			
		||||
    template <class ReaderClass> DistilParameters(Reader<ReaderClass>& Reader){read(Reader,"Distil",*this);}
 | 
			
		||||
    
 | 
			
		||||
    // Numeric parameter is allowed to be empty (in which case it = Default),
 | 
			
		||||
    // but assert during setup() if specified but not numeric
 | 
			
		||||
    
 | 
			
		||||
    static int ParameterDefault( const std::string & s, int Default, bool bCalledFromSetup )
 | 
			
		||||
    {
 | 
			
		||||
        int i = Default;
 | 
			
		||||
        if( s.length() > 0 ) {
 | 
			
		||||
            std::istringstream ss( s );
 | 
			
		||||
            ss >> i;
 | 
			
		||||
            if( bCalledFromSetup )
 | 
			
		||||
                assert( !ss.fail() && "Parameter should either be empty or integer" );
 | 
			
		||||
        }
 | 
			
		||||
        return i;
 | 
			
		||||
    }
 | 
			
		||||
                                    int, TI,
 | 
			
		||||
                                    int, LI,
 | 
			
		||||
                                    int, SI )
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
#define DISTIL_PARAMETERS_DEFINE( inSetup ) \
 | 
			
		||||
const int Nt{env().getDim(Tdir)}; \
 | 
			
		||||
const int nvec{par().nvec}; \
 | 
			
		||||
const int nnoise{par().Distil.nnoise}; \
 | 
			
		||||
const int tsrc{par().Distil.tsrc}; \
 | 
			
		||||
const int TI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.TI, Nt,   inSetup)}; \
 | 
			
		||||
const int LI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.LI, nvec, inSetup)}; \
 | 
			
		||||
const int SI{Hadrons::MDistil::DistilParameters::ParameterDefault(par().Distil.SI, Ns,   inSetup)}; \
 | 
			
		||||
const bool full_tdil{ TI == Nt }; \
 | 
			
		||||
const bool exact_distillation{ full_tdil && LI == nvec }; \
 | 
			
		||||
const int Nt_inv{ full_tdil ? 1 : TI }
 | 
			
		||||
 | 
			
		||||
/******************************************************************************
 | 
			
		||||
 Make a lower dimensional grid in preparation for local slice operations
 | 
			
		||||
 ******************************************************************************/
 | 
			
		||||
 
 | 
			
		||||
							
								
								
									
										7
									
								
								Hadrons/Modules/MDistil/DistilPar.cc
									
									
									
									
									
										Normal file
									
								
							
							
						
						
									
										7
									
								
								Hadrons/Modules/MDistil/DistilPar.cc
									
									
									
									
									
										Normal file
									
								
							@@ -0,0 +1,7 @@
 | 
			
		||||
#include <Hadrons/Modules/MDistil/DistilPar.hpp>
 | 
			
		||||
 | 
			
		||||
using namespace Grid;
 | 
			
		||||
using namespace Hadrons;
 | 
			
		||||
using namespace MDistil;
 | 
			
		||||
 | 
			
		||||
template class Grid::Hadrons::MDistil::TDistilPar<FIMPL>;
 | 
			
		||||
							
								
								
									
										96
									
								
								Hadrons/Modules/MDistil/DistilPar.hpp
									
									
									
									
									
										Normal file
									
								
							
							
						
						
									
										96
									
								
								Hadrons/Modules/MDistil/DistilPar.hpp
									
									
									
									
									
										Normal file
									
								
							@@ -0,0 +1,96 @@
 | 
			
		||||
#ifndef Hadrons_MDistil_DistilPar_hpp_
 | 
			
		||||
#define Hadrons_MDistil_DistilPar_hpp_
 | 
			
		||||
 | 
			
		||||
#include <Hadrons/Global.hpp>
 | 
			
		||||
#include <Hadrons/Module.hpp>
 | 
			
		||||
#include <Hadrons/ModuleFactory.hpp>
 | 
			
		||||
#include <Hadrons/Modules/MDistil/DistilCommon.hpp>
 | 
			
		||||
 | 
			
		||||
BEGIN_HADRONS_NAMESPACE
 | 
			
		||||
 | 
			
		||||
/******************************************************************************
 | 
			
		||||
 *                         DistilPar                                 *
 | 
			
		||||
 ******************************************************************************/
 | 
			
		||||
BEGIN_MODULE_NAMESPACE(MDistil)
 | 
			
		||||
 | 
			
		||||
 | 
			
		||||
class DistilParPar: Serializable
 | 
			
		||||
{
 | 
			
		||||
public:
 | 
			
		||||
    GRID_SERIALIZABLE_CLASS_MEMBERS(DistilParPar,
 | 
			
		||||
                                    int, nnoise,
 | 
			
		||||
                                    int, tsrc,
 | 
			
		||||
                                    int, TI,
 | 
			
		||||
                                    int, LI,
 | 
			
		||||
                                    int, SI );
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
class TDistilPar: public Module<DistilParPar>
 | 
			
		||||
{
 | 
			
		||||
public:
 | 
			
		||||
    // constructor
 | 
			
		||||
    TDistilPar(const std::string name);
 | 
			
		||||
    // destructor
 | 
			
		||||
    virtual ~TDistilPar(void) {};
 | 
			
		||||
    // dependency relation
 | 
			
		||||
    virtual std::vector<std::string> getInput(void);
 | 
			
		||||
    virtual std::vector<std::string> getOutput(void);
 | 
			
		||||
    // setup
 | 
			
		||||
    virtual void setup(void);
 | 
			
		||||
    // execution
 | 
			
		||||
    virtual void execute(void);
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
MODULE_REGISTER_TMP(DistilPar, TDistilPar<FIMPL>, MDistil);
 | 
			
		||||
 | 
			
		||||
/******************************************************************************
 | 
			
		||||
 *                 TDistilPar implementation                             *
 | 
			
		||||
 ******************************************************************************/
 | 
			
		||||
// constructor /////////////////////////////////////////////////////////////////
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
TDistilPar<FImpl>::TDistilPar(const std::string name)
 | 
			
		||||
: Module<DistilParPar>(name)
 | 
			
		||||
{}
 | 
			
		||||
 | 
			
		||||
// dependencies/products ///////////////////////////////////////////////////////
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
std::vector<std::string> TDistilPar<FImpl>::getInput(void)
 | 
			
		||||
{
 | 
			
		||||
    std::vector<std::string> in;
 | 
			
		||||
    
 | 
			
		||||
    return in;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
std::vector<std::string> TDistilPar<FImpl>::getOutput(void)
 | 
			
		||||
{
 | 
			
		||||
    std::vector<std::string> out = {getName()};
 | 
			
		||||
    
 | 
			
		||||
    return out;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
// setup ///////////////////////////////////////////////////////////////////////
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
void TDistilPar<FImpl>::setup(void)
 | 
			
		||||
{
 | 
			
		||||
    
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
// execution ///////////////////////////////////////////////////////////////////
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
void TDistilPar<FImpl>::execute(void)
 | 
			
		||||
{
 | 
			
		||||
    Hadrons::MDistil::DistilParameters &out = envGet(Hadrons::MDistil::DistilParameters, getName());
 | 
			
		||||
    out.nnoise=par().nnoise;
 | 
			
		||||
    out.tsrc=par().tsrc;
 | 
			
		||||
    out.TI=par().TI;
 | 
			
		||||
    out.LI=par().LI;
 | 
			
		||||
    out.SI=par().SI;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
END_MODULE_NAMESPACE
 | 
			
		||||
 | 
			
		||||
END_HADRONS_NAMESPACE
 | 
			
		||||
 | 
			
		||||
#endif // Hadrons_MDistil_DistilPar_hpp_
 | 
			
		||||
@@ -47,11 +47,9 @@ public:
 | 
			
		||||
                                    std::string, noise,
 | 
			
		||||
                                    std::string, perambulator,
 | 
			
		||||
                                    std::string, lapevec,
 | 
			
		||||
                                    std::string, source,
 | 
			
		||||
                                    std::string, sink,
 | 
			
		||||
                                    int, tsrc,
 | 
			
		||||
                                    std::string, nvec,
 | 
			
		||||
                                    std::string, TI)
 | 
			
		||||
                                    std::string, rho,
 | 
			
		||||
                                    std::string, phi,
 | 
			
		||||
                                    MDistil::DistilParameters, DistilPar);
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
@@ -80,10 +78,10 @@ public:
 | 
			
		||||
    std::string PerambulatorName;
 | 
			
		||||
    std::string NoiseVectorName;
 | 
			
		||||
    std::string LapEvecName;
 | 
			
		||||
    bool bMakeSource;
 | 
			
		||||
    bool bMakeSink;
 | 
			
		||||
    std::string SourceName;
 | 
			
		||||
    std::string SinkName;
 | 
			
		||||
    bool bMakeRho;
 | 
			
		||||
    bool bMakePhi;
 | 
			
		||||
    std::string RhoName;
 | 
			
		||||
    std::string PhiName;
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
MODULE_REGISTER_TMP(DistilVectors, TDistilVectors<FIMPL>, MDistil);
 | 
			
		||||
@@ -131,24 +129,24 @@ std::vector<std::string> TDistilVectors<FImpl>::getInput(void)
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
std::vector<std::string> TDistilVectors<FImpl>::getOutput(void)
 | 
			
		||||
{
 | 
			
		||||
    SourceName  = par().source;
 | 
			
		||||
    SinkName    = par().sink;
 | 
			
		||||
    bMakeSource = ( SourceName.size() > 0 );
 | 
			
		||||
    bMakeSink   = (   SinkName.size() > 0 );
 | 
			
		||||
    if (!bMakeSource && !bMakeSink)
 | 
			
		||||
    RhoName  = par().rho;
 | 
			
		||||
    PhiName  = par().phi;
 | 
			
		||||
    bMakeRho = ( RhoName.size() > 0 );
 | 
			
		||||
    bMakePhi = ( PhiName.size() > 0 );
 | 
			
		||||
    if (!bMakeRho && !bMakePhi)
 | 
			
		||||
    {
 | 
			
		||||
        SourceName = getName();
 | 
			
		||||
        SinkName   = SourceName;
 | 
			
		||||
        SourceName.append("_rho");
 | 
			
		||||
        SinkName.append("_phi");
 | 
			
		||||
        bMakeSource = true;
 | 
			
		||||
        bMakeSink   = true;
 | 
			
		||||
        RhoName = getName();
 | 
			
		||||
        PhiName = RhoName;
 | 
			
		||||
        RhoName.append("_rho");
 | 
			
		||||
        PhiName.append("_phi");
 | 
			
		||||
        bMakeRho = true;
 | 
			
		||||
        bMakePhi = true;
 | 
			
		||||
    }
 | 
			
		||||
    std::vector<std::string> out;
 | 
			
		||||
    if (bMakeSource)
 | 
			
		||||
        out.push_back(SourceName);
 | 
			
		||||
    if (bMakeSink)
 | 
			
		||||
        out.push_back(SinkName);
 | 
			
		||||
    if (bMakeRho)
 | 
			
		||||
        out.push_back(RhoName);
 | 
			
		||||
    if (bMakePhi)
 | 
			
		||||
        out.push_back(PhiName);
 | 
			
		||||
    return out;
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
@@ -157,28 +155,25 @@ template <typename FImpl>
 | 
			
		||||
void TDistilVectors<FImpl>::setup(void)
 | 
			
		||||
{
 | 
			
		||||
    Cleanup();
 | 
			
		||||
    auto &noise        = envGet(NoiseTensor, NoiseVectorName);
 | 
			
		||||
    auto &noise        = envGet(NoiseTensor,  NoiseVectorName);
 | 
			
		||||
    auto &perambulator = envGet(PerambTensor, PerambulatorName);
 | 
			
		||||
    
 | 
			
		||||
    // We expect the perambulator to have been created with these indices
 | 
			
		||||
    assert( perambulator.ValidateIndexNames() && "Perambulator index names bad" );
 | 
			
		||||
    
 | 
			
		||||
    const int Nt{ env().getDim(Tdir) };
 | 
			
		||||
    assert( Nt == static_cast<int>( perambulator.tensor.dimension(0) ) && "PerambTensor time dimensionality bad" );
 | 
			
		||||
    const int LI{ static_cast<int>( perambulator.tensor.dimension(2) ) };
 | 
			
		||||
    const int SI{ static_cast<int>( perambulator.tensor.dimension(5) ) };
 | 
			
		||||
    const int Nt_inv{ static_cast<int>( perambulator.tensor.dimension(4) ) };
 | 
			
		||||
    const int nnoise{ static_cast<int>( perambulator.tensor.dimension(3) ) };
 | 
			
		||||
    assert( nnoise >= static_cast<int>( noise.tensor.dimension(0) ) && "Not enough noise vectors for perambulator" );
 | 
			
		||||
    // Nvec defaults to what's in the perambulator unless overriden
 | 
			
		||||
    const int nvec_per{ static_cast<int>( perambulator.tensor.dimension(1) ) };
 | 
			
		||||
    const int nvec{Hadrons::MDistil::DistilParameters::ParameterDefault(par().nvec, nvec_per, true) };
 | 
			
		||||
    assert( nvec <= nvec_per && "Not enough distillation sub-space vectors" );
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
    
 | 
			
		||||
    if (bMakeSource)
 | 
			
		||||
        envCreate(std::vector<FermionField>, SourceName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
 | 
			
		||||
    if (bMakeSink)
 | 
			
		||||
        envCreate(std::vector<FermionField>,   SinkName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
 | 
			
		||||
    if (bMakeRho)
 | 
			
		||||
        envCreate(std::vector<FermionField>, RhoName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
 | 
			
		||||
    if (bMakePhi)
 | 
			
		||||
        envCreate(std::vector<FermionField>,   PhiName, 1, nnoise*LI*SI*Nt_inv, envGetGrid(FermionField));
 | 
			
		||||
    
 | 
			
		||||
    grid4d = env().getGrid();
 | 
			
		||||
    Coordinate latt_size   = GridDefaultLatt();
 | 
			
		||||
@@ -189,10 +184,10 @@ void TDistilVectors<FImpl>::setup(void)
 | 
			
		||||
    mpi_layout[Nd-1] = 1;
 | 
			
		||||
    grid3d = MakeLowerDimGrid(grid4d);
 | 
			
		||||
    
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "tmp2",1,LatticeSpinColourVector(grid4d));
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "sink_tslice",1,LatticeSpinColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "source4d",1,LatticeSpinColourVector(grid4d));
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "source3d",1,LatticeSpinColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "source3d_nospin",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "sink3d",1,LatticeSpinColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
@@ -216,50 +211,49 @@ void TDistilVectors<FImpl>::execute(void)
 | 
			
		||||
    auto &perambulator = envGet(PerambTensor, PerambulatorName);
 | 
			
		||||
    auto &epack        = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, LapEvecName);
 | 
			
		||||
    
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, tmp2);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, tmp3d);
 | 
			
		||||
    envGetTmp(LatticeColourVector,     tmp3d_nospin);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, sink_tslice);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, source4d);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, source3d);
 | 
			
		||||
    envGetTmp(LatticeColourVector,     source3d_nospin);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, sink3d);
 | 
			
		||||
    envGetTmp(LatticeColourVector,     evec3d);
 | 
			
		||||
    
 | 
			
		||||
    const int Ntlocal{ grid4d->LocalDimensions()[3] };
 | 
			
		||||
    const int Ntfirst{ grid4d->LocalStarts()[3] };
 | 
			
		||||
    
 | 
			
		||||
    const int Nt{ env().getDim(Tdir) };
 | 
			
		||||
    const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, false ) };
 | 
			
		||||
    const int LI{ static_cast<int>( perambulator.tensor.dimension(2) ) };
 | 
			
		||||
    const int SI{ static_cast<int>( perambulator.tensor.dimension(5) ) };
 | 
			
		||||
    const int Nt_inv{ static_cast<int>( perambulator.tensor.dimension(4) ) };
 | 
			
		||||
    const int nnoise{ static_cast<int>( perambulator.tensor.dimension(3) ) };
 | 
			
		||||
    // Nvec defaults to what's in the perambulator unless overriden
 | 
			
		||||
    const int nvec{Hadrons::MDistil::DistilParameters::ParameterDefault(par().nvec, static_cast<int>( perambulator.tensor.dimension(1) ), false)};
 | 
			
		||||
    const int tsrc{ par().tsrc };
 | 
			
		||||
    const bool full_tdil{ TI==Nt };
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int tsrc{par().DistilPar.tsrc}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
    
 | 
			
		||||
    int vecindex;
 | 
			
		||||
    int t_inv;
 | 
			
		||||
    if (bMakeSource)
 | 
			
		||||
    if (bMakeRho)
 | 
			
		||||
    {
 | 
			
		||||
        auto &rho = envGet(std::vector<FermionField>, SourceName);
 | 
			
		||||
        auto &rho = envGet(std::vector<FermionField>, RhoName);
 | 
			
		||||
        for (int inoise = 0; inoise < nnoise; inoise++) {
 | 
			
		||||
            for (int dk = 0; dk < LI; dk++) {
 | 
			
		||||
                for (int dt = 0; dt < Nt_inv; dt++) {
 | 
			
		||||
                    for (int ds = 0; ds < SI; ds++) {
 | 
			
		||||
                        vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt;
 | 
			
		||||
                        rho[vecindex] = 0;
 | 
			
		||||
                        tmp3d_nospin = 0;
 | 
			
		||||
                        source3d_nospin = 0;
 | 
			
		||||
                        for (int it = dt; it < Nt; it += TI){
 | 
			
		||||
                            if (full_tdil) t_inv = tsrc; else t_inv = it;
 | 
			
		||||
                            if (t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal) {
 | 
			
		||||
                                for (int ik = dk; ik < nvec; ik += LI){
 | 
			
		||||
                                    for (int is = ds; is < Ns; is += SI){
 | 
			
		||||
                                        ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
 | 
			
		||||
                                        tmp3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
 | 
			
		||||
                                        tmp3d=0;
 | 
			
		||||
                                        pokeSpin(tmp3d,tmp3d_nospin,is);
 | 
			
		||||
                                        tmp2=0;
 | 
			
		||||
                                        InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Tdir);
 | 
			
		||||
                                        rho[vecindex] += tmp2;
 | 
			
		||||
                                        source3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
 | 
			
		||||
                                        source3d=0;
 | 
			
		||||
                                        pokeSpin(source3d,source3d_nospin,is);
 | 
			
		||||
                                        source4d=0;
 | 
			
		||||
                                        InsertSliceLocal(source3d,source4d,0,t_inv-Ntfirst,Tdir);
 | 
			
		||||
                                        rho[vecindex] += source4d;
 | 
			
		||||
                                    }
 | 
			
		||||
                                }
 | 
			
		||||
                            }
 | 
			
		||||
@@ -269,8 +263,8 @@ void TDistilVectors<FImpl>::execute(void)
 | 
			
		||||
            }
 | 
			
		||||
        }
 | 
			
		||||
    }
 | 
			
		||||
    if (bMakeSink) {
 | 
			
		||||
        auto &phi = envGet(std::vector<FermionField>, SinkName);
 | 
			
		||||
    if (bMakePhi) {
 | 
			
		||||
        auto &phi = envGet(std::vector<FermionField>, PhiName);
 | 
			
		||||
        for (int inoise = 0; inoise < nnoise; inoise++) {
 | 
			
		||||
            for (int dk = 0; dk < LI; dk++) {
 | 
			
		||||
                for (int dt = 0; dt < Nt_inv; dt++) {
 | 
			
		||||
@@ -278,12 +272,12 @@ void TDistilVectors<FImpl>::execute(void)
 | 
			
		||||
                        vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * SI*dt;
 | 
			
		||||
                        phi[vecindex] = 0;
 | 
			
		||||
                        for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) {
 | 
			
		||||
                            sink_tslice=0;
 | 
			
		||||
                            sink3d=0;
 | 
			
		||||
                            for (int ivec = 0; ivec < nvec; ivec++) {
 | 
			
		||||
                                ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
 | 
			
		||||
                                sink_tslice += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds);
 | 
			
		||||
                                sink3d += evec3d * perambulator.tensor(t, ivec, dk, inoise,dt,ds);
 | 
			
		||||
                            }
 | 
			
		||||
                            InsertSliceLocal(sink_tslice,phi[vecindex],0,t-Ntfirst,Tdir);
 | 
			
		||||
                            InsertSliceLocal(sink3d,phi[vecindex],0,t-Ntfirst,Tdir);
 | 
			
		||||
                        }
 | 
			
		||||
                    }
 | 
			
		||||
                }
 | 
			
		||||
 
 | 
			
		||||
@@ -43,10 +43,7 @@ class NoisesPar: Serializable
 | 
			
		||||
{
 | 
			
		||||
public:
 | 
			
		||||
    GRID_SERIALIZABLE_CLASS_MEMBERS(NoisesPar,
 | 
			
		||||
                                    int, nnoise,
 | 
			
		||||
                                    int, nvec,
 | 
			
		||||
                                    std::string, TI,
 | 
			
		||||
                                    std::string, LI,
 | 
			
		||||
                                    MDistil::DistilParameters, DistilPar,
 | 
			
		||||
                                    std::string, NoiseFileName)
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
@@ -96,9 +93,9 @@ std::vector<std::string> TNoises<FImpl>::getOutput(void)
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
void TNoises<FImpl>::setup(void)
 | 
			
		||||
{
 | 
			
		||||
    const int Nt{env().getDim(Tdir)};
 | 
			
		||||
    const int nnoise{par().nnoise};
 | 
			
		||||
    const int nvec{par().nvec};
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    envCreate(NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
@@ -107,12 +104,12 @@ template <typename FImpl>
 | 
			
		||||
void TNoises<FImpl>::execute(void)
 | 
			
		||||
{
 | 
			
		||||
    const int Nt{env().getDim(Tdir)};
 | 
			
		||||
    const int nnoise{par().nnoise};
 | 
			
		||||
    const int nvec{par().nvec};
 | 
			
		||||
    const int TI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().TI, Nt, false) };
 | 
			
		||||
    const int LI{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI, nvec, false) };
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; \
 | 
			
		||||
    const bool exact_distillation{ full_tdil && LI == nvec }; \
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise};
 | 
			
		||||
    const int nvec{par().DistilPar.nvec};
 | 
			
		||||
    const int TI{par().DistilPar.TI};
 | 
			
		||||
    const int LI{par().DistilPar.LI};
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const bool exact_distillation{ full_tdil && LI == nvec }; 
 | 
			
		||||
    
 | 
			
		||||
    // We use our own seeds so we can specify different noises per quark
 | 
			
		||||
    Real rn;
 | 
			
		||||
 
 | 
			
		||||
@@ -51,10 +51,9 @@ public:
 | 
			
		||||
                                    std::string, eigenPack,
 | 
			
		||||
                                    std::string, PerambFileName,
 | 
			
		||||
                                    std::string, solve,
 | 
			
		||||
                                    int, nvec,
 | 
			
		||||
                                    std::string, nvec_reduced,
 | 
			
		||||
                                    std::string, LI_reduced,
 | 
			
		||||
                                    DistilParameters, Distil);
 | 
			
		||||
                                    MDistil::DistilParameters, DistilPar);
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
@@ -116,16 +115,23 @@ template <typename FImpl>
 | 
			
		||||
void TPerambFromSolve<FImpl>::setup(void)
 | 
			
		||||
{
 | 
			
		||||
    Cleanup();
 | 
			
		||||
    DISTIL_PARAMETERS_DEFINE( true );
 | 
			
		||||
    const int nvec_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().nvec_reduced, nvec, true) };
 | 
			
		||||
    const int LI_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI_reduced, LI, true) };
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
    const int nvec_reduced{ par().nvec_reduced.empty() ? nvec:std::stoi(par().nvec_reduced)};
 | 
			
		||||
    const int LI_reduced{ par().LI_reduced.empty() ? LI:std::stoi(par().LI_reduced)};
 | 
			
		||||
    grid4d = env().getGrid();
 | 
			
		||||
    grid3d = MakeLowerDimGrid(grid4d);
 | 
			
		||||
    envCreate(PerambTensor, getName(), 1, Nt,nvec_reduced,LI_reduced,nnoise,Nt_inv,SI);
 | 
			
		||||
    envCreate(NoiseTensor, getName() + "_noise", 1, nnoise, Nt, nvec, Ns );
 | 
			
		||||
    envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "result3d_nospin",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmpLat(LatticeColourVector, "result_nospin");
 | 
			
		||||
    envTmpLat(LatticeColourVector, "result4d_nospin");
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
@@ -146,28 +152,42 @@ void TPerambFromSolve<FImpl>::execute(void)
 | 
			
		||||
    GridCartesian * grid4d = env().getGrid();
 | 
			
		||||
    const int Ntlocal{grid4d->LocalDimensions()[3]};
 | 
			
		||||
    const int Ntfirst{grid4d->LocalStarts()[3]};
 | 
			
		||||
    DISTIL_PARAMETERS_DEFINE( false );
 | 
			
		||||
    const int nvec_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().nvec_reduced, nvec, false) };
 | 
			
		||||
    const int LI_reduced{ Hadrons::MDistil::DistilParameters::ParameterDefault( par().LI_reduced, LI, false) };
 | 
			
		||||
    auto &perambulator = envGet(PerambTensor, getName());
 | 
			
		||||
    auto &solve       = envGet(std::vector<FermionField>, par().solve);
 | 
			
		||||
    auto &epack   = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
    const int nvec_reduced{ par().nvec_reduced.empty() ? nvec:std::stoi(par().nvec_reduced)};
 | 
			
		||||
    const int LI_reduced{ par().LI_reduced.empty() ? LI:std::stoi(par().LI_reduced)};
 | 
			
		||||
    auto &perambulator  = envGet(PerambTensor, getName());
 | 
			
		||||
    auto &solve         = envGet(std::vector<FermionField>, par().solve);
 | 
			
		||||
    auto &epack         = envGet(Grid::Hadrons::EigenPack<LatticeColourVector>, par().eigenPack);
 | 
			
		||||
    
 | 
			
		||||
    envGetTmp(LatticeColourVector, result_nospin);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result_3d);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result4d_nospin);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result3d_nospin);
 | 
			
		||||
    envGetTmp(LatticeColourVector, evec3d);
 | 
			
		||||
    
 | 
			
		||||
    for (int inoise = 0; inoise < nnoise; inoise++) {
 | 
			
		||||
        for (int dk = 0; dk < LI_reduced; dk++) {
 | 
			
		||||
            for (int dt = 0; dt < Nt_inv; dt++) {
 | 
			
		||||
                for (int ds = 0; ds < SI; ds++) {
 | 
			
		||||
                    for (int is = 0; is < Ns; is++) {
 | 
			
		||||
                        result_nospin = peekSpin(solve[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))],is);
 | 
			
		||||
                        for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++) {
 | 
			
		||||
                            ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Tdir);
 | 
			
		||||
                            for (int ivec = 0; ivec < nvec_reduced; ivec++) {
 | 
			
		||||
    for (int inoise = 0; inoise < nnoise; inoise++)
 | 
			
		||||
    {
 | 
			
		||||
        for (int dk = 0; dk < LI_reduced; dk++)
 | 
			
		||||
       	{
 | 
			
		||||
            for (int dt = 0; dt < Nt_inv; dt++)
 | 
			
		||||
	    {
 | 
			
		||||
                for (int ds = 0; ds < SI; ds++)
 | 
			
		||||
	       	{
 | 
			
		||||
                    for (int is = 0; is < Ns; is++)
 | 
			
		||||
		    {
 | 
			
		||||
                        result4d_nospin = peekSpin(solve[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))],is);
 | 
			
		||||
                        for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++)
 | 
			
		||||
		       	{
 | 
			
		||||
                            ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir);
 | 
			
		||||
                            for (int ivec = 0; ivec < nvec_reduced; ivec++)
 | 
			
		||||
			    {
 | 
			
		||||
                                ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
 | 
			
		||||
                                pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result_3d)),is);
 | 
			
		||||
                                pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result3d_nospin)),is);
 | 
			
		||||
                                LOG(Message) <<  "perambulator(t, ivec, dk, inoise,dt,ds)(is) = (" << t << "," << ivec << "," << dk << "," << inoise << "," << dt << "," << ds << ")(" << is << ") = " <<  perambulator.tensor(t, ivec, dk, inoise,dt,ds)()(is)() << std::endl;
 | 
			
		||||
                            }
 | 
			
		||||
                        }
 | 
			
		||||
@@ -179,8 +199,6 @@ void TPerambFromSolve<FImpl>::execute(void)
 | 
			
		||||
    if(grid4d->IsBoss())
 | 
			
		||||
    {
 | 
			
		||||
        std::string sPerambName{par().PerambFileName};
 | 
			
		||||
        if (sPerambName.empty())
 | 
			
		||||
            sPerambName = getName();
 | 
			
		||||
        sPerambName.append( "." );
 | 
			
		||||
        sPerambName.append( std::to_string(vm().getTrajectory()));
 | 
			
		||||
        perambulator.write(sPerambName.c_str());
 | 
			
		||||
 
 | 
			
		||||
@@ -49,8 +49,7 @@ public:
 | 
			
		||||
                                    std::string, PerambFileName,
 | 
			
		||||
                                    std::string, UnsmearedSinkFileName,
 | 
			
		||||
                                    std::string, UnsmearedSinkMultiFile,
 | 
			
		||||
                                    int, nvec,
 | 
			
		||||
                                    DistilParameters, Distil);
 | 
			
		||||
                                    MDistil::DistilParameters, DistilPar);
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
@@ -106,8 +105,6 @@ std::vector<std::string> TPerambulator<FImpl>::getInput(void)
 | 
			
		||||
{
 | 
			
		||||
    sLapEvecName = par().lapevec;
 | 
			
		||||
    sNoiseName = par().noise;
 | 
			
		||||
    if( sNoiseName.length() == 0 )
 | 
			
		||||
        sNoiseName = getName() + "_noise";
 | 
			
		||||
    return {sLapEvecName, par().solver, sNoiseName };
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
@@ -124,23 +121,32 @@ void TPerambulator<FImpl>::setup(void)
 | 
			
		||||
    Cleanup();
 | 
			
		||||
    grid4d = env().getGrid();
 | 
			
		||||
    grid3d = MakeLowerDimGrid(grid4d);
 | 
			
		||||
    DISTIL_PARAMETERS_DEFINE( true );
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int tsrc{par().DistilPar.tsrc}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
    const std::string UnsmearedSinkFileName{ par().UnsmearedSinkFileName };
 | 
			
		||||
    if (!UnsmearedSinkFileName.empty())
 | 
			
		||||
        bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, true ) != 0 );
 | 
			
		||||
        //bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, true ) != 0 );
 | 
			
		||||
        bool bMulti = 0;
 | 
			
		||||
    
 | 
			
		||||
    envCreate(PerambTensor, getName(), 1, Nt,nvec,LI,nnoise,Nt_inv,SI);
 | 
			
		||||
    envCreate(std::vector<FermionField>, getName() + "_unsmeared_sink", 1,
 | 
			
		||||
              nnoise*LI*Ns*Nt_inv, envGetGrid(FermionField));
 | 
			
		||||
    
 | 
			
		||||
    envTmpLat(LatticeSpinColourVector, "dist_source");
 | 
			
		||||
    envTmpLat(LatticeSpinColourVector, "tmp2");
 | 
			
		||||
    envTmpLat(LatticeSpinColourVector, "result");
 | 
			
		||||
    envTmpLat(LatticeColourVector, "result_nospin");
 | 
			
		||||
    envTmp(LatticeSpinColourVector, "tmp3d",1,LatticeSpinColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "tmp3d_nospin",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "result_3d",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector, "evec3d",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmpLat(LatticeSpinColourVector,   "dist_source");
 | 
			
		||||
    envTmpLat(LatticeSpinColourVector,   "source4d");
 | 
			
		||||
    envTmp(LatticeSpinColourVector,      "source3d",1,LatticeSpinColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector,          "source3d_nospin",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmpLat(LatticeSpinColourVector,   "result4d");
 | 
			
		||||
    envTmpLat(LatticeColourVector,       "result4d_nospin");
 | 
			
		||||
    envTmp(LatticeColourVector,          "result3d_nospin",1,LatticeColourVector(grid3d));
 | 
			
		||||
    envTmp(LatticeColourVector,          "evec3d",1,LatticeColourVector(grid3d));
 | 
			
		||||
    
 | 
			
		||||
    Ls_ = env().getObjectLs(par().solver);
 | 
			
		||||
    envTmpLat(FermionField, "v4dtmp");
 | 
			
		||||
@@ -164,7 +170,16 @@ void TPerambulator<FImpl>::Cleanup(void)
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
void TPerambulator<FImpl>::execute(void)
 | 
			
		||||
{
 | 
			
		||||
    DISTIL_PARAMETERS_DEFINE( false );
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int tsrc{par().DistilPar.tsrc}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
 | 
			
		||||
    auto &solver=envGet(Solver, par().solver);
 | 
			
		||||
    auto &mat = solver.getFMat();
 | 
			
		||||
    envGetTmp(FermionField, v4dtmp);
 | 
			
		||||
@@ -175,12 +190,12 @@ void TPerambulator<FImpl>::execute(void)
 | 
			
		||||
    auto &epack = envGet(LapEvecs, sLapEvecName);
 | 
			
		||||
    auto &unsmeared_sink = envGet(std::vector<FermionField>, getName() + "_unsmeared_sink");
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, dist_source);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, tmp2);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, result);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result_nospin);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, tmp3d);
 | 
			
		||||
    envGetTmp(LatticeColourVector, tmp3d_nospin);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result_3d);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, source4d);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, source3d);
 | 
			
		||||
    envGetTmp(LatticeColourVector, source3d_nospin);
 | 
			
		||||
    envGetTmp(LatticeSpinColourVector, result4d);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result4d_nospin);
 | 
			
		||||
    envGetTmp(LatticeColourVector, result3d_nospin);
 | 
			
		||||
    envGetTmp(LatticeColourVector, evec3d);
 | 
			
		||||
    const int Ntlocal{grid4d->LocalDimensions()[3]};
 | 
			
		||||
    const int Ntfirst{grid4d->LocalStarts()[3]};
 | 
			
		||||
@@ -196,7 +211,7 @@ void TPerambulator<FImpl>::execute(void)
 | 
			
		||||
                {
 | 
			
		||||
                    LOG(Message) <<  "LapH source vector from noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl;
 | 
			
		||||
                    dist_source = 0;
 | 
			
		||||
                    tmp3d_nospin = 0;
 | 
			
		||||
                    source3d_nospin = 0;
 | 
			
		||||
                    evec3d = 0;
 | 
			
		||||
                    for (int it = dt; it < Nt; it += TI)
 | 
			
		||||
                    {
 | 
			
		||||
@@ -208,39 +223,39 @@ void TPerambulator<FImpl>::execute(void)
 | 
			
		||||
                                for (int is = ds; is < Ns; is += SI)
 | 
			
		||||
                                {
 | 
			
		||||
                                    ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv-Ntfirst,Tdir);
 | 
			
		||||
                                    tmp3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
 | 
			
		||||
                                    tmp3d=0;
 | 
			
		||||
                                    pokeSpin(tmp3d,tmp3d_nospin,is);
 | 
			
		||||
                                    tmp2=0;
 | 
			
		||||
                                    InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Tdir);
 | 
			
		||||
                                    dist_source += tmp2;
 | 
			
		||||
                                    source3d_nospin = evec3d * noise.tensor(inoise, t_inv, ik, is);
 | 
			
		||||
                                    source3d=0;
 | 
			
		||||
                                    pokeSpin(source3d,source3d_nospin,is);
 | 
			
		||||
                                    source4d=0;
 | 
			
		||||
                                    InsertSliceLocal(source3d,source4d,0,t_inv-Ntfirst,Tdir);
 | 
			
		||||
                                    dist_source += source4d;
 | 
			
		||||
                                }
 | 
			
		||||
                            }
 | 
			
		||||
                        }
 | 
			
		||||
                    }
 | 
			
		||||
                    result=0;
 | 
			
		||||
                    result4d=0;
 | 
			
		||||
                    v4dtmp = dist_source;
 | 
			
		||||
                    if (Ls_ == 1)
 | 
			
		||||
                        solver(result, v4dtmp);
 | 
			
		||||
                        solver(result4d, v4dtmp);
 | 
			
		||||
                    else
 | 
			
		||||
                    {
 | 
			
		||||
                        mat.ImportPhysicalFermionSource(v4dtmp, v5dtmp);
 | 
			
		||||
                        solver(v5dtmp_sol, v5dtmp);
 | 
			
		||||
                        mat.ExportPhysicalFermionSolution(v5dtmp_sol, v4dtmp);
 | 
			
		||||
                        result = v4dtmp;
 | 
			
		||||
                        result4d = v4dtmp;
 | 
			
		||||
                    }
 | 
			
		||||
                    if (!UnsmearedSinkFileName.empty())
 | 
			
		||||
                        unsmeared_sink[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))] = result;
 | 
			
		||||
                        unsmeared_sink[inoise+nnoise*(dk+LI*(dt+Nt_inv*ds))] = result4d;
 | 
			
		||||
                    for (int is = 0; is < Ns; is++)
 | 
			
		||||
                    {
 | 
			
		||||
                        result_nospin = peekSpin(result,is);
 | 
			
		||||
                        result4d_nospin = peekSpin(result4d,is);
 | 
			
		||||
                        for (int t = Ntfirst; t < Ntfirst + Ntlocal; t++)
 | 
			
		||||
                        {
 | 
			
		||||
                            ExtractSliceLocal(result_3d,result_nospin,0,t-Ntfirst,Tdir);
 | 
			
		||||
                            for (int ivec = 0; ivec < nvec; ivec++)
 | 
			
		||||
                            ExtractSliceLocal(result3d_nospin,result4d_nospin,0,t-Ntfirst,Tdir); 
 | 
			
		||||
			    for (int ivec = 0; ivec < nvec; ivec++)
 | 
			
		||||
                            {
 | 
			
		||||
                                ExtractSliceLocal(evec3d,epack.evec[ivec],0,t-Ntfirst,Tdir);
 | 
			
		||||
                                pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result_3d)),is);
 | 
			
		||||
                                pokeSpin(perambulator.tensor(t, ivec, dk, inoise,dt,ds),static_cast<Complex>(innerProduct(evec3d, result3d_nospin)),is);
 | 
			
		||||
                            }
 | 
			
		||||
                        }
 | 
			
		||||
                    }
 | 
			
		||||
@@ -284,7 +299,8 @@ void TPerambulator<FImpl>::execute(void)
 | 
			
		||||
    //Save the unsmeared sinks if filename specified
 | 
			
		||||
    if (!UnsmearedSinkFileName.empty())
 | 
			
		||||
    {
 | 
			
		||||
        bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, false ) != 0 );
 | 
			
		||||
       // bool bMulti = ( Hadrons::MDistil::DistilParameters::ParameterDefault( par().UnsmearedSinkMultiFile, 1, false ) != 0 );
 | 
			
		||||
        bool bMulti =0;
 | 
			
		||||
        LOG(Message) << "Writing unsmeared sink to " << UnsmearedSinkFileName << std::endl;
 | 
			
		||||
        A2AVectorsIo::write(UnsmearedSinkFileName, unsmeared_sink, bMulti, vm().getTrajectory());
 | 
			
		||||
    }
 | 
			
		||||
 
 | 
			
		||||
@@ -44,8 +44,7 @@ class LoadDistilNoisePar: Serializable
 | 
			
		||||
public:
 | 
			
		||||
    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadDistilNoisePar,
 | 
			
		||||
                                        std::string, NoiseFileName,
 | 
			
		||||
	                                int, nvec,
 | 
			
		||||
	                                MDistil::DistilParameters, Distil);
 | 
			
		||||
                                        MDistil::DistilParameters, DistilPar);
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
@@ -97,8 +96,10 @@ std::vector<std::string> TLoadDistilNoise<FImpl>::getOutput(void)
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
void TLoadDistilNoise<FImpl>::setup(void)
 | 
			
		||||
{
 | 
			
		||||
  DISTIL_PARAMETERS_DEFINE( true );
 | 
			
		||||
  envCreate(MDistil::NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns);
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    envCreate(MDistil::NoiseTensor, getName(), 1, nnoise, Nt, nvec, Ns);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
// execution ///////////////////////////////////////////////////////////////////
 | 
			
		||||
 
 | 
			
		||||
@@ -43,9 +43,8 @@ class LoadPerambulatorPar: Serializable
 | 
			
		||||
{
 | 
			
		||||
public:
 | 
			
		||||
    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadPerambulatorPar,
 | 
			
		||||
                                        std::string, PerambFileName, //stem!!!
 | 
			
		||||
	                                int, nvec,
 | 
			
		||||
	                                MDistil::DistilParameters, Distil);
 | 
			
		||||
                                        std::string, PerambFileName,
 | 
			
		||||
                                        MDistil::DistilParameters, DistilPar);
 | 
			
		||||
};
 | 
			
		||||
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
@@ -97,8 +96,16 @@ std::vector<std::string> TLoadPerambulator<FImpl>::getOutput(void)
 | 
			
		||||
template <typename FImpl>
 | 
			
		||||
void TLoadPerambulator<FImpl>::setup(void)
 | 
			
		||||
{
 | 
			
		||||
  DISTIL_PARAMETERS_DEFINE( true );
 | 
			
		||||
  //std::array<std::string,6> sIndexNames{"Nt", "nvec", "LI", "nnoise", "Nt_inv", "SI"};
 | 
			
		||||
    const int Nt{env().getDim(Tdir)}; 
 | 
			
		||||
    const int nvec{par().DistilPar.nvec}; 
 | 
			
		||||
    const int nnoise{par().DistilPar.nnoise}; 
 | 
			
		||||
    const int tsrc{par().DistilPar.tsrc}; 
 | 
			
		||||
    const int TI{par().DistilPar.TI}; 
 | 
			
		||||
    const int LI{par().DistilPar.LI}; 
 | 
			
		||||
    const int SI{par().DistilPar.SI}; 
 | 
			
		||||
    const bool full_tdil{ TI == Nt }; 
 | 
			
		||||
    const bool exact_distillation{ full_tdil && LI == nvec }; 
 | 
			
		||||
    const int Nt_inv{ full_tdil ? 1 : TI };
 | 
			
		||||
  envCreate(MDistil::PerambTensor, getName(), 1, Nt,nvec,LI,nnoise,Nt_inv,SI);
 | 
			
		||||
}
 | 
			
		||||
 | 
			
		||||
 
 | 
			
		||||
@@ -65,7 +65,9 @@ modules_cc =\
 | 
			
		||||
  Modules/MScalarSUN/TrMag.cc \
 | 
			
		||||
  Modules/MDistil/PerambFromSolve.cc \
 | 
			
		||||
  Modules/MDistil/DistilVectors.cc \
 | 
			
		||||
  Modules/MDistil/DistilPar.cc \
 | 
			
		||||
  Modules/MDistil/LapEvec.cc \
 | 
			
		||||
  Modules/MDistil/DistilParameters.cc \
 | 
			
		||||
  Modules/MDistil/Perambulator.cc \
 | 
			
		||||
  Modules/MDistil/Noises.cc \
 | 
			
		||||
  Modules/MAction/MobiusDWF.cc \
 | 
			
		||||
@@ -151,6 +153,7 @@ modules_hpp =\
 | 
			
		||||
  Modules/MDistil/Noises.hpp \
 | 
			
		||||
  Modules/MDistil/LapEvec.hpp \
 | 
			
		||||
  Modules/MDistil/Perambulator.hpp \
 | 
			
		||||
  Modules/MDistil/DistilPar.hpp \
 | 
			
		||||
  Modules/MDistil/DistilCommon.hpp \
 | 
			
		||||
  Modules/MAction/ZMobiusDWF.hpp \
 | 
			
		||||
  Modules/MAction/ScaledDWF.hpp \
 | 
			
		||||
 
 | 
			
		||||
		Reference in New Issue
	
	Block a user