Field

Field is the most important Datatype offered by HILA. The Field is a container of lattice fields, and is the general object we evolve and iterate over. The Field can be comprised of either Standard types or Basic types listed on the datatypes page.

To see all the possible methods of a Field see the class page which lists comprehensive documentation.

Constructors and assignments

By default Field variables are constructed uninitialized. Fields can also be constructed with a constant or another Field:

Field<Complex<double>> f;       // Default constructor assigns the fields to `nullptr`
Field<Complex<double>> g(1);    // Each elements becomes (1,0)
Field<Complex<double>> h = 1;   // Equivalent to above
Field<Complex<double>> m(g);    // Copy constructor: content of g is copied
m = f;                          // ERROR: trying to use the value of the uninitialized Field f
f = 2 + g;                      // This initializes f

Assginment to Field variables is possible either from another Field or from a single element which are assignment compatible. Assignment from 0 is always possible.

Uninitialized Field does not take any storage, it is automatically allocated on first use. Fields are destructed when they run out of scope.

Access and Traversal

The principal traversal of the lattice is with site loops onsites(Parity), and a special location identifier X (effectively a new keyword). Within the onsites loop X location identifier represents the current location of a point that is currently being indexed. Note that the X identifier is only defined within onsites loops. Access operation f[X] can be applied only to field variables, and has the type of the field element. X is of type X_index_type. All available methods can be seen in the class documentation. Note that the X_index_type class is only a dummy class with decelerations, yet hilapp handles defining the contents of the methods, so the source is not available.

To illustrate how looping over a Field object works we will first define a few fields:

using mytype = Matrix<3,3,Complex<double>>;   // use type alias
Field<mytype> f,g,h; // Default constructor assigns the fields to `nullptr`
g = 2 // Assigning g to be 2*I throughout the field

For a field comprised of square-matrix elements, real numbers are algebraically interpreted as \( 2 = 2\cdot\mathbb{1}\), multiples of identity matrix.

We can now iterate over the fields with the onsites loop:

onsites(ALL) f[X] = 2 + g[X];          // 2 acts as 2* unit matrix for square matrices

Above we linearly add each element at X from g to each element at X in f with an additional \(2\cdot\mathbb{1}\) at each site. The ALL statement is a Parity which indicates that we will iterate over the whole Field. Other options are EVEN and ODD which indicate that we only iterate over the even or odd sites.

Similarly we can write this same statement in the following short form:

f[ALL] = 2 + g[X];                     // equivalent shorter form for simple 1-line assignments
f = 2 + g;                             // this is also equivalent!

Above you can also notice the simplest algebraic form, which allows for applying linear operations of the fields. The main difference is in sequencing: the first form goes through the lattice sites in one site loop, whereas the second stores the result of 2 + g to a temporary field variable which is copied to f (in this case std::moved). The site loop form is faster since it minimizes temporaries and memory accesses.

Now to demonstrate a more complicated onsites loop we will apply neighboring effects.

parity p = EVEN;
Direction d = e_x;
 
onsites(p) {
    auto t = g[X + d];                 // fetch from neighboring site in the e_x direction
    f[X] += t + t*t;                   // can define variables in the loop   
 
    h[X] = g[X + e_x - 2*e_y];         // non-nearest neighbour fetch
 
    if (X.coordinate(e_t) == 0) {      // Do this on 1st timeslice only
        h[X] *= 0.5;
    }
}

On the first line in the onsites loop we define a variable which we assign the e_x neighbor to. As we can see, variables can be defined within the scope of the loop.

Non nearest neighboring indexing also works, which is illustrated on the fourth line of the onsites loop.

On line 6 of the onsites loop we can also see that if statements can be used to apply limitations, in the above case we use it to index a slice of the field.

Because f[X] is of type field element (in this case mytype), the methods defined for the element type can be used. Within onsites loop f[X].dagger() is ok, f.dagger() is not. f[X] also serves as a visual identifier for a field variable access.

External non-Field variables cannot be changed inside onsites loops (except in reductions, see below)

double d = 2;
onsites(ALL) {
    d = f[X];              // ERROR: cannot change the value of variables defined outside onsites
    double a = f[X] + d;   // OK, a defined inside loop and d is constant
}

External variables

Assignment and manipulation of variables defined outside of onsites loops are illustrated below:

double a = 3, b = 5;
Field<double> f, g=0;
 
onsites(ALL) {
    f[X] = (a + b);            // ok, loop extern variables a,b do not change within the loop
    b = f[X];                  // ERROR: cannot change a loop external non-field variable (except reductions)
    double c = sin(f[X]);      // ok, variable c defined within the loop
    f[X] = c + g;              // ERROR: using field variable g without [X]
}

Loop external field access

There are a few ways of accessing Field elements which exclusively are to be used outside of onsites loops.

Access field at a single point:

CoordinateVector v = {2,3,4,5};
auto value = f[v];              // "value" is broadcast to all nodes!
f[v] = 1;                       // In assignment, values are not broadcast: the node which
                                // owns site v must have correct rhs.

Accessing multiple elements:

Accessing multiple elements is done by using Field::get_elements and Field::set_elements which take in a std::vector of coordinates.

std::vector<CoordinateVector> coordinates = {{2,3,4,5}, {1,3,4,5}}; // 4 Dimensional field
std::vector<MyType> values = f.get_elements(coordinates); // Returns std::vector of elements that were fetched
g.set_elements(values, coordinates); // Sets values fetched from f

Additionally Field::get_elements takes in an additional argument broadcast get_elements(coordinates, broadcast). If set to true then the fetched elements are broadcasted to all MPI ranks, by default it is set to false.

For getting a sub volume there are the following methods:

• Field::get_slice
• Field::get_subvolume

Reductions

Field reduction methods

Reduction on a single Field variable can be done with functions:

• T Field::sum() : sum of all field elements.
• T Field::product() : product of field elements.
• T Field::min() : minimum element
• T Field::max() : maximum

sum() and product() take optional arguments: sum(Parity par = ALL) or sum(Parity par = ALL ,bool allreduce = true). Parity gives the parity of sites over which the reduction is done, and if allreduce = false the reduction result is only sent to rank 0. Defaults are parity=ALL and allreduce = true.

Field<Complex<double>> f;
f.random();
hila::out0 << "Sum of elements of f = " << f.sum() << '\n';
hila::out0 << "Sum over even sites = " << f.sum(EVEN) << '\n';
// following result goes only to rank 0, but it is OK because
// the result is printed from that rank
hila::out0 << "Sum over odd sites = " << f.sum(ODD,false) << '\n';

min() and max() have signatures

• T Field<T>::min(Parity p) : minimum over Parity
• T Field<T>::min(CoordinateVector &loc) : also location of the min value
• T Field<T>::min(Parity p, CoordinateVector &loc) : location over parity

Field<double> f;
f.random();
CoordinateVector cv;
hila::out0 << "Max of f is " << f.max(cv) << " at location " << cv << '\n';

‍NOTE: sum() can be used for any allowed Field variable, product(), min() and max() only over integer or floating point Fields (C++ arithmetic types).

Standard reduction in onsites()

Standard reduction is done by using compound operators += or *= and regular variables in onsites()-loops:

Field<Complex<double>> f;
...
Complex<double> s2 = 0, s4 = 0;
onsites(ALL) {
    auto sqn = f[X].squarenorm();
    s2 += sqn;
    s4 += sqr(sqn);
}
hila::out0 << "Average of |f|^2 = " << s2 / lattice.volume() 
           << " and |f|^4 = " << s4 / lattice.volume() << '\n';

By default the reduction is sent to all MPI ranks (allreduce, implemented with MPI_Allreduce()). If reduction only to rank 0 is needed, hila::set_allreduce(false) can be inserted before onsites() (uses MPI_Reduce()). This may give higher performance.

hila::set_allreduce(false);
onsites(ALL) s += f[X];
// at this point only rank 0 contains valid reduction

Allreduce remains off only for the immediately following onsites()-loop, thus hila::set_allreduce(false) has to be inserted on each reduction loop where one wishes to have allreduce off.

Sum reduction += can use any hila Field element type, product reduction *= only floating point or integer variables (C++ arithmetic types).

Field<double> f;
...
double p = 1;  // note initialization to 1
onsites(ALL) p *= f[X];
// now p contains the product of all elements of f

Reductions can be accumulated using the same reduction variable:

double s = 5;        // initializing, for some reason, to 5
onsites(ALL) {
    s += f[X];
    if (g[X] > 0) 
        s += g[X];
}
onsites(EVEN) {
    s += h[X];
}
// s will now contain 5 + sum of operations above

As an implementation detail, the initial value of the reduction variable is taken from rank 0, i.e. it is zeroed on other ranks.

Individual data elements of a reduction variable can be accumulated independently:

Field<int> intf;
intf[ALL] = 1111*hila::random();
 
// Do the histogram of the last decimal digit in Field intf
Vector<10,int> dvec = 0;
onsites(ALL) {
    dvec[intf[X] % 10] += 1;
}
hila::out0 << "Histogram is " << dvec << '\n';

‍Using Vector<T> -types for histograms (as above) is inefficient for longer vectors. Use ReductionVector instead.

Reduction<T> -type

Using type Reduction<T>, where T can be any hila Field element type, enables delayed and non-blocking reductions which are not possible with the standard reductions. These are optimisations, the results of the reductions are not affected. This type implements only the sum reduction (+=).

The class Reduction<T> implements the following methods (see examples below):

• Option setting, called before reduction loop:

  • Reduction<T> & Reduction<T>::allreduce(bool on = true): set allreduce on/off (default: on)
  • Reduction<T> & Reduction<T>::delayed(bool on = true): set delayed on/off (default: off)
  • Reduction<T> & Reduction<T>::nonblocking(bool on = true): set nonblocking on/off (default: off)

  These return a reference to the reduction variable, so that these can be chained.

• Option values can be inquired with functions
  • bool Reduction<T>::is_allreduce()
  • bool Reduction<T>::is_delayed()
  • bool Reduction<T>::is_nonblocking()
• T Reduction<T>::value() : get the final value of the reduction.
• void Reduction<T>::reduce() : complete the reduction for non-blocking and/or delayed reduction. For other reductions does nothing. Is implicitly called by value() if not called before.
• void Reduction<T>::start_reduce() : start the non-blocking reduction (if non-blocking is on). Must be called by all MPI ranks.
• Assignment operator = : initializes or resets the value.
• Compound add operator += : reduction operator inside onsites().

The first call of value() or reduce() must be called by all ranks, otherwise the program deadlocks. Subsequent calls to value() just return the previously reduced value.

Delayed reduction means that the MPI reduction is delayed until the final result is desired:

Field<Complex<double>> f[NDIM];
...
 
Reduction<Complex<double>> r;  // definition of reduction variable
r = 0; 
r.allreduce(false).delayed();  // set delay, and allreduce = off
foralldir(d) {
    onsites(ALL) {
        r += exp(f[d][X]);
    }  
    // MPI reduction is not yet done here, as with the std. reduction
}
// r.value() gives the result. Could also call r.reduce() before r.value().
hila::out0 << "Sum is " << r.value() << '\n';

Reduction across MPI ranks is done in call r.value(). Thus, there is only one (MPI) reduction operation, whereas the standard method (using Complex<double> r) would perform MPI reduction after each onsites(). A The function r.value() must be called by all ranks or program deadlocks. The result of the reduction is not affected.

Non-blocking reduction enables overlapping MPI communication and computation (implemented using MPI_Ireduce()). This potentially gives higher performance.

Reduction<Complex<double>> r(0);
r.nonblocking();
onsites(ALL) {
    r += f[X] + g[X];
}
// reduction starts automatically after onsites()
 
// do some computations here not involving r
...
 
// reduction is completed when value()-method is called
hila::out0 << "Reduction is " << r.value() << '\n';

Non-blocking and delayed reduction can be combined. In this case, the reduction should be started with start_reduce()-method:

Reduction<Complex<double>> r(0);
r.allreduce(false).nonblocking().delayed();
foralldir(d) {
    onsites(ALL) {
        r += f[d][X];
    }  
}
// start reduction across ranks
r.start_reduce();
 
// do something else
... 
 
// reduction is completed with value()-method. Could also call .reduce() 
hila::out0 << "Result is " << r.value() << '\n';

If r.start_reduce() is omitted, it is implicitly done at r.value(), ensuring correct results. However, the possible benefit of overlapping communication and computation is lost.

ReductionVector<T>

ReductionVector enables (weighted) histogramming or general reduction to arrays. As an example, measuring "wall-to-wall" correlation function

Field<Complex<double>> F;
...
// allocate ReductionVector of length lattice.size(e_z)
ReductionVector<Complex<double>> wall(lattice.size(e_z));
wall.allreduce(false);
wall = 0;
 
// average values of F on each z-plane
onsites(ALL) {
    wall[X.coordinate(e_z)] += F[X];
}
 
// reduction done, calculate periodic correlation function
if (hila::myrank() == 0) {
    for (int dist = 0; dist <= lattice.size(e_z)/2; dist++) {
        Complex<double> corr = 0;
        for (int z = 0; z < lattice.size(e_z); z++) {
            corr += wall[z].dagger() * wall[(z + dist) % lattice.size(e_z)];
        }
        hila::out0 << dist << ' ' << corr << '\n';
    }
}

The type implements the methods (for clarity, using placeholder variable name a)

• Constructor ReductionVector<T> a(int size) : initializes and allocates vector of size elements
• Option setting, called before reduction loop:
  • a.allreduce(bool on = true): set allreduce on/off (default: on)
  • a.delayed(bool on = true): set delayed on/off (default: off)
• Boolean query methods a.is_allreduce(), a.is_delayed().
• a[i] where i is int : access the element number i.
• a.reduce() : in delayed reduction completes the reduction. This must be called before the the result is used.
• Assignment a = <value> : fills in a with the value
• std::vector<T> a.vector() : returns std::vector of the data
• Methods directly analogous to std::vector<T>:
  • size_t a.size() : size of a
  • a.clear() : zero all elements (equivalent to a = 0)
  • a.resize(new_size) : resize vector a
  • a.push_back(T val) : insert value val, increasing size of a
  • a.pop_back() : remove last element, decreasing size of a
  • T & front() : first element
  • T & back() : last element

Example of delayed ReductionVector use: make a histogram with 100 bins of N fields G[N], where we know the field values are between 0 and fmax:

Field<double> G[N];
double fmax = ...
...
 
ReductionVector<long> rv(100);     // 100 bins
rv = 0;
rv.allreduce(false).delayed(true);
for (int i = 0; i < N; i++) {
    onsites(ALL) {
        rv[ G[i][X] / fmax * 100] += 1;
    }
}
// do the reduction - without .reduce() result is undefined!
rv.reduce();
for (int i = 0; i < rv.size(); i++) {
    hila::out0 << i << ' ' << rv[i] << '\n';
}

‍NOTE: ReductionVector indexing has no protection. Using index past the ends of the array will cause trouble!

Shift and reflect operators

Field::shift operations allow shifting all elements by a certain displacement vector v. Even and Odd elements cannot be shifted with Field::shift method

CoordinateVector v = {0,1,1,0};
 
f = g.shift(v);                // these three
g.shift(v,f);                  //
f[ALL] = g[X + v];             // are equivalent
 
f[EVEN] = g[X + v];            // Cannot be done with g.shift() alone

Field::reflect operations reflect the field along a certain axis. Let g be an existing field of MyType, then we have:

Field<MyType> f;
f = g.reflect() // reflects the Filed g along all axes
f = g.reflect(e_x) // reflect the Field g along the x-axis
f = g.reflect({1,0,0,0}) // also reflect the Field g along the x-axis

with this in mind g.reflect() == g.reflect(e_x + e_y + e_z + e_t) == g.reflect({1,1,1,1}).

Random number generators

Random number generators are used with onsites loops. There are two random number generators available.

Uniform distribution:

Field<MyType> f;
onsites(ALL) f[X].random();

Gaussian distribution:

Field<MyType> f;
onsites(ALL) f[X].gaussian_random();

The used gaussian_random or random function is defined based on the element type MyType.

FFT

HILA has built-in Fast Fourier Transform for Field variables, using method Field::FFT. Complex-to-complex, complex-to-real and real-to-complex transformations are implemented.

Complex-to-complex FFT

The complex-to-complex FFT is defined for Field<T> -variables where T is a type which contains complex numbers, for example Complex<float>, Vector<4,Complex<double>>, SU<3,float> (SU<n> matrix is built with complex numbers). The result of the FFT is a Field<T> of the type of the input.

Forward and inverse FFTs are defined as

‍ \( \tilde f[\mathbf{n}] = \sum_{\mathbf{x}} f[\mathbf{x}] \exp(i \mathbf{k}\cdot\mathbf{x}), ~~~~~ f'[\mathbf{x}] = \sum_{\mathbf{n}} \tilde f[\mathbf{n}] \exp(-i \mathbf{k}\cdot\mathbf{x}) = \mbox{Volume} \times f[\mathbf{x}] \)

where \( n_i = 0 \ldots (N_i-1)\), \( ~ N_i\) is the lattice size to direction \(i\) and \( k_i \equiv 2\pi n_i/N_i\).

‍NOTE: as often in numerical FFTs, FFT followed by inverse FFT gives the original field multiplied by the lattice volume.

Example:

Field<Vector<3,Complex<double>>> f;
f = ...
auto res = f.FFT()         // Forward transform
auto res_2 = res.FFT(fft_direction::back) / lattice.volume();  // res_2 is same as f
f = f.FFT();               

Above the argument fft_direction::back indicates to perform the inverse transform \(\mathcal{F}^{-1}\). By default the arugment fftdir is fft_direction::forward.

HILA provides a convenience function Vector<NDIM, double > CoordinateVector::convert_to_k() which converts the momentum lattice coordinate n to wave vector k so that \(~ k_i = 2\pi n_i / N_i ~\) if \(~ 0 \le n_i \le N_i/2 ~\), and \(~ k_i = 2\pi (n_i - N_i)/N_i ~\) if \(~ n_i > N_i / 2~ \). Now \( ~ \frac \pi 2 < k_i \le \frac \pi 2 ~ \).

For example, if we want to generate a field where the wave numbers are random but with magnitude \(~ \exp(-k^2/\sigma^2)~\):

Field<Complex<double>> f;
onsites(ALL) {
    // sigma2 is the variance in k-space
    double width = exp(-X.coordinates().convert_to_k().squarenorm() / sigma2 );
    f[X].gaussian_random(width);
}
f = f.fft(fft_direction::back);  // convert to x-space

One can also specify the coordinate directions along which FFT is done with a CoordinateVector or Direction. FFT is done only to directions where the CoordinateVector is non-zero:

auto res = f.FFT(e_x);                                              //Forward transform only to x-direction
auto res_2 = res.FFT(e_x,fft_direction::back) / lattice.size(e_x);  // res_2 is same as f
 
CoordinateVector d;
d.fill(1);
d[e_y] = 0;
res = f.FFT(d);     // FFT is done to other directions except to y

With this in mind f.FFT(e_x + e_y + e_z) == f.FFT() if the lattice is 3-dimensional.

Real-to-complex / complex-to-real FFT

Real-to-complex FFT can only be done for Fields of type Field<float> or Field<double>. As opposed to complex-to-complex FFT the element types cannot be composites.

Real-to-complex FFT is done with function Field<T>::FFT_real_to_complex([coordinate directions],[FFT direction]), where the coordinate directions and FFT direction are optional. The result is a Field<Complex<float>> or Field<Complex<double>>, depending on the input Field.

Field<float> g;
g = ...
auto gt = g.FFT_real_to_complex();   // gt is type Field<Complex<float>>

The result obeys gt[-n] == gt[n].conj(), where negative n is understood to be shifted to the opposite side of the lattice \( -n_i = N_i - n_i\) (i.e mod to interval \(N_i\)).

In complex-to-real FFT only half of the degrees of freedom of the complex Field are used. The routine implicitly symmetrizes the field (f[-n] = f[n].conj()) so that the user does not have to do it.

Field<float> g;
g = ...
 
auto gt = g.FFT_real_to_complex();       // gt is type Field<Complex<float>>
 
// transform back, now g_new = g * lattice.volume()
auto g_new = gt.FFT_complex_to_real(fft_direction::back);

Helper routine hila::FFT_complex_to_real_site(CoordinateVector &cv) can be used to find out which sites are significant in complex-to-real transform:

hila::FFT_complex_to_real_site(cv) =  1 :  both real and imaginary parts significant at site cv
                                   =  0 :  only real part is significant
                                   = -1 :  value ignored on this site

As an example, sites on 2 dimensions on a 8x8 lattice: (+ and - refer to values +1 and -1, and * is origin (0,0), which has value 0)

- + + + - - - -                          - - - 0 + + + 0
- + + + - - - -                          - - - + + + + +
- + + + - - - -    after centering       - - - + + + + +
0 + + + 0 - - -    (0,0) to center       - - - + + + + +
+ + + + + - - -    ----------------->    - - - * + + + 0
+ + + + + - - -                          - - - - + + + -
+ + + + + - - -                          - - - - + + + -
  + + + 0 - - -                          - - - - + + + -

The value on sites with - are reconstructed from + -sites using f[-n] = f[n].conj().

Mathematical methods

Standard arithmetic methods:

• Field::operator+
• Field::operator-
• Field::operator*
• Field::operator/

The standard arithmetic methods above operate element wise. Both fields must have compatible types. All operators are also defined for scalars. Let @ be one of the four operators above, then the operator works in the following ways:

Field<MyType> f,g;
MyType a;
.
.
.
f = f @ g; // Field operated with other field
f = f @ a; // Field operated with scalar 

Arithmetic assignment methods:

• Field::operator+=
• Field::operator-=
• Field::operator*=
• Field::operator/=

The same which holds for the arithmetic operators hold also for the assignment operators.

Field<MyType> f,g;
MyType a;
.
.
.
f @= g; // Field operated with other field
f @= a; // Field operated with scalar 

Mathematical functions:

Element wise functions:

• Field::exp
• Field::log
• Field::sin
• Field::cos
• Field::tan
• Field::asin
• Field::acos
• Field::atan
• Field::abs
• Field::pow
• Field::conj
• Field::dagger
• Field::real
• Field::imag

Let func be a mathematical method defined above, then usage is:

Field<MyType> f;
f = func(f)

An exception to this is pow which takes in additional arguments, see Field::pow documentation for details. Additionally norm, conj, dagger, real and imag can be used as member functions f.fun().

All mathematical functions are operated element wise on the Field using the definition of the element type T of the Field<T>:

Field<R> res;
onsites(ALL) {
    res[X] = func(arg[X]);
}
return res;

if func is not defined for type T then an error is thrown.

Scalar returning functions

• Field::norm
• Field::squarenorm
• Field::squarenorm_relative

All of these mathematical methods return a scalar metric of the Field. Let func be a mathematical method defined above, then usage is:

Field<MyType> f;
double a;
a = func(f)

An exception to this is pow which takes in additional arguments, see Field::squarenorm_relative documentation for details.