Adding metric and the implicit steps

lib/qcd/action/Actions.h

@@ -50,6 +50,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 ////////////////////////////////////////////
 #include <Grid/qcd/action/gauge/GaugeImplementations.h>
 #include <Grid/qcd/utils/WilsonLoops.h>
+#include <Grid/qcd/utils/Metric.h>
 #include <Grid/qcd/utils/CovariantLaplacian.h>
 
 #include <Grid/qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions

lib/qcd/action/gauge/GaugeImplTypes.h

@@ -80,7 +80,7 @@ public:
 
   ///////////////////////////////////////////////////////////
   // Move these to another class
-  // HMC auxiliary functions
+  // HMC auxiliary functions 
   static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) {
     // specific for SU gauge fields
     LinkField Pmu(P._grid);

lib/qcd/hmc/HMC_GridModules.h

@@ -57,7 +57,7 @@ public:
   }    
 
   template <class ReaderClass>
-  GridModuleParameters(Reader<ReaderClass>& Reader) {
+  GridModuleParameters(Reader<ReaderClass>& Reader, std::string n = "LatticeGrid"):name(n) {
     read(Reader, name, *this);
     check();
   }
@@ -69,7 +69,7 @@ public:
     write(Writer, name, *this);
   }
 private:
-    std::string name = "LatticeGrid";
+    std::string name;
 };
 
 // Lower level class
@@ -94,7 +94,7 @@ class GridModule {
 ////////////////////////////////////
 // Classes for the user
 ////////////////////////////////////
-// Note: the space time grid must be out of the QCD namespace
+// Note: the space time grid should be out of the QCD namespace
 template< class vector_type>
 class GridFourDimModule : public GridModule {
  public:

lib/qcd/hmc/integrators/Integrator.h

@@ -77,7 +77,8 @@ class Integrator {
   double t_U;  // Track time passing on each level and for U and for P
   std::vector<double> t_P;  
 
-  MomentaField P;
+  //MomentaField P;
+  GeneralisedMomenta<FieldImplementation, TrivialMetric<MomentaField>> P;
   SmearingPolicy& Smearer;
   RepresentationPolicy Representations;
   IntegratorParameters Params;
@@ -86,7 +87,7 @@ class Integrator {
 
   void update_P(Field& U, int level, double ep) {
     t_P[level] += ep;
-    update_P(P, U, level, ep);
+    update_P(P.Mom, U, level, ep);
 
     std::cout << GridLogIntegrator << "[" << level << "] P "
               << " dt " << ep << " : t_P " << t_P[level] << std::endl;
@@ -131,51 +132,56 @@ class Integrator {
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
   }
 
-  void implicit_update_P(MomentaField& Mom, Field& U, int level, double ep) {
+  void implicit_update_P(Field& U, int level, double ep) {
+    t_P[level] += ep;
+
+    std::cout << GridLogIntegrator << "[" << level << "] P "
+              << " dt " << ep << " : t_P " << t_P[level] << std::endl;
     // Fundamental updates, include smearing
-    MomentaField Msum(Mom._grid);
+    MomentaField Msum(P.Mom._grid);
     Msum = zero;
     for (int a = 0; a < as[level].actions.size(); ++a) {
-      // Compute the force 
+      // Compute the force
       // We need to compute the derivative of the actions
       // only once
       Field force(U._grid);
-      conformable(U._grid, Mom._grid);
+      conformable(U._grid, P.Mom._grid);
       Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
       as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta
 
       std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
       if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
-      force = FieldImplementation::projectForce(force); // Ta for gauge fields
+      force = FieldImplementation::projectForce(force);  // Ta for gauge fields
-      Real force_abs = std::sqrt(norm2(force)/U._grid->gSites());
+      Real force_abs = std::sqrt(norm2(force) / U._grid->gSites());
-      std::cout << GridLogIntegrator << "Force average: " << force_abs << std::endl;
+      std::cout << GridLogIntegrator << "Force average: " << force_abs
+                << std::endl;
       Msum += force;
     }
 
-    MomentaField NewMom = Mom;
+    MomentaField NewMom = P.Mom;
-    MomentaField OldMom = Mom;
+    MomentaField OldMom = P.Mom;
     double threshold = 1e-6;
+    P.M.ImportGauge(U);
     // Here run recursively
-    do{
+    do {
-      MomentaField MomDer(Mom._grid);
+      MomentaField MomDer(P.Mom._grid);
+      MomentaField X(P.Mom._grid);
       OldMom = NewMom;
       // Compute the derivative of the kinetic term
       // with respect to the gauge field
-
+      P.DerivativeU(NewMom, MomDer);
-      // Laplacian.Mder(NewMom, MomDer);
+      NewMom = P.Mom - ep * (MomDer + Msum);
-      // NewMom = Mom - ep*(MomDer + Msum);
 
     } while (norm2(NewMom - OldMom) > threshold);
 
-    Mom = NewMom;
+    P.Mom = NewMom;
-
 
     // update the auxiliary fields momenta
   }
 
 
   void update_U(Field& U, double ep) {
-    update_U(P, U, ep);
+    update_U(P.Mom, U, ep);
 
     t_U += ep;
     int fl = levels - 1;
@@ -193,6 +199,27 @@ class Integrator {
     Representations.update(U);  // void functions if fundamental representation
   }
 
+  void implicit_update_U(Field&U, double ep){
+    t_U += ep;
+    int fl = levels - 1;
+    std::cout << GridLogIntegrator << "   " << "[" << fl << "] U " << " dt " << ep << " : t_U " << t_U << std::endl;
+
+    Real threshold = 1e-6;
+    P.M.ImportGauge(U);
+    MomentaField Mom1(P.Mom._grid);
+    MomentaField Mom2(P.Mom._grid);
+    Field OldU = U;
+    Field NewU = U;
+    P.DerivativeP(Mom1); // first term in the derivative 
+    do {
+      OldU = NewU; // some redundancy to be eliminated
+      P.DerivativeP(Mom2); // second term in the derivative, on the updated U
+      FieldImplementation::update_field(Mom1 + Mom2, NewU, ep);
+      P.M.ImportGauge(NewU);
+    } while (norm2(NewU - OldU) > threshold);
+  }
+
+
   virtual void step(Field& U, int level, int first, int last) = 0;
 
  public:
@@ -249,10 +276,11 @@ class Integrator {
 
   // Initialization of momenta and actions
   void refresh(Field& U, GridParallelRNG& pRNG) {
-    assert(P._grid == U._grid);
+    assert(P.Mom._grid == U._grid);
     std::cout << GridLogIntegrator << "Integrator refresh\n";
-    FieldImplementation::generate_momenta(P, pRNG);
+    //FieldImplementation::generate_momenta(P, pRNG);
-     
+    P.MomentaDistribution(pRNG);
+
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
     // of the Metropolis
@@ -297,7 +325,7 @@ class Integrator {
   // Calculate action
   RealD S(Field& U) {  // here also U not used
 
-    RealD H = - FieldImplementation::FieldSquareNorm(P); // - trace (P*P)
+    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom); // - trace (P*P)
     RealD Hterm;
     std::cout << GridLogMessage << "Momentum action H_p = " << H << "\n";
 

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -335,7 +335,7 @@ class ImplicitLeapFrog : public Integrator<FieldImplementation, SmearingPolicy,
       }
 
       if (level == fl) {  // lowest level
-        this->update_U(U, eps);
+        this->implicit_update_U(U, eps / 2.0);
       } else {  // recursive function call
         this->step(U, level + 1, first_step, last_step);
       }

lib/qcd/utils/CovariantLaplacian.h

@@ -33,6 +33,32 @@ directory
 namespace Grid {
 namespace QCD {
 
+struct LaplacianParams : Serializable {
+  GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianParams, 
+                                  RealD, lo, 
+                                  RealD, hi, 
+                                  int,   MaxIter, 
+                                  RealD, tolerance, 
+                                  int,   degree, 
+                                  int,   precision);
+  
+  // constructor 
+  LaplacianParams(RealD lo      = 0.0, 
+                  RealD hi      = 1.0, 
+                  int maxit     = 1000,
+                  RealD tol     = 1.0e-8, 
+                  int degree    = 10,
+                  int precision = 64)
+    : lo(lo),
+      hi(hi),
+      MaxIter(maxit),
+      tolerance(tol),
+      degree(degree),
+      precision(precision){};
+};
+
+
+
 ////////////////////////////////////////////////////////////
 // Laplacian operator L on adjoint fields
 //
@@ -48,12 +74,22 @@ namespace QCD {
 ////////////////////////////////////////////////////////////
 
 template <class Impl>
-class LaplacianAdjointField {
+class LaplacianAdjointField: public Metric<typename Impl::Field> {
+  OperatorFunction<typename Impl::Field> &Solver;
+  LaplacianParams param;
+  MultiShiftFunction PowerNegHalf;    
+
  public:
   INHERIT_GIMPL_TYPES(Impl);
-  
+
-  LaplacianAdjointField(GridBase* grid, const RealD k = 1.0) : 
+  LaplacianAdjointField(GridBase* grid, OperatorFunction<GaugeField>& S, LaplacianParams& p, const RealD k = 1.0)
-    U(Nd, grid), kappa(k){};
+      : U(Nd, grid), Solver(S), param(p), kappa(k){
+        AlgRemez remez(param.lo,param.hi,param.precision);
+        std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
+        remez.generateApprox(param.degree,1,2);
+        PowerNegHalf.Init(remez,param.tolerance,true);
+
+      };
 
   void ImportGauge(const GaugeField& _U) {
     for (int mu = 0; mu < Nd; mu++) {
@@ -61,41 +97,63 @@ class LaplacianAdjointField {
     }
   }
 
-  void Mdiag(const GaugeLinkField& in, GaugeLinkField& out) { assert(0); }
+  void M(const GaugeField& in, GaugeField& out) {
-
-  void Mdir(const GaugeLinkField& in, GaugeLinkField& out, int dir, int disp) {
-    assert(0);
-  }
-
-  void M(const GaugeLinkField& in, GaugeLinkField& out) {
     GaugeLinkField tmp(in._grid);
     GaugeLinkField tmp2(in._grid);
     GaugeLinkField sum(in._grid);
-    sum = zero;
+
-    for (int mu = 0; mu < Nd; mu++) {
+    for (int nu = 0; nu < Nd; nu++) {
-      tmp = U[mu] * Cshift(in, mu, +1) * adj(U[mu]);
+      sum = zero;
-      tmp2 = adj(U[mu]) * in * U[mu];
+      GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
-      sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in;
+      GaugeLinkField out_nu(out._grid);
+      for (int mu = 0; mu < Nd; mu++) {
+        tmp = U[mu] * Cshift(in_nu, mu, +1) * adj(U[mu]);
+        tmp2 = adj(U[mu]) * in_nu * U[mu];
+        sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_nu;
+      }
+      out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum;
+      PokeIndex<LorentzIndex>(out, out_nu, nu);
     }
-    out = (1.0 - kappa) * in - kappa / (double(4 * Nd)) * sum;
   }
 
-  void MDeriv(const GaugeLinkField& in, GaugeLinkField& out, bool dag){
+  void MDeriv(const GaugeField& in, GaugeField& der, bool dag) {
-    RealD factor = - kappa / (double(4 * Nd))
+    RealD factor = -kappa / (double(4 * Nd));
-    if (!dag)
+    for (int mu = 0; mu < Nd; mu++) {
-      out = factor * Cshift(in, mu, +1) * adj(U[mu]) + adj(U[mu]) * in;
+      GaugeLinkField in_mu = PeekIndex<LorentzIndex>(in, mu);
-    else
+      GaugeLinkField der_mu(der._grid);
-      out = factor * U[mu] * Cshift(in, mu, +1) + in * U[mu];
+      if (!dag)
+        der_mu =
+            factor * Cshift(in_mu, mu, +1) * adj(U[mu]) + adj(U[mu]) * in_mu;
+      else
+        der_mu = factor * U[mu] * Cshift(in_mu, mu, +1) + in_mu * U[mu];
+    }
   }
 
+  void Minv(const GaugeField& in, GaugeField& inverted){
+    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
+    Solver(HermOp, in, inverted);
+  }
+
+  void MInvSquareRoot(GaugeField& P){
+    // Takes a gaussian gauge field and multiplies by the metric
+    // need the rational approximation for the square root 
+    GaugeField Gp(P._grid);
+    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
+    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerNegHalf);
+    msCG(HermOp,P,Gp);
+    P = Gp; // now P has the correct distribution
+  }
+
+
+
  private:
   RealD kappa;
   std::vector<GaugeLinkField> U;
 };
 
 
-// This is just a debug tests
+// This is just for debuggin purposes
-// not meant to be used 
+// not meant to be used by the final users
 
 template <class Impl>
 class LaplacianAlgebraField {

lib/qcd/utils/Metric.h

@@ -0,0 +1,121 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/integrators/Integrator.h
+
+Copyright (C) 2015
+
+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 */
+//--------------------------------------------------------------------
+#ifndef METRIC_H
+#define METRIC_H
+
+namespace Grid{
+namespace QCD{
+
+template <typename Field> 
+class Metric{
+public:
+  virtual void ImportGauge(const Field&)   = 0;
+  virtual void M(const Field&, Field&)     = 0;
+  virtual void Minv(const Field&, Field&)  = 0;
+  virtual void MInvSquareRoot(Field&) = 0;
+};
+
+
+// Need a trivial operator
+template <typename Field>
+class TrivialMetric : public Metric<Field>{
+public:
+  virtual void ImportGauge(const Field&){};
+  virtual void M(const Field& in, Field& out){
+    out = in;
+  }
+  virtual void Minv(const Field& in, Field& out){
+    out = in;
+  }
+  virtual void MInvSquareRoot(Field& P){
+    // do nothing
+  }
+
+};
+
+///////////////////////////////
+// Generalised momenta
+///////////////////////////////
+
+template <typename Implementation, typename Metric>
+class GeneralisedMomenta{
+public:
+  typedef typename Implementation::Field MomentaField;  //for readability
+
+  Metric M;
+  MomentaField Mom;
+
+  GeneralisedMomenta(GridBase* grid): Mom(grid){}
+
+
+  void MomentaDistribution(GridParallelRNG& pRNG){
+    // Generate a distribution for
+    // 1/2 P^dag G P
+    // where G = M^-1
+
+    // Generate gaussian momenta
+    Implementation::generate_momenta(Mom, pRNG);
+    // Modify the distribution with the metric
+    M.MInvSquareRoot(Mom);
+  }
+
+  void Derivative(MomentaField& in, MomentaField& der){
+    // Compute the derivative of the kinetic term
+    // with respect to the gauge field
+    MomentaField MomDer(in._grid);
+    MomentaField X(in._grid);
+
+    M.Minv(in, X); // X = G in
+    M.MDeriv(X, MomDer, DaggerNo); // MomDer = dM/dU X
+    // MomDer is just the derivative
+    MomDer = adj(X)* MomDer;
+    // Traceless Antihermitian
+    // assuming we are in the algebra
+    der = Implementation::projectForce(MomDer);
+  }
+
+  void DerivativeP(MomentaField& der){
+    M.Minv(Mom, der);
+  }
+
+};
+
+
+
+
+
+
+
+
+}
+}
+
+
+#endif //METRIC_H

tests/solver/Test_laplacian.cc

@@ -50,8 +50,11 @@ int main (int argc, char ** argv)
 
   double Kappa = 0.9999;
 
+  std::cout << GridLogMessage << "Running with kappa: " << Kappa << std::endl;
+
   typedef SU<Nc>::LatticeAlgebraVector AVector;
   // Source and result in the algebra
+  // needed for the second test
   AVector src_vec(&Grid); random(pRNG, src_vec);
   AVector result_vec(&Grid); result_vec = zero;
   
@@ -59,16 +62,33 @@ int main (int argc, char ** argv)
   SU<Nc>::FundamentalLieAlgebraMatrix(src_vec, src);
   LatticeColourMatrix result(&Grid); result=zero;
 
-  LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid, Kappa);
-  Laplacian.ImportGauge(Umu);
 
-  HermitianLinearOperator<LaplacianAdjointField<PeriodicGimplR>,LatticeColourMatrix> HermOp(Laplacian);
+  // Generate a field of adjoint matrices
-  ConjugateGradient<LatticeColourMatrix> CG(1.0e-8,10000);
+  LatticeGaugeField src_f(&Grid);
+
+  // A matrix in the adjoint
+  LatticeColourMatrix src_mu(&Grid);
+  for (int mu = 0; mu < Nd; mu++) {
+    SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, src_mu);
+    PokeIndex<LorentzIndex>(src_f, src_mu, mu);
+  }
+  LatticeGaugeField result_f(&Grid);
+
+  // Definition of the Laplacian operator
+  ConjugateGradient<LatticeGaugeField> CG(1.0e-8,10000);
+  LaplacianParams LapPar(0.001, 1.0, 1000, 1e-8, 10, 64);
+  LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid, CG, LapPar, Kappa);
+  Laplacian.ImportGauge(Umu);
   std::cout << GridLogMessage << "Testing the Laplacian using the full matrix" <<std::endl;
-  CG(HermOp,src,result); // fastest
+  Laplacian.Minv(src_f, result_f);
   
 
+
+  Laplacian.MomentaDistribution(src_f);
+
+
   // Tests also the version using the algebra decomposition
+  /*
   LaplacianAlgebraField<PeriodicGimplR> LaplacianAlgebra(&Grid, Kappa);
   LaplacianAlgebra.ImportGauge(Umu);
 
@@ -82,7 +102,7 @@ int main (int argc, char ** argv)
 
   result2 -= result;
   std::cout << GridLogMessage << "Results difference " << norm2(result2) << std::endl;
-
+  */
 
   Grid_finalize();
 }