diff --git a/extras/Hadrons/AllToAllVectors.hpp b/extras/Hadrons/AllToAllVectors.hpp new file mode 100644 index 00000000..7d9aa62e --- /dev/null +++ b/extras/Hadrons/AllToAllVectors.hpp @@ -0,0 +1,217 @@ +#ifndef A2A_Vectors_hpp_ +#define A2A_Vectors_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +//////////////////////////////// +// A2A Modes +//////////////////////////////// + +template +class A2AModesSchurDiagTwo +{ + private: + const std::vector *evec; + const std::vector *eval; + Matrix &action; + Solver &solver; + const int Nl, Nh; + const bool return_5d; + std::vector w_high_5d, v_high_5d, w_high_4d, v_high_4d; + + public: + A2AModesSchurDiagTwo(const std::vector *_evec, const std::vector *_eval, + Matrix &_action, + Solver &_solver, + const int _Nl, const int _Nh, + const bool _return_5d) + : evec(_evec), eval(_eval), + action(_action), + solver(_solver), + Nl(_Nl), Nh(_Nh), + return_5d(_return_5d) + { + init_resize(1, Nh); + if (return_5d) init_resize(Nh, Nh); + }; + + void init_resize(const size_t size_5d, const size_t size_4d) + { + GridBase *grid_5d = action.Grid(); + GridBase *grid_4d = action.GaugeGrid(); + + w_high_5d.resize(size_5d, grid_5d); + v_high_5d.resize(size_5d, grid_5d); + + w_high_4d.resize(size_4d, grid_4d); + v_high_4d.resize(size_4d, grid_4d); + } + + void high_modes(Field &source_5d, Field &w_source_5d, Field &source_4d, int i) + { + int i5d; + LOG(Message) << "A2A high modes for i = " << i << std::endl; + i5d = 0; + if (return_5d) i5d = i; + this->high_mode_v(action, solver, source_5d, v_high_5d[i5d], v_high_4d[i]); + this->high_mode_w(w_source_5d, source_4d, w_high_5d[i5d], w_high_4d[i]); + } + + void return_v(int i, Field &vout_5d, Field &vout_4d) + { + if (i < Nl) + { + this->low_mode_v(action, evec->at(i), eval->at(i), vout_5d, vout_4d); + } + else + { + vout_4d = v_high_4d[i - Nl]; + if (!(return_5d)) i = Nl; + vout_5d = v_high_5d[i - Nl]; + } + } + void return_w(int i, Field &wout_5d, Field &wout_4d) + { + if (i < Nl) + { + this->low_mode_w(action, evec->at(i), eval->at(i), wout_5d, wout_4d); + } + else + { + wout_4d = w_high_4d[i - Nl]; + if (!(return_5d)) i = Nl; + wout_5d = w_high_5d[i - Nl]; + } + } + + void low_mode_v(Matrix &action, const Field &evec, const RealD &eval, Field &vout_5d, Field &vout_4d) + { + GridBase *grid = action.RedBlackGrid(); + Field src_o(grid); + Field sol_e(grid); + Field sol_o(grid); + Field tmp(grid); + + src_o = evec; + src_o.checkerboard = Odd; + pickCheckerboard(Even, sol_e, vout_5d); + pickCheckerboard(Odd, sol_o, vout_5d); + + ///////////////////////////////////////////////////// + // v_ie = -(1/eval_i) * MeeInv Meo MooInv evec_i + ///////////////////////////////////////////////////// + action.MooeeInv(src_o, tmp); + assert(tmp.checkerboard == Odd); + action.Meooe(tmp, sol_e); + assert(sol_e.checkerboard == Even); + action.MooeeInv(sol_e, tmp); + assert(tmp.checkerboard == Even); + sol_e = (-1.0 / eval) * tmp; + assert(sol_e.checkerboard == Even); + + ///////////////////////////////////////////////////// + // v_io = (1/eval_i) * MooInv evec_i + ///////////////////////////////////////////////////// + action.MooeeInv(src_o, tmp); + assert(tmp.checkerboard == Odd); + sol_o = (1.0 / eval) * tmp; + assert(sol_o.checkerboard == Odd); + + setCheckerboard(vout_5d, sol_e); + assert(sol_e.checkerboard == Even); + setCheckerboard(vout_5d, sol_o); + assert(sol_o.checkerboard == Odd); + + action.ExportPhysicalFermionSolution(vout_5d, vout_4d); + } + + void low_mode_w(Matrix &action, const Field &evec, const RealD &eval, Field &wout_5d, Field &wout_4d) + { + GridBase *grid = action.RedBlackGrid(); + SchurDiagTwoOperator _HermOpEO(action); + + Field src_o(grid); + Field sol_e(grid); + Field sol_o(grid); + Field tmp(grid); + + GridBase *fgrid = action.Grid(); + Field tmp_wout(fgrid); + + src_o = evec; + src_o.checkerboard = Odd; + pickCheckerboard(Even, sol_e, tmp_wout); + pickCheckerboard(Odd, sol_o, tmp_wout); + + ///////////////////////////////////////////////////// + // w_ie = - MeeInvDag MoeDag Doo evec_i + ///////////////////////////////////////////////////// + _HermOpEO.Mpc(src_o, tmp); + assert(tmp.checkerboard == Odd); + action.MeooeDag(tmp, sol_e); + assert(sol_e.checkerboard == Even); + action.MooeeInvDag(sol_e, tmp); + assert(tmp.checkerboard == Even); + sol_e = (-1.0) * tmp; + + ///////////////////////////////////////////////////// + // w_io = Doo evec_i + ///////////////////////////////////////////////////// + _HermOpEO.Mpc(src_o, sol_o); + assert(sol_o.checkerboard == Odd); + + setCheckerboard(tmp_wout, sol_e); + assert(sol_e.checkerboard == Even); + setCheckerboard(tmp_wout, sol_o); + assert(sol_o.checkerboard == Odd); + + action.DminusDag(tmp_wout, wout_5d); + + action.ExportPhysicalFermionSource(wout_5d, wout_4d); + } + + void high_mode_v(Matrix &action, Solver &solver, const Field &source, Field &vout_5d, Field &vout_4d) + { + GridBase *fgrid = action.Grid(); + solver(vout_5d, source); // Note: solver is solver(out, in) + action.ExportPhysicalFermionSolution(vout_5d, vout_4d); + } + + void high_mode_w(const Field &w_source_5d, const Field &source_4d, Field &wout_5d, Field &wout_4d) + { + wout_5d = w_source_5d; + wout_4d = source_4d; + } +}; + +// TODO: A2A for coarse eigenvectors + +// template +// class A2ALMSchurDiagTwoCoarse : public A2AModesSchurDiagTwo +// { +// private: +// const std::vector &subspace; +// const std::vector &evec_coarse; +// const std::vector &eval_coarse; +// Matrix &action; + +// public: +// A2ALMSchurDiagTwoCoarse(const std::vector &_subspace, const std::vector &_evec_coarse, const std::vector &_eval_coarse, Matrix &_action) +// : subspace(_subspace), evec_coarse(_evec_coarse), eval_coarse(_eval_coarse), action(_action){}; + +// void operator()(int i, FineField &vout, FineField &wout) +// { +// FineField prom_evec(subspace[0]._grid); +// blockPromote(evec_coarse[i], prom_evec, subspace); +// this->low_mode_v(action, prom_evec, eval_coarse[i], vout); +// this->low_mode_w(action, prom_evec, eval_coarse[i], wout); +// } +// }; + +END_HADRONS_NAMESPACE + +#endif // A2A_Vectors_hpp_ \ No newline at end of file diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 6bf5d3ed..06cd804a 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -1,57 +1,58 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include #include +#include #include #include -#include -#include #include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include #include +#include +#include +#include +#include #include -#include #include +#include #include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include diff --git a/extras/Hadrons/Modules/MSolver/A2AVectors.cc b/extras/Hadrons/Modules/MSolver/A2AVectors.cc new file mode 100644 index 00000000..f72f405d --- /dev/null +++ b/extras/Hadrons/Modules/MSolver/A2AVectors.cc @@ -0,0 +1,8 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MSolver; + +template class Grid::Hadrons::MSolver::TA2AVectors; +template class Grid::Hadrons::MSolver::TA2AVectors; diff --git a/extras/Hadrons/Modules/MSolver/A2AVectors.hpp b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp new file mode 100644 index 00000000..d7e4cfe4 --- /dev/null +++ b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp @@ -0,0 +1,235 @@ +#ifndef Hadrons_MSolver_A2AVectors_hpp_ +#define Hadrons_MSolver_A2AVectors_hpp_ + +#include +#include +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * A2AVectors * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSolver) + +class A2AVectorsPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(A2AVectorsPar, + bool, return_5d, + int, Nl, + int, N, + std::vector, sources, + std::string, action, + std::string, eigenPack, + std::string, solver); +}; + +template +class TA2AVectors : public Module +{ + public: + FERM_TYPE_ALIASES(FImpl,); + SOLVER_TYPE_ALIASES(FImpl,); + + typedef FermionEigenPack EPack; + typedef CoarseFermionEigenPack CoarseEPack; + + typedef A2AModesSchurDiagTwo A2ABase; + + public: + // constructor + TA2AVectors(const std::string name); + // destructor + virtual ~TA2AVectors(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getReference(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); + + private: + unsigned int Ls_; + std::string className_; +}; + +MODULE_REGISTER_TMP(A2AVectors, ARG(TA2AVectors), MSolver); +MODULE_REGISTER_TMP(ZA2AVectors, ARG(TA2AVectors), MSolver); + +/****************************************************************************** + * TA2AVectors implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TA2AVectors::TA2AVectors(const std::string name) +: Module(name) +, className_ (name + "_class") +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TA2AVectors::getInput(void) +{ + int Nl = par().Nl; + std::string sub_string = ""; + if (Nl > 0) sub_string = "_subtract"; + std::vector in = {par().solver + sub_string}; + + int n = par().sources.size(); + + for (unsigned int t = 0; t < n; t += 1) + { + in.push_back(par().sources[t]); + } + + return in; +} + +template +std::vector TA2AVectors::getReference(void) +{ + std::vector ref = {par().action}; + + if (!par().eigenPack.empty()) + { + ref.push_back(par().eigenPack); + } + + return ref; +} + +template +std::vector TA2AVectors::getOutput(void) +{ + std::vector out = {getName(), className_}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TA2AVectors::setup(void) +{ + int N = par().N; + int Nl = par().Nl; + int Nh = N - Nl; + bool return_5d = par().return_5d; + int Ls; + + std::string sub_string = ""; + if (Nl > 0) sub_string = "_subtract"; + auto &solver = envGet(Solver, par().solver + sub_string); + Ls = env().getObjectLs(par().solver + sub_string); + + auto &action = envGet(FMat, par().action); + + envTmpLat(FermionField, "ferm_src", Ls); + envTmpLat(FermionField, "unphys_ferm", Ls); + envTmpLat(FermionField, "tmp"); + + std::vector *evec; + const std::vector *eval; + + if (Nl > 0) + { + // Low modes + auto &epack = envGet(EPack, par().eigenPack); + + LOG(Message) << "Creating a2a vectors " << getName() << + " using eigenpack '" << par().eigenPack << "' (" + << epack.evec.size() << " modes)" << + " and " << Nh << " high modes." << std::endl; + evec = &epack.evec; + eval = &epack.eval; + } + else + { + LOG(Message) << "Creating a2a vectors " << getName() << + " using " << Nh << " high modes only." << std::endl; + } + + envCreate(A2ABase, className_, Ls, + evec, eval, + action, + solver, + Nl, Nh, + return_5d); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TA2AVectors::execute(void) +{ + auto &action = envGet(FMat, par().action); + + int Nt = env().getDim(Tp); + int Nc = FImpl::Dimension; + int Ls_; + int Nl = par().Nl; + + std::string sub_string = ""; + if (Nl > 0) sub_string = "_subtract"; + Ls_ = env().getObjectLs(par().solver + sub_string); + + auto &a2areturn = envGet(A2ABase, className_); + + // High modes + auto sources = par().sources; + int Nsrc = par().sources.size(); + + envGetTmp(FermionField, ferm_src); + envGetTmp(FermionField, unphys_ferm); + envGetTmp(FermionField, tmp); + + int N_count = 0; + for (unsigned int s = 0; s < Ns; ++s) + for (unsigned int c = 0; c < Nc; ++c) + for (unsigned int T = 0; T < Nsrc; T++) + { + auto &prop_src = envGet(PropagatorField, sources[T]); + LOG(Message) << "A2A src for s = " << s << " , c = " << c << ", T = " << T << std::endl; + // source conversion for 4D sources + if (!env().isObject5d(sources[T])) + { + if (Ls_ == 1) + { + PropToFerm(ferm_src, prop_src, s, c); + tmp = ferm_src; + } + else + { + PropToFerm(tmp, prop_src, s, c); + action.ImportPhysicalFermionSource(tmp, ferm_src); + action.ImportUnphysicalFermion(tmp, unphys_ferm); + } + } + // source conversion for 5D sources + else + { + if (Ls_ != env().getObjectLs(sources[T])) + { + HADRONS_ERROR(Size, "Ls mismatch between quark action and source"); + } + else + { + PropToFerm(ferm_src, prop_src, s, c); + action.ExportPhysicalFermionSolution(ferm_src, tmp); + unphys_ferm = ferm_src; + } + } + LOG(Message) << "a2areturn.high_modes Ncount = " << N_count << std::endl; + a2areturn.high_modes(ferm_src, unphys_ferm, tmp, N_count); + N_count++; + } +} +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MSolver_A2AVectors_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 58b082d0..ca454086 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,115 +1,117 @@ modules_cc =\ - Modules/MContraction/WeakHamiltonianEye.cc \ - Modules/MContraction/Baryon.cc \ - Modules/MContraction/Meson.cc \ - Modules/MContraction/WeakNeutral4ptDisc.cc \ - Modules/MContraction/WeakHamiltonianNonEye.cc \ - Modules/MContraction/WardIdentity.cc \ - Modules/MContraction/DiscLoop.cc \ - Modules/MContraction/Gamma3pt.cc \ - Modules/MFermion/FreeProp.cc \ - Modules/MFermion/GaugeProp.cc \ - Modules/MSource/Point.cc \ - Modules/MSource/Wall.cc \ Modules/MSource/SeqConserved.cc \ Modules/MSource/SeqGamma.cc \ + Modules/MSource/Point.cc \ Modules/MSource/Z2.cc \ + Modules/MSource/Wall.cc \ Modules/MSink/Point.cc \ Modules/MSink/Smear.cc \ - Modules/MSolver/RBPrecCG.cc \ - Modules/MSolver/LocalCoherenceLanczos.cc \ - Modules/MGauge/StoutSmearing.cc \ - Modules/MGauge/Unit.cc \ - Modules/MGauge/UnitEm.cc \ - Modules/MGauge/StochEm.cc \ - Modules/MGauge/Random.cc \ - Modules/MGauge/FundtoHirep.cc \ - Modules/MUtilities/TestSeqGamma.cc \ - Modules/MUtilities/TestSeqConserved.cc \ - Modules/MLoop/NoiseLoop.cc \ - Modules/MScalar/FreeProp.cc \ - Modules/MScalar/VPCounterTerms.cc \ - Modules/MScalar/ChargedProp.cc \ - Modules/MScalar/ScalarVP.cc \ - Modules/MAction/Wilson.cc \ + Modules/MIO/LoadNersc.cc \ + Modules/MIO/LoadEigenPack.cc \ + Modules/MIO/LoadCoarseEigenPack.cc \ + Modules/MIO/LoadBinary.cc \ + Modules/MScalarSUN/TwoPointNPR.cc \ + Modules/MScalarSUN/TrPhi.cc \ + Modules/MScalarSUN/StochFreeField.cc \ + Modules/MScalarSUN/TrMag.cc \ + Modules/MScalarSUN/Div.cc \ + Modules/MScalarSUN/ShiftProbe.cc \ + Modules/MScalarSUN/Grad.cc \ + Modules/MScalarSUN/TwoPoint.cc \ + Modules/MScalarSUN/TimeMomProbe.cc \ + Modules/MScalarSUN/EMT.cc \ + Modules/MScalarSUN/TransProj.cc \ + Modules/MScalarSUN/TrKinetic.cc \ + Modules/MAction/DWF.cc \ Modules/MAction/MobiusDWF.cc \ Modules/MAction/ZMobiusDWF.cc \ - Modules/MAction/WilsonClover.cc \ - Modules/MAction/DWF.cc \ + Modules/MAction/Wilson.cc \ Modules/MAction/ScaledDWF.cc \ - Modules/MScalarSUN/TrPhi.cc \ - Modules/MScalarSUN/Grad.cc \ - Modules/MScalarSUN/TimeMomProbe.cc \ - Modules/MScalarSUN/TrMag.cc \ - Modules/MScalarSUN/TrKinetic.cc \ - Modules/MScalarSUN/EMT.cc \ - Modules/MScalarSUN/ShiftProbe.cc \ - Modules/MScalarSUN/TransProj.cc \ - Modules/MScalarSUN/StochFreeField.cc \ - Modules/MScalarSUN/TwoPoint.cc \ - Modules/MScalarSUN/TwoPointNPR.cc \ - Modules/MScalarSUN/Div.cc \ - Modules/MIO/LoadEigenPack.cc \ - Modules/MIO/LoadBinary.cc \ - Modules/MIO/LoadNersc.cc \ - Modules/MIO/LoadCoarseEigenPack.cc + Modules/MAction/WilsonClover.cc \ + Modules/MContraction/WeakHamiltonianNonEye.cc \ + Modules/MContraction/WardIdentity.cc \ + Modules/MContraction/WeakHamiltonianEye.cc \ + Modules/MContraction/DiscLoop.cc \ + Modules/MContraction/Baryon.cc \ + Modules/MContraction/Gamma3pt.cc \ + Modules/MContraction/WeakNeutral4ptDisc.cc \ + Modules/MContraction/Meson.cc \ + Modules/MScalar/VPCounterTerms.cc \ + Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/FreeProp.cc \ + Modules/MScalar/ScalarVP.cc \ + Modules/MUtilities/TestSeqGamma.cc \ + Modules/MUtilities/TestSeqConserved.cc \ + Modules/MFermion/FreeProp.cc \ + Modules/MFermion/GaugeProp.cc \ + Modules/MSolver/RBPrecCG.cc \ + Modules/MSolver/LocalCoherenceLanczos.cc \ + Modules/MSolver/A2AVectors.cc \ + Modules/MLoop/NoiseLoop.cc \ + Modules/MGauge/Unit.cc \ + Modules/MGauge/Random.cc \ + Modules/MGauge/UnitEm.cc \ + Modules/MGauge/StochEm.cc \ + Modules/MGauge/FundtoHirep.cc \ + Modules/MGauge/StoutSmearing.cc modules_hpp =\ - Modules/MContraction/Baryon.hpp \ - Modules/MContraction/Meson.hpp \ - Modules/MContraction/WeakHamiltonian.hpp \ - Modules/MContraction/WeakHamiltonianNonEye.hpp \ - Modules/MContraction/DiscLoop.hpp \ - Modules/MContraction/WeakNeutral4ptDisc.hpp \ - Modules/MContraction/Gamma3pt.hpp \ - Modules/MContraction/WardIdentity.hpp \ - Modules/MContraction/WeakHamiltonianEye.hpp \ - Modules/MFermion/FreeProp.hpp \ - Modules/MFermion/GaugeProp.hpp \ + Modules/MSource/SeqConserved.hpp \ Modules/MSource/SeqGamma.hpp \ + Modules/MSource/Z2.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/Wall.hpp \ - Modules/MSource/Z2.hpp \ - Modules/MSource/SeqConserved.hpp \ Modules/MSink/Smear.hpp \ Modules/MSink/Point.hpp \ - Modules/MSolver/LocalCoherenceLanczos.hpp \ - Modules/MSolver/RBPrecCG.hpp \ - Modules/MGauge/UnitEm.hpp \ - Modules/MGauge/StoutSmearing.hpp \ - Modules/MGauge/Unit.hpp \ - Modules/MGauge/Random.hpp \ - Modules/MGauge/FundtoHirep.hpp \ - Modules/MGauge/StochEm.hpp \ - Modules/MUtilities/TestSeqGamma.hpp \ - Modules/MUtilities/TestSeqConserved.hpp \ - Modules/MLoop/NoiseLoop.hpp \ - Modules/MScalar/FreeProp.hpp \ - Modules/MScalar/VPCounterTerms.hpp \ - Modules/MScalar/ScalarVP.hpp \ - Modules/MScalar/Scalar.hpp \ - Modules/MScalar/ChargedProp.hpp \ - Modules/MAction/DWF.hpp \ - Modules/MAction/MobiusDWF.hpp \ - Modules/MAction/Wilson.hpp \ - Modules/MAction/WilsonClover.hpp \ - Modules/MAction/ZMobiusDWF.hpp \ - Modules/MAction/ScaledDWF.hpp \ - Modules/MScalarSUN/StochFreeField.hpp \ + Modules/MIO/LoadBinary.hpp \ + Modules/MIO/LoadEigenPack.hpp \ + Modules/MIO/LoadCoarseEigenPack.hpp \ + Modules/MIO/LoadNersc.hpp \ + Modules/MScalarSUN/Utils.hpp \ + Modules/MScalarSUN/Grad.hpp \ + Modules/MScalarSUN/TrPhi.hpp \ Modules/MScalarSUN/TwoPointNPR.hpp \ + Modules/MScalarSUN/TwoPoint.hpp \ + Modules/MScalarSUN/TransProj.hpp \ + Modules/MScalarSUN/TrKinetic.hpp \ + Modules/MScalarSUN/StochFreeField.hpp \ Modules/MScalarSUN/ShiftProbe.hpp \ - Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TimeMomProbe.hpp \ + Modules/MScalarSUN/Div.hpp \ Modules/MScalarSUN/TrMag.hpp \ Modules/MScalarSUN/EMT.hpp \ - Modules/MScalarSUN/TwoPoint.hpp \ - Modules/MScalarSUN/TrPhi.hpp \ - Modules/MScalarSUN/Utils.hpp \ - Modules/MScalarSUN/TransProj.hpp \ - Modules/MScalarSUN/Grad.hpp \ - Modules/MScalarSUN/TrKinetic.hpp \ - Modules/MIO/LoadEigenPack.hpp \ - Modules/MIO/LoadNersc.hpp \ - Modules/MIO/LoadCoarseEigenPack.hpp \ - Modules/MIO/LoadBinary.hpp + Modules/MAction/ZMobiusDWF.hpp \ + Modules/MAction/ScaledDWF.hpp \ + Modules/MAction/Wilson.hpp \ + Modules/MAction/WilsonClover.hpp \ + Modules/MAction/MobiusDWF.hpp \ + Modules/MAction/DWF.hpp \ + Modules/MContraction/WeakHamiltonian.hpp \ + Modules/MContraction/DiscLoop.hpp \ + Modules/MContraction/Meson.hpp \ + Modules/MContraction/WardIdentity.hpp \ + Modules/MContraction/WeakHamiltonianEye.hpp \ + Modules/MContraction/Gamma3pt.hpp \ + Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MContraction/Baryon.hpp \ + Modules/MContraction/WeakNeutral4ptDisc.hpp \ + Modules/MScalar/ScalarVP.hpp \ + Modules/MScalar/Scalar.hpp \ + Modules/MScalar/FreeProp.hpp \ + Modules/MScalar/ChargedProp.hpp \ + Modules/MScalar/VPCounterTerms.hpp \ + Modules/MUtilities/TestSeqConserved.hpp \ + Modules/MUtilities/TestSeqGamma.hpp \ + Modules/MFermion/FreeProp.hpp \ + Modules/MFermion/GaugeProp.hpp \ + Modules/MSolver/A2AVectors.hpp \ + Modules/MSolver/RBPrecCG.hpp \ + Modules/MSolver/LocalCoherenceLanczos.hpp \ + Modules/MLoop/NoiseLoop.hpp \ + Modules/MGauge/StoutSmearing.hpp \ + Modules/MGauge/StochEm.hpp \ + Modules/MGauge/FundtoHirep.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MGauge/Random.hpp \ + Modules/MGauge/UnitEm.hpp diff --git a/lib/algorithms/iterative/Deflation.h b/lib/algorithms/iterative/Deflation.h index 316afe90..7c3e8e57 100644 --- a/lib/algorithms/iterative/Deflation.h +++ b/lib/algorithms/iterative/Deflation.h @@ -63,7 +63,7 @@ public: DeflatedGuesser(const std::vector & _evec,const std::vector & _eval) : evec(_evec), eval(_eval) {}; - virtual void operator()(const Field &src,Field &guess) { + virtual void operator()(const Field &src,Field &guess) { guess = zero; assert(evec.size()==eval.size()); auto N = evec.size(); @@ -71,6 +71,7 @@ public: const Field& tmp = evec[i]; axpy(guess,TensorRemove(innerProduct(tmp,src)) / eval[i],tmp,guess); } + guess.checkerboard = src.checkerboard; } }; @@ -101,6 +102,7 @@ public: axpy(guess_coarse,TensorRemove(innerProduct(tmp,src_coarse)) / eval_coarse[i],tmp,guess_coarse); } blockPromote(guess_coarse,guess,subspace); + guess.checkerboard = src.checkerboard; }; }; diff --git a/lib/algorithms/iterative/SchurRedBlack.h b/lib/algorithms/iterative/SchurRedBlack.h index 091330b2..5abc4d9f 100644 --- a/lib/algorithms/iterative/SchurRedBlack.h +++ b/lib/algorithms/iterative/SchurRedBlack.h @@ -95,16 +95,26 @@ namespace Grid { private: OperatorFunction & _HermitianRBSolver; int CBfactorise; + bool subGuess; public: ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver) : + SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver) { CBfactorise=0; + subtractGuess(initSubGuess); }; + void subtractGuess(const bool initSubGuess) + { + subGuess = initSubGuess; + } + bool isSubtractGuess(void) + { + return subGuess; + } template void operator() (Matrix & _Matrix,const Field &in, Field &out){ @@ -150,9 +160,12 @@ namespace Grid { // Call the red-black solver ////////////////////////////////////////////////////////////// std::cout< using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve; @@ -184,15 +201,25 @@ namespace Grid { private: OperatorFunction & _HermitianRBSolver; int CBfactorise; + bool subGuess; public: ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - SchurRedBlackDiagMooeeSolve(OperatorFunction &HermitianRBSolver,int cb=0) : _HermitianRBSolver(HermitianRBSolver) + SchurRedBlackDiagMooeeSolve(OperatorFunction &HermitianRBSolver,int cb=0, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver) { CBfactorise=cb; + subtractGuess(initSubGuess); }; + void subtractGuess(const bool initSubGuess) + { + subGuess = initSubGuess; + } + bool isSubtractGuess(void) + { + return subGuess; + } template void operator() (Matrix & _Matrix,const Field &in, Field &out){ ZeroGuesser guess; @@ -236,7 +263,10 @@ namespace Grid { ////////////////////////////////////////////////////////////// std::cout< & _HermitianRBSolver; int CBfactorise; + bool subGuess; public: ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver) : + SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver) { - CBfactorise=0; + CBfactorise = 0; + subtractGuess(initSubGuess); }; + void subtractGuess(const bool initSubGuess) + { + subGuess = initSubGuess; + } + bool isSubtractGuess(void) + { + return subGuess; + } template void operator() (Matrix & _Matrix,const Field &in, Field &out){ @@ -322,8 +366,11 @@ namespace Grid { std::cout< & _HermitianRBSolver; int CBfactorise; + bool subGuess; public: ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - SchurRedBlackDiagTwoMixed(LinearFunction &HermitianRBSolver) : + SchurRedBlackDiagTwoMixed(LinearFunction &HermitianRBSolver, const bool initSubGuess = false) : _HermitianRBSolver(HermitianRBSolver) { CBfactorise=0; + subtractGuess(initSubGuess); }; + void subtractGuess(const bool initSubGuess) + { + subGuess = initSubGuess; + } + bool isSubtractGuess(void) + { + return subGuess; + } template void operator() (Matrix & _Matrix,const Field &in, Field &out){ @@ -408,7 +469,10 @@ namespace Grid { // _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); // _HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd); guess(src_o,tmp); - _HermitianRBSolver(src_o,tmp); assert(tmp.checkerboard==Odd); + Mtmp = tmp; + _HermitianRBSolver(_HermOpEO,src_o,tmp); assert(tmp.checkerboard==Odd); + // Fionn A2A boolean behavioural control + if (subGuess) tmp = tmp-Mtmp; _Matrix.MooeeInv(tmp,sol_o); assert( sol_o.checkerboard ==Odd); /////////////////////////////////////////////////// @@ -422,12 +486,16 @@ namespace Grid { setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); // Verify the unprec residual - _Matrix.M(out,resid); - resid = resid-in; - RealD ns = norm2(in); - RealD nr = norm2(resid); + if ( ! subGuess ) { + _Matrix.M(out,resid); + resid = resid-in; + RealD ns = norm2(in); + RealD nr = norm2(resid); - std::cout<::ExportPhysicalFermionSolution(const FermionField &so axpby_ssp_pplus (tmp, 1., tmp , 1., solution5d, 0, Ls-1); ExtractSlice(exported4d, tmp, 0, 0); } +template +void CayleyFermion5D::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d) +{ + int Ls = this->Ls; + FermionField tmp(this->FermionGrid()); + tmp = solution5d; + conformable(solution5d._grid,this->FermionGrid()); + conformable(exported4d._grid,this->GaugeGrid()); + axpby_ssp_pplus (tmp, 0., solution5d, 1., solution5d, 0, 0); + axpby_ssp_pminus(tmp, 1., tmp , 1., solution5d, 0, Ls-1); + ExtractSlice(exported4d, tmp, 0, 0); +} +template +void CayleyFermion5D::ImportUnphysicalFermion(const FermionField &input4d,FermionField &imported5d) +{ + int Ls = this->Ls; + FermionField tmp(this->FermionGrid()); + conformable(imported5d._grid,this->FermionGrid()); + conformable(input4d._grid ,this->GaugeGrid()); + tmp = zero; + InsertSlice(input4d, tmp, 0 , 0); + InsertSlice(input4d, tmp, Ls-1, 0); + axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, 0, 0); + axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1); + imported5d=tmp; +} + template void CayleyFermion5D::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d) { diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h index b370b09d..3bf2a8b6 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.h +++ b/lib/qcd/action/fermion/CayleyFermion5D.h @@ -89,7 +89,9 @@ namespace Grid { virtual void Dminus(const FermionField &psi, FermionField &chi); virtual void DminusDag(const FermionField &psi, FermionField &chi); virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d); + virtual void ExportPhysicalFermionSource(const FermionField &solution5d, FermionField &exported4d); virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d); + virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d); ///////////////////////////////////////////////////// // Instantiate different versions depending on Impl diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 77fdbec5..dc1ab924 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -162,10 +162,18 @@ namespace Grid { { imported = input; }; + virtual void ImportUnphysicalFermion(const FermionField &input,FermionField &imported) + { + imported=input; + }; virtual void ExportPhysicalFermionSolution(const FermionField &solution,FermionField &exported) { exported=solution; }; + virtual void ExportPhysicalFermionSource(const FermionField &solution,FermionField &exported) + { + exported=solution; + }; }; }