mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-12 20:27:06 +01:00
Merge branch 'develop' into feature/hadrons
This commit is contained in:
@ -425,7 +425,7 @@ namespace Grid {
|
||||
A[p]=zero;
|
||||
}
|
||||
|
||||
GridParallelRNG RNG(Grid()); RNG.SeedRandomDevice();
|
||||
GridParallelRNG RNG(Grid()); RNG.SeedFixedIntegers(std::vector<int>({55,72,19,17,34}));
|
||||
Lattice<iScalar<CComplex> > val(Grid()); random(RNG,val);
|
||||
|
||||
Complex one(1.0);
|
||||
|
@ -177,9 +177,11 @@ public:
|
||||
// Global addressing
|
||||
////////////////////////////////////////////////////////////////
|
||||
void GlobalIndexToGlobalCoor(int gidx,std::vector<int> &gcoor){
|
||||
assert(gidx< gSites());
|
||||
Lexicographic::CoorFromIndex(gcoor,gidx,_gdimensions);
|
||||
}
|
||||
void LocalIndexToLocalCoor(int lidx,std::vector<int> &lcoor){
|
||||
assert(lidx<lSites());
|
||||
Lexicographic::CoorFromIndex(lcoor,lidx,_ldimensions);
|
||||
}
|
||||
void GlobalCoorToGlobalIndex(const std::vector<int> & gcoor,int & gidx){
|
||||
|
@ -206,7 +206,7 @@ void CartesianCommunicator::Init(int *argc, char ***argv) {
|
||||
sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",GroupRank,r);
|
||||
|
||||
shm_unlink(shm_name);
|
||||
int fd=shm_open(shm_name,O_RDWR|O_CREAT,0660);
|
||||
int fd=shm_open(shm_name,O_RDWR|O_CREAT,0666);
|
||||
if ( fd < 0 ) { perror("failed shm_open"); assert(0); }
|
||||
ftruncate(fd, size);
|
||||
|
||||
@ -226,7 +226,7 @@ void CartesianCommunicator::Init(int *argc, char ***argv) {
|
||||
|
||||
sprintf(shm_name,"/Grid_mpi3_shm_%d_%d",GroupRank,r);
|
||||
|
||||
int fd=shm_open(shm_name,O_RDWR,0660);
|
||||
int fd=shm_open(shm_name,O_RDWR,0666);
|
||||
if ( fd<0 ) { perror("failed shm_open"); assert(0); }
|
||||
|
||||
void * ptr = mmap(NULL,size, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0);
|
||||
|
@ -30,12 +30,19 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
#define GRID_LATTICE_RNG_H
|
||||
|
||||
#include <random>
|
||||
|
||||
#ifdef RNG_SITMO
|
||||
#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
|
||||
#endif
|
||||
|
||||
#if defined(RNG_SITMO)
|
||||
#define RNG_FAST_DISCARD
|
||||
#else
|
||||
#undef RNG_FAST_DISCARD
|
||||
#endif
|
||||
|
||||
namespace Grid {
|
||||
|
||||
//http://nvlpubs.nist.gov/nistpubs/SpecialPublications/NIST.SP.800-90Ar1.pdf ?
|
||||
|
||||
//////////////////////////////////////////////////////////////
|
||||
// Allow the RNG state to be less dense than the fine grid
|
||||
//////////////////////////////////////////////////////////////
|
||||
@ -65,120 +72,139 @@ namespace Grid {
|
||||
|
||||
multiplicity = multiplicity *fine->_rdimensions[fd] / coarse->_rdimensions[d];
|
||||
}
|
||||
|
||||
return multiplicity;
|
||||
}
|
||||
|
||||
// Wrap seed_seq to give common interface with random_device
|
||||
// Should rather wrap random_device and have a generate
|
||||
class fixedSeed {
|
||||
public:
|
||||
|
||||
typedef std::seed_seq::result_type result_type;
|
||||
|
||||
std::seed_seq src;
|
||||
|
||||
template<class int_type> fixedSeed(const std::vector<int_type> &seeds) : src(seeds.begin(),seeds.end()) {};
|
||||
|
||||
template< class RandomIt > void generate( RandomIt begin, RandomIt end ) {
|
||||
src.generate(begin,end);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
class deviceSeed {
|
||||
public:
|
||||
|
||||
std::random_device rd;
|
||||
|
||||
typedef std::random_device::result_type result_type;
|
||||
|
||||
deviceSeed(void) : rd(){};
|
||||
|
||||
template< class RandomIt > void generate( RandomIt begin, RandomIt end ) {
|
||||
for(RandomIt it=begin; it!=end;it++){
|
||||
*it = rd();
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
// real scalars are one component
|
||||
template<class scalar,class distribution,class generator> void fillScalar(scalar &s,distribution &dist,generator & gen)
|
||||
template<class scalar,class distribution,class generator>
|
||||
void fillScalar(scalar &s,distribution &dist,generator & gen)
|
||||
{
|
||||
s=dist(gen);
|
||||
}
|
||||
template<class distribution,class generator> void fillScalar(ComplexF &s,distribution &dist, generator &gen)
|
||||
template<class distribution,class generator>
|
||||
void fillScalar(ComplexF &s,distribution &dist, generator &gen)
|
||||
{
|
||||
s=ComplexF(dist(gen),dist(gen));
|
||||
}
|
||||
template<class distribution,class generator> void fillScalar(ComplexD &s,distribution &dist,generator &gen)
|
||||
template<class distribution,class generator>
|
||||
void fillScalar(ComplexD &s,distribution &dist,generator &gen)
|
||||
{
|
||||
s=ComplexD(dist(gen),dist(gen));
|
||||
}
|
||||
|
||||
class GridRNGbase {
|
||||
|
||||
public:
|
||||
|
||||
int _seeded;
|
||||
// One generator per site.
|
||||
// Uniform and Gaussian distributions from these generators.
|
||||
#ifdef RNG_RANLUX
|
||||
typedef uint64_t RngStateType;
|
||||
typedef std::ranlux48 RngEngine;
|
||||
typedef uint64_t RngStateType;
|
||||
static const int RngStateCount = 15;
|
||||
#elif RNG_MT19937
|
||||
#endif
|
||||
#ifdef RNG_MT19937
|
||||
typedef std::mt19937 RngEngine;
|
||||
typedef uint32_t RngStateType;
|
||||
static const int RngStateCount = std::mt19937::state_size;
|
||||
#elif RNG_SITMO
|
||||
#endif
|
||||
#ifdef RNG_SITMO
|
||||
typedef sitmo::prng_engine RngEngine;
|
||||
typedef uint64_t RngStateType;
|
||||
static const int RngStateCount = 4;
|
||||
#endif
|
||||
std::vector<RngEngine> _generators;
|
||||
std::vector<std::uniform_real_distribution<RealD>> _uniform;
|
||||
std::vector<std::normal_distribution<RealD>> _gaussian;
|
||||
std::vector<std::discrete_distribution<int32_t>> _bernoulli;
|
||||
|
||||
void GetState(std::vector<RngStateType> & saved,int gen) {
|
||||
std::vector<RngEngine> _generators;
|
||||
std::vector<std::uniform_real_distribution<RealD> > _uniform;
|
||||
std::vector<std::normal_distribution<RealD> > _gaussian;
|
||||
std::vector<std::discrete_distribution<int32_t> > _bernoulli;
|
||||
std::vector<std::uniform_int_distribution<uint32_t> > _uid;
|
||||
|
||||
///////////////////////
|
||||
// support for parallel init
|
||||
///////////////////////
|
||||
#ifdef RNG_FAST_DISCARD
|
||||
static void Skip(RngEngine &eng)
|
||||
{
|
||||
/////////////////////////////////////////////////////////////////////////////////////
|
||||
// Skip by 2^40 elements between successive lattice sites
|
||||
// This goes by 10^12.
|
||||
// Consider quenched updating; likely never exceeding rate of 1000 sweeps
|
||||
// per second on any machine. This gives us of order 10^9 seconds, or 100 years
|
||||
// skip ahead.
|
||||
// For HMC unlikely to go at faster than a solve per second, and
|
||||
// tens of seconds per trajectory so this is clean in all reasonable cases,
|
||||
// and margin of safety is orders of magnitude.
|
||||
// We could hack Sitmo to skip in the higher order words of state if necessary
|
||||
/////////////////////////////////////////////////////////////////////////////////////
|
||||
uint64_t skip = 0x1; skip = skip<<40;
|
||||
eng.discard(skip);
|
||||
}
|
||||
#endif
|
||||
static RngEngine Reseed(RngEngine &eng)
|
||||
{
|
||||
std::vector<uint32_t> newseed;
|
||||
std::uniform_int_distribution<uint32_t> uid;
|
||||
return Reseed(eng,newseed,uid);
|
||||
}
|
||||
static RngEngine Reseed(RngEngine &eng,std::vector<uint32_t> & newseed,
|
||||
std::uniform_int_distribution<uint32_t> &uid)
|
||||
{
|
||||
const int reseeds=4;
|
||||
|
||||
newseed.resize(reseeds);
|
||||
for(int i=0;i<reseeds;i++){
|
||||
newseed[i] = uid(eng);
|
||||
}
|
||||
std::seed_seq sseq(newseed.begin(),newseed.end());
|
||||
return RngEngine(sseq);
|
||||
}
|
||||
|
||||
void GetState(std::vector<RngStateType> & saved,RngEngine &eng) {
|
||||
saved.resize(RngStateCount);
|
||||
std::stringstream ss;
|
||||
ss<<_generators[gen];
|
||||
ss<<eng;
|
||||
ss.seekg(0,ss.beg);
|
||||
for(int i=0;i<RngStateCount;i++){
|
||||
ss>>saved[i];
|
||||
}
|
||||
}
|
||||
void SetState(std::vector<RngStateType> & saved,int gen){
|
||||
void GetState(std::vector<RngStateType> & saved,int gen) {
|
||||
GetState(saved,_generators[gen]);
|
||||
}
|
||||
void SetState(std::vector<RngStateType> & saved,RngEngine &eng){
|
||||
assert(saved.size()==RngStateCount);
|
||||
std::stringstream ss;
|
||||
for(int i=0;i<RngStateCount;i++){
|
||||
ss<< saved[i]<<" ";
|
||||
}
|
||||
ss.seekg(0,ss.beg);
|
||||
ss>>_generators[gen];
|
||||
ss>>eng;
|
||||
}
|
||||
void SetState(std::vector<RngStateType> & saved,int gen){
|
||||
SetState(saved,_generators[gen]);
|
||||
}
|
||||
void SetEngine(RngEngine &Eng, int gen){
|
||||
_generators[gen]=Eng;
|
||||
}
|
||||
void GetEngine(RngEngine &Eng, int gen){
|
||||
Eng=_generators[gen];
|
||||
}
|
||||
template<class source> void Seed(source &src, int gen)
|
||||
{
|
||||
_generators[gen] = RngEngine(src);
|
||||
}
|
||||
};
|
||||
|
||||
class GridSerialRNG : public GridRNGbase {
|
||||
public:
|
||||
|
||||
// FIXME ... do we require lockstep draws of randoms
|
||||
// from all nodes keeping seeds consistent.
|
||||
// place a barrier/broadcast in the fill routine
|
||||
|
||||
GridSerialRNG() : GridRNGbase() {
|
||||
_generators.resize(1);
|
||||
_uniform.resize(1,std::uniform_real_distribution<RealD>{0,1});
|
||||
_gaussian.resize(1,std::normal_distribution<RealD>(0.0,1.0) );
|
||||
_bernoulli.resize(1,std::discrete_distribution<int32_t>{1,1});
|
||||
_seeded=0;
|
||||
_uid.resize(1,std::uniform_int_distribution<uint32_t>() );
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class sobj,class distribution> inline void fill(sobj &l,std::vector<distribution> &dist){
|
||||
|
||||
typedef typename sobj::scalar_type scalar_type;
|
||||
@ -191,7 +217,7 @@ namespace Grid {
|
||||
for(int idx=0;idx<words;idx++){
|
||||
fillScalar(buf[idx],dist[0],_generators[0]);
|
||||
}
|
||||
|
||||
|
||||
CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
|
||||
|
||||
};
|
||||
@ -250,28 +276,18 @@ namespace Grid {
|
||||
CartesianCommunicator::BroadcastWorld(0,(void *)&l,sizeof(l));
|
||||
}
|
||||
|
||||
template<class source> void Seed(source &src)
|
||||
{
|
||||
_generators[0] = RngEngine(src);
|
||||
_seeded=1;
|
||||
}
|
||||
void SeedRandomDevice(void){
|
||||
deviceSeed src;
|
||||
Seed(src);
|
||||
}
|
||||
void SeedFixedIntegers(const std::vector<int> &seeds){
|
||||
CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());
|
||||
fixedSeed src(seeds);
|
||||
Seed(src);
|
||||
std::seed_seq src(seeds.begin(),seeds.end());
|
||||
Seed(src,0);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
class GridParallelRNG : public GridRNGbase {
|
||||
public:
|
||||
|
||||
GridBase *_grid;
|
||||
int _vol;
|
||||
public:
|
||||
|
||||
int generator_idx(int os,int is){
|
||||
return is*_grid->oSites()+os;
|
||||
@ -285,15 +301,9 @@ namespace Grid {
|
||||
_uniform.resize(_vol,std::uniform_real_distribution<RealD>{0,1});
|
||||
_gaussian.resize(_vol,std::normal_distribution<RealD>(0.0,1.0) );
|
||||
_bernoulli.resize(_vol,std::discrete_distribution<int32_t>{1,1});
|
||||
_seeded=0;
|
||||
_uid.resize(_vol,std::uniform_int_distribution<uint32_t>() );
|
||||
}
|
||||
|
||||
|
||||
|
||||
//FIXME implement generic IO and create state save/restore
|
||||
//void SaveState(const std::string<char> &file);
|
||||
//void LoadState(const std::string<char> &file);
|
||||
|
||||
template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist){
|
||||
|
||||
typedef typename vobj::scalar_object scalar_object;
|
||||
@ -306,7 +316,6 @@ namespace Grid {
|
||||
int osites=_grid->oSites();
|
||||
int words=sizeof(scalar_object)/sizeof(scalar_type);
|
||||
|
||||
|
||||
parallel_for(int ss=0;ss<osites;ss++){
|
||||
|
||||
std::vector<scalar_object> buf(Nsimd);
|
||||
@ -329,104 +338,114 @@ namespace Grid {
|
||||
}
|
||||
};
|
||||
|
||||
// This loop could be made faster to avoid the Ahmdahl by
|
||||
// i) seed generators on each timeslice, for x=y=z=0;
|
||||
// ii) seed generators on each z for x=y=0
|
||||
// iii)seed generators on each y,z for x=0
|
||||
// iv) seed generators on each y,z,x
|
||||
// made possible by physical indexing.
|
||||
template<class source> void Seed(source &src)
|
||||
{
|
||||
void SeedFixedIntegers(const std::vector<int> &seeds){
|
||||
|
||||
typedef typename source::result_type seed_t;
|
||||
std::uniform_int_distribution<seed_t> uid;
|
||||
// Everyone generates the same seed_seq based on input seeds
|
||||
CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());
|
||||
|
||||
int numseed=4;
|
||||
int gsites = _grid->_gsites;
|
||||
std::vector<seed_t> site_init(numseed);
|
||||
std::seed_seq source(seeds.begin(),seeds.end());
|
||||
|
||||
RngEngine master_engine(source);
|
||||
|
||||
#ifdef RNG_FAST_DISCARD
|
||||
////////////////////////////////////////////////
|
||||
// Skip ahead through a single stream.
|
||||
// Applicable to SITMO and other has based/crypto RNGs
|
||||
// Should be applicable to Mersenne Twister, but the C++11
|
||||
// MT implementation does not implement fast discard even though
|
||||
// in principle this is possible
|
||||
////////////////////////////////////////////////
|
||||
std::vector<int> gcoor;
|
||||
int rank,o_idx,i_idx;
|
||||
|
||||
// Everybody loops over global volume.
|
||||
for(int gidx=0;gidx<_grid->_gsites;gidx++){
|
||||
|
||||
// Master RngEngine
|
||||
std::vector<seed_t> master_init(numseed); src.generate(master_init.begin(),master_init.end());
|
||||
_grid->Broadcast(0,(void *)&master_init[0],sizeof(seed_t)*numseed);
|
||||
fixedSeed master_seed(master_init);
|
||||
RngEngine master_engine(master_seed);
|
||||
|
||||
// Per node RngEngine
|
||||
std::vector<seed_t> node_init(numseed);
|
||||
for(int r=0;r<_grid->ProcessorCount();r++) {
|
||||
|
||||
std::vector<seed_t> rank_init(numseed);
|
||||
for(int i=0;i<numseed;i++) rank_init[i] = uid(master_engine);
|
||||
|
||||
std::cout << GridLogMessage << "SeedSeq for rank "<<r;
|
||||
for(int i=0;i<numseed;i++) std::cout<<" "<<rank_init[i];
|
||||
std::cout <<std::endl;
|
||||
|
||||
if ( r==_grid->ThisRank() ) {
|
||||
for(int i=0;i<numseed;i++) node_init[i] = rank_init[i];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////
|
||||
// Set up a seed_seq wrapper with these 8 words
|
||||
// and draw for each site within node.
|
||||
////////////////////////////////////////////////////
|
||||
fixedSeed node_seed(node_init);
|
||||
RngEngine node_engine(node_seed);
|
||||
|
||||
for(int gidx=0;gidx<gsites;gidx++){
|
||||
int rank,o_idx,i_idx;
|
||||
Skip(master_engine); // Skip to next RNG sequence
|
||||
|
||||
// Where is it?
|
||||
_grid->GlobalIndexToGlobalCoor(gidx,gcoor);
|
||||
_grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
|
||||
|
||||
// If this is one of mine we take it
|
||||
if( rank == _grid->ThisRank() ){
|
||||
int l_idx=generator_idx(o_idx,i_idx);
|
||||
for(int i=0;i<numseed;i++) site_init[i] = uid(node_engine);
|
||||
fixedSeed site_seed(site_init);
|
||||
_generators[l_idx] = RngEngine(site_seed);
|
||||
_generators[l_idx] = master_engine;
|
||||
}
|
||||
|
||||
}
|
||||
#else
|
||||
////////////////////////////////////////////////////////////////
|
||||
// Machine and thread decomposition dependent seeding is efficient
|
||||
// and maximally parallel; but NOT reproducible from machine to machine.
|
||||
// Not ideal, but fastest way to reseed all nodes.
|
||||
////////////////////////////////////////////////////////////////
|
||||
{
|
||||
// Obtain one Reseed per processor
|
||||
int Nproc = _grid->ProcessorCount();
|
||||
std::vector<RngEngine> seeders(Nproc);
|
||||
int me= _grid->ThisRank();
|
||||
for(int p=0;p<Nproc;p++){
|
||||
seeders[p] = Reseed(master_engine);
|
||||
}
|
||||
master_engine = seeders[me];
|
||||
}
|
||||
|
||||
{
|
||||
// Obtain one reseeded generator per thread
|
||||
int Nthread = GridThread::GetThreads();
|
||||
std::vector<RngEngine> seeders(Nthread);
|
||||
for(int t=0;t<Nthread;t++){
|
||||
seeders[t] = Reseed(master_engine);
|
||||
}
|
||||
|
||||
parallel_for(int t=0;t<Nthread;t++) {
|
||||
// set up one per local site in threaded fashion
|
||||
std::vector<uint32_t> newseeds;
|
||||
std::uniform_int_distribution<uint32_t> uid;
|
||||
for(int l=0;l<_grid->lSites();l++) {
|
||||
if ( (l%Nthread)==t ) {
|
||||
_generators[l] = Reseed(seeders[t],newseeds,uid);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
_seeded=1;
|
||||
}
|
||||
void SeedRandomDevice(void){
|
||||
deviceSeed src;
|
||||
Seed(src);
|
||||
#endif
|
||||
}
|
||||
void SeedFixedIntegers(const std::vector<int> &seeds){
|
||||
CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());
|
||||
fixedSeed src(seeds);
|
||||
Seed(src);
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Support for rigorous test of RNG's
|
||||
// Return uniform random uint32_t from requested site generator
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
uint32_t GlobalU01(int gsite){
|
||||
|
||||
uint32_t the_number;
|
||||
|
||||
// who
|
||||
std::vector<int> gcoor;
|
||||
int rank,o_idx,i_idx;
|
||||
_grid->GlobalIndexToGlobalCoor(gsite,gcoor);
|
||||
_grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
|
||||
|
||||
// draw
|
||||
int l_idx=generator_idx(o_idx,i_idx);
|
||||
if( rank == _grid->ThisRank() ){
|
||||
the_number = _uid[l_idx](_generators[l_idx]);
|
||||
}
|
||||
|
||||
// share & return
|
||||
_grid->Broadcast(rank,(void *)&the_number,sizeof(the_number));
|
||||
return the_number;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
template <class vobj> inline void random(GridParallelRNG &rng,Lattice<vobj> &l){
|
||||
rng.fill(l,rng._uniform);
|
||||
}
|
||||
template <class vobj> inline void random(GridParallelRNG &rng,Lattice<vobj> &l) { rng.fill(l,rng._uniform); }
|
||||
template <class vobj> inline void gaussian(GridParallelRNG &rng,Lattice<vobj> &l) { rng.fill(l,rng._gaussian); }
|
||||
template <class vobj> inline void bernoulli(GridParallelRNG &rng,Lattice<vobj> &l){ rng.fill(l,rng._bernoulli);}
|
||||
|
||||
template <class vobj> inline void gaussian(GridParallelRNG &rng,Lattice<vobj> &l){
|
||||
rng.fill(l,rng._gaussian);
|
||||
}
|
||||
|
||||
template <class vobj> inline void bernoulli(GridParallelRNG &rng,Lattice<vobj> &l){
|
||||
rng.fill(l,rng._bernoulli);
|
||||
}
|
||||
|
||||
template <class sobj> inline void random(GridSerialRNG &rng,sobj &l){
|
||||
rng.fill(l,rng._uniform);
|
||||
}
|
||||
|
||||
template <class sobj> inline void gaussian(GridSerialRNG &rng,sobj &l){
|
||||
rng.fill(l,rng._gaussian);
|
||||
}
|
||||
|
||||
template <class sobj> inline void bernoulli(GridSerialRNG &rng,sobj &l){
|
||||
rng.fill(l,rng._bernoulli);
|
||||
}
|
||||
template <class sobj> inline void random(GridSerialRNG &rng,sobj &l) { rng.fill(l,rng._uniform ); }
|
||||
template <class sobj> inline void gaussian(GridSerialRNG &rng,sobj &l) { rng.fill(l,rng._gaussian ); }
|
||||
template <class sobj> inline void bernoulli(GridSerialRNG &rng,sobj &l){ rng.fill(l,rng._bernoulli); }
|
||||
|
||||
}
|
||||
#endif
|
||||
|
@ -491,10 +491,15 @@ static inline void writeRNGState(GridSerialRNG &serial,GridParallelRNG ¶llel
|
||||
#ifdef RNG_RANLUX
|
||||
header.floating_point = std::string("UINT64");
|
||||
header.data_type = std::string("RANLUX48");
|
||||
#else
|
||||
#endif
|
||||
#ifdef RNG_MT19937
|
||||
header.floating_point = std::string("UINT32");
|
||||
header.data_type = std::string("MT19937");
|
||||
#endif
|
||||
#ifdef RNG_SITMO
|
||||
header.floating_point = std::string("UINT64");
|
||||
header.data_type = std::string("SITMO");
|
||||
#endif
|
||||
|
||||
truncate(file);
|
||||
offset = writeHeader(header,file);
|
||||
@ -522,10 +527,15 @@ static inline void readRNGState(GridSerialRNG &serial,GridParallelRNG & parallel
|
||||
#ifdef RNG_RANLUX
|
||||
assert(format == std::string("UINT64"));
|
||||
assert(data_type == std::string("RANLUX48"));
|
||||
#else
|
||||
#endif
|
||||
#ifdef RNG_MT19937
|
||||
assert(format == std::string("UINT32"));
|
||||
assert(data_type == std::string("MT19937"));
|
||||
#endif
|
||||
#ifdef RNG_SITMO
|
||||
assert(format == std::string("UINT64"));
|
||||
assert(data_type == std::string("SITMO"));
|
||||
#endif
|
||||
|
||||
// depending on datatype, set up munger;
|
||||
// munger is a function of <floating point, Real, data_type>
|
||||
|
@ -58,6 +58,7 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
|
||||
#include <Grid/qcd/action/fermion/DomainWallFermion.h>
|
||||
#include <Grid/qcd/action/fermion/MobiusFermion.h>
|
||||
#include <Grid/qcd/action/fermion/ZMobiusFermion.h>
|
||||
#include <Grid/qcd/action/fermion/SchurDiagTwoKappa.h>
|
||||
#include <Grid/qcd/action/fermion/ScaledShamirFermion.h>
|
||||
#include <Grid/qcd/action/fermion/MobiusZolotarevFermion.h>
|
||||
#include <Grid/qcd/action/fermion/ShamirZolotarevFermion.h>
|
||||
|
102
lib/qcd/action/fermion/SchurDiagTwoKappa.h
Normal file
102
lib/qcd/action/fermion/SchurDiagTwoKappa.h
Normal file
@ -0,0 +1,102 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: SchurDiagTwoKappa.h
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Christoph Lehner
|
||||
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 */
|
||||
#ifndef _SCHUR_DIAG_TWO_KAPPA_H
|
||||
#define _SCHUR_DIAG_TWO_KAPPA_H
|
||||
|
||||
namespace Grid {
|
||||
|
||||
// This is specific to (Z)mobius fermions
|
||||
template<class Matrix, class Field>
|
||||
class KappaSimilarityTransform {
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Matrix);
|
||||
std::vector<Coeff_t> kappa, kappaDag, kappaInv, kappaInvDag;
|
||||
|
||||
KappaSimilarityTransform (Matrix &zmob) {
|
||||
for (int i=0;i<(int)zmob.bs.size();i++) {
|
||||
Coeff_t k = 1.0 / ( 2.0 * (zmob.bs[i] *(4 - zmob.M5) + 1.0) );
|
||||
kappa.push_back( k );
|
||||
kappaDag.push_back( conj(k) );
|
||||
kappaInv.push_back( 1.0 / k );
|
||||
kappaInvDag.push_back( 1.0 / conj(k) );
|
||||
}
|
||||
}
|
||||
|
||||
template<typename vobj>
|
||||
void sscale(const Lattice<vobj>& in, Lattice<vobj>& out, Coeff_t* s) {
|
||||
GridBase *grid=out._grid;
|
||||
out.checkerboard = in.checkerboard;
|
||||
assert(grid->_simd_layout[0] == 1); // should be fine for ZMobius for now
|
||||
int Ls = grid->_rdimensions[0];
|
||||
parallel_for(int ss=0;ss<grid->oSites();ss++){
|
||||
vobj tmp = s[ss % Ls]*in._odata[ss];
|
||||
vstream(out._odata[ss],tmp);
|
||||
}
|
||||
}
|
||||
|
||||
RealD sscale_norm(const Field& in, Field& out, Coeff_t* s) {
|
||||
sscale(in,out,s);
|
||||
return norm2(out);
|
||||
}
|
||||
|
||||
virtual RealD M (const Field& in, Field& out) { return sscale_norm(in,out,&kappa[0]); }
|
||||
virtual RealD MDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaDag[0]);}
|
||||
virtual RealD MInv (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInv[0]);}
|
||||
virtual RealD MInvDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInvDag[0]);}
|
||||
|
||||
};
|
||||
|
||||
template<class Matrix,class Field>
|
||||
class SchurDiagTwoKappaOperator : public SchurOperatorBase<Field> {
|
||||
public:
|
||||
KappaSimilarityTransform<Matrix, Field> _S;
|
||||
SchurDiagTwoOperator<Matrix, Field> _Mat;
|
||||
|
||||
SchurDiagTwoKappaOperator (Matrix &Mat): _S(Mat), _Mat(Mat) {};
|
||||
|
||||
virtual RealD Mpc (const Field &in, Field &out) {
|
||||
Field tmp(in._grid);
|
||||
|
||||
_S.MInv(in,out);
|
||||
_Mat.Mpc(out,tmp);
|
||||
return _S.M(tmp,out);
|
||||
|
||||
}
|
||||
virtual RealD MpcDag (const Field &in, Field &out){
|
||||
Field tmp(in._grid);
|
||||
|
||||
_S.MDag(in,out);
|
||||
_Mat.MpcDag(out,tmp);
|
||||
return _S.MInvDag(tmp,out);
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
@ -54,7 +54,7 @@ THE SOFTWARE.
|
||||
|
||||
#define GRID_MACRO_EMPTY()
|
||||
|
||||
#define GRID_MACRO_EVAL(...) GRID_MACRO_EVAL1024(__VA_ARGS__)
|
||||
#define GRID_MACRO_EVAL(...) GRID_MACRO_EVAL64(__VA_ARGS__)
|
||||
#define GRID_MACRO_EVAL1024(...) GRID_MACRO_EVAL512(GRID_MACRO_EVAL512(__VA_ARGS__))
|
||||
#define GRID_MACRO_EVAL512(...) GRID_MACRO_EVAL256(GRID_MACRO_EVAL256(__VA_ARGS__))
|
||||
#define GRID_MACRO_EVAL256(...) GRID_MACRO_EVAL128(GRID_MACRO_EVAL128(__VA_ARGS__))
|
||||
|
Reference in New Issue
Block a user