template<class floatT,size_t HaloDepth,CompressionType comp>
__host__ __device__ SU3<floatT> inline hypThreeLinkStaple_first_level(SU3Accessor<floatT,comp> gAcc, gSiteMu siteMu, int excluded_dir1, int excluded_dir2) {
    typedef GIndexer<All,HaloDepth> GInd;
    SU3<floatT> temp = su3_zero<floatT>();
    /* printf("ex_dir 1 %d , ex_dir2 %d  \n",excluded_dir1,excluded_dir2); */
    gSite site = GInd::getSite(siteMu.isite);
    gSite site_save = site;
    int mu = siteMu.mu;
    if(mu == excluded_dir1 || mu == excluded_dir2) return temp*(floatT)NAN;
    if(excluded_dir1 == excluded_dir2) return temp*(floatT)NAN;

    int check_count = 0;

    for (int nu_h = 1; nu_h < 4; nu_h++) {
        int nu = (mu+nu_h)%4;
        if(nu == excluded_dir1 || nu == excluded_dir2) continue;
        check_count += 1;
        temp += gAcc.template getLinkPath<All, HaloDepth>(site, nu, mu, Back(nu));
        site = site_save;
        temp += gAcc.template getLinkPath<All, HaloDepth>(site, Back(nu), mu, nu);
        site = site_save;
    }

    if(excluded_dir1 >= 0 || excluded_dir2 >= 0) assert(check_count == 1); //loop length is 1 in the case of use in hyp smearing
    if(excluded_dir1 >= 0 || excluded_dir2 >= 0) assert(excluded_dir1 >= 0 && excluded_dir2 >= 0); //only one excluded dir not supported

    return temp;
}


//It is assumed that the exluded directions in gaAcc_0, gaAcc_1, gaAcc_2 go in asending order (10 < 20 < 21 < 30 < 31 < 32)
// example:  suppose the exluded directions are 20, 21, 32, (for Vtilde_{i, mu; 2}), then gAcc_0 = _gauge_lvl1_20, gAcc_1 = _gauge_lvl1_21, gAcc_2 = _gauge_lvl1_32
// in other words, the primary excluded direction in this example is 2, so the extra excluded directions go in asending order:  0, 1, 3.
template<class floatT,size_t HaloDepth,CompressionType comp>
__host__ __device__ SU3<floatT> inline hypThreeLinkStaple_second_level(SU3Accessor<floatT,comp> gAcc_0, SU3Accessor<floatT,comp> gAcc_1, SU3Accessor<floatT,comp> gAcc_2, gSiteMu siteMu, int excluded_dir, SU3Accessor<floatT,comp> temp_gAcc_mu_excluded_dir, SU3Accessor<floatT,comp> temp_gAcc_nu_excluded_dir) {

    typedef GIndexer<All,HaloDepth> GInd;
    SU3<floatT> temp = su3_zero<floatT>();
    int mu = siteMu.mu;
    if(mu==excluded_dir) return temp*(floatT)NAN;
    gSite origin = GInd::getSite(siteMu.isite);
    gSite upMu = GInd::site_up(origin,mu);
    //
//gAcc_idx = mu > excluded_dir ? mu - 1 : mu;

    //ex = 0, then mu=1,2,3, => 10, 20, 30, then gAcc_idx= (0, 1, 2) = (mu-1, mu-1, mu-1)
    //ex = 1, then mu=0,2,3, => 10, 21, 31, then gAcc_idx= (0, 1, 2) = (mu, mu-1, mu-1)
    //ex = 2, then mu=0,1,3, => 20, 21, 32, then gAcc_idx= (0, 1, 2) = (mu, mu, mu-1)
    //ex = 3, then mu=0,1,2, => 30, 31, 32, then gAcc_idx= (0, 1, 2) = (mu, mu, mu)
    //=>gAcc_idx = mu > excluded_dir ? mu - 1 : mu

    int gAcc_idx = mu > excluded_dir ? mu - 1 : mu;
    assert(gAcc_idx < 3);
    if(gAcc_idx == 0) temp_gAcc_mu_excluded_dir=gAcc_0; // mu is temp excluded dir
    else if(gAcc_idx == 1)temp_gAcc_mu_excluded_dir=gAcc_1;
    else if(gAcc_idx == 2)temp_gAcc_mu_excluded_dir=gAcc_2;
    else assert(0);

    int nu_add = 0;
    for (int nu_h = 0; nu_h < 3; nu_h++) {

        if(nu_h == excluded_dir) nu_add = 1;
        int nu = nu_h + nu_add; // we want nu here to go in ascending order; if excluded_dir == 1, then nu should go 0, 2, 3
        if(nu == mu) continue;

        gSite downNu = GInd::site_dn(origin,nu);
        gSite upNu = GInd::site_up(origin,nu);
        gSite upMudownNu = GInd::site_dn(upMu, nu);
        assert(gAcc_idx != nu_h);

        if(nu_h == 0)temp_gAcc_nu_excluded_dir=gAcc_0; //nu is temp excluded dir
        else if(nu_h == 1)temp_gAcc_nu_excluded_dir=gAcc_1;
        else if(nu_h == 2)temp_gAcc_nu_excluded_dir=gAcc_2;

        // nu > 0
        temp += temp_gAcc_mu_excluded_dir.getLink(GInd::getSiteMu(origin,nu))*temp_gAcc_nu_excluded_dir.getLink(GInd::getSiteMu(upNu,mu))*temp_gAcc_mu_excluded_dir.getLinkDagger(GInd::getSiteMu(upMu,nu));
        // nu < 0
        temp += temp_gAcc_mu_excluded_dir.getLinkDagger(GInd::getSiteMu(downNu,nu))*temp_gAcc_nu_excluded_dir.getLink(GInd::getSiteMu(downNu,mu))*temp_gAcc_mu_excluded_dir.getLink(GInd::getSiteMu(upMudownNu,nu));
    }
    return temp;
}


template<class floatT,size_t HaloDepth,CompressionType comp>
__host__ __device__ SU3<floatT> inline hypThreeLinkStaple_third_level(SU3Accessor<floatT,comp> gAcc_0, SU3Accessor<floatT,comp> gAcc_1, SU3Accessor<floatT,comp> gAcc_2, SU3Accessor<floatT,comp> gAcc_3, gSiteMu siteMu, SU3Accessor<floatT,comp> temp_gAcc_mu, SU3Accessor<floatT,comp> temp_gAcc_nu) {

    typedef GIndexer<All,HaloDepth> GInd;
    SU3<floatT> temp = su3_zero<floatT>();
    SU3<floatT> temp_chk = su3_zero<floatT>();
    int mu = siteMu.mu;
    gSite origin = GInd::getSite(siteMu.isite);
    gSite upMu = GInd::site_up(origin,mu);

    if(mu == 0)temp_gAcc_mu=gAcc_0;
    else if(mu == 1)temp_gAcc_mu=gAcc_1;
    else if(mu == 2)temp_gAcc_mu=gAcc_2;
    else if(mu == 3)temp_gAcc_mu=gAcc_3;

    for (int nu_h = 1; nu_h < 4; nu_h++) {
        int nu = (mu+nu_h)%4;
        gSite downNu = GInd::site_dn(origin,nu);
        gSite upNu = GInd::site_up(origin,nu);
        gSite upMudownNu = GInd::site_dn(upMu, nu);

        if(nu == 0)temp_gAcc_nu=gAcc_0;
        else if(nu == 1)temp_gAcc_nu=gAcc_1;
        else if(nu == 2)temp_gAcc_nu=gAcc_2;
        else if(nu == 3)temp_gAcc_nu=gAcc_3;

        // nu > 0
        temp += temp_gAcc_mu.getLink(GInd::getSiteMu(origin,nu))*temp_gAcc_nu.getLink(GInd::getSiteMu(upNu,mu))*temp_gAcc_mu.getLinkDagger(GInd::getSiteMu(upMu,nu));
        // nu < 0
        temp += temp_gAcc_mu.getLinkDagger(GInd::getSiteMu(downNu,nu))*temp_gAcc_nu.getLink(GInd::getSiteMu(downNu,mu))*temp_gAcc_mu.getLink(GInd::getSiteMu(upMudownNu,nu));
    }
    assert(!(temp==temp_chk));
    return temp;
}