From 8c3cc3236447b4b8eef95a29da1b48166b5eb03d Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Thu, 29 Dec 2016 22:42:58 +0100
Subject: [PATCH] Scalar action

---
 lib/qcd/action/Actions.h       |   5 +
 lib/qcd/action/scalar/Scalar.h | 211 +++++++++++++++++++++++++++++++++
 2 files changed, 216 insertions(+)
 create mode 100644 lib/qcd/action/scalar/Scalar.h

diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h
index fea75f8a..efd6a5bc 100644
--- a/lib/qcd/action/Actions.h
+++ b/lib/qcd/action/Actions.h
@@ -292,4 +292,9 @@ typedef MobiusFermion<GparityWilsonImplD> GparityMobiusFermionD;
 #include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h>
 #include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h>
 
+////////////////////
+// Scalar actions
+////////////////////
+#include <Grid/qcd/action/scalar/Scalar.h>
+
 #endif
diff --git a/lib/qcd/action/scalar/Scalar.h b/lib/qcd/action/scalar/Scalar.h
new file mode 100644
index 00000000..194f6767
--- /dev/null
+++ b/lib/qcd/action/scalar/Scalar.h
@@ -0,0 +1,211 @@
+#ifndef QCD_SCALAR_ACTION_H
+#define QCD_SCALAR_ACTION_H
+
+#define INHERIT_SIMPL_TYPES(Impl)\
+typedef typename Impl::SiteScalar      SiteScalar;	 \
+typedef typename Impl::SiteSpinor      SiteSpinor;	 \
+typedef typename Impl::SitePropagator  SitePropagator;   \
+typedef typename Impl::ScalarField     ScalarField;	 \
+typedef typename Impl::FermionField    FermionField;	 \
+typedef typename Impl::PropagatorField PropagatorField;  \
+typedef typename Impl::StencilImpl     StencilImpl;
+
+namespace Grid{
+namespace QCD{
+  // Scalar implementation class ///////////////////////////////////////////////
+  // FIXME: it is not very nice to have the FImpl aliases
+  template <class S,
+            class Representation = FundamentalRep<1>,
+            class _Coeff_t = RealD>
+  class ScalarImpl:
+    public PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension>>
+  {
+  public:
+    static constexpr unsigned int rDim = Representation::Dimension;
+  public:
+    // gauge types
+    typedef PeriodicGaugeImpl<GaugeImplTypes<S, rDim>> Gimpl;
+    INHERIT_GIMPL_TYPES(Gimpl);
+    // site types
+    // (using classes instead of aliases to allow for partial specialisation)
+    template <typename vtype, unsigned int d>
+    class iImplScalar
+    {
+    public:
+      typedef iScalar<iScalar<iVector<vtype, d>>> type;
+    };
+    template <typename vtype>
+    class iImplScalar<vtype, 1>
+    {
+    public:
+      typedef iScalar<iScalar<iScalar<vtype>>> type;
+    };
+    template <typename vtype, unsigned int d>
+    class iImplPropagator
+    {
+    public:
+      typedef iScalar<iScalar<iMatrix<vtype, d>>> type;
+    };
+    template <typename vtype>
+    class iImplPropagator<vtype, 1>
+    {
+    public:
+      typedef iScalar<iScalar<iScalar<vtype>>> type;
+    };
+    // type aliases
+    typedef typename iImplScalar<S, rDim>::type      SiteScalar;
+    typedef SiteScalar                               SiteSpinor;
+    typedef typename iImplPropagator<S, rDim>::type  SitePropagator;
+    typedef Lattice<SiteScalar>                      ScalarField;
+    typedef ScalarField                              FermionField;
+    typedef Lattice<SitePropagator>                  PropagatorField;
+    typedef CartesianStencil<SiteScalar, SiteScalar> StencilImpl;
+  };
+  
+  // single scalar implementation
+  typedef ScalarImpl<vComplex> ScalarImplR;
+  
+  // Scalar action /////////////////////////////////////////////////////////////
+  template <typename SImpl>
+  class Scalar:
+    public CheckerBoardedSparseMatrixBase<typename SImpl::ScalarField>,
+    public SImpl
+  {
+  public:
+    INHERIT_GIMPL_TYPES(SImpl);
+    INHERIT_SIMPL_TYPES(SImpl);
+  public:
+    // constructor
+    Scalar(GaugeField &_Umu, GridCartesian &Sgrid, GridRedBlackCartesian &Hgrid,
+           RealD _mass)
+    : _grid(&Sgrid)
+    , _cbgrid(&Hgrid)
+    , mass(_mass)
+    , Lebesgue(_grid)
+    , LebesgueEvenOdd(_cbgrid)
+    , Umu(&Sgrid)
+    , UmuEven(&Hgrid)
+    , UmuOdd(&Hgrid)
+    {
+      Umu = _Umu;
+      pickCheckerboard(Even, UmuEven, Umu);
+      pickCheckerboard(Odd, UmuOdd, Umu);
+    }
+    // grid access
+    virtual GridBase *RedBlackGrid(void) {return _grid;}
+    // half checkerboard operations
+    // FIXME: do implementation
+    virtual void Meooe(const ScalarField &in, ScalarField &out)
+    {
+      assert(0);
+    }
+    virtual void Mooee(const ScalarField &in, ScalarField &out)
+    {
+      assert(0);
+    }
+    virtual void MooeeInv(const ScalarField &in, ScalarField &out)
+    {
+      assert(0);
+    }
+    virtual void MeooeDag(const ScalarField &in, ScalarField &out)
+    {
+      assert(0);
+    }
+    virtual void MooeeDag(const ScalarField &in, ScalarField &out)
+    {
+      assert(0);
+    }
+    virtual void MooeeInvDag(const ScalarField &in, ScalarField &out)
+    {
+      assert(0);
+    }
+    // free propagators
+    static void MomentumSpacePropagator(ScalarField &out, RealD m);
+    static void FreePropagator(const ScalarField &in, ScalarField &out,
+                               const ScalarField &momKernel);
+    static void FreePropagator(const ScalarField &in, ScalarField &out, RealD m);
+  public:
+    RealD mass;
+    
+    GridBase *_grid;
+    GridBase *_cbgrid;
+    
+    // Defines the stencils for even and odd
+    StencilImpl Stencil;
+    StencilImpl StencilEven;
+    StencilImpl StencilOdd;
+    
+    // Copy of the gauge field, with even and odd subsets
+    GaugeField Umu;
+    GaugeField UmuEven;
+    GaugeField UmuOdd;
+    
+    LebesgueOrder Lebesgue;
+    LebesgueOrder LebesgueEvenOdd;
+  };
+  
+  template <typename SImpl>
+  void Scalar<SImpl>::MomentumSpacePropagator(ScalarField &out, RealD m)
+  {
+    GridBase           *grid = out._grid;
+    ScalarField        kmu(grid);
+    const unsigned int nd    = grid->_ndimension;
+    std::vector<int>   &l    = grid->_fdimensions;
+    
+    out = m*m;
+    for(int mu = 0; mu < nd; mu++)
+    {
+      Real twoPiL = M_PI*2./l[mu];
+      
+      LatticeCoordinate(kmu,mu);
+      kmu = 2.*sin(.5*twoPiL*kmu);
+      out = out + kmu*kmu;
+    }
+  }
+  
+  template <typename SImpl>
+  void Scalar<SImpl>::FreePropagator(const ScalarField &in, ScalarField &out,
+                                     const ScalarField &FTKernel)
+  {
+    FFT         fft((GridCartesian *)in._grid);
+    ScalarField inFT(in._grid);
+    
+    fft.FFT_all_dim(inFT, in, FFT::forward);
+    inFT = inFT*FTKernel;
+    fft.FFT_all_dim(out, inFT, FFT::backward);
+  }
+  
+  template <typename SImpl>
+  void Scalar<SImpl>::FreePropagator(const ScalarField &in, ScalarField &out,
+                                     RealD m)
+  {
+    ScalarField FTKernel(in._grid);
+    
+    MomentumSpacePropagator(FTKernel, m);
+    FreePropagator(in, out, FTKernel);
+  }
+  
+  template <class SImpl>
+  void ScalarToProp(typename SImpl::PropagatorField &p,
+                    const typename SImpl::ScalarField &s,
+                    const int c)
+  {
+    for(int i = 0; i < SImpl::rDim; ++i)
+    {
+      pokeColour(p, peekColour(s, i), i);
+    }
+  }
+  
+  template <class SImpl>
+  void PropToScalar(typename SImpl::ScalarField &s,
+                    const typename SImpl::PropagatorField &p,
+                    const int c)
+  {
+    for(int i = 0; i < SImpl::rDim; ++i)
+    {
+      pokeColour(s, peekColour(p, i), i);
+    }
+  }
+}}
+
+#endif