clover + test (valence)

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -86,8 +86,8 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   CloverTerm += fillCloverXY(Bz);
   CloverTerm += fillCloverXT(Ex);
   CloverTerm += fillCloverYT(Ey);
-  CloverTerm += fillCloverZT(Ez) ;
+  CloverTerm += fillCloverZT(Ez);
-  CloverTerm *= csw;
+  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<Impl>::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<Impl>::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<Impl>::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<Impl>::Mooee(const FermionField &in, FermionField &out)
 template <class Impl>
 void WilsonCloverFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out)
 {
-  this->MooeeInternal(in, out, DaggerNo, InverseYes);
+  this->MooeeInternal(in, out, DaggerYes, InverseNo);
 }
 
 template <class Impl>
@@ -154,7 +164,7 @@ void WilsonCloverFermion<Impl>::MooeeInv(const FermionField &in, FermionField &o
 template <class Impl>
 void WilsonCloverFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
 {
-  this->MooeeInternal(in, out, DaggerNo, InverseYes);
+  this->MooeeInternal(in, out, DaggerYes, InverseYes);
 }
 
 template <class Impl>
@@ -164,26 +174,98 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
   CloverFieldType *Clover;
   assert(in.checkerboard == Odd || in.checkerboard == Even);
 
-  if (in._grid->_isCheckerBoarded)
+
-  {
+
-    if (in.checkerboard == Odd)
+
-    {
+ if (dag){ 
-      std::cout << "Calling clover term Odd" << std::endl;
+  if (in._grid->_isCheckerBoarded){
-      Clover = (inv) ? &CloverTermInvOdd : &CloverTermOdd;
+    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<int> 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<int> 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 << GridLogMessage << "*Clover.checkerboard "  << (*Clover).checkerboard << std::endl;
-    {
+    out = *Clover * in;
-      std::cout << "Calling clover term Even" << std::endl;
+  } else { 
-      Clover = (inv) ? &CloverTermInvEven : &CloverTermEven;
+  Clover = (inv) ? &CloverTermInv : &CloverTerm; 
-    }
+  out = adj(*Clover) * in;
-  }
-  else
-  {
-    Clover = (inv) ? &CloverTermInv : &CloverTerm;
   }
 
-  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

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]()(0, 1) = timesI(F._odata[i]()());           //fixed
-      T._odata[i]()(1, 0) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(1, 0) = timesI(F._odata[i]()());           //fixed
-      T._odata[i]()(2, 3) = timesI(F._odata[i]()());
+      T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()());      //fixed
-      T._odata[i]()(3, 2) = timesI(F._odata[i]()());
+      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]()(0, 1) = -(F._odata[i]()());         //fixed
-      T._odata[i]()(1, 0) = -(F._odata[i]()());
+      T._odata[i]()(1, 0) = (F._odata[i]()());          //fixed
-      T._odata[i]()(2, 3) = -(F._odata[i]()());
+      T._odata[i]()(2, 3) = (F._odata[i]()());          //fixed
-      T._odata[i]()(3, 2) = (F._odata[i]()());
+      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]()(0, 0) = timesI(F._odata[i]()());           //fixed
-      T._odata[i]()(1, 1) = timesI(F._odata[i]()());
+      T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()());      //fixed
-      T._odata[i]()(2, 2) = timesI(F._odata[i]()());
+      T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()());      //fixed
-      T._odata[i]()(3, 3) = timesMinusI(F._odata[i]()());
+      T._odata[i]()(3, 3) = timesI(F._odata[i]()());           //fixed
     }
   
   return T;

tests/core/Test_wilson_clover.cc

@@ -172,11 +172,6 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
   std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
 
-  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
-  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
-  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
-  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
-
   std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
   std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
 
@@ -196,11 +191,12 @@ int main (int argc, char ** argv)
   setCheckerboard(phi,phi_e);
   setCheckerboard(phi,phi_o);
 
+  std::cout<<GridLogMessage << "err   "<< norm2(err)<< std::endl;
   err = phi-chi;
   std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
 
   std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1                                   "<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MeeDag MeeInvDag = 1                                   "<<std::endl;
   std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
 
   pickCheckerboard(Even,chi_e,chi);
@@ -215,37 +211,29 @@ int main (int argc, char ** argv)
   setCheckerboard(phi,phi_e);
   setCheckerboard(phi,phi_o);
 
+ std::cout<<GridLogMessage << "err   "<< std::endl;
   err = phi-chi;
   std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
 
   std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
-  std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian              "<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MeeInv MeeDag = 1                                      "<<std::endl;
   std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
   
-  random(pRNG,phi);
+ pickCheckerboard(Even,chi_e,chi);
-  random(pRNG,chi);
+ pickCheckerboard(Odd ,chi_o,chi);
-  pickCheckerboard(Even,chi_e,chi);
-  pickCheckerboard(Odd ,chi_o,chi);
-  pickCheckerboard(Even,phi_e,phi);
-  pickCheckerboard(Odd ,phi_o,phi);
 
-  SchurDiagMooeeOperator<WilsonCloverFermionR,FermionField> HermOpEO(Dwc);
+  Dwc.MooeeDag(chi_e,src_e);
-  HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
+  Dwc.MooeeInv(src_e,phi_e);
-  HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
 
-  HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
+  Dwc.MooeeDag(chi_o,src_o);
-  HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
+  Dwc.MooeeInv(src_o,phi_o);
 
-  pDce = innerProduct(phi_e,dchi_e);
+  setCheckerboard(phi,phi_e);
-  pDco = innerProduct(phi_o,dchi_o);
+  setCheckerboard(phi,phi_o);
-  cDpe = innerProduct(chi_e,dphi_e);
-  cDpo = innerProduct(chi_o,dphi_o);
-
-  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
-  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
-
-  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
-  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
 
+ std::cout<<GridLogMessage << "err   "<< std::endl;
+  err = phi-chi;
+  std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
+  
   Grid_finalize();
 }