Merge branch 'feature/scalar_adjointFT' into feature/hadrons

# Conflicts:
#	lib/qcd/action/scalar/ScalarImpl.h

lib/qcd/action/scalar/Scalar.h

@@ -31,6 +31,7 @@ directory
 
 #include <Grid/qcd/action/scalar/ScalarImpl.h>
 #include <Grid/qcd/action/scalar/ScalarAction.h>
+#include <Grid/qcd/action/scalar/ScalarInteractionAction.h>
 
 namespace Grid {
 namespace QCD {
@@ -39,6 +40,10 @@ namespace QCD {
   typedef ScalarAction<ScalarImplF>                 ScalarActionF;
   typedef ScalarAction<ScalarImplD>                 ScalarActionD;
 
+  template <int Colours, int Dimensions> using ScalarAdjActionR = ScalarInteractionAction<ScalarNxNAdjImplR<Colours>, Dimensions>;
+  template <int Colours, int Dimensions> using ScalarAdjActionF = ScalarInteractionAction<ScalarNxNAdjImplF<Colours>, Dimensions>;
+  template <int Colours, int Dimensions> using ScalarAdjActionD = ScalarInteractionAction<ScalarNxNAdjImplD<Colours>, Dimensions>;
+  
 }
 }
 

lib/qcd/action/scalar/ScalarAction.h

@@ -46,20 +46,17 @@ namespace Grid {
     RealD lambda;
 
  public:
-    ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){};
+    ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l) {}
 
     virtual std::string LogParameters() {
       std::stringstream sstream;
       sstream << GridLogMessage << "[ScalarAction] lambda      : " << lambda      << std::endl;
       sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl;
       return sstream.str();
-      
     }
-    
     virtual std::string action_name() {return "ScalarAction";}
 
-    virtual void refresh(const Field &U,
+    virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}  // noop as no pseudoferms
-			 GridParallelRNG &pRNG){};  // noop as no pseudoferms
 
     virtual RealD S(const Field &p) {
       return (mass_square * 0.5 + QCD::Nd) * ScalarObs<Impl>::sumphisquared(p) +
@@ -76,9 +73,11 @@ namespace Grid {
       for (int mu = 1; mu < QCD::Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
 
       force =+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp;
-    };
+    }
 };
 
-} // Grid
+
+
+}  // namespace Grid
 
 #endif // SCALAR_ACTION_H

lib/qcd/action/scalar/ScalarImpl.h

@@ -88,20 +88,18 @@ namespace Grid {
   };
 
   template <class S, unsigned int N>
-  class ScalarMatrixImplTypes {
+  class ScalarAdjMatrixImplTypes {
   public:
     typedef S Simd;
-    
     template <typename vtype>
     using iImplField = iScalar<iScalar<iMatrix<vtype, N> > >;
 
     typedef iImplField<Simd> SiteField;
 
-    
     typedef Lattice<SiteField> Field;
 
     static inline void generate_momenta(Field& P, GridParallelRNG& pRNG) {
-      gaussian(pRNG, P);
+      QCD::SU<N>::GaussianFundamentalLieAlgebraMatrix(pRNG, P);
     }
 
     static inline Field projectForce(Field& P) {return P;}
@@ -111,19 +109,19 @@ namespace Grid {
     }
 
     static inline RealD FieldSquareNorm(Field& U) {
-      return (TensorRemove(- sum(trace(U*U))*0.5).real());
+      return (TensorRemove(sum(trace(U*U))).real());
     }
 
     static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
-      gaussian(pRNG, U);
+      QCD::SU<N>::LieRandomize(pRNG, U);
     }
 
     static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
-      gaussian(pRNG, U);
+      QCD::SU<N>::LieRandomize(pRNG, U, 0.01);
     }
 
     static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
-      U = 1.0;
+      U = zero;
     }
 
   };
@@ -138,6 +136,15 @@ namespace Grid {
   typedef ScalarImplTypes<vComplexF> ScalarImplCF;
   typedef ScalarImplTypes<vComplexD> ScalarImplCD;
     
+  // Hardcoding here the size of the matrices
+  typedef ScalarAdjMatrixImplTypes<vComplex,  QCD::Nc> ScalarAdjImplR;
+  typedef ScalarAdjMatrixImplTypes<vComplexF, QCD::Nc> ScalarAdjImplF;
+  typedef ScalarAdjMatrixImplTypes<vComplexD, QCD::Nc> ScalarAdjImplD;
+
+  template <int Colours > using ScalarNxNAdjImplR = ScalarAdjMatrixImplTypes<vComplex,   Colours >;
+  template <int Colours > using ScalarNxNAdjImplF = ScalarAdjMatrixImplTypes<vComplexF,  Colours >;
+  template <int Colours > using ScalarNxNAdjImplD = ScalarAdjMatrixImplTypes<vComplexD,  Colours >;
+  
   //}
 }
 

lib/qcd/action/scalar/ScalarInteractionAction.h

@@ -6,10 +6,7 @@
 
   Copyright (C) 2015
 
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+  Author: Guido Cossu <guido,cossu@ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: neo <cossu@post.kek.jp>
-Author: paboyle <paboyle@ph.ed.ac.uk>
 
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
@@ -30,55 +27,122 @@ directory
   *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef SCALAR_ACTION_H
+#ifndef SCALAR_INT_ACTION_H
-#define SCALAR_ACTION_H
+#define SCALAR_INT_ACTION_H
+
+
+// Note: this action can completely absorb the ScalarAction for real float fields
+// use the scalarObjs to generalise the structure
 
 namespace Grid {
   // FIXME drop the QCD namespace everywhere here
 
-  template <class Impl>
+  template <class Impl, int Ndim >
   class ScalarInteractionAction : public QCD::Action<typename Impl::Field> {
   public:
     INHERIT_FIELD_TYPES(Impl);
-    
   private:
     RealD mass_square;
     RealD lambda;
 
+
+    typedef typename Field::vector_object vobj;
+    typedef CartesianStencil<vobj,vobj> Stencil;
+
+    SimpleCompressor<vobj> compressor;
+    int npoint = 2*Ndim;
+    std::vector<int> directions;//    = {0,1,2,3,0,1,2,3};  // forcing 4 dimensions
+    std::vector<int> displacements;//  = {1,1,1,1, -1,-1,-1,-1};
+
+
   public:
-    ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){};
+
+    ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l), displacements(2*Ndim,0), directions(2*Ndim,0){
+      for (int mu = 0 ; mu < Ndim; mu++){
+		directions[mu]         = mu; directions[mu+Ndim]    = mu;
+		displacements[mu]      =  1; displacements[mu+Ndim] = -1;
+      }
+    }
 
     virtual std::string LogParameters() {
       std::stringstream sstream;
       sstream << GridLogMessage << "[ScalarAction] lambda      : " << lambda      << std::endl;
       sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl;
       return sstream.str();
-      
     }
 
     virtual std::string action_name() {return "ScalarAction";}
 
-    virtual void refresh(const Field &U,
+    virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}
-			 GridParallelRNG &pRNG){};  // noop as no pseudoferms
 
     virtual RealD S(const Field &p) {
-      return (mass_square * 0.5 + QCD::Nd) * ScalarObs<Impl>::sumphisquared(p) +
+      assert(p._grid->Nd() == Ndim);
-	(lambda / 24.) * ScalarObs<Impl>::sumphifourth(p) +
+      static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
-	ScalarObs<Impl>::sumphider(p);
+      phiStencil.HaloExchange(p, compressor);
+      Field action(p._grid), pshift(p._grid), phisquared(p._grid);
+      phisquared = p*p;
+      action = (2.0*Ndim + mass_square)*phisquared + lambda*phisquared*phisquared;
+      for (int mu = 0; mu < Ndim; mu++) {
+	//  pshift = Cshift(p, mu, +1);  // not efficient, implement with stencils
+	parallel_for (int i = 0; i < p._grid->oSites(); i++) {
+	  int permute_type;
+	  StencilEntry *SE;
+	  vobj temp2;
+	  const vobj *temp, *t_p;
+	    
+	  SE = phiStencil.GetEntry(permute_type, mu, i);
+	  t_p  = &p._odata[i];
+	  if ( SE->_is_local ) {
+	    temp = &p._odata[SE->_offset];
+	    if ( SE->_permute ) {
+	      permute(temp2, *temp, permute_type);
+	      action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2;
+	    } else {
+	      action._odata[i] -= *temp*(*t_p) + (*t_p)*(*temp);
+	    }
+	  } else {
+	    action._odata[i] -= phiStencil.CommBuf()[SE->_offset]*(*t_p) + (*t_p)*phiStencil.CommBuf()[SE->_offset];
+	  }
+	}
+	//  action -= pshift*p + p*pshift;
+      }
+      // NB the trace in the algebra is normalised to 1/2
+      // minus sign coming from the antihermitian fields
+      return -(TensorRemove(sum(trace(action)))).real();
     };
 
-    virtual void deriv(const Field &p,
+    virtual void deriv(const Field &p, Field &force) {
-		       Field &force) {
+      assert(p._grid->Nd() == Ndim);
-      Field tmp(p._grid);
+      force = (2.0*Ndim + mass_square)*p + 2.0*lambda*p*p*p;
-      Field p2(p._grid);
+      // move this outside
-      ScalarObs<Impl>::phisquared(p2, p);
+      static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
-      tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1));
+      phiStencil.HaloExchange(p, compressor);
-      for (int mu = 1; mu < QCD::Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
       
-      force=+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp;
+      //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
-    };
+      for (int point = 0; point < npoint; point++) {
+	parallel_for (int i = 0; i < p._grid->oSites(); i++) {
+	  const vobj *temp;
+	  vobj temp2;
+	  int permute_type;
+	  StencilEntry *SE;
+	  SE = phiStencil.GetEntry(permute_type, point, i);
+	  
+	  if ( SE->_is_local ) {
+	    temp = &p._odata[SE->_offset];
+	    if ( SE->_permute ) {
+	      permute(temp2, *temp, permute_type);
+	      force._odata[i] -= temp2;
+	    } else {
+	      force._odata[i] -= *temp;
+	    }
+	  } else {
+	    force._odata[i] -= phiStencil.CommBuf()[SE->_offset];
+	  }
+	}
+      }
+    }
   };
   
-} // Grid
+}  // namespace Grid
 
-#endif // SCALAR_ACTION_H
+#endif  // SCALAR_INT_ACTION_H

lib/qcd/hmc/GenericHMCrunner.h

@@ -207,6 +207,12 @@ using GenericHMCRunnerTemplate = HMCWrapperTemplate<Implementation, Integrator,
 typedef HMCWrapperTemplate<ScalarImplR, MinimumNorm2, ScalarFields>
     ScalarGenericHMCRunner;
 
+typedef HMCWrapperTemplate<ScalarAdjImplR, MinimumNorm2, ScalarMatrixFields>
+    ScalarAdjGenericHMCRunner;
+
+template <int Colours> 
+using ScalarNxNAdjGenericHMCRunner = HMCWrapperTemplate < ScalarNxNAdjImplR<Colours>, MinimumNorm2, ScalarNxNMatrixFields<Colours> >;
+
 }  // namespace QCD
 }  // namespace Grid
 

lib/qcd/representations/hmc_types.h

@@ -62,7 +62,10 @@ class Representations {
 
 typedef Representations<FundamentalRepresentation> NoHirep;
 typedef Representations<EmptyRep<typename ScalarImplR::Field> > ScalarFields;
-  //typedef Representations<EmptyRep<typename ScalarMatrixImplR::Field> > ScalarMatrixFields;
+typedef Representations<EmptyRep<typename ScalarAdjImplR::Field> > ScalarMatrixFields;
+
+template < int Colours> 
+using ScalarNxNMatrixFields = Representations<EmptyRep<typename ScalarNxNAdjImplR<Colours>::Field> >;
 
 // Helper classes to access the elements
 // Strips the first N parameters from the tuple

tests/Test_stencil.cc

@@ -33,8 +33,7 @@ using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
-int main (int argc, char ** argv)
+int main(int argc, char ** argv) {
-{
   Grid_init(&argc, &argv);
 
   //  typedef LatticeColourMatrix Field;
@@ -190,7 +189,6 @@ int main (int argc, char ** argv)
 
 	SimpleCompressor<vobj> compress;
 
-	
 	Bar = Cshift(Foo,dir,disp);
 
 	if ( disp & 0x1 ) {

tests/hmc/Test_hmc_ScalarActionNxN.cc

@@ -0,0 +1,193 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2016
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Grid.h>
+namespace Grid {
+class ScalarActionParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarActionParameters,
+    double, mass_squared,
+    double, lambda);
+
+    template <class ReaderClass >
+  ScalarActionParameters(Reader<ReaderClass>& Reader){
+    read(Reader, "ScalarAction", *this);
+  }
+
+};
+}
+
+using namespace Grid;
+using namespace Grid::QCD;
+
+template <class Impl>
+class MagLogger : public HmcObservable<typename Impl::Field> {
+public:
+  typedef typename Impl::Field Field;
+  typedef typename Impl::Simd::scalar_type Trace;
+  
+  void TrajectoryComplete(int traj,
+                          Field &U,
+                          GridSerialRNG &sRNG,
+                          GridParallelRNG &pRNG) {
+    
+    int def_prec = std::cout.precision();
+    
+    std::cout << std::setprecision(std::numeric_limits<Real>::digits10 + 1);
+    std::cout << GridLogMessage
+              << "m= " << TensorRemove(trace(sum(U))) << std::endl;
+    std::cout << GridLogMessage
+              << "m^2= " << TensorRemove(trace(sum(U)*sum(U))) << std::endl;
+    std::cout << GridLogMessage
+    << "phi^2= " << TensorRemove(sum(trace(U*U))) << std::endl;
+    std::cout.precision(def_prec);
+    
+  }
+private:
+  
+};
+
+template <class Impl>
+class MagMod: public ObservableModule<MagLogger<Impl>, NoParameters>{
+  typedef ObservableModule<MagLogger<Impl>, NoParameters> ObsBase;
+  using ObsBase::ObsBase; // for constructors
+  
+  // acquire resource
+  virtual void initialize(){
+    this->ObservablePtr.reset(new MagLogger<Impl>());
+  }
+public:
+  MagMod(): ObsBase(NoParameters()){}
+};
+
+int main(int argc, char **argv) {
+  typedef Grid::JSONReader       Serialiser;
+  
+  Grid_init(&argc, &argv);
+  int threads = GridThread::GetThreads();
+  // here make a routine to print all the relevant information on the run
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  // Typedefs to simplify notation
+  constexpr int Ncolours    = 2;
+  constexpr int Ndimensions = 3;
+  typedef ScalarNxNAdjGenericHMCRunner<Ncolours> HMCWrapper;  // Uses the default minimum norm, real scalar fields
+  typedef ScalarAdjActionR<Ncolours, Ndimensions> ScalarAction;
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+  TheHMC.ReadCommandLine(argc, argv);
+
+  if (TheHMC.ParameterFile.empty()){
+    std::cout << "Input file not specified."
+              << "Use --ParameterFile option in the command line.\nAborting" 
+              << std::endl;
+    exit(1);
+  }
+  Serialiser Reader(TheHMC.ParameterFile);
+
+  // Grid from the command line
+  GridModule ScalarGrid;
+  if (GridDefaultLatt().size() != Ndimensions){
+    std::cout << "Incorrect dimension of the grid\n. Expected dim="<< Ndimensions << std::endl;
+    exit(1);
+  }
+  if (GridDefaultMpi().size() != Ndimensions){
+    std::cout << "Incorrect dimension of the mpi grid\n. Expected dim="<< Ndimensions << std::endl;
+    exit(1);
+  }
+  ScalarGrid.set_full(new GridCartesian(GridDefaultLatt(),GridDefaultSimd(Ndimensions, vComplex::Nsimd()),GridDefaultMpi()));
+  ScalarGrid.set_rb(new GridRedBlackCartesian(ScalarGrid.get_full()));
+  TheHMC.Resources.AddGrid("scalar", ScalarGrid);
+  std::cout << "Lattice size : " << GridDefaultLatt() << std::endl;
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams(Reader);
+  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar(Reader);
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+  
+  // Construct observables
+  typedef MagMod<HMCWrapper::ImplPolicy> MagObs;
+  TheHMC.Resources.AddObservable<MagObs>();
+  
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+
+  // Scalar action in adjoint representation
+  ScalarActionParameters SPar(Reader);
+  ScalarAction Saction(SPar.mass_squared, SPar.lambda);
+
+  // Collect actions
+  ActionLevel<ScalarAction::Field, ScalarNxNMatrixFields<Ncolours>> Level1(1);
+  Level1.push_back(&Saction);
+  TheHMC.TheAction.push_back(Level1);
+  /////////////////////////////////////////////////////////////
+  TheHMC.Parameters.initialize(Reader);
+
+  TheHMC.Run();
+
+  Grid_finalize();
+}  // main
+
+/* Examples for input files
+
+JSON
+
+{
+    "Checkpointer": {
+    "config_prefix": "ckpoint_scalar_lat",
+    "rng_prefix": "ckpoint_scalar_rng",
+    "saveInterval": 1,
+    "format": "IEEE64BIG"
+    },
+    "RandomNumberGenerator": {
+    "serial_seeds": "1 2 3 4 6",
+    "parallel_seeds": "6 7 8 9 11"
+    },
+    "ScalarAction":{
+      "mass_squared": 0.5,
+      "lambda": 0.1
+    },
+    "HMC":{
+    "StartTrajectory": 0,
+    "Trajectories": 100,
+    "MetropolisTest": true,
+    "NoMetropolisUntil": 10,
+    "StartingType": "HotStart",
+    "MD":{
+        "name": "MinimumNorm2",
+	      "MDsteps": 15,
+	      "trajL": 2.0
+	    }
+    }
+}
+
+
+XML example not provided yet
+
+*/