Measure the trace of the product of two wilson lines going in opposite time directins separated by distance r (here 'r' is the physical sepration between two temporal transporters(wilson lines in time directions).).
/*
* main_wilson_lines_correlator_stacked.cu
*
* Rasmus Larsen, 25 Feb 2021
*
*/
#include "../simulateqcd.h" /* This is the principle header file which is stored in ../ directory relative to this source code. This header file is primary header file for the applications, profiling, and test source codes. This header file contains the inclusion of many header files from base, gauge, spinnor folders.For more look at [[Header_files]]. The header files basically contains veriable declaration, function declaration, or function defination. Here #include is preprocessor which tells compiler to load the header file at given address or relative address */
#include "../modules/observables/wilsonLineCorrelatorMultiGPU.h"
#define PREC double //this is called macro, for more look at [[scalar_cpp_macro]], this is an object like macro, here we have defined the precision as double which is widely used in scientific calculations.
#define STACKS 96 // this is object like macro
template<class floatT> //this is an function template used for generic programming(i.e. for writing a single function for different data types, here then compiler creates code accordingly as one calls these functions based on different data types(which is also known as instantiation)), for more look at [[Generics in c++ (Template Introduction)]]
struct WLParam : LatticeParameters { /* This is structures in c++, here we have defined a structure named WLParam which is inherited from another structure (LatticeParameters), for more look at [[https://www.programiz.com/cpp-programming/structure|c++ structures_programizz]] */
Parameter<floatT> gtolerance; // these are member template variables of the structure WLParam, with typename Parameter which is defined elsewhere which can be found in simulateqcd.h header file. The definition of this Parameter class can be found at SIMULATeQCD/src/base/IO/parameterManagement.h.
Parameter<int,1> maxgfsteps; /* In short all these variables of structure are named as structure members */
Parameter<int,1> numunit;
Parameter<int> load_conf;
Parameter <std::string> gauge_file;
Parameter <std::string> directory;
Parameter <std::string> file_type;
WLParam() { /* This is Constructor of the structure named WLParam also known as Member function, here the default values are provided*/
addDefault (gtolerance,"gtolerance",1e-6);
addDefault (maxgfsteps,"maxgfsteps",9000);
addDefault (numunit ,"numunit" ,20);
addDefault(load_conf, "load_conf", 0);
addOptional(gauge_file, "gauge_file");
add(directory, "directory");
add(file_type, "file_type");
}
};
template<class floatT> /* Here is again some structure inherited from class ParameterList, just below are structure members and member functions/constructor (milcInfo) defined. */
struct milcInfo : ParameterList {
Parameter<floatT> ssplaq;
Parameter<floatT> stplaq;
Parameter<floatT> linktr;
milcInfo() {
add(ssplaq, "gauge.ssplaq");
add(stplaq, "gauge.stplaq");
add(linktr, "gauge.nersc_linktr");
}
};
/* Below is the start of main() function which will contain the most part of the program. Here it main takes two arguments of type integer and char which are variable(named argc) and character pointer named argv[](), In C and C++ programming, the `int main(int argc, char *argv[])` is the signature of the `main` function, which serves as the entry point for the program. Let's break down the components of this signature:
- **`int main`:** Indicates that the `main` function returns an integer value. The returned integer is often used as an exit status, where 0 typically signifies successful execution, and nonzero values indicate errors.
- **`int argc`:** Stands for "argument count" and represents the number of command-line arguments passed to the program, including the name of the program itself.
- **`char *argv[]`:** Stands for "argument vector" and is an array of strings (`char *`). It holds the actual command-line arguments. `argv[0]` is the name of the program, and subsequent elements (`argv[1]`, `argv[2]`, etc.) contain the additional arguments.
The `argv` parameter in the `main` function is declared as a pointer to an array of strings. The declaration `char *argv[]` indicates an array of pointers to characters, and in the context of the `main` function, it is often used to represent an array of strings.
In C and C++, arrays are often represented as pointers to their first elements, and strings are represented as arrays of characters. So, `char *argv[]` indicates an array of pointers to strings.
Each element of the `argv` array is a pointer to the first character of a null-terminated string. The first element (`argv[0]`) is typically the name of the program itself, and subsequent elements (`argv[1]`, `argv[2]`, and so on) are the command-line arguments passed to the program.
It's important to note that `argv` is a pointer to an array, and you can use array indexing (`argv[i]`) or pointer arithmetic to access individual elements of the array. The `argc` parameter provides the count of elements in the `argv` array.
Here's a brief explanation of how this works:
- `argc` is the count of command-line arguments, so it tells you how many items are in the `argv` array.
- `argv` is an array of strings where each element is a command-line argument. The first element (`argv[0]`) is the name of the program itself.
A simple example of using command-line arguments might look like this:
```cpp
#include <iostream>
int main(int argc, char *argv[]) {
// Print the program name
std::cout << "Program name: " << argv[0] << std::endl;
// Print command-line arguments
for (int i = 1; i < argc; ++i) {
std::cout << "Argument " << i << ": " << argv[i] << std::endl;
}
return 0;
}
If you run this program with command-line arguments like `./program arg1 arg2`, it will output:
Program name: ./program
Argument 1: arg1
Argument 2: arg2
This mechanism allows programs to accept inputs from the command line, making them more flexible and configurable.*/
int main(int argc, char *argv[]) {
/// Controls whether DEBUG statements are shown as it runs; could also set to INFO, which is less verbose.
stdLogger.setVerbosity(INFO); /* This line suggests that there is a logger or logging system in use, and it's setting the verbosity level to INFO. Logging is a way to output information about the program's behavior to aid in debugging or monitoring. The verbosity level determines the amount of detail in the log messages, and INFO is a standard level indicating informational messages. For more refer to [[https://latticeqcd.github.io/SIMULATeQCD/02_contributions/terminalIO.html|SIMULATeQCD documentation]]*/
/// Initialize parameter class.
WLParam<PREC> param;/* It declares an instance of the WLParam template structure with the template parameter PREC */
/// Initialize the CommunicationBase.
CommunicationBase commBase(&argc, &argv);/* This line initializes an instance of the CommunicationBase class (i.e. creating an object named commBase), passing the addresses of argc and argv to its constructor. This suggests that CommunicationBase might be involved in handling communication-related tasks and could be designed to work with command-line arguments. */
param.readfile(commBase, "../parameter/test.param", argc, argv); /* This line calls a readfile method on the param object, passing several arguments including the commBase object, a file path ("../parameter/test.param"), and the command-line arguments argc and argv. This implies that the readfile method is responsible for reading configuration or parameter information from a file, and the commBase object might be used in this process. */
commBase.init(param.nodeDim()); /* This line initializes the commBase object by calling its init method and passing the result of param.nodeDim() as an argument. This implies that the nodeDim method of the param object(which is instantiation of WLParam<template>class) returns some value related to the dimensions of nodes, and this information is used to initialize the communication base. */
// cout << param.nodeDim[0] << " param 0 " << param.nodeDim[1] << " param 1 " << param.nodeDim[2] << " param 2 " << param.nodeDim[3] << " param 3 " <<endl;
// Set the HaloDepth.
const size_t HaloDepth = 2;
rootLogger.info("Initialize Lattice"); /* For more click here, [[Terminal output and terminating the program]] */
/// Initialize the Lattice class.
initIndexer(HaloDepth,param,commBase); /* this function is defined in SIMULATeQCD/src/base/indexer/initCPUIndexer.cpp ,
The line `initIndexer(HaloDepth, param, commBase);` appears to be a function call to a function named `initIndexer`. Let's break down the components:
1. **`initIndexer`**:
- This is the name of the function being called. The function likely performs some initialization related to an indexer. The exact behavior and purpose of this function depend on its definition elsewhere in your codebase.
2. **`(HaloDepth, param, commBase)`**:
- This is the list of arguments being passed to the `initIndexer` function. It appears to be passing three arguments:
- `HaloDepth`: This is likely a variable or constant that indicates the depth of a halo region. The purpose of this parameter depends on the context and the design of your program.
- `param`: This is an instance of the `WLParam` template class, which presumably holds parameters relevant to the initialization process.
- `commBase`: This is an instance of the `CommunicationBase` class, suggesting that communication-related functionality may be involved in the initialization.
The purpose of the `initIndexer` function would be defined in the codebase, possibly performing tasks such as setting up indices or configurations related to indexing, and using the provided parameters for the initialization process. Without the definition of the `initIndexer` function and additional context from the surrounding code, it's not possible to provide more specific details about its behavior. */
/// Initialize the Gaugefield.
rootLogger.info("Initialize Gaugefield");
Gaugefield<PREC,true,HaloDepth> gauge(commBase); /* concept of halo: [[Multi-GPU: Distribution of local lattices on the individual GPUs]] */
/// Initialize gaugefield with unit-matrices.
gauge.one();
std::string gauge_file;
// load gauge file, 0 start from 1, 1 and 2 load file, 2 will also gauge fix
/* i.e. If we don't have any configuration availalbe, it will use the unit matrices at each links. If we have configuration availalbe then it will load the configuration. */
if (param.load_conf() == 0)
{
rootLogger.info("Starting from unit configuration");
gauge.one();
}
else if(param.load_conf() == 2 || param.load_conf() == 1)
{
std::string file_path = param.directory();
file_path.append(param.gauge_file()); //adding full path to gague_file configuration.
rootLogger.info("Starting from configuration: " , file_path);
// rootLogger.info(param.gauge_file() , endl);
if(param.file_type() == "nersc"){
gauge.readconf_nersc(file_path);
}
else if(param.file_type() == "milc"){
gauge.readconf_milc(file_path);
//Exchange Halos
gauge.updateAll(); /* This line likely invokes a method updateAll on the object gauge. The comment //Exchange Halos suggests that this method may be involved in updating or exchanging information related to halo regions in the gauge field. Halo regions often represent additional boundary layers beyond the core computational domain and are commonly used in parallel computing to facilitate communication between neighboring processes */
GaugeAction<PREC,true,HaloDepth> enDensity(gauge); /* This line declares and initializes an object enDensity of the GaugeAction class template. It appears to be associated with actions or computations related to the gauge field. The template parameters include the precision (PREC), a boolean value (true), and a depth parameter for halo regions (HaloDepth). The gauge object is passed to the constructor, suggesting that the GaugeAction class is designed to operate on gauge fields.*/
PREC SpatialPlaq = enDensity.plaquetteSS();/* This line calculates the spatial plaquette (plaquette on spatial links) using the plaquetteSS method of the enDensity object. The result is stored in a variable SpatialPlaq. */
PREC TemporalPlaq = enDensity.plaquette()*2.0-SpatialPlaq;/* This line calculates the temporal plaquette (plaquette on temporal links) using the plaquette method of the enDensity object. It then subtracts twice the spatial plaquette (SpatialPlaq). The result is stored in a variable TemporalPlaq. */
rootLogger.info("plaquetteST: " , TemporalPlaq);/* This line logs an informational message using the info method of the rootLogger. The message includes the label "plaquetteST: " and the value of TemporalPlaq. */
rootLogger.info("plaquetteSS: " , SpatialPlaq);/* Similarly, this line logs another informational message including the label "plaquetteSS: " and the value of SpatialPlaq. */
std::string info_path = file_path;
info_path.append(".info");
milcInfo<PREC> paramMilc; /* This declares an instance of the milcInfo template class with the template parameter PREC and names the instance paramMilc. */
paramMilc.readfile(commBase,info_path);
rootLogger.info("plaquette SS info file: " , (paramMilc.ssplaq())/3.0);
rootLogger.info("plaquette ST info file: " , (paramMilc.stplaq())/3.0);
rootLogger.info("linktr info file: " , paramMilc.linktr());
if(abs((paramMilc.ssplaq())/3.0-SpatialPlaq) > 1e-5){
throw std::runtime_error(stdLogger.fatal("Error ssplaq!"));
}
if(abs((paramMilc.stplaq())/3.0-TemporalPlaq) > 1e-5){
throw std::runtime_error(stdLogger.fatal("Error stplaq!"));
}
}
}
// Exchange Halos
gauge.updateAll();
// Initialize ReductionBase.
LatticeContainer<true,PREC> redBase(commBase);
// We need to tell the Reductionbase how large our array will be. Again it runs on the spacelike volume only,
/// so make sure you adjust this parameter accordingly, so that you don't waste memory.
typedef GIndexer<All,HaloDepth> GInd; /* The line `typedef GIndexer<All,HaloDepth> GInd;` is creating a typedef in C++. Let's break down the components:
- `GIndexer<All, HaloDepth>`: This part is a template instantiation. It creates a type using the `GIndexer` template, where `All` and `HaloDepth` are template parameters.
- `typedef GIndexer<All, HaloDepth> GInd;`: This line declares a new type named `GInd` using the `typedef` keyword. It essentially creates an alias for the type `GIndexer<All, HaloDepth>`. After this line, you can use `GInd` as a shorthand for `GIndexer<All, HaloDepth>`.
In this specific case, it seems like `GIndexer` is a template class, and `GInd` is an alias for a specific instantiation of that template with `All` and `HaloDepth` as template arguments.
This typedef might be used to simplify the code and improve readability, especially if the `GIndexer<All, HaloDepth>` type is used frequently in the program. Instead of writing the full template instantiation each time, developers can use the shorter and more convenient `GInd`. */
redBase.adjustSize(GInd::getLatData().vol4)
rootLogger.info("volume size " , GInd::getLatData().globvol4);
///////////// gauge fixing
if(param.load_conf() ==2){
GaugeFixing<PREC,true,HaloDepth> GFixing(gauge);
int ngfstep=0;
PREC gftheta=1e10;
const PREC gtol = param.gtolerance(); //1e-6; /// When theta falls below this number, stop...
const int ngfstepMAX = param.maxgfsteps() ; //9000; /// ...or stop after a fixed number of steps; this way the program doesn't get stuck.
const int nunit= param.numunit(); //20; /// Re-unitarize every 20 steps.
while ( (ngfstep<ngfstepMAX) && (gftheta>gtol) ) {
/// Compute starting GF functional and update the lattice.
GFixing.gaugefixOR();
/// Due to the nature of the update, we have to re-unitarize every so often.
if ( (ngfstep%nunit) == 0 ) {
gauge.su3latunitarize();
}
/// Re-calculate theta to determine whether we are sufficiently fixed.
gftheta=GFixing.getTheta();
ngfstep+=1;
}
gauge.su3latunitarize(); /// One final re-unitarization.
rootLogger.info("Gauge fixing finished in " , ngfstep , " steps, with gftheta = " , gftheta);
}
//end of gauge fixing snippet
std::string Name = "WLine_";
std::string Name_r2 = "WLine_r2_";
if(param.load_conf() == 2 || param.load_conf() == 1){
Name.append(param.gauge_file());
Name_r2.append(param.gauge_file());
}
else{
Name.append("one");
Name_r2.append("one");
}
FileWriter file(gauge.getComm(), param, Name);
FileWriter file_r2(gauge.getComm(), param, Name_r2);
/* Initialization of Name and Name_r2:
1. Two std::string variables, Name and Name_r2, are initialized with base names ("WLine_" and "WLine_r2_" respectively).
Condition Based on param.load_conf():
2. Checks if the load_conf parameter in the param object is equal to 2 or 1.
If true, it appends the value of param.gauge_file() to both Name and Name_r2. This suggests that the file names are being constructed based on the value of gauge_file.
If false, it appends the string "one" to both Name and Name_r2.
FileWriter Initialization:
3. Two instances of the FileWriter class are initialized:
file: Initialized with the communication object (gauge.getComm()), the param object, and the constructed Name.
file_r2: Similar to file, but initialized with the constructed Name_r2. */
StopWatch<true> timer;
/// Start timer.
timer.start();
/// Exchange Halos
gauge.updateAll();
PREC dot;
WilsonLineCorrelatorMultiGPU<PREC,HaloDepth,STACKS> WilsonClass;
std::vector<PREC> dotVector;
PREC * results;
results = new PREC[(GInd::getLatData().globvol3/2+GInd::getLatData().globLX*GInd::getLatData().globLY)*GInd::getLatData().globLT];
///
timer.start();
//// loop over length of wilson lines
for(int length = 1; length<GInd::getLatData().globLT+1;length++){
/// calculate the wilson line starting from any spacetime point save in mu=0 direction
WilsonClass.gWilson(gauge, length);
/// copy from mu=0 to mu=1
gauge.template iterateOverBulkAtMu<1,256>(CopyFromMu<PREC,HaloDepth,All,0>(gauge));
gauge.updateAll();
// initial position x0=-1 due to adding dx in first line
int x0 = -STACKS;
int y0 = 0;
int z0 = 0;
int dx = STACKS;
int dy = 1;
int dz = 1;
int points = (GInd::getLatData().globvol3/2+GInd::getLatData().globLX*GInd::getLatData().globLY);
for(size_t i = 0; i<GInd::getLatData().globvol3/2+GInd::getLatData().globLX*GInd::getLatData().globLY;i+=STACKS){
x0 += dx;
if(x0 >= (int)GInd::getLatData().globLX || x0 <0){
dx *= -1;
x0 += dx;
y0 += dy;
if(y0 >= (int)GInd::getLatData().globLY|| y0 <0){
dy *= -1;
y0 += dy;
z0 += dz;
/// move mu=1 direction by dz
WilsonClass.gMoveOne(gauge,2,dz);
gauge.updateAll();
}
else if(param.nodeDim[1]>1){
/// move mu=1 direction by dy
WilsonClass.gMoveOne(gauge,1,dy);
gauge.updateAll();
}
}
// A(x).A(x+r)^dag along direction x and also y if not split on different GPU's
if(param.nodeDim[1] == 1){
dotVector = WilsonClass.gDotAlongXYStackedShared(gauge,y0,redBase);
}
else{
dotVector = WilsonClass.gDotAlongXYStackedShared(gauge,0,redBase);
}
if(dx>0){
for(int j = 0;j < STACKS ; j++){
results[i+j+points*(length-1)] = dotVector[j];
}
}
else{
for(int j = 0;j < STACKS ; j++){
results[i+STACKS-1-j+points*(length-1)] = dotVector[j];
}
}
for(int j = 0;j < STACKS ; j++){
// rootLogger.info(x0+j , " " , y0 , " ", z0 , " " , length , " " , dotVector[j]);
file << x0+j << " " << y0 << " "<< z0 << " " << length << " " << dotVector[j] << "\n";
}
}
}
/////////////// Save into radius squared file
int r2max = 0;
r2max += GInd::getLatData().globLX*GInd::getLatData().globLX;
r2max += GInd::getLatData().globLY*GInd::getLatData().globLY;
r2max += GInd::getLatData().globLZ*GInd::getLatData().globLZ;
PREC * results_r2 = new PREC[r2max+1];
PREC * norm_r2 = new PREC[r2max+1];
for(int length = 1; length<GInd::getLatData().globLT+1;length++){
for(int ir2=0; ir2<r2max+1; ir2++) {
results_r2[ir2] = 0.0;
norm_r2[ir2] = 0.0;
}
int x0 = -1;
int y0 = 0;
int z0 = 0;
int dx = 1;
int dy = 1;
int dz = 1;
size_t entries = GInd::getLatData().globvol3/2+GInd::getLatData().globLX*GInd::getLatData().globLY;
for(size_t i = 0; i<entries;i++){
x0 += dx;
if(x0 >= (int)GInd::getLatData().globLX || x0 <0){
dx *= -1;
x0 += dx;
y0 += dy;
if(y0 >= (int)GInd::getLatData().globLY|| y0 <0){
dy *= -1;
y0 += dy;
z0 += dz;
}
}
dot = results[i+entries*(length-1)];
// save results
int ir2 = 0;
if(x0 > (int)GInd::getLatData().globLX/2){
ir2 += (x0-(int)GInd::getLatData().globLX)*(x0-(int)GInd::getLatData().globLX);
}
else{
ir2 += x0*x0;
}
if(y0 > (int)GInd::getLatData().globLY/2){
ir2 += (y0-(int)GInd::getLatData().globLY)*(y0-(int)GInd::getLatData().globLY);
}
else{
ir2 += y0*y0;
}
if(z0 > (int)GInd::getLatData().globLZ/2){
ir2 += (z0-(int)GInd::getLatData().globLZ)*(z0-(int)GInd::getLatData().globLZ);
}
else{
ir2 += z0*z0;
}
// factor for counting contributions
// Initial factor 2 for symmetry between z and -z
// double the factor if x or y = l/2 due to periodicity
double factor = 1.0;
if(z0 == (int)GInd::getLatData().globLZ/2 || z0 == 0){
factor = 0.5*factor;
if((y0 == (int)GInd::getLatData().globLY/2 || y0 == 0) && (x0 == (int)GInd::getLatData().globLX/2 || x0==0) ){
factor = 2.0*factor;
}
}
if( (z0 == (int)GInd::getLatData().globLZ/2 || z0==0) && (y0 == (int)GInd::getLatData().globLY/2 || y0 == 0) && (x0 == (int)GInd::getLatData().globLX/2 || x0==0) ){
factor = 0.5*factor;
}
results_r2[ir2] += factor*dot;
norm_r2[ir2] += factor;
}
/////// write to file
for(int ir2=0; ir2<r2max+1; ir2++) {
if(norm_r2[ir2] > 0.1){
file_r2 << length << " " << ir2 << " " << results_r2[ir2]/norm_r2[ir2] << "\n";
}
}
}
/////////////// final check
PREC sum = 0.0;
for(size_t i = 0; i<(GInd::getLatData().globvol3/2+GInd::getLatData().globLX*GInd::getLatData().globLY)*GInd::getLatData().globLT;i+=1){
sum += abs(results[i]);
}
timer.stop();
rootLogger.info("Time for operators: " , timer);
rootLogger.info("abs(sum) = " , sum);
delete [] results;
delete [] results_r2;
delete [] norm_r2;
return 0;
}