Covariant laplacian and implicit integration

lib/qcd/hmc/HMCModules.h

@@ -37,25 +37,19 @@ namespace Grid {
 namespace QCD {
 
 ////////////////////////////////////////////////////////////////////
-class RNGModuleParameters: Serializable {
+struct RNGModuleParameters: Serializable {
-
   GRID_SERIALIZABLE_CLASS_MEMBERS(RNGModuleParameters,
   std::string, serial_seeds,
   std::string, parallel_seeds,);
-public:
-  std::vector<int> SerialSeed;
-  std::vector<int> ParallelSeed;
 
-  RNGModuleParameters(const std::vector<int> S = std::vector<int>(),
+  std::vector<int> getSerialSeeds(){return strToVec<int>(serial_seeds);}
-                      const std::vector<int> P = std::vector<int>())
+  std::vector<int> getParallelSeeds(){return strToVec<int>(parallel_seeds);}
-      : SerialSeed(S), ParallelSeed(P) {}
 
+  RNGModuleParameters(): serial_seeds("1"), parallel_seeds("1"){}
 
   template <class ReaderClass >
   RNGModuleParameters(Reader<ReaderClass>& Reader){
     read(Reader, "RandomNumberGenerator", *this); 
-    SerialSeed = strToVec<int>(serial_seeds);
-    ParallelSeed = strToVec<int>(parallel_seeds);
   }
   
 };
@@ -82,12 +76,14 @@ public:
   GridParallelRNG& get_pRNG() { return *pRNG_.get(); }
 
   void seed() {
-    if (Params_.SerialSeed.size() == 0 && Params_.ParallelSeed.size() == 0) {
+    auto SerialSeeds   = Params_.getSerialSeeds();
-      std::cout << "Seeds not initialized" << std::endl;
+    auto ParallelSeeds = Params_.getParallelSeeds();
+    if (SerialSeeds.size() == 0 && ParallelSeeds.size() == 0) {
+      std::cout << GridLogError << "Seeds not initialized" << std::endl;
       exit(1);
     }
-    sRNG_.SeedFixedIntegers(Params_.SerialSeed);
+    sRNG_.SeedFixedIntegers(SerialSeeds);
-    pRNG_->SeedFixedIntegers(Params_.ParallelSeed);
+    pRNG_->SeedFixedIntegers(ParallelSeeds);
   }
 };
 

lib/qcd/hmc/HMCResourceManager.h

@@ -75,6 +75,8 @@ class HMCResourceManager {
   bool have_RNG;
   bool have_CheckPointer;
 
+  // NOTE: operator << is not overloaded for std::vector<string> 
+  // so thsi function is necessary
   void output_vector_string(const std::vector<std::string> &vs){
     for (auto &i: vs)
       std::cout << i << " ";
@@ -85,13 +87,13 @@ class HMCResourceManager {
  public:
   HMCResourceManager() : have_RNG(false), have_CheckPointer(false) {}
 
-  template <class ReaderClass >
+  template <class ReaderClass, class vector_type = vComplex >
   void initialize(ReaderClass &Read){
     // assumes we are starting from the main node
 
     // Geometry
     GridModuleParameters GridPar(Read);
-    GridFourDimModule GridMod( GridPar) ;
+    GridFourDimModule<vector_type> GridMod( GridPar) ;
     AddGrid("gauge", GridMod);
 
     // Checkpointer
@@ -100,9 +102,6 @@ class HMCResourceManager {
     std::string cp_type;
     read(Read,"name", cp_type);
     std::cout << "Registered types " << std::endl;
-    // NOTE: operator << is not overloaded for std::vector<string> 
-    // so it complains here
-    //std::cout << CPfactory.getBuilderList() << std::endl;
     output_vector_string(CPfactory.getBuilderList());
 
 
@@ -178,7 +177,7 @@ class HMCResourceManager {
 
   // Add a named grid set, 4d shortcut
   void AddFourDimGrid(std::string s) {
-    GridFourDimModule Mod;
+    GridFourDimModule<vComplex> Mod;
     AddGrid(s, Mod);
   }
 

lib/qcd/hmc/HMC_GridModules.h

@@ -31,41 +31,48 @@ directory
 #define HMC_GRID_MODULES
 
 namespace Grid {
-namespace QCD {
 
 // Resources
 // Modules for grids 
 
-class GridModuleParameters: Serializable{
+// Introduce another namespace HMCModules?
-  
+
+class GridModuleParameters: Serializable{   
+public:
   GRID_SERIALIZABLE_CLASS_MEMBERS(GridModuleParameters,
   std::string, lattice,
-  std::string,  mpi);
+  std::string, mpi);
 
-public: 
+  std::vector<int> getLattice(){return strToVec<int>(lattice);}
-  // these namings are ugly
+  std::vector<int> getMpi()    {return strToVec<int>(mpi);}
-  // also ugly the distinction between the serializable members
-  // and this
-  std::vector<int> lattice_v;
-  std::vector<int> mpi_v;
 
-  GridModuleParameters(const std::vector<int> l_ = std::vector<int>(),
+  void check(){
-    const std::vector<int> mpi_ = std::vector<int>()):lattice_v(l_), mpi_v(mpi_){}
+    if (getLattice().size() != getMpi().size()) {
-
+      std::cout << GridLogError 
-  template <class ReaderClass>
+                << "Error in GridModuleParameters: lattice and mpi dimensions "
-  GridModuleParameters(Reader<ReaderClass>& Reader) {
-    read(Reader, "LatticeGrid", *this);
-    lattice_v = strToVec<int>(lattice);
-    mpi_v = strToVec<int>(mpi);
-    if (mpi_v.size() != lattice_v.size()) {
-      std::cout << "Error in GridModuleParameters: lattice and mpi dimensions "
                    "do not match"
                 << std::endl;
       exit(1);
     }
+  }    
+
+  template <class ReaderClass>
+  GridModuleParameters(Reader<ReaderClass>& Reader) {
+    read(Reader, name, *this);
+    check();
   }
+
+  // Save on file
+  template< class WriterClass>
+  void save(Writer<WriterClass>& Writer){
+    check();
+    write(Writer, name, *this);
+  }
+private:
+    std::string name = "LatticeGrid";
 };
 
+// Lower level class
 class GridModule {
  public:
   GridCartesian* get_full() { 
@@ -84,27 +91,33 @@ class GridModule {
   
 };
 
-// helpers
+////////////////////////////////////
-// FIXME define a class accepting also real vtypes
+// Classes for the user
+////////////////////////////////////
+// Note: the space time grid must be out of the QCD namespace
+template< class vector_type>
 class GridFourDimModule : public GridModule {
  public:
-  // add a function to create the module from a Reader
   GridFourDimModule() {
+    using namespace QCD;
     set_full(SpaceTimeGrid::makeFourDimGrid(
-        GridDefaultLatt(), GridDefaultSimd(4, vComplex::Nsimd()),
+        GridDefaultLatt(), GridDefaultSimd(4, vector_type::Nsimd()),
         GridDefaultMpi()));
     set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
   }
 
-  template <class vector_type = vComplex>
   GridFourDimModule(GridModuleParameters Params) {
-    if (Params.lattice_v.size() == 4) {
+    using namespace QCD;
+    Params.check();
+    std::vector<int> lattice_v = Params.getLattice();
+    std::vector<int> mpi_v = Params.getMpi();
+    if (lattice_v.size() == 4) {
       set_full(SpaceTimeGrid::makeFourDimGrid(
-          Params.lattice_v, GridDefaultSimd(4, vector_type::Nsimd()),
+          lattice_v, GridDefaultSimd(4, vector_type::Nsimd()),
-          Params.mpi_v));
+          mpi_v));
       set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(grid_.get()));
     } else {
-      std::cout
+      std::cout << GridLogError 
           << "Error in GridFourDimModule: lattice dimension different from 4"
           << std::endl;
       exit(1);
@@ -112,10 +125,9 @@ class GridFourDimModule : public GridModule {
   }
 };
 
+typedef GridFourDimModule<vComplex> GridDefaultFourDimModule;
 
 
-
-}  // namespace QCD
 }  // namespace Grid
 
 #endif  // HMC_GRID_MODULES

lib/qcd/hmc/integrators/Integrator.h

@@ -131,6 +131,49 @@ class Integrator {
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
   }
 
+  void implicit_update_P(MomentaField& Mom, Field& U, int level, double ep) {
+    // Fundamental updates, include smearing
+    MomentaField Msum(Mom._grid);
+    Msum = zero;
+    for (int a = 0; a < as[level].actions.size(); ++a) {
+      // Compute the force 
+      // We need to compute the derivative of the actions
+      // only once
+      Field force(U._grid);
+      conformable(U._grid, 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 average: " << force_abs << std::endl;
+      Msum += force;
+    }
+
+    MomentaField NewMom = Mom;
+    MomentaField OldMom = Mom;
+    double threshold = 1e-6;
+    // Here run recursively
+    do{
+      MomentaField MomDer(Mom._grid);
+      OldMom = NewMom;
+      // Compute the derivative of the kinetic term
+      // with respect to the gauge field
+
+      // Laplacian.Mder(NewMom, MomDer);
+      // NewMom = Mom - ep*(MomDer + Msum);
+
+    } while (norm2(NewMom - OldMom) > threshold);
+
+    Mom = NewMom;
+
+
+    // update the auxiliary fields momenta
+  }
+
+
   void update_U(Field& U, double ep) {
     update_U(P, U, ep);
 

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -285,6 +285,74 @@ class ForceGradient : public Integrator<FieldImplementation, SmearingPolicy,
     }
   }
 };
+
+
+////////////////////////////////
+// Riemannian Manifold HMC
+// Girolami et al
+////////////////////////////////
+
+
+
+
+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)
+      : Integrator<FieldImplementation, SmearingPolicy, RepresentationPolicy>(
+            grid, Par, Aset, Sm){};
+
+  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->update_U(U, eps);
+      } else {  // recursive function call
+        this->step(U, level + 1, first_step, last_step);
+      }
+
+      int mm = last_step ? 1 : 2;
+      this->update_P(U, level, mm * eps / 2.0);
+    }
+  }
+};
+
+
+
+
+
+
+
+
 }
 }
 

lib/qcd/utils/CovariantLaplacian.h

@@ -33,13 +33,27 @@ directory
 namespace Grid {
 namespace QCD {
 
+////////////////////////////////////////////////////////////
+// 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 LaplacianAdjointField {
  public:
   INHERIT_GIMPL_TYPES(Impl);
-  typedef SU<Nc>::LatticeAlgebraVector AVector;
+  
-
+  LaplacianAdjointField(GridBase* grid, const RealD k = 1.0) : 
-  LaplacianAdjointField(GridBase* grid) : U(Nd, grid){};
+    U(Nd, grid), kappa(k){};
 
   void ImportGauge(const GaugeField& _U) {
     for (int mu = 0; mu < Nd; mu++) {
@@ -49,61 +63,48 @@ class LaplacianAdjointField {
 
   void Mdiag(const GaugeLinkField& in, GaugeLinkField& out) { assert(0); }
 
-  void Mdir(const GaugeLinkField& in, GaugeLinkField& out, int dir, int disp) { assert(0); }
+  void Mdir(const GaugeLinkField& in, GaugeLinkField& out, int dir, int disp) {
-
+    assert(0);
-/*
-  // Operator with algebra vector inputs and outputs
-  void M2(const AVector& in, AVector& out) {
-    double kappa = 0.9;
-    //Reconstruct matrix
-
-    GaugeLinkField tmp(in._grid);
-    GaugeLinkField tmp2(in._grid);
-    GaugeLinkField sum(in._grid);
-    GaugeLinkField out_mat(in._grid);
-    GaugeLinkField in_mat(in._grid);
-    SU<Nc>::FundamentalLieAlgebraMatrix(in, in_mat);
-
-    sum = zero;
-    for (int mu = 0; mu < Nd; mu++) {
-      tmp  = U[mu] * Cshift(in_mat, mu, +1) * adj(U[mu]);
-      tmp2 = adj(U[mu]) * in_mat * U[mu];
-      sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_mat;
-    }
-    out_mat = (1.0 - kappa) * in_mat - kappa/(double(4*Nd)) * sum;
-    // Project
-    SU<Nc>::projectOnAlgebra(out, out_mat);
   }
-*/
-
-    void M(const GaugeLinkField& in, GaugeLinkField& out) {
-    double kappa = 0.999;
-    //Reconstruct matrix
 
+  void M(const GaugeLinkField& in, GaugeLinkField& out) {
     GaugeLinkField tmp(in._grid);
     GaugeLinkField tmp2(in._grid);
     GaugeLinkField sum(in._grid);
     sum = zero;
     for (int mu = 0; mu < Nd; mu++) {
-      tmp  = U[mu] * Cshift(in, mu, +1) * adj(U[mu]);
+      tmp = U[mu] * Cshift(in, mu, +1) * adj(U[mu]);
       tmp2 = adj(U[mu]) * in * U[mu];
       sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in;
     }
-    out = (1.0 - kappa) * in - kappa/(double(4*Nd)) * sum;
+    out = (1.0 - kappa) * in - kappa / (double(4 * Nd)) * sum;
+  }
+
+  void MDeriv(const GaugeLinkField& in, GaugeLinkField& out, bool dag){
+    RealD factor = - kappa / (double(4 * Nd))
+    if (!dag)
+      out = factor * Cshift(in, mu, +1) * adj(U[mu]) + adj(U[mu]) * in;
+    else
+      out = factor * U[mu] * Cshift(in, mu, +1) + in * U[mu];
   }
 
  private:
+  RealD kappa;
   std::vector<GaugeLinkField> U;
 };
 
 
+// This is just a debug tests
+// not meant to be used 
+
 template <class Impl>
 class LaplacianAlgebraField {
  public:
   INHERIT_GIMPL_TYPES(Impl);
   typedef SU<Nc>::LatticeAlgebraVector AVector;
 
-  LaplacianAlgebraField(GridBase* grid) : U(Nd, grid){};
+  LaplacianAlgebraField(GridBase* grid, const RealD k) : 
+    U(Nd, grid), kappa(k){};
 
   void ImportGauge(const GaugeField& _U) {
     for (int mu = 0; mu < Nd; mu++) {
@@ -117,28 +118,28 @@ class LaplacianAlgebraField {
 
   // Operator with algebra vector inputs and outputs
   void M(const AVector& in, AVector& out) {
-    double kappa = 0.999;
-    //Reconstruct matrix
-
     GaugeLinkField tmp(in._grid);
     GaugeLinkField tmp2(in._grid);
     GaugeLinkField sum(in._grid);
     GaugeLinkField out_mat(in._grid);
     GaugeLinkField in_mat(in._grid);
+
+    // Reconstruct matrix
     SU<Nc>::FundamentalLieAlgebraMatrix(in, in_mat);
 
     sum = zero;
     for (int mu = 0; mu < Nd; mu++) {
-      tmp  = U[mu] * Cshift(in_mat, mu, +1) * adj(U[mu]);
+      tmp = U[mu] * Cshift(in_mat, mu, +1) * adj(U[mu]);
       tmp2 = adj(U[mu]) * in_mat * U[mu];
       sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_mat;
     }
-    out_mat = (1.0 - kappa) * in_mat - kappa/(double(4*Nd)) * sum;
+    out_mat = (1.0 - kappa) * in_mat - kappa / (double(4 * Nd)) * sum;
     // Project
     SU<Nc>::projectOnAlgebra(out, out_mat);
   }
 
  private:
+  RealD kappa;
   std::vector<GaugeLinkField> U;
 };
 

tests/hmc/Test_hmc_WilsonGauge.cc

@@ -59,8 +59,8 @@ int main(int argc, char **argv) {
   TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
 
   RNGModuleParameters RNGpar;
-  RNGpar.SerialSeed = {1,2,3,4,5};
+  RNGpar.serial_seeds = "1 2 3 4 5";
-  RNGpar.ParallelSeed = {6,7,8,9,10};
+  RNGpar.parallel_seeds = "6 7 8 9 10";
   TheHMC.Resources.SetRNGSeeds(RNGpar);
 
   // Construct observables

tests/hmc/Test_hmc_WilsonTMFermionGauge.cc

@@ -0,0 +1,147 @@
+/*************************************************************************************
+
+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>
+
+int main(int argc, char **argv) {
+  using namespace Grid;
+  using namespace Grid::QCD;
+
+  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 GenericHMCRunner<MinimumNorm2> HMCWrapper;  // Uses the default minimum norm
+  typedef WilsonImplR FermionImplPolicy;
+  typedef WilsonTMFermionR FermionAction;
+  typedef typename FermionAction::FermionField FermionField;
+
+
+  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+  HMCWrapper TheHMC;
+
+  // Grid from the command line
+  TheHMC.Resources.AddFourDimGrid("gauge");
+  // 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 = 5;
+  CPparams.format = "IEEE64BIG";
+  
+  TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
+
+  RNGModuleParameters RNGpar;
+  RNGpar.SerialSeed = {1,2,3,4,5};
+  RNGpar.ParallelSeed = {6,7,8,9,10};
+  TheHMC.Resources.SetRNGSeeds(RNGpar);
+
+  // Construct observables
+  // here there is too much indirection 
+  PlaquetteObsParameters PlPar;
+  PlPar.output_prefix = "Plaquette";
+  PlaquetteMod<HMCWrapper::ImplPolicy> PlaqModule(PlPar);
+  TheHMC.Resources.AddObservable(&PlaqModule);
+  //////////////////////////////////////////////
+
+  /////////////////////////////////////////////////////////////
+  // Collect actions, here use more encapsulation
+  // need wrappers of the fermionic classes 
+  // that have a complex construction
+  // standard
+  RealD beta = 3.9 ;
+  SymanzikGaugeActionR Waction(beta);
+    
+  auto GridPtr = TheHMC.Resources.GetCartesian();
+  auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
+
+  // temporarily need a gauge field
+  LatticeGaugeField U(GridPtr);
+
+  Real mass = -0.89163;
+  Real mu   = 0.01;
+
+  // Can we define an overloaded operator that does not need U and initialises
+  // it with zeroes?
+  FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, mu);
+
+  ConjugateGradient<FermionField> CG(1.0e-8, 2000);
+
+  TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
+
+  // With modules
+  /*
+
+  TwoFlavourFmodule<FermionImplPolicy> TwoFMod(Reader);
+  
+  */
+
+    // Set smearing (true/false), default: false
+  Nf2.is_smeared = false;
+
+
+    // Collect actions
+  ActionLevel<HMCWrapper::Field> Level1(1);
+  Level1.push_back(&Nf2);
+
+  ActionLevel<HMCWrapper::Field> Level2(4);
+  Level2.push_back(&Waction);
+
+  TheHMC.TheAction.push_back(Level1);
+  TheHMC.TheAction.push_back(Level2);
+  /////////////////////////////////////////////////////////////
+
+  /*
+    double rho = 0.1;  // smearing parameter
+    int Nsmear = 2;    // number of smearing levels
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(
+        UGrid, Nsmear, Stout);
+  */
+
+  // HMC parameters are serialisable 
+  TheHMC.Parameters.MD.MDsteps = 20;
+  TheHMC.Parameters.MD.trajL   = 1.0;
+
+  TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
+  TheHMC.Run();  // no smearing
+  // TheHMC.Run(SmearingPolicy); // for smearing
+
+  Grid_finalize();
+
+} // main
+
+

tests/solver/Test_laplacian.cc

@@ -48,28 +48,33 @@ int main (int argc, char ** argv)
   LatticeGaugeField Umu(&Grid); 
   SU<Nc>::HotConfiguration(pRNG,Umu);
 
+  double Kappa = 0.9999;
+
   typedef SU<Nc>::LatticeAlgebraVector AVector;
   // Source and result in the algebra
   AVector src_vec(&Grid); random(pRNG, src_vec);
   AVector result_vec(&Grid); result_vec = zero;
+  
   LatticeColourMatrix src(&Grid); 
   SU<Nc>::FundamentalLieAlgebraMatrix(src_vec, src);
   LatticeColourMatrix result(&Grid); result=zero;
 
-  LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid);
+  LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid, Kappa);
   Laplacian.ImportGauge(Umu);
 
   HermitianLinearOperator<LaplacianAdjointField<PeriodicGimplR>,LatticeColourMatrix> HermOp(Laplacian);
   ConjugateGradient<LatticeColourMatrix> CG(1.0e-8,10000);
+  std::cout << GridLogMessage << "Testing the Laplacian using the full matrix" <<std::endl;
   CG(HermOp,src,result); // fastest
   
 
   // Tests also the version using the algebra decomposition
-  LaplacianAlgebraField<PeriodicGimplR> LaplacianAlgebra(&Grid);
+  LaplacianAlgebraField<PeriodicGimplR> LaplacianAlgebra(&Grid, Kappa);
   LaplacianAlgebra.ImportGauge(Umu);
 
   HermitianLinearOperator<LaplacianAlgebraField<PeriodicGimplR>,AVector> HermOpAlg(LaplacianAlgebra);
   ConjugateGradient<AVector> CG_Algebra(1.0e-8,10000);
+  std::cout << GridLogMessage << "Testing the Laplacian using the algebra vectors" <<std::endl;
   CG_Algebra(HermOpAlg,src_vec,result_vec);
   
   LatticeColourMatrix result2(&Grid);