portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Merge branch 'develop' of https://github.com/paboyle/Grid into develop

Browse Source
This commit is contained in:
fionnoh 2019-03-19 16:35:48 +00:00
commit 9790926cc5
14 changed files with 624 additions and 192 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Grid/Grid.h
View File

@ -42,6 +42,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#include <Grid/GridQCDcore.h>
#include <Grid/qcd/action/Action.h>
#include <Grid/qcd/utils/GaugeFix.h>
#include <Grid/qcd/utils/CovariantSmearing.h>
#include <Grid/qcd/smearing/Smearing.h>
#include <Grid/parallelIO/MetaData.h>
#include <Grid/qcd/hmc/HMC_aggregate.h>

87
Grid/qcd/utils/CovariantSmearing.h Normal file
View File

@ -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 Gimpl> class CovariantSmearing : public Gimpl
{
public:
INHERIT_GIMPL_TYPES(Gimpl);
typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz;
template<typename T>
static void GaussianSmear(const std::vector<LatticeColourMatrix>& 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<Nd;mu++) {
if ( mu != orthog ) {
psi = psi + Gimpl::CovShiftForward(U[mu],mu,chi);
psi = psi + Gimpl::CovShiftBackward(U[mu],mu,chi);
}
}
chi = chi + coeff*psi;
}
}
};
}}

23
Grid/simd/Grid_vector_types.h
View File

@ -10,6 +10,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Guido Cossu <cossu@iroiro-pc.kek.jp>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: neo <cossu@post.kek.jp>
Author: Michael Marshall <michael.marshall@ed.ac.au>
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 <typename Condition, typename ReturnType> using NotEnableIf = Invoke<st
////////////////////////////////////////////////////////
// Check for complexity with type traits
template <typename T> struct is_complex : public std::false_type {};
template <> struct is_complex<std::complex<double> > : public std::true_type {};
template <> struct is_complex<std::complex<float> > : public std::true_type {};
template <> struct is_complex<ComplexD> : public std::true_type {};
template <> struct is_complex<ComplexF> : public std::true_type {};
template <typename T> using IfReal = Invoke<std::enable_if<std::is_floating_point<T>::value, int> >;
template<typename T, typename V=void> struct is_real : public std::false_type {};
template<typename T> struct is_real<T, typename std::enable_if<std::is_floating_point<T>::value,
void>::type> : public std::true_type {};
template<typename T, typename V=void> struct is_integer : public std::false_type {};
template<typename T> struct is_integer<T, typename std::enable_if<std::is_integral<T>::value,
void>::type> : public std::true_type {};
template <typename T> using IfReal = Invoke<std::enable_if<is_real<T>::value, int> >;
template <typename T> using IfComplex = Invoke<std::enable_if<is_complex<T>::value, int> >;
template <typename T> using IfInteger = Invoke<std::enable_if<std::is_integral<T>::value, int> >;
template <typename T> using IfInteger = Invoke<std::enable_if<is_integer<T>::value, int> >;
template <typename T1,typename T2> using IfSame = Invoke<std::enable_if<std::is_same<T1,T2>::value, int> >;
template <typename T> using IfNotReal = Invoke<std::enable_if<!std::is_floating_point<T>::value, int> >;
template <typename T> using IfNotReal = Invoke<std::enable_if<!is_real<T>::value, int> >;
template <typename T> using IfNotComplex = Invoke<std::enable_if<!is_complex<T>::value, int> >;
template <typename T> using IfNotInteger = Invoke<std::enable_if<!std::is_integral<T>::value, int> >;
template <typename T> using IfNotInteger = Invoke<std::enable_if<!is_integer<T>::value, int> >;
template <typename T1,typename T2> using IfNotSame = Invoke<std::enable_if<!std::is_same<T1,T2>::value, int> >;
////////////////////////////////////////////////////////
@ -857,8 +866,10 @@ template <typename T>
struct is_simd : public std::false_type {};
template <> struct is_simd<vRealF> : public std::true_type {};
template <> struct is_simd<vRealD> : public std::true_type {};
template <> struct is_simd<vRealH> : public std::true_type {};
template <> struct is_simd<vComplexF> : public std::true_type {};
template <> struct is_simd<vComplexD> : public std::true_type {};
template <> struct is_simd<vComplexH> : public std::true_type {};
template <> struct is_simd<vInteger> : public std::true_type {};
template <typename T> using IfSimd = Invoke<std::enable_if<is_simd<T>::value, int> >;

89
Grid/tensors/Tensor_class.h
View File

@ -5,6 +5,7 @@ Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Michael Marshall <michael.marshall@ed.ac.au>
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 vtype>
class iScalar {
public:
vtype _internal;
typedef vtype element;
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
typedef iScalar<tensor_reduced_v> tensor_reduced;
typedef iScalar<recurse_scalar_object> scalar_object;
// substitutes a real or complex version with same tensor structure
typedef iScalar<typename GridTypeMapper<vtype>::Complexified> Complexified;
typedef iScalar<typename GridTypeMapper<vtype>::Realified> Realified;
// get double precision version
typedef iScalar<typename GridTypeMapper<vtype>::DoublePrecision> DoublePrecision;
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
using Traits = GridTypeMapper<iScalar<vtype> >;
GridVector_CopyTraits;
// Scalar no action
// template<int Level> using tensor_reduce_level = typename
@ -173,7 +173,10 @@ class iScalar {
return stream;
};
strong_inline const scalar_type * begin() const { return reinterpret_cast<const scalar_type *>(&_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<scalar<scalar<double>>>> back to double.
@ -194,21 +197,8 @@ class iVector {
public:
vtype _internal[N];
typedef vtype element;
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
typedef iScalar<tensor_reduced_v> tensor_reduced;
typedef iVector<recurse_scalar_object, N> scalar_object;
// substitutes a real or complex version with same tensor structure
typedef iVector<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
typedef iVector<typename GridTypeMapper<vtype>::Realified, N> Realified;
// get double precision version
typedef iVector<typename GridTypeMapper<vtype>::DoublePrecision, N> DoublePrecision;
using Traits = GridTypeMapper<iVector<vtype, N> >;
GridVector_CopyTraits;
template <class T, typename std::enable_if<!isGridTensor<T>::value, T>::type
* = nullptr>
@ -218,7 +208,6 @@ class iVector {
return *this;
}
enum { TensorLevel = GridTypeMapper<vtype>::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<const scalar_type *>(_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 <class vtype, int N>
@ -310,25 +304,8 @@ class iMatrix {
public:
vtype _internal[N][N];
typedef vtype element;
typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
typedef typename GridTypeMapper<vtype>::vector_type vector_type;
typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
// substitutes a real or complex version with same tensor structure
typedef iMatrix<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
typedef iMatrix<typename GridTypeMapper<vtype>::Realified, N> Realified;
// get double precision version
typedef iMatrix<typename GridTypeMapper<vtype>::DoublePrecision, N> DoublePrecision;
// Tensor removal
typedef iScalar<tensor_reduced_v> tensor_reduced;
typedef iMatrix<recurse_scalar_object, N> scalar_object;
enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
using Traits = GridTypeMapper<iMatrix<vtype, N> >;
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<const scalar_type *>(_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 <class v>
@ -480,6 +462,3 @@ void vprefetch(const iMatrix<v, N> &vv) {
}
}
#endif

197
Grid/tensors/Tensor_traits.h
View File

@ -5,6 +5,7 @@
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Christopher Kelly <ckelly@phys.columbia.edu>
Author: Michael Marshall <michael.marshall@ed.ac.au>
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 <ckelly@phys.columbia.edu>
namespace Grid {
// Forward declarations
template<class T> class iScalar;
template<class T, int N> class iVector;
template<class T, int N> class iMatrix;
// These are the Grid tensors
template<typename T> struct isGridTensor : public std::false_type { static constexpr bool notvalue = true; };
template<class T> struct isGridTensor<iScalar<T>> : public std::true_type { static constexpr bool notvalue = false; };
template<class T, int N> struct isGridTensor<iVector<T, N>> : public std::true_type { static constexpr bool notvalue = false; };
template<class T, int N> struct isGridTensor<iMatrix<T, N>> : public std::true_type { static constexpr bool notvalue = false; };
//////////////////////////////////////////////////////////////////////////////////
// Want to recurse: GridTypeMapper<Matrix<vComplexD> >::scalar_type == ComplexD.
// Use of a helper class like this allows us to template specialise and "dress"
@ -41,24 +53,25 @@ namespace Grid {
//
//////////////////////////////////////////////////////////////////////////////////
template <class T> 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<RealF> {
public:
template<typename T> struct GridTypeMapper {};
template<> struct GridTypeMapper<RealF> : 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<RealD> {
public:
template<> struct GridTypeMapper<RealD> : 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<ComplexF> {
public:
template<> struct GridTypeMapper<ComplexF> : 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<ComplexD> {
public:
template<> struct GridTypeMapper<ComplexD> : 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<Integer> {
public:
template<> struct GridTypeMapper<Integer> : 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<vRealF> {
public:
template<> struct GridTypeMapper<vRealF> : 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<vRealD> {
public:
template<> struct GridTypeMapper<vRealD> : 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<vComplexH> {
public:
template<> struct GridTypeMapper<vRealH> : 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<vComplexH> : 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<vComplexF> {
public:
template<> struct GridTypeMapper<vComplexF> : 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<vComplexD> {
public:
template<> struct GridTypeMapper<vComplexD> : 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<vInteger> {
public:
template<> struct GridTypeMapper<vInteger> : 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<typename T> struct isGridTensor {
static const bool value = true;
static const bool notvalue = false;
#define GridTypeMapper_RepeatedTypes \
using BaseTraits = GridTypeMapper<T>; \
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<typename T> struct GridTypeMapper<iScalar<T>> {
GridTypeMapper_RepeatedTypes;
using tensor_reduced = iScalar<typename BaseTraits::tensor_reduced>;
using scalar_object = iScalar<typename BaseTraits::scalar_object>;
using Complexified = iScalar<typename BaseTraits::Complexified>;
using Realified = iScalar<typename BaseTraits::Realified>;
using DoublePrecision = iScalar<typename BaseTraits::DoublePrecision>;
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<int > {
static const bool value = false;
static const bool notvalue = true;
template<typename T, int N> struct GridTypeMapper<iVector<T, N>> {
GridTypeMapper_RepeatedTypes;
using tensor_reduced = iScalar<typename BaseTraits::tensor_reduced>;
using scalar_object = iVector<typename BaseTraits::scalar_object, N>;
using Complexified = iVector<typename BaseTraits::Complexified, N>;
using Realified = iVector<typename BaseTraits::Realified, N>;
using DoublePrecision = iVector<typename BaseTraits::DoublePrecision, N>;
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<RealD > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<RealF > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<ComplexD > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<ComplexF > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<Integer > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<vRealD > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<vRealF > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<vComplexD > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<vComplexF > {
static const bool value = false;
static const bool notvalue = true;
};
template<> struct isGridTensor<vInteger > {
static const bool value = false;
static const bool notvalue = true;
template<typename T, int N> struct GridTypeMapper<iMatrix<T, N>> {
GridTypeMapper_RepeatedTypes;
using tensor_reduced = iScalar<typename BaseTraits::tensor_reduced>;
using scalar_object = iMatrix<typename BaseTraits::scalar_object, N>;
using Complexified = iMatrix<typename BaseTraits::Complexified, N>;
using Realified = iMatrix<typename BaseTraits::Realified, N>;
using DoublePrecision = iMatrix<typename BaseTraits::DoublePrecision, N>;
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,19 +261,12 @@ namespace Grid {
typedef T type;
};
//Query if a tensor or Lattice<Tensor> is SIMD vector or scalar
template<typename T>
class isSIMDvectorized{
template<typename U>
static typename std::enable_if< !std::is_same< typename GridTypeMapper<typename getVectorType<U>::type>::scalar_type,
typename GridTypeMapper<typename getVectorType<U>::type>::vector_type>::value, char>::type test(void *);
template<typename U>
static double test(...);
public:
enum {value = sizeof(test<T>(0)) == sizeof(char) };
};
//Query whether a tensor or Lattice<Tensor> is SIMD vector or scalar
template<typename T, typename V=void> struct isSIMDvectorized : public std::false_type {};
template<typename U> struct isSIMDvectorized<U, typename std::enable_if< !std::is_same<
typename GridTypeMapper<typename getVectorType<U>::type>::scalar_type,
typename GridTypeMapper<typename getVectorType<U>::type>::vector_type>::value, void>::type>
: public std::true_type {};
//Get the precision of a Lattice, tensor or scalar type in units of sizeof(float)
template<typename T>

2
Hadrons/Environment.hpp
View File

@ -182,7 +182,7 @@ private:
std::map<CoarseGridKey, GridPt> gridCoarse5d_;
unsigned int nd_;
// random number generator
RngPt rng4d_;
RngPt rng4d_{nullptr};
// object store
std::vector<ObjInfo> object_;
std::map<std::string, unsigned int> objectAddress_;

20
Hadrons/Modules/MAction/DWF.hpp
View File

@ -112,8 +112,6 @@ void TDWF<FImpl>::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<FImpl>::setup(void)
auto &g5 = *envGetGrid(FermionField, par().Ls);
auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
typename DomainWallFermion<FImpl>::ImplParams implParams;
if (!par().boundary.empty())
{
implParams.boundary_phases = strToVec<Complex>(par().boundary);
}
if (!par().twist.empty())
{
implParams.twist_n_2pi_L = strToVec<Real>(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<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, implParams);
}

20
Hadrons/Modules/MAction/MobiusDWF.hpp
View File

@ -112,8 +112,6 @@ void TMobiusDWF<FImpl>::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);
@ -121,8 +119,26 @@ void TMobiusDWF<FImpl>::setup(void)
auto &g5 = *envGetGrid(FermionField, par().Ls);
auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
typename MobiusFermion<FImpl>::ImplParams implParams;
if (!par().boundary.empty())
{
implParams.boundary_phases = strToVec<Complex>(par().boundary);
}
if (!par().twist.empty())
{
implParams.twist_n_2pi_L = strToVec<Real>(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<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
implParams);

20
Hadrons/Modules/MAction/ScaledDWF.hpp
View File

@ -111,8 +111,6 @@ void TScaledDWF<FImpl>::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<FImpl>::setup(void)
auto &g5 = *envGetGrid(FermionField, par().Ls);
auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
typename ScaledShamirFermion<FImpl>::ImplParams implParams;
if (!par().boundary.empty())
{
implParams.boundary_phases = strToVec<Complex>(par().boundary);
}
if (!par().twist.empty())
{
implParams.twist_n_2pi_L = strToVec<Real>(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<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, par().scale,
implParams);

20
Hadrons/Modules/MAction/Wilson.hpp
View File

@ -109,15 +109,31 @@ void TWilson<FImpl>::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<FImpl>::ImplParams implParams;
if (!par().boundary.empty())
{
implParams.boundary_phases = strToVec<Complex>(par().boundary);
}
if (!par().twist.empty())
{
implParams.twist_n_2pi_L = strToVec<Real>(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<FImpl>, getName(), 1, U, grid, gridRb,
par().mass, implParams);
}

21
Hadrons/Modules/MAction/WilsonClover.hpp
View File

@ -112,17 +112,34 @@ void TWilsonClover<FImpl>::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<FImpl>::ImplParams implParams;
if (!par().boundary.empty())
{
implParams.boundary_phases = strToVec<Complex>(par().boundary);
}
if (!par().twist.empty())
{
implParams.twist_n_2pi_L = strToVec<Real>(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<FImpl>, getName(), 1, U, grid,
gridRb, par().mass, par().csw_r, par().csw_t,
par().clover_anisotropy, implParams);

21
Hadrons/Modules/MAction/ZMobiusDWF.hpp
View File

@ -118,10 +118,7 @@ void TZMobiusDWF<FImpl>::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<FImpl>::setup(void)
auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
auto omega = par().omega;
typename ZMobiusFermion<FImpl>::ImplParams implParams;
if (!par().boundary.empty())
{
implParams.boundary_phases = strToVec<Complex>(par().boundary);
}
if (!par().twist.empty())
{
implParams.twist_n_2pi_L = strToVec<Real>(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<FImpl>, getName(), par().Ls, U, g5,
grb5, g4, grb4, par().mass, par().M5, omega,
par().b, par().c, implParams);

73
tests/smearing/Test_smearing.cc Normal file
View File

@ -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 <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);
std::vector<int> 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<LatticeColourMatrix> U(4,&Grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
std::vector<int> site({4,4,0,0});
src=zero;
SpinColourVector scv;
scv=zero;
scv()(0)(0) = 1.0;
pokeSite(scv,src,site);
CovariantSmearing<PeriodicGimplR>::GaussianSmear(U, src, 2.0, 50, Tdir);
std::cout << src <<std::endl;
}

194
tests/solver/Test_wilsonclover_mg_lime.cc Normal file
View File

@ -0,0 +1,194 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/solver/Test_wilsonclover_mg_mp.cc
Copyright (C) 2015-2018
Author: Daniel Richtmann <daniel.richtmann@ur.de>
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 <Grid/Grid.h>
#include <Test_multigrid_common.h>
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<int> 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 <<GridLogMessage<<"**************************************"<<std::endl;
std::cout <<GridLogMessage<<"** Reading back ILDG conf *********"<<std::endl;
std::cout <<GridLogMessage<<"**************************************"<<std::endl;
IldgReader _IldgReader;
_IldgReader.open(file);
_IldgReader.readConfiguration(Umu_d,header);
_IldgReader.close();
}
#else
{
FieldMetaData header;
std::string file("./ckpoint_lat.IEEE64BIG.1100");
NerscIO::readConfiguration(Umu_d,header,file);
}
#endif
// SU3::HotConfiguration(fPRNG, Umu_d);
LatticeGaugeFieldF Umu_f(FGrid_f); precisionChange(Umu_f, Umu_d);
// clang-format on
RealD mass = -0.25;
RealD csw_r = 1.0;
RealD csw_t = 1.0;
MultiGridParams mgParams;
std::string inputXml{"./mg_params.xml"};
if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {
inputXml = GridCmdOptionPayload(argv, argv + argc, "--inputxml");
assert(inputXml.length() != 0);
}
{
XmlWriter writer("mg_params_template.xml");
write(writer, "Params", mgParams);
std::cout << GridLogMessage << "Written mg_params_template.xml" << std::endl;
XmlReader reader(inputXml);
read(reader, "Params", mgParams);
std::cout << GridLogMessage << "Read in " << inputXml << std::endl;
}
checkParameterValidity(mgParams);
std::cout << mgParams << std::endl;
LevelInfo levelInfo_d(FGrid_d, mgParams);
LevelInfo levelInfo_f(FGrid_f, mgParams);
// Note: We do chiral doubling, so actually only nbasis/2 full basis vectors are used
const int nbasis = 40;
WilsonCloverFermionD Dwc_d(Umu_d, *FGrid_d, *FrbGrid_d, mass, csw_r, csw_t);
WilsonCloverFermionF Dwc_f(Umu_f, *FGrid_f, *FrbGrid_f, mass, csw_r, csw_t);
MdagMLinearOperator<WilsonCloverFermionD, LatticeFermionD> MdagMOpDwc_d(Dwc_d);
MdagMLinearOperator<WilsonCloverFermionF, LatticeFermionF> 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<vSpinColourVectorF, vTComplexF, nbasis, WilsonCloverFermionF>(mgParams, levelInfo_f, Dwc_f, Dwc_f);
MGPreconDwc_f->setup();
if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) {
MGPreconDwc_f->runChecks(1e-6);
}
MixedPrecisionFlexibleGeneralisedMinimalResidual<LatticeFermionD, LatticeFermionF> 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<vSpinColourVectorD, vTComplexD, nbasis, WilsonCloverFermionD>(mgParams, levelInfo_d, Dwc_d, Dwc_d);
MGPreconDwc_d->setup();
if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) {
MGPreconDwc_d->runChecks(1e-13);
}
FlexibleGeneralisedMinimalResidual<LatticeFermionD> 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();
}