Merge branch 'feature/bugfix-ck-cj' into develop

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -452,10 +452,10 @@ PARALLEL_FOR_LOOP
 	// DhopDir provides U or Uconj depending on coor/flavour.
 	GaugeLinkField link(mat._grid);
 	// use lorentz for flavour as hack.
-	auto tmp = TraceIndex<SpinIndex>(outerProduct(Btilde,A));  
 PARALLEL_FOR_LOOP
-        for(auto ss=tmp.begin();ss<tmp.end();ss++){
+        for(auto ss=link.begin();ss<link.end();ss++){
-	  link[ss]() = tmp[ss](0,0) - conjugate(tmp[ss](1,1)) ; // IS THIS SIGN RIGHT?
+	  auto ttmp = traceIndex<SpinIndex>(outerProduct(Btilde[ss],A[ss]));  
+	  link[ss]() = ttmp(0,0) + conjugate(ttmp(1,1)) ; 
 	}
 	PokeIndex<LorentzIndex>(mat,link,mu);
 	return;

tests/forces/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS += Test_contfrac_force Test_dwf_force Test_dwf_gpforce Test_gp_rect_force Test_gpdwf_force Test_partfrac_force Test_rect_force Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
+bin_PROGRAMS += Test_contfrac_force Test_dwf_force Test_dwf_gpforce Test_gp_rect_force Test_gpdwf_force Test_gpwilson_force Test_partfrac_force Test_rect_force Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi
 
 
 Test_contfrac_force_SOURCES=Test_contfrac_force.cc
@@ -22,6 +22,10 @@ Test_gpdwf_force_SOURCES=Test_gpdwf_force.cc
 Test_gpdwf_force_LDADD=-lGrid
 
 
+Test_gpwilson_force_SOURCES=Test_gpwilson_force.cc
+Test_gpwilson_force_LDADD=-lGrid
+
+
 Test_partfrac_force_SOURCES=Test_partfrac_force.cc
 Test_partfrac_force_LDADD=-lGrid
 

tests/forces/Test_gpwilson_force.cc

@@ -0,0 +1,169 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_gpdwf_force.cc
+
+    Copyright (C) 2015
+
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+#define parallel_for PARALLEL_FOR_LOOP for
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+
+  const int Ls=8;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = UGrid;
+  GridRedBlackCartesian * FrbGrid = UrbGrid;
+
+  std::vector<int> seeds4({1,2,3,4});
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+  typedef typename GparityDomainWallFermionR::FermionField FermionField;
+  FermionField phi        (FGrid); gaussian(RNG4,phi);
+  FermionField Mphi       (FGrid); 
+  FermionField MphiPrime  (FGrid); 
+
+  LatticeGaugeField U(UGrid);
+
+  SU3::HotConfiguration(RNG4,U);
+  
+  ////////////////////////////////////
+  // Unmodified matrix element
+  ////////////////////////////////////
+  RealD mass=0.01; 
+  RealD M5=1.8; 
+
+  const int nu = 3;
+  std::vector<int> twists(Nd,0);  twists[nu] = 1;
+  GparityWilsonFermionR::ImplParams params;  params.twists = twists;
+  GparityWilsonFermionR Wil(U,*UGrid,*UrbGrid,mass,params);
+  Wil.M   (phi,Mphi);
+
+  ComplexD S    = innerProduct(Mphi,Mphi); // pdag MdagM p
+
+  // get the deriv of phidag MdagM phi with respect to "U"
+  LatticeGaugeField UdSdU(UGrid);
+  LatticeGaugeField tmp(UGrid);
+
+  Wil.MDeriv(tmp , Mphi,  phi,DaggerNo );  UdSdU=tmp;
+  Wil.MDeriv(tmp , phi,  Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);  
+  
+  FermionField Ftmp      (FGrid);
+
+  ////////////////////////////////////
+  // Modify the gauge field a little 
+  ////////////////////////////////////
+  RealD dt = 0.00001;
+
+  LatticeColourMatrix mommu(UGrid); 
+  LatticeColourMatrix forcemu(UGrid); 
+  LatticeGaugeField mom(UGrid); 
+  LatticeGaugeField Uprime(UGrid); 
+
+  for(int mu=0;mu<Nd;mu++){
+
+    SU3::GaussianLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+
+    PokeIndex<LorentzIndex>(mom,mommu,mu);
+
+    // fourth order exponential approx
+    parallel_for(auto i=mom.begin();i<mom.end();i++){
+      Uprime[i](mu) =
+	  U[i](mu)
+	+ mom[i](mu)*U[i](mu)*dt 
+	+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
+	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
+	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
+	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
+	+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
+	;
+    }
+
+  }
+  
+  Wil.ImportGauge(Uprime);
+  Wil.M          (phi,MphiPrime);
+
+  ComplexD Sprime    = innerProduct(MphiPrime   ,MphiPrime);
+
+  //////////////////////////////////////////////
+  // Use derivative to estimate dS
+  //////////////////////////////////////////////
+
+
+  LatticeComplex dS(UGrid); dS = zero;
+  for(int mu=0;mu<Nd;mu++){
+    mommu   = PeekIndex<LorentzIndex>(UdSdU,mu);
+    mommu=Ta(mommu)*2.0;
+    PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
+  }
+
+  for(int mu=0;mu<Nd;mu++){
+    forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
+    mommu   = PeekIndex<LorentzIndex>(mom,mu);
+
+    // Update PF action density
+    dS = dS+trace(mommu*forcemu)*dt;
+
+  }
+
+  Complex dSpred    = sum(dS);
+
+  // From TwoFlavourPseudoFermion:
+  //////////////////////////////////////////////////////
+  // dS/du = - phi^dag  (Mdag M)^-1 [ Mdag dM + dMdag M ]  (Mdag M)^-1 phi
+  //       = - phi^dag M^-1 dM (MdagM)^-1 phi -  phi^dag (MdagM)^-1 dMdag dM (Mdag)^-1 phi 
+  //
+  //       = - Ydag dM X  - Xdag dMdag Y
+  //
+  //////////////////////////////////////////////////////
+  // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
+  //
+  // 
+
+  std::cout << GridLogMessage << " S      "<<S<<std::endl;
+  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
+  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+  std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
+
+  std::cout<< GridLogMessage << "Done" <<std::endl;
+  Grid_finalize();
+}