diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index f8b62ba4..54e72353 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -86,8 +86,8 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) CloverTerm += fillCloverXY(Bz); CloverTerm += fillCloverXT(Ex); CloverTerm += fillCloverYT(Ey); - CloverTerm += fillCloverZT(Ez) ; - CloverTerm *= csw; + CloverTerm += fillCloverZT(Ez); + CloverTerm *= 0.5 * csw; // FieldStrength normalization? should be ( -i/8 ). Is it the anti-symmetric combination? int lvol = _Umu._grid->lSites(); int DimRep = Impl::Dimension; @@ -109,7 +109,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) for (int a = 0; a < DimRep; a++) for (int b = 0; b < DimRep; b++) EigenCloverOp(a + j * DimRep, b + k * DimRep) = Qx()(j, k)(a, b); - //std::cout << EigenCloverOp << std::endl; + // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; EigenInvCloverOp = EigenCloverOp.inverse(); @@ -119,6 +119,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) for (int a = 0; a < DimRep; a++) for (int b = 0; b < DimRep; b++) Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); + // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; pokeLocalSite(Qxinv, CloverTermInv, lcoor); } @@ -127,8 +128,17 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) pickCheckerboard(Even, CloverTermEven, CloverTerm); pickCheckerboard( Odd, CloverTermOdd, CloverTerm); + + pickCheckerboard(Even, CloverTermDagEven, adj(CloverTerm)); + pickCheckerboard( Odd, CloverTermDagOdd, adj(CloverTerm)); + + pickCheckerboard(Even, CloverTermInvEven, CloverTermInv); pickCheckerboard( Odd, CloverTermInvOdd, CloverTermInv); + + + pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); + pickCheckerboard( Odd, CloverTermInvDagOdd, adj(CloverTermInv)); } @@ -142,7 +152,7 @@ void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) template void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out) { - this->MooeeInternal(in, out, DaggerNo, InverseYes); + this->MooeeInternal(in, out, DaggerYes, InverseNo); } template @@ -154,7 +164,7 @@ void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &o template void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) { - this->MooeeInternal(in, out, DaggerNo, InverseYes); + this->MooeeInternal(in, out, DaggerYes, InverseYes); } template @@ -164,26 +174,98 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie CloverFieldType *Clover; assert(in.checkerboard == Odd || in.checkerboard == Even); - if (in._grid->_isCheckerBoarded) - { - if (in.checkerboard == Odd) - { - std::cout << "Calling clover term Odd" << std::endl; - Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd; + + + + if (dag){ + if (in._grid->_isCheckerBoarded){ + if (in.checkerboard == Odd){ + std::cout << "Calling clover term adj Odd" << std::endl; + Clover = (inv) ? &CloverTermInvDagOdd : &CloverTermDagOdd; + +/* test + int DimRep = Impl::Dimension; + Eigen::MatrixXcd A = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + std::vector lcoor; + typename SiteCloverType::scalar_object Qx2 = zero; + GridBase *grid = in._grid; + int site = 0 ; + grid->LocalIndexToLocalCoor(site, lcoor); + peekLocalSite(Qx2, *Clover, lcoor); + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + A(a + j * DimRep, b + k * DimRep) = Qx2()(j, k)(a, b); + std::cout << "adj Odd =" << site << "\n" << A << std::endl; + end test */ + + + + } else { + std::cout << "Calling clover term adj Even" << std::endl; + Clover = (inv) ? &CloverTermInvDagEven : &CloverTermDagEven; + +/* test + int DimRep = Impl::Dimension; + Eigen::MatrixXcd A = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + std::vector lcoor; + typename SiteCloverType::scalar_object Qx2 = zero; + GridBase *grid = in._grid; + int site = 0 ; + grid->LocalIndexToLocalCoor(site, lcoor); + peekLocalSite(Qx2, *Clover, lcoor); + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + A(a + j * DimRep, b + k * DimRep) = Qx2()(j, k)(a, b); + std::cout << "adj Odd =" << site << "\n" << A << std::endl; + end test */ + + } - else - { - std::cout << "Calling clover term Even" << std::endl; - Clover = (inv) ? &CloverTermInvEven : &CloverTermEven; - } - } - else - { - Clover = (inv) ? &CloverTermInv : &CloverTerm; + std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; + out = *Clover * in; + } else { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = adj(*Clover) * in; } - std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; - if (dag){ out = adj(*Clover) * in;} else { out = *Clover * in;} + + + + } else { + if (in._grid->_isCheckerBoarded){ + + if (in.checkerboard == Odd){ + std::cout << "Calling clover term Odd" << std::endl; + Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd; + } else { + std::cout << "Calling clover term Even" << std::endl; + Clover = (inv) ? &CloverTermInvEven : &CloverTermEven; + } + out = *Clover * in; + std::cout << GridLogMessage << "*Clover.checkerboard " << (*Clover).checkerboard << std::endl; + } else { + Clover = (inv) ? &CloverTermInv : &CloverTerm; + out = *Clover * in; + } + } + + + + + + + +/* + } else { + out = *Clover * in; + } + */ + + } // MooeeInternal // Derivative parts diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index fd9d1f60..e8654513 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -63,7 +63,11 @@ public: CloverTermEven(&Hgrid), CloverTermOdd(&Hgrid), CloverTermInvEven(&Hgrid), - CloverTermInvOdd(&Hgrid) + CloverTermInvOdd(&Hgrid), + CloverTermDagEven(&Hgrid), //test + CloverTermDagOdd(&Hgrid), //test + CloverTermInvDagEven(&Hgrid), //test + CloverTermInvDagOdd(&Hgrid) //test { csw = _csw; assert(Nd == 4); // require 4 dimensions @@ -91,6 +95,11 @@ private: CloverFieldType CloverTerm, CloverTermInv; // Clover term CloverFieldType CloverTermEven, CloverTermOdd; CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term + + CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; //test + CloverFieldType CloverTermDagEven, CloverTermDagOdd; //test + + // eventually these two can be compressed into 6x6 blocks instead of the 12x12 // using the DeGrand-Rossi basis for the gamma matrices @@ -149,10 +158,10 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { - T._odata[i]()(0, 1) = timesMinusI(F._odata[i]()()); - T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()()); - T._odata[i]()(2, 3) = timesI(F._odata[i]()()); - T._odata[i]()(3, 2) = timesI(F._odata[i]()()); + T._odata[i]()(0, 1) = timesI(F._odata[i]()()); //fixed + T._odata[i]()(1, 0) = timesI(F._odata[i]()()); //fixed + T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); //fixed + T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); //fixed } return T; @@ -165,10 +174,10 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { - T._odata[i]()(0, 1) = (F._odata[i]()()); - T._odata[i]()(1, 0) = -(F._odata[i]()()); - T._odata[i]()(2, 3) = -(F._odata[i]()()); - T._odata[i]()(3, 2) = (F._odata[i]()()); + T._odata[i]()(0, 1) = -(F._odata[i]()()); //fixed + T._odata[i]()(1, 0) = (F._odata[i]()()); //fixed + T._odata[i]()(2, 3) = (F._odata[i]()()); //fixed + T._odata[i]()(3, 2) = -(F._odata[i]()()); //fixed } return T; @@ -181,10 +190,10 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { - T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()()); - T._odata[i]()(1, 1) = timesI(F._odata[i]()()); - T._odata[i]()(2, 2) = timesI(F._odata[i]()()); - T._odata[i]()(3, 3) = timesMinusI(F._odata[i]()()); + T._odata[i]()(0, 0) = timesI(F._odata[i]()()); //fixed + T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()()); //fixed + T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); //fixed + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); //fixed } return T; diff --git a/tests/core/Test_wilson_clover.cc b/tests/core/Test_wilson_clover.cc index 1b208e2d..9e5b246e 100644 --- a/tests/core/Test_wilson_clover.cc +++ b/tests/core/Test_wilson_clover.cc @@ -172,11 +172,6 @@ int main (int argc, char ** argv) std::cout< HermOpEO(Dwc); - HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); - HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); + Dwc.MooeeDag(chi_e,src_e); + Dwc.MooeeInv(src_e,phi_e); - HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); - HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); + Dwc.MooeeDag(chi_o,src_o); + Dwc.MooeeInv(src_o,phi_o); - pDce = innerProduct(phi_e,dchi_e); - pDco = innerProduct(phi_o,dchi_o); - cDpe = innerProduct(chi_e,dphi_e); - cDpo = innerProduct(chi_o,dphi_o); - - std::cout<