Contents

Codes

code1


SU3<floatT> temp;
gSite site1;
site1=GInd::site_up(site,mu);
temp = SU3Accessor.template getLinkPath<All, HaloDepth>(site1,mu, nu,Back(mu), Back(mu), Back(nu));
temp += SU3Accessor.template getLinkPath<All, HaloDepth>(site1,nu, nu,Back(mu), Back(nu), Back(nu));
result += tr_d(temp,SU3Accessor.getLink(GInd::getSiteMu(site, mu)));

result1

Reduced Rectangle is: 0.139255

code2


SU3<floatT> temp;
temp = SU3Accessor.template getLinkPath<All, HaloDepth>(site,mu,mu, nu,Back(mu), Back(mu), Back(nu));
temp += SU3Accessor.template getLinkPath<All, HaloDepth>(site,mu,nu, nu,Back(mu), Back(nu), Back(nu));
result += tr_d(temp);

result2

Reduced Rectangle is :269532

code3


SU3<floatT> temp;
SU3<floatT> link;
gSite site1;
site1=GInd::site_up(site,mu);
temp = SU3Accessor.template getLinkPath<All, HaloDepth>(site1,mu, nu,Back(mu), Back(mu), Back(nu));
temp += SU3Accessor.template getLinkPath<All, HaloDepth>(site1,nu, nu,Back(mu), Back(nu), Back(nu));
link=SU3Accessor.getLink(GInd::getSiteMu(site, mu));
result += tr_d(temp,link);

result3

Reduced Rectangle is: 0.139255

code4


SU3<floatT> temp;
SU3<floatT> link;
gSite site1;
site1=GInd::site_up(site,mu);
temp = SU3Accessor.template getLinkPath<All, HaloDepth>(site1,mu, nu,Back(mu), Back(mu), Back(nu));
temp += SU3Accessor.template getLinkPath<All, HaloDepth>(site1,nu, nu,Back(mu), Back(nu), Back(nu));
link=SU3Accessor.getLink(GInd::getSiteMu(site, mu));
link=dagger(link);
result += tr_d(temp,link);

result4

Reduced Rectangle is: 0.000282382

code5


SU3<floatT> temp;
gSite site1;
gSite site2;
gSite site3;
gSite site4;
gSite site5;
gSite site6;
gSite site8;
site1=site;
site2=GInd::site_up(site1,mu);
site3=GInd::site_up(site2,mu);
site4=GInd::site_up(site3,nu);
site5=GInd::site_up(site2,nu);
site6=GInd::site_up(site1,nu);
site8=GInd::site_up(site6,nu);
//temp=SU3Accessor.getLink(GInd::getSiteMu(site1,mu))*SU3Accessor.getLink(GInd::getSiteMu(site2,mu))*SU3Accessor.getLink(GInd::getSiteMu(site3,nu))*dagger(SU3Accessor.getLink(GInd::getSiteMu(site5,mu)))*dagger(SU3Accessor.getLink(GInd::getSiteMu(site6,mu)))*dagger(SU3Accessor.getLink(GInd::getSiteMu(site,nu)));

temp=SU3Accessor.getLink(GInd::getSiteMu(site1,mu))
    *SU3Accessor.getLink(GInd::getSiteMu(site2,mu))
    *SU3Accessor.getLink(GInd::getSiteMu(site3,nu))
    *dagger(SU3Accessor.getLink(GInd::getSiteMu(site5,mu)))
    *dagger(SU3Accessor.getLink(GInd::getSiteMu(site6,mu)))
    *dagger(SU3Accessor.getLink(GInd::getSiteMu(site,nu)));

temp += SU3Accessor.getLink(GInd::getSiteMu(site1,mu))*SU3Accessor.getLink(GInd::getSiteMu(site2,nu))*SU3Accessor.getLink(GInd::getSiteMu(site5,nu))*dagger(SU3Accessor.getLink(GInd::getSiteMu(site8,mu)))*dagger(SU3Accessor.getLink(GInd::getSiteMu(site6,nu)))*dagger(SU3Accessor.getLink(GInd::getSiteMu(site,nu)));
result +=tr_d(temp);

result5

Reduced Rectangle is: 0.269532

code6


SU3<floatT> temp;
gSite site1;
gSite site2;
gSite site3;
gSite site4;
gSite site5;
gSite site6;
gSite site8;
site1=site;
site2=GInd::site_up(site1,mu);
site3=GInd::site_up(site2,mu);
site4=GInd::site_up(site3,nu);
site5=GInd::site_up(site2,nu);
site6=GInd::site_up(site1,nu);
site8=GInd::site_up(site6,nu);
temp=SU3Accessor.getLink(GInd::getSiteMu(site1,mu))*SU3Accessor.getLink(GInd::getSiteMu(site2,mu))*SU3Accessor.getLink(GInd::getSiteMu(site3,nu))*(SU3Accessor.getLinkDagger(GInd::getSiteMu(site5,mu)))*(SU3Accessor.getLinkDagger(GInd::getSiteMu(site6,mu)))*SU3Accessor.getLinkDagger(GInd::getSiteMu(site,nu));

temp += SU3Accessor.getLink(GInd::getSiteMu(site1,mu))*SU3Accessor.getLink(GInd::getSiteMu(site2,nu))*SU3Accessor.getLink(GInd::getSiteMu(site5,nu))*(SU3Accessor.getLinkDagger(GInd::getSiteMu(site8,mu)))*(SU3Accessor.getLinkDagger(GInd::getSiteMu(site6,nu)))*(SU3Accessor.getLinkDagger(GInd::getSiteMu(site,nu)));
result +=tr_d(temp);

result6

Reduced Rectangle is: 0.269532

code7

I am copying the code of calculating ractangle from rhmc.cpp 


temp = SU3Accessor.getLink(GInd::getSiteMu(GInd::site_up(site, mu), mu) )
                        * SU3Accessor.getLink(GInd::getSiteMu(GInd::site_2up(site, mu), nu) )
                        * SU3Accessor.getLinkDagger(GInd::getSiteMu(GInd::site_up_up(site, mu, nu), mu))
                        * SU3Accessor.getLinkDagger(GInd::getSiteMu(GInd::site_up(site, nu), mu))
                        * SU3Accessor.getLinkDagger(GInd::getSiteMu(site, nu));
                temp += SU3Accessor.getLink(GInd::getSiteMu(GInd::site_up(site, mu), nu) )
                         * SU3Accessor.getLink(GInd::getSiteMu(GInd::site_up_up(site, mu, nu), nu) )
                         * SU3Accessor.getLinkDagger(GInd::getSiteMu(GInd::site_2up(site, nu), mu))
                         * SU3Accessor.getLinkDagger(GInd::getSiteMu(GInd::site_up(site, nu), nu))
                         * SU3Accessor.getLinkDagger(GInd::getSiteMu(site, nu));

                result += tr_d(SU3Accessor.getLink(GInd::getSiteMu(site, mu)), temp);

result6

Reduced Rectangle is: 0.269532

Remarks

  1. In rhmc calcuations they inherited frm gAcc, rather than from SU3Accessor class.

code8

I am copying the code of calculating ractangle from rhmc.cpp 


temp = gAcc.getLink(GInd::getSiteMu(GInd::site_up(site, mu), mu) )
                        * gAcc.getLink(GInd::getSiteMu(GInd::site_2up(site, mu), nu) )
                        * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_up_up(site, mu, nu), mu))
                        * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_up(site, nu), mu))
                        * gAcc.getLinkDagger(GInd::getSiteMu(site, nu));
                temp += gAcc.getLink(GInd::getSiteMu(GInd::site_up(site, mu), nu) )
                         * gAcc.getLink(GInd::getSiteMu(GInd::site_up_up(site, mu, nu), nu) )
                         * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_2up(site, nu), mu))
                         * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_up(site, nu), nu))
                         * gAcc.getLinkDagger(GInd::getSiteMu(site, nu));

                result += tr_d(gAcc.getLink(GInd::getSiteMu(site, mu)), temp);

I have tried adding more parameter while defining template structure. 


#include "../simulateqcd.h"
#include <string>

#define PREC double

template<class floatT, bool onDevice,size_t HaloDepth, CompressionType comp>
struct CalcPlaq{

    //Gauge accessor to access the gauge field
    /* SU3Accessor<floatT> SU3Accessor; */
    SU3Accessor<floatT,comp> gAcc;

    //Constructor to initialize all necessary members.
    CalcPlaq(Gaugefield<floatT,true,HaloDepth,comp> &gauge) : gAcc(gauge.getAccessor()){
    }

    //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;

        /// Define a SU(3) matrix
        SU3<floatT> temp;
        SU3<floatT> link;

        floatT result = 0;
        /// loop through all directions which are needed by the four link variables in the plaquette
        for (int nu = 1; nu < 4; nu++) {
            for (int mu = 0; mu < nu; mu++) {
              temp = gAcc.getLink(GInd::getSiteMu(GInd::site_up(site, mu), mu) )
                        * gAcc.getLink(GInd::getSiteMu(GInd::site_2up(site, mu), nu) )
                        * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_up_up(site, mu, nu), mu))
                        * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_up(site, nu), mu))
                        * gAcc.getLinkDagger(GInd::getSiteMu(site, nu));
                temp += gAcc.getLink(GInd::getSiteMu(GInd::site_up(site, mu), nu) )
                         * gAcc.getLink(GInd::getSiteMu(GInd::site_up_up(site, mu, nu), nu) )
                         * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_2up(site, nu), mu))
                         * gAcc.getLinkDagger(GInd::getSiteMu(GInd::site_up(site, nu), nu))
                         * gAcc.getLinkDagger(GInd::getSiteMu(site, nu));

                result += tr_d(gAcc.getLink(GInd::getSiteMu(site, mu)), temp);
            }
        }

        return result*0.5;
    }
};


//Function to compute the plaquette using the above struct CalcPlaq.
template<class floatT, size_t HaloDepth,CompressionType comp>
floatT gPlaq(Gaugefield<floatT,true, HaloDepth,comp> &gauge, LatticeContainer<true,floatT> &redBase){

    typedef GIndexer<All,HaloDepth> GInd;
    const size_t elems = GInd::getLatData().vol4;
    std::cout << "Total number of sites are : "<<elems<< std::endl;
    //Make sure, redBase is large enough
    redBase.adjustSize(elems);

   redBase.template iterateOverBulk<All, HaloDepth>(CalcPlaq<floatT,true,HaloDepth,comp>(gauge));

    //Do the final reduction
    floatT plaq;
    redBase.reduce(plaq, elems);

    //Normalize the result
    plaq /= (GInd::getLatData().globalLattice().mult()*18);
    return plaq;
}


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[] = {24, 24, 20, 20}; */
    const int LatDim[] = {24, 24, 24, 6};

    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()); //it will print node dimension 

    /// Set the HaloDepth. It should be a constant values, since this value should be passed as an non-type template
    /// parameter to each kernel.
    const size_t HaloDepth = 1;


    // Initialize the Indexer on GPU and CPU.
    initIndexer(HaloDepth,param,commBase);
    typedef GIndexer<All,HaloDepth> GInd;


    /// Initialize the Gaugefield. Basically, this object holds all SU(3)-matrices of the gaugefield.
    /// The second template parameter determines whether the gaugefield should be stored on GPU or CPU but this will
    /// be changed in the future ...
    rootLogger.info("Initialize Gaugefield");
    Gaugefield<PREC, true,HaloDepth> gauge(commBase);

    /// 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);
    grnd_state<false> h_rand;
    grnd_state<true> d_rand;

    h_rand.make_rng_state(1337); // 1337 is random number generator seed 
    d_rand = h_rand;

//let me write a program to find it for all gague files:stored at "/root/try_16_output_dir/after_eqm/node0/l246f21b6285m0039185m0783706a_0.") starting from 200. in seperation of 10.
    /// Read a configuration from hard drive.
    /// I am incluiding the for loop to do it for all files
    for (int i=250; i<=250;i+=10)
    {
        rootLogger.info("Reading Gauge Configuration");
        gauge.readconf_nersc("/root/try_16_output_dir/after_eqm/node0/l246f21b6285m0039185m0783706a_0."+std::to_string(i));
        /* gauge.random(d_rand.state); */
        /// start timer...
        timer.start();
        /// define variable where the plaquette should be stored.
        PREC plaq = 0;

        /// Exchange Halos before calculating the plaquette! , it gives nan values if we don't exchange halos.
        gauge.updateAll();

        timer.start();
        /// compute plaquette
        plaq = gPlaq<PREC,HaloDepth>(gauge, redBase);
        /// stop timer and print time
        timer.stop();
        rootLogger.info("Time for operators " ,  timer);
        rootLogger.info("Reduced Rectangle is: " ,  plaq);
        /* rootLogger.info("Reduced Plaquette is: " ,  plaq); */
}
    return 0;
}

result8

Reduced Rectangle is: 0.269532