Add WIP SplitEven modules

Hadrons/Modules/MFermion/SplitEvenProp.cpp

@@ -0,0 +1,7 @@
+#include <Hadrons/Modules/MFermion/SplitEvenProp.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MFermion;
+
+template class Grid::Hadrons::MFermion::TSplitEvenProp<FIMPL>;

Hadrons/Modules/MFermion/SplitEvenProp.hpp

@@ -0,0 +1,148 @@
+#ifndef Hadrons_MFermion_SplitEvenProp_hpp_
+#define Hadrons_MFermion_SplitEvenProp_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         SplitEvenProp                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MFermion)
+
+class SplitEvenPropPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SplitEvenPropPar,
+                                    std::string,    q1,
+                                    // Must be a propagator solved on g5*src
+                                    std::string,    q2g5,
+                                    double,         m1,
+                                    double,         m2,
+                                    // TODO: Remove 'rotated'
+                                    bool,           rotated,
+                                    Gamma::Algebra, gamma);
+};
+
+template <typename FImpl>
+class TSplitEvenProp: public Module<SplitEvenPropPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+
+    // constructor
+    TSplitEvenProp(const std::string name);
+    // destructor
+    virtual ~TSplitEvenProp(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);
+private:
+    enum class Mode
+    {
+        SINGLE,
+        STDVECTOR
+    };
+
+    Mode mode;
+};
+
+MODULE_REGISTER_TMP(SplitEvenProp, TSplitEvenProp<FIMPL>, MFermion);
+
+/******************************************************************************
+ *                 TSplitEvenProp implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSplitEvenProp<FImpl>::TSplitEvenProp(const std::string name)
+: Module<SplitEvenPropPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSplitEvenProp<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q1, par().q2g5};
+
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSplitEvenProp<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSplitEvenProp<FImpl>::setup(void)
+{
+    if (envHasType(PropagatorField, par().q1))
+    {
+        if (!envHasType(PropagatorField, par().q2g5))
+            {HADRONS_ERROR(Definition, "q1 is a 'PropagatorField' and g5q2 is not");}
+        envCreateLat(PropagatorField, getName());
+        this->mode = Mode::SINGLE;
+    }
+    else if (envHasType(std::vector<PropagatorField>, par().q1))
+    {
+        if (!envHasType(std::vector<PropagatorField>, par().q2g5))
+            {HADRONS_ERROR(Definition, "q1 is a 'std::vector<PropagatorField>' and g5q2 is not");}
+        std::vector<PropagatorField>& q1 = envGet(std::vector<PropagatorField>, par().q1);
+        envCreate(std::vector<PropagatorField>, getName(), 1, q1.size(), envGetGrid(PropagatorField));
+        this->mode = Mode::STDVECTOR;
+    }
+    else
+    {
+        HADRONS_ERROR(Definition, "q1 and q2g5 are neither 'PropagatorField's nor 'std::vector<PropagatorField>'s");
+    }
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+#define SplitEven(q1, q2g5, g5, gamma, massfactor) massfactor*adj(g5*q2g5)*gamma*q1
+#define SplitEvenRotated(q1, q2g5, g5, gamma, massfactor) massfactor*q1*adj(g5*q2g5)*gamma
+template <typename FImpl>
+void TSplitEvenProp<FImpl>::execute(void)
+{
+    Gamma g5(Gamma::Algebra::Gamma5);
+    Gamma gamma(par().gamma);
+    double massfactor = par().m2 - par().m1;
+    if (this->mode == Mode::SINGLE)
+    {
+        PropagatorField& q1   = envGet(PropagatorField, par().q1);
+        PropagatorField& q2g5 = envGet(PropagatorField, par().q2g5);
+        PropagatorField& out  = envGet(PropagatorField, getName());
+
+        // TODO: Remove 'rotated' - hack for regression
+        if (par().rotated) out = SplitEvenRotated(q1, q2g5, g5, gamma, massfactor);
+        else               out = SplitEven       (q1, q2g5, g5, gamma, massfactor);
+
+    }
+    else if (this->mode == Mode::STDVECTOR)
+    {
+        std::vector<PropagatorField>& q1   = envGet(std::vector<PropagatorField>, par().q1);
+        std::vector<PropagatorField>& q2g5 = envGet(std::vector<PropagatorField>, par().q2g5);
+        std::vector<PropagatorField>& out  = envGet(std::vector<PropagatorField>, getName());
+        for (int i=0; i < q1.size(); ++i)
+        {
+            // TODO: Remove 'rotated' - hack for regression
+            if (par().rotated) out[i] = SplitEvenRotated(q1[i], q2g5[i], g5, gamma, massfactor);
+            else               out[i] = SplitEven       (q1[i], q2g5[i], g5, gamma, massfactor);
+        }
+    }
+}
+#undef SplitEven
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MFermion_SplitEvenProp_hpp_

Hadrons/Modules/MUtilities/MulGamma.cpp

@@ -0,0 +1,7 @@
+#include <Hadrons/Modules/MUtilities/MulGamma.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MUtilities;
+
+template class HADRONS_NAMESPACE::MUtilities::TMulGamma<FIMPL>;

Hadrons/Modules/MUtilities/MulGamma.hpp

@@ -0,0 +1,140 @@
+#ifndef Hadrons_MUtilities_MulGamma_hpp_
+#define Hadrons_MUtilities_MulGamma_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         MulGamma                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MUtilities)
+
+class MulGammaPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(MulGammaPar,
+                                    Gamma::Algebra, gammaleft,
+                                    std::string,    q,
+                                    Gamma::Algebra, gammaright);
+};
+
+template <typename FImpl>
+class TMulGamma: public Module<MulGammaPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    // constructor
+    TMulGamma(const std::string name);
+    // destructor
+    virtual ~TMulGamma(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);
+private:
+    enum class Mode
+    {
+        SINGLE,
+        STDVECTOR
+    };
+
+    Mode mode;
+};
+
+MODULE_REGISTER_TMP(MulGamma, TMulGamma<FIMPL>, MUtilities);
+
+/******************************************************************************
+ *                 TMulGamma implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TMulGamma<FImpl>::TMulGamma(const std::string name)
+: Module<MulGammaPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TMulGamma<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q};
+
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TMulGamma<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMulGamma<FImpl>::setup(void)
+{
+    if (envHasType(PropagatorField, par().q))
+    {
+        LOG(Message) << "SINGLE PROP FIELD" << std::endl;
+        envCreateLat(PropagatorField, getName());
+        this->mode = Mode::SINGLE;
+    }
+    else if (envHasType(std::vector<PropagatorField>, par().q))
+    {
+        LOG(Message) << "VECTOR PROP FIELD" << std::endl;
+        std::vector<PropagatorField>& q = envGet(std::vector<PropagatorField>, par().q);
+        envCreate(std::vector<PropagatorField>, getName(), 1, q.size(), envGetGrid(PropagatorField));
+        this->mode = Mode::STDVECTOR;
+    }
+    else
+    {
+        HADRONS_ERROR(Definition, "q is neither a 'PropagatorField' nor a 'std::vector<PropagatorField>'");
+    }
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMulGamma<FImpl>::execute(void)
+{
+    Gamma gammaleft (par().gammaleft);
+    Gamma gammaright(par().gammaright);
+
+    if (this->mode == Mode::SINGLE)
+    {
+        PropagatorField& q   = envGet(PropagatorField, par().q);
+        PropagatorField& out  = envGet(PropagatorField, getName());
+
+        // Optimise out identity multiplications
+        if (par().gammaleft == Gamma::Algebra::Identity)
+            out = q*gammaright;
+        else if (par().gammaright == Gamma::Algebra::Identity)
+            out = gammaleft*q;
+        else
+            out = gammaleft*q*gammaright;
+    }
+    else if (this->mode == Mode::STDVECTOR)
+    {
+        std::vector<PropagatorField>& q   = envGet(std::vector<PropagatorField>, par().q);
+        std::vector<PropagatorField>& out = envGet(std::vector<PropagatorField>, getName());
+
+        // Optimise out identity multiplications
+        if (par().gammaleft == Gamma::Algebra::Identity)
+            for (int i=0; i < q.size(); ++i) out[i] = q[i]*gammaright;
+        else if (par().gammaright == Gamma::Algebra::Identity)
+            for (int i=0; i < q.size(); ++i) out[i] = gammaleft*q[i];
+        else
+            for (int i=0; i < q.size(); ++i) out[i] = gammaleft*q[i]*gammaright;
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MUtilities_MulGamma_hpp_