Merge branch 'develop' into hisq_fat_links

Grid/qcd/action/ActionBase.h

@@ -34,10 +34,24 @@ directory
 
 NAMESPACE_BEGIN(Grid);
 
+///////////////////////////////////
+// Smart configuration base class
+///////////////////////////////////
+template< class Field >
+class ConfigurationBase
+{
+public:
+  ConfigurationBase() {}
+  virtual ~ConfigurationBase() {}
+  virtual void set_Field(Field& U) =0;
+  virtual void smeared_force(Field&) const = 0;
+  virtual Field& get_SmearedU() =0;
+  virtual Field &get_U(bool smeared = false) = 0;
+};
+
 template <class GaugeField >
 class Action 
 {
-
 public:
   bool is_smeared = false;
   RealD deriv_norm_sum;
@@ -77,11 +91,39 @@ public:
   void refresh_timer_stop(void)  { refresh_us+=usecond(); }
   void S_timer_start(void)       { S_us-=usecond(); }
   void S_timer_stop(void)        { S_us+=usecond(); }
+  /////////////////////////////
   // Heatbath?
+  /////////////////////////////
   virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
   virtual RealD S(const GaugeField& U) = 0;                             // evaluate the action
   virtual RealD Sinitial(const GaugeField& U) { return this->S(U); } ;  // if the refresh computes the action, can cache it. Alternately refreshAndAction() ?
   virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0;        // evaluate the action derivative
+
+  /////////////////////////////////////////////////////////////
+  // virtual smeared interface through configuration container
+  /////////////////////////////////////////////////////////////
+  virtual void refresh(ConfigurationBase<GaugeField> & U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
+  {
+    refresh(U.get_U(is_smeared),sRNG,pRNG);
+  }
+  virtual RealD S(ConfigurationBase<GaugeField>& U)
+  {
+    return S(U.get_U(is_smeared));
+  }
+  virtual RealD Sinitial(ConfigurationBase<GaugeField>& U) 
+  {
+    return Sinitial(U.get_U(is_smeared));
+  }
+  virtual void deriv(ConfigurationBase<GaugeField>& U, GaugeField& dSdU)
+  {
+    deriv(U.get_U(is_smeared),dSdU); 
+    if ( is_smeared ) {
+      U.smeared_force(dSdU);
+    }
+  }
+  ///////////////////////////////
+  // Logging
+  ///////////////////////////////
   virtual std::string action_name()    = 0;                             // return the action name
   virtual std::string LogParameters()  = 0;                             // prints action parameters
   virtual ~Action(){}

Grid/qcd/action/ActionCore.h

@@ -30,6 +30,8 @@ directory
 #ifndef QCD_ACTION_CORE
 #define QCD_ACTION_CORE
 
+#include <Grid/qcd/action/gauge/GaugeImplementations.h>
+
 #include <Grid/qcd/action/ActionBase.h>
 NAMESPACE_CHECK(ActionBase);
 #include <Grid/qcd/action/ActionSet.h>

Grid/qcd/smearing/GaugeConfiguration.h

@@ -7,23 +7,10 @@
 
 NAMESPACE_BEGIN(Grid);
 
-template< class Impl >
-class ConfigurationBase
-{
-public:
-  INHERIT_FIELD_TYPES(Impl);
-
-  ConfigurationBase() {}
-  virtual ~ConfigurationBase() {}
-  virtual void set_Field(Field& U) =0;
-  virtual void smeared_force(Field&) const = 0;
-  virtual Field& get_SmearedU() =0;
-  virtual Field &get_U(bool smeared = false) = 0;
-};
 
 //trivial class for no smearing
 template< class Impl >
-class NoSmearing : public ConfigurationBase<Impl>
+class NoSmearing : public ConfigurationBase<typename Impl::Field>
 {
 public:
   INHERIT_FIELD_TYPES(Impl);
@@ -56,7 +43,7 @@ public:
   It stores a list of smeared configurations.
 */
 template <class Gimpl>
-class SmearedConfiguration : public ConfigurationBase<Gimpl>
+class SmearedConfiguration : public ConfigurationBase<typename Gimpl::Field>
 {
 public:
   INHERIT_GIMPL_TYPES(Gimpl);

Grid/qcd/smearing/GaugeConfigurationMasked.h

@@ -422,8 +422,7 @@ public:
     InsertForce(Fdet1,Fdet1_mu,mu);
     InsertForce(Fdet2,Fdet2_mu,mu);
 
-    force = Fdet1 + Fdet2;
+    force= (-0.5)*( Fdet1 + Fdet2);
-    
   }
   RealD logDetJacobianLevel(const GaugeField &U,int smr)
   {
@@ -534,27 +533,82 @@ public:
     RealD ln_det = 0;
     if (this->smearingLevels > 0)
     {
+      double start = usecond();
       for (int ismr = this->smearingLevels - 1; ismr > 0; --ismr) {
 	ln_det+= logDetJacobianLevel(this->get_smeared_conf(ismr-1),ismr);
       }
       ln_det +=logDetJacobianLevel(*(this->ThinLinks),0);
+
+      double end = usecond();
+      double time = (end - start)/ 1e3;
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: logDetJacobian took " << time << " ms" << std::endl;  
     }
     return ln_det;
   }
   void logDetJacobianForce(GaugeField &force)
   {
-    RealD ln_det = 0;
+    force =Zero();
+    GaugeField force_det(force.Grid());
+
     if (this->smearingLevels > 0)
     {
+      double start = usecond();
+
+      GaugeLinkField tmp_mu(force.Grid());
+
       for (int ismr = this->smearingLevels - 1; ismr > 0; --ismr) {
-	ln_det+= logDetJacobianForceLevel(this->get_smeared_conf(ismr-1),force,ismr);
+
+	// remove U in UdSdU...
+	for (int mu = 0; mu < Nd; mu++) {
+	  tmp_mu = adj(peekLorentz(this->get_smeared_conf(ismr), mu)) * peekLorentz(force, mu);
+	  pokeLorentz(force, tmp_mu, mu);
+	}
+	
+      	// Propagate existing force
+        force = this->AnalyticSmearedForce(force, this->get_smeared_conf(ismr - 1), ismr);
+
+	// Add back U in UdSdU...
+	for (int mu = 0; mu < Nd; mu++) {
+	  tmp_mu = peekLorentz(this->get_smeared_conf(ismr - 1), mu) * peekLorentz(force, mu);
+	  pokeLorentz(force, tmp_mu, mu);
+	}
+    	
+	// Get this levels determinant force
+	force_det = Zero();
+	logDetJacobianForceLevel(this->get_smeared_conf(ismr-1),force_det,ismr);
+
+	// Sum the contributions
+	force = force + force_det;
       }
-      ln_det +=logDetJacobianForeceLevel(*(this->ThinLinks),force,0);
+    
-    }
+      // remove U in UdSdU...
+      for (int mu = 0; mu < Nd; mu++) {
+	tmp_mu = adj(peekLorentz(this->get_smeared_conf(0), mu)) * peekLorentz(force, mu);
+	pokeLorentz(force, tmp_mu, mu);
+      }
+
+      force = this->AnalyticSmearedForce(force, *this->ThinLinks,0);
+
+      for (int mu = 0; mu < Nd; mu++) {
+	tmp_mu = peekLorentz(*this->ThinLinks, mu) * peekLorentz(force, mu);
+	pokeLorentz(force, tmp_mu, mu);
+      }
+
+      force_det = Zero();
+
+      logDetJacobianForceLevel(*this->ThinLinks,force_det,0);
+
+      force = force + force_det;
+
+      force=Ta(force); // Ta
+      
+      double end = usecond();
+      double time = (end - start)/ 1e3;
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: lnDetJacobianForce took " << time << " ms" << std::endl;  
+    }  // if smearingLevels = 0 do nothing
   }
 
 private:
-  // Member functions
   //====================================================================
   // Override base clas here to mask it
   virtual void fill_smearedSet(GaugeField &U)
@@ -574,6 +628,7 @@ private:
       GaugeField smeared_B(this->ThinLinks->Grid());
 
       previous_u = *this->ThinLinks;
+      double start = usecond();
       for (int smearLvl = 0; smearLvl < this->smearingLevels; ++smearLvl)
       {
         this->StoutSmearing->smear(smeared_A, previous_u);
@@ -586,8 +641,11 @@ private:
 
         // For debug purposes
         RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(previous_u);
-        std::cout << GridLogMessage << "[SmearedConfigurationMasked] Plaq: " << impl_plaq << std::endl;
+        std::cout << GridLogMessage << "[SmearedConfigurationMasked] smeared Plaq: " << impl_plaq << std::endl;
       }
+      double end = usecond();
+      double time = (end - start)/ 1e3;
+      std::cout << GridLogMessage << "GaugeConfigurationMasked: Link smearing took " << time << " ms" << std::endl;  
     }
   }
   //====================================================================
@@ -603,13 +661,11 @@ private:
     GaugeLinkField iQ(grid), e_iQ(grid);
     GaugeLinkField SigmaKPrime_mu(grid);
     GaugeLinkField GaugeKmu(grid), Cmu(grid);
-
+    
     this->StoutSmearing->BaseSmear(C, GaugeK);
     SigmaK = Zero();
     iLambda = Zero();
 
-    SigmaK = Zero();
-
     SigmaKPrimeA = SigmaKPrime;
     ApplyMask(SigmaKPrimeA,level);
     SigmaKPrimeB = SigmaKPrime - SigmaKPrimeA;
@@ -633,136 +689,28 @@ private:
     // propagate the rest of the force as identity map, just add back
     ////////////////////////////////////////////////////////////////////////////////////
     SigmaK = SigmaK+SigmaKPrimeB;
+
     return SigmaK;
   }
 
-  ////////////////////////////////////////
-  // INHERIT THESE
-  ////////////////////////////////////////
-  
-  /*! @brief Returns smeared configuration at level 'Level' */
-  /*
-  const GaugeField &get_smeared_conf(int Level) const
-  {
-    return SmearedSet[Level];
-  }
-  */
-  
-  // Duplicates code that is in GaugeConfiguration.h
-  // Should inherit or share.
-  //====================================================================
-  /*
-  void set_iLambda(GaugeLinkField& iLambda, GaugeLinkField& e_iQ,
-                   const GaugeLinkField& iQ, const GaugeLinkField& Sigmap,
-                   const GaugeLinkField& GaugeK) const 
-  {
-    GridBase* grid = iQ.Grid();
-    GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid);
-    GaugeLinkField unity(grid);
-    unity = 1.0;
-
-    LatticeComplex u(grid), w(grid);
-    LatticeComplex f0(grid), f1(grid), f2(grid);
-    LatticeComplex xi0(grid), xi1(grid), tmp(grid);
-    LatticeComplex u2(grid), w2(grid), cosw(grid);
-    LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(grid);
-    LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid);
-    LatticeComplex r22(grid), tr1(grid), tr2(grid);
-    LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid),
-      b22(grid);
-    LatticeComplex LatticeUnitComplex(grid);
-
-    LatticeUnitComplex = 1.0;
-
-    // Exponential
-    iQ2 = iQ * iQ;
-    iQ3 = iQ * iQ2;
-    StoutSmearing->set_uw(u, w, iQ2, iQ3);
-    StoutSmearing->set_fj(f0, f1, f2, u, w);
-    e_iQ = f0 * unity + timesMinusI(f1) * iQ - f2 * iQ2;
-
-    // Getting B1, B2, Gamma and Lambda
-    // simplify this part, reduntant calculations in set_fj
-    xi0 = StoutSmearing->func_xi0(w);
-    xi1 = StoutSmearing->func_xi1(w);
-    u2 = u * u;
-    w2 = w * w;
-    cosw = cos(w);
-
-    emiu = cos(u) - timesI(sin(u));
-    e2iu = cos(2.0 * u) + timesI(sin(2.0 * u));
-
-    r01 = (2.0 * u + timesI(2.0 * (u2 - w2))) * e2iu +
-      emiu * ((16.0 * u * cosw + 2.0 * u * (3.0 * u2 + w2) * xi0) +
-	      timesI(-8.0 * u2 * cosw + 2.0 * (9.0 * u2 + w2) * xi0));
-
-    r11 = (2.0 * LatticeUnitComplex + timesI(4.0 * u)) * e2iu +
-      emiu * ((-2.0 * cosw + (3.0 * u2 - w2) * xi0) +
-	      timesI((2.0 * u * cosw + 6.0 * u * xi0)));
-
-    r21 =
-      2.0 * timesI(e2iu) + emiu * (-3.0 * u * xi0 + timesI(cosw - 3.0 * xi0));
-
-    r02 = -2.0 * e2iu +
-      emiu * (-8.0 * u2 * xi0 +
-	      timesI(2.0 * u * (cosw + xi0 + 3.0 * u2 * xi1)));
-
-    r12 = emiu * (2.0 * u * xi0 + timesI(-cosw - xi0 + 3.0 * u2 * xi1));
-
-    r22 = emiu * (xi0 - timesI(3.0 * u * xi1));
-
-    fden = LatticeUnitComplex / (2.0 * (9.0 * u2 - w2) * (9.0 * u2 - w2));
-
-    b10 = 2.0 * u * r01 + (3.0 * u2 - w2) * r02 - (30.0 * u2 + 2.0 * w2) * f0;
-    b11 = 2.0 * u * r11 + (3.0 * u2 - w2) * r12 - (30.0 * u2 + 2.0 * w2) * f1;
-    b12 = 2.0 * u * r21 + (3.0 * u2 - w2) * r22 - (30.0 * u2 + 2.0 * w2) * f2;
-
-    b20 = r01 - (3.0 * u) * r02 - (24.0 * u) * f0;
-    b21 = r11 - (3.0 * u) * r12 - (24.0 * u) * f1;
-    b22 = r21 - (3.0 * u) * r22 - (24.0 * u) * f2;
-
-    b10 *= fden;
-    b11 *= fden;
-    b12 *= fden;
-    b20 *= fden;
-    b21 *= fden;
-    b22 *= fden;
-
-    B1 = b10 * unity + timesMinusI(b11) * iQ - b12 * iQ2;
-    B2 = b20 * unity + timesMinusI(b21) * iQ - b22 * iQ2;
-    USigmap = GaugeK * Sigmap;
-
-    tr1 = trace(USigmap * B1);
-    tr2 = trace(USigmap * B2);
-
-    GaugeLinkField QUS = iQ * USigmap;
-    GaugeLinkField USQ = USigmap * iQ;
-
-    GaugeLinkField iGamma = tr1 * iQ - timesI(tr2) * iQ2 +
-      timesI(f1) * USigmap + f2 * QUS + f2 * USQ;
-
-    iLambda = Ta(iGamma);
-  }
-  */
-  //====================================================================
 public:
-  //  GaugeField* ThinLinks; /* Pointer to the thin links configuration -- base class*/ 
+
-  ////////////////////////
-  // Derived class
-  ////////////////////////
   /* Standard constructor */
   SmearedConfigurationMasked(GridCartesian* _UGrid, unsigned int Nsmear, Smear_Stout<Gimpl>& Stout,bool domask=false)
     : SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout)
   {
     if(domask) assert(Nsmear%(2*Nd)==0); // Or multiply by 8??
 
+    // was resized in base class
+    assert(this->SmearedSet.size()==Nsmear);
+    
     GridRedBlackCartesian * UrbGrid;
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(_UGrid);
     LatticeComplex one(_UGrid); one = ComplexD(1.0,0.0);
     LatticeComplex tmp(_UGrid);
 
     for (unsigned int i = 0; i < this->smearingLevels; ++i) {
-      this->SmearedSet.push_back(*(new GaugeField(_UGrid)));
+
       masks.push_back(*(new LatticeLorentzComplex(_UGrid)));
       if (domask) {
 
@@ -787,97 +735,41 @@ public:
     }
     delete UrbGrid;
   }
-
-  //////////////////////////////////////////////////////////////
-  //Base functionality:
-  //////////////////////////////////////////////////////////////
   
-  /*! For just thin links */
+  virtual void smeared_force(GaugeField &SigmaTilde) 
-  //  SmearedConfigurationMasked()
-  //    : smearingLevels(0), StoutSmearing(nullptr), SmearedSet(), ThinLinks(NULL), UGrid(NULL), UrbGrid(NULL), masks() {}
-
-  // attach the smeared routines to the thin links U and fill the smeared set
-  /*
-  void set_Field(GaugeField &U)
   {
-    double start = usecond();
+    if (this->smearingLevels > 0)
-    fill_smearedSet(U);
-    double end = usecond();
-    double time = (end - start)/ 1e3;
-    std::cout << GridLogMessage << "Smearing in " << time << " ms" << std::endl;  
-  }
-  */
-  
-  //====================================================================
-  /*
-  void smeared_force(GaugeField &SigmaTilde) 
-  {
-    if (smearingLevels > 0)
     {
       double start = usecond();
       GaugeField force = SigmaTilde; // actually = U*SigmaTilde
       GaugeLinkField tmp_mu(SigmaTilde.Grid());
 
+      // Remove U from UdSdU
       for (int mu = 0; mu < Nd; mu++)
       {
         // to get just SigmaTilde
-        tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels - 1], mu)) * peekLorentz(force, mu);
+        tmp_mu = adj(peekLorentz(this->SmearedSet[this->smearingLevels - 1], mu)) * peekLorentz(force, mu);
         pokeLorentz(force, tmp_mu, mu);
       }
 
-      for (int ismr = smearingLevels - 1; ismr > 0; --ismr) {
+      for (int ismr = this->smearingLevels - 1; ismr > 0; --ismr) {
-        force = AnalyticSmearedForce(force, get_smeared_conf(ismr - 1),ismr);
+        force = this->AnalyticSmearedForce(force, this->get_smeared_conf(ismr - 1),ismr);
       }
       
-      force = AnalyticSmearedForce(force, *ThinLinks,0);
+      force = this->AnalyticSmearedForce(force, *this->ThinLinks,0);
 
+      // Add U to UdSdU
       for (int mu = 0; mu < Nd; mu++)
       {
-        tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
+        tmp_mu = peekLorentz(*this->ThinLinks, mu) * peekLorentz(force, mu);
         pokeLorentz(SigmaTilde, tmp_mu, mu);
       }
       double end = usecond();
       double time = (end - start)/ 1e3;
-      std::cout << GridLogMessage << "Smearing force in " << time << " ms" << std::endl;  
+      std::cout << GridLogMessage << " GaugeConfigurationMasked: Smeared Force chain rule took " << time << " ms" << std::endl;  
     }  // if smearingLevels = 0 do nothing
   }
-  */
-  //====================================================================
 
-  //  GaugeField& get_SmearedU() { return SmearedSet[smearingLevels - 1]; }
-  //  GaugeField& get_SmearedU(int n) { return this->SmearedSet[n]; }
-
-  /*
-  GaugeField &get_U(bool smeared = false)
-  {
-    // get the config, thin links by default
-    if (smeared)
-    {
-      if (smearingLevels)
-      {
-        RealD impl_plaq =
-	  WilsonLoops<Gimpl>::avgPlaquette(SmearedSet[smearingLevels - 1]);
-        std::cout << GridLogDebug << "getting Usmr Plaq: " << impl_plaq
-                  << std::endl;
-        return get_SmearedU();
-      }
-      else
-      {
-        RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(*ThinLinks);
-        std::cout << GridLogDebug << "getting Thin Plaq: " << impl_plaq
-                  << std::endl;
-        return *ThinLinks;
-      }
-    }
-    else
-    {
-      RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(*ThinLinks);
-      std::cout << GridLogDebug << "getting Thin Plaq: " << impl_plaq
-                << std::endl;
-      return *ThinLinks;
-    }
-  }
-  */
 };
 
 NAMESPACE_END(Grid);

Grid/qcd/smearing/JacobianAction.h

@@ -0,0 +1,91 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+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 */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+////////////////////////////////////////////////////////////////////////
+// Jacobian Action .. 
+////////////////////////////////////////////////////////////////////////
+template <class Gimpl>
+class JacobianAction : public Action<typename Gimpl::GaugeField> {
+public:  
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  SmearedConfigurationMasked<Gimpl> * smearer;
+  /////////////////////////// constructors
+  explicit JacobianAction(SmearedConfigurationMasked<Gimpl> * _smearer ) { smearer=_smearer;};
+
+  virtual std::string action_name() {return "JacobianAction";}
+
+  virtual std::string LogParameters(){
+    std::stringstream sstream;
+    sstream << GridLogMessage << "[JacobianAction] " << std::endl;
+    return sstream.str();
+  }
+
+  //////////////////////////////////
+  // Usual cases are not used
+  //////////////////////////////////
+  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){ assert(0);};
+  virtual RealD S(const GaugeField &U) { assert(0); }
+  virtual void deriv(const GaugeField &U, GaugeField &dSdU) { assert(0);  }
+
+  //////////////////////////////////
+  // Functions of smart configs only
+  //////////////////////////////////
+  virtual void refresh(ConfigurationBase<GaugeField> & U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
+  {
+    return;
+  }
+  virtual RealD S(ConfigurationBase<GaugeField>& U)
+  {
+    // det M = e^{ - ( - logDetM) }
+    assert( &U == smearer );
+    return -smearer->logDetJacobian();
+  }
+  virtual RealD Sinitial(ConfigurationBase<GaugeField>& U) 
+  {
+    return S(U);
+  }
+  virtual void deriv(ConfigurationBase<GaugeField>& U, GaugeField& dSdU)
+  {
+    assert( &U == smearer );
+    smearer->logDetJacobianForce(dSdU);
+  }
+
+private:
+ };
+
+NAMESPACE_END(Grid);
+

Grid/stencil/Stencil.h

@@ -451,7 +451,6 @@ public:
     else if ( this->fullDirichlet ) DslashLogDirichlet();
     else DslashLogFull();
     acceleratorCopySynchronise();
-    // Everyone agrees we are all done
     _grid->StencilBarrier(); 
   }
   ////////////////////////////////////////////////////////////////////////
@@ -540,6 +539,7 @@ public:
       compress.Point(point);
       HaloGatherDir(source,compress,point,face_idx);
     }
+    accelerator_barrier();
     face_table_computed=1;
     assert(u_comm_offset==_unified_buffer_size);
 

grid-config.in

@@ -60,7 +60,7 @@ while test $# -gt 0; do
     ;;
     
     --cxxflags)
-      echo @GRID_CXXFLAGS@
+      echo @GRID_CXXFLAGS@ -I@prefix@/include
     ;;
     
     --cxx)
@@ -72,11 +72,11 @@ while test $# -gt 0; do
     ;;
     
     --ldflags)
-      echo @GRID_LDFLAGS@
+      echo @GRID_LDFLAGS@ -L@prefix@/lib
     ;;
     
     --libs)
-      echo @GRID_LIBS@
+      echo @GRID_LIBS@ -lGrid
     ;;
     
     --summary)

systems/Lumi/config-command

@@ -0,0 +1,32 @@
+spack load c-lime
+spack load gmp
+spack load mpfr
+CLIME=`spack find --paths c-lime | grep c-lime| cut -c 15-`
+GMP=`spack find --paths gmp | grep gmp | cut -c 12-`
+MPFR=`spack find --paths mpfr | grep mpfr | cut -c 12-`
+echo clime $CLIME
+echo gmp $GMP
+echo mpfr $MPFR
+
+../../configure --enable-comms=mpi-auto \
+--with-lime=$CLIME \
+--enable-unified=no \
+--enable-shm=nvlink \
+--enable-tracing=timer \
+--enable-accelerator=hip \
+--enable-gen-simd-width=64 \
+--enable-simd=GPU \
+--disable-accelerator-cshift \
+--with-gmp=$OLCF_GMP_ROOT \
+--with-fftw=$FFTW_DIR/.. \
+--with-mpfr=/opt/cray/pe/gcc/mpfr/3.1.4/ \
+--disable-fermion-reps \
+--disable-gparity \
+CXX=hipcc MPICXX=mpicxx \
+CXXFLAGS="-fPIC -I{$ROCM_PATH}/include/ -std=c++14 -I${MPICH_DIR}/include -L/lib64 --amdgpu-target=gfx90a" \
+ LDFLAGS="-L/lib64 -L/opt/rocm-5.2.0/lib/ -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lamdhip64 --amdgpu-target=gfx90a "
+
+
+#--enable-simd=GPU-RRII \
+
+

systems/Lumi/sourceme.sh

@@ -0,0 +1 @@
+module load CrayEnv LUMI/22.12 partition/G  cray-fftw/3.3.10.1

systems/mac-arm/config-command-mpi

@@ -1,2 +1,4 @@
 BREW=/opt/local/
-MPICXX=mpicxx CXX=c++-12 ../../configure --enable-simd=GEN --enable-comms=mpi-auto --enable-unified=yes --prefix $HOME/QCD/GridInstallOpt --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW
+MPICXX=mpicxx CXX=c++-12 ../../configure --enable-simd=GEN --enable-comms=mpi-auto --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-fermion-reps --disable-gparity 
+
+

tests/forces/Test_fthmc.cc

@@ -0,0 +1,170 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_fthmc.cc
+
+    Copyright (C) 2022
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    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>
+#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
+#include <Grid/qcd/smearing/JacobianAction.h>
+
+using namespace std;
+using namespace Grid;
+
+
+template<class Gimpl>
+void ForceTest(Action<LatticeGaugeField> &action,SmearedConfigurationMasked<PeriodicGimplR> & smU,MomentumFilterBase<LatticeGaugeField> &Filter)
+{
+  LatticeGaugeField U = smU.get_U(false); // unsmeared config
+  GridBase *UGrid = U.Grid();
+
+  std::vector<int> seeds({1,2,3,5});
+  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
+
+  LatticeColourMatrix Pmu(UGrid); 
+  LatticeGaugeField P(UGrid); 
+  LatticeGaugeField UdSdU(UGrid); 
+
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+  
+  RealD eps=0.005;
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Refresh "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  
+  Gimpl::generate_momenta(P,sRNG,RNG4);
+  Filter.applyFilter(P);
+
+  action.refresh(smU,sRNG,RNG4);
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+
+  RealD S1 = action.S(smU);
+
+  Gimpl::update_field(P,U,eps);
+  smU.set_Field(U);
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Derivative "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  action.deriv(smU,UdSdU);
+  UdSdU = Ta(UdSdU);
+  Filter.applyFilter(UdSdU);
+
+  DumpSliceNorm("Force",UdSdU,Nd-1);
+  
+  Gimpl::update_field(P,U,eps);
+  smU.set_Field(U);
+
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
+  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  
+  RealD S2 = action.S(smU);
+
+  // Use the derivative
+  LatticeComplex dS(UGrid); dS = Zero();
+  for(int mu=0;mu<Nd;mu++){
+    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
+    Pmu= PeekIndex<LorentzIndex>(P,mu);
+    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
+  }
+  ComplexD dSpred    = sum(dS);
+  RealD diff =  S2-S1-dSpred.real();
+
+  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
+  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
+  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
+  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
+  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
+  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
+  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
+  //  assert(diff<1.0);
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::cout << std::setprecision(14);
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate mpi_layout  = GridDefaultMpi();
+  Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  Coordinate shm;
+  GlobalSharedMemory::GetShmDims(mpi_layout,shm);
+
+  const int Ls=12;
+  const int Nt = latt_size[3];
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+
+  ///////////////////// Gauge Field and Gauge Forces ////////////////////////////
+  LatticeGaugeField U(UGrid);
+
+#if  0
+  FieldMetaData header;
+  std::string file("./ckpoint_lat.2000");
+  NerscIO::readConfiguration(U,header,file);
+#else
+  std::vector<int> seeds({1,2,3,4,5,6,7,8});
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
+  SU<Nc>::HotConfiguration(RNG4,U);
+#endif
+
+  
+  RealD beta=6.0;
+  WilsonGaugeActionR  PlaqAction(beta);
+  IwasakiGaugeActionR RectAction(beta);
+
+
+  ////////////////////////////////////////////////
+  // Plaquette only FTHMC smearer
+  ////////////////////////////////////////////////
+  double rho = 0.1;
+  Smear_Stout<PeriodicGimplR> Smearer(rho);
+  SmearedConfigurationMasked<PeriodicGimplR> SmartConfig(UGrid,2*Nd,Smearer,true);
+
+  JacobianAction<PeriodicGimplR> Jacobian(&SmartConfig);
+  
+  ////////////////////////////////////////////////
+  // Run some tests
+  ////////////////////////////////////////////////
+  MomentumFilterNone<LatticeGaugeField> FilterNone;
+  SmartConfig.set_Field(U);
+  ForceTest<GimplTypesR>(PlaqAction,SmartConfig,FilterNone);
+  SmartConfig.set_Field(U);
+  ForceTest<GimplTypesR>(RectAction,SmartConfig,FilterNone);
+  SmartConfig.set_Field(U);
+  ForceTest<GimplTypesR>(Jacobian,SmartConfig,FilterNone);
+
+  Grid_finalize();
+}