Added laplacian operator for smearing sources

tests/core/Test_laplacian.cc

@@ -0,0 +1,105 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/Test_laplacian.cc
+
+    Copyright (C) 2017
+
+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>
+
+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(&Grid);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+
+  std::vector<int> point({0,0,0,0});
+
+  LatticeFermion src   (&Grid); //random(pRNG,src);
+  SpinColourVectorD Sp;
+  for (unsigned int s = 0; s < Ns; ++s)
+      for (unsigned int c = 0; c < Nc; ++c)
+        Sp()(s)(c) = 1;
+
+  src = zero;
+  pokeSite(Sp,src,point);
+
+  LatticeFermion result(&Grid); result=zero;
+  LatticeFermion tmp(&Grid); tmp=zero;
+
+  // Gauge configuration
+  LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Testing the laplacian operator on a point source            "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  Laplacian<WilsonImplR> LaplaceOperator(src._grid);
+  LaplaceOperator.ImportGauge(Umu);
+  LaplaceOperator.M(src, result);
+
+  std::cout << "Source vector" << std::endl;
+  std::cout << src << std::endl;
+
+  std::cout << "Result vector" << std::endl;
+  std::cout << result << std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Testing the laplacian smearing operator on a point source   "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  LatticeFermion smeared  (&Grid); smeared = src;
+  for (int smr = 0; smr < 10; ++smr)
+  {
+      LaplaceOperator.M(smeared, tmp);
+      smeared += 0.1*tmp;
+  }
+
+  std::cout << "Smeared vector" << std::endl;
+  std::cout << smeared << std::endl;
+
+  // Norm of vector
+  LatticeComplex smr_norm(&Grid);
+  smr_norm = localNorm2(smeared);
+  std::cout << "Smeared vector norm" << std::endl;
+  std::cout << smr_norm << std::endl;
+
+
+
+  Grid_finalize();
+}