If Field<T> a belongs to original (parent) lattice and Field<T> belongs to blocked lattice, the operation
copies data from the sparse (blocked) sites of a to b.
The sites which are blocked can be accessed with the CoordinateVector function is_divisible:
#ifndef HILA_FIELD_H
#define HILA_FIELD_H
#include <sstream>
#include <iostream>
#include <string>
#include <cstring>
#include <math.h>
#include <type_traits>
#include "plumbing/lattice.h"
#include "plumbing/backend_vector/vector_types.h"
#include "plumbing/com_mpi.h"
#define onsites(p) for (Parity par_dummy__(p); par_dummy__ == EVEN; par_dummy__ = ODD)
template <typename T>
template <typename T>
void ensure_field_operators_exist();
#include "plumbing/ensure_loop_functions.h"
template <typename T>
public:
enum class gather_status_t : unsigned { NOT_DONE, STARTED, DONE };
private:
class field_struct {
public:
#ifdef VECTORIZED
#endif
unsigned assigned_to;
gather_status_t gather_status_arr[3][
NDIRS];
MPI_Request receive_request[3][
NDIRS];
MPI_Request send_request[3][
NDIRS];
#ifndef VANILLA
T *receive_buffer[
NDIRS];
#endif
void initialize_communication() {
for (
int d = 0; d <
NDIRS; d++) {
for (int p = 0; p < 3; p++)
gather_status_arr[p][d] = gather_status_t::NOT_DONE;
send_buffer[d] = nullptr;
#ifndef VANILLA
receive_buffer[d] = nullptr;
#endif
}
}
void free_communication() {
for (
int d = 0; d <
NDIRS; d++) {
if (send_buffer[d] != nullptr)
payload.
free_mpi_buffer(send_buffer[d]);
#ifndef VANILLA
if (receive_buffer[d] != nullptr)
payload.free_mpi_buffer(receive_buffer[d]);
#endif
}
}
void allocate_payload() {
payload.
allocate_field(lattice);
}
void free_payload() {
}
#ifndef VECTORIZED
inline auto get(const unsigned i) const {
return payload.
get(i, lattice->
mynode.
field_alloc_size);
}
template <typename A>
inline void set(const A &value, const unsigned i) {
payload.
set(value, i, lattice->mynode.field_alloc_size);
}
inline auto get_element(const unsigned i) const {
}
template <typename A>
inline void set_element(const A &value, const unsigned i) {
payload.
set_element(value, i, lattice);
}
#else
template <typename vecT>
inline vecT get_vector(const unsigned i) const {
return payload.template get_vector<vecT>(i);
}
inline T get_element(const unsigned i) const {
}
template <typename vecT>
inline void set_vector(const vecT &val, const unsigned i) {
return payload.set_vector(val, i);
}
inline void set_element(const T &val, const unsigned i) {
return payload.set_element(val, i);
}
#endif
void gather_elements(T *buffer, const std::vector<CoordinateVector> &coord_list,
int root = 0) const;
void scatter_elements(T *buffer, const std::vector<CoordinateVector> &coord_list,
int root = 0);
};
static_assert(std::is_trivial<T>::value && std::is_standard_layout<T>::value,
"Field expects only pod-type elements (plain data): default "
"constructor, copy and delete");
public:
#ifdef VECTORIZED
static_assert(sizeof(hila::arithmetic_type<T>) == 4 ||
sizeof(hila::arithmetic_type<T>) == 8,
"In vectorized arch (e.g. AVX2), only 4 or 8 byte (32 or 64 bit) numbers for "
"Field<> implemented, sorry!");
#endif
#if defined(CUDA) || defined(HIP)
static_assert(!std::is_same<hila::arithmetic_type<T>, long double>::value,
"Type 'long double' numbers in Field<> not supported by cuda/hip");
#endif
}
}
template <typename A, std::enable_if_t<std::is_convertible<A, T>::value, int> = 0>
}
template <typename A,
std::enable_if_t<
hila::is_assignable<T &, A>::value || std::is_convertible<A, T>::value, int> = 0>
}
}
assert(rhs.fs->mylattice.ptr() == lattice.
ptr());
rhs.fs = nullptr;
}
#ifdef HILAPP
ensure_field_operators_exist<T>();
#endif
}
void allocate() {
assert(lattice.
is_initialized() &&
"Fields cannot be used before lattice.setup()");
fs = (field_struct *)memalloc(
sizeof(field_struct));
fs->initialize_communication();
#if !defined(CUDA) && !defined(HIP)
#else
fs->payload.neighbours[d] = lattice->
backend_lattice->
d_neighb[d];
#endif
}
#ifdef SPECIAL_BOUNDARY_CONDITIONS
fs->boundary_condition[dir] = hila::bc::PERIODIC;
fs->boundary_condition[-dir] = hila::bc::PERIODIC;
}
#endif
#ifdef VECTORIZED
fs->vector_lattice = lattice->backend_lattice
} else {
fs->vector_lattice =
nullptr;
}
#endif
}
if (
fs !=
nullptr && !hila::about_to_finish) {
fs->free_communication();
}
}
}
return fs !=
nullptr && ((
fs->assigned_to & parity_bits(p)) != 0);
}
gather_status_t gather_status(
Parity p,
int d)
const {
assert(parity_bits(p) && d >= 0 && d <
NDIRS);
return fs->gather_status_arr[(int)p - 1][d];
}
void set_gather_status(
Parity p,
int d, gather_status_t stat)
const {
assert(parity_bits(p) && d >= 0 && d <
NDIRS);
fs->gather_status_arr[(int)p - 1][d] = stat;
}
allocate();
}
}
void mark_changed(
const Parity p)
const {
drop_comms(i, opp_parity(p));
set_gather_status(opp_parity(p), i, gather_status_t::NOT_DONE);
set_gather_status(
ALL, i, gather_status_t::NOT_DONE);
} else {
set_gather_status(
EVEN, i, gather_status_t::NOT_DONE);
set_gather_status(
ODD, i, gather_status_t::NOT_DONE);
}
}
fs->assigned_to |= parity_bits(p);
}
void mark_gathered(
int dir,
const Parity p)
const {
set_gather_status(p, dir, gather_status_t::DONE);
}
bool is_gathered(
int dir,
Parity par)
const {
return gather_status(par, dir) == gather_status_t::DONE ||
gather_status(
ALL, dir) == gather_status_t::DONE;
} else {
return gather_status(
ALL, dir) == gather_status_t::DONE ||
(gather_status(
EVEN, dir) == gather_status_t::DONE &&
gather_status(
ODD, dir) == gather_status_t::DONE);
}
}
void mark_gather_started(
int dir,
Parity p)
const {
set_gather_status(p, dir, gather_status_t::STARTED);
}
bool is_gather_started(
int dir,
Parity par)
const {
return gather_status(par, dir) == gather_status_t::STARTED;
}
bool gather_not_done(
int dir,
Parity par)
const {
return gather_status(par, dir) == gather_status_t::NOT_DONE;
}
bool boundary_need_to_communicate(
const Direction dir)
const {
#ifdef SPECIAL_BOUNDARY_CONDITIONS
!
fs->mylattice->mynode.is_on_edge(dir);
#else
return true;
#endif
}
#ifdef SPECIAL_BOUNDARY_CONDITIONS
#ifdef HILAPP
"BC possible only for types which implement unary "
"minus (-) -operator and are not unsigned");
ensure_unary_minus_is_loop_function<T>();
#endif
fs->boundary_condition[dir] =
bc;
fs->boundary_condition[-dir] =
bc;
#if !defined(CUDA) && !defined(HIP)
fs->neighbours[dir] = mylat->
get_neighbour_array(dir, bc);
fs->neighbours[-dir] = mylat->get_neighbour_array(-dir, bc);
#else
if (bc == hila::bc::PERIODIC) {
fs->payload.neighbours[dir] = mylat->backend_lattice->
d_neighb[dir];
fs->payload.neighbours[-dir] = mylat->backend_lattice->
d_neighb[-dir];
} else {
fs->payload.neighbours[dir] = mylat->backend_lattice->d_neighb_special[dir];
fs->payload.neighbours[-dir] = mylat->backend_lattice->d_neighb_special[-dir];
}
#endif
#else
assert(bc == hila::bc::PERIODIC &&
"Only periodic bondary conditions when SPECIAL_BOUNDARY_CONDITIONS is undefined");
#endif
}
#ifdef SPECIAL_BOUNDARY_CONDITIONS
return fs->boundary_condition[dir];
#else
return hila::bc::PERIODIC;
#endif
}
void print_boundary_condition() {
for (
int dir = 0; dir <
NDIRS; dir++) {
}
}
template <typename A>
}
}
T operator[](
const Parity p)
const;
T &operator[](
const Parity p);
inline auto field_buffer() const {
return fs->payload.get_buffer();
}
#ifndef VECTORIZED
return fs->get_element(i);
}
#else
return fs->get_element(i);
}
template <typename vecT>
inline auto get_vector_at(unsigned i) const {
return fs->template get_vector<vecT>(i);
}
inline auto get_value_at_nb_site(
Direction d,
unsigned i)
const {
return fs->get_element(
fs->vector_lattice->site_neighbour(d, i));
}
#endif
#ifndef VECTORIZED
template <typename A>
fs->set_element(value, i);
}
#else
template <typename vecT>
inline void set_vector_at(const vecT &value, unsigned i) {
fs->set_vector(value, i);
}
template <typename A>
fs->set_element(value, i);
}
#endif
return *this;
}
template <typename A,
std::enable_if_t<
hila::is_assignable<T &, A>::value || std::is_convertible<A, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<
hila::is_assignable<T &, A>::value || std::is_convertible<A, T>::value, int> = 0>
return *this;
}
return *this;
}
if (this != &rhs) {
assert((
fs ==
nullptr || (rhs.fs->mylattice.ptr() ==
fs->mylattice.ptr())) &&
"Field = Field assignment possible only if fields belong to the same lattice");
rhs.fs = nullptr;
}
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_plus<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_minus<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_mul<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_div<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_plus<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_minus<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_mul<T, A>, T>::value, int> = 0>
return *this;
}
template <typename A,
std::enable_if_t<std::is_convertible<hila::type_div<T, A>, T>::value, int> = 0>
return *this;
}
return *this;
}
return f;
}
template <typename S>
double s = 0;
onsites(
ALL) s += !((*this)[X] == rhs[X]);
return (s == 0);
}
template <typename S>
bool operator!=(
const Field<S> &rhs)
const {
return !(*this == rhs);
}
double n = 0;
}
return n;
}
}
return f;
}
template <typename R = T, typename A = decltype(::dagger(std::declval<R>()))>
return f;
}
template <typename R = T, typename A = decltype(::real(std::declval<R>()))>
return f;
}
template <typename R = T, typename A = decltype(::imag(std::declval<R>()))>
return f;
}
}
template <typename A, std::enable_if_t<std::is_assignable<T &, A>::value, int> = 0>
T element;
element = value;
if (
fs->mylattice->is_on_mynode(coord)) {
}
mark_changed(coord.
parity());
return element;
}
T element;
int owner =
fs->mylattice->node_rank(coord);
}
return element;
}
void set_elements(const std::vector<T> &elements,
const std::vector<CoordinateVector> &coord_list);
std::vector<T>
get_elements(
const std::vector<CoordinateVector> &coord_list,
bool broadcast = false) const;
bool broadcast = false) const;
void copy_local_data(std::vector<T> &buffer) const;
void set_local_data(const std::vector<T> &buffer);
void copy_local_data_with_halo(std::vector<T> &buffer) const;
template <typename A,
std::enable_if_t<std::is_assignable<T &, hila::type_plus<T, A>>::value, int> = 0>
if (
fs->mylattice->is_on_mynode(coord)) {
auto i =
fs->mylattice->site_index(coord);
v += av;
}
mark_changed(coord.parity());
}
template <typename A,
std::enable_if_t<std::is_assignable<T &, hila::type_minus<T, A>>::value, int> = 0>
if (
fs->mylattice->is_on_mynode(coord)) {
auto i =
fs->mylattice->site_index(coord);
v -= av;
}
mark_changed(coord.parity());
}
template <typename A,
std::enable_if_t<std::is_assignable<T &, hila::type_mul<T, A>>::value, int> = 0>
if (
fs->mylattice->is_on_mynode(coord)) {
auto i =
fs->mylattice->site_index(coord);
v *= av;
}
mark_changed(coord.parity());
}
template <typename A,
std::enable_if_t<std::is_assignable<T &, hila::type_div<T, A>>::value, int> = 0>
if (
fs->mylattice->is_on_mynode(coord)) {
auto i =
fs->mylattice->site_index(coord);
v /= av;
}
mark_changed(coord.parity());
}
T r, v;
r = v;
v++;
return r;
}
++v;
return v;
}
T r, v;
r = v;
v--;
return r;
}
--v;
return v;
}
FFT_complex_to_real(
fft_direction fdir = fft_direction::forward)
const;
void write(std::ofstream &outputfile,
bool binary =
true,
int precision = 8)
const;
void write(
const std::string &filename,
bool binary =
true,
int precision = 8)
const;
void read(std::ifstream &inputfile);
void read(
const std::string &filename);
template <typename Out>
void write_slice(Out &outputfile,
const CoordinateVector &slice,
int precision = 6)
const;
T
sum(
Parity par = Parity::all,
bool allreduce =
true)
const;
T
product(
Parity par = Parity::all,
bool allreduce =
true)
const;
void random();
void gaussian_random(double width = 1.0);
}
};
template <typename A, typename B>
tmp[
ALL] = lhs[X] + rhs[X];
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
return rhs + lhs;
}
template <typename A, typename B>
tmp[
ALL] = lhs[X] - rhs[X];
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
tmp[
ALL] = lhs[X] * rhs[X];
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename A, typename B>
return tmp;
}
template <typename T>
}
}
template <typename T, typename R = decltype(exp(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(log(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(sin(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(cos(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(tan(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(asin(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(acos(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(atan(std::declval<T>()))>
}
return res;
}
template <typename T, typename R = decltype(abs(std::declval<T>()))>
}
return res;
}
template <typename T, typename P, typename R = decltype(pow(std::declval<T>()), std::declval<P>())>
}
return res;
}
template <typename T>
double r = 0;
}
return r;
}
template <typename T>
}
template <typename T>
}
template <typename T, typename A = decltype(::dagger(std::declval<T>()))>
}
template <typename T, typename A = decltype(::real(std::declval<T>()))>
}
template <typename T, typename A = decltype(::imag(std::declval<T>()))>
}
template <typename A, typename B, typename R = decltype(std::declval<A>() - std::declval<B>())>
double res = 0;
}
return res;
}
template <typename T>
return res;
}
template <typename T>
if (hila::is_comm_initialized()) {
if (is_gather_started(d,
ALL)) {
drop_comms_timer.
start();
}
if (is_gather_started(d, p)) {
drop_comms_timer.start();
wait_gather(d, p);
drop_comms_timer.stop();
}
} else {
if (is_gather_started(d,
EVEN)) {
drop_comms_timer.start();
drop_comms_timer.stop();
}
if (is_gather_started(d,
ODD)) {
drop_comms_timer.start();
drop_comms_timer.stop();
}
}
}
}
template <typename T>
start_gather(d, p);
wait_gather(d, p);
}
template <typename T>
#if defined(CUDA) || defined(HIP)
std::vector<T> rng_buffer(lattice->mynode.
volume);
for (auto &element : rng_buffer)
(*this).set_local_data(rng_buffer);
} else {
}
}
#else
}
#endif
}
template <typename T>
#if defined(CUDA) || defined(HIP)
std::vector<T> rng_buffer(lattice->mynode.volume);
for (auto &element : rng_buffer)
(*this).set_local_data(rng_buffer);
} else {
}
}
#else
}
#endif
}
#ifdef HILAPP
inline void dummy_X_f() {
d = +d1;
d = -d1;
vec = +v1;
vec = -v1;
vec = d + d1;
vec = d - d1;
vec = v1 + d1;
vec = v1 - d1;
vec = d1 + v1;
vec = d1 - v1;
vec = vec + v1;
vec = vec - v1;
vec[e_x] = vec[0] + vec[e_y] + vec.
e(e_y);
}
}
template <typename T>
inline void ensure_field_operators_exist() {
ensure_unary_minus_is_loop_function<T>();
ensure_assign_zero_is_loop_function<T>();
}
#endif
#endif
Array< n, m, T > conj(const Array< n, m, T > &arg)
Return conjugate Array.
Array< n, m, hila::arithmetic_type< T > > imag(const Array< n, m, T > &arg)
Return imaginary part of Array.
hila::arithmetic_type< T > squarenorm(const Array< n, m, T > &rhs)
Return square norm of Array.
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Field< T > & operator=(const Field< T > &rhs)
Assignment operator.
Field()
Field constructor.
Field< A > real() const
Returns real part of Field.
std::vector< T > get_subvolume(const CoordinateVector &cmin, const CoordinateVector &cmax, bool broadcast=false) const
Get a subvolume of the field elements to all nodes.
Field< T > & shift(const CoordinateVector &v, Field< T > &r, Parity par) const
Create a periodically shifted copy of the field.
void set_boundary_condition(Direction dir, hila::bc bc)
Set the boundary condition in a given Direction (periodic or antiperiodic)
void write_subvolume(std::ofstream &outputfile, const CoordinateVector &cmin, const CoordinateVector &cmax, int precision=6) const
const T get_element(const CoordinateVector &coord) const
Get singular element which will be broadcast to all nodes.
hila::bc get_boundary_condition(Direction dir) const
Get the boundary condition of the Field.
Field< A > imag() const
Returns imaginary part of Field.
Lattice mylattice() const
Return the lattice to which this field instance belongs to.
T gpu_minmax(bool min_or_max, Parity par, CoordinateVector &loc) const
Declare gpu_reduce here, defined only for GPU targets.
Field< A > dagger() const
Returns dagger or Hermitian conjugate of Field depending on how it is defined for Field type T.
const T set_element(const CoordinateVector &coord, const A &value)
Get singular element which will be broadcast to all nodes.
T max(Parity par=ALL) const
Find maximum value from Field.
Field< T > operator+() const
Summation operator.
bool is_initialized(Parity p) const
Returns true if the Field has been assigned to.
auto get_value_at(const unsigned i) const
Get an individual element outside a loop. This is also used as a getter in the vanilla code.
field_struct *__restrict__ fs
Field::fs holds all of the field content in Field::field_struct.
void clear()
Destroys field data.
void set_value_at(const A &value, unsigned i)
Set an individual element outside a loop. This is also used as a getter in the vanilla code.
Field< T > operator-() const
Subtraction operator.
Field< Complex< hila::arithmetic_type< T > > > FFT_real_to_complex(fft_direction fdir=fft_direction::forward) const
Field< T > & operator*=(const Field< A > &rhs)
Product assignment operator.
bool is_allocated() const
Returns true if the Field data has been allocated.
void copy_boundary_condition(const Field< A > &rhs)
Copy the boundary condition from another field.
Field< T > conj() const
Returns field with all elements conjugated depending how conjugate is defined for type.
T product(Parity par=Parity::all, bool allreduce=true) const
Product reduction of Field.
void check_alloc()
Allocate Field if it is not already allocated.
bool operator==(const Field< S > &rhs) const
Field comparison operator.
double squarenorm() const
Squarenorm.
dir_mask_t start_gather(Direction d, Parity p=ALL) const
Field< T > & operator-=(const Field< A > &rhs)
Subtraction assignment operator.
void write(std::ofstream &outputfile, bool binary=true, int precision=8) const
Write the field to a file stream.
std::vector< T > get_elements(const std::vector< CoordinateVector > &coord_list, bool broadcast=false) const
Get a list of elements and store them into a vector.
T minmax(bool is_min, Parity par, CoordinateVector &loc) const
Function to perform min or max operations.
T min(Parity par=ALL) const
Find minimum value from Field.
T sum(Parity par=Parity::all, bool allreduce=true) const
Sum reduction of Field.
void read(std::ifstream &inputfile)
Read the Field from a stream.
std::vector< T > get_slice(const CoordinateVector &c, bool broadcast=false) const
Get a slice (subvolume)
void unblock_to(Field< T > &target) const
Copy content to the argument Field on blocked (sparse) sites. a.unblock_to(b) is the inverse of a....
Field< T > & operator/=(const Field< A > &rhs)
Division assignment operator.
Field< T > & operator+=(const Field< A > &rhs)
Addition assignment operator.
main interface class to lattices.
lattice_struct * ptr() const
get non-const pointer to lattice_struct (cf. operator ->)
T e(const int i, const int j) const
Standard array indexing operation for matrices.
X-coordinate type - "dummy" class.
The field_storage struct contains minimal information for using the field in a loop....
auto get_element(const unsigned i, const Lattice lattice) const
Conditionally reture bool type false if type T does't have unary - operator.
unsigned *__restrict__ neighb[NDIRS]
Main neighbour index array.
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
T arg(const Complex< T > &a)
Return argument of Complex number.
This header file defines:
Parity
Parity enum with values EVEN, ODD, ALL; refers to parity of the site. Parity of site (x,...
constexpr Parity EVEN
bit pattern: 001
#define foralldir(d)
Macro to loop over (all) Direction(s)
constexpr unsigned NDIRS
Number of directions.
constexpr Parity ODD
bit pattern: 010
Direction
Enumerator for direction that assigns integer to direction to be interpreted as unit vector.
This file defines all includes for HILA.
fft_direction
define a class for FFT direction
Implement hila::swap for gauge fields.
T gaussian_random()
Template function T hila::gaussian_random<T>(),generates gaussian random value of type T,...
logger_class log
Now declare the logger.
double random()
Real valued uniform random number generator.
int myrank()
rank of this node
bool is_device_rng_on()
Check if the RNG on GPU is allocated and ready to use.
std::ostream out0
This writes output only from main process (node 0)
bc
list of field boundary conditions - used only if SPECIAL_BOUNDARY_CONDITIONS defined
T broadcast(T &var, int rank=0)
Broadcast the value of var to all MPI ranks from rank (default=0).
bool bc_need_communication(hila::bc bc)
False if we have b.c. which does not require communication.
T gaussian_random(out_only T &val, double w=1.0)
Template function const T & hila::gaussian_random(T & variable,double width=1)
X + coordinate offset, used in f[X+CoordinateVector] or f[X+dir1+dir2] etc.
vectorized_lattice_struct< vector_size > * get_vectorized_lattice()
Returns a vectorized lattice with given vector size.
unsigned * d_neighb[NDIRS]
Storage for the neighbour indexes. Stored on device.
is_vectorizable_type<T>::value is always false if the target is not vectorizable
Information necessary to communicate with a node.