diff --git a/examples/Example_Laplacian_inverter.cc b/examples/Example_Laplacian_inverter.cc deleted file mode 100644 index 1235d2b8..00000000 --- a/examples/Example_Laplacian_inverter.cc +++ /dev/null @@ -1,122 +0,0 @@ -#include -using namespace Grid; - -// Function used for Chebyshev smearing -// -Real MomentumSmearing(Real p2) -{ - return (1 - 4.0*p2) * exp(-p2/4); -} -Real DistillationSmearing(Real p2) -{ - if ( p2 > 0.5 ) return 0.0; - else return 1.0; -} - -// Flip sign to make prop to p^2, not -p^2 relative to last example -template class CovariantLaplacianCshift : public SparseMatrixBase -{ -public: - INHERIT_GIMPL_TYPES(Gimpl); - - GridBase *grid; - GaugeField U; - - CovariantLaplacianCshift(GaugeField &_U) : - grid(_U.Grid()), - U(_U) { }; - - virtual GridBase *Grid(void) { return grid; }; - - virtual void M (const Field &in, Field &out) - { - out=Zero(); - for(int mu=0;mu(U, mu); // NB: Inefficent - out = out - Gimpl::CovShiftForward(Umu,mu,in); - out = out - Gimpl::CovShiftBackward(Umu,mu,in); - out = out + 2.0*in; - } - }; - 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 &out) {assert(0);}; // Unimplemented need only for multigrid -}; - - - -int main(int argc, char ** argv) -{ - Grid_init(&argc, &argv); - - typedef LatticeColourVector Field; - - auto latt_size = GridDefaultLatt(); - auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); - auto mpi_layout = GridDefaultMpi(); - - GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector({45,12,81,9})); - - - LatticeGaugeField U(&Grid); - - SU::ColdConfiguration(RNG,U); - - typedef CovariantLaplacianCshift Laplacian_t; - Laplacian_t Laplacian(U); - - ColourVector ColourKronecker; - ColourKronecker = Zero(); - ColourKronecker()()(0) = 1.0; - - Coordinate site({latt_size[0]/2, - latt_size[1]/2, - latt_size[2]/2, - 0}); - - Field kronecker(&Grid); - kronecker = Zero(); - pokeSite(ColourKronecker,kronecker,site); - - - Field psi(&Grid), chi(&Grid); - - ////////////////////////////////////// - // Classic Wuppertal smearing - ////////////////////////////////////// - - Integer Iterations = 80; - Real width = 2.0; - Real coeff = (width*width) / Real(4*Iterations); - - chi=kronecker; - // chi = (1-p^2/2N)^N kronecker - for(int n = 0; n < Iterations; ++n) { - Laplacian.M(chi,psi); - chi = chi - coeff*psi; - } - - std::cout << " Wuppertal smeared operator is chi = \n" << chi < HermOp(Laplacian); - - Chebyshev ChebySmear(lo,hi,20,DistillationSmearing); - // Chebyshev ChebySmear(lo,hi,20,MomentumSmearing); - { - std::ofstream of("chebysmear"); - ChebySmear.csv(of); - } - - ChebySmear(HermOp,kronecker,chi); - - std::cout << " Chebyshev smeared operator is chi = \n" << chi <