Remarks: Please avoid this code, it was written to enhance the Statistics but it is highly highly higly compute intensive.
/*
* main_plaquette.cpp
*
* Lukas Mazur, 10 Apr 2018 (original Author for Plaquette)
* Vishal Rao , 06 Jan 2024 (generalisation from plaquette to wilson loop)
*
* This is just an example how a very very basic program works. Look at src/testing/main_GeneralOperatorTest.cpp to see
* how to write more advanced GPU code.
*
*/
#include "../simulateqcd.h"
/* #include "../modules/hyp/hypParameters.h" was already included in ../modules/hypSmearing.h */
/* #include "../modules/hyp/smearParameters.h" */
#include "../modules/hyp/hypSmearing.h"
/* #include "../modules/hyp/hypSmearing.cpp" */
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <strings.h>
#include <vector>
#include <cmath>
#define PREC double
template<class floatT,size_t HaloDepth>
struct CalcWilson{
//Gauge accessor to access the gauge field
SU3Accessor<floatT> SU3Accessor;
int _length=1;
int _height=1;
//Constructor to initialize all necessary members.
CalcWilson(Gaugefield<floatT,true,HaloDepth> &gauge, int &length, int &height) : SU3Accessor(gauge.getAccessor()), _length(length),_height(height){
}
__device__ __host__ SU3<floatT> HalfLoop(int &length, int &height,gSite site,int &mu, int &nu)
{
// we have not taken site variable as reference as if it were take as reference, it will update the values at same memory location
SU3<floatT> temp;
typedef GIndexer<All, HaloDepth> GInd;
temp = SU3Accessor.getLink(GInd::getSiteMu(site, mu));
for (int len=1;len<length;len++){
site=GInd::site_up(site, mu);
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, mu));
}
site=GInd::site_up(site,mu);
for (int hei=1;hei<=height;hei++){
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, nu));
site=GInd::site_up(site,nu);
}
return temp;
}
//This is the operator that is called inside the Kernel
__device__ __host__ floatT operator()(gSite site) {
/// We need to choose the type of indexer. The first template is the layout of the lattice.
typedef GIndexer<All, HaloDepth> GInd;
int l=_length;
int h=_height;
/// Define a SU(3) matrix
SU3<floatT> temp1,temp2;
floatT result = 0;
for (int nu = 1; nu < 3; nu++)
{
for (int mu = 0; mu < nu; mu++)
{
temp1=HalfLoop(l,h,site,mu,nu);
temp2=HalfLoop(h,l,site,nu,mu);
result+=tr_d(temp1*dagger(temp2));
/* temp1=HalfLoop(h,l,site,mu,nu); */
/* temp2=HalfLoop(l,h,site,nu,mu); */
/* result+=tr_d(temp1*dagger(temp2)); */
//Return the result
//The return value will be stored in the array of the reductionbase at index site.isite.
}
}
return result;
}
};
template<class floatT,size_t HaloDepth>
struct CalcWilsonBresenham{
int _z;
int _columns;
int _rows;
int *_d_path_xy;
SU3Accessor<floatT> SU3Accessor;
CalcWilsonBresenham(Gaugefield<floatT,true,HaloDepth> &gauge, int *d_path_xy,int &rows,int &columns, int &z) : SU3Accessor(gauge.getAccessor()), _d_path_xy(d_path_xy),_rows(rows),_columns(columns),_z(z){}
//This is the operator that is called inside the Kernel
__device__ __host__ SU3<floatT> HalfLoopBresenhamLower(int *d_path_xy,int &rows,int &columns, int &z, gSite site, int &mu, int &nu,int &index)
{
/* int n_diagonal_coordinates=0; //could be useful for later improvements in code */
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
typedef GIndexer<All, HaloDepth> GInd;
for (int iy=0;iy<columns-1;iy++) //accessing all coordinates, one minus for last link in z direction
{
if(d_path_xy[iy+index*columns]==0)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,mu));
site=GInd::site_up(site,mu);
}
else if (d_path_xy[iy+index*columns]==1)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,nu));
site=GInd::site_up(site,nu);
}
}
for (int iz=0;iz<z;iz++) //iz is iterator in z direction
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //3-mu-nu is the fictitious time direction
site=GInd::site_up(site,3-mu-nu); //half loop will complete at end of loop.
}
return temp;
}
__device__ __host__ SU3<floatT> HalfLoopBresenhamUpper(int *d_path_xy,int &rows,int &columns, int &z, gSite site, int &mu, int &nu,int &index)
{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1);
typedef GIndexer<All, HaloDepth> GInd;
for (int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //3-mu-nu is the fictitious time direction
site=GInd::site_up(site,3-mu-nu); //half loop will complete at end of loop.
}
for (int iy=0;iy<columns-1;iy++) //accessing all coordinates, one minus for last link in z direction
{
if(d_path_xy[iy+index*columns]==0)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,mu));
site=GInd::site_up(site,mu);
}
else if (d_path_xy[iy+index*columns]==1)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,nu));
site=GInd::site_up(site,nu);
}
}
return temp;
}
__device__ __host__ floatT operator()(gSite site) {
/// We need to choose the type of indexer. The first template is the layout of the lattice.
typedef GIndexer<All, HaloDepth> GInd;
/// Define a SU(3) matrix
SU3<floatT> temp1,temp2;
floatT result = 0;
for (int nu = 1; nu < 3; nu++)
{
for (int mu = 0; mu < nu; mu++)
{
for(int i=0;i<_rows;i++)
{
temp1=HalfLoopBresenhamLower(_d_path_xy,_rows,_columns,_z,site,mu,nu,i); //passing 'i' to construct ith row in the function
for(int j=0;j<_rows;j++)
{
temp2=HalfLoopBresenhamUpper(_d_path_xy,_rows,_columns,_z,site,mu,nu,j);
result+=tr_d(temp1*dagger(temp2));
}
}
}
}
return result/(_rows*_rows);
}
};
//Function to compute the wilson loop using the above struct CalcWilson.
template<class floatT, size_t HaloDepth>
floatT WilsonLoop(Gaugefield<floatT,true, HaloDepth> &gauge, LatticeContainer<true,floatT> &redBase, int &l, int &h){
typedef GIndexer<All,HaloDepth> GInd;
const size_t elems = GInd::getLatData().vol4;
//Make sure, redBase is large enough
redBase.adjustSize(elems);
// we will iterate the process of finding wilson loop over whole bulk.
redBase.template iterateOverBulk<All, HaloDepth>(CalcWilson<floatT, HaloDepth>(gauge, l,h));
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
Wloop /= (GInd::getLatData().globalLattice().mult()*9); //3 loops(one averaged loop per plane) times 3 colors.
return Wloop;
}
//defining function for bresenham
template<class floatT, size_t HaloDepth>
floatT WilsonLoopBresenham(Gaugefield<floatT,true, HaloDepth> &gauge, LatticeContainer<true,floatT> &redBase, int *d_path, int &rows,int &columns, int &z){
typedef GIndexer<All,HaloDepth> GInd;
const size_t elems = GInd::getLatData().vol4;
//Make sure, redBase is large enough
redBase.adjustSize(elems);
// we will iterate the process of finding wilson loop over whole bulk.
redBase.template iterateOverBulk<All, HaloDepth>(CalcWilsonBresenham<floatT, HaloDepth>(gauge,d_path,rows,columns,z));
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
Wloop /= (GInd::getLatData().globalLattice().mult()*9); //3 loops(one averaged loop per plane) times 3 colors=9.
return Wloop;
}
//
// I am going to paste the hypSmearing.cpp file content here for the moment.
// also call updateAll()
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::Su3Unitarize(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out, Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_base){
HypStaple<floatT, HaloDepth, comp, 4> su_3_unitarize(gauge_out.getAccessor(), gauge_base.getAccessor(), gauge_base.getAccessor(), gauge_base.getAccessor());
gauge_out.iterateOverBulkAllMu(su_3_unitarize);
if(update_all)gauge_out.updateAll();
}
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::SmearAll(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out) {
// create level 1 fields
_dummy.iterateOverBulkAllMu(staple3_lvl1_10);
_gauge_lvl1_10 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_10, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_20);
_gauge_lvl1_20 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_20, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_30);
_gauge_lvl1_30 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_30, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_21);
_gauge_lvl1_21 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_21, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_31);
_gauge_lvl1_31 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_31, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_32);
_gauge_lvl1_32 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_32, _gauge_base);
// now that we have level 1 fields, create level 2 staples
// note: the order of the gauge fields goes in ascending order (10 < 20 < 30, 10 < 21 < 31, 20 < 21 < 32, 30 < 31 < 32)
// this is ASSUMED by HypStaple<floatT, HaloDepth, comp, 2>; DO NOT change this order without also modifying HypStaple<floatT, HaloDepth, comp, 2> and threeLinkStaple_second_level
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_0(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_20.getAccessor(), _gauge_lvl1_30.getAccessor(), _dummy.getAccessor(), 0);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_1(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_31.getAccessor(), _dummy.getAccessor(), 1);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_2(_gauge_lvl1_20.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 2);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_3(_gauge_lvl1_30.getAccessor(), _gauge_lvl1_31.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 3);
//second level fields
_dummy.iterateOverBulkAllMu(staple3_lvl2_0);
_gauge_lvl2_0 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_0, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_1);
_gauge_lvl2_1 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_1, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_2);
_gauge_lvl2_2 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_2, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_3);
_gauge_lvl2_3 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_3, _gauge_base);
// now that we have level 2 fields, create level 3 staple
HypStaple<floatT, HaloDepth, comp, 1> staple3_lvl3(_gauge_lvl2_0.getAccessor(), _gauge_lvl2_1.getAccessor(), _gauge_lvl2_2.getAccessor(), _gauge_lvl2_3.getAccessor());
_dummy.iterateOverBulkAllMu(staple3_lvl3);
// OLD VERSION, (MAYBE) DOES NOT WORK FOR SOME REASON
//_gauge_lvl2_0 = (1-params.alpha_1) * _gauge_base + params.alpha_1/6 * _dummy; //reused _gauge_lvl2_0
//Su3Unitarize(_gauge_lvl2_0);
// NEW VERSION, USES EXTRA FIELD RATHER THAN REUSE _gauge_lvl2_0
gauge_out = (1-params.alpha_1) * _gauge_base + params.alpha_1/6 * _dummy; //reused _gauge_lvl2_0
Su3Unitarize(gauge_out, _gauge_base);
}
#define CLASS_INIT(floatT,HALO) \
template class HypSmearing<floatT,true,HALO,R18>;
INIT_PH(CLASS_INIT)
std::vector<std::vector<int>> Bresenham(int x0,int y0, int x1, int y1)
{
std::vector<std::vector<int>> v;
std::vector<int> x_values;
std::vector<int> y_values;
int dx=x1-x0;
int dy=y1-y0;
int swap_flag=0;
if (dy>dx){
int c=x1;
x1=y1;
y1=c;
//swapping x, y for both the points.
int d=x0;
x0=y0;
y0=d;
dx=x1-x0;
dy=y1-y0;
swap_flag=1;
}
int p=2*dy-dx;
int x=x0;
int y=y0;
x_values.push_back(x);
y_values.push_back(y);
int i=0;
do {
if (p<0) {
x+=1;
p+=2*dy;
}
else{
x+=1;
y+=1;
p+=2*(dy-dx);
}
/* cout<<x<<" "<<y<<std::endl; */
x_values.push_back(x);
y_values.push_back(y);
i+=1;
} while (i<dx);
if(swap_flag==1){
v.push_back(y_values);
v.push_back(x_values);
}
else{
v.push_back(x_values);
v.push_back(y_values);
}
return v;
}
std::vector<std::vector<int>> path_index(std::vector<std::vector<int>> xy_bresenham_vec)
{
std::vector<std::vector<int>> path_vec;
std::vector<int> path_temp;
std::vector<int> index; //index where dx=dy
// 0 means step in x, 1 means step in y, 2 means step in z
//first we shall find the number of points where delta_x=delta_y, and the number of path shall be 2^number(where delta_x=delta_y)
for (int j=0;j<xy_bresenham_vec[0].size()-1;j++)
{
int delta_x=xy_bresenham_vec[0][j+1]-xy_bresenham_vec[0][j];
int delta_y=xy_bresenham_vec[1][j+1]-xy_bresenham_vec[1][j];
if(delta_x==1 and delta_y==0)
{
path_temp.push_back(0);
}
else if(delta_x==0 and delta_y==1)
{
path_temp.push_back(1);
}
else if(delta_x==1 and delta_y==1)
{
path_temp.push_back(0);
path_temp.push_back(1);
index.push_back(j);
}
}
path_vec.push_back(path_temp);
path_vec.push_back(index);
return path_vec;
}
int invert(int num)
{
int result;
if(num==0)
{
result=1;
}
else if(num==1)
{
result=0;
}
return result;
}
std::vector<int> int_to_binary(int num,int size)
{
//we size=size binary vector
std::vector<int> binary_vec;
std::vector<int> binary_vec_rev;
int result=num;
while(result!=1 && result!=0)
{
binary_vec.push_back(result%2);
result/=2;
}
if(result==1 or result==0) binary_vec.push_back(result);
int present_size=binary_vec.size();
if (present_size<size)
{
for(int i=0;i<size-present_size;i++)
{
binary_vec.push_back(0);
}
}
int correct_size=binary_vec.size();
for(int i=0;i<size;i++)
{
binary_vec_rev.push_back(binary_vec[size-1-i]);
}
return binary_vec_rev;
}
int main(int argc, char *argv[]) {
stdLogger.setVerbosity(DEBUG);
/// Initialize parameter class. This class can also read parameter from textfiles!
LatticeParameters param;
/// Initialize the Lattice dimension
const int LatDim[] = {32, 32, 32, 8};
const int NodeDim[] = {1, 1, 1, 1};
/// Just pass these dimensions to the parameter class
param.latDim.set(LatDim);
param.nodeDim.set(NodeDim);
/// Initialize a timer
StopWatch<true> timer;
/// Initialize the CommunicationBase. This class handles the communitation between different Cores/GPU's.
CommunicationBase commBase(&argc, &argv, true);
commBase.init(param.nodeDim());
const size_t HaloDepth = 1; //since it is not multi-gpu code
/// highlight the output differently.
rootLogger.info("Initialize Lattice");
/// Initialize the Indexer on GPU and CPU.
initIndexer(HaloDepth,param,commBase);
typedef GIndexer<All,HaloDepth> GInd;
rootLogger.info("Initialize Gaugefield");
Gaugefield<PREC, true,HaloDepth> gauge(commBase);
typedef double floatT;
GaugeAction<floatT, true, HaloDepth, R18> gaugeaction(gauge);
/// Initialize gaugefield with unity-matrices.
gauge.one();
/// Initialize LatticeContainer. This is in principle the "array", where the values of the plaquette are
/// stored which are summed up in the end
LatticeContainer<true,PREC> redBase(commBase);
/// We need to tell the Reductionbase how large our Array will be
redBase.adjustSize(GInd::getLatData().vol4);
int eqm_point=155;
int last_point=1150;
int node;
sscanf(argv[1], "%d", &node);
//sscanf(argv[2], "%d", &height);
for (int i=eqm_point; i<=last_point;i+=5)
{
rootLogger.info("Read configuration");
gauge.readconf_nersc("/root/project1/build_SIMULATeQCD/applications/try_20_output_dir/after_eqm/node"+std::to_string(node)+"/l328f21b6285m0039185m0783706a_"+std::to_string(node)+"."+std::to_string(i));
gauge.updateAll();
Gaugefield<PREC, true,HaloDepth> gauge_out(commBase);
gauge_out.one(); //initialize the gauge field variable
for(int no_smearing_steps=0;no_smearing_steps<=60;no_smearing_steps+=30)
{
//we shall make a variable which stores gauge configuration
//doing smearing 5,10,15.... times, (loop will run just 5 times where we left of using gauge=gauge_out)
if (no_smearing_steps!=0){
for(int j=1;j<=30;j++)
{
HypSmearing<floatT, true, HaloDepth, R18> hypsmearing(gauge); //smearing of gauge
hypsmearing.SmearAll(gauge_out); //and storing at gauge_out
gauge=gauge_out;
}
}
const int Ns= LatDim[0];
for (int length=1; length<Ns; length++)
{
for (int height=1; height<Ns; height++)
{
PREC Wloop = 0;
/* timer.start(); */
/// compute wilson loop with smeared gauge_out
Wloop = WilsonLoop<PREC,HaloDepth>(gauge, redBase, length,height );
printf("Node:%d,Smear_count:%d,Length=%d,Height=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,length, height,Wloop);
/* rootLogger.info("Reduced Rectangle from rhmc : " , gaugeaction.rectangle()); */
/* rootLogger.info("Reduced plaquette from rhmc : " , gaugeaction.plaquette()); */
}
}
// we shall now implement bresenham algorithm to calculate the new wilson loops.int x_coord=length;
for (int x=1;x<Ns;x++)
{
for(int y=1;y<Ns;y++)
{
std::vector<std::vector<int>> xy_bresenham_vec=Bresenham(0,0,x,y);
std::vector<int> path_xy= path_index(xy_bresenham_vec)[0];
std::vector<std::vector<int>> final_paths_vec;
std::vector<int> index_xy = path_index(xy_bresenham_vec)[1];
int size_dx_eq_dy= index_xy.size();
int n_paths=pow(2,size_dx_eq_dy);
//we shall find the binary eqvivalent of the 2^size_dx_eq_dy and replace the corresponding value by the binary equivalent value at that index
//construct a function which returns the binary equivalent array of any integer
for (int ip=0;ip<n_paths;ip++)
{
std::vector<int> binary_rep = int_to_binary(ip,size_dx_eq_dy);
std::vector<int> path_xy_duplicate=path_xy;
for(int iindex=0;iindex<size_dx_eq_dy;iindex++)
{
int index=iindex+index_xy[iindex];
path_xy_duplicate[index]=binary_rep[iindex];
path_xy_duplicate[index+1]=invert(binary_rep[iindex]);
}
final_paths_vec.push_back(path_xy_duplicate);
}
int rows=n_paths;
int columns=final_paths_vec[0].size();
//allocating memory for paths on gpu
//we can not copy the 2d array to GPU, we have to convert it to 1d array and then decode it to our requirements on GPU
//let us convert this 2d array to 1d,
int *host_final_paths=new int[rows*columns];
for(int ir=0;ir<rows;ir++)
{
for(int ic=0;ic<columns;ic++)
{
host_final_paths[ic+ir*columns]=final_paths_vec[ir][ic];
}
}
//allocating memory for paths on gpu
int* d_path; //declaring pointer to pointer
cudaMalloc((void**)&d_path,rows*columns*sizeof(int));
cudaMemcpy(d_path,host_final_paths,columns*rows*sizeof(int),cudaMemcpyHostToDevice);
for(int z=1;z<Ns;z++)
{
PREC Wloop=WilsonLoopBresenham<PREC,HaloDepth>(gauge, redBase,d_path,rows,columns,z );
printf("Node:%d,Smear_count:%d,x=%d,y=%d,z=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,x,y,z,Wloop);
}
cudaFree(d_path);
delete [] host_final_paths;
}
}
}
}
return 0;
}
//