Added and tested the covariant laplacian + CG solver

lib/qcd/action/Actions.h

@@ -48,9 +48,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 ////////////////////////////////////////////
 // Utility functions
 ////////////////////////////////////////////
-//#include <Grid/qcd/action/gauge/GaugeImplTypes.h>
 #include <Grid/qcd/action/gauge/GaugeImplementations.h>
 #include <Grid/qcd/utils/WilsonLoops.h>
+#include <Grid/qcd/utils/CovariantLaplacian.h>
 
 #include <Grid/qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions
 #include <Grid/qcd/action/fermion/FermionOperatorImpl.h>

lib/qcd/utils/CovariantLaplacian.h

@@ -0,0 +1,82 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h
+
+Copyright (C) 2016
+
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef COVARIANT_LAPLACIAN_H
+#define COVARIANT_LAPLACIAN_H
+
+namespace Grid {
+namespace QCD {
+
+template <class Impl>
+class LaplacianAdjointField {
+ public:
+  INHERIT_GIMPL_TYPES(Impl);
+  typedef SU<Nc>::LatticeAlgebraVector AVector;
+
+  LaplacianAdjointField(GridBase* grid) : U(Nd, grid){};
+
+  void ImportGauge(const GaugeField& _U) {
+    for (int mu = 0; mu < Nd; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(_U, mu);
+    }
+  }
+
+  void Mdiag(const AVector& in, AVector& out) { assert(0); }
+
+  void Mdir(const AVector& in, AVector& out, int dir, int disp) { assert(0); }
+
+  void M(const AVector& in, AVector& out) {
+    double kappa = 0.99;
+    //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);
+  }
+
+ private:
+  std::vector<GaugeLinkField> U;
+};
+}
+}
+
+#endif

tests/solver/Test_laplacian.cc

@@ -0,0 +1,64 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_wilson_cg_prec.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    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>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({1,2,3,4,5});
+  GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);
+
+  LatticeGaugeField Umu(&Grid); 
+  SU<Nc>::HotConfiguration(pRNG,Umu);
+
+  typedef SU<Nc>::LatticeAlgebraVector AVector;
+  // Source and result in the algebra
+  AVector src(&Grid); random(pRNG,src);
+  AVector result(&Grid); result=zero;
+
+  LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid);
+  Laplacian.ImportGauge(Umu);
+
+  HermitianLinearOperator<LaplacianAdjointField<PeriodicGimplR>,AVector> HermOp(Laplacian);
+  ConjugateGradient<AVector> CG(1.0e-8,10000);
+  CG(HermOp,src,result);
+  
+  Grid_finalize();
+}