diff --git a/Grid/Grid.h b/Grid/Grid.h index 9dcc207b..70859c81 100644 --- a/Grid/Grid.h +++ b/Grid/Grid.h @@ -42,6 +42,7 @@ Author: paboyle #include #include #include +#include #include #include #include diff --git a/Grid/qcd/utils/CovariantSmearing.h b/Grid/qcd/utils/CovariantSmearing.h new file mode 100644 index 00000000..eae5ed71 --- /dev/null +++ b/Grid/qcd/utils/CovariantSmearing.h @@ -0,0 +1,87 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h + +Copyright (C) 2016 + +Author: Azusa Yamaguchi + +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 +*************************************************************************************/ +#pragma once + +namespace Grid { +namespace QCD { + +template class CovariantSmearing : public Gimpl +{ +public: + INHERIT_GIMPL_TYPES(Gimpl); + + typedef typename Gimpl::GaugeLinkField GaugeMat; + typedef typename Gimpl::GaugeField GaugeLorentz; + + template + static void GaussianSmear(const std::vector& U, + T& chi, + const Real& width, int Iterations, int orthog) + { + GridBase *grid = chi._grid; + T psi(grid); + + //////////////////////////////////////////////////////////////////////////////////// + // Follow Chroma conventions for width to keep compatibility with previous data + // Free field iterates + // chi = (1 - w^2/4N p^2)^N chi + // + // ~ (e^(-w^2/4N p^2)^N chi + // ~ (e^(-w^2/4 p^2) chi + // ~ (e^(-w'^2/2 p^2) chi [ w' = w/sqrt(2) ] + // + // Which in coordinate space is proportional to + // + // e^(-x^2/w^2) = e^(-x^2/2w'^2) + // + // The 4 is a bit unconventional from Gaussian width perspective, but... it's Chroma convention. + // 2nd derivative approx d^2/dx^2 = x+mu + x-mu - 2x + // + // d^2/dx^2 = - p^2 + // + // chi = ( 1 + w^2/4N d^2/dx^2 )^N chi + // + //////////////////////////////////////////////////////////////////////////////////// + Real coeff = (width*width) / Real(4*Iterations); + + int dims = Nd; + if( orthog < Nd ) dims=Nd-1; + + for(int n = 0; n < Iterations; ++n) { + psi = (-2.0*dims)*chi; + for(int mu=0;mu Author: Guido Cossu Author: Peter Boyle Author: neo +Author: Michael Marshall 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 @@ -89,17 +90,25 @@ template using NotEnableIf = Invoke struct is_complex : public std::false_type {}; -template <> struct is_complex > : public std::true_type {}; -template <> struct is_complex > : public std::true_type {}; +template <> struct is_complex : public std::true_type {}; +template <> struct is_complex : public std::true_type {}; -template using IfReal = Invoke::value, int> >; +template struct is_real : public std::false_type {}; +template struct is_real::value, + void>::type> : public std::true_type {}; + +template struct is_integer : public std::false_type {}; +template struct is_integer::value, + void>::type> : public std::true_type {}; + +template using IfReal = Invoke::value, int> >; template using IfComplex = Invoke::value, int> >; -template using IfInteger = Invoke::value, int> >; +template using IfInteger = Invoke::value, int> >; template using IfSame = Invoke::value, int> >; -template using IfNotReal = Invoke::value, int> >; +template using IfNotReal = Invoke::value, int> >; template using IfNotComplex = Invoke::value, int> >; -template using IfNotInteger = Invoke::value, int> >; +template using IfNotInteger = Invoke::value, int> >; template using IfNotSame = Invoke::value, int> >; //////////////////////////////////////////////////////// @@ -857,8 +866,10 @@ template struct is_simd : public std::false_type {}; template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; +template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; +template <> struct is_simd : public std::true_type {}; template <> struct is_simd : public std::true_type {}; template using IfSimd = Invoke::value, int> >; diff --git a/Grid/tensors/Tensor_class.h b/Grid/tensors/Tensor_class.h index c7f868db..d59640df 100644 --- a/Grid/tensors/Tensor_class.h +++ b/Grid/tensors/Tensor_class.h @@ -5,6 +5,7 @@ Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle +Author: Michael Marshall 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 @@ -42,27 +43,26 @@ namespace Grid { // class GridTensorBase {}; +// Too late to remove these traits from Grid Tensors, so inherit from GridTypeMapper +#define GridVector_CopyTraits \ + using element = vtype; \ + using scalar_type = typename Traits::scalar_type; \ + using vector_type = typename Traits::vector_type; \ + using vector_typeD = typename Traits::vector_typeD; \ + using tensor_reduced = typename Traits::tensor_reduced; \ + using scalar_object = typename Traits::scalar_object; \ + using Complexified = typename Traits::Complexified; \ + using Realified = typename Traits::Realified; \ + using DoublePrecision = typename Traits::DoublePrecision; \ + static constexpr int TensorLevel = Traits::TensorLevel + template class iScalar { public: vtype _internal; - typedef vtype element; - typedef typename GridTypeMapper::scalar_type scalar_type; - typedef typename GridTypeMapper::vector_type vector_type; - typedef typename GridTypeMapper::vector_typeD vector_typeD; - typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; - typedef typename GridTypeMapper::scalar_object recurse_scalar_object; - typedef iScalar tensor_reduced; - typedef iScalar scalar_object; - // substitutes a real or complex version with same tensor structure - typedef iScalar::Complexified> Complexified; - typedef iScalar::Realified> Realified; - - // get double precision version - typedef iScalar::DoublePrecision> DoublePrecision; - - enum { TensorLevel = GridTypeMapper::TensorLevel + 1 }; + using Traits = GridTypeMapper >; + GridVector_CopyTraits; // Scalar no action // template using tensor_reduce_level = typename @@ -173,7 +173,10 @@ class iScalar { return stream; }; - + strong_inline const scalar_type * begin() const { return reinterpret_cast(&_internal); } + strong_inline scalar_type * begin() { return reinterpret_cast< scalar_type *>(&_internal); } + strong_inline const scalar_type * end() const { return begin() + Traits::count; } + strong_inline scalar_type * end() { return begin() + Traits::count; } }; /////////////////////////////////////////////////////////// // Allows to turn scalar>>> back to double. @@ -194,22 +197,9 @@ class iVector { public: vtype _internal[N]; - typedef vtype element; - typedef typename GridTypeMapper::scalar_type scalar_type; - typedef typename GridTypeMapper::vector_type vector_type; - typedef typename GridTypeMapper::vector_typeD vector_typeD; - typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; - typedef typename GridTypeMapper::scalar_object recurse_scalar_object; - typedef iScalar tensor_reduced; - typedef iVector scalar_object; + using Traits = GridTypeMapper >; + GridVector_CopyTraits; - // substitutes a real or complex version with same tensor structure - typedef iVector::Complexified, N> Complexified; - typedef iVector::Realified, N> Realified; - - // get double precision version - typedef iVector::DoublePrecision, N> DoublePrecision; - template ::value, T>::type * = nullptr> strong_inline auto operator=(T arg) -> iVector { @@ -218,7 +208,6 @@ class iVector { return *this; } - enum { TensorLevel = GridTypeMapper::TensorLevel + 1 }; iVector(const Zero &z) { *this = zero; }; iVector() = default; /* @@ -303,6 +292,11 @@ class iVector { // strong_inline vtype && operator ()(int i) { // return _internal[i]; // } + + strong_inline const scalar_type * begin() const { return reinterpret_cast(_internal); } + strong_inline scalar_type * begin() { return reinterpret_cast< scalar_type *>(_internal); } + strong_inline const scalar_type * end() const { return begin() + Traits::count; } + strong_inline scalar_type * end() { return begin() + Traits::count; } }; template @@ -310,25 +304,8 @@ class iMatrix { public: vtype _internal[N][N]; - typedef vtype element; - typedef typename GridTypeMapper::scalar_type scalar_type; - typedef typename GridTypeMapper::vector_type vector_type; - typedef typename GridTypeMapper::vector_typeD vector_typeD; - typedef typename GridTypeMapper::tensor_reduced tensor_reduced_v; - typedef typename GridTypeMapper::scalar_object recurse_scalar_object; - - // substitutes a real or complex version with same tensor structure - typedef iMatrix::Complexified, N> Complexified; - typedef iMatrix::Realified, N> Realified; - - // get double precision version - typedef iMatrix::DoublePrecision, N> DoublePrecision; - - // Tensor removal - typedef iScalar tensor_reduced; - typedef iMatrix scalar_object; - - enum { TensorLevel = GridTypeMapper::TensorLevel + 1 }; + using Traits = GridTypeMapper >; + GridVector_CopyTraits; iMatrix(const Zero &z) { *this = zero; }; iMatrix() = default; @@ -458,6 +435,11 @@ class iMatrix { // strong_inline vtype && operator ()(int i,int j) { // return _internal[i][j]; // } + + strong_inline const scalar_type * begin() const { return reinterpret_cast(_internal[0]); } + strong_inline scalar_type * begin() { return reinterpret_cast< scalar_type *>(_internal[0]); } + strong_inline const scalar_type * end() const { return begin() + Traits::count; } + strong_inline scalar_type * end() { return begin() + Traits::count; } }; template @@ -480,6 +462,3 @@ void vprefetch(const iMatrix &vv) { } } #endif - - - diff --git a/Grid/tensors/Tensor_traits.h b/Grid/tensors/Tensor_traits.h index c1ef397a..9cb93e17 100644 --- a/Grid/tensors/Tensor_traits.h +++ b/Grid/tensors/Tensor_traits.h @@ -5,6 +5,7 @@ Author: Azusa Yamaguchi Author: Peter Boyle Author: Christopher Kelly +Author: Michael Marshall 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 @@ -26,6 +27,17 @@ Author: Christopher Kelly namespace Grid { + // Forward declarations + template class iScalar; + template class iVector; + template class iMatrix; + + // These are the Grid tensors + template struct isGridTensor : public std::false_type { static constexpr bool notvalue = true; }; + template struct isGridTensor> : public std::true_type { static constexpr bool notvalue = false; }; + template struct isGridTensor> : public std::true_type { static constexpr bool notvalue = false; }; + template struct isGridTensor> : public std::true_type { static constexpr bool notvalue = false; }; + ////////////////////////////////////////////////////////////////////////////////// // Want to recurse: GridTypeMapper >::scalar_type == ComplexD. // Use of a helper class like this allows us to template specialise and "dress" @@ -40,25 +52,26 @@ namespace Grid { // to study C++11's type_traits.h file. (std::enable_if >) // ////////////////////////////////////////////////////////////////////////////////// - - template class GridTypeMapper { - public: - typedef typename T::scalar_type scalar_type; - typedef typename T::vector_type vector_type; - typedef typename T::vector_typeD vector_typeD; - typedef typename T::tensor_reduced tensor_reduced; - typedef typename T::scalar_object scalar_object; - typedef typename T::Complexified Complexified; - typedef typename T::Realified Realified; - typedef typename T::DoublePrecision DoublePrecision; - enum { TensorLevel = T::TensorLevel }; + + // This saves repeating common properties for supported Grid Scalar types + // TensorLevel How many nested grid tensors + // Rank Rank of the grid tensor + // count Total number of elements, i.e. product of dimensions + // Dimension(dim) Size of dimension dim + struct GridTypeMapper_Base { + static constexpr int TensorLevel = 0; + static constexpr int Rank = 0; + static constexpr std::size_t count = 1; + static constexpr int Dimension(int dim) { return 0; } }; ////////////////////////////////////////////////////////////////////////////////// // Recursion stops with these template specialisations ////////////////////////////////////////////////////////////////////////////////// - template<> class GridTypeMapper { - public: + + template struct GridTypeMapper {}; + + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealF scalar_type; typedef RealF vector_type; typedef RealD vector_typeD; @@ -67,10 +80,8 @@ namespace Grid { typedef ComplexF Complexified; typedef RealF Realified; typedef RealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; typedef RealD vector_type; typedef RealD vector_typeD; @@ -79,10 +90,8 @@ namespace Grid { typedef ComplexD Complexified; typedef RealD Realified; typedef RealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; typedef ComplexF vector_type; typedef ComplexD vector_typeD; @@ -91,10 +100,8 @@ namespace Grid { typedef ComplexF Complexified; typedef RealF Realified; typedef ComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; typedef ComplexD vector_type; typedef ComplexD vector_typeD; @@ -103,10 +110,8 @@ namespace Grid { typedef ComplexD Complexified; typedef RealD Realified; typedef ComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef Integer scalar_type; typedef Integer vector_type; typedef Integer vector_typeD; @@ -115,11 +120,9 @@ namespace Grid { typedef void Complexified; typedef void Realified; typedef void DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealF scalar_type; typedef vRealF vector_type; typedef vRealD vector_typeD; @@ -128,10 +131,8 @@ namespace Grid { typedef vComplexF Complexified; typedef vRealF Realified; typedef vRealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef RealD scalar_type; typedef vRealD vector_type; typedef vRealD vector_typeD; @@ -140,10 +141,18 @@ namespace Grid { typedef vComplexD Complexified; typedef vRealD Realified; typedef vRealD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { + typedef RealF scalar_type; + typedef vRealH vector_type; + typedef vRealD vector_typeD; + typedef vRealH tensor_reduced; + typedef RealF scalar_object; + typedef vComplexH Complexified; + typedef vRealH Realified; + typedef vRealD DoublePrecision; + }; + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; typedef vComplexH vector_type; typedef vComplexD vector_typeD; @@ -152,10 +161,8 @@ namespace Grid { typedef vComplexH Complexified; typedef vRealH Realified; typedef vComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexF scalar_type; typedef vComplexF vector_type; typedef vComplexD vector_typeD; @@ -164,10 +171,8 @@ namespace Grid { typedef vComplexF Complexified; typedef vRealF Realified; typedef vComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef ComplexD scalar_type; typedef vComplexD vector_type; typedef vComplexD vector_typeD; @@ -176,10 +181,8 @@ namespace Grid { typedef vComplexD Complexified; typedef vRealD Realified; typedef vComplexD DoublePrecision; - enum { TensorLevel = 0 }; }; - template<> class GridTypeMapper { - public: + template<> struct GridTypeMapper : public GridTypeMapper_Base { typedef Integer scalar_type; typedef vInteger vector_type; typedef vInteger vector_typeD; @@ -188,57 +191,52 @@ namespace Grid { typedef void Complexified; typedef void Realified; typedef void DoublePrecision; - enum { TensorLevel = 0 }; }; - // First some of my own traits - template struct isGridTensor { - static const bool value = true; - static const bool notvalue = false; +#define GridTypeMapper_RepeatedTypes \ + using BaseTraits = GridTypeMapper; \ + using scalar_type = typename BaseTraits::scalar_type; \ + using vector_type = typename BaseTraits::vector_type; \ + using vector_typeD = typename BaseTraits::vector_typeD; \ + static constexpr int TensorLevel = BaseTraits::TensorLevel + 1 + + template struct GridTypeMapper> { + GridTypeMapper_RepeatedTypes; + using tensor_reduced = iScalar; + using scalar_object = iScalar; + using Complexified = iScalar; + using Realified = iScalar; + using DoublePrecision = iScalar; + static constexpr int Rank = BaseTraits::Rank + 1; + static constexpr std::size_t count = BaseTraits::count; + static constexpr int Dimension(int dim) { + return ( dim == 0 ) ? 1 : BaseTraits::Dimension(dim - 1); } }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; + + template struct GridTypeMapper> { + GridTypeMapper_RepeatedTypes; + using tensor_reduced = iScalar; + using scalar_object = iVector; + using Complexified = iVector; + using Realified = iVector; + using DoublePrecision = iVector; + static constexpr int Rank = BaseTraits::Rank + 1; + static constexpr std::size_t count = BaseTraits::count * N; + static constexpr int Dimension(int dim) { + return ( dim == 0 ) ? N : BaseTraits::Dimension(dim - 1); } }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; - }; - template<> struct isGridTensor { - static const bool value = false; - static const bool notvalue = true; + + template struct GridTypeMapper> { + GridTypeMapper_RepeatedTypes; + using tensor_reduced = iScalar; + using scalar_object = iMatrix; + using Complexified = iMatrix; + using Realified = iMatrix; + using DoublePrecision = iMatrix; + static constexpr int Rank = BaseTraits::Rank + 2; + static constexpr std::size_t count = BaseTraits::count * N * N; + static constexpr int Dimension(int dim) { + return ( dim == 0 || dim == 1 ) ? N : BaseTraits::Dimension(dim - 2); } }; // Match the index @@ -263,20 +261,13 @@ namespace Grid { typedef T type; }; - //Query if a tensor or Lattice is SIMD vector or scalar - template - class isSIMDvectorized{ - template - static typename std::enable_if< !std::is_same< typename GridTypeMapper::type>::scalar_type, - typename GridTypeMapper::type>::vector_type>::value, char>::type test(void *); + //Query whether a tensor or Lattice is SIMD vector or scalar + template struct isSIMDvectorized : public std::false_type {}; + template struct isSIMDvectorized::type>::scalar_type, + typename GridTypeMapper::type>::vector_type>::value, void>::type> + : public std::true_type {}; - template - static double test(...); - - public: - enum {value = sizeof(test(0)) == sizeof(char) }; - }; - //Get the precision of a Lattice, tensor or scalar type in units of sizeof(float) template class getPrecision{ diff --git a/Hadrons/Environment.hpp b/Hadrons/Environment.hpp index 05d3d8dd..24f6b31e 100644 --- a/Hadrons/Environment.hpp +++ b/Hadrons/Environment.hpp @@ -182,7 +182,7 @@ private: std::map gridCoarse5d_; unsigned int nd_; // random number generator - RngPt rng4d_; + RngPt rng4d_{nullptr}; // object store std::vector object_; std::map objectAddress_; diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp index bac70ae8..21024a4f 100644 --- a/Hadrons/Modules/MAction/DWF.hpp +++ b/Hadrons/Modules/MAction/DWF.hpp @@ -112,8 +112,6 @@ void TDWF::setup(void) << par().mass << ", M5= " << par().M5 << " and Ls= " << par().Ls << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); @@ -121,8 +119,26 @@ void TDWF::setup(void) auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); typename DomainWallFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, DomainWallFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, implParams); } diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp index 6b960e2e..3ce02b91 100644 --- a/Hadrons/Modules/MAction/MobiusDWF.hpp +++ b/Hadrons/Modules/MAction/MobiusDWF.hpp @@ -112,17 +112,33 @@ void TMobiusDWF::setup(void) << ", b= " << par().b << ", c= " << par().c << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; - + auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); auto &grb4 = *envGetRbGrid(FermionField); auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); typename MobiusFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, MobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().b, par().c, implParams); diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp index 8742a820..3e052225 100644 --- a/Hadrons/Modules/MAction/ScaledDWF.hpp +++ b/Hadrons/Modules/MAction/ScaledDWF.hpp @@ -111,8 +111,6 @@ void TScaledDWF::setup(void) << ", scale= " << par().scale << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); @@ -120,8 +118,26 @@ void TScaledDWF::setup(void) auto &g5 = *envGetGrid(FermionField, par().Ls); auto &grb5 = *envGetRbGrid(FermionField, par().Ls); typename ScaledShamirFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, ScaledShamirFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, par().scale, implParams); diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp index 3ce5569e..11df89c3 100644 --- a/Hadrons/Modules/MAction/Wilson.hpp +++ b/Hadrons/Modules/MAction/Wilson.hpp @@ -109,15 +109,31 @@ void TWilson::setup(void) { LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); typename WilsonFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, WilsonFermion, getName(), 1, U, grid, gridRb, par().mass, implParams); } diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp index a9ea8de5..d80bccdf 100644 --- a/Hadrons/Modules/MAction/WilsonClover.hpp +++ b/Hadrons/Modules/MAction/WilsonClover.hpp @@ -112,17 +112,34 @@ void TWilsonClover::setup(void) { LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass << " using gauge field '" << par().gauge << "'" << std::endl; - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; LOG(Message) << "Clover term csw_r: " << par().csw_r << " csw_t: " << par().csw_t << std::endl; + auto &U = envGet(GaugeField, par().gauge); auto &grid = *envGetGrid(FermionField); auto &gridRb = *envGetRbGrid(FermionField); typename WilsonCloverFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, WilsonCloverFermion, getName(), 1, U, grid, gridRb, par().mass, par().csw_r, par().csw_t, par().clover_anisotropy, implParams); diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp index e780bd73..56efe1a2 100644 --- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp +++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp @@ -118,10 +118,7 @@ void TZMobiusDWF::setup(void) { LOG(Message) << " omega[" << i << "]= " << par().omega[i] << std::endl; } - LOG(Message) << "Fermion boundary conditions: " << par().boundary - << std::endl; - env().createGrid(par().Ls); auto &U = envGet(GaugeField, par().gauge); auto &g4 = *envGetGrid(FermionField); auto &grb4 = *envGetRbGrid(FermionField); @@ -129,8 +126,26 @@ void TZMobiusDWF::setup(void) auto &grb5 = *envGetRbGrid(FermionField, par().Ls); auto omega = par().omega; typename ZMobiusFermion::ImplParams implParams; - implParams.boundary_phases = strToVec(par().boundary); - implParams.twist_n_2pi_L = strToVec(par().twist); + if (!par().boundary.empty()) + { + implParams.boundary_phases = strToVec(par().boundary); + } + if (!par().twist.empty()) + { + implParams.twist_n_2pi_L = strToVec(par().twist); + } + LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases + << std::endl; + LOG(Message) << "Twists: " << implParams.twist_n_2pi_L + << std::endl; + if (implParams.boundary_phases.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of boundary phase"); + } + if (implParams.twist_n_2pi_L.size() != env().getNd()) + { + HADRONS_ERROR(Size, "Wrong number of twist"); + } envCreateDerived(FMat, ZMobiusFermion, getName(), par().Ls, U, g5, grb5, g4, grb4, par().mass, par().M5, omega, par().b, par().c, implParams); diff --git a/tests/smearing/Test_smearing.cc b/tests/smearing/Test_smearing.cc new file mode 100644 index 00000000..996671f0 --- /dev/null +++ b/tests/smearing/Test_smearing.cc @@ -0,0 +1,73 @@ +/* + * + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_prec.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=zero; + LatticeGaugeField Umu(&Grid); + // SU3::HotConfiguration(pRNG,Umu); + SU3::ColdConfiguration(Umu); + std::vector U(4,&Grid); + + for(int mu=0;mu(Umu,mu); + } + + std::vector site({4,4,0,0}); + src=zero; + SpinColourVector scv; + scv=zero; + scv()(0)(0) = 1.0; + pokeSite(scv,src,site); + + CovariantSmearing::GaussianSmear(U, src, 2.0, 50, Tdir); + + std::cout << src < + + 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 + *************************************************************************************/ +/* */ + + +#include +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main(int argc, char **argv) { + + + Grid_init(&argc, &argv); + + // clang-format off + GridCartesian *FGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridCartesian *FGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_d); + GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_f); + // clang-format on + + std::vector fSeeds({1, 2, 3, 4}); + GridParallelRNG fPRNG(FGrid_d); + fPRNG.SeedFixedIntegers(fSeeds); + + // clang-format off + LatticeFermionD src_d(FGrid_d); gaussian(fPRNG, src_d); + LatticeFermionD resultMGD_d(FGrid_d); resultMGD_d = zero; + LatticeFermionD resultMGF_d(FGrid_d); resultMGF_d = zero; + LatticeGaugeFieldD Umu_d(FGrid_d); + +#if 0 + { + FieldMetaData header; + std::string file("./qcdsf.769.00399.lime"); + std::cout < MdagMOpDwc_d(Dwc_d); + MdagMLinearOperator MdagMOpDwc_f(Dwc_f); + + std::cout << GridLogMessage << "**************************************************" << std::endl; + std::cout << GridLogMessage << "Testing single-precision Multigrid for Wilson Clover" << std::endl; + std::cout << GridLogMessage << "**************************************************" << std::endl; + + auto MGPreconDwc_f = createMGInstance(mgParams, levelInfo_f, Dwc_f, Dwc_f); + + MGPreconDwc_f->setup(); + + if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) { + MGPreconDwc_f->runChecks(1e-6); + } + + MixedPrecisionFlexibleGeneralisedMinimalResidual MPFGMRESPREC( + 1.0e-12, 50000, FGrid_f, *MGPreconDwc_f, 100, false); + + std::cout << std::endl << "Starting with a new solver" << std::endl; + MPFGMRESPREC(MdagMOpDwc_d, src_d, resultMGF_d); + + MGPreconDwc_f->reportTimings(); + + if(GridCmdOptionExists(argv, argv + argc, "--docomparison")) { + + std::cout << GridLogMessage << "**************************************************" << std::endl; + std::cout << GridLogMessage << "Testing double-precision Multigrid for Wilson Clover" << std::endl; + std::cout << GridLogMessage << "**************************************************" << std::endl; + + auto MGPreconDwc_d = createMGInstance(mgParams, levelInfo_d, Dwc_d, Dwc_d); + + MGPreconDwc_d->setup(); + + if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) { + MGPreconDwc_d->runChecks(1e-13); + } + + FlexibleGeneralisedMinimalResidual FGMRESPREC(1.0e-12, 50000, *MGPreconDwc_d, 100, false); + + std::cout << std::endl << "Starting with a new solver" << std::endl; + FGMRESPREC(MdagMOpDwc_d, src_d, resultMGD_d); + + MGPreconDwc_d->reportTimings(); + + std::cout << GridLogMessage << "**************************************************" << std::endl; + std::cout << GridLogMessage << "Comparing single-precision Multigrid with double-precision one for Wilson Clover" << std::endl; + std::cout << GridLogMessage << "**************************************************" << std::endl; + + LatticeFermionD diffFullSolver(FGrid_d); + + RealD deviationFullSolver = axpy_norm(diffFullSolver, -1.0, resultMGF_d, resultMGD_d); + + // clang-format off + LatticeFermionF src_f(FGrid_f); precisionChange(src_f, src_d); + LatticeFermionF resMGF_f(FGrid_f); resMGF_f = zero; + LatticeFermionD resMGD_d(FGrid_d); resMGD_d = zero; + // clang-format on + + (*MGPreconDwc_f)(src_f, resMGF_f); + (*MGPreconDwc_d)(src_d, resMGD_d); + + LatticeFermionD diffOnlyMG(FGrid_d); + LatticeFermionD resMGF_d(FGrid_d); + precisionChange(resMGF_d, resMGF_f); + + RealD deviationOnlyPrec = axpy_norm(diffOnlyMG, -1.0, resMGF_d, resMGD_d); + + // clang-format off + std::cout << GridLogMessage << "Absolute difference between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver << std::endl; + std::cout << GridLogMessage << "Relative deviation between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver / norm2(resultMGD_d) << std::endl; + std::cout << GridLogMessage << "Absolute difference between one iteration of MG Prec in double and single precision: " << deviationOnlyPrec << std::endl; + std::cout << GridLogMessage << "Relative deviation between one iteration of MG Prec in double and single precision: " << deviationOnlyPrec / norm2(resMGD_d) << std::endl; + // clang-format on + } + + Grid_finalize(); +}