The Taming of QCD by Fortran 90

TL;DR

This paper presents a comprehensive implementation of lattice QCD in Fortran 90, featuring modular, machine-independent code, efficient data storage, and parallel computing capabilities, enabling effective simulations on supercomputers.

Contribution

The authors developed a modular, machine-independent Fortran 90 framework for lattice QCD with optimized data handling and parallel processing, enhancing simulation efficiency and portability.

Findings

Successful implementation on supercomputers

High performance in parallel and single processor modes

Ease of writing QCD programs with new modules

Abstract

We implement lattice QCD using the Fortran 90 language. We have designed machine independent modules that define fields (gauge, fermions, scalars, etc...) and have defined overloaded operators for all possible operations between fields, matrices and numbers. With these modules it is very simple to write QCD programs. We have also created a useful compression standard for storing the lattice configurations, a parallel implementation of the random generators, an assignment that does not require temporaries, and a machine independent precision definition. We have tested our program on parallel and single processor supercomputers obtaining excellent performances.

Equations

$$f_{2,\mathbf{x}}=(1-{\gamma }_{\mu })U_{\mathbf{x}}^{\mu }{f}_{1,\mathbf{x}+\hat{\mu }}$$

$$f_{2,\mathbf{x}}=(1+{\gamma }_{\mu })U_{\mathbf{x}-\hat{\mu }}^{†\mu }{f}_{1,\mathbf{x}-\hat{\mu }}$$

$$s_{i+1}=a_1\ast s_i+b_1\mathrm{mod}{2}^{48},$$

Full text

1 Introduction

With the twofold goal of facilitating the development of algorithms and applications for lattice QCD, and of maintaining good code performance, we have taken advantage of the possibilities offered by Fortran 90 to write a set of modules for a high-level, yet efficient implementation of QCD simulations. In particular Fortran 90 offers the possibility to define both types and overloaded operators. These two key features make Fortran 90 particularly suitable for QCD simulations.

Our end product is a package, “QCDF90”, which is fully described in a long documentation [1], where we provide all the information needed to use our package. In the following we will highlight the main characteristics of QCDF90.

2 Geometry and field definitions

The set of all lattice sites can be subdivided into “even” and “odd” sites according to whether the sum of the integer valued coordinates $\tt x+y+z+t$ is even or odd (checkerboard subdivision). There are many algorithms which demand, especially in the context of a parallel implementation, that even and odd variables be treated separately. In QCDF90 we have implemented such checkerboarded separation of the lattice in two sublattices. Correspondingly all field variables are divided into even and odd variables. The first component of the type definition is an integer variable $\tt parity$ which will take values $\tt 0$ and $\tt 1$ for variables defined over even and odd sites respectively.

For the gauge variables it is convenient to include in the type a single direction $\mu$ component (of definite parity, of course). The second component of the gauge field type is an integer variable $\tt dir$ which takes values from $\tt 1$ to $\tt 4$ for variables defined in the corresponding direction.

In a vectorized or superscalar architecture pipelined instructions and longer arrays give origin to better performance. Therefore the lattice is most efficiently indexed by a single lattice volume index ranging $\tt xyzt$ from [math] to $NX*NY*NZ*NT/2-1$ for each sublattices (where $Ni$ is the lattice size in direction $i$).

We have defined the following field types: $\tt gauge\_field$ (a $3*3$ complex matrix in a given direction); $\tt fermi\_field$ (a $3$ component complex vector times $4$ spinor indices); $\tt complex\_field$ (a complex scalar); $\tt real\_field$ (a real scalar); $\tt generator\_field$ ($8$ real variables in a given direction, for the $SU(3)$ generators). These are defined on an even or odd sublattice. In addition we define the type $\tt full\_gauge\_field$ (a collection of $8$ gauge_field) and the type $\tt matrix$ as a $3*3$ complex matrix.

3 Overloaded operators

The above type definitions are easily manipulated once one defines overloaded operators for all possible operations between fields, matrices, complex and real numbers. The overloaded operator set includes multiplication, multiplication with the adjoint, division, addition, subtraction, lambda matrices algebraic manipulations, gamma matrices algebraic manipulations, adjoin-ing, conjugation, real and complex traces, exponentiation, square root, contraction, etc…

For example, if $\tt g_{i}$ are gauge fields, the use of overloaded operators allows instructions to be as simple as $\tt g_{1}=g_{2}*g_{3}$.

To implement some useful algebraic operations on the $SU(3)$ generators we have overloaded some further operators which perform special operations involving $\tt generator\_fields$ and $\tt gauge\_fields$. In particular it is very important to have an efficient algorithm for the exponentiation of a matrix, since this operation can be a time consuming component of several QCD calculations. The algorithm that we have used takes advantage of the properties of the $3*3$ Hermitian traceless matrices to perform the exponentiation with a minimal number of arithmetic operations.

4 Shifts

Shifts are also implemented as overloaded operators. For each field types C-shift implements an ordinary shift of the field with respect to the Cartesian geometry of the lattice.

Moreover, because gauge theories are characterized by the property of local gauge invariance, we find it very useful to directly define a U-shift operator that consists of the shift with the appropriate parallel transport factor. In a gauge theory the U-shift is the relevant shift operator and is the natural building block for programming.

In QCD the efficient manipulation of the Dirac operator is a very critical issue. Therefore for the Fermi fields we have also defined other shift operators which incorporate fundamental features of the Dirac operator. First we define a W-shift that combines the shift with the parallel transport and the appropriate gamma matrix manipulation, i.e. acting on a Fermi field $\tt f_{1}$ W-shift produces a Fermi field $\tt f_{2}$, given by

[TABLE]

for positive W-shift, and

[TABLE]

for negative W-shift. The direct definition of this combined operator entails advantages of efficiency because, from the properties of the gamma matrices, it follows that only one half of the spin components undergo the transport.

We overload the (Wilson) lattice Dirac operator and, finally, we overload the operator Xdirac$=\gamma_{5}$ Dirac $\gamma_{5}$.

5 Assignments

The use of overloaded operators may imply the creation of more temporaries and, consequently, more motion of data than a straightforward implementation of operations among arrays. Consider for example the following operation among variables of type fermi_field: $\tt f1=f1+f2+f3$. As far as we know, Fortran 90 does not specify how the variables should be passed in function calls. As a consequence, the above instruction requires as many as four temporaries. The procedure could be drastically simplified through the use of an overloaded assignment $\tt+=$. The above instruction could be written $\tt f1~{}+=f2+f3$ which the compiler would implement by issuing first a call to a function that adds $\tt f2$ and $\tt f3$ returning the result in $\tt t1$. The addresses of $\tt f1$ and $\tt t1$ would then be passed to a subroutine that implements the operation $\tt f1=f1+t1$ among the components of the data types. The required number of copies to memory would be only two, instead of four.

In order to allow for these possible gains in efficiency, we have defined a large set of overloaded assignments. Since Fortran 90 permits only the use of the $\tt=$ symbol for the assignment, we defined two global variables: a character variable $\tt assign\_type$ and an integer variable $\tt assign\_spec$ (for assign specification, introduced to accommodate assignments of a more elaborate nature). The default values of these variables are “$\tt=$” and “$\tt 0$”. Overloaded assignments are obtained by setting $\tt assign\_type$ (and possibly $\tt assign\_spec$) to the appropriate value immediately before the assignment. (Our example become $\tt assign\_type=^{\prime}+^{\prime};f1=f2+f3$).

In the module “assign_mixed”, assignments are also defined between variables of different types. For example if $\tt c_{1}$ is a complex_field and $\tt complex$ is a complex variable, the instruction $\tt complex=c_{1}$, is interpreted as setting the variable $\tt complex$ to the sum over all the lattice of the components of $\tt c_{1}$.

6 Random numbers

We have implemented a parallelizable version of the unix pseudorandom number generator erand48, which also provides added functionality. Erand48 is a congruential pseudorandom number generator based on the iterative formula

[TABLE]

where $a_{1}=\tt 0x5DEECE66D$, $b_{1}=\tt 0xB$, $s_{i}$ and $s_{i+1}$ are integers of at least 48 bits of precision. The “seeds” $s_{i}$ are converted to real pseudorandom numbers $r_{i}$ with uniform distribution between [math] and $1$ by $r_{i}=2^{-48}\,s_{i}$.

As presented above, the algorithm is intrinsically serial. However it follows from Eq. (3) that the $\rm N^{th}$ iterate $s_{i+N}$ is still of the form $s_{i+N}=a_{N}*s_{i}+b_{N}\;{\rm mod}\;2^{48}$ with integers $a_{N}$ and $b_{N}$ which are uniquely determined by $a_{1}$, $b_{1}$. The module takes advantage of this fact and of the definitions of a global variable $\tt seeds$ to generate pseudorandom numbers in a parallelizable fashion. The module “random_numbers” overloads operators that generate Gaussian or uniformly distributed fields and all the necessary seed manipulations.

7 Conditionals

The module “conditional” defines six overloaded relational operators, $\tt>$, $\tt>=$, $\tt<$, $\tt<=$, $\tt==$, $\tt/=$, and the $\tt.Xor.$ operator.

The relational operators return a dummy logical variable and at the same time the global variable $\tt context$ is set to $\tt.TRUE.$ at all sites where the relation is satisfied, to $\tt.FALSE.$ at all other sites.

The operator $\tt.Xor.$ admits as operands a pair of fields of the same type and returns a field, also of the same type, having as elements the corresponding elements of the first operand at the sites where the global variable $\tt context$ is $\tt.TRUE.$, the elements of the second operand at the sites where $\tt context$ is $\tt.FALSE.$. This can be used to select elements out of two fields according to some local condition, an operation which lies at the foundation of stochastic simulation techniques.

8 Precision

To render the precision definitions machine independent, the module “precision” defines two $\tt kind$ parameters, $\tt REAL8$ and $\tt LONG$. These parameters store the $\tt kind$ of an 8-byte floating point variable and of an 8-byte integer variable.

9 Write and read configurations

To store and retrieve an entire SU(3) gauge field configuration, we developed a portable, compressed ASCII format. Only the first two columns of the gauge field matrices are stored, thanks to unitarity and unimodularity. Our subroutines takes advantage of the fact that all of the elements of the gauge field matrices have magnitude smaller or equal to 1 to re-express their real and imaginary parts as 48bit integers. These integers are then written in base 64, with the digits being given by the ASCII collating sequence starting from character “0” (to avoid unwanted ASCII characters). Thus, an entire gauge field matrix is represented by 96 ASCII characters, without loss of numerical information. A detailed description of the compressing procedure is written, as a header, at the beginning of the stored configuration file itself.