diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp index 40df28e0..272380d1 100644 --- a/Hadrons/Modules/MDistil/DistilVectors.hpp +++ b/Hadrons/Modules/MDistil/DistilVectors.hpp @@ -8,59 +8,60 @@ BEGIN_HADRONS_NAMESPACE /****************************************************************************** - * * Class to generate Distillation $\varrho$ and $\varphi$ vectors * - * ******************************************************************************/ + * Class to generate Distillation $\varrho$ and $\varphi$ vectors * + ******************************************************************************/ /* - * template - * class DistilVectorsFromPerambulator - * { - * public: - * FERM_TYPE_ALIASES(FImpl,); - * SOLVER_TYPE_ALIASES(FImpl,); - * public: - * DistilVectorsFromPerambulator(FMat &action); - * virtual ~DistilVectorsFromPerambulator(void) = default; - * void makeRho(FermionField &rhoOut, - * const FermionField &evec, const Complex &noises); - * void makePhi(FermionField &phiOut, - * const FermionField &evec, const SpinVector &Perambulators); - * private: - * GridBase *fGrid_, *frbGrid_, *gGrid_; - * bool is5d_; - * FermionField src_o_, sol_e_, sol_o_, tmp_, tmp5_; - * }; - * - * */ + template + class DistilVectorsFromPerambulator + { + public: + FERM_TYPE_ALIASES(FImpl,); + SOLVER_TYPE_ALIASES(FImpl,); + public: + DistilVectorsFromPerambulator(FMat &action); + virtual ~DistilVectorsFromPerambulator(void) = default; + void makeRho(FermionField &rhoOut, + const FermionField &evec, const Complex &noises); + void makePhi(FermionField &phiOut, + const FermionField &evec, const SpinVector &Perambulators); + private: + GridBase *fGrid_, *frbGrid_, *gGrid_; + bool is5d_; + FermionField src_o_, sol_e_, sol_o_, tmp_, tmp5_; + }; + +*/ /****************************************************************************** - * * * A2AVectorsSchurDiagTwo template implementation * - * * ******************************************************************************/ + * A2AVectorsSchurDiagTwo template implementation * + ******************************************************************************/ /* - * - * template - * A2AVectorsSchurDiagTwo::A2AVectorsSchurDiagTwo(FMat &action) - * : action_(action) - * , fGrid_(action_.FermionGrid()) - * , frbGrid_(action_.FermionRedBlackGrid()) - * , gGrid_(action_.GaugeGrid()) - * , src_o_(frbGrid_) - * , sol_e_(frbGrid_) - * , sol_o_(frbGrid_) - * , tmp_(frbGrid_) - * , tmp5_(fGrid_) - * , op_(action_) - * {} - * - * template - * void A2AVectorsSchurDiagTwo::makeRho(FermionField &rhoOut, const FermionField &evec, const Complex &noises) - * { - * - * LatticeSpinColourVector tmp2(grid4d); - * LatticeSpinColourVector tmp3d(grid3d); - * LatticeColourVector tmp3d_nospin(grid3d); - * LatticeColourVector tmp_nospin(grid4d); - * + +template +A2AVectorsSchurDiagTwo::A2AVectorsSchurDiagTwo(FMat &action) +: action_(action) +, fGrid_(action_.FermionGrid()) +, frbGrid_(action_.FermionRedBlackGrid()) +, gGrid_(action_.GaugeGrid()) +, src_o_(frbGrid_) +, sol_e_(frbGrid_) +, sol_o_(frbGrid_) +, tmp_(frbGrid_) +, tmp5_(fGrid_) +, op_(action_) +{} + +template +void A2AVectorsSchurDiagTwo::makeRho(FermionField &rhoOut, const FermionField &evec, const Complex &noises) +{ + + LatticeSpinColourVector tmp2(grid4d); + LatticeSpinColourVector tmp3d(grid3d); + LatticeColourVector tmp3d_nospin(grid3d); + LatticeColourVector tmp_nospin(grid4d); + LatticeColourVector evec3d(grid3d); + for (int inoise = 0; inoise < nnoise; inoise++) { for (int dk = 0; dk < LI; dk++) { @@ -74,7 +75,8 @@ BEGIN_HADRONS_NAMESPACE if( it >= Ntfirst && it < Ntfirst + Ntlocal ) { for (int ik = dk; ik < nvec; ik += LI){ for (int is = ds; is < Ns; is += Ns){ //at the moment, full spin dilution is enforced - tmp3d_nospin = eig[it].evec[ik] * noises[inoise][it][ik]()(is)(); //noises do not have to be a spin vector + ExtractSliceLocal(evec3d,eig4d.evec[ik],0,it,3); + tmp3d_nospin = evec3d * noises[inoise][it][ik]()(is)(); //noises do not have to be a spin vector tmp3d=zero; pokeSpin(tmp3d,tmp3d_nospin,is); tmp2=zero;