Getting rid of one more non-auto View, comms overlap in Laplace operator

Fixing Laplace flopcount Minor cleanup
Laplace benchmark added
2025-10-31 03:54:33 +00:00 · 2024-02-25 22:37:48 -05:00 · 2024-02-13 12:06:08 -05:00 · 2024-02-12 21:23:36 -05:00 · 2024-02-12 21:10:21 -05:00 · 2024-02-08 17:13:10 -05:00
19 changed files with 2598 additions and 104 deletions
--- a/Grid/algorithms/LinearOperator.h
+++ b/Grid/algorithms/LinearOperator.h
@@ -460,6 +460,53 @@ class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase<Fi
  }
 };

+template<class Matrix,class Field>
+class QuadLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat;
+public:
+  RealD a0,a1,a2;
+  QuadLinearOperator(Matrix &Mat): _Mat(Mat),a0(0.),a1(0.),a2(1.) {};
+  QuadLinearOperator(Matrix &Mat, RealD _a0,RealD _a1,RealD _a2): _Mat(Mat),a0(_a0),a1(_a1),a2(_a2) {};
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {
+    assert(0);
+    _Mat.Mdiag(in,out);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    assert(0);
+    _Mat.Mdir(in,out,dir,disp);
+  }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){
+    assert(0);
+    _Mat.MdirAll(in,out);
+  }
+  void HermOp (const Field &in, Field &out){
+//    _Mat.M(in,out);
+    Field tmp1(in.Grid());
+//    Linop.HermOpAndNorm(psi, mmp, d, b);
+    _Mat.M(in,tmp1);
+    _Mat.M(tmp1,out);
+    out *= a2;
+    axpy(out, a1, tmp1, out);
+    axpy(out, a0, in, out);
+//    d=real(innerProduct(psi,mmp));
+//    b=norm2(mmp);
+  }
+  void AdjOp     (const Field &in, Field &out){
+    assert(0);
+    _Mat.M(in,out);
+  }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    HermOp(in,out);
+    ComplexD dot= innerProduct(in,out); n1=real(dot);
+    n2=norm2(out);
+  }
+  void Op(const Field &in, Field &out){
+    assert(0);
+    _Mat.M(in,out);
+  }
+};
+
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 // Left  handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta  -->  ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta
 // Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta  -->  ( 1 - Moe Mee^-1 Meo Moo^-1) phi=eta ; psi = Moo^-1 phi
--- a/Grid/algorithms/approx/Forecast.h
+++ b/Grid/algorithms/approx/Forecast.h
@@ -36,11 +36,12 @@ NAMESPACE_BEGIN(Grid);
 // Abstract base class.
 // Takes a matrix (Mat), a source (phi), and a vector of Fields (chi)
 // and returns a forecasted solution to the system D*psi = phi (psi).
-template<class Matrix, class Field>
+// Changing to operator
+template<class LinearOperatorBase, class Field>
 class Forecast
 {
 public:
-  virtual Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& chi) = 0;
+  virtual Field operator()(LinearOperatorBase &Mat, const Field& phi, const std::vector<Field>& chi) = 0;
 };

 // Implementation of Brower et al.'s chronological inverter (arXiv:hep-lat/9509012),
@@ -54,13 +55,13 @@ public:
  Field operator()(Matrix &Mat, const Field& phi, const std::vector<Field>& prev_solns)
  {
    int degree = prev_solns.size();
+    std::cout << GridLogMessage << "ChronoForecast: degree= " << degree << std::endl;
    Field chi(phi); // forecasted solution

    // Trivial cases
    if(degree == 0){ chi = Zero(); return chi; }
    else if(degree == 1){ return prev_solns[0]; }

-    //    RealD dot;
    ComplexD xp;
    Field r(phi); // residual
    Field Mv(phi);
@@ -83,8 +84,9 @@ public:
    // Perform sparse matrix multiplication and construct rhs
    for(int i=0; i<degree; i++){
      b[i] = innerProduct(v[i],phi);
-      Mat.M(v[i],Mv);
-      Mat.Mdag(Mv,MdagMv[i]);
+//      Mat.M(v[i],Mv);
+//      Mat.Mdag(Mv,MdagMv[i]);
+      Mat.HermOp(v[i],MdagMv[i]);
      G[i][i] = innerProduct(v[i],MdagMv[i]);
    }

--- a/Grid/qcd/action/ActionCore.h
+++ b/Grid/qcd/action/ActionCore.h
@@ -67,6 +67,7 @@ NAMESPACE_CHECK(Scalar);
 #include <Grid/qcd/utils/Metric.h>
 NAMESPACE_CHECK(Metric);
 #include <Grid/qcd/utils/CovariantLaplacian.h>
+#include <Grid/qcd/utils/CovariantLaplacianRat.h>
 NAMESPACE_CHECK(CovariantLaplacian);


--- a/Grid/qcd/action/ActionParams.h
+++ b/Grid/qcd/action/ActionParams.h
@@ -65,6 +65,19 @@ struct WilsonImplParams {
  }
 };

+struct GaugeImplParams {
+//  bool overlapCommsCompute;
+//  AcceleratorVector<Real,Nd> twist_n_2pi_L;
+  AcceleratorVector<Complex,Nd> boundary_phases;
+  GaugeImplParams()  {
+    boundary_phases.resize(Nd, 1.0);
+//      twist_n_2pi_L.resize(Nd, 0.0);
+  };
+  GaugeImplParams(const AcceleratorVector<Complex,Nd> phi) : boundary_phases(phi) {
+//    twist_n_2pi_L.resize(Nd, 0.0);
+  }
+};
+
 struct StaggeredImplParams {
  Coordinate dirichlet; // Blocksize of dirichlet BCs
  int  partialDirichlet;
--- a/Grid/qcd/action/gauge/GaugeImplTypes.h
+++ b/Grid/qcd/action/gauge/GaugeImplTypes.h
@@ -32,7 +32,7 @@ directory

 NAMESPACE_BEGIN(Grid);

-#define CPS_MD_TIME
+#undef CPS_MD_TIME

 #ifdef CPS_MD_TIME
 #define HMC_MOMENTUM_DENOMINATOR (2.0)
--- a/Grid/qcd/action/gauge/WilsonGaugeAction.h
+++ b/Grid/qcd/action/gauge/WilsonGaugeAction.h
@@ -42,9 +42,13 @@ template <class Gimpl>
 class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
 public:  
  INHERIT_GIMPL_TYPES(Gimpl);
+  typedef GaugeImplParams ImplParams;
+  ImplParams Params;

  /////////////////////////// constructors
-  explicit WilsonGaugeAction(RealD beta_):beta(beta_){};
+  explicit WilsonGaugeAction(RealD beta_,
+		  const ImplParams &p = ImplParams()
+		  ):beta(beta_),Params(p){};

  virtual std::string action_name() {return "WilsonGaugeAction";}

@@ -56,14 +60,53 @@ public:

  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){};  // noop as no pseudoferms

+// Umu<->U maximally confusing
+  virtual void boundary(const GaugeField &Umu, GaugeField &Ub){
+    typedef typename Simd::scalar_type scalar_type;
+    assert(Params.boundary_phases.size() == Nd);
+    GridBase *GaugeGrid=Umu.Grid();
+    GaugeLinkField U(GaugeGrid);
+    GaugeLinkField tmp(GaugeGrid);
+
+    Lattice<iScalar<vInteger> > coor(GaugeGrid);
+    for (int mu = 0; mu < Nd; mu++) {
+	////////// boundary phase /////////////
+      auto pha = Params.boundary_phases[mu];
+      scalar_type phase( real(pha),imag(pha) );
+      std::cout<< GridLogIterative << "[WilsonGaugeAction] boundary "<<mu<<" "<<phase<< std::endl; 
+
+	int L   = GaugeGrid->GlobalDimensions()[mu];
+        int Lmu = L - 1;
+
+      LatticeCoordinate(coor, mu);
+
+      U = PeekIndex<LorentzIndex>(Umu, mu);
+      tmp = where(coor == Lmu, phase * U, U);
+      PokeIndex<LorentzIndex>(Ub, tmp, mu);
+//      PokeIndex<LorentzIndex>(Ub, U, mu);
+//      PokeIndex<LorentzIndex>(Umu, tmp, mu);
+
+    }
+  };
+
  virtual RealD S(const GaugeField &U) {
-    RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
-    RealD vol = U.Grid()->gSites();
+    GaugeField Ub(U.Grid());
+    this->boundary(U,Ub);
+    static RealD lastG=0.;
+    RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(Ub);
+    RealD vol = Ub.Grid()->gSites();
    RealD action = beta * (1.0 - plaq) * (Nd * (Nd - 1.0)) * vol * 0.5;
+    std::cout << GridLogMessage << "[WilsonGaugeAction] dH: " << action-lastG << std::endl;
+    RealD plaq_o = WilsonLoops<Gimpl>::avgPlaquette(U);
+    RealD action_o = beta * (1.0 - plaq_o) * (Nd * (Nd - 1.0)) * vol * 0.5;
+    std::cout << GridLogMessage << "[WilsonGaugeAction] U: " << action_o <<" Ub: "<< action  << std::endl;
+    lastG=action;
    return action;
  };

  virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
+    GaugeField Ub(U.Grid());
+    this->boundary(U,Ub);
    // not optimal implementation FIXME
    // extend Ta to include Lorentz indexes

@@ -73,10 +116,9 @@ public:
    GaugeLinkField dSdU_mu(U.Grid());
    for (int mu = 0; mu < Nd; mu++) {

-      Umu = PeekIndex<LorentzIndex>(U, mu);
-      
+      Umu = PeekIndex<LorentzIndex>(Ub, mu);
      // Staple in direction mu
-      WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
+      WilsonLoops<Gimpl>::Staple(dSdU_mu, Ub, mu);
      dSdU_mu = Ta(Umu * dSdU_mu) * factor;
      
      PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
--- a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
+++ b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h
@@ -178,7 +178,10 @@ NAMESPACE_BEGIN(Grid);
        // Use chronological inverter to forecast solutions across poles
        std::vector<FermionField> prev_solns;
        if(use_heatbath_forecasting){ prev_solns.reserve(param.degree); }
-        ChronoForecast<AbstractEOFAFermion<Impl>, FermionField> Forecast;
+	MdagMLinearOperator<AbstractEOFAFermion<Impl> ,FermionField> MdagML(Lop);
+	MdagMLinearOperator<AbstractEOFAFermion<Impl> ,FermionField> MdagMR(Rop);
+//        ChronoForecast<AbstractEOFAFermion<Impl>, FermionField> Forecast;
+	ChronoForecast<MdagMLinearOperator<AbstractEOFAFermion<Impl>, FermionField> , FermionField> Forecast;

        // \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta
        RealD N(PowerNegHalf.norm);
@@ -198,7 +201,7 @@ NAMESPACE_BEGIN(Grid);
          heatbathRefreshShiftCoefficients(0, -gamma_l);
          if(use_heatbath_forecasting){ // Forecast CG guess using solutions from previous poles
            Lop.Mdag(CG_src, Forecast_src);
-            CG_soln = Forecast(Lop, Forecast_src, prev_solns);
+            CG_soln = Forecast(MdagML, Forecast_src, prev_solns);
            SolverHBL(Lop, CG_src, CG_soln);
            prev_solns.push_back(CG_soln);
          } else {
@@ -225,7 +228,7 @@ NAMESPACE_BEGIN(Grid);
 	  heatbathRefreshShiftCoefficients(1, -gamma_l*PowerNegHalf.poles[k]);
          if(use_heatbath_forecasting){
            Rop.Mdag(CG_src, Forecast_src);
-            CG_soln = Forecast(Rop, Forecast_src, prev_solns);
+            CG_soln = Forecast(MdagMR, Forecast_src, prev_solns);
            SolverHBR(Rop, CG_src, CG_soln);
            prev_solns.push_back(CG_soln);
          } else {
--- a/Grid/qcd/action/scalar/ScalarImpl.h
+++ b/Grid/qcd/action/scalar/ScalarImpl.h
@@ -1,6 +1,6 @@
 #pragma once

-#define CPS_MD_TIME 
+#undef CPS_MD_TIME 

 #ifdef CPS_MD_TIME
 #define HMC_MOMENTUM_DENOMINATOR (2.0)
--- a/Grid/qcd/hmc/GenericHMCrunner.h
+++ b/Grid/qcd/hmc/GenericHMCrunner.h
@@ -121,12 +121,19 @@ public:

  template <class SmearingPolicy>
  void Run(SmearingPolicy &S) {
-    Runner(S);
+    TrivialMetric<typename Implementation::Field> Mtr;
+    Runner(S,Mtr);
+  }
+
+  template <class SmearingPolicy, class Metric>
+  void Run(SmearingPolicy &S, Metric &Mtr) {
+    Runner(S,Mtr);
  }

  void Run(){
    NoSmearing<Implementation> S;
-    Runner(S);
+    TrivialMetric<typename Implementation::Field> Mtr;
+    Runner(S,Mtr);
  }

  //Use the checkpointer to initialize the RNGs and the gauge field, writing the resulting gauge field into U.
@@ -176,15 +183,15 @@ public:
  //////////////////////////////////////////////////////////////////

 private:
-  template <class SmearingPolicy>
-  void Runner(SmearingPolicy &Smearing) {
+  template <class SmearingPolicy, class Metric>
+  void Runner(SmearingPolicy &Smearing, Metric &Mtr) {
    auto UGrid = Resources.GetCartesian();
    Field U(UGrid);

    initializeGaugeFieldAndRNGs(U);

    typedef IntegratorType<SmearingPolicy> TheIntegrator;
-    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing);
+    TheIntegrator MDynamics(UGrid, Parameters.MD, TheAction, Smearing,Mtr);

    // Sets the momentum filter
    MDynamics.setMomentumFilter(*(Resources.GetMomentumFilter()));
--- a/Grid/qcd/hmc/HMC.h
+++ b/Grid/qcd/hmc/HMC.h
@@ -55,6 +55,8 @@ struct HMCparameters: Serializable {
                                  Integer, NoMetropolisUntil,
 				  bool, PerformRandomShift, /* @brief Randomly shift the gauge configuration at the start of a trajectory */
                                  std::string, StartingType,
+				  Integer, SW,
+                                  RealD, Kappa,
                                  IntegratorParameters, MD)

  HMCparameters() {
@@ -110,6 +112,8 @@ private:
  IntegratorType &TheIntegrator;
  ObsListType Observables;

+  int traj_num;
+
  /////////////////////////////////////////////////////////
  // Metropolis step
  /////////////////////////////////////////////////////////
@@ -200,14 +204,14 @@ private:

    std::cout << GridLogMessage << "--------------------------------------------------\n";
    std::cout << GridLogMessage << " Molecular Dynamics evolution ";
-    TheIntegrator.integrate(U);
+    TheIntegrator.integrate(U,traj_num);
    std::cout << GridLogMessage << "--------------------------------------------------\n";

    //////////////////////////////////////////////////////////////////////////////////////////////////////
    // updated state action
    //////////////////////////////////////////////////////////////////////////////////////////////////////
    std::cout << GridLogMessage << "--------------------------------------------------\n";
-    std::cout << GridLogMessage << "Compute final action";
+    std::cout << GridLogMessage << "Compute final action" <<std::endl;
    RealD H1 = TheIntegrator.S(U);  
    std::cout << GridLogMessage << "--------------------------------------------------\n";

@@ -242,7 +246,7 @@ public:
  HybridMonteCarlo(HMCparameters _Pams, IntegratorType &_Int,
                   GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, 
                   ObsListType _Obs, Field &_U)
-    : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U) {}
+    : Params(_Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Observables(_Obs), Ucur(_U),traj_num(0) {}
  ~HybridMonteCarlo(){};

  void evolve(void) {
@@ -257,9 +261,10 @@ public:
    unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory;

    for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) {
+    

      std::cout << GridLogHMC << "-- # Trajectory = " << traj << "\n";
-
+      traj_num=traj;
      if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
      	std::cout << GridLogHMC << "-- Thermalization" << std::endl;
      }
--- a/Grid/qcd/hmc/integrators/Integrator.h
+++ b/Grid/qcd/hmc/integrators/Integrator.h
@@ -9,6 +9,7 @@ Copyright (C) 2015
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Guido Cossu <cossu@post.kek.jp>
+Author: Chulwoo Jung <chulwoo@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
@@ -33,6 +34,7 @@ directory
 #define INTEGRATOR_INCLUDED

 #include <memory>
+#include <Grid/parallelIO/NerscIO.h>

 NAMESPACE_BEGIN(Grid);

@@ -41,10 +43,19 @@ public:
  GRID_SERIALIZABLE_CLASS_MEMBERS(IntegratorParameters,
 				  std::string, name,      // name of the integrator
 				  unsigned int, MDsteps,  // number of outer steps
+				  RealD, RMHMCTol,
+                                  RealD, RMHMCCGTol,
+                                  RealD, lambda0,
+                                  RealD, lambda1,
+                                  RealD, lambda2,
 				  RealD, trajL)           // trajectory length

  IntegratorParameters(int MDsteps_ = 10, RealD trajL_ = 1.0)
  : MDsteps(MDsteps_),
+   lambda0(0.1931833275037836),
+   lambda1(0.1931833275037836),
+   lambda2(0.1931833275037836),
+   RMHMCTol(1e-8),RMHMCCGTol(1e-8),
    trajL(trajL_) {};

  template <class ReaderClass, typename std::enable_if<isReader<ReaderClass>::value, int >::type = 0 >
@@ -75,11 +86,14 @@ public:
  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 > P;
  SmearingPolicy& Smearer;
  RepresentationPolicy Representations;
  IntegratorParameters Params;

+  RealD Saux,Smom,Sg;
+
  //Filters allow the user to manipulate the conjugate momentum, for example to freeze links in DDHMC
  //It is applied whenever the momentum is updated / refreshed
  //The default filter does nothing
@@ -96,7 +110,16 @@ public:
  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;
+  }
+
+  void update_P2(Field& U, int level, double ep) 
+  {
+    t_P[level] += ep;
+    update_P2(P.Mom, U, level, ep);
+
    std::cout << GridLogIntegrator << "[" << level << "] P " << " dt " << ep << " : t_P " << t_P[level] << std::endl;
  }

@@ -119,62 +142,174 @@ public:
    }
  } update_P_hireps{};

- 
  void update_P(MomentaField& Mom, Field& U, int level, double ep) {
    // input U actually not used in the fundamental case
    // Fundamental updates, include smearing

    for (int a = 0; a < as[level].actions.size(); ++a) {
-
      double start_full = usecond();
      Field force(U.Grid());
      conformable(U.Grid(), Mom.Grid());

+      Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
      double start_force = usecond();
+      as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta

-      as[level].actions.at(a)->deriv_timer_start();
-      as[level].actions.at(a)->deriv(Smearer, force);  // deriv should NOT include Ta
-      as[level].actions.at(a)->deriv_timer_stop();
-
-      auto name = as[level].actions.at(a)->action_name();
-
+      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
      double end_force = usecond();
-      
-      MomFilter->applyFilter(force);
-
-      std::cout << GridLogIntegrator << " update_P : Level [" << level <<"]["<<a <<"] "<<name<<" dt "<<ep<<  std::endl;
-      
-      Real force_abs   = std::sqrt(norm2(force)/U.Grid()->gSites()); //average per-site norm.  nb. norm2(latt) = \sum_x norm2(latt[x]) 
-      Real impulse_abs = force_abs * ep * HMC_MOMENTUM_DENOMINATOR;    
-
-      Real force_max   = std::sqrt(maxLocalNorm2(force));
-      Real impulse_max = force_max * ep * HMC_MOMENTUM_DENOMINATOR;    
-
-      as[level].actions.at(a)->deriv_log(force_abs,force_max,impulse_abs,impulse_max);
-      
-      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] dt           : " << ep <<" "<<name<<std::endl;
-      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Force average: " << force_abs <<" "<<name<<std::endl;
-      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Force max    : " << force_max <<" "<<name<<std::endl;
-      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Fdt average  : " << impulse_abs <<" "<<name<<std::endl;
-      std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Fdt max      : " << impulse_max <<" "<<name<<std::endl;
-
+      Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites());
+      std::cout << GridLogIntegrator << "["<<level<<"]["<<a<<"] Force average: " << force_abs << std::endl;
      Mom -= force * ep* HMC_MOMENTUM_DENOMINATOR;; 
      double end_full = usecond();
      double time_full  = (end_full - start_full) / 1e3;
      double time_force = (end_force - start_force) / 1e3;
      std::cout << GridLogMessage << "["<<level<<"]["<<a<<"] P update elapsed time: " << time_full << " ms (force: " << time_force << " ms)"  << std::endl;
-
    }

    // Force from the other representations
    as[level].apply(update_P_hireps, Representations, Mom, U, ep);
+  }

+  void update_P2(MomentaField& Mom, Field& U, int level, double ep) {
+    // input U actually not used in the fundamental case
+    // Fundamental updates, include smearing
+
+    std::cout << GridLogIntegrator << "U before update_P2: " << std::sqrt(norm2(U)) << std::endl;
+    // Generalised momenta  
+    // Derivative of the kinetic term must be computed before
+    // Mom is the momenta and gets updated by the 
+    // actions derivatives
+    MomentaField MomDer(P.Mom.Grid());
+    P.M.ImportGauge(U);
+    P.DerivativeU(P.Mom, MomDer);
+    std::cout << GridLogIntegrator << "MomDer update_P2: " << std::sqrt(norm2(MomDer)) << std::endl;
+//    Mom -= MomDer * ep;
+    Mom -= MomDer * ep * HMC_MOMENTUM_DENOMINATOR;
+    std::cout << GridLogIntegrator << "Mom update_P2: " << std::sqrt(norm2(Mom)) << std::endl;
+
+    // Auxiliary fields
+    P.update_auxiliary_momenta(ep*0.5 );
+    P.AuxiliaryFieldsDerivative(MomDer);
+    std::cout << GridLogIntegrator << "MomDer(Aux) update_P2: " << std::sqrt(norm2(Mom)) << std::endl;
+//    Mom -= MomDer * ep;
+    Mom -= MomDer * ep * HMC_MOMENTUM_DENOMINATOR;
+    P.update_auxiliary_momenta(ep*0.5 );
+
+    for (int a = 0; a < as[level].actions.size(); ++a) {
+      double start_full = usecond();
+      Field force(U.Grid());
+      conformable(U.Grid(), Mom.Grid());
+
+      Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
+      double start_force = usecond();
+      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
+      double end_force = usecond();
+      Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites());
+      std::cout << GridLogIntegrator << "["<<level<<"]["<<a<<"] Force average: " << force_abs << std::endl;
+      Mom -= force * ep* HMC_MOMENTUM_DENOMINATOR;; 
+      double end_full = usecond();
+      double time_full  = (end_full - start_full) / 1e3;
+      double time_force = (end_force - start_force) / 1e3;
+      std::cout << GridLogMessage << "["<<level<<"]["<<a<<"] P update elapsed time: " << time_full << " ms (force: " << time_force << " ms)"  << std::endl;
+    }
+
+    // Force from the other representations
+    as[level].apply(update_P_hireps, Representations, Mom, U, ep);
+  }
+
+  void implicit_update_P(Field& U, int level, double ep, double ep1, bool intermediate = false) {
+    t_P[level] += ep;
+
+    double ep2= ep-ep1;
+
+    std::cout << GridLogIntegrator << "[" << level << "] P "
+              << " dt " << ep << " : t_P " << t_P[level] << std::endl;
+    std::cout << GridLogIntegrator << "U before implicit_update_P: " << std::sqrt(norm2(U)) << std::endl;
+    // Fundamental updates, include smearing
+    MomentaField Msum(P.Mom.Grid());
+    Msum = Zero();
+    for (int a = 0; a < as[level].actions.size(); ++a) {
+      // Compute the force terms for the lagrangian part
+      // We need to compute the derivative of the actions
+      // only once
+      Field force(U.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
+      Real force_abs = std::sqrt(norm2(force) / U.Grid()->gSites());
+      std::cout << GridLogIntegrator << "|Force| site average: " << force_abs
+                << std::endl;
+      Msum += force;
+    }
+
+    MomentaField NewMom = P.Mom;
+    MomentaField OldMom = P.Mom;
+    double threshold = Params.RMHMCTol;
+    P.M.ImportGauge(U);
+    MomentaField MomDer(P.Mom.Grid());
+    MomentaField MomDer1(P.Mom.Grid());
+    MomentaField AuxDer(P.Mom.Grid());
+    MomDer1 = Zero();
+    MomentaField diff(P.Mom.Grid());
+    double factor = 2.0;
+    if (intermediate){
+      P.DerivativeU(P.Mom, MomDer1);
+      factor = 1.0;
+    }
+//    std::cout << GridLogIntegrator << "MomDer1 implicit_update_P: " << std::sqrt(norm2(MomDer1)) << std::endl;
+
+    // Auxiliary fields
+    P.update_auxiliary_momenta(ep1);
+    P.AuxiliaryFieldsDerivative(AuxDer);
+    Msum += AuxDer;
+    
+
+    // Here run recursively
+    int counter = 1;
+    RealD RelativeError;
+    do {
+      std::cout << GridLogIntegrator << "UpdateP implicit step "<< counter << std::endl;
+
+      // Compute the derivative of the kinetic term
+      // with respect to the gauge field
+      P.DerivativeU(NewMom, MomDer);
+      Real force_abs = std::sqrt(norm2(MomDer) / U.Grid()->gSites());
+      std::cout << GridLogIntegrator << "|Force| laplacian site average: " << force_abs
+                << std::endl;
+
+//      NewMom = P.Mom - ep* 0.5 * HMC_MOMENTUM_DENOMINATOR * (2.0*Msum + factor*MomDer + MomDer1);// simplify
+      NewMom = P.Mom -  HMC_MOMENTUM_DENOMINATOR * (ep*Msum + ep1* factor*MomDer + ep2* MomDer1);// simplify
+      diff = NewMom - OldMom;
+      counter++;
+      RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewMom));
+      std::cout << GridLogIntegrator << "UpdateP RelativeError: " << RelativeError << std::endl;
+      OldMom = NewMom;
+    } while (RelativeError > threshold);
+
+    P.Mom = NewMom;
+    std::cout << GridLogIntegrator << "NewMom implicit_update_P: " << std::sqrt(norm2(NewMom)) << std::endl;
+
+    // update the auxiliary fields momenta    
+    P.update_auxiliary_momenta(ep2);
+  }
+
+  void implicit_update_P(Field& U, int level, double ep, bool intermediate = false) {
+      implicit_update_P( U, level, ep, ep*0.5, intermediate ); 
  }

  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;
@@ -183,12 +318,8 @@ public:
  
  void update_U(MomentaField& Mom, Field& U, double ep) 
  {
-    MomentaField MomFiltered(Mom.Grid());
-    MomFiltered = Mom;
-    MomFilter->applyFilter(MomFiltered);
-
    // exponential of Mom*U in the gauge fields case
-    FieldImplementation::update_field(MomFiltered, U, ep);
+    FieldImplementation::update_field(Mom, U, ep);

    // Update the smeared fields, can be implemented as observer
    Smearer.set_Field(U);
@@ -197,18 +328,74 @@ public:
    Representations.update(U);  // void functions if fundamental representation
  }

+  void implicit_update_U(Field&U, double ep, double ep1 ){
+    double ep2=ep-ep1;
+    t_U += ep;
+    int fl = levels - 1;
+    std::cout << GridLogIntegrator << "   " << "[" << fl << "] U " << " dt " << ep << " : t_U " << t_U << std::endl;
+    std::cout << GridLogIntegrator << "U before implicit_update_U: " << std::sqrt(norm2(U)) << std::endl;
+
+    MomentaField Mom1(P.Mom.Grid());
+    MomentaField Mom2(P.Mom.Grid());
+    RealD RelativeError;
+    Field diff(U.Grid());
+    Real threshold =  Params.RMHMCTol;
+    int counter = 1;
+    int MaxCounter = 100;
+
+    Field OldU = U;
+    Field NewU = U;
+
+    P.M.ImportGauge(U);
+    P.DerivativeP(Mom1); // first term in the derivative 
+    std::cout << GridLogIntegrator << "implicit_update_U: Mom1: " << std::sqrt(norm2(Mom1)) << std::endl;
+
+    P.update_auxiliary_fields(ep1);
+
+
+    MomentaField sum=Mom1;
+    do {
+      std::cout << GridLogIntegrator << "UpdateU implicit step "<< counter << std::endl;
+      
+      P.DerivativeP(Mom2); // second term in the derivative, on the updated U
+      std::cout << GridLogIntegrator << "implicit_update_U: Mom1: " << std::sqrt(norm2(Mom1)) << std::endl;
+      sum = (Mom1*ep1 + Mom2*ep2);
+
+      for (int mu = 0; mu < Nd; mu++) {
+        auto Umu = PeekIndex<LorentzIndex>(U, mu);
+        auto Pmu = PeekIndex<LorentzIndex>(sum, mu);
+        Umu = expMat(Pmu, 1, 12) * Umu;
+        PokeIndex<LorentzIndex>(NewU, ProjectOnGroup(Umu), mu);
+      }
+
+      diff = NewU - OldU;
+      RelativeError = std::sqrt(norm2(diff))/std::sqrt(norm2(NewU));
+      std::cout << GridLogIntegrator << "UpdateU RelativeError: " << RelativeError << std::endl;
+      
+      P.M.ImportGauge(NewU);
+      OldU = NewU; // some redundancy to be eliminated
+      counter++;
+    } while (RelativeError > threshold && counter < MaxCounter);
+
+    U = NewU;
+    std::cout << GridLogIntegrator << "NewU implicit_update_U: " << std::sqrt(norm2(U)) << std::endl;
+    P.update_auxiliary_fields(ep2);
+  }
+
+
  virtual void step(Field& U, int level, int first, int last) = 0;

 public:
  Integrator(GridBase* grid, IntegratorParameters Par,
             ActionSet<Field, RepresentationPolicy>& Aset,
-             SmearingPolicy& Sm)
+             SmearingPolicy& Sm, Metric<MomentaField>& M)
    : Params(Par),
      as(Aset),
-      P(grid),
+      P(grid, M),
      levels(Aset.size()),
      Smearer(Sm),
-      Representations(grid) 
+      Representations(grid),
+      Saux(0.),Smom(0.),Sg(0.)
  {
    t_P.resize(levels, 0.0);
    t_U = 0.0;
@@ -324,7 +511,8 @@ public:

  void reverse_momenta()
  {
-    P *= -1.0;
+    P.Mom *= -1.0;
+    P.AuxMom *= -1.0;
  }

  // to be used by the actionlevel class to iterate
@@ -343,11 +531,14 @@ public:
  // Initialization of momenta and actions
  void refresh(Field& U,  GridSerialRNG & sRNG, GridParallelRNG& pRNG) 
  {
-    assert(P.Grid() == U.Grid());
+    assert(P.Mom.Grid() == U.Grid());
    std::cout << GridLogIntegrator << "Integrator refresh" << std::endl;

    std::cout << GridLogIntegrator << "Generating momentum" << std::endl;
-    FieldImplementation::generate_momenta(P, sRNG, pRNG);
+//    FieldImplementation::generate_momenta(P.Mom, sRNG, pRNG);
+    P.M.ImportGauge(U);
+    P.MomentaDistribution(sRNG,pRNG);
+

    // Update the smeared fields, can be implemented as observer
    // necessary to keep the fields updated even after a reject
@@ -402,9 +593,22 @@ public:

    std::cout << GridLogIntegrator << "Integrator action\n";

-    RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+//    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+//    RealD Hterm;
+
+//    static RealD Saux=0.,Smom=0.,Sg=0.;
+
+    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+    std::cout << GridLogMessage << "S:FieldSquareNorm H_p = " << H << "\n";
+    std::cout << GridLogMessage << "S:dSField = " << H-Smom << "\n";
+    Smom=H;
+    P.M.ImportGauge(U);
+    RealD Hterm = - P.MomentaAction();
+    std::cout << GridLogMessage << "S:Momentum action H_p = " << Hterm << "\n";
+    std::cout << GridLogMessage << "S:dSMom = " << Hterm-Saux << "\n";
+    Saux=Hterm;
+    H = Hterm;

-    RealD Hterm;

    // Actions
    for (int level = 0; level < as.size(); ++level) {
@@ -446,9 +650,18 @@ public:

    std::cout << GridLogIntegrator << "Integrator initial action\n";

-    RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
-
-    RealD Hterm;
+//    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+//    RealD Hterm;
+    RealD H = - FieldImplementation::FieldSquareNorm(P.Mom)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
+    std::cout << GridLogMessage << "S:FieldSquareNorm H_p = " << H << "\n";
+    std::cout << GridLogMessage << "S:dSField = " << H-Smom << "\n";
+    Smom=H;
+    P.M.ImportGauge(U);
+    RealD Hterm = - P.MomentaAction();
+    std::cout << GridLogMessage << "S:Momentum action H_p = " << Hterm << "\n";
+    std::cout << GridLogMessage << "S:dSMom = " << Hterm-Saux << "\n";
+    Saux=Hterm;
+    H = Hterm;

    // Actions
    for (int level = 0; level < as.size(); ++level) {
@@ -471,7 +684,7 @@ public:
  }

  
-  void integrate(Field& U) 
+  void integrate(Field& U, int traj=-1 ) 
  {
    // reset the clocks
    t_U = 0;
@@ -483,6 +696,12 @@ public:
      int first_step = (stp == 0);
      int last_step = (stp == Params.MDsteps - 1);
      this->step(U, 0, first_step, last_step);
+      if (traj>=0){
+        std::string file("./config."+std::to_string(traj)+"_"+std::to_string(stp+1) );
+        int precision32 = 0;
+        int tworow      = 0;
+        NerscIO::writeConfiguration(U,file,tworow,precision32);
+      }
    }

    // Check the clocks all match on all levels
@@ -492,7 +711,6 @@ public:
    }

    FieldImplementation::Project(U);
-
    // and that we indeed got to the end of the trajectory
    assert(fabs(t_U - Params.trajL) < 1.0e-6);

--- a/Grid/qcd/hmc/integrators/Integrator_algorithm.h
+++ b/Grid/qcd/hmc/integrators/Integrator_algorithm.h
@@ -102,8 +102,8 @@ public:

  std::string integrator_name(){return "LeapFrog";}

-  LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
-    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm){};
+  LeapFrog(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm,M){};

  void step(Field& U, int level, int _first, int _last) {
    int fl = this->as.size() - 1;
@@ -140,14 +140,14 @@ template <class FieldImplementation_, class SmearingPolicy, class Representation
 class MinimumNorm2 : public Integrator<FieldImplementation_, SmearingPolicy, RepresentationPolicy> 
 {
 private:
-  const RealD lambda = 0.1931833275037836;
+//  const RealD lambda = 0.1931833275037836;

 public:
  typedef FieldImplementation_ FieldImplementation;
  INHERIT_FIELD_TYPES(FieldImplementation);

-  MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm)
-    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm){};
+  MinimumNorm2(GridBase* grid, IntegratorParameters Par, ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(grid, Par, Aset, Sm,M){};

  std::string integrator_name(){return "MininumNorm2";}

@@ -155,6 +155,11 @@ public:
    // level  : current level
    // fl     : final level
    // eps    : current step size
+    assert(level<3);
+    RealD lambda= this->Params.lambda0;
+    if (level>0) lambda= this->Params.lambda1;
+    if (level>1) lambda= this->Params.lambda2;
+    std::cout << GridLogMessage << "level: "<<level<< "lambda: "<<lambda<<std::endl;

    int fl = this->as.size() - 1;

@@ -210,9 +215,9 @@ public:
  // Looks like dH scales as dt^4. tested wilson/wilson 2 level.
  ForceGradient(GridBase* grid, IntegratorParameters Par,
                ActionSet<Field, RepresentationPolicy>& Aset,
-                SmearingPolicy& Sm)
+                SmearingPolicy& Sm, Metric<Field>& M)
    : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
-									    grid, Par, Aset, Sm){};
+									    grid, Par, Aset, Sm,M){};

  std::string integrator_name(){return "ForceGradient";}
  
@@ -275,6 +280,255 @@ public:
  }
 };

+////////////////////////////////
+// Riemannian Manifold HMC
+// Girolami et al
+////////////////////////////////
+
+
+
+// correct
+template <class FieldImplementation, class SmearingPolicy,
+          class RepresentationPolicy =
+              Representations<FundamentalRepresentation> >
+class ImplicitLeapFrog : public Integrator<FieldImplementation, SmearingPolicy,
+                                           RepresentationPolicy> {
+ public:
+  typedef ImplicitLeapFrog<FieldImplementation, SmearingPolicy, RepresentationPolicy>
+      Algorithm;
+  INHERIT_FIELD_TYPES(FieldImplementation);
+
+  // Riemannian manifold metric operator
+  // Hermitian operator Fisher
+
+  std::string integrator_name(){return "ImplicitLeapFrog";}
+
+  ImplicitLeapFrog(GridBase* grid, IntegratorParameters Par,
+           ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm, M){};
+
+  void step(Field& U, int level, int _first, int _last) {
+    int fl = this->as.size() - 1;
+    // level  : current level
+    // fl     : final level
+    // eps    : current step size
+
+    // Get current level step size
+    RealD eps = this->Params.trajL/this->Params.MDsteps;
+    for (int l = 0; l <= level; ++l) eps /= this->as[l].multiplier;
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+       this->implicit_update_P(U, level, eps / 2.0);
+      }
+
+      if (level == fl) {  // lowest level
+        this->implicit_update_U(U, eps,eps/2.);
+      } else {  // recursive function call
+        this->step(U, level + 1, first_step, last_step);
+      }
+
+      //int mm = last_step ? 1 : 2;
+      if (last_step){
+        this->update_P2(U, level, eps / 2.0);
+      } else {
+      this->implicit_update_P(U, level, eps, true);// works intermediate step
+      }
+    }
+  }
+};
+
+
+template <class FieldImplementation, class SmearingPolicy,
+          class RepresentationPolicy =
+              Representations<FundamentalRepresentation> >
+class ImplicitMinimumNorm2 : public Integrator<FieldImplementation, SmearingPolicy,
+                                       RepresentationPolicy> {
+ private:
+//  const RealD lambda = 0.1931833275037836;
+
+ public:
+  INHERIT_FIELD_TYPES(FieldImplementation);
+
+  ImplicitMinimumNorm2(GridBase* grid, IntegratorParameters Par,
+               ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm, M){};
+
+  std::string integrator_name(){return "ImplicitMininumNorm2";}
+
+  void step(Field& U, int level, int _first, int _last) {
+    // level  : current level
+    // fl     : final level
+    // eps    : current step size
+
+    int fl = this->as.size() - 1;
+//    assert(Params.lambda.size()>level);
+//    RealD lambda= Params.lambda[level];
+    assert(level<3);
+    RealD lambda= this->Params.lambda0;
+    if (level>0) lambda= this->Params.lambda1;
+    if (level>1) lambda= this->Params.lambda2;
+    std::cout << GridLogMessage << "level: "<<level<< "lambda: "<<lambda<<std::endl;
+
+  if(level<fl){
+
+    RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
+    for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
+
+    // Nesting:  2xupdate_U of size eps/2
+    // Next level is eps/2/multiplier
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {  // steps per step
+
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+        this->update_P(U, level, lambda * eps);
+      }
+
+        this->step(U, level + 1, first_step, 0);
+
+      this->update_P(U, level, (1.0 - 2.0 * lambda) * eps);
+
+        this->step(U, level + 1, 0, last_step);
+
+      int mm = (last_step) ? 1 : 2;
+      this->update_P(U, level, lambda * eps * mm);
+    }
+  } 
+  else 
+  { // last level
+    RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
+    for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
+
+    // Nesting:  2xupdate_U of size eps/2
+    // Next level is eps/2/multiplier
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {  // steps per step
+
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+        this->implicit_update_P(U, level, lambda * eps);
+      }
+
+      this->implicit_update_U(U, 0.5 * eps,lambda*eps);
+
+      this->implicit_update_P(U, level, (1.0 - 2.0 * lambda) * eps, true);
+
+      this->implicit_update_U(U, 0.5 * eps, (0.5-lambda)*eps);
+
+      if (last_step) {
+        this->update_P2(U, level, eps * lambda);
+      } else {
+        this->implicit_update_P(U, level, lambda * eps*2.0, true);
+      }
+    }
+  }
+
+  }
+};
+
+template <class FieldImplementation, class SmearingPolicy,
+          class RepresentationPolicy =
+              Representations<FundamentalRepresentation> >
+class ImplicitCampostrini : public Integrator<FieldImplementation, SmearingPolicy,
+                                       RepresentationPolicy> {
+ private:
+//  const RealD lambda = 0.1931833275037836;
+
+ public:
+  INHERIT_FIELD_TYPES(FieldImplementation);
+
+  ImplicitCampostrini(GridBase* grid, IntegratorParameters Par,
+               ActionSet<Field, RepresentationPolicy>& Aset, SmearingPolicy& Sm, Metric<Field>& M)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm, M){};
+
+  std::string integrator_name(){return "ImplicitCampostrini";}
+
+  void step(Field& U, int level, int _first, int _last) {
+    // level  : current level
+    // fl     : final level
+    // eps    : current step size
+
+    int fl = this->as.size() - 1;
+//    assert(Params.lambda.size()>level);
+//    RealD lambda= Params.lambda[level];
+    assert(level<3);
+    RealD lambda= this->Params.lambda0;
+    if (level>0) lambda= this->Params.lambda1;
+    if (level>1) lambda= this->Params.lambda2;
+    std::cout << GridLogMessage << "level: "<<level<< "lambda: "<<lambda<<std::endl;
+    
+    RealD sigma=pow(2.0,1./3.);
+
+  if(level<fl){
+//Still Omelyan. Needs to change step() to accept variable stepsize
+    RealD eps = this->Params.trajL/this->Params.MDsteps * 2.0;
+    for (int l = 0; l <= level; ++l) eps /= 2.0 * this->as[l].multiplier;
+
+    // Nesting:  2xupdate_U of size eps/2
+    // Next level is eps/2/multiplier
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {  // steps per step
+
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+
+      if (first_step) {  // initial half step
+        this->update_P(U, level, lambda * eps);
+      }
+
+        this->step(U, level + 1, first_step, 0);
+
+      this->update_P(U, level, (1.0 - 2.0 * lambda) * eps);
+
+        this->step(U, level + 1, 0, last_step);
+
+      int mm = (last_step) ? 1 : 2;
+      this->update_P(U, level, lambda * eps * mm);
+    }
+  } 
+  else 
+  { // last level
+    RealD dt = this->Params.trajL/this->Params.MDsteps * 2.0;
+    for (int l = 0; l <= level; ++l) dt /= 2.0 * this->as[l].multiplier;
+
+    RealD epsilon = dt/(2.0 - sigma);
+
+    int multiplier = this->as[level].multiplier;
+    for (int e = 0; e < multiplier; ++e) {  // steps per step
+
+      int first_step = _first && (e == 0);
+      int last_step = _last && (e == multiplier - 1);
+      // initial half step
+      if (first_step) {  this->implicit_update_P(U, level, epsilon*0.5); }
+      this->implicit_update_U(U, epsilon,epsilon*0.5);
+      this->implicit_update_P(U, level, (1.0 - sigma) * epsilon *0.5, epsilon*0.5, true);
+      this->implicit_update_U(U, -epsilon*sigma, -epsilon*sigma*0.5);
+      this->implicit_update_P(U, level, (1.0 - sigma) * epsilon *0.5, -epsilon*sigma*0.5, true);
+      this->implicit_update_U(U, epsilon,epsilon*0.5);
+      if (last_step) { this->update_P2(U, level, epsilon*0.5 ); } 
+      else
+      this->implicit_update_P(U, level, epsilon,epsilon*0.5);
+    }
+  }
+
+  }
+};
+
 NAMESPACE_END(Grid);

 #endif  // INTEGRATOR_INCLUDED
--- a/Grid/qcd/utils/CovariantLaplacian.h
+++ b/Grid/qcd/utils/CovariantLaplacian.h
@@ -54,7 +54,361 @@ struct LaplacianParams : Serializable {
      precision(precision){};
 };

+#define LEG_LOAD(Dir)						 \
+  SE = st.GetEntry(ptype, Dir, ss);				 \
+  if (SE->_is_local ) {						 \
+    int perm= SE->_permute;					 \
+    chi = coalescedReadPermute(in[SE->_offset],ptype,perm,lane); \
+  } else {							 \
+    chi = coalescedRead(buf[SE->_offset],lane);			 \
+  }								 \
+  acceleratorSynchronise();

+const std::vector<int> directions4D   ({Xdir,Ydir,Zdir,Tdir,Xdir,Ydir,Zdir,Tdir});
+const std::vector<int> displacements4D({1,1,1,1,-1,-1,-1,-1});
+
+template<class Gimpl,class Field> class CovariantAdjointLaplacianStencil : public SparseMatrixBase<Field>
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+//  RealD kappa;
+
+  typedef typename Field::vector_object siteObject;
+
+  template <typename vtype> using iImplDoubledGaugeField = iVector<iScalar<iMatrix<vtype, Nc> >, Nds>;
+  typedef iImplDoubledGaugeField<Simd> SiteDoubledGaugeField;
+  typedef Lattice<SiteDoubledGaugeField> DoubledGaugeField;
+  typedef CartesianStencil<siteObject, siteObject, DefaultImplParams> StencilImpl;
+
+  GridBase *grid;
+  StencilImpl Stencil;
+  SimpleCompressor<siteObject> Compressor;
+  DoubledGaugeField Uds;
+
+  CovariantAdjointLaplacianStencil( GridBase *_grid)
+    : grid(_grid),
+      Stencil    (grid,8,Even,directions4D,displacements4D),
+      Uds(grid){}
+
+  CovariantAdjointLaplacianStencil(GaugeField &Umu)
+    :
+      grid(Umu.Grid()),
+      Stencil    (grid,8,Even,directions4D,displacements4D),
+      Uds(grid)
+  { GaugeImport(Umu); }
+
+  void GaugeImport (const GaugeField &Umu)
+  {
+    assert(grid == Umu.Grid());
+    for (int mu = 0; mu < Nd; mu++) {
+      auto U = PeekIndex<LorentzIndex>(Umu, mu);
+      PokeIndex<LorentzIndex>(Uds, U, mu );
+      U = adj(Cshift(U, mu, -1));
+      PokeIndex<LorentzIndex>(Uds, U, mu + 4);
+    }
+  };
+  
+  virtual GridBase *Grid(void) { return grid; };
+//broken
+#if 0
+  virtual void  MDeriv(const Field &_left, Field &_right,Field &_der, int mu)
+  {
+    ///////////////////////////////////////////////
+    // Halo exchange for this geometry of stencil
+    ///////////////////////////////////////////////
+    Stencil.HaloExchange(_lef, Compressor);
+
+    ///////////////////////////////////
+    // Arithmetic expressions
+    ///////////////////////////////////
+    autoView( st     , Stencil    , AcceleratorRead);
+    auto buf = st.CommBuf();
+
+    autoView( in     , _left    , AcceleratorRead);
+    autoView( right    , _right   , AcceleratorRead);
+    autoView( der    , _der   , AcceleratorWrite);
+    autoView( U     , Uds    , AcceleratorRead);
+
+    typedef typename Field::vector_object        vobj;
+    typedef decltype(coalescedRead(left[0]))    calcObj;
+    typedef decltype(coalescedRead(U[0](0))) calcLink;
+
+    const int      Nsimd = vobj::Nsimd();
+    const uint64_t NN = grid->oSites();
+
+    accelerator_for( ss, NN, Nsimd, {
+
+	StencilEntry *SE;
+	
+	const int lane=acceleratorSIMTlane(Nsimd);
+
+	calcObj chi;
+	calcObj phi;
+	calcObj res;
+	calcObj Uchi;
+	calcObj Utmp;
+	calcObj Utmp2;
+	calcLink UU;
+	calcLink Udag;
+	int ptype;
+
+	res                 = coalescedRead(def[ss]);
+	phi                 = coalescedRead(right[ss]);
+
+#define LEG_LOAD_MULT_LINK(leg,polarisation)			\
+	UU = coalescedRead(U[ss](polarisation));	\
+	Udag = adj(UU);					\
+	LEG_LOAD(leg);					\
+	mult(&Utmp(), &UU, &chi());			\
+	Utmp2 = adj(Utmp);				\
+	mult(&Utmp(), &UU, &Utmp2());			\
+	Utmp2 = adj(Utmp);				\
+	mult(&Uchi(), &phi(), &Utmp2());			\
+	res = res + Uchi;
+	
+	LEG_LOAD_MULT_LINK(0,Xp);
+	LEG_LOAD_MULT_LINK(1,Yp);
+	LEG_LOAD_MULT_LINK(2,Zp);
+	LEG_LOAD_MULT_LINK(3,Tp);
+
+	coalescedWrite(der[ss], res,lane);
+    });
+
+  };
+#endif
+
+  virtual void  Morig(const Field &_in, Field &_out)
+  {
+    ///////////////////////////////////////////////
+    // Halo exchange for this geometry of stencil
+    ///////////////////////////////////////////////
+    Stencil.HaloExchange(_in, Compressor);
+
+    ///////////////////////////////////
+    // Arithmetic expressions
+    ///////////////////////////////////
+//    auto st = Stencil.View(AcceleratorRead);
+    autoView( st     , Stencil    , AcceleratorRead);
+    auto buf = st.CommBuf();
+
+    autoView( in     , _in    , AcceleratorRead);
+    autoView( out    , _out   , AcceleratorWrite);
+    autoView( U     , Uds    , AcceleratorRead);
+
+    typedef typename Field::vector_object        vobj;
+    typedef decltype(coalescedRead(in[0]))    calcObj;
+    typedef decltype(coalescedRead(U[0](0))) calcLink;
+
+    const int      Nsimd = vobj::Nsimd();
+    const uint64_t NN = grid->oSites();
+
+    accelerator_for( ss, NN, Nsimd, {
+
+	StencilEntry *SE;
+	
+	const int lane=acceleratorSIMTlane(Nsimd);
+
+	calcObj chi;
+	calcObj res;
+	calcObj Uchi;
+	calcObj Utmp;
+	calcObj Utmp2;
+	calcLink UU;
+	calcLink Udag;
+	int ptype;
+
+	res                 = coalescedRead(in[ss])*(-8.0);
+
+#define LEG_LOAD_MULT(leg,polarisation)			\
+	UU = coalescedRead(U[ss](polarisation));	\
+	Udag = adj(UU);					\
+	LEG_LOAD(leg);					\
+	mult(&Utmp(), &UU, &chi());			\
+	Utmp2 = adj(Utmp);				\
+	mult(&Utmp(), &UU, &Utmp2());			\
+	Uchi = adj(Utmp);				\
+	res = res + Uchi;
+	
+	LEG_LOAD_MULT(0,Xp);
+	LEG_LOAD_MULT(1,Yp);
+	LEG_LOAD_MULT(2,Zp);
+	LEG_LOAD_MULT(3,Tp);
+	LEG_LOAD_MULT(4,Xm);
+	LEG_LOAD_MULT(5,Ym);
+	LEG_LOAD_MULT(6,Zm);
+	LEG_LOAD_MULT(7,Tm);
+
+	coalescedWrite(out[ss], res,lane);
+    });
+
+  };
+  virtual void  Mnew (const Field &_in, Field &_out)
+  {
+    ///////////////////////////////////////////////
+    // Halo exchange for this geometry of stencil
+    ///////////////////////////////////////////////
+//    Stencil.HaloExchange(_in, Compressor);
+      std::vector<std::vector<CommsRequest_t> > requests;
+      Stencil.Prepare();
+  {
+    GRID_TRACE("Laplace Gather");
+    Stencil.HaloGather(_in,Compressor);
+  }
+
+  tracePush("Laplace Communication");
+  Stencil.CommunicateBegin(requests);
+  {
+    GRID_TRACE("MergeSHM");
+    Stencil.CommsMergeSHM(Compressor);
+  }
+    
+
+    ///////////////////////////////////
+    // Arithmetic expressions
+    ///////////////////////////////////
+//    auto st = Stencil.View(AcceleratorRead);
+    autoView( st     , Stencil    , AcceleratorRead);
+    auto buf = st.CommBuf();
+
+    autoView( in     , _in    , AcceleratorRead);
+    autoView( out    , _out   , AcceleratorWrite);
+    autoView( U     , Uds    , AcceleratorRead);
+
+    typedef typename Field::vector_object        vobj;
+    typedef decltype(coalescedRead(in[0]))    calcObj;
+    typedef decltype(coalescedRead(U[0](0))) calcLink;
+
+    const int      Nsimd = vobj::Nsimd();
+    const uint64_t NN = grid->oSites();
+
+    accelerator_for( ss, NN, Nsimd, {
+
+	StencilEntry *SE;
+	
+	const int lane=acceleratorSIMTlane(Nsimd);
+
+	calcObj chi;
+	calcObj res;
+	calcObj Uchi;
+	calcObj Utmp;
+	calcObj Utmp2;
+	calcLink UU;
+	calcLink Udag;
+	int ptype;
+
+	res                 = coalescedRead(in[ss])*(-8.0);
+
+
+        SE = st.GetEntry(ptype, 0, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(0,Xp);
+	}
+        SE = st.GetEntry(ptype, 1, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(1,Yp);
+	}
+        SE = st.GetEntry(ptype, 2, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(2,Zp);
+	}
+        SE = st.GetEntry(ptype, 3, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(3,Tp);
+	}
+        SE = st.GetEntry(ptype, 4, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(4,Xm);
+	}
+        SE = st.GetEntry(ptype, 5, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(5,Ym);
+	}
+        SE = st.GetEntry(ptype, 6, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(6,Zm);
+	}
+        SE = st.GetEntry(ptype, 7, ss);				 
+        if (SE->_is_local ) {
+	LEG_LOAD_MULT(7,Tm);
+	}
+
+	coalescedWrite(out[ss], res,lane);
+    });
+
+    Stencil.CommunicateComplete(requests);
+  tracePop("Communication");
+
+  {
+    GRID_TRACE("Merge");
+    Stencil.CommsMerge(Compressor);
+  }
+
+
+    accelerator_for( ss, NN, Nsimd, {
+
+	StencilEntry *SE;
+	
+	const int lane=acceleratorSIMTlane(Nsimd);
+
+	calcObj chi;
+	calcObj res;
+	calcObj Uchi;
+	calcObj Utmp;
+	calcObj Utmp2;
+	calcLink UU;
+	calcLink Udag;
+	int ptype;
+
+//	res                 = coalescedRead(in[ss])*(-8.0);
+	res                 = coalescedRead(out[ss]);
+
+        SE = st.GetEntry(ptype, 0, ss);				 
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(0,Xp);
+	}
+        SE = st.GetEntry(ptype, 1, ss);				 
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(1,Yp);
+	}
+        SE = st.GetEntry(ptype, 2, ss);				 
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(2,Zp);
+	}
+        SE = st.GetEntry(ptype, 3, ss);
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(3,Tp);
+	}
+        SE = st.GetEntry(ptype, 4, ss);
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(4,Xm);
+	}
+        SE = st.GetEntry(ptype, 5, ss);
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(5,Ym);
+	}
+        SE = st.GetEntry(ptype, 6, ss);
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(6,Zm);
+	}
+        SE = st.GetEntry(ptype, 7, ss);
+        if ((SE->_is_local )==0){
+	LEG_LOAD_MULT(7,Tm);
+	}
+
+	coalescedWrite(out[ss], res,lane);
+    });
+  };
+
+  virtual void  M(const Field &in, Field &out) {Mnew(in,out);};
+  virtual void  Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
+  virtual  void Mdiag    (const Field &in, Field &out)                  {assert(0);}; // Unimplemented need only for multigrid
+  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
+  virtual  void MdirAll  (const Field &in, std::vector<Field> &out)     {assert(0);}; // Unimplemented need only for multigrid
+};
+
+#undef LEG_LOAD_MULT
+#undef LEG_LOAD_MULT_LINK
+#undef LEG_LOAD

 ////////////////////////////////////////////////////////////
 // Laplacian operator L on adjoint fields
@@ -76,29 +430,40 @@ class LaplacianAdjointField: public Metric<typename Impl::Field> {
  LaplacianParams param;
  MultiShiftFunction PowerHalf;    
  MultiShiftFunction PowerInvHalf;    
+//template<class Gimpl,class Field> class CovariantAdjointLaplacianStencil : public SparseMatrixBase<Field>
+  CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField> LapStencil;

 public:
  INHERIT_GIMPL_TYPES(Impl);

-  LaplacianAdjointField(GridBase* grid, OperatorFunction<GaugeField>& S, LaplacianParams& p, const RealD k = 1.0)
-    : U(Nd, grid), Solver(S), param(p), kappa(k){
+  LaplacianAdjointField(GridBase* grid, OperatorFunction<GaugeField>& S, LaplacianParams& p, const RealD k = 1.0, bool if_remez=true)
+    : U(Nd, grid), Solver(S), param(p), kappa(k)
+	,LapStencil(grid){
    AlgRemez remez(param.lo,param.hi,param.precision);
    std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
+    if(if_remez){
    remez.generateApprox(param.degree,1,2);
    PowerHalf.Init(remez,param.tolerance,false);
    PowerInvHalf.Init(remez,param.tolerance,true);
+    }
+    this->triv=0;
        

  };
-
+  LaplacianAdjointField(){this->triv=0; printf("triv=%d\n",this->Trivial());}
  void Mdir(const GaugeField&, GaugeField&, int, int){ assert(0);}
  void MdirAll(const GaugeField&, std::vector<GaugeField> &){ assert(0);}
  void Mdiag(const GaugeField&, GaugeField&){ assert(0);}

  void ImportGauge(const GaugeField& _U) {
+    RealD total=0.;
    for (int mu = 0; mu < Nd; mu++) {
      U[mu] = PeekIndex<LorentzIndex>(_U, mu);
+      total += norm2(U[mu]);
    }
+    LapStencil.GaugeImport (_U);
+
+    std::cout << GridLogDebug <<"ImportGauge:norm2(U _U) = "<<total<<std::endl;
  }

  void M(const GaugeField& in, GaugeField& out) {
@@ -106,10 +471,12 @@ public:
    // test
    //GaugeField herm = in + adj(in);
    //std::cout << "AHermiticity: " << norm2(herm) << std::endl;
+//    std::cout << GridLogDebug <<"M:Kappa = "<<kappa<<std::endl;

+    GaugeLinkField sum(in.Grid());
+#if 0
    GaugeLinkField tmp(in.Grid());
    GaugeLinkField tmp2(in.Grid());
-    GaugeLinkField sum(in.Grid());

    for (int nu = 0; nu < Nd; nu++) {
      sum = Zero();
@@ -123,10 +490,22 @@ public:
      out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum;
      PokeIndex<LorentzIndex>(out, out_nu, nu);
    }
+#else
+    for (int nu = 0; nu < Nd; nu++) {
+      GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
+      GaugeLinkField out_nu(out.Grid());
+      LapStencil.M(in_nu,sum);
+      out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum;
+      PokeIndex<LorentzIndex>(out, out_nu, nu);
+    }
+#endif
+//    std::cout << GridLogDebug <<"M:norm2(out) = "<<norm2(out)<<std::endl;
  }

+
  void MDeriv(const GaugeField& in, GaugeField& der) {
    // in is anti-hermitian
+//    std::cout << GridLogDebug <<"MDeriv:Kappa = "<<kappa<<std::endl;
    RealD factor = -kappa / (double(4 * Nd));
    
    for (int mu = 0; mu < Nd; mu++){
@@ -140,6 +519,7 @@ public:
      // adjoint in the last multiplication
      PokeIndex<LorentzIndex>(der,  -2.0 * factor * der_mu, mu);
    } 
+    std::cout << GridLogDebug <<"MDeriv: Kappa= "<< kappa << " norm2(der) = "<<norm2(der)<<std::endl;
  }

  // separating this temporarily
@@ -159,11 +539,22 @@ public:
      }
      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
    }
+    std::cout << GridLogDebug <<"MDeriv: Kappa= "<< kappa << " norm2(der) = "<<norm2(der)<<std::endl;
  }

  void Minv(const GaugeField& in, GaugeField& inverted){
    HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
    Solver(HermOp, in, inverted);
+    std::cout << GridLogDebug <<"Minv:norm2(inverted) = "<<norm2(inverted)<<std::endl;
+  }
+
+
+  void MinvDeriv(const GaugeField& in, GaugeField& der) {
+    GaugeField X(in.Grid());
+    Minv(in,X);
+    MDeriv(X,der);
+    der *=-1.0;
+    std::cout << GridLogDebug <<"MinvDeriv:norm2(der) = "<<norm2(der)<<std::endl;
  }

  void MSquareRoot(GaugeField& P){
@@ -172,6 +563,7 @@ public:
    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerHalf);
    msCG(HermOp,P,Gp);
    P = Gp; 
+    std::cout << GridLogDebug <<"MSquareRoot:norm2(P) = "<<norm2(P)<<std::endl;
  }

  void MInvSquareRoot(GaugeField& P){
@@ -180,6 +572,7 @@ public:
    ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerInvHalf);
    msCG(HermOp,P,Gp);
    P = Gp; 
+    std::cout << GridLogDebug <<"MInvSquareRoot:norm2(P) = "<<norm2(P)<<std::endl;
  }


--- a/Grid/qcd/utils/CovariantLaplacianRat.h
+++ b/Grid/qcd/utils/CovariantLaplacianRat.h
@@ -0,0 +1,403 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/scalar/CovariantLaplacianRat.h
+
+Copyright (C) 2021
+
+Author: Chulwoo Jung <chulwoo@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 */
+#pragma once 
+#define MIXED_CG
+//enable/disable push_back
+#undef USE_CHRONO 
+
+//#include <roctracer/roctx.h>
+
+NAMESPACE_BEGIN(Grid);
+
+struct LaplacianRatParams {
+
+  RealD offset;
+  int order;
+  std::vector<RealD> a0;
+  std::vector<RealD> a1;
+  std::vector<RealD> b0;
+  std::vector<RealD> b1;
+  RealD b2; //for debugging
+  int   MaxIter;
+  RealD tolerance;
+  int   precision;
+  
+  // constructor 
+  LaplacianRatParams(int ord = 1,
+                  int maxit     = 1000,
+                  RealD tol     = 1.0e-8, 
+                  int precision = 64)
+    : offset(1.), order(ord),b2(1.),
+      MaxIter(maxit),
+      tolerance(tol),
+      precision(precision){ 
+      a0.resize(ord,0.);
+      a1.resize(ord,0.);
+      b0.resize(ord,0.);
+      b1.resize(ord,0.);
+      };
+};
+
+
+
+////////////////////////////////////////////////////////////
+// Laplacian operator L on adjoint fields
+//
+// phi: adjoint field
+// L: D_mu^dag D_mu
+//
+// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag + 
+//                     U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu)
+//                     -2phi(x)]
+//
+// Operator designed to be encapsulated by
+// an HermitianLinearOperator<.. , ..>
+////////////////////////////////////////////////////////////
+
+template <class Impl, class ImplF>
+class LaplacianAdjointRat: public Metric<typename Impl::Field> {
+  OperatorFunction<typename Impl::Field> &Solver;
+  LaplacianRatParams Gparam;
+  LaplacianRatParams Mparam;
+  GridBase *grid;
+  GridBase *grid_f;
+  CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField> LapStencil;
+  CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField> LapStencilF;
+public:
+  INHERIT_GIMPL_TYPES(Impl);
+//   typedef typename GImpl::LinkField GaugeLinkField; \
+//  typedef typename GImpl::Field GaugeField;         
+  typedef typename ImplF::Field GaugeFieldF;
+  typedef typename ImplF::LinkField GaugeLinkFieldF; \
+  GaugeField Usav;
+  GaugeFieldF UsavF;
+  std::vector< std::vector<GaugeLinkField> > prev_solnsM;
+  std::vector< std::vector<GaugeLinkField> > prev_solnsMinv;
+  std::vector< std::vector<GaugeLinkField> > prev_solnsMDeriv;
+  std::vector< std::vector<GaugeLinkField> > prev_solnsMinvDeriv;
+
+	  LaplacianAdjointRat(GridBase* _grid, GridBase* _grid_f, OperatorFunction<GaugeField>& S, LaplacianRatParams& gpar, LaplacianRatParams& mpar)
+    : grid(_grid),grid_f(_grid_f), LapStencil(_grid), LapStencilF(_grid_f), U(Nd, _grid), Solver(S), Gparam(gpar), Mparam(mpar),Usav(_grid), UsavF(_grid_f),
+      prev_solnsM(4),prev_solnsMinv(4),prev_solnsMDeriv(4),prev_solnsMinvDeriv(4) {
+//    std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
+    this->triv=0;
+        
+
+  };
+  LaplacianAdjointRat(){this->triv=0; printf("triv=%d\n",this->Trivial());}
+  void Mdir(const GaugeField&, GaugeField&, int, int){ assert(0);}
+  void MdirAll(const GaugeField&, std::vector<GaugeField> &){ assert(0);}
+  void Mdiag(const GaugeField&, GaugeField&){ assert(0);}
+
+  void ImportGauge(const GaugeField& _U) {
+    RealD total=0.;
+    for (int mu = 0; mu < Nd; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(_U, mu);
+      total += norm2(U[mu]);
+    }
+    Usav = _U;
+    precisionChange(UsavF,Usav);
+    std::cout <<GridLogDebug << "ImportGauge:norm2(_U) = "<<" "<<total<<std::endl;
+  }
+
+  void MDerivLink(const GaugeLinkField& left, const GaugeLinkField& right,
+              GaugeField& der) {
+    std::cout<<GridLogMessage << "MDerivLink start "<< std::endl;
+    RealD factor = -1. / (double(4 * Nd));
+    for (int mu = 0; mu < Nd; mu++) {
+      GaugeLinkField der_mu(der.Grid());
+      der_mu = Zero();
+//      for (int nu = 0; nu < Nd; nu++) {
+//        GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
+//        GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
+        der_mu += U[mu] * Cshift(left, mu, 1) * adj(U[mu]) * right;
+        der_mu += U[mu] * Cshift(right, mu, 1) * adj(U[mu]) * left;
+//      }
+      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
+    }
+//    std::cout << GridLogDebug <<"MDerivLink:  norm2(der) = "<<norm2(der)<<std::endl;
+    std::cout<<GridLogMessage << "MDerivLink end "<< std::endl;
+  }
+
+  void MDerivLink(const GaugeLinkField& left, const GaugeLinkField& right,
+              std::vector<GaugeLinkField> & der) {
+//    std::cout<<GridLogMessage << "MDerivLink "<< std::endl;
+    RealD factor = -1. / (double(4 * Nd));
+
+    for (int mu = 0; mu < Nd; mu++) {
+      GaugeLinkField der_mu(left.Grid());
+      der_mu = Zero();
+        der_mu += U[mu] * Cshift(left, mu, 1) * adj(U[mu]) * right;
+        der_mu += U[mu] * Cshift(right, mu, 1) * adj(U[mu]) * left;
+//      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
+      der[mu] = -factor*der_mu;
+//      std::cout << GridLogDebug <<"MDerivLink:  norm2(der) = "<<norm2(der[mu])<<std::endl;
+        
+    }
+//    std::cout<<GridLogMessage << "MDerivLink end "<< std::endl;
+  }
+
+  void MDerivInt(LaplacianRatParams &par, const GaugeField& left, const GaugeField& right,
+              GaugeField& der ,  std::vector< std::vector<GaugeLinkField> >& prev_solns ) {
+
+// get rid of this please
+    std::cout<<GridLogMessage << "LaplaceStart " <<std::endl;
+    RealD fac =  - 1. / (double(4 * Nd)) ;
+    RealD coef=0.5;
+    LapStencil.GaugeImport(Usav);
+    LapStencilF.GaugeImport(UsavF);
+
+
+    for (int nu=0;nu<Nd;nu++){
+        GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
+        GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
+        GaugeLinkField LMinvMom(left.Grid());
+    
+        GaugeLinkField GMom(left.Grid());
+        GaugeLinkField LMinvGMom(left.Grid());
+    
+        GaugeLinkField AGMom(left.Grid());
+        GaugeLinkField MinvAGMom(left.Grid());
+        GaugeLinkField LMinvAGMom(left.Grid());
+    
+        GaugeLinkField AMinvMom(left.Grid());
+        GaugeLinkField LMinvAMom(left.Grid());
+        GaugeLinkField temp(left.Grid());
+        GaugeLinkField temp2(left.Grid());
+    
+        std::vector<GaugeLinkField> MinvMom(par.order,left.Grid());
+    
+        GaugeLinkField MinvGMom(left.Grid());
+        GaugeLinkField Gtemp(left.Grid());
+        GaugeLinkField Gtemp2(left.Grid());
+    
+    
+        ConjugateGradient<GaugeLinkField> CG(par.tolerance,10000,false);
+    //    ConjugateGradient<GaugeFieldF> CG_f(par.tolerance,10000,false);
+        LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64);
+    
+        ChronoForecast< QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,GaugeLinkField>,GaugeLinkField> , GaugeLinkField> Forecast;
+    
+        GMom = par.offset * right_nu;
+    
+        for(int i =0;i<par.order;i++){
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
+#if USE_CHRONO
+        MinvMom[i] = Forecast(QuadOp, right_nu, prev_solns[nu]);
+#endif
+#ifndef MIXED_CG
+        CG(QuadOp,right_nu,MinvMom[i]);
+#else
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
+    //    QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
+        MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
+        MixedCG.InnerTolerance=par.tolerance;
+        MixedCG(right_nu,MinvMom[i]);
+    #endif
+    #if USE_CHRONO
+        prev_solns[nu].push_back(MinvMom[i]);
+    #endif
+        
+        GMom += par.a0[i]*MinvMom[i]; 
+        LapStencil.M(MinvMom[i],Gtemp2);
+        GMom += par.a1[i]*fac*Gtemp2; 
+        }
+        for(int i =0;i<par.order;i++){
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
+    
+        MinvGMom = Forecast(QuadOp, GMom, prev_solns[nu]);
+    #ifndef MIXED_CG
+        CG(QuadOp,GMom,MinvGMom);
+        LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
+        CG(QuadOp,right_nu,MinvMom[i]);
+    #else
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
+    //    QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
+        MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
+        MixedCG.InnerTolerance=par.tolerance;
+        MixedCG(GMom,MinvGMom);
+        LapStencil.M(MinvGMom, Gtemp2); LMinvGMom=fac*Gtemp2;
+    //    Laplacian.M(MinvGMom, LMinvGMom);
+        MixedCG(right_nu,MinvMom[i]);
+    #endif
+#if USE_CHRONO
+        prev_solns[nu].push_back(MinvGMom);
+#endif
+    
+        LapStencil.M(MinvMom[i], Gtemp2); LMinvMom=fac*Gtemp2;
+        AMinvMom = par.a1[i]*LMinvMom;
+        AMinvMom += par.a0[i]*MinvMom[i];
+    
+        LapStencil.M(AMinvMom, Gtemp2); LMinvAMom=fac*Gtemp2;
+        LapStencil.M(MinvGMom, Gtemp2); temp=fac*Gtemp2;
+        MinvAGMom = par.a1[i]*temp;
+        MinvAGMom += par.a0[i]*MinvGMom;
+        LapStencil.M(MinvAGMom, Gtemp2); LMinvAGMom=fac*Gtemp2;
+    
+    
+        GaugeField tempDer(left.Grid());
+        std::vector<GaugeLinkField> DerLink(Nd,left.Grid());
+        std::vector<GaugeLinkField> tempDerLink(Nd,left.Grid());
+
+        std::cout<<GridLogMessage << "force contraction "<< i <<std::endl;
+    //    roctxRangePushA("RMHMC force contraction");
+ #if 0
+        MDerivLink(GMom,MinvMom[i],tempDer); der += coef*2*par.a1[i]*tempDer;
+        MDerivLink(left_nu,MinvGMom,tempDer); der += coef*2*par.a1[i]*tempDer;
+        MDerivLink(LMinvAGMom,MinvMom[i],tempDer); der += coef*-2.*par.b2*tempDer;
+        MDerivLink(LMinvAMom,MinvGMom,tempDer); der += coef*-2.*par.b2*tempDer;
+        MDerivLink(MinvAGMom,LMinvMom,tempDer); der += coef*-2.*par.b2*tempDer;
+        MDerivLink(AMinvMom,LMinvGMom,tempDer); der += coef*-2.*par.b2*tempDer;
+        MDerivLink(MinvAGMom,MinvMom[i],tempDer); der += coef*-2.*par.b1[i]*tempDer;
+        MDerivLink(AMinvMom,MinvGMom,tempDer); der += coef*-2.*par.b1[i]*tempDer;
+#else
+	for (int mu=0;mu<Nd;mu++) DerLink[mu]=Zero();
+        MDerivLink(GMom,MinvMom[i],tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*2*par.a1[i]*tempDerLink[mu];
+        MDerivLink(left_nu,MinvGMom,tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*2*par.a1[i]*tempDerLink[mu];
+        MDerivLink(LMinvAGMom,MinvMom[i],tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*-2.*par.b2*tempDerLink[mu];
+        MDerivLink(LMinvAMom,MinvGMom,tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*-2.*par.b2*tempDerLink[mu];
+        MDerivLink(MinvAGMom,LMinvMom,tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*-2.*par.b2*tempDerLink[mu];
+        MDerivLink(AMinvMom,LMinvGMom,tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*-2.*par.b2*tempDerLink[mu];
+        MDerivLink(MinvAGMom,MinvMom[i],tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*-2.*par.b1[i]*tempDerLink[mu];
+        MDerivLink(AMinvMom,MinvGMom,tempDerLink); 	for (int mu=0;mu<Nd;mu++) DerLink[mu] += coef*-2.*par.b1[i]*tempDerLink[mu];
+//      PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
+        for (int mu=0;mu<Nd;mu++) PokeIndex<LorentzIndex>(tempDer, tempDerLink[mu], mu);
+
+	der += tempDer;
+#endif
+        std::cout<<GridLogMessage << "coef =  force contraction "<< i << "done "<< coef <<std::endl;
+    //    roctxRangePop();
+    
+        }
+    }
+    std::cout<<GridLogMessage << "LaplaceEnd " <<std::endl;
+//  exit(-42);
+  }
+
+  void MDeriv(const GaugeField& in, GaugeField& der) {
+    MDeriv(in,in, der);
+  }
+
+  void MDeriv(const GaugeField& left, const GaugeField& right,
+              GaugeField& der) {
+
+    der=Zero();
+    MDerivInt(Mparam, left, right, der,prev_solnsMDeriv );
+    std::cout <<GridLogDebug << "MDeriv:norm2(der) = "<<norm2(der)<<std::endl;
+  }
+
+  void MinvDeriv(const GaugeField& in, GaugeField& der) {
+    std::vector< std::vector<GaugeLinkField> > prev_solns(4);
+    der=Zero();
+    MDerivInt(Gparam, in, in, der,prev_solnsMinvDeriv);
+    std::cout <<GridLogDebug << "MinvDeriv:norm2(der) = "<<norm2(der)<<std::endl;
+  }
+
+
+  void MSquareRootInt(LaplacianRatParams &par, GaugeField& P, std::vector< std::vector<GaugeLinkField> > & prev_solns ){
+
+    std::cout<<GridLogMessage << "LaplaceStart " <<std::endl;
+    RealD fac = -1. / (double(4 * Nd));
+    LapStencil.GaugeImport(Usav);
+    LapStencilF.GaugeImport(UsavF);
+    for(int nu=0; nu<Nd;nu++){
+        GaugeLinkField P_nu = PeekIndex<LorentzIndex>(P, nu);
+        GaugeLinkField Gp(P.Grid());
+        Gp = par.offset * P_nu;
+        ConjugateGradient<GaugeLinkField> CG(par.tolerance,10000);
+    //    ConjugateGradient<GaugeLinkFieldF> CG_f(1.0e-8,10000);
+    
+        ChronoForecast< QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> , GaugeLinkField> Forecast;
+    
+        GaugeLinkField Gtemp(P.Grid());
+        GaugeLinkField Gtemp2(P.Grid());
+    
+    
+        for(int i =0;i<par.order;i++){
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<Impl,typename Impl::LinkField>,GaugeLinkField> QuadOp(LapStencil,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
+    
+        Gtemp = Forecast(QuadOp, P_nu, prev_solns[nu]);
+    #ifndef MIXED_CG
+        CG(QuadOp,P_nu,Gtemp);
+    #else
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<ImplF,typename ImplF::LinkField>,GaugeLinkFieldF> QuadOpF(LapStencilF,par.b0[i],fac*par.b1[i],fac*fac*par.b2);
+    //    QuadLinearOperator<LaplacianAdjointField<ImplF>,GaugeFieldF> QuadOpF(LapStencilF,par.b0[i],par.b1[i],par.b2);
+        MixedPrecisionConjugateGradient<GaugeLinkField,GaugeLinkFieldF> MixedCG(par.tolerance,10000,10000,grid_f,QuadOpF,QuadOp);
+        MixedCG.InnerTolerance=par.tolerance;
+        MixedCG(P_nu,Gtemp);
+    #endif
+    #if USE_CHRONO
+        prev_solns[nu].push_back(Gtemp);
+    #endif
+    
+        Gp += par.a0[i]*Gtemp; 
+        LapStencil.M(Gtemp,Gtemp2);
+        Gp += par.a1[i]*fac*Gtemp2; 
+        }
+        PokeIndex<LorentzIndex>(P, Gp, nu);
+    }
+    std::cout<<GridLogMessage << "LaplaceEnd " <<std::endl;
+  }
+
+  void MSquareRoot(GaugeField& P){
+    std::vector< std::vector<GaugeLinkField> > prev_solns(4);
+    MSquareRootInt(Mparam,P,prev_solns);
+    std::cout <<GridLogDebug << "MSquareRoot:norm2(P) = "<<norm2(P)<<std::endl;
+  }
+
+  void MInvSquareRoot(GaugeField& P){
+    std::vector< std::vector<GaugeLinkField> > prev_solns(4);
+    MSquareRootInt(Gparam,P,prev_solns);
+    std::cout <<GridLogDebug << "MInvSquareRoot:norm2(P) = "<<norm2(P)<<std::endl;
+  }
+
+  void M(const GaugeField& in, GaugeField& out) {
+      out = in;
+      std::vector< std::vector<GaugeLinkField> > prev_solns(4);
+      MSquareRootInt(Mparam,out,prev_solns);
+      MSquareRootInt(Mparam,out,prev_solns);
+      std::cout <<GridLogDebug << "M:norm2(out) = "<<norm2(out)<<std::endl;
+  }
+
+  void Minv(const GaugeField& in, GaugeField& inverted){
+      inverted = in;
+      std::vector< std::vector<GaugeLinkField> > prev_solns(4);
+      MSquareRootInt(Gparam,inverted,prev_solns);
+      MSquareRootInt(Gparam,inverted,prev_solns);
+      std::cout <<GridLogDebug << "Minv:norm2(inverted) = "<<norm2(inverted)<<std::endl;
+  }
+
+
+
+private:
+  std::vector<GaugeLinkField> U;
+};
+#undef MIXED_CG
+
+NAMESPACE_END(Grid);
--- a/Grid/qcd/utils/Metric.h
+++ b/Grid/qcd/utils/Metric.h
@@ -7,6 +7,7 @@ Source file: ./lib/qcd/hmc/integrators/Integrator.h
 Copyright (C) 2015

 Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Chulwoo Jung <chulwoo@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
@@ -33,7 +34,12 @@ NAMESPACE_BEGIN(Grid);

 template <typename Field> 
 class Metric{
+protected:
+  int triv;
 public:
+  Metric(){this->triv=1;}
+  int Trivial(){ return triv;}
+//printf("Metric::Trivial=%d\n",triv); ;
  virtual void ImportGauge(const Field&)   = 0;
  virtual void M(const Field&, Field&)     = 0;
  virtual void Minv(const Field&, Field&)  = 0;
@@ -41,6 +47,8 @@ public:
  virtual void MInvSquareRoot(Field&) = 0;
  virtual void MDeriv(const Field&, Field&) = 0;
  virtual void MDeriv(const Field&, const Field&, Field&) = 0;
+  virtual void MinvDeriv(const Field&, Field&) = 0;
+//  virtual void MinvDeriv(const Field&, const Field&, Field&) = 0;
 };


@@ -48,23 +56,36 @@ public:
 template <typename Field>
 class TrivialMetric : public Metric<Field>{
 public:
+//  TrivialMetric(){this->triv=1;printf("TrivialMetric::triv=%d\n",this->Trivial());}
  virtual void ImportGauge(const Field&){};
  virtual void M(const Field& in, Field& out){
+//    printf("M:norm=%0.15e\n",norm2(in));
+    std::cout << GridLogIntegrator << " M:norm(in)= " << std::sqrt(norm2(in)) << std::endl;
    out = in;
  }
  virtual void Minv(const Field& in, Field& out){
+    std::cout << GridLogIntegrator << " Minv:norm(in)= " << std::sqrt(norm2(in)) << std::endl;
    out = in;
  }
  virtual void MSquareRoot(Field& P){
+    std::cout << GridLogIntegrator << " MSquareRoot:norm(P)= " << std::sqrt(norm2(P)) << std::endl;
    // do nothing
  }
  virtual void MInvSquareRoot(Field& P){
+    std::cout << GridLogIntegrator << " MInvSquareRoot:norm(P)= " << std::sqrt(norm2(P)) << std::endl;
    // do nothing
  }
  virtual void MDeriv(const Field& in, Field& out){
+    std::cout << GridLogIntegrator << " MDeriv:norm(in)= " << std::sqrt(norm2(in)) << std::endl;
+    out = Zero();
+  }
+  virtual void MinvDeriv(const Field& in, Field& out){
+    std::cout << GridLogIntegrator << " MinvDeriv:norm(in)= " << std::sqrt(norm2(in)) << std::endl;
    out = Zero();
  }
  virtual void MDeriv(const Field& left, const Field& right, Field& out){
+    std::cout << GridLogIntegrator << " MDeriv:norm(left)= " << std::sqrt(norm2(left)) << std::endl;
+    std::cout << GridLogIntegrator << " MDeriv:norm(right)= " << std::sqrt(norm2(right)) << std::endl;
    out = Zero();
  }

@@ -101,14 +122,15 @@ public:
    // Generate gaussian momenta
    Implementation::generate_momenta(Mom, sRNG, pRNG);
    // Modify the distribution with the metric
+//    if(M.Trivial()) return;
    M.MSquareRoot(Mom);

    if (1) {
      // Auxiliary momenta
      // do nothing if trivial, so hide in the metric
      MomentaField AuxMomTemp(Mom.Grid());
-      Implementation::generate_momenta(AuxMom, sRNG, pRNG);
-      Implementation::generate_momenta(AuxField, sRNG, pRNG);
+      Implementation::generate_momenta(AuxMom, sRNG,pRNG);
+      Implementation::generate_momenta(AuxField, sRNG,pRNG);
      // Modify the distribution with the metric
      // Aux^dag M Aux
      M.MInvSquareRoot(AuxMom);  // AuxMom = M^{-1/2} AuxMomTemp
@@ -117,11 +139,12 @@ public:

  // Correct
  RealD MomentaAction(){
+    static RealD Saux=0.,Smom=0.;
    MomentaField inv(Mom.Grid());
    inv = Zero();
    M.Minv(Mom, inv);
-    LatticeComplex Hloc(Mom.Grid());
-    Hloc = Zero();
+    LatticeComplex Hloc(Mom.Grid()); Hloc = Zero();
+    LatticeComplex Hloc2(Mom.Grid()); Hloc2 = Zero();
    for (int mu = 0; mu < Nd; mu++) {
      // This is not very general
      // hide in the metric
@@ -129,8 +152,15 @@ public:
      auto inv_mu = PeekIndex<LorentzIndex>(inv, mu);
      Hloc += trace(Mom_mu * inv_mu);
    }
+    auto Htmp1 = TensorRemove(sum(Hloc));
+    std::cout << GridLogMessage << "S:dSmom = " << Htmp1.real()-Smom << "\n";
+    Smom=Htmp1.real()/HMC_MOMENTUM_DENOMINATOR;
+    

-    if (1) {
+    
+
+//    if(!M.Trivial()) 
+    {
      // Auxiliary Fields
      // hide in the metric
      M.M(AuxMom, inv);
@@ -140,13 +170,18 @@ public:
        auto inv_mu = PeekIndex<LorentzIndex>(inv, mu);
        auto am_mu = PeekIndex<LorentzIndex>(AuxMom, mu);
        auto af_mu = PeekIndex<LorentzIndex>(AuxField, mu);
-        Hloc += trace(am_mu * inv_mu);// p M p
-        Hloc += trace(af_mu * af_mu);
+        Hloc += trace(am_mu * inv_mu);
+        Hloc2 += trace(af_mu * af_mu);
      }
    }
+    auto Htmp2 = TensorRemove(sum(Hloc))-Htmp1;
+    std::cout << GridLogMessage << "S:dSaux = " << Htmp2.real()-Saux << "\n";
+    Saux=Htmp2.real();

-    auto Hsum = TensorRemove(sum(Hloc));
-    return Hsum.real();
+    auto Hsum = TensorRemove(sum(Hloc))/HMC_MOMENTUM_DENOMINATOR;
+    auto Hsum2 = TensorRemove(sum(Hloc2));
+    std::cout << GridLogIntegrator << "MomentaAction: " <<  Hsum.real()+Hsum2.real() << std::endl;
+    return Hsum.real()+Hsum2.real();
  }

  // Correct
@@ -157,15 +192,17 @@ public:
    MomentaField MDer(in.Grid());
    MomentaField X(in.Grid());
    X = Zero();
-    M.Minv(in, X);  // X = G in
-    M.MDeriv(X, MDer);  // MDer = U * dS/dU
-    der = Implementation::projectForce(MDer);  // Ta if gauge fields
+    M.MinvDeriv(in, MDer);  // MDer = U * dS/dU
+    der = -1.0* Implementation::projectForce(MDer);  // Ta if gauge fields
+//    std::cout << GridLogIntegrator << " DerivativeU: norm(in)= " << std::sqrt(norm2(in)) << std::endl;
+//    std::cout << GridLogIntegrator << " DerivativeU: norm(der)= " << std::sqrt(norm2(der)) << std::endl;
    
  }

  void AuxiliaryFieldsDerivative(MomentaField& der){
    der = Zero();
-    if (1){
+//    if(!M.Trivial()) 
+    {
      // Auxiliary fields
      MomentaField der_temp(der.Grid());
      MomentaField X(der.Grid());
@@ -173,6 +210,7 @@ public:
      //M.M(AuxMom, X); // X = M Aux
      // Two derivative terms
      // the Mderiv need separation of left and right terms
+    std::cout << GridLogIntegrator << " AuxiliaryFieldsDerivative:norm(AuxMom)= " << std::sqrt(norm2(AuxMom)) << std::endl;
      M.MDeriv(AuxMom, der); 


@@ -180,6 +218,7 @@ public:
      //M.MDeriv(X, AuxMom, der_temp); der += der_temp;

      der = -1.0*Implementation::projectForce(der);
+      std::cout << GridLogIntegrator << " AuxiliaryFieldsDerivative:norm(der)= " << std::sqrt(norm2(der)) << std::endl;
    }
  }

@@ -189,22 +228,28 @@ public:
    // is the projection necessary here?
    // no for fields in the algebra
    der = Implementation::projectForce(der); 
+    std::cout << GridLogIntegrator << " DerivativeP:norm(der)= " << std::sqrt(norm2(der)) << std::endl;
  }

  void update_auxiliary_momenta(RealD ep){
-    if(1){
-      AuxMom -= ep * AuxField;
+      std::cout << GridLogIntegrator << "AuxMom update_auxiliary_fields: " << std::sqrt(norm2(AuxMom)) << std::endl;
+      std::cout << GridLogIntegrator << "AuxField update_auxiliary_fields: " << std::sqrt(norm2(AuxField)) << std::endl;
+    {
+      AuxMom -= ep * AuxField * HMC_MOMENTUM_DENOMINATOR;
+      std::cout << GridLogIntegrator << "AuxMom update_auxiliary_fields: " << std::sqrt(norm2(AuxMom)) << std::endl;
    }
  }

  void update_auxiliary_fields(RealD ep){
-    if (1) {
+//    if(!M.Trivial()) 
+    {
      MomentaField tmp(AuxMom.Grid());
      MomentaField tmp2(AuxMom.Grid());
      M.M(AuxMom, tmp);
      // M.M(tmp, tmp2);
      AuxField += ep * tmp;  // M^2 AuxMom
      // factor of 2?
+      std::cout << GridLogIntegrator << "AuxField update_auxiliary_fields: " << std::sqrt(norm2(AuxField)) << std::endl;
    }
  }

--- a/HMC/Mobius2p1p1fEOFA_4Gev.cc
+++ b/HMC/Mobius2p1p1fEOFA_4Gev.cc
@@ -0,0 +1,637 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: 
+
+Copyright (C) 2015-2016
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: Guido Cossu
+Author: David Murphy
+Author: Chulwoo Jung <chulwoo@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>
+
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
+#define MIXED_PRECISION
+#endif
+// second level EOFA
+#undef EOFA_H
+#undef USE_OBC
+#define DO_IMPLICIT
+
+NAMESPACE_BEGIN(Grid);
+
+  /*
+   * Need a plan for gauge field update for mixed precision in HMC                      (2x speed up)
+   *    -- Store the single prec action operator.
+   *    -- Clone the gauge field from the operator function argument.
+   *    -- Build the mixed precision operator dynamically from the passed operator and single prec clone.
+   */
+
+  template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class  SchurOperatorF> 
+  class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> {
+  public:
+    typedef typename FermionOperatorD::FermionField FieldD;
+    typedef typename FermionOperatorF::FermionField FieldF;
+
+    using OperatorFunction<FieldD>::operator();
+
+    RealD   Tolerance;
+    RealD   InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
+    Integer MaxInnerIterations;
+    Integer MaxOuterIterations;
+    GridBase* SinglePrecGrid4; //Grid for single-precision fields
+    GridBase* SinglePrecGrid5; //Grid for single-precision fields
+    RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance
+
+    FermionOperatorF &FermOpF;
+    FermionOperatorD &FermOpD;;
+    SchurOperatorF &LinOpF;
+    SchurOperatorD &LinOpD;
+
+    Integer TotalInnerIterations; //Number of inner CG iterations
+    Integer TotalOuterIterations; //Number of restarts
+    Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
+
+    MixedPrecisionConjugateGradientOperatorFunction(RealD tol, 
+						    Integer maxinnerit, 
+						    Integer maxouterit, 
+						    GridBase* _sp_grid4, 
+						    GridBase* _sp_grid5, 
+						    FermionOperatorF &_FermOpF,
+						    FermionOperatorD &_FermOpD,
+						    SchurOperatorF   &_LinOpF,
+						    SchurOperatorD   &_LinOpD): 
+      LinOpF(_LinOpF),
+      LinOpD(_LinOpD),
+      FermOpF(_FermOpF),
+      FermOpD(_FermOpD),
+      Tolerance(tol), 
+      InnerTolerance(tol), 
+      MaxInnerIterations(maxinnerit), 
+      MaxOuterIterations(maxouterit), 
+      SinglePrecGrid4(_sp_grid4),
+      SinglePrecGrid5(_sp_grid5),
+      OuterLoopNormMult(100.) 
+    { 
+      /* Debugging instances of objects; references are stored
+      std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl;
+      std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl;
+      std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl;
+      std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl;
+      */
+    };
+
+    void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
+
+      std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
+
+      SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
+      
+      //      std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
+      //      std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
+      // Assumption made in code to extract gauge field
+      // We could avoid storing LinopD reference alltogether ?
+      assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
+
+      ////////////////////////////////////////////////////////////////////////////////////
+      // Must snarf a single precision copy of the gauge field in Linop_d argument
+      ////////////////////////////////////////////////////////////////////////////////////
+      typedef typename FermionOperatorF::GaugeField GaugeFieldF;
+      typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
+      typedef typename FermionOperatorD::GaugeField GaugeFieldD;
+      typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
+
+      GridBase * GridPtrF = SinglePrecGrid4;
+      GridBase * GridPtrD = FermOpD.Umu.Grid();
+      GaugeFieldF     U_f  (GridPtrF);
+      GaugeLinkFieldF Umu_f(GridPtrF);
+      //      std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
+      //      std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
+
+      ////////////////////////////////////////////////////////////////////////////////////
+      // Moving this to a Clone method of fermion operator would allow to duplicate the 
+      // physics parameters and decrease gauge field copies
+      ////////////////////////////////////////////////////////////////////////////////////
+      GaugeLinkFieldD Umu_d(GridPtrD);
+      for(int mu=0;mu<Nd*2;mu++){ 
+	Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
+	precisionChange(Umu_f,Umu_d);
+	PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
+      }
+      pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
+      pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
+
+      ////////////////////////////////////////////////////////////////////////////////////
+      // Make a mixed precision conjugate gradient
+      ////////////////////////////////////////////////////////////////////////////////////
+      MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
+      std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
+      MPCG(src,psi);
+    }
+  };
+
+NAMESPACE_END(Grid);
+
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+
+  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
+  typedef WilsonImplR FermionImplPolicy;
+  typedef MobiusFermionD FermionAction;
+  typedef MobiusFermionF FermionActionF;
+  typedef MobiusEOFAFermionD FermionEOFAAction;
+  typedef MobiusEOFAFermionF FermionEOFAActionF;
+  typedef typename FermionAction::FermionField FermionField;
+  typedef typename FermionActionF::FermionField FermionFieldF;
+
+  typedef Grid::XmlReader       Serialiser;
+  
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+  HMCparameters HMCparams;
+#if 1
+  {
+    XmlReader  HMCrd("HMCparameters.xml");
+    read(HMCrd,"HMCparameters",HMCparams);
+  }
+#else
+  {
+//    HMCparameters HMCparams;
+  //  "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
+  //  HMCparams.StartingType     =std::string("ColdStart");
+    HMCparams.StartingType     =std::string("CheckpointStart");
+    HMCparams.StartTrajectory  =7;
+    HMCparams.SW  =4;
+    HMCparams.Trajectories     =1000;
+    HMCparams.NoMetropolisUntil=0;
+    HMCparams.MD.name          =std::string("Force Gradient");
+    HMCparams.MD.MDsteps       = 10;
+    HMCparams.MD.trajL         = 1.0;
+  }
+#endif
+
+#ifdef DO_IMPLICIT
+//    typedef GenericHMCRunner<ImplicitLeapFrog> HMCWrapper; 
+  typedef GenericHMCRunner<ImplicitMinimumNorm2> HMCWrapper; 
+  HMCparams.MD.name          =std::string("ImplicitMinimumNorm2");
+#else
+//  typedef GenericHMCRunner<LeapFrog> HMCWrapper; 
+  typedef GenericHMCRunner<ForceGradient> HMCWrapper; 
+//  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; 
+  HMCparams.MD.name          =std::string("ForceGradient");
+#endif
+
+  std::cout << GridLogMessage<< HMCparams <<std::endl;
+  HMCWrapper TheHMC(HMCparams);
+  TheHMC.ReadCommandLine(argc, argv);
+  { 
+    XmlWriter HMCwr("HMCparameters.xml.out");
+    write(HMCwr,"HMCparameters",TheHMC.Parameters);
+  }
+
+  // Grid from the command line arguments --grid and --mpi
+  TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
+  
+  CheckpointerParameters CPparams;
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix    = "ckpoint_rng";
+  CPparams.saveInterval  = 1;
+  CPparams.format        = "IEEE64BIG";
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  // here there is too much indirection 
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+  //////////////////////////////////////////////
+
+  const int Ls      = 12;
+  Real beta         = 5.983;
+  std::cout << GridLogMessage << " beta  "<< beta << std::endl;
+  Real light_mass   = 0.00049;
+  Real strange_mass = 0.0158;
+  Real charm_mass = 0.191;
+  Real pv_mass    = 1.0;
+  RealD M5  = 1.4;
+  RealD b   = 2.0; 
+  RealD c   = 1.0;
+
+  // Copied from paper
+//  std::vector<Real> hasenbusch({ 0.045 }); // Paper values from F1 incorrect run
+  std::vector<Real> hasenbusch({ 0.0038, 0.0145, 0.045, 0.108 , 0.25, 0.51 }); // Paper values from F1 incorrect run
+  std::vector<Real> hasenbusch2({ 0.4 }); // Paper values from F1 incorrect run
+
+//  RealD eofa_mass=0.05 ;
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////
+  //Bad choices with large dH. Equalising force L2 norm was not wise.
+  ///////////////////////////////////////////////////////////////////////////////////////////////
+  //std::vector<Real> hasenbusch({ 0.03, 0.2, 0.3, 0.5, 0.8 }); 
+
+  auto GridPtr   = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+  auto FGrid     = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
+  auto FrbGrid   = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
+
+  Coordinate latt  = GridDefaultLatt();
+  Coordinate mpi   = GridDefaultMpi();
+  Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd());
+  Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd());
+//  auto GridPtrF   = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
+  auto UGrid_f    = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
+  auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
+  auto FGridF     = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_f);
+  auto FrbGridF   = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_f);
+
+
+#ifndef USE_OBC
+//  IwasakiGaugeActionR GaugeAction(beta);
+  WilsonGaugeActionR GaugeAction(beta);
+#else
+  std::vector<Complex> boundaryG = {1,1,1,0};
+  WilsonGaugeActionR::ImplParams ParamsG(boundaryG);
+  WilsonGaugeActionR GaugeAction(beta,ParamsG);
+#endif
+
+  // temporarily need a gauge field
+  LatticeGaugeField U(GridPtr);
+  LatticeGaugeFieldF UF(UGrid_f);
+
+  // These lines are unecessary if BC are all periodic
+#ifndef USE_OBC
+  std::vector<Complex> boundary = {1,1,1,-1};
+#else
+  std::vector<Complex> boundary = {1,1,1,0};
+#endif
+  FermionAction::ImplParams Params(boundary);
+  FermionActionF::ImplParams ParamsF(boundary);
+  
+  double ActionStoppingCondition     = 1e-8;
+  double DerivativeStoppingCondition = 1e-8;
+  double MaxCGIterations =  100000;
+
+  ////////////////////////////////////
+  // Collect actions
+  ////////////////////////////////////
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  ActionLevel<HMCWrapper::Field> Level2(HMCparams.SW);
+
+  ////////////////////////////////////
+  // Strange action
+  ////////////////////////////////////
+  typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> LinearOperatorF;
+  typedef SchurDiagMooeeOperator<FermionAction ,FermionField > LinearOperatorD;
+  typedef SchurDiagMooeeOperator<FermionEOFAActionF,FermionFieldF> LinearOperatorEOFAF;
+  typedef SchurDiagMooeeOperator<FermionEOFAAction ,FermionField > LinearOperatorEOFAD;
+
+  typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusFermionD,MobiusFermionF,LinearOperatorD,LinearOperatorF> MxPCG;
+  typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusEOFAFermionD,MobiusEOFAFermionF,LinearOperatorEOFAD,LinearOperatorEOFAF> MxPCG_EOFA;
+
+  // DJM: setup for EOFA ratio (Mobius)
+  OneFlavourRationalParams OFRp;
+  OFRp.lo       = 0.99; // How do I know this on F1?
+  OFRp.hi       = 20;
+  OFRp.MaxIter  = 100000;
+  OFRp.tolerance= 1.0e-12;
+  OFRp.degree   = 12;
+  OFRp.precision= 50;
+
+  
+  MobiusEOFAFermionD Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, charm_mass, 0.0, -1, M5, b, c);
+  MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, strange_mass, strange_mass, charm_mass, 0.0, -1, M5, b, c);
+  MobiusEOFAFermionD Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , charm_mass, strange_mass,      charm_mass, -1.0, 1, M5, b, c);
+  MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, charm_mass, strange_mass,      charm_mass, -1.0, 1, M5, b, c);
+  
+#ifdef EOFA_H
+  MobiusEOFAFermionD Strange2_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , eofa_mass, eofa_mass, charm_mass , 0.0, -1, M5, b, c);
+  MobiusEOFAFermionF Strange2_Op_LF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, eofa_mass, eofa_mass, charm_mass , 0.0, -1, M5, b, c);
+  MobiusEOFAFermionD Strange2_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , charm_mass , eofa_mass,      charm_mass , -1.0, 1, M5, b, c);
+  MobiusEOFAFermionF Strange2_Op_RF(UF, *FGridF, *FrbGridF, *UGrid_f, *GridRBPtrF, charm_mass , eofa_mass,      charm_mass , -1.0, 1, M5, b, c);
+#endif
+
+  ConjugateGradient<FermionField>      ActionCG(ActionStoppingCondition,MaxCGIterations);
+  ConjugateGradient<FermionField>  DerivativeCG(DerivativeStoppingCondition,MaxCGIterations);
+#ifdef MIXED_PRECISION
+  const int MX_inner = 50000;
+
+  // Mixed precision EOFA
+  LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L);
+  LinearOperatorEOFAD Strange_LinOp_R (Strange_Op_R);
+  LinearOperatorEOFAF Strange_LinOp_LF(Strange_Op_LF);
+  LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF);
+
+#ifdef EOFA_H
+  // Mixed precision EOFA
+  LinearOperatorEOFAD Strange2_LinOp_L (Strange2_Op_L);
+  LinearOperatorEOFAD Strange2_LinOp_R (Strange2_Op_R);
+  LinearOperatorEOFAF Strange2_LinOp_LF(Strange2_Op_LF);
+  LinearOperatorEOFAF Strange2_LinOp_RF(Strange2_Op_RF);
+#endif
+
+  MxPCG_EOFA ActionCGL(ActionStoppingCondition,
+		       MX_inner,
+		       MaxCGIterations,
+		       UGrid_f,
+		       FrbGridF,
+		       Strange_Op_LF,Strange_Op_L,
+		       Strange_LinOp_LF,Strange_LinOp_L);
+
+#ifdef EOFA_H
+  MxPCG_EOFA ActionCGL2(ActionStoppingCondition,
+		       MX_inner,
+		       MaxCGIterations,
+		       UGrid_f,
+		       FrbGridF,
+		       Strange2_Op_LF,Strange2_Op_L,
+		       Strange2_LinOp_LF,Strange2_LinOp_L);
+#endif
+
+  MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition,
+			   MX_inner,
+			   MaxCGIterations,
+			   UGrid_f,
+			   FrbGridF,
+			   Strange_Op_LF,Strange_Op_L,
+			   Strange_LinOp_LF,Strange_LinOp_L);
+
+#ifdef EOFA_H
+  MxPCG_EOFA DerivativeCGL2(DerivativeStoppingCondition,
+			   MX_inner,
+			   MaxCGIterations,
+			   UGrid_f,
+			   FrbGridF,
+			   Strange2_Op_LF,Strange2_Op_L,
+			   Strange2_LinOp_LF,Strange2_LinOp_L);
+#endif
+  
+  MxPCG_EOFA ActionCGR(ActionStoppingCondition,
+		       MX_inner,
+		       MaxCGIterations,
+		       UGrid_f,
+		       FrbGridF,
+		       Strange_Op_RF,Strange_Op_R,
+		       Strange_LinOp_RF,Strange_LinOp_R);
+  
+#ifdef EOFA_H
+  MxPCG_EOFA ActionCGR2(ActionStoppingCondition,
+		       MX_inner,
+		       MaxCGIterations,
+		       UGrid_f,
+		       FrbGridF,
+		       Strange2_Op_RF,Strange2_Op_R,
+		       Strange2_LinOp_RF,Strange2_LinOp_R);
+#endif
+  
+  MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
+			   MX_inner,
+			   MaxCGIterations,
+			   UGrid_f,
+			   FrbGridF,
+			   Strange_Op_RF,Strange_Op_R,
+			   Strange_LinOp_RF,Strange_LinOp_R);
+  
+#ifdef EOFA_H
+  MxPCG_EOFA DerivativeCGR2(DerivativeStoppingCondition,
+			   MX_inner,
+			   MaxCGIterations,
+			   UGrid_f,
+			   FrbGridF,
+			   Strange2_Op_RF,Strange2_Op_R,
+			   Strange2_LinOp_RF,Strange2_LinOp_R);
+#endif
+  
+  ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> 
+    EOFA(Strange_Op_L, Strange_Op_R, 
+	 ActionCG, 
+	 ActionCGL, ActionCGR,
+	 DerivativeCGL, DerivativeCGR,
+	 OFRp, true);
+  
+#ifdef EOFA_H
+  ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> 
+    EOFA2(Strange2_Op_L, Strange2_Op_R, 
+	 ActionCG, 
+	 ActionCGL2, ActionCGR2,
+	 DerivativeCGL2, DerivativeCGR2,
+	 OFRp, true);
+#endif
+
+  Level1.push_back(&EOFA);
+#ifdef EOFA_H
+  Level1.push_back(&EOFA2);
+#endif
+
+#else
+  ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> 
+    EOFA(Strange_Op_L, Strange_Op_R, 
+	 ActionCG, 
+	 ActionCG, ActionCG,
+	 ActionCG, ActionCG,
+	 //         DerivativeCG, DerivativeCG,
+	 OFRp, true);
+  Level1.push_back(&EOFA);
+#endif
+
+  ////////////////////////////////////
+  // up down action
+  ////////////////////////////////////
+  std::vector<Real> light_den;
+  std::vector<Real> light_num;
+
+  int n_hasenbusch = hasenbusch.size();
+  light_den.push_back(light_mass);
+  for(int h=0;h<n_hasenbusch;h++){
+    light_den.push_back(hasenbusch[h]);
+    light_num.push_back(hasenbusch[h]);
+  }
+  light_num.push_back(pv_mass);
+
+  int n_hasenbusch2 = hasenbusch2.size();
+  light_den.push_back(charm_mass);
+  for(int h=0;h<n_hasenbusch2;h++){
+    light_den.push_back(hasenbusch2[h]);
+    light_num.push_back(hasenbusch2[h]);
+  }
+  light_num.push_back(pv_mass);
+
+
+  //////////////////////////////////////////////////////////////
+  // Forced to replicate the MxPCG and DenominatorsF etc.. because
+  // there is no convenient way to "Clone" physics params from double op
+  // into single op for any operator pair.
+  // Same issue prevents using MxPCG in the Heatbath step
+  //////////////////////////////////////////////////////////////
+  std::vector<FermionAction *> Numerators;
+  std::vector<FermionAction *> Denominators;
+  std::vector<TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy> *> Quotients;
+  std::vector<MxPCG *> ActionMPCG;
+  std::vector<MxPCG *> MPCG;
+  std::vector<FermionActionF *> DenominatorsF;
+  std::vector<LinearOperatorD *> LinOpD;
+  std::vector<LinearOperatorF *> LinOpF; 
+
+  for(int h=0;h<light_den.size();h++){
+
+    std::cout << GridLogMessage << " 2f quotient Action  "<< light_num[h] << " / " << light_den[h]<< std::endl;
+
+    Numerators.push_back  (new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_num[h],M5,b,c, Params));
+    Denominators.push_back(new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_den[h],M5,b,c, Params));
+
+#ifdef MIXED_PRECISION
+    ////////////////////////////////////////////////////////////////////////////
+    // Mixed precision CG for 2f force
+    ////////////////////////////////////////////////////////////////////////////
+    double DerivativeStoppingConditionLoose = 1e-8;
+
+    DenominatorsF.push_back(new FermionActionF(UF,*FGridF,*FrbGridF,*UGrid_f,*GridRBPtrF,light_den[h],M5,b,c, ParamsF));
+    LinOpD.push_back(new LinearOperatorD(*Denominators[h]));
+    LinOpF.push_back(new LinearOperatorF(*DenominatorsF[h]));
+
+    double conv  = DerivativeStoppingCondition;
+    if (h<3) conv= DerivativeStoppingConditionLoose; // Relax on first two hasenbusch factors
+    MPCG.push_back(new MxPCG(conv,
+			     MX_inner,
+			     MaxCGIterations,
+			     UGrid_f,
+			     FrbGridF,
+			     *DenominatorsF[h],*Denominators[h],
+			     *LinOpF[h], *LinOpD[h]) );
+
+    ActionMPCG.push_back(new MxPCG(ActionStoppingCondition,
+				   MX_inner,
+				   MaxCGIterations,
+				   UGrid_f,
+				   FrbGridF,
+				   *DenominatorsF[h],*Denominators[h],
+				   *LinOpF[h], *LinOpD[h]) );
+
+    // Heatbath not mixed yet. As inverts numerators not so important as raised mass.
+    Quotients.push_back (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*ActionMPCG[h],ActionCG));
+#else
+    ////////////////////////////////////////////////////////////////////////////
+    // Standard CG for 2f force
+    ////////////////////////////////////////////////////////////////////////////
+    Quotients.push_back   (new TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],DerivativeCG,ActionCG));
+#endif
+
+  }
+
+  for(int h=0;h<n_hasenbusch+1;h++){
+    Level1.push_back(Quotients[h]);
+  }
+
+  /////////////////////////////////////////////////////////////
+  // Gauge action
+  /////////////////////////////////////////////////////////////
+  Level2.push_back(&GaugeAction);
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  std::cout << GridLogMessage << " Action complete "<< std::endl;
+
+  /////////////////////////////////////////////////////////////
+  // HMC parameters are serialisable
+
+  NoSmearing<HMCWrapper::ImplPolicy> S;
+#ifndef DO_IMPLICIT
+  TrivialMetric<HMCWrapper::ImplPolicy::Field> Mtr;
+#else
+    LaplacianRatParams gpar(2),mpar(2);
+    gpar.offset = 1.;
+    gpar.a0[0] = 500.;
+    gpar.a1[0] = 0.;
+    gpar.b0[0] = 0.25;
+    gpar.b1[0] = 1.;
+    gpar.a0[1] = -500.;
+    gpar.a1[1] = 0.;
+    gpar.b0[1] = 0.36;
+    gpar.b1[1] = 1.2;
+    gpar.b2=1.;
+
+    mpar.offset = 1.;
+    mpar.a0[0] =  -0.850891906532;
+    mpar.a1[0] = -1.54707654538;
+    mpar. b0[0] = 2.85557166137;
+    mpar. b1[0] = 5.74194794773;
+    mpar.a0[1] = -13.5120056831218384729709214298;
+    mpar.a1[1] = 1.54707654538396877086370295729;
+    mpar.b0[1] = 19.2921090880640520026645390317;
+    mpar.b1[1] = -3.54194794773029020262811172870;
+    mpar.b2=1.;
+    for(int i=0;i<2;i++){
+       gpar.a1[i] *=16.;
+       gpar.b1[i] *=16.;
+       mpar.a1[i] *=16.;
+       mpar.b1[i] *=16.;
+    }
+    gpar.b2 *= 16.*16.;
+    mpar.b2 *= 16.*16.;
+
+    ConjugateGradient<LatticeGaugeField> CG(1.0e-8,10000);
+    LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64);
+
+    std::cout << GridLogMessage << "LaplacianRat " << std::endl;
+    gpar.tolerance=HMCparams.MD.RMHMCCGTol;
+    mpar.tolerance=HMCparams.MD.RMHMCCGTol;
+    std::cout << GridLogMessage << "gpar offset= " << gpar.offset <<std::endl;
+    std::cout << GridLogMessage << " a0= " << gpar.a0 <<std::endl;
+    std::cout << GridLogMessage << " a1= " << gpar.a1 <<std::endl;
+    std::cout << GridLogMessage << " b0= " << gpar.b0 <<std::endl;
+    std::cout << GridLogMessage << " b1= " << gpar.b1 <<std::endl;
+    std::cout << GridLogMessage << " b2= " << gpar.b2 <<std::endl ;;
+
+    std::cout << GridLogMessage << "mpar offset= " << mpar.offset <<std::endl;
+    std::cout << GridLogMessage << " a0= " << mpar.a0 <<std::endl;
+    std::cout << GridLogMessage << " a1= " << mpar.a1 <<std::endl;
+    std::cout << GridLogMessage << " b0= " << mpar.b0 <<std::endl;
+    std::cout << GridLogMessage << " b1= " << mpar.b1 <<std::endl;
+    std::cout << GridLogMessage << " b2= " << mpar.b2 <<std::endl;
+//  Assumes PeriodicGimplR or D at the moment
+    auto UGrid = TheHMC.Resources.GetCartesian("gauge");
+//    auto UGrid_f   = GridPtrF;
+//  auto GridPtrF   = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
+//    std::cout << GridLogMessage << " UGrid= " << UGrid <<std::endl;
+//    std::cout << GridLogMessage << " UGrid_f= " << UGrid_f <<std::endl;
+
+    LaplacianAdjointRat<HMCWrapper::ImplPolicy, PeriodicGimplF> Mtr(UGrid, UGrid_f ,CG, gpar, mpar);
+#endif
+
+  std::cout << GridLogMessage << " Running the HMC "<< std::endl;
+  TheHMC.Run(S,Mtr);  // no smearing
+
+  Grid_finalize();
+} // main
+
+
+
--- a/benchmarks/Benchmark_ITT.cc
+++ b/benchmarks/Benchmark_ITT.cc
@@ -365,9 +365,15 @@ public:
    GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
    std::cout << GridLogMessage << "Initialised RNGs" << std::endl;

+#if 1
    typedef DomainWallFermionF Action;
    typedef typename Action::FermionField Fermion;
    typedef LatticeGaugeFieldF Gauge;
+#else
+    typedef GparityDomainWallFermionF Action;
+    typedef typename Action::FermionField Fermion;
+    typedef LatticeGaugeFieldF Gauge;
+#endif
    
    ///////// Source preparation ////////////
    Gauge Umu(UGrid);  SU<Nc>::HotConfiguration(RNG4,Umu); 
@@ -635,6 +641,170 @@ public:
    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
    return mflops_best;
  }
+
+  static double Laplace(int L)
+  {
+    double mflops;
+    double mflops_best = 0;
+    double mflops_worst= 0;
+    std::vector<double> mflops_all;
+
+    ///////////////////////////////////////////////////////
+    // Set/Get the layout & grid size
+    ///////////////////////////////////////////////////////
+    int threads = GridThread::GetThreads();
+    Coordinate mpi = GridDefaultMpi(); assert(mpi.size()==4);
+    Coordinate local({L,L,L,L});
+    Coordinate latt4({local[0]*mpi[0],local[1]*mpi[1],local[2]*mpi[2],local[3]*mpi[3]});
+    
+    GridCartesian         * TmpGrid   = SpaceTimeGrid::makeFourDimGrid(latt4,
+								       GridDefaultSimd(Nd,vComplex::Nsimd()),
+								       GridDefaultMpi());
+    uint64_t NP = TmpGrid->RankCount();
+    uint64_t NN = TmpGrid->NodeCount();
+    NN_global=NN;
+    uint64_t SHM=NP/NN;
+
+
+    ///////// Welcome message ////////////
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+    std::cout<<GridLogMessage << "Benchmark Laplace on "<<L<<"^4 local volume "<<std::endl;
+    std::cout<<GridLogMessage << "* Global volume  : "<<GridCmdVectorIntToString(latt4)<<std::endl;
+    std::cout<<GridLogMessage << "* ranks          : "<<NP  <<std::endl;
+    std::cout<<GridLogMessage << "* nodes          : "<<NN  <<std::endl;
+    std::cout<<GridLogMessage << "* ranks/node     : "<<SHM <<std::endl;
+    std::cout<<GridLogMessage << "* ranks geom     : "<<GridCmdVectorIntToString(mpi)<<std::endl;
+    std::cout<<GridLogMessage << "* Using "<<threads<<" threads"<<std::endl;
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+
+    ///////// Lattice Init ////////////
+    GridCartesian         * FGrid   = SpaceTimeGrid::makeFourDimGrid(latt4, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
+    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid);
+    
+    ///////// RNG Init ////////////
+    std::vector<int> seeds4({1,2,3,4});
+    GridParallelRNG          RNG4(FGrid);  RNG4.SeedFixedIntegers(seeds4);
+    std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
+
+    RealD mass=0.1;
+    RealD c1=9.0/8.0;
+    RealD c2=-1.0/24.0;
+    RealD u0=1.0;
+
+//    typedef ImprovedStaggeredFermionF Action;
+//    typedef typename Action::FermionField Fermion; 
+    typedef LatticeGaugeFieldF Gauge;
+    
+    Gauge Umu(FGrid);  SU<Nc>::HotConfiguration(RNG4,Umu); 
+
+//    typename Action::ImplParams params;
+//    Action Ds(Umu,Umu,*FGrid,*FrbGrid,mass,c1,c2,u0,params);
+
+//  PeriodicGimplF
+    typedef typename PeriodicGimplF::LinkField GaugeLinkFieldF;
+
+    ///////// Source preparation ////////////
+    GaugeLinkFieldF src   (FGrid); random(RNG4,src);
+//    GaugeLinkFieldF src_e (FrbGrid);
+//    GaugeLinkFieldF src_o (FrbGrid);
+//    GaugeLinkFieldF r_e   (FrbGrid);
+//    GaugeLinkFieldF r_o   (FrbGrid);
+    GaugeLinkFieldF r_eo  (FGrid);
+  
+    {
+
+ //     pickCheckerboard(Even,src_e,src);
+ //     pickCheckerboard(Odd,src_o,src);
+    
+      const int num_cases = 1;
+      std::string fmt("G/O/C  ");
+      
+      controls Cases [] = {
+	{  StaggeredKernelsStatic::OptGeneric   ,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  },
+      }; 
+
+      for(int c=0;c<num_cases;c++) {
+        CovariantAdjointLaplacianStencil<PeriodicGimplF,typename PeriodicGimplF::LinkField> LapStencilF(FGrid);
+        QuadLinearOperator<CovariantAdjointLaplacianStencil<PeriodicGimplF,typename PeriodicGimplF::LinkField>,PeriodicGimplF::LinkField> QuadOpF(LapStencilF,c2,c1,1.);
+        LapStencilF.GaugeImport(Umu);
+	
+
+	StaggeredKernelsStatic::Comms = Cases[c].CommsOverlap;
+	StaggeredKernelsStatic::Opt   = Cases[c].Opt;
+	CartesianCommunicator::SetCommunicatorPolicy(Cases[c].CommsAsynch);
+      
+	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+	if ( StaggeredKernelsStatic::Opt == StaggeredKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using Stencil Nc Laplace" <<std::endl;
+	if ( StaggeredKernelsStatic::Comms == StaggeredKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
+	if ( StaggeredKernelsStatic::Comms == StaggeredKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential Comms/Compute" <<std::endl;
+	std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
+	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+	
+	int nwarm = 10;
+	double t0=usecond();
+	FGrid->Barrier();
+	for(int i=0;i<nwarm;i++){
+//	  Ds.DhopEO(src_o,r_e,DaggerNo);
+          QuadOpF.HermOp(src,r_eo);
+	}
+	FGrid->Barrier();
+	double t1=usecond();
+	uint64_t ncall = 500;
+
+	FGrid->Broadcast(0,&ncall,sizeof(ncall));
+
+	//	std::cout << GridLogMessage << " Estimate " << ncall << " calls per second"<<std::endl;
+
+	time_statistics timestat;
+	std::vector<double> t_time(ncall);
+	for(uint64_t i=0;i<ncall;i++){
+	  t0=usecond();
+//	  Ds.DhopEO(src_o,r_e,DaggerNo);
+          QuadOpF.HermOp(src,r_eo);
+	  t1=usecond();
+	  t_time[i] = t1-t0;
+	}
+	FGrid->Barrier();
+	
+	double volume=1;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
+//	double flops=(1146.0*volume)/2;
+	double flops=(2*2*8*216.0*volume);
+	double mf_hi, mf_lo, mf_err;
+	
+	timestat.statistics(t_time);
+	mf_hi = flops/timestat.min;
+	mf_lo = flops/timestat.max;
+	mf_err= flops/timestat.min * timestat.err/timestat.mean;
+
+	mflops = flops/timestat.mean;
+	mflops_all.push_back(mflops);
+	if ( mflops_best == 0   ) mflops_best = mflops;
+	if ( mflops_worst== 0   ) mflops_worst= mflops;
+	if ( mflops>mflops_best ) mflops_best = mflops;
+	if ( mflops<mflops_worst) mflops_worst= mflops;
+	
+	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Quad mflop/s =   "<< mflops << " ("<<mf_err<<") " << mf_lo<<"-"<<mf_hi <<std::endl;
+	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Quad mflop/s per rank   "<< mflops/NP<<std::endl;
+	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Quad mflop/s per node   "<< mflops/NN<<std::endl;
+	FGrid->Barrier();
+      
+      }
+
+      std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+      std::cout<<GridLogMessage << L<<"^4  Quad Best  mflop/s        =   "<< mflops_best << " ; " << mflops_best/NN<<" per node " <<std::endl;
+      std::cout<<GridLogMessage << L<<"^4  Quad Worst mflop/s        =   "<< mflops_worst<< " ; " << mflops_worst/NN<<" per node " <<std::endl;
+      std::cout<<GridLogMessage <<fmt << std::endl;
+      std::cout<<GridLogMessage ;
+	FGrid->Barrier();
+
+      for(int i=0;i<mflops_all.size();i++){
+	std::cout<<mflops_all[i]/NN<<" ; " ;
+      }
+      std::cout<<std::endl;
+    }
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+    return mflops_best;
+  }
 };


@@ -662,6 +832,7 @@ int main (int argc, char ** argv)
  std::vector<double> wilson;
  std::vector<double> dwf4;
  std::vector<double> staggered;
+  std::vector<double> lap;

  int Ls=1;
  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
@@ -688,12 +859,20 @@ int main (int argc, char ** argv)
    staggered.push_back(result);
  }

+  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << " Laplace QuadOp 4D " <<std::endl;
+  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+  for(int l=0;l<L_list.size();l++){
+    double result = Benchmark::Laplace(L_list[l]) ;
+    lap.push_back(result);
+  }
+
  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
  std::cout<<GridLogMessage << " Summary table Ls="<<Ls <<std::endl;
  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << "L \t\t Wilson \t\t DWF4 \t\t Staggered" <<std::endl;
+  std::cout<<GridLogMessage << "L \t\t Wilson \t\t DWF4 \t\t Staggered \t\t Quad Laplace" <<std::endl;
  for(int l=0;l<L_list.size();l++){
-    std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< wilson[l]<<" \t\t "<<dwf4[l] << " \t\t "<< staggered[l]<<std::endl;
+    std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< wilson[l]<<" \t\t "<<dwf4[l] << " \t\t "<< staggered[l]<< " \t\t "<< lap[l]<< std::endl;
  }
  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;

--- a/tests/hmc/Test_hmc_WilsonGauge.cc
+++ b/tests/hmc/Test_hmc_WilsonGauge.cc
@@ -83,8 +83,15 @@ int main(int argc, char **argv)
  // need wrappers of the fermionic classes 
  // that have a complex construction
  // standard
-  RealD beta = 5.6 ;
+  RealD beta = 6.6 ; 
+
+#if 0
  WilsonGaugeActionR Waction(beta);
+#else
+  std::vector<Complex> boundaryG = {1,1,1,0};
+  WilsonGaugeActionR::ImplParams ParamsG(boundaryG);
+  WilsonGaugeActionR Waction(beta,ParamsG);
+#endif
  
  ActionLevel<HMCWrapper::Field> Level1(1);
  Level1.push_back(&Waction);
--- a/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc
+++ b/tests/hmc/Test_hmc_WilsonGauge_Implicit.cc
@@ -0,0 +1,238 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+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>
+
+#undef USE_OBC
+#define DO_IMPLICIT
+
+
+int main(int argc, char **argv) 
+{
+  using namespace Grid;
+
+  Grid_init(&argc, &argv);
+  GridLogLayout();
+  
+  std::string arg;
+ 
+  HMCparameters HMCparams;
+#if 1
+  {
+    XmlReader  HMCrd("HMCparameters.xml");
+    read(HMCrd,"HMCparameters",HMCparams);
+  }
+#else
+//IntegratorParameters MD;
+  std::vector<int> steps(0);
+  if( GridCmdOptionExists(argv,argv+argc,"--MDsteps") ){
+    arg= GridCmdOptionPayload(argv,argv+argc,"--MDsteps");
+    GridCmdOptionIntVector(arg,steps);
+    assert(steps.size()==1);
+  }
+  MD.trajL   = 0.001*std::sqrt(2.);
+  MD.MDsteps = 1;
+  if (steps.size()>0) MD.MDsteps = steps[0];
+  if( GridCmdOptionExists(argv,argv+argc,"--trajL") ){
+    arg= GridCmdOptionPayload(argv,argv+argc,"--trajL");
+    std::vector<int> traj(0);
+    GridCmdOptionIntVector(arg,traj);
+    assert(traj.size()==1);
+    MD.trajL *= double(traj[0]);
+  }
+  MD.RMHMCTol=1e-8;
+  MD.RMHMCCGTol=1e-8;
+  std::cout << "RMHMCTol= "<<  MD.RMHMCTol<<" RMHMCCGTol= "<<MD.RMHMCCGTol<<std::endl;
+
+  HMCparameters HMCparams;
+  HMCparams.StartTrajectory  = 0;
+  HMCparams.Trajectories     = 1;
+  HMCparams.NoMetropolisUntil=  100;
+  // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
+  HMCparams.StartingType     =std::string("ColdStart");
+  HMCparams.Kappa=0.01; //checking against trivial. Pathetic.
+  HMCparams.MD = MD;
+#endif
+
+
+
+   // Typedefs to simplify notation
+#ifdef DO_IMPLICIT
+  typedef GenericHMCRunner<ImplicitMinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+//  typedef GenericHMCRunner<ImplicitCampostrini> HMCWrapper;  // 4th order
+  HMCparams.MD.name    = std::string("ImplicitMinimumNorm2");
+#else
+  typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+  HMCparams.MD.name    = std::string("MinimumNorm2");
+#endif
+
+
+
+  // Possibile to create the module by hand 
+  // hardcoding parameters or using a Reader
+
+
+  // Checkpointer definition
+  CheckpointerParameters CPparams;  
+  CPparams.config_prefix = "ckpoint_lat";
+  CPparams.rng_prefix = "ckpoint_rng";
+  CPparams.saveInterval = 1;
+  CPparams.format = "IEEE64BIG";
+  
+  HMCWrapper TheHMC(HMCparams);
+  // Grid from the command line
+  TheHMC.Resources.AddFourDimGrid("gauge");
+  TheHMC.Resources.LoadNerscCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.serial_seeds = "1 2 3 4 5";
+  RNGpar.parallel_seeds = "6 7 8 9 10";
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  // here there is too much indirection 
+  typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
+  typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
+  TheHMC.Resources.AddObservable<PlaqObs>();
+  TopologyObsParameters TopParams;
+  TopParams.interval = 1;
+  TopParams.do_smearing = true;
+//  TopParams.Smearing.steps = 1600;
+//  TopParams.Smearing.step_size = 0.01;
+  TopParams.Smearing.init_step_size = 0.01;
+  TopParams.Smearing.meas_interval = 10;
+  TopParams.Smearing.maxTau = 16.0; 
+//  TheHMC.Resources.AddObservable<QObs>(TopParams);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes 
+  // that have a complex construction
+  // standard
+
+  RealD beta = 6.6;
+  std::cout << "Wilson Gauge beta= " <<beta <<std::endl;
+#ifndef USE_OBC
+  WilsonGaugeActionR Waction(beta);
+#else
+  std::vector<Complex> boundaryG = {1,1,1,0};
+  WilsonGaugeActionR::ImplParams ParamsG(boundaryG);
+  WilsonGaugeActionR Waction(beta,ParamsG);
+  std::cout << "boundaryG = " <<boundaryG  <<std::endl;
+#endif
+
+  
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Waction);
+  TheHMC.TheAction.push_back(Level1);
+
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  std::cout << "trajL= " <<TheHMC.Parameters.MD.trajL <<" steps= "<<TheHMC.Parameters.MD.MDsteps << " integrator= "<<TheHMC.Parameters.MD.name<<std::endl;
+
+  NoSmearing<HMCWrapper::ImplPolicy> S;
+#ifndef DO_IMPLICIT
+  TrivialMetric<HMCWrapper::ImplPolicy::Field> Mtr;
+#else
+// g_x3_2
+    LaplacianRatParams gpar(2),mpar(2);
+    gpar.offset = 1.;
+    gpar.a0[0] = 500.;
+    gpar.a1[0] = 0.;
+    gpar.b0[0] = 0.25;
+    gpar.b1[0] = 1.;
+    gpar.a0[1] = -500.;
+    gpar.a1[1] = 0.;
+    gpar.b0[1] = 0.36;
+    gpar.b1[1] = 1.2;
+    gpar.b2=1.;
+
+    mpar.offset = 1.;
+    mpar.a0[0] =  -0.850891906532;
+    mpar.a1[0] = -1.54707654538;
+    mpar. b0[0] = 2.85557166137;
+    mpar. b1[0] = 5.74194794773;
+    mpar.a0[1] = -13.5120056831218384729709214298;
+    mpar.a1[1] = 1.54707654538396877086370295729;
+    mpar.b0[1] = 19.2921090880640520026645390317;
+    mpar.b1[1] = -3.54194794773029020262811172870;
+    mpar.b2=1.;
+    for(int i=0;i<2;i++){
+       gpar.a1[i] *=16.;
+       gpar.b1[i] *=16.;
+       mpar.a1[i] *=16.;
+       mpar.b1[i] *=16.;
+    }
+    gpar.b2 *= 16.*16.;
+    mpar.b2 *= 16.*16.;
+
+    ConjugateGradient<LatticeGaugeField> CG(1.0e-8,10000);
+    LaplacianParams LapPar(0.0001, 1.0, 10000, 1e-8, 12, 64);
+
+    std::cout << GridLogMessage << "LaplacianRat " << std::endl;
+
+    gpar.tolerance=HMCparams.MD.RMHMCCGTol;
+    mpar.tolerance=HMCparams.MD.RMHMCCGTol;
+
+    std::cout << GridLogMessage << "gpar offset= " << gpar.offset <<std::endl;
+    std::cout << GridLogMessage << " a0= " << gpar.a0 <<std::endl;
+    std::cout << GridLogMessage << " a1= " << gpar.a1 <<std::endl;
+    std::cout << GridLogMessage << " b0= " << gpar.b0 <<std::endl;
+    std::cout << GridLogMessage << " b1= " << gpar.b1 <<std::endl;
+    std::cout << GridLogMessage << " b2= " << gpar.b2 <<std::endl ;;
+
+    std::cout << GridLogMessage << "mpar offset= " << mpar.offset <<std::endl;
+    std::cout << GridLogMessage << " a0= " << mpar.a0 <<std::endl;
+    std::cout << GridLogMessage << " a1= " << mpar.a1 <<std::endl;
+    std::cout << GridLogMessage << " b0= " << mpar.b0 <<std::endl;
+    std::cout << GridLogMessage << " b1= " << mpar.b1 <<std::endl;
+    std::cout << GridLogMessage << " b2= " << mpar.b2 <<std::endl;
+//  Assumes PeriodicGimplR or D at the moment
+    Coordinate latt  = GridDefaultLatt();
+    Coordinate mpi   = GridDefaultMpi();
+    auto UGrid = TheHMC.Resources.GetCartesian("gauge");
+    Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd());
+    auto UGrid_f   = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
+    std::cout << GridLogMessage << " UGrid= " << UGrid <<std::endl;
+    std::cout << GridLogMessage << " UGrid_f= " << UGrid_f <<std::endl;
+
+    LaplacianAdjointRat<HMCWrapper::ImplPolicy, PeriodicGimplF> Mtr(UGrid, UGrid_f,CG, gpar, mpar);
+#endif
+ 
+  {
+    XmlWriter HMCwr("HMCparameters.xml.out");
+    write(HMCwr,"HMCparameters",TheHMC.Parameters);
+  }
+
+  TheHMC.Run(S,Mtr);  // no smearing
+
+  Grid_finalize();
+
+} // main
Author	SHA1	Message	Date
Chulwoo Jung	cfa0576ffd	Getting rid of one more non-auto View, comms overlap in Laplace operator	2024-02-25 22:37:48 -05:00
Chulwoo Jung	fe98e9f555	Fixing Laplace flopcount Minor cleanup	2024-02-13 12:06:08 -05:00
Chulwoo Jung	948d16fb06	Laplace benchmark added	2024-02-12 21:23:36 -05:00
Chulwoo Jung	58fbcaa399	Checking in before cleaning up	2024-02-12 21:10:21 -05:00
Chulwoo Jung	9ad6836b0f	Mixed precision for Laplace. Main program with Metric	2024-02-08 17:13:10 -05:00
Chulwoo Jung	026eb8a695	Wilson RMHMC main program	2023-12-12 15:34:03 -05:00
Chulwoo Jung	076580c232	Recovering mixed precision CG for Laplace Checking in to move to aurora	2023-12-12 15:32:00 -05:00
Chulwoo Jung	7af6022a2a	Added midMD checkpointing (for lattice only for now)	2023-12-04 20:05:41 -05:00
Chulwoo Jung	982a60536c	Checking in before forking	2023-11-22 16:33:15 -05:00
Chulwoo Jung	dc36d272ce	Gauge RMHMC conserving dH	2023-11-21 13:48:51 -05:00
Chulwoo Jung	515ff6bf62	Added Laplacian metric, Gauge OpenBC	2023-11-09 21:42:46 -05:00