diff --git a/Grid/Makefile.am b/Grid/Makefile.am index 8472dd71a2..0148ba3e15 100644 --- a/Grid/Makefile.am +++ b/Grid/Makefile.am @@ -56,6 +56,8 @@ include Eigen.inc extra_sources+=$(WILS_FERMION_FILES) extra_sources+=$(STAG_FERMION_FILES) +extra_sources+=$(XCONJ_FERMION_FILES) + if BUILD_ZMOBIUS extra_sources+=$(ZWILS_FERMION_FILES) endif diff --git a/Grid/qcd/action/ActionParams.h b/Grid/qcd/action/ActionParams.h index 122dfb9cba..3be5ae0ca5 100644 --- a/Grid/qcd/action/ActionParams.h +++ b/Grid/qcd/action/ActionParams.h @@ -46,6 +46,18 @@ struct GparityWilsonImplParams { }; }; +struct XconjWilsonImplParams { + Coordinate dirichlet; // Blocksize of dirichlet BCs + int partialDirichlet; + + Coordinate twists; //Here the first Nd-1 directions are treated as "spatial", and a twist value of 1 indicates G-parity BCs in that direction. + //mu=Nd-1 is assumed to be the time direction and a twist value of 1 indicates antiperiodic BCs + ComplexD boundary_phase; //+1 for X-conjugate, -1 for Xbar-conjugate, or other + XconjWilsonImplParams() : twists(Nd, 0), boundary_phase(1.0) { dirichlet.resize(0); partialDirichlet=0; }; +}; + + + struct WilsonImplParams { bool overlapCommsCompute; Coordinate dirichlet; // Blocksize of dirichlet BCs diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index 8d7b3fdc30..0a37e93501 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -271,6 +271,19 @@ typedef NaiveStaggeredFermion NaiveStaggeredFermionD; typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermion5DF; typedef ImprovedStaggeredFermion5D ImprovedStaggeredFermion5DD; + +typedef DomainWallFermion XconjugateDomainWallFermionD2; +typedef DomainWallFermion XconjugateDomainWallFermionF; +typedef DomainWallFermion XconjugateDomainWallFermionD; + +typedef MobiusFermion XconjugateMobiusFermionD2; +typedef MobiusFermion XconjugateMobiusFermionF; +typedef MobiusFermion XconjugateMobiusFermionD; + +typedef MobiusEOFAFermion XconjugateMobiusEOFAFermionD2; +typedef MobiusEOFAFermion XconjugateMobiusEOFAFermionF; +typedef MobiusEOFAFermion XconjugateMobiusEOFAFermionD; + NAMESPACE_END(Grid); //////////////////// diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index 56aaca12cf..8dbb292763 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -168,6 +168,10 @@ NAMESPACE_CHECK(ImplBase); #include NAMESPACE_CHECK(ImplWilson); +#include +NAMESPACE_CHECK(ImplXconj); + + //////////////////////////////////////////////////////////////////////////////////////// // Flavour doubled spinors; is Gparity the only? what about C*? //////////////////////////////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/WilsonCompressor.h b/Grid/qcd/action/fermion/WilsonCompressor.h index 186fa2786f..a129037492 100644 --- a/Grid/qcd/action/fermion/WilsonCompressor.h +++ b/Grid/qcd/action/fermion/WilsonCompressor.h @@ -391,6 +391,23 @@ class WilsonProjector { default: assert(0); break; } } + + template + static accelerator void Recon(fsp &result,const hsp &in,int mu,int dag){ + int mudag=dag? mu : (mu+Nd)%(2*Nd); + switch(mudag) { + case Xp: spReconXp(result,in); break; + case Yp: spReconYp(result,in); break; + case Zp: spReconZp(result,in); break; + case Tp: spReconTp(result,in); break; + case Xm: spReconXm(result,in); break; + case Ym: spReconYm(result,in); break; + case Zm: spReconZm(result,in); break; + case Tm: spReconTm(result,in); break; + default: assert(0); break; + } + } + }; template using WilsonCompressor = WilsonCompressorTemplate; diff --git a/Grid/qcd/action/fermion/XconjImpl.h b/Grid/qcd/action/fermion/XconjImpl.h new file mode 100644 index 0000000000..032beb3123 --- /dev/null +++ b/Grid/qcd/action/fermion/XconjImpl.h @@ -0,0 +1,276 @@ +#pragma once + +NAMESPACE_BEGIN(Grid); + + +///////////////////////////////////////////////////////////////////////////// +// Single flavour four spinors with colour index +///////////////////////////////////////////////////////////////////////////// +template +class XconjugateWilsonImpl : public ConjugateGaugeImpl > { +public: + + static const int Dimension = Representation::Dimension; + static const bool isFundamental = Representation::isFundamental; + static const bool LsVectorised=false; + static const bool isGparity=false; + static const int Nhcs = Options::Nhcs; + + typedef ConjugateGaugeImpl< GaugeImplTypes > Gimpl; + INHERIT_GIMPL_TYPES(Gimpl); + + //Necessary? + constexpr bool is_fundamental() const{return Dimension == Nc ? 1 : 0;} + + typedef typename Options::_Coeff_t Coeff_t; + typedef typename Options::template PrecisionMapper::LowerPrecVector SimdL; + + template using iImplSpinor = iScalar, Ns> >; + template using iImplPropagator = iScalar, Ns> >; + template using iImplHalfSpinor = iScalar, Nhs> >; + template using iImplHalfCommSpinor = iScalar, Nhcs> >; + template using iImplDoubledGaugeField = iVector >, Nds>; + + typedef iImplSpinor SiteSpinor; + typedef iImplPropagator SitePropagator; + typedef iImplHalfSpinor SiteHalfSpinor; + typedef iImplHalfCommSpinor SiteHalfCommSpinor; + typedef iImplDoubledGaugeField SiteDoubledGaugeField; + + typedef Lattice FermionField; + typedef Lattice PropagatorField; + typedef Lattice DoubledGaugeField; + + typedef WilsonCompressor Compressor; + typedef XconjWilsonImplParams ImplParams; + typedef WilsonStencil StencilImpl; + typedef const typename StencilImpl::View_type StencilView; + + ImplParams Params; + + XconjugateWilsonImpl(const ImplParams &p = ImplParams()) : Params(p){ + }; + + template + static accelerator_inline void multLink(_Spinor &phi, + const SiteDoubledGaugeField &U, + const _Spinor &chi, + int mu) + { + assert(0); + } + + + template + static accelerator_inline void multLink(_Spinor &phi, + const SiteDoubledGaugeField &U, + const _Spinor &chi, + int mu, + StencilEntry *SE, + StencilView &St) + { + Gamma X = Gamma(Gamma::Algebra::MinusSigmaXZ); + typedef decltype(coalescedRead( *( (SiteSpinor*)NULL ) )) _FullSpinor; + + int direction = St._directions[mu]; + int distance = St._distances[mu]; + int ptype = St._permute_type[mu]; + int sl = St._simd_layout[direction]; + Coordinate icoor; + + const int Nsimd =SiteDoubledGaugeField::Nsimd(); + int s = acceleratorSIMTlane(Nsimd); + St.iCoorFromIindex(icoor,s); + + int mmu = mu % Nd; + int dag = mu / Nd; + + auto UU=coalescedRead(U(mu)); + + //Decide whether we do a G-parity flavor twist + //Note: this assumes (but does not check) that sl==1 || sl==2 i.e. max 2 SIMD lanes in G-parity dir + //It also assumes (but does not check) that abs(distance) == 1 + int permute_lane = (sl==1) + || ((distance== 1)&&(icoor[direction]==1)) + || ((distance==-1)&&(icoor[direction]==0)); + + permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu] && mmu < Nd-1; //only if we are going around the world in a spatial direction + + //Implement the X-conjugate BC + /* + Note + C g^mu C = g^mu,T + g^mu C = -C g^mu,T + + We need + (1+g^mu) Cg^5 psi^* + = g^5(1-g^mu) C psi^* + = g^5 (C - g^mu C) psi^* + = g^5 (C + C g^mu,T ) psi^* + = g^5 C (1 + g^mu,T ) psi^* + = [ g^5 C (1 + g^mu,dag ) psi ]^* + = [ g^5 C (1 + g^mu ) psi ]^* + + We also need to return a half spinor but we can apply the projection again once we have + (1+g^mu) Cg^5 psi^* + + Note: + Applying Proj after Recon picks up an extra factor of 2 + e.g. (from TwoSpinor.h) + spProjXp = [0]+i[3] + [1]+i[2] + + spReconXp = [0] + [1] + -i[1] + -i[0] + + spProjXp spReconXp = [0] + [0] + [1] + [1] + */ + if(permute_lane){ + _FullSpinor fs; + WilsonProjector::Recon(fs,chi,mmu,dag); + fs = distance == 1 ? -0.5*(X*conjugate(fs)) : 0.5*(X*conjugate(fs)); //extra 0.5 due to Recon/Proj cycle + fs = fs * St.parameters.boundary_phase; + _Spinor XchiStar; + WilsonProjector::Proj(XchiStar,fs,mmu,dag); + mult(&phi(),&UU,&XchiStar()); + }else{ + mult(&phi(),&UU,&chi()); + } + } + + template + inline void multLinkField(_SpinorField & out, + const DoubledGaugeField &Umu, + const _SpinorField & phi, + int mu) + { + assert(0); + } + + template + static accelerator_inline void loadLinkElement(Simd ®, ref &memory) + { + reg = memory; + } + + inline void DoubleStore(GridBase *GaugeGrid, + DoubledGaugeField &Uds, + const GaugeField &Umu) + { + typedef typename Simd::scalar_type scalar_type; + + conformable(Uds.Grid(), GaugeGrid); + conformable(Umu.Grid(), GaugeGrid); + + GaugeLinkField U(GaugeGrid); + GaugeLinkField tmp(GaugeGrid); + + Lattice > coor(GaugeGrid); + for (int mu = 0; mu <= Nd-1; mu++) { //mu=Nd-1 is either periodic or antiperiodic + int L = GaugeGrid->GlobalDimensions()[mu]; + int Lmu = L - 1; + + LatticeCoordinate(coor, mu); + + U = PeekIndex(Umu, mu); + + if(mu == Nd-1 && Params.twists[mu] ){ //antiperiodic BCs + tmp = where(coor == Lmu, -U, U); + PokeIndex(Uds, tmp, mu); + + U = adj(Cshift(U, mu, -1)); + U = where(coor == 0, -U, U); + PokeIndex(Uds, U, mu + 4); + }else if(mu < Nd-1 && Params.twists[mu] ){ //charge conj BCs + PokeIndex(Uds, U, mu); + + U = adj(Cshift(U, mu, -1)); + tmp = conjugate(U); + U = where(coor == 0, tmp, U); + + PokeIndex(Uds, U, mu + 4); + }else{ //periodic BCs + PokeIndex(Uds, U, mu); + + U = adj(Cshift(U, mu, -1)); + PokeIndex(Uds, U, mu + 4); + } + } + } + + inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A,int mu){ + GaugeLinkField link(mat.Grid()); + link = TraceIndex(outerProduct(Btilde,A)); + PokeIndex(mat,link,mu); + } + + inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){ + mat = outerProduct(B,A); + } + + inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) { + mat = TraceIndex(P); + } + + inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds) + { + for (int mu = 0; mu < Nd; mu++) + mat[mu] = PeekIndex(Uds, mu); + } + + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã,int mu) + { +#undef USE_OLD_INSERT_FORCE + int Ls=Btilde.Grid()->_fdimensions[0]; + autoView( mat_v , mat, AcceleratorWrite); +#ifdef USE_OLD_INSERT_FORCE + GaugeLinkField tmp(mat.Grid()); + tmp = Zero(); + { + const int Nsimd = SiteSpinor::Nsimd(); + autoView( tmp_v , tmp, AcceleratorWrite); + autoView( Btilde_v , Btilde, AcceleratorRead); + autoView( Atilde_v , Atilde, AcceleratorRead); + accelerator_for(sss,tmp.Grid()->oSites(),1,{ + int sU=sss; + for(int s=0;s(outerProduct(Btilde_v[sF],Atilde_v[sF])); // ordering here + } + }); + } + PokeIndex(mat,tmp,mu); +#else + { + const int Nsimd = SiteSpinor::Nsimd(); + autoView( Btilde_v , Btilde, AcceleratorRead); + autoView( Atilde_v , Atilde, AcceleratorRead); + accelerator_for(sss,mat.Grid()->oSites(),Nsimd,{ + int sU=sss; + typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType; + ColorMatrixType sum; + zeroit(sum); + for(int s=0;s XconjugateWilsonImplR; // Real.. whichever prec +typedef XconjugateWilsonImpl XconjugateWilsonImplF; // Float +typedef XconjugateWilsonImpl XconjugateWilsonImplD; // Double + +NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/CayleyFermion5DInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/CayleyFermion5DInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..cb1db62588 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/CayleyFermion5DInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../CayleyFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/ContinuedFractionFermion5DInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/ContinuedFractionFermion5DInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..c2d4b8fccc --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/ContinuedFractionFermion5DInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../ContinuedFractionFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/DomainWallEOFAFermionInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/DomainWallEOFAFermionInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..2f550a2b19 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/DomainWallEOFAFermionInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../DomainWallEOFAFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/MobiusEOFAFermionInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/MobiusEOFAFermionInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..7a8f1172e7 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/MobiusEOFAFermionInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../MobiusEOFAFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/PartialFractionFermion5DInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/PartialFractionFermion5DInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..7f4cea7185 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/PartialFractionFermion5DInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../PartialFractionFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/WilsonFermion5DInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/WilsonFermion5DInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..804d088449 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/WilsonFermion5DInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/WilsonKernelsInstantiationXconjugateWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/WilsonKernelsInstantiationXconjugateWilsonImplD.cc new file mode 120000 index 0000000000..01c35e7bc1 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/WilsonKernelsInstantiationXconjugateWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/impl.h b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/impl.h new file mode 100644 index 0000000000..5ee8ad4d04 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplD/impl.h @@ -0,0 +1 @@ +#define IMPLEMENTATION XconjugateWilsonImplD diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/CayleyFermion5DInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/CayleyFermion5DInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..cb1db62588 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/CayleyFermion5DInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../CayleyFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/ContinuedFractionFermion5DInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/ContinuedFractionFermion5DInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..c2d4b8fccc --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/ContinuedFractionFermion5DInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../ContinuedFractionFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/DomainWallEOFAFermionInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/DomainWallEOFAFermionInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..2f550a2b19 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/DomainWallEOFAFermionInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../DomainWallEOFAFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/MobiusEOFAFermionInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/MobiusEOFAFermionInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..7a8f1172e7 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/MobiusEOFAFermionInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../MobiusEOFAFermionInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/PartialFractionFermion5DInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/PartialFractionFermion5DInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..7f4cea7185 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/PartialFractionFermion5DInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../PartialFractionFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/WilsonFermion5DInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/WilsonFermion5DInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..804d088449 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/WilsonFermion5DInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonFermion5DInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/WilsonKernelsInstantiationXconjugateWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/WilsonKernelsInstantiationXconjugateWilsonImplF.cc new file mode 120000 index 0000000000..01c35e7bc1 --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/WilsonKernelsInstantiationXconjugateWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/impl.h b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/impl.h new file mode 100644 index 0000000000..d50b56e15e --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/XconjugateWilsonImplF/impl.h @@ -0,0 +1 @@ +#define IMPLEMENTATION XconjugateWilsonImplF diff --git a/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh b/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh index 728dc5e783..4cc3a3fe56 100755 --- a/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh +++ b/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh @@ -34,7 +34,10 @@ DWF_IMPL_LIST=" \ WilsonImplD \ WilsonImplD2 \ ZWilsonImplF \ - ZWilsonImplD2 " + ZWilsonImplD \ + ZWilsonImplD2 \ + XconjugateWilsonImplF \ + XconjugateWilsonImplD " GDWF_IMPL_LIST=" \ GparityWilsonImplF \ diff --git a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h index 9d4df1d35e..343e050d04 100644 --- a/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h +++ b/Grid/qcd/action/pseudofermion/ExactOneFlavourRatio.h @@ -126,6 +126,11 @@ NAMESPACE_BEGIN(Grid); const FermionField &getPhi() const{ return Phi; } + //Set the pseudofermion Phi field for testing or other purposes + void setPhi(const FermionField &phi_in){ + Phi = phi_in; + } + virtual std::string action_name() { return "ExactOneFlavourRatioPseudoFermionAction"; } virtual std::string LogParameters() { diff --git a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h index f237dee44a..a9a40c1ab3 100644 --- a/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h +++ b/Grid/qcd/action/pseudofermion/GeneralEvenOddRationalRatio.h @@ -226,6 +226,13 @@ NAMESPACE_BEGIN(Grid); //Access the fermion field const FermionField &getPhiOdd() const{ return PhiOdd; } + //Set the the pseudofermion Phi field for testing or other purposes + //Note: this action ignores the contribution from PhiEven + void setPhi(const FermionField &phi_in){ + pickCheckerboard(Even,PhiEven,phi_in); + pickCheckerboard(Odd,PhiOdd,phi_in); + } + virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) { std::cout< + Author: Peter Boyle + Author: paboyle + + 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 +#include +using namespace std; +using namespace Grid; + +typedef typename XconjugateDomainWallFermionF::FermionField LatticeFermionF; +typedef typename XconjugateDomainWallFermionD::FermionField LatticeFermionD; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); +#ifdef ENABLE_GPARITY + int Ls=16; + bool do_fp32 = true; + bool do_fp64 = false; + + for(int i=0;i> Ls; + }else if(std::string(argv[i]) == "-precision"){ + std::string p = argv[i+1]; + if(p == "both"){ do_fp32 = true; do_fp64 = true; } + else if(p == "single"){ do_fp32 = true; do_fp64 = false; } //default + else if(p == "double"){ do_fp32 = false; do_fp64 = true; } + else{ + assert(0 && "Invalid precision argument"); + } + } + } + + int threads = GridThread::GetThreads(); + std::cout<::HotConfiguration(RNG4,Umu); + std::cout << GridLogMessage << "Random gauge initialised " << std::endl; + + RealD mass=0.1; + RealD M5 =1.8; + + RealD NP = UGrid->_Nprocessors; + RealD NN = UGrid->NodeCount(); + + std::cout << GridLogMessage<< "*****************************************************************" < twists({1,1,1,0}); + XconjugateDomainWallFermionF::ImplParams xparams; + xparams.twists = twists; + xparams.boundary_phase = 1.0; + + int ncall =1000; + + if(do_fp32){ + std::cout << GridLogMessage<< "* SINGLE/SINGLE"<Barrier(); + Dw.Dhop(src,result,0); + std::cout<Barrier(); + + double volume=Ls; for(int mu=0;muBarrier(); + DwD.Dhop(src_d,result_d,0); + std::cout<Barrier(); + + double volume=Ls; for(int mu=0;mu + +using namespace std; +using namespace Grid; + +//A test implementation of the X-conjugate action as a wrapper around the regular 2f G-parity action +template +class XconjugateDWF{ +public: + typedef typename GPaction::FermionField GPfermion; + typedef LatticeFermionD FermionField; + GPaction *action; + bool Xbar; + XconjugateDWF(GPaction *action, bool Xbar = false): action(action), Xbar(Xbar){} + + template + LatticeFermionD op11(const LatticeFermionD &in, const Op &op){ + GPfermion tmp1(in.Grid()), tmp2(in.Grid()); + tmp1.Checkerboard() = in.Checkerboard(); + tmp1 = Zero(); + PokeIndex(tmp1, in, 0); + op(tmp2, tmp1); + return PeekIndex(tmp2,0); + } + template + LatticeFermionD op12(const LatticeFermionD &in, const Op &op){ + GPfermion tmp1(in.Grid()), tmp2(in.Grid()); + tmp1.Checkerboard() = in.Checkerboard(); + tmp1 = Zero(); + PokeIndex(tmp1, in, 1); + op(tmp2, tmp1); + return PeekIndex(tmp2,0); + } + + template + LatticeFermionD opFull(const LatticeFermionD &in, const Op &op){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + LatticeFermionD tmp1(in.Grid()); + tmp1.Checkerboard() = in.Checkerboard(); + tmp1 = -(X*conjugate(in)); + if(Xbar) tmp1 = -tmp1; + + LatticeFermionD out11 = op11(in, op); + LatticeFermionD out12 = op12(tmp1, op); + LatticeFermionD out = out11 + out12; + return out; + } + +#define DEFOP(OP) void OP(const FermionField &in, FermionField &out){ out=opFull(in, [&](GPfermion &out, const GPfermion &in){ return action->OP(in,out); }); } + DEFOP(M); + DEFOP(Mdag); + DEFOP(Meooe); + DEFOP(MeooeDag); + DEFOP(Mooee); + DEFOP(MooeeDag); + DEFOP(MooeeInv); + DEFOP(MooeeInvDag); +#undef DEFOP; +}; + +//A 2-flavor representation of the X-conjugate action acting on X-conjugate fermion fields +template +class Xconjugate2fWrapper{ +public: + typedef typename XconjAction::FermionField OneFlavorFermionField; + typedef Lattice > FermionField; + + XconjAction *action; + GPAction *gaction; + + Xconjugate2fWrapper(XconjAction *action, GPAction *gaction): action(action), gaction(gaction){} + + template + void opFull(FermionField &out, const FermionField &in, const Op &op){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + + OneFlavorFermionField tmp(in.Grid()); + op(tmp, PeekIndex(in,0)); + + out.Checkerboard() = tmp.Checkerboard(); + PokeIndex(out, tmp, 0); + tmp = -(X*conjugate(tmp)); + tmp = tmp * action->Params.boundary_phase; + PokeIndex(out, tmp, 1); + } + +#define DEFOP(OP) void OP(const FermionField &in, FermionField &out){ opFull(out, in, \ + [&](OneFlavorFermionField &out1f, \ + const OneFlavorFermionField &in1f){ \ + return action->OP(in1f,out1f); \ + } \ + ); \ + } + DEFOP(M); + DEFOP(Mdag); + DEFOP(Meooe); + DEFOP(MeooeDag); + DEFOP(Mooee); + DEFOP(MooeeDag); + DEFOP(MooeeInv); + DEFOP(MooeeInvDag); +#undef DEFOP; +}; + +const Gamma & Xmatrix(){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + return X; +} + +typedef typename GparityMobiusFermionD::FermionField FermionField2f; +typedef typename XconjugateMobiusFermionD::FermionField FermionField1f; + + +void boostXconjToGparity(FermionField2f &out, const FermionField1f &in){ + PokeIndex(out, in, 0); + FermionField1f tmp = -(Xmatrix()*conjugate(in)); + PokeIndex(out, tmp, 1); +} +void boostXbarConjToGparity(FermionField2f &out, const FermionField1f &in){ + PokeIndex(out, in, 0); + FermionField1f tmp = Xmatrix()*conjugate(in); + PokeIndex(out, tmp, 1); +} + +template +inline RealD norm2Diff(const Field &a, const Field &b){ + return norm2(Field(a-b)); +} + +void norm2DiffXconj(double &diff_f0, double &diff_f1, const FermionField2f &gp, const FermionField1f &xconj){ + FermionField1f tmp = PeekIndex(gp,0); + diff_f0 = norm2(FermionField1f(tmp-xconj)); + tmp = PeekIndex(gp,1); + FermionField1f xconj_f1 = -(Xmatrix()*conjugate(xconj)); + diff_f1 = norm2(FermionField1f(tmp-xconj_f1)); +} +void norm2DiffXbarConj(double &diff_f0, double &diff_f1, const FermionField2f &gp, const FermionField1f &xconj){ + FermionField1f tmp = PeekIndex(gp,0); + diff_f0 = norm2(FermionField1f(tmp-xconj)); + tmp = PeekIndex(gp,1); + FermionField1f xconj_f1 = Xmatrix()*conjugate(xconj); + diff_f1 = norm2(FermionField1f(tmp-xconj_f1)); +} + +//1/sqrt(2) * +//| X -1 | +//| -i iX | +void applyUdag(FermionField2f &out, const FermionField2f &in){ + FermionField1f u = PeekIndex(in,0); + FermionField1f l = PeekIndex(in,1); + FermionField1f tmp = ComplexD(1./sqrt(2),0) * (Xmatrix()*u - l); + PokeIndex(out, tmp, 0); + tmp = ComplexD(0,-1./sqrt(2))*(u - Xmatrix()*l); + PokeIndex(out, tmp, 1); +} +//1/sqrt(2) * +//| -X i | +//| -1 iX | +void applyU(FermionField2f &out, const FermionField2f &in){ + FermionField1f u = PeekIndex(in,0); + FermionField1f l = PeekIndex(in,1); + FermionField1f tmp = ComplexD(-1./sqrt(2)) * (Xmatrix()*u + ComplexD(0,-1)*l); + PokeIndex(out, tmp, 0); + tmp = ComplexD(-1./sqrt(2))*(u + ComplexD(0,-1)*(Xmatrix()*l)); + PokeIndex(out, tmp, 1); +} +template +void applyR(FermionField2f &out, const FermionField2f &in, GPAction &action){ + FermionField2f tmp(in.Grid()), tmp2(in.Grid()); + applyU(tmp, in); + action.M(tmp, tmp2); + applyUdag(out,tmp2); //Rv +} +template +void applyRij(FermionField1f &out, const FermionField1f &in, const int i, const int j, GPAction &action){ + FermionField2f tmp2f(in.Grid()); + tmp2f = Zero(); + PokeIndex(tmp2f, in, j); + FermionField2f tmp2f_2(in.Grid()); + applyR(tmp2f_2, tmp2f, action); + out = PeekIndex(tmp2f_2,i); +} +template +void applyMij(FermionField1f &out, const FermionField1f &in, const int i, const int j, GPAction &action){ + FermionField2f tmp2f(in.Grid()); + tmp2f = Zero(); + PokeIndex(tmp2f, in, j); + FermionField2f tmp2f_2(in.Grid()); + action.M(tmp2f,tmp2f_2); + out = PeekIndex(tmp2f_2,i); +} +template +void applyMijstar(FermionField1f &out, const FermionField1f &in, const int i, const int j, GPAction &action){ + FermionField1f tmp(in.Grid()), tmp2(in.Grid()); //Mij* v = (Mij v*)* + tmp = conjugate(in); + applyMij(tmp2,tmp,i,j,action); + out = conjugate(tmp2); +} + + +//A reference implementation of the real, rotated GP Dirac operator +template +class RotatedGPrefAction{ +public: + typedef typename GPAction::FermionField TwoFlavorFermionField; + typedef Lattice > OneFlavorFermionField; + + GPAction *gaction; + + RotatedGPrefAction(GPAction *gaction): gaction(gaction){} + + //out = D_ij in + template + void DD(OneFlavorFermionField &out, const OneFlavorFermionField &in, const int i, const int j, const Op &op){ + TwoFlavorFermionField tmp2f(in.Grid()); + tmp2f = Zero(); + PokeIndex(tmp2f, in, j); + TwoFlavorFermionField tmp2f_2(in.Grid()); + op(tmp2f_2, tmp2f); + out = PeekIndex(tmp2f_2,i); + } + //out = D*_ij in + template + void DDstar(OneFlavorFermionField &out, const OneFlavorFermionField &in, const int i, const int j, const Op &op){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + tmp = conjugate(in); + DD(tmp2,tmp,i,j,op); + out = conjugate(tmp2); + } + + //Use Re(U) V = 0.5 ( U + U^*) V = 0.5 ( UV + [U V*]* ) + template + void do_real(OneFlavorFermionField &out, const OneFlavorFermionField &in, const U &u){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + u(out,in); + tmp = conjugate(in); + u(tmp2,tmp); + tmp = conjugate(tmp2); + out = out + tmp; + out = out * 0.5; + } + //Use Im(U) V = -0.5i ( U - U^*) V = -0.5i ( UV - [U V*]* ) + template + void do_imag(OneFlavorFermionField &out, const OneFlavorFermionField &in, const U &u){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + u(out,in); + + tmp = conjugate(in); + u(tmp2,tmp); + tmp = conjugate(tmp2); + out = out - tmp; + out = out * ComplexD(0,-0.5); + } + + //-(X D_11 X + X D_12) in + template + void opAparen(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + tmp = Xmatrix()*in; + DD(tmp2,tmp,0,0,op); + out = Xmatrix()*tmp2; + + DD(tmp,in,0,1,op); + tmp2 = Xmatrix() * tmp; + out = out + tmp2; + out = -out; + } + template + void opA(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + do_real(out,in, [&](OneFlavorFermionField &out, const OneFlavorFermionField &in){ opAparen(out,in,op); } ); + } + + + //(-X D11 -X D12 X) in + template + void opBparen(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + DD(tmp,in,0,0,op); + out = Xmatrix()*tmp; + + tmp = Xmatrix()*in; + DD(tmp2,tmp,0,1,op); + tmp = Xmatrix() * tmp2; + out = out + tmp; + + out = -out; + } + template + void opB(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + do_imag(out,in, [&](OneFlavorFermionField &out, const OneFlavorFermionField &in){ opBparen(out,in,op); } ); + } + + //-(D11 X + D12 ) in + template + void opCparen(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + tmp = Xmatrix()*in; + DD(out,tmp,0,0,op); + + DD(tmp,in,0,1,op); + out = out + tmp; + + out = -out; + } + template + void opC(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + do_imag(out,in, [&](OneFlavorFermionField &out, const OneFlavorFermionField &in){ opCparen(out,in,op); } ); + } + + //(D11 + D12 X ) in + template + void opDparen(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + OneFlavorFermionField tmp(in.Grid()), tmp2(in.Grid()); + DD(out,in,0,0,op); + + tmp = Xmatrix()*in; + DD(tmp2,tmp,0,1,op); + out = out + tmp2; + } + template + void opD(OneFlavorFermionField &out, const OneFlavorFermionField &in, const Op &op){ + do_real(out,in, [&](OneFlavorFermionField &out, const OneFlavorFermionField &in){ opDparen(out,in,op); } ); + } + + template + void opFull(TwoFlavorFermionField &out, const TwoFlavorFermionField &in, const Op &op){ + OneFlavorFermionField u(in.Grid()), l(in.Grid()), tmp(in.Grid()), tmp2(in.Grid()); + u = PeekIndex(in,0); + l = PeekIndex(in,1); + opA(tmp,u,op); + opB(tmp2,l,op); + tmp = tmp + tmp2; + PokeIndex(out, tmp, 0); + + opC(tmp,u,op); + opD(tmp2,l,op); + tmp = tmp + tmp2; + PokeIndex(out, tmp, 1); + } + +#define DEFOP(OP) void OP(const TwoFlavorFermionField &in, TwoFlavorFermionField &out){ opFull(out, in, \ + [&](TwoFlavorFermionField &out, \ + const TwoFlavorFermionField &in){ \ + return gaction->OP(in,out); \ + } \ + ); \ + } + DEFOP(M); + DEFOP(Mdag); + DEFOP(Meooe); + DEFOP(MeooeDag); + DEFOP(Mooee); + DEFOP(MooeeDag); + DEFOP(MooeeInv); + DEFOP(MooeeInvDag); +#undef DEFOP; +}; + + +template +void applyRijRef(FermionField1f &out, const FermionField1f &in, const int i, const int j, RefAction &action){ + FermionField2f tmp2f(in.Grid()); + tmp2f = Zero(); + PokeIndex(tmp2f, in, j); + FermionField2f tmp2f_2(in.Grid()); + action.M(tmp2f,tmp2f_2); + out = PeekIndex(tmp2f_2,i); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeGaugeField Umu(UGrid); + SU::HotConfiguration(RNG4, Umu); + + Gamma X = Xmatrix(); + + //Set up a regular MDWF action instance as well as X-conj and Xbar-conj versions + RealD mass=0.01; + RealD M5=1.8; + RealD mob_b=1.5; + GparityMobiusFermionD ::ImplParams params; + std::vector twists({1,1,1,0}); + params.twists = twists; + + GparityMobiusFermionD reg_action(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); + + XconjugateMobiusFermionD::ImplParams xparams; + xparams.twists = twists; + xparams.boundary_phase = 1.0; + + XconjugateMobiusFermionD xconj_action(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,xparams); + + xparams.boundary_phase = -1.0; + XconjugateMobiusFermionD xbarconj_action(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,xparams); + + + //####################################################################################################################################### + { + std::cout << "Testing M on full grid" << std::endl; + + FermionField1f rand_sc(FGrid), rand_sc2(FGrid); + gaussian(RNG5,rand_sc); + gaussian(RNG5,rand_sc2); + + FermionField2f rand_f0(FGrid); //v \delta_f,0 + rand_f0 = Zero(); + PokeIndex(rand_f0, rand_sc, 0); + + FermionField2f rand_f1(FGrid); //v \delta_f,1 + rand_f1 = Zero(); + PokeIndex(rand_f1, rand_sc, 1); + + FermionField2f tmp(FGrid), tmp2(FGrid), tmp3(FGrid), tmp4(FGrid), tmp5(FGrid); + FermionField1f tmpsc(FGrid), tmpsc2(FGrid), tmpsc3(FGrid), tmpsc4(FGrid); + RealD nrm; + + std::cout << "Test the relationship between the upper and lower rows in flavor space of the G-parity Dirac operator" << std::endl; + + //Test M00 v = -X M11* X v = -X [M11 X v*]* + reg_action.M(rand_f0,tmp); + FermionField1f M00v = PeekIndex(tmp,0); + FermionField1f M10v = PeekIndex(tmp,1); + + tmp = X*conjugate(rand_f1); + reg_action.M(tmp,tmp2); + FermionField1f M11Xvconj = PeekIndex(tmp2,1); + FermionField1f M01Xvconj = PeekIndex(tmp2,0); + + tmpsc = -(X*conjugate(M11Xvconj)); + nrm = norm2Diff(tmpsc, M00v); + std::cout << "Test of M00 v = -X M11* X v = -X [M11 X v*]* (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test M10 v = X M01* X v = X [M01 X v*]* + tmpsc = X*conjugate(M01Xvconj); + nrm = norm2Diff(tmpsc, M10v); + std::cout << "Test of M10 v = X M11* X v = X [M11 X v*]* (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test the X-conjugate implementation + std::cout << "Test the implementation of the X-conjugate action against the reference"<< std::endl; + XconjugateDWF xconj_action_ref(®_action); + XconjugateDWF xbarconj_action_ref(®_action,true); + + //Test upper boundary + std::vector L(4); + for(int mu=0;mu<4;mu++) + L[mu] = UGrid->GlobalDimensions()[mu]; + Coordinate site(5,0); + typedef FermionField1f::vector_object SpinorV; + typedef typename SpinorV::scalar_object SpinorS; + + tmpsc3 = Zero(); + for(int i=2;i<5;i++) site[i] = L[i-1]/2; //midpoint in y,z,t + site[1] = 0; //field only on site on lower boundary + + SpinorS v; + peekSite(v, rand_sc, site); + pokeSite(v, tmpsc3, site); + + xconj_action_ref.M(tmpsc3, tmpsc); + xconj_action.M(tmpsc3, tmpsc2); + nrm = norm2Diff(tmpsc, tmpsc2); + std::cout << "Test of Xconjugate action M upper boundary only (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + site[1] = L[0]-1; + + SpinorS r1,r2; + peekSite(r1, tmpsc, site); + peekSite(r2, tmpsc2, site); + nrm = norm2Diff(r1,r2); + std::cout << "Results L-1\nref: " << r1 << "\nimpl: " << r2 << "\ndiff (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + site[1] = 1; + peekSite(r1, tmpsc, site); + peekSite(r2, tmpsc2, site); + nrm = norm2Diff(r1,r2); + std::cout << "Results 1\nref: " << r1 << "\nimpl: " << r2 << "\ndiff (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test lower boundary + site[1] = L[0]-1; + tmpsc3 = Zero(); + pokeSite(v, tmpsc3, site); + + + xconj_action_ref.M(tmpsc3, tmpsc); + xconj_action.M(tmpsc3, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of Xconjugate action M lower boundary only (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + site[1] = 0; + + peekSite(r1, tmpsc, site); + peekSite(r2, tmpsc2, site); + nrm = norm2Diff(r1,r2); + std::cout << "Results 0\nref: " << r1 << "\nimpl: " << r2 << "\ndiff (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + site[1] = L[0]-2; + + peekSite(r1, tmpsc, site); + peekSite(r2, tmpsc2, site); + nrm = norm2Diff(r1,r2); + std::cout << "Results L-2\nref: " << r1 << "\nimpl: " << r2 << "\ndiff (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test full + xconj_action_ref.M(rand_sc, tmpsc); + xconj_action.M(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of Xconjugate action M against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xconj_action_ref.Mdag(rand_sc, tmpsc); + xconj_action.Mdag(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of Xconjugate action Mdag against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xbarconj_action_ref.M(rand_sc, tmpsc); + xbarconj_action.M(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of Xbar-conjugate action M against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xbarconj_action_ref.Mdag(rand_sc, tmpsc); + xbarconj_action.Mdag(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of Xbar-conjugate action Mdag against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + RealD u, l; + + //Test the X-conjugate Dirac op acting on a field is the same as the G-parity Dirac op acting on the equivalent 2f field + xconj_action.M(rand_sc, tmpsc); + boostXconjToGparity(tmp, rand_sc); + reg_action.M(tmp, tmp2); + norm2DiffXconj(u,l,tmp2,tmpsc); + std::cout << "Test X-conj Dop reproduces G-parity Dop acting on X-conjugate field, f=0 (expect 0): " << u << " f=1 (expect 0): " << l << std::endl; + assert(l < 1e-10); + assert(u < 1e-10); + + //Test the Xbar-conjugate Dirac op acting on a field is the same as the G-parity Dirac op acting on the equivalent 2f field + xbarconj_action.M(rand_sc, tmpsc); + boostXbarConjToGparity(tmp, rand_sc); + reg_action.M(tmp, tmp2); + norm2DiffXbarConj(u,l,tmp2,tmpsc); + std::cout << "Test Xbar-conj Dop reproduces G-parity Dop acting on Xbar-conjugate field, f=0 (expect 0): " << u << " f=1 (expect 0): " << l << std::endl; + assert(l < 1e-10); + assert(u < 1e-10); + + //Test the X-conj 2f wrapper reproduces G-parity Dop acting on X-conjugate field + Xconjugate2fWrapper xconj_2f_wrapper(&xconj_action,®_action); + boostXconjToGparity(tmp, rand_sc); + xconj_2f_wrapper.M(tmp, tmp3); + reg_action.M(tmp, tmp2); + nrm = norm2Diff(tmp3,tmp2); + std::cout << "Test the X-conj 2f wrapper reproduces G-parity Dop acting on X-conjugate field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Show that U^dag * v with X-conj v is real + boostXconjToGparity(tmp, rand_sc); + applyUdag(tmp2,tmp); + tmp3 = tmp2 - conjugate(tmp2); + nrm = norm2(tmp3); + std::cout << "Test U^dag v with v an X-conj vector results in a real vector (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Check its form + tmpsc = PeekIndex(tmp2,0); + tmpsc2 = sqrt(2.)*(Xmatrix()*real(rand_sc)); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test U^dag v upper component has expected form (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + tmpsc = PeekIndex(tmp2,1); + tmpsc2 = sqrt(2.)*imag(rand_sc); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test U^dag v lower component has expected form (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Show that U is unitary + gaussian(RNG5,tmp); + applyUdag(tmp2,tmp); + applyU(tmp3,tmp2); + nrm = norm2Diff(tmp3,tmp); + std::cout << "Test U U^dag = 1 (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + applyU(tmp2,tmp); + applyUdag(tmp3,tmp2); + nrm = norm2Diff(tmp3,tmp); + std::cout << "Test U^dag U = 1 (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + { + //Show -U^* U^dag = -U U^T = Xi + gaussian(RNG5,tmp); + + tmp2 = conjugate(tmp); // U^T v = (U^dag v^*)^* + applyUdag(tmp3,tmp2); + tmp2 = conjugate(tmp3); + applyU(tmp3,tmp2); //U U^T v + tmp3 = -tmp3; + + //Xi = -i sigma2 X + GparityFlavour sigma2 = GparityFlavour(GparityFlavour::Algebra::SigmaY); + tmp2 = ComplexD(0,-1)*(sigma2*(Xmatrix()*tmp)); + nrm = norm2Diff(tmp3,tmp2); + std::cout << "Test Xi = -U U^T (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + } + + + { + //Show R = U^dag M U is a real matrix through + // (R v)^* = R v^* + FermionField2f Rv_allstar(FGrid), R_vstar(FGrid), v(FGrid); + gaussian(RNG5,v); + applyU(tmp, v); + reg_action.M(tmp, tmp2); + applyUdag(Rv_allstar,tmp2); //Rv + Rv_allstar = conjugate(Rv_allstar); + + tmp = conjugate(v); + applyU(tmp2,tmp); + reg_action.M(tmp2, tmp); + applyUdag(R_vstar,tmp); //Rv + + nrm = norm2Diff(R_vstar,Rv_allstar); + std::cout << "Test U^dag M U is a real matrix (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + } + + { + //Test relations between elements of M + //M11 = -X M00^* X + FermionField1f tmp(FGrid),tmp2(FGrid),tmp3(FGrid), v(FGrid), lhs(FGrid), rhs(FGrid); + gaussian(RNG5,v); + tmp = Xmatrix()*v; + applyMijstar(tmp2,tmp,0,0,reg_action); + rhs = -(Xmatrix()*tmp2); + + applyMij(lhs,v,1,1,reg_action); + nrm = norm2Diff(lhs, rhs); + std::cout << "Test M11 = -X M00* X (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //M10 = X M01* X + tmp = Xmatrix()*v; //X M01^* X v = (X M01 X v^* )^* + applyMijstar(tmp2,tmp,0,1,reg_action); + rhs = Xmatrix()*tmp2; + + applyMij(lhs,v,1,0,reg_action); + nrm = norm2Diff(lhs, rhs); + std::cout << "Test M10 = X M01* X (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + } + + { + //Test expressions for elements of R + //R11 = 0.5(M11 + M12 X + M12* X + M11*) + FermionField1f tmp(FGrid),tmp2(FGrid),tmp3(FGrid), v(FGrid), lhs(FGrid), rhs(FGrid); + gaussian(RNG5,v); + applyMij(rhs,v,0,0,reg_action); + + tmp = Xmatrix()*v; + applyMij(tmp2,tmp,0,1,reg_action); + rhs = rhs + tmp2; + + tmp = Xmatrix()*v; + applyMijstar(tmp2,tmp,0,1,reg_action); + rhs = rhs + tmp2; + + applyMijstar(tmp2,v,0,0,reg_action); + rhs = rhs + tmp2; + + rhs = rhs*0.5; + + applyRij(lhs,v,1,1,reg_action); + + nrm = norm2Diff(lhs, rhs); + std::cout << "Test R11 = 0.5(M11 + M12 X + M12* X + M11*) (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //R00 = 0.5( -X M00 X -X M01 -X M01* -X M00* X ) + tmp = Xmatrix()*v; + applyMij(tmp2,tmp,0,0,reg_action); + rhs = -(Xmatrix()*tmp2); + + applyMij(tmp,v,0,1,reg_action); + tmp2 = -(Xmatrix()*tmp); + rhs = rhs + tmp2; + + applyMijstar(tmp,v,0,1,reg_action); + tmp2 = -(Xmatrix()*tmp); + rhs = rhs + tmp2; + + tmp = Xmatrix()*v; + applyMijstar(tmp2,tmp,0,0,reg_action); + tmp = -(Xmatrix()*tmp2); + rhs = rhs + tmp; + + rhs = rhs * 0.5; + + applyRij(lhs,v,0,0,reg_action); + + nrm = norm2Diff(lhs, rhs); + std::cout << "Test R00 = 0.5( -X M00 X -X M01 -X M01* -X M00* X ) (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //R01 = 0.5*( iX M00 + iX M01 X -iX M01*X - iX M00* ) + applyMij(tmp,v,0,0,reg_action); + rhs = Xmatrix()*tmp; + + tmp = Xmatrix()*v; + applyMij(tmp2,tmp,0,1,reg_action); + tmp = Xmatrix()*tmp2; + rhs = rhs + tmp; + + tmp = Xmatrix()*v; + applyMijstar(tmp2,tmp,0,1,reg_action); + tmp = Xmatrix()*tmp2; + rhs = rhs - tmp; + + applyMijstar(tmp,v,0,0,reg_action); + tmp2 = Xmatrix()*tmp; + rhs = rhs - tmp2; + + rhs = rhs * ComplexD(0,0.5); + + applyRij(lhs,v,0,1,reg_action); + + nrm = norm2Diff(lhs, rhs); + std::cout << "Test R01 = 0.5*( iX M00 + iX M01 X -iX M01*X - iX M00* ) (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //R10 = 0.5( i M00 X + i M01 -i M01* - iM00* X ) + tmp = Xmatrix() * v; + applyMij(rhs,tmp,0,0,reg_action); + + applyMij(tmp,v,0,1,reg_action); + rhs = rhs + tmp; + + applyMijstar(tmp,v,0,1,reg_action); + rhs = rhs - tmp; + + tmp = Xmatrix() * v; + applyMijstar(tmp2,tmp,0,0,reg_action); + rhs = rhs - tmp2; + + rhs = rhs * ComplexD(0,0.5); + + applyRij(lhs,v,1,0,reg_action); + + nrm = norm2Diff(lhs, rhs); + std::cout << "Test R10 = 0.5( i M00 X + i M01 -i M01* - iM00* X ) (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + } + + { + //Test reference implementation + RotatedGPrefAction real_action_ref(®_action); + FermionField1f tmp(FGrid),tmp2(FGrid),tmp3(FGrid), v(FGrid), lhs(FGrid), rhs(FGrid); + gaussian(RNG5,v); + applyRijRef(lhs,v,0,0,real_action_ref); + applyRij(rhs,v,0,0,reg_action); + nrm = norm2Diff(lhs, rhs); + std::cout << "Test reference impl R00 (expect 0): " << nrm << std::endl; + //assert(nrm < 1e-10); + + applyRijRef(lhs,v,0,1,real_action_ref); + applyRij(rhs,v,0,1,reg_action); + nrm = norm2Diff(lhs, rhs); + std::cout << "Test reference impl R01 (expect 0): " << nrm << std::endl; + //assert(nrm < 1e-10); + + applyRijRef(lhs,v,1,0,real_action_ref); + applyRij(rhs,v,1,0,reg_action); + nrm = norm2Diff(lhs, rhs); + std::cout << "Test reference impl R10 (expect 0): " << nrm << std::endl; + //assert(nrm < 1e-10); + + applyRijRef(lhs,v,1,1,real_action_ref); + applyRij(rhs,v,1,1,reg_action); + nrm = norm2Diff(lhs, rhs); + std::cout << "Test reference impl R11 (expect 0): " << nrm << std::endl; + //assert(nrm < 1e-10); + } + + { + //Test reference implementation of U^dag M U + FermionField2f Rv_test(FGrid), Rv_ref(FGrid), v(FGrid); + gaussian(RNG5,v); + applyU(tmp, v); + reg_action.M(tmp, tmp2); + applyUdag(Rv_test,tmp2); //Rv + + RotatedGPrefAction real_action_ref(®_action); + real_action_ref.M(v,Rv_ref); + + nrm = norm2Diff(Rv_test,Rv_ref); + std::cout << "Test reference implementation of real matrix M (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + } + } + + //######################################################################## + { + std::cout << "Test of red-black preconditioned operator" << std::endl; + SchurDiagMooeeOperator reg_schurop(reg_action); + SchurDiagMooeeOperator xconj_schurop(xconj_action); + SchurDiagMooeeOperator xbarconj_schurop(xbarconj_action); + XconjugateDWF xconj_action_ref(®_action); + + FermionField1f rand_sc(FrbGrid), rand_sc2(FrbGrid); //v + gaussian(RNG5rb,rand_sc); + gaussian(RNG5rb,rand_sc2); + + FermionField2f rand_f0(FrbGrid); //v \delta_f,0 + rand_f0 = Zero(); + PokeIndex(rand_f0, rand_sc, 0); + + FermionField2f rand_f1(FrbGrid); //v \delta_f,1 + rand_f1 = Zero(); + PokeIndex(rand_f1, rand_sc, 1); + + FermionField2f tmp(FrbGrid), tmp2(FrbGrid), tmp3(FrbGrid), tmp4(FrbGrid); + FermionField1f tmpsc(FrbGrid), tmpsc2(FrbGrid), tmpsc3(FrbGrid); + RealD nrm; + + std::cout << "Test the relationship between the upper and lower rows in flavor space of the G-parity preconditioned Dirac operator" << std::endl; + + reg_schurop.Mpc(rand_f0,tmp); + FermionField1f M00v = PeekIndex(tmp,0); + FermionField1f M10v = PeekIndex(tmp,1); + + //Test M00 v = -X M11* X v = -X [M11 X v*]* + tmp = X*conjugate(rand_f1); + + reg_schurop.Mpc(tmp,tmp2); + FermionField1f M11Xvconj = PeekIndex(tmp2,1); + FermionField1f M01Xvconj = PeekIndex(tmp2,0); + + tmpsc = -(X*conjugate(M11Xvconj)); + nrm = norm2Diff(tmpsc, M00v); + std::cout << "Test of M00 v = -X M11* X v = -X [M11 X v*]* (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test M10 v = X M01* X v = X [M01 X v*]* + tmpsc = X*conjugate(M01Xvconj); + nrm = norm2Diff(tmpsc,M10v); + std::cout << "Test of M10 v = X M11* X v = X [M11 X v*]* (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + std::cout << "Test the relationship between the upper and lower rows in flavor space of the G-parity preconditioned squared Dirac operator" << std::endl; + + reg_schurop.HermOp(rand_f0,tmp); + M00v = PeekIndex(tmp,0); + M10v = PeekIndex(tmp,1); + + //Test M00 v = -X M11* X v = -X [M11 X v*]* + tmp = X*conjugate(rand_f1); + + reg_schurop.HermOp(tmp,tmp2); + M11Xvconj = PeekIndex(tmp2,1); + M01Xvconj = PeekIndex(tmp2,0); + + tmpsc = -(X*conjugate(M11Xvconj)); + nrm = norm2Diff(tmpsc,M00v); + std::cout << "Test of M00 v = -X M11* X v = -X [M11 X v*]* (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test M10 v = X M01* X v = X [M01 X v*]* + tmpsc = X*conjugate(M01Xvconj); + nrm = norm2Diff(tmpsc,M10v); + std::cout << "Test of M10 v = X M11* X v = X [M11 X v*]* (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test the X-conjugate implementation + std::cout << "Test the implementation of the X-conjugate preconditioned action against the reference"<< std::endl; + xconj_action_ref.Meooe(rand_sc, tmpsc); + xconj_action.Meooe(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of X-conjugate action Meooe against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xconj_action_ref.MeooeDag(rand_sc, tmpsc); + xconj_action.MeooeDag(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of X-conjugate action MeooeDag against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xconj_action_ref.Mooee(rand_sc, tmpsc); + xconj_action.Mooee(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of X-conjugate action Mooee against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xconj_action_ref.MooeeDag(rand_sc, tmpsc); + xconj_action.MooeeDag(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of X-conjugate action MooeeDag against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xconj_action_ref.MooeeInv(rand_sc, tmpsc); + xconj_action.MooeeInv(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of X-conjugate action MooeeInv against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + xconj_action_ref.MooeeInvDag(rand_sc, tmpsc); + xconj_action.MooeeInvDag(rand_sc, tmpsc2); + nrm = norm2Diff(tmpsc,tmpsc2); + std::cout << "Test of X-conjugate action MooeeInvDag against reference, full field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + RealD u,l; + + //Test the X-conjugate Dirac op acting on a field is the same as the G-parity Dirac op acting on a field with explicit X-conjugate BCs + xconj_schurop.HermOp(rand_sc, tmpsc); + boostXconjToGparity(tmp, rand_sc); + reg_schurop.HermOp(tmp, tmp2); + norm2DiffXconj(u,l,tmp2,tmpsc); + std::cout << "Test X-conj HermOp reproduces G-parity HermOp acting on X-conjugate field, f=0 (expect 0): " << u << " f=1 (expect 0): " << l << std::endl; + assert(u < 1e-10); + assert(l < 1e-10); + + //Test the X-conj 2f wrapper reproduces G-parity Dop acting on X-conjugate field + Xconjugate2fWrapper xconj_2f_wrapper(&xconj_action,®_action); + SchurDiagMooeeOperator, FermionField2f> xconj_2f_schurop(xconj_2f_wrapper); + boostXconjToGparity(tmp, rand_sc); + reg_schurop.HermOp(tmp, tmp2); + xconj_2f_schurop.HermOp(tmp,tmp3); + nrm = norm2Diff(tmp3,tmp2); + std::cout << "Test the X-conj 2f wrapper reproduces G-parity HermOp acting on X-conjugate field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test the Xbar-conjugate Dirac op acting on a field is the same as the G-parity Dirac op acting on a field with explicit Xbar-conjugate BCs + xbarconj_schurop.HermOp(rand_sc, tmpsc); + boostXbarConjToGparity(tmp,rand_sc); + reg_schurop.HermOp(tmp, tmp2); + norm2DiffXbarConj(u,l,tmp2,tmpsc); + std::cout << "Test Xbar-conj HermOp reproduces G-parity HermOp acting on Xbar-conjugate field, f=0 (expect 0): " << u << " f=1 (expect 0): " << l << std::endl; + assert(u < 1e-10); + assert(l < 1e-10); + + //Test the Xbar-conj 2f wrapper reproduces G-parity Dop acting on Xbar-conjugate field + Xconjugate2fWrapper xbarconj_2f_wrapper(&xbarconj_action,®_action); + SchurDiagMooeeOperator, FermionField2f> xbarconj_2f_schurop(xbarconj_2f_wrapper); + boostXbarConjToGparity(tmp, rand_sc); + reg_schurop.HermOp(tmp, tmp2); + xbarconj_2f_schurop.HermOp(tmp,tmp3); + nrm = norm2Diff(tmp3,tmp2); + std::cout << "Test the Xbar-conj 2f wrapper reproduces G-parity HermOp acting on Xbar-conjugate field (expect 0): " << nrm << std::endl; + assert(nrm < 1e-10); + + //Test reconstruction of G-parity Dop on arbitrary flavor vector using Xconj action + PokeIndex(tmp, rand_sc, 0); + PokeIndex(tmp, rand_sc2, 1); + reg_schurop.HermOp(tmp, tmp2); + + FermionField1f rho(FrbGrid), xi(FrbGrid); + rho = 0.5 * ( rand_sc + (X*conjugate(rand_sc2)) ); + xi = 0.5 * ( rand_sc - (X*conjugate(rand_sc2)) ); + xconj_schurop.HermOp(rho, tmpsc); + xbarconj_schurop.HermOp(xi, tmpsc2); + + tmpsc3 = PeekIndex(tmp2,0) - tmpsc - tmpsc2; + u = norm2(tmpsc3); + + tmpsc = -(X*conjugate(tmpsc)); + tmpsc2 = X*conjugate(tmpsc2); + tmpsc3 = PeekIndex(tmp2,1) - tmpsc - tmpsc2; + l = norm2(tmpsc3); + + std::cout << "Test reconstruction of GP HermOp on random 2f field from Xconj ops f=0 (expect 0): " << u << " f=1 (expect 0): " << l << std::endl; + assert(u < 1e-10); + assert(l < 1e-10); + } + std::cout << "All tests passed" << std::endl; + + Grid_finalize(); +} + diff --git a/tests/core/Test_gpdwf_Xconj_doublelatt.cc b/tests/core/Test_gpdwf_Xconj_doublelatt.cc new file mode 100644 index 0000000000..d44ede7064 --- /dev/null +++ b/tests/core/Test_gpdwf_Xconj_doublelatt.cc @@ -0,0 +1,400 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/core/Test_gpdwf_Xconj_doublelatt.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle + + 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 */ + +/** + This test uses the relationship between G-parity BCs, X-conjugate BCs and antiperiodic BCs on a doubled lattice + to perform a complete test of the implementations + **/ + +#include + +using namespace std; +using namespace Grid; + ; + +//typedef GparityDomainWallFermionD GparityDiracOp; +//typedef DomainWallFermionD StandardDiracOp; +//typedef XconjugateDomainWallFermionD XconjDiracOp; +//#define DOP_PARAMS + +typedef GparityMobiusFermionD GparityDiracOp; +typedef MobiusFermionD StandardDiracOp; +typedef XconjugateMobiusFermionD XconjDiracOp; +#define DOP_PARAMS ,1.5, 0.5 + + +typedef typename GparityDiracOp::FermionField GparityFermionField; +typedef typename GparityDiracOp::GaugeField GparityGaugeField; +typedef typename GparityFermionField::vector_type vComplexType; + +typedef typename StandardDiracOp::FermionField StandardFermionField; +typedef typename StandardDiracOp::GaugeField StandardGaugeField; + +enum{ same_vComplex = std::is_same::value }; +static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIMD complex type"); + +const Gamma & Xmatrix(){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + return X; +} + +void boostXconjToGparity(GparityFermionField &out, const StandardFermionField &in){ + PokeIndex(out, in, 0); + StandardFermionField tmp = -(Xmatrix()*conjugate(in)); + PokeIndex(out, tmp, 1); +} +void boostXconjToDoubleLatt(StandardFermionField &out, const StandardFermionField &in, const int nu, const int L){ + assert(out.Grid() != in.Grid()); + assert(!out.Grid()->_isCheckerBoarded); + assert(!in.Grid()->_isCheckerBoarded); + + LatticeInteger xcoor_2L_5(out.Grid()); + LatticeCoordinate(xcoor_2L_5,1+nu); //note '1+nu'! This is because for 5D fields the s-direction is direction 0 + + StandardFermionField tmp_2L(out.Grid()); + Replicate(in, out); + tmp_2L = -(Xmatrix()*conjugate(out)); + out = where( xcoor_2L_5 >= Integer(L), tmp_2L, out ); +} +void boostGparityToDoubleLatt(StandardFermionField &out, const GparityFermionField &in, const int nu, const int L){ + assert(out.Grid() != in.Grid()); + assert(!out.Grid()->_isCheckerBoarded); + assert(!in.Grid()->_isCheckerBoarded); + + LatticeInteger xcoor_2L_5(out.Grid()); + LatticeCoordinate(xcoor_2L_5,1+nu); //note '1+nu'! This is because for 5D fields the s-direction is direction 0 + + StandardFermionField tmp_L(in.Grid()); + tmp_L = PeekIndex(in,0); + Replicate(tmp_L, out); + + tmp_L = PeekIndex(in,1); + StandardFermionField tmp_2L(out.Grid()); + Replicate(tmp_L, tmp_2L); + + out = where( xcoor_2L_5 >= Integer(L), tmp_2L, out ); +} + +//Break up a larger lattice into smaller lattices +//output vector is resized to the number of subdivisions, with lexicographic block indexing x+rx*(y+ry*(z+rz*t)) +//where rx,ry,rz,rt are the grid size ratios +template +void SubDivide(std::vector > &out, const Lattice &in, GridBase* out_grid) +{ + out.resize(0,out_grid); + typedef typename vobj::scalar_object sobj; + + GridBase *og = out_grid; + GridBase *ig = in.Grid(); + + int nd = ig->_ndimension; + assert(og->_ndimension==nd); + + subdivides(og,ig); + + int nout = 1; + Coordinate ratio(nd); + for(int d=0;d_fdimensions[d]/og->_fdimensions[d]; + nout *= ratio[d]; + } + out.resize(nout, out_grid); + + Coordinate ocoor(nd); + Coordinate icoor(nd); + for(int g=0;ggSites();g++){ + + int oidx = 0; + ig->GlobalIndexToGlobalCoor(g,icoor); + for(int d=nd-1;d>=0;d--){ + ocoor[d] = icoor[d]%og->_gdimensions[d]; + oidx = icoor[d] / og->_gdimensions[d] + ratio[d]*oidx; //lex mapping x+rx*(y+ry*(z+rz*t)) + } + + sobj tmp; + peekSite(tmp,in,icoor); + pokeSite(tmp,out[oidx],ocoor); + } + +} + + + +int main (int argc, char ** argv) +{ + int nu = 0; + int tbc_aprd = 0; //use antiperiodic BCs in the time direction? + + Grid_init(&argc,&argv); + + for(int i=1;i> nu; + std::cout << GridLogMessage << "Set Gparity direction to " << nu << std::endl; + }else if(std::string(argv[i]) == "--Tbc-APRD"){ + tbc_aprd = 1; + std::cout << GridLogMessage << "Using antiperiodic BCs in the time direction" << std::endl; + } + } + + std::cout << GridLogMessage<< "*****************************************************************" < seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5_L(FGrid_L); RNG5_L.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4_L(UGrid_L); RNG4_L.SeedFixedIntegers(seeds4); + + GparityGaugeField Umu_L(UGrid_L); //note type is same as StandardGaugeField + SU::HotConfiguration(RNG4_L,Umu_L); + + //////////////////////// + //Copy-conjugate the gauge field + //////////////////////// + StandardGaugeField Umu_2L(UGrid_2L); + Replicate(Umu_L,Umu_2L); //output on right, grr! + + LatticeInteger xcoor_2L_4(UGrid_2L); + LatticeCoordinate(xcoor_2L_4,nu); + StandardGaugeField Umu_2L_conj = conjugate(Umu_2L); + Umu_2L = where( xcoor_2L_4 >= Integer(L), Umu_2L_conj, Umu_2L ); //don't need Cshift as replicate already duplicates U onto the doubled lattice + + //////////////////////// + //Generate two one-flavor fields with implicit X-conjugate BCs + //////////////////////// + + StandardFermionField src_X_L(FGrid_L), src2_X_L(FGrid_L); + random(RNG5_L,src_X_L); + random(RNG5_L,src2_X_L); + + //////////////////////// + //Create the corresponding two-flavor G-parity fields + //////////////////////// + GparityFermionField src_GP_L(FGrid_L), src2_GP_L(FGrid_L); + boostXconjToGparity(src_GP_L,src_X_L); + boostXconjToGparity(src2_GP_L,src2_X_L); + + //////////////////////// + //Create the corresponding doubled-lattice fields + //////////////////////// + StandardFermionField src_std_2L(FGrid_2L), src2_std_2L(FGrid_2L); + boostXconjToDoubleLatt(src_std_2L, src_X_L, nu, L); + boostXconjToDoubleLatt(src2_std_2L, src2_X_L, nu, L); + + //////////////////////// + //Create the Dirac operators + //////////////////////// + RealD mass=0.01; + RealD M5=1.8; + + //Standard Dirac op on doubled lattice + AcceleratorVector bc_std(Nd, 1.0); + if(tbc_aprd) bc_std[Nd-1] = -1.; //antiperiodic time BC + bc_std[nu] = -1; //antiperiodic in G-parity direction + StandardDiracOp::ImplParams std_params(bc_std); + StandardDiracOp Dstd_2L(Umu_2L,*FGrid_2L,*FrbGrid_2L,*UGrid_2L,*UrbGrid_2L,mass,M5 DOP_PARAMS, std_params); + + SchurDiagMooeeOperator HermOp_std_2L(Dstd_2L); + ConjugateGradient CG_std(1.0e-8,10000); + + //G-parity Dirac op on single lattice + std::vector twists(Nd,0); + twists[nu] = 1; + if(tbc_aprd) twists[Nd-1] = 1; + + GparityDiracOp::ImplParams gp_params; + gp_params.twists = twists; + GparityDiracOp DGP_L(Umu_L,*FGrid_L,*FrbGrid_L,*UGrid_L,*UrbGrid_L,mass,M5 DOP_PARAMS,gp_params); + + SchurDiagMooeeOperator HermOp_GP_L(DGP_L); + ConjugateGradient CG_GP(1.0e-8,10000); + + //X-conj Dirac op on single lattice + XconjDiracOp::ImplParams x_params; + x_params.twists = twists; + XconjDiracOp DX_L(Umu_L,*FGrid_L,*FrbGrid_L,*UGrid_L,*UrbGrid_L,mass,M5 DOP_PARAMS,x_params); + + SchurDiagMooeeOperator HermOp_X_L(DX_L); + + ///////////////////// + //Solve inverse with 1st source + ///////////////////// + StandardFermionField src_o_std_2L(FrbGrid_2L); + StandardFermionField result_o_std_2L(FrbGrid_2L); + pickCheckerboard(Odd,src_o_std_2L,src_std_2L); + result_o_std_2L=Zero(); + CG_std(HermOp_std_2L,src_o_std_2L,result_o_std_2L); + StandardFermionField result_std_2L(FGrid_2L); + result_std_2L = Zero(); + setCheckerboard(result_std_2L,result_o_std_2L); + + StandardFermionField src_o_X_L(FrbGrid_L); + StandardFermionField result_o_X_L(FrbGrid_L); + pickCheckerboard(Odd,src_o_X_L,src_X_L); + result_o_X_L=Zero(); + CG_std(HermOp_X_L,src_o_X_L,result_o_X_L); + StandardFermionField result_X_L(FGrid_L); + result_X_L = Zero(); + setCheckerboard(result_X_L,result_o_X_L); + + GparityFermionField src_o_GP_L(FrbGrid_L); + GparityFermionField result_o_GP_L(FrbGrid_L); + pickCheckerboard(Odd,src_o_GP_L,src_GP_L); + result_o_GP_L=Zero(); + CG_GP(HermOp_GP_L,src_o_GP_L,result_o_GP_L); + GparityFermionField result_GP_L(FGrid_L); + result_GP_L = Zero(); + setCheckerboard(result_GP_L,result_o_GP_L); + + std::cout << "Norms: std:" << norm2(result_o_std_2L) << " X:" << 2*norm2(result_o_X_L) << " GP:" << norm2(result_o_GP_L) << std::endl; + + ///////////////////// + //Convert X-conj and GP solutions to 2L formulation for direct comparison + //////////////////// + std::cout << "Computing differences between solutions in 2L formulation" << std::endl; + StandardFermionField result_X_2L(FGrid_2L); + boostXconjToDoubleLatt(result_X_2L, result_X_L, nu, L); + + StandardFermionField result_GP_2L(FGrid_2L); + boostGparityToDoubleLatt(result_GP_2L, result_GP_L, nu, L); + + StandardFermionField diff_2L(FGrid_2L); + diff_2L = result_X_2L - result_std_2L; + std::cout << "X vs std: " << norm2(diff_2L) << std::endl; + assert(norm2(diff_2L) < 1e-8); + diff_2L = result_GP_2L - result_std_2L; + std::cout << "GP vs std: " << norm2(diff_2L) << std::endl; + assert(norm2(diff_2L) < 1e-8); + + //////////////////////////////// + //Repeat experiments with a G-parity source, back convert to X-conjugate via doubled lattice as in doc + //////////////////////////////// + std::cout << "Using a random two-flavor G-parity field" << std::endl; + + random(RNG5_L,src_GP_L); + boostGparityToDoubleLatt(src_std_2L, src_GP_L, nu, L); + std::vector src_std_2L_split; + SubDivide(src_std_2L_split, src_std_2L, FGrid_L); + assert(src_std_2L_split.size() == 2); + + //v[0] = a+ib + //v[1] = -X[a-ib]* + StandardFermionField tmp_X_L(FGrid_L), tmp2_X_L(FGrid_L); + tmp_X_L = Xmatrix()*conjugate(src_std_2L_split[1]); + src_X_L = 0.5*( src_std_2L_split[0] + tmp_X_L ); //a + src2_X_L = ComplexD(0,-0.5)*( src_std_2L_split[0] - tmp_X_L ); //b + + ///////////////////// + //Solve inverse with new sources + ///////////////////// + pickCheckerboard(Odd,src_o_std_2L,src_std_2L); + result_o_std_2L=Zero(); + CG_std(HermOp_std_2L,src_o_std_2L,result_o_std_2L); + result_std_2L = Zero(); + setCheckerboard(result_std_2L,result_o_std_2L); + + //note 2 inversions for Xconj + pickCheckerboard(Odd,src_o_X_L,src_X_L); + result_o_X_L=Zero(); + CG_std(HermOp_X_L,src_o_X_L,result_o_X_L); + result_X_L = Zero(); + setCheckerboard(result_X_L,result_o_X_L); + + StandardFermionField src2_o_X_L(FrbGrid_L); + StandardFermionField result2_o_X_L(FrbGrid_L); + pickCheckerboard(Odd,src2_o_X_L,src2_X_L); + result2_o_X_L=Zero(); + CG_std(HermOp_X_L,src2_o_X_L,result2_o_X_L); + StandardFermionField result2_X_L(FGrid_L); + result2_X_L = Zero(); + setCheckerboard(result2_X_L,result2_o_X_L); + + //Gparity + pickCheckerboard(Odd,src_o_GP_L,src_GP_L); + result_o_GP_L=Zero(); + CG_GP(HermOp_GP_L,src_o_GP_L,result_o_GP_L); + result_GP_L = Zero(); + setCheckerboard(result_GP_L,result_o_GP_L); + + /////////////////// + //Convert to doubled latt + /////////////////// + std::cout << "Computing differences between solutions in 2L formulation" << std::endl; + boostXconjToDoubleLatt(result_X_2L, result_X_L, nu, L); + //Note we need to reconstruct the full solution + StandardFermionField tmp_X_2L(FGrid_2L); + boostXconjToDoubleLatt(tmp_X_2L, result2_X_L, nu, L); + result_X_2L = result_X_2L + ComplexD(0,1)*tmp_X_2L; + + boostGparityToDoubleLatt(result_GP_2L, result_GP_L, nu, L); + + diff_2L = result_X_2L - result_std_2L; + std::cout << "X vs std: " << norm2(diff_2L) << std::endl; + assert(norm2(diff_2L) < 1e-8); + diff_2L = result_GP_2L - result_std_2L; + std::cout << "GP vs std: " << norm2(diff_2L) << std::endl; + assert(norm2(diff_2L) < 1e-8); + + Grid_finalize(); +} diff --git a/tests/forces/Test_gpdwf_Xconj_force.cc b/tests/forces/Test_gpdwf_Xconj_force.cc new file mode 100644 index 0000000000..2ab1e89227 --- /dev/null +++ b/tests/forces/Test_gpdwf_Xconj_force.cc @@ -0,0 +1,221 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/forces/Test_gpdwf_Xconj_force.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + +typedef GparityDomainWallFermionD::FermionField FermionField2f; +typedef XconjugateDomainWallFermionD::FermionField FermionField1f; + +const Gamma & Xmatrix(){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + return X; +} + +void XconjBoost(FermionField2f &to, const FermionField1f &from){ + FermionField1f tmp = -(Xmatrix()*conjugate(from)); + PokeIndex(to, from, 0); + PokeIndex(to, tmp, 1); + to.Checkerboard() = from.Checkerboard(); +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + std::vector seeds4({1,2,3,5}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + int threads = GridThread::GetThreads(); + std::cout<::HotConfiguration(RNG4,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD mass=0.01; + RealD M5=1.8; + + XconjugateDomainWallFermionD::ImplParams xparams; + std::vector twists({1,1,1,1}); //GPBC in 3 dirs, antiperiodic in time + xparams.twists = twists; + xparams.boundary_phase = 1.0; + + XconjugateDomainWallFermionD Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,xparams); + + GparityDomainWallFermionD::ImplParams params; params.twists = twists; + GparityDomainWallFermionD Ddwf_gp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); + typedef GparityDomainWallFermionD::FermionField FermionField2f; + + //Action of X-conjugate operator + LatticeFermion phi (FGrid); gaussian(RNG5,phi); + LatticeFermion Mphi (FGrid); + LatticeFermion MphiPrime (FGrid); + Ddwf.M (phi,Mphi); + ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p + + //Check against action of G-parity operator + FermionField2f phi_gp (FGrid); + FermionField2f Mphi_gp (FGrid); + FermionField2f MphiPrime_gp (FGrid); + XconjBoost(phi_gp, phi); + phi_gp = phi_gp * sqrt(0.5); //rho from rho' + + Ddwf_gp.M (phi_gp,Mphi_gp); + ComplexD S_gp = innerProduct(Mphi_gp,Mphi_gp); // pdag MdagM p + + std::cout << GridLogMessage << "Initial action Xconjugate: " << S << " Gparity: " << S_gp << " diff: " << S_gp - S << std::endl; + + + // get the deriv of phidag MdagM phi with respect to "U" + LatticeGaugeField UdSdU(UGrid); + LatticeGaugeField tmp(UGrid); + Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; + Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + + //and the same for the G-parity operator + LatticeGaugeField UdSdU_gp(UGrid); + LatticeGaugeField tmp_gp(UGrid); + Ddwf_gp.MDeriv(tmp_gp , Mphi_gp, phi_gp,DaggerNo ); UdSdU_gp=tmp_gp; + Ddwf_gp.MDeriv(tmp_gp , phi_gp, Mphi_gp,DaggerYes ); UdSdU_gp=(UdSdU_gp+tmp_gp); + + LatticeGaugeField force_diff = UdSdU_gp - UdSdU; + std::cout << GridLogMessage << "Force Xconjugate: " << norm2(UdSdU) << " G-parity: " << norm2(UdSdU_gp) << " diff: " << norm2(force_diff) << std::endl; + + + + LatticeFermion Ftmp (FGrid); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.0001; + + LatticeColourMatrix mommu(UGrid); + LatticeColourMatrix forcemu(UGrid); + LatticeGaugeField mom(UGrid); + LatticeGaugeField Uprime(UGrid); + + for(int mu=0;mu::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg + + PokeIndex(mom,mommu,mu); + + // fourth order exponential approx + + autoView( mom_v, mom, CpuRead); + autoView( U_v , U, CpuRead); + autoView(Uprime_v, Uprime, CpuWrite); + + thread_foreach( i,mom_v,{ + Uprime_v[i](mu) = U_v[i](mu) + + mom_v[i](mu)*U_v[i](mu)*dt + + mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt/2.0) + + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt/6.0) + + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt/24.0) + + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt/120.0) + + mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *mom_v[i](mu) *U_v[i](mu)*(dt*dt*dt*dt*dt*dt/720.0) + ; + }); + } + + Ddwf.ImportGauge(Uprime); + Ddwf.M (phi,MphiPrime); + ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime); + + Ddwf_gp.ImportGauge(Uprime); + Ddwf_gp.M (phi_gp,MphiPrime_gp); + ComplexD Sprime_gp = innerProduct(MphiPrime_gp ,MphiPrime_gp); + + std::cout << GridLogMessage << "Final action Xconjugate: " << Sprime << " Gparity: " << Sprime_gp << " diff: " << Sprime_gp - Sprime << std::endl; + + + ////////////////////////////////////////////// + // Use derivative to estimate dS + ////////////////////////////////////////////// + + LatticeComplex dS(UGrid); dS = Zero(); + LatticeComplex dS_gp(UGrid); dS_gp = Zero(); + + for(int mu=0;mu(UdSdU,mu); + mommu=Ta(mommu)*2.0; + PokeIndex(UdSdU,mommu,mu); + + mommu = PeekIndex(UdSdU_gp,mu); + mommu=Ta(mommu)*2.0; + PokeIndex(UdSdU_gp,mommu,mu); + } + + for(int mu=0;mu(mom,mu); + + // Update PF action density + forcemu = PeekIndex(UdSdU,mu); + dS = dS+trace(mommu*forcemu)*dt; + + forcemu = PeekIndex(UdSdU_gp,mu); + dS_gp = dS_gp+trace(mommu*forcemu)*dt; + } + + ComplexD dSpred = sum(dS); + ComplexD dSpred_gp = sum(dS_gp); + std::cout << std::setprecision(14); + + std::cout << GridLogMessage << " Xconj Gparity Diff" << std::endl; + std::cout << GridLogMessage << "S "<< S<< " " << S_gp << " " << S_gp-S << std::endl; + std::cout << GridLogMessage << "Sprime "<< Sprime << " " << Sprime_gp << " " << Sprime_gp - Sprime << std::endl; + std::cout << GridLogMessage << "dS "<< Sprime-S << " " << Sprime_gp - S_gp << " " << (Sprime-S)-(Sprime_gp-S_gp) << std::endl; + std::cout << GridLogMessage << "predict dS "<< dSpred << " " << dSpred_gp << " " << dSpred_gp - dSpred << std::endl; + + assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ; + + std::cout<< GridLogMessage << "Done" < +Author: Dennis Bollweg +Author: Peter Boyle + + 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 */ + +/** + * Tests for the action of X-conjugate Dirac operator + * + */ + +#include + +using namespace std; +using namespace Grid; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=4; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + auto const & latt_size = GridDefaultLatt(); + size_t vol4d = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3]; + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeGaugeField Umu(UGrid); + SU::HotConfiguration(RNG4, Umu); + + Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + Gamma X = C*g5; + + //Set up a regular MDWF action instance as well as X-conj and Xbar-conj versions + RealD m1=0.1; + RealD m2=1.0; + + RealD M5=1.8; + RealD mob_b=1.5; + + typedef typename GparityMobiusFermionD::FermionField FermionField2f; + typedef typename XconjugateMobiusFermionD::FermionField FermionField1f; + + GparityMobiusFermionD ::ImplParams params; + std::vector twists({1,1,0,1}); //G-parity in x,y periodic in z, antiperiodic in t + params.twists = twists; + GparityMobiusFermionD reg_action_m1(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,m1,M5,mob_b,mob_b-1.,params); + GparityMobiusFermionD reg_action_m2(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,m2,M5,mob_b,mob_b-1.,params); + + XconjugateMobiusFermionD::ImplParams xparams; + xparams.twists = twists; + xparams.boundary_phase = 1.0; + XconjugateMobiusFermionD xconj_action_m1(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,m1,M5,mob_b,mob_b-1.,xparams); + XconjugateMobiusFermionD xconj_action_m2(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,m2,M5,mob_b,mob_b-1.,xparams); + + xparams.boundary_phase = -1.0; + XconjugateMobiusFermionD xbarconj_action_m1(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,m1,M5,mob_b,mob_b-1.,xparams); + XconjugateMobiusFermionD xbarconj_action_m2(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,m2,M5,mob_b,mob_b-1.,xparams); + + //Gauge BCs + typedef ConjugateGimplD GimplD; + std::vector gauge_twists = twists; + gauge_twists[3] = 0; //periodic in time + GimplD::setDirections(twists); + + typedef typename GparityMobiusFermionD::Impl_t FImplD; + typedef typename XconjugateMobiusFermionD::Impl_t XFImplD; + + //First demonstrate that the X-conjugate action is computed correctly + //It should be the same as the G-parity Dirac operator applied to an X-conjugate vector + { + typedef TwoFlavourEvenOddRatioPseudoFermionAction XconjTwoFlavorHMC; + typedef TwoFlavourEvenOddRatioPseudoFermionAction GparityTwoFlavorHMC; + + ConjugateGradient CG2f(1e-10,10000); + ConjugateGradient CG1f(1e-10,10000); + + XconjTwoFlavorHMC hmc_xconj(xconj_action_m2, xconj_action_m1, CG1f, CG1f); + GparityTwoFlavorHMC hmc_gp(reg_action_m2, reg_action_m1, CG2f, CG2f); + + RealD scale = std::sqrt(0.5); + FermionField1f eta1f(FGrid); + gaussian(RNG5,eta1f); + + FermionField2f eta2f(FGrid); + + FermionField1f tmp = -(X*conjugate(eta1f)); + PokeIndex(eta2f, eta1f, 0); + PokeIndex(eta2f, tmp, 1); + + hmc_xconj.setPhi(eta1f); + RealD S_xconj = 2. * hmc_xconj.S(Umu); + + hmc_gp.setPhi(eta2f); + RealD S_gp = hmc_gp.S(Umu); + + std::cout << "Test 2 rho^dag [M_X M_X^dag]^-1 rho = \\hat rho^dag [M_GP M_GP^dag] \\hat rho (expect 0): " << S_gp - S_xconj << std::endl; + assert(fabs(S_gp - S_xconj) < 1e-8); + } + //Next demonstrate the G-parity Dirac operator acting on a general 2f vector can be computed in terms of the X-conjugate action + { + typedef TwoFlavourEvenOddRatioPseudoFermionAction XconjTwoFlavorHMC; + typedef TwoFlavourEvenOddRatioPseudoFermionAction GparityTwoFlavorHMC; + + ConjugateGradient CG2f(1e-10,10000); + ConjugateGradient CG1f(1e-10,10000); + + XconjTwoFlavorHMC hmc_xconj(xconj_action_m2, xconj_action_m1, CG1f, CG1f); + XconjTwoFlavorHMC hmc_xbarconj(xbarconj_action_m2, xbarconj_action_m1, CG1f, CG1f); + GparityTwoFlavorHMC hmc_gp(reg_action_m2, reg_action_m1, CG2f, CG2f); + + RealD scale = std::sqrt(0.5); + FermionField2f eta2f(FGrid); + gaussian(RNG5,eta2f); + + FermionField1f eta2f_0 = PeekIndex(eta2f,0); + FermionField1f eta2f_1 = PeekIndex(eta2f,1); + + FermionField1f rho1f = 0.5*( eta2f_0 + X*conjugate(eta2f_1) ); + FermionField1f tau1f = Complex(0,0.5)*( eta2f_0 - X*conjugate(eta2f_1) ); + + hmc_gp.setPhi(eta2f); + RealD S_gp = hmc_gp.S(Umu); + + hmc_xconj.setPhi(rho1f); + RealD S_xconj = 2.*hmc_xconj.S(Umu); + hmc_xconj.setPhi(tau1f); + RealD S_xconj2 = 2.*hmc_xconj.S(Umu); + + RealD S_xtotal = S_xconj + S_xconj2; + + std::cout << "Test eta^dag [M_GP M_GP^dag]^-1 eta = 2 rho^dag [M_X M_X^dag]^-1 rho + 2 tau^dag [M_X M_X^dag]^-1 tau (expect 0): " << S_xtotal-S_gp << std::endl; + assert(fabs(S_xtotal-S_gp) < 1e-8); + + //The above can also be computed with the Xbar-conjugate action after transforming the source + FermionField1f xi1f = Complex(0,-1)*tau1f; + hmc_xbarconj.setPhi(xi1f); + RealD S_xbarconj = 2.*hmc_xbarconj.S(Umu); + + S_xtotal = S_xconj + S_xbarconj; + + std::cout << "Test eta^dag [M_GP M_GP^dag]^-1 eta = 2 rho^dag [M_X M_X^dag]^-1 rho + 2 xi^dag [M_Xbar M_Xbar^dag]^-1 xi (expect 0): " << S_xtotal-S_gp << std::endl; + assert(fabs(S_xtotal-S_gp) < 1e-8); + } + //Do the same but with a rational approximation to (M M^dag)^-1/2 + { + typedef GeneralEvenOddRatioRationalPseudoFermionAction GparityRHMC; + typedef GeneralEvenOddRatioRationalPseudoFermionAction XconjRHMC; + + RationalActionParams rat_act_params; + rat_act_params.inv_pow = 2; // (M^dag M)^{1/2} + rat_act_params.precision= 60; + rat_act_params.MaxIter = 50000; + rat_act_params.action_degree = 15; + rat_act_params.action_tolerance = 1.0e-10; + rat_act_params.md_degree = 12; + rat_act_params.md_tolerance = 1.0e-5; + rat_act_params.lo = 0.01; + rat_act_params.hi = 80.0; + + GparityRHMC rhmc_gp(reg_action_m2, reg_action_m1, rat_act_params); + XconjRHMC rhmc_xconj(xconj_action_m2, xconj_action_m1, rat_act_params); + XconjRHMC rhmc_xbarconj(xbarconj_action_m2, xbarconj_action_m1, rat_act_params); + + RealD scale = std::sqrt(0.5); + FermionField2f eta2f(FGrid); + gaussian(RNG5,eta2f); + + FermionField1f eta2f_0 = PeekIndex(eta2f,0); + FermionField1f eta2f_1 = PeekIndex(eta2f,1); + + FermionField1f rho1f = 0.5*( eta2f_0 + X*conjugate(eta2f_1) ); + FermionField1f tau1f = Complex(0,0.5)*( eta2f_0 - X*conjugate(eta2f_1) ); + + rhmc_gp.setPhi(eta2f); + RealD S_gp = rhmc_gp.S(Umu); + + rhmc_xconj.setPhi(rho1f); + RealD S_xconj = 2.*rhmc_xconj.S(Umu); + rhmc_xconj.setPhi(tau1f); + RealD S_xconj2 = 2.*rhmc_xconj.S(Umu); + + RealD S_xtotal = S_xconj + S_xconj2; + + std::cout << "Test eta^dag [M_GP M_GP^dag]^-1/2 eta = 2 rho^dag [M_X M_X^dag]^-1/2 rho + 2 tau^dag [M_X M_X^dag]^-1/2 tau using RHMC (expect 0): " << S_xtotal-S_gp << std::endl; + assert(fabs(S_xtotal-S_gp) < 1e-8); + + //The above can also be computed with the Xbar-conjugate action after transforming the source + FermionField1f xi1f = Complex(0,-1)*tau1f; + rhmc_xbarconj.setPhi(xi1f); + RealD S_xbarconj = 2.*rhmc_xbarconj.S(Umu); + + S_xtotal = S_xconj + S_xbarconj; + + std::cout << "Test eta^dag [M_GP M_GP^dag]^-1/2 eta = 2 rho^dag [M_X M_X^dag]^-1/2 rho + 2 xi^dag [M_Xbar M_Xbar^dag]^-1/2 xi using RHMC (expect 0): " << S_xtotal-S_gp << std::endl; + assert(fabs(S_xtotal-S_gp) < 1e-8); + } + + //Examine the EOFA Dirac operator + { + typedef XconjugateMobiusEOFAFermionD XconjEOFAaction; + typedef ExactOneFlavourRatioPseudoFermionAction XconjEOFA; + typedef GparityMobiusEOFAFermionD GparityEOFAaction; + typedef ExactOneFlavourRatioPseudoFermionAction GparityEOFA; + + ConjugateGradient CG2f(1e-10,10000); + ConjugateGradient CG1f(1e-10,10000); + + RealD iml = m1; + RealD ipv = m2; + + OneFlavourRationalParams eofa_params; //not needed because we are not doing the refresh + eofa_params.lo = 1.0; + eofa_params.hi = 2.0; + eofa_params.degree = 12; + eofa_params.tolerance = 1e-4; //this is just for doing the Remez + + xparams.boundary_phase = 1.0; + XconjEOFAaction Lop_xconj(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, iml, iml, ipv, 0.0, -1, M5, mob_b, mob_b-1.,xparams); + XconjEOFAaction Rop_xconj(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, ipv, iml, ipv, -1.0, 1, M5, mob_b, mob_b-1.,xparams); + XconjEOFA eofa_xconj(Lop_xconj, Rop_xconj, CG1f, eofa_params); + + xparams.boundary_phase = -1.0; + XconjEOFAaction Lop_xbarconj(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, iml, iml, ipv, 0.0, -1, M5, mob_b, mob_b-1.,xparams); + XconjEOFAaction Rop_xbarconj(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, ipv, iml, ipv, -1.0, 1, M5, mob_b, mob_b-1.,xparams); + XconjEOFA eofa_xbarconj(Lop_xbarconj, Rop_xbarconj, CG1f, eofa_params); + + GparityEOFAaction Lop_gp(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, iml, iml, ipv, 0.0, -1, M5, mob_b, mob_b-1.,params); + GparityEOFAaction Rop_gp(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, ipv, iml, ipv, -1.0, 1, M5, mob_b, mob_b-1.,params); + GparityEOFA eofa_gp(Lop_gp, Rop_gp, CG2f, eofa_params); + + RealD scale = std::sqrt(0.5); + FermionField2f eta2f(FGrid); + gaussian(RNG5,eta2f); + + FermionField1f eta2f_0 = PeekIndex(eta2f,0); + FermionField1f eta2f_1 = PeekIndex(eta2f,1); + + FermionField1f rho1f = 0.5*( eta2f_0 + X*conjugate(eta2f_1) ); + FermionField1f tau1f = Complex(0,0.5)*( eta2f_0 - X*conjugate(eta2f_1) ); + + eofa_gp.setPhi(eta2f); + RealD S_gp = eofa_gp.S(Umu); + + eofa_xconj.setPhi(rho1f); + RealD S_xconj = 2.*eofa_xconj.S(Umu); + eofa_xconj.setPhi(tau1f); + RealD S_xconj2 = 2.*eofa_xconj.S(Umu); + + RealD S_xtotal = S_xconj + S_xconj2; + + std::cout << "Test eta^dag [M_GP M_GP^dag]^-1/2 eta = 2 rho^dag [M_X M_X^dag]^-1/2 rho + 2 tau^dag [M_X M_X^dag]^-1/2 tau using EOFA (expect 0): " << S_xtotal-S_gp << std::endl; + assert(fabs(S_xtotal-S_gp) < 1e-8); + + //The above can also be computed with the Xbar-conjugate action after transforming the source + FermionField1f xi1f = Complex(0,-1)*tau1f; + eofa_xbarconj.setPhi(xi1f); + RealD S_xbarconj = 2.*eofa_xbarconj.S(Umu); + + S_xtotal = S_xconj + S_xbarconj; + + std::cout << "Test eta^dag [M_GP M_GP^dag]^-1/2 eta = 2 rho^dag [M_X M_X^dag]^-1/2 rho + 2 xi^dag [M_Xbar M_Xbar^dag]^-1/2 xi using EOFA (expect 0): " << S_xtotal-S_gp << std::endl; + assert(fabs(S_xtotal-S_gp) < 1e-8); + } + + + + + Grid_finalize(); +} diff --git a/tests/solver/Test_gpdwf_Xconj_inverse.cc b/tests/solver/Test_gpdwf_Xconj_inverse.cc new file mode 100644 index 0000000000..d963a783c3 --- /dev/null +++ b/tests/solver/Test_gpdwf_Xconj_inverse.cc @@ -0,0 +1,702 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/solver/Test_gpdwf_Xconj_inverse.cc + + Copyright (C) 2015 + +Author: Christopher Kelly +Author: Peter Boyle + + 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 */ + +/** + * Tests for the implementation of X-conjugate BCs + * + */ + +#include + +using namespace std; +using namespace Grid; + +typedef typename GparityMobiusFermionD::FermionField FermionField2f; +typedef typename XconjugateMobiusFermionD::FermionField FermionField1f; + +void invertGparity(std::vector &out, GparityMobiusFermionD &action){ + GridBase* UGrid = action.GaugeGrid(); + GridBase* FGrid = action.FermionGrid(); + GridBase* FrbGrid = action.FermionRedBlackGrid(); + + ConjugateGradient cg(1e-08,10000); + SchurRedBlackDiagMooeeSolve solver(cg); + + typedef FermionField2f::scalar_object Spinor; + + LatticeInteger tcoor4d(UGrid); + LatticeCoordinate(tcoor4d,3); + + FermionField2f tmp4d(UGrid); + FermionField2f src4d(UGrid); + FermionField2f src5d(FGrid), sol5d(FGrid); + FermionField2f src5d_e(FrbGrid), src5d_o(FrbGrid), sol5d_o(FrbGrid); + + out = std::vector(24, UGrid); + + for(int f=0;f<2;f++){ + for(int s=0;s<4;s++){ + for(int c=0;c<3;c++){ + Spinor v = Zero(); + v(f)(s)(c) = 1.; + + tmp4d = v; + src4d = Zero(); + src4d = where( tcoor4d == Integer(0), tmp4d, src4d); //unit vector on every site on timeslice 0, zero elsewhere + + action.ImportPhysicalFermionSource(src4d, src5d); + solver.RedBlackSource(action, src5d, src5d_e, src5d_o); + solver.RedBlackSolve(action, src5d_o, sol5d_o); + solver.RedBlackSolution(action, sol5d_o, src5d_e, sol5d); + action.ExportPhysicalFermionSolution(sol5d, out[c+3*(s + 4*f)]); + } + } + } +} + +//0.5*(1+iX) in +template +T mulPplusLeft(const T &in){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static ComplexD _I(0,1); + T out = 0.5*(in + _I * (X*in)); + return out; +} +//0.5*(1+iX) in +template +T mulPminusLeft(const T &in){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static ComplexD _mI(0,-1); + T out = 0.5*(in + _mI * (X*in)); + return out; +} + +//0.5*(1+iX) in +template +T mulPplusRight(const T &in){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static ComplexD _I(0,1); + T out = 0.5*(in + _I * (in*X)); + return out; +} +//0.5*(1+iX) in +template +T mulPminusRight(const T &in){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static ComplexD _mI(0,-1); + T out = 0.5*(in + _mI * (in*X)); + return out; +} + + + +typedef iScalar, Ns> > vSCmatrix; +typedef iScalar, Ns> > vSCvector; +typedef Lattice SCmatrixField; + +//Poke the spin-color vector field 'from' onto the sc,cc spin/color column of 'into' +void pokeSpinColorColumn(SCmatrixField &into, const FermionField1f &from, const int sc, const int cc){ + size_t Nsimd = FermionField1f::vector_type::Nsimd(); + autoView( from_v, from, AcceleratorRead); + autoView( into_v, into, AcceleratorWrite); + accelerator_for( i, from_v.size(), Nsimd, { + auto site_from = from_v(i); + auto site_into = into_v(i); + for(int sr=0;sr, Ns>, Ngp> vSCFmatrix; +typedef iVector, Ns>, Ngp> vSCFvector; +typedef Lattice SCFmatrixField; + +//Poke the spin-color-flavor vector field 'from' onto the sc,cc,ff spin/color/flavor column of 'into' +void pokeSpinColorFlavorColumn(SCFmatrixField &into, const FermionField2f &from, const int sc, const int cc, const int fc){ + size_t Nsimd = FermionField2f::vector_type::Nsimd(); + autoView( from_v, from, AcceleratorRead); + autoView( into_v, into, AcceleratorWrite); + accelerator_for( i, from_v.size(), Nsimd, { + auto site_from = from_v(i); + auto site_into = into_v(i); + for(int fr=0;fr &out, GparityMobiusFermionD &action){ + GridBase* UGrid = action.GaugeGrid(); + GridBase* FGrid = action.FermionGrid(); + GridBase* FrbGrid = action.FermionRedBlackGrid(); + + ConjugateGradient cg(1e-08,10000); + SchurRedBlackDiagMooeeSolve solver(cg); + + typedef FermionField2f::scalar_object Spinor; + typedef iScalar, Ns> > FSpinor; + + LatticeInteger tcoor4d(UGrid); + LatticeCoordinate(tcoor4d,3); + + FermionField2f tmp4d(UGrid); + FermionField2f src4d(UGrid); + FermionField2f src5d(FGrid), sol5d(FGrid); + FermionField2f src5d_e(FrbGrid), src5d_o(FrbGrid), sol5d_o(FrbGrid); + + SCmatrixField tmp_vplus_f0_4d(UGrid); + SCmatrixField tmp_vminus_f0_4d(UGrid); + + FermionField1f tmp1f(UGrid); + + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + + for(int s=0;s<4;s++){ + for(int c=0;c<3;c++){ + FSpinor vb = Zero(); + vb()(s)(c) = 1.; + + Spinor vplus; + vplus(0) = ( mulPplusLeft(vb) )(); + vplus(1) = ( -(X*mulPminusLeft(vb)) )(); + + Spinor vminus; + vminus(0) = ( mulPminusLeft(vb) )(); + vminus(1) = ( -(X*mulPplusLeft(vb)) )(); + + Spinor* vpm[2] = {&vplus, &vminus}; + SCmatrixField* solpm[2] = {&tmp_vplus_f0_4d, &tmp_vminus_f0_4d}; + + for(int pm=0;pm<2;pm++){ + tmp4d = *vpm[pm]; + src4d = Zero(); + src4d = where( tcoor4d == Integer(0), tmp4d, src4d); //vector on every site on timeslice 0, zero elsewhere + + action.ImportPhysicalFermionSource(src4d, src5d); + solver.RedBlackSource(action, src5d, src5d_e, src5d_o); + solver.RedBlackSolve(action, src5d_o, sol5d_o); + solver.RedBlackSolution(action, sol5d_o, src5d_e, sol5d); + action.ExportPhysicalFermionSolution(sol5d, tmp4d); + tmp1f = PeekIndex(tmp4d,0); //only need f0 component + pokeSpinColorColumn(*solpm[pm], tmp1f, s,c); + } + } + } + + SCmatrixField vplus_f0_Pplus = mulPplusRight(tmp_vplus_f0_4d); + SCmatrixField vplus_f0_Pminus = mulPminusRight(tmp_vplus_f0_4d); + + SCmatrixField vminus_f0_Pplus = mulPplusRight(tmp_vminus_f0_4d); + SCmatrixField vminus_f0_Pminus = mulPminusRight(tmp_vminus_f0_4d); + + SCmatrixField Minv11 = vplus_f0_Pplus + vminus_f0_Pminus; + SCmatrixField Minv12 = vplus_f0_Pminus*X + vminus_f0_Pplus*X; + SCmatrixField Minv22 = -(X*(conjugate(Minv11)*X)); + SCmatrixField Minv21 = X*(conjugate(Minv12)*X); + + out = std::vector(24, UGrid); + + SCmatrixField* Minv[2][2] = { {&Minv11,&Minv12}, {&Minv21,&Minv22} }; + + for(int fc=0;fc<2;fc++){ + for(int sc=0;sc<4;sc++){ + for(int cc=0;cc<3;cc++){ + int out_idx = cc+3*(sc + 4*fc); + + for(int fr=0;fr<2;fr++){ + peekSpinColorColumn(tmp1f, *Minv[fr][fc] , sc, cc); + + PokeIndex(out[out_idx], tmp1f, fr); + } + } + } + } +} + + + + +void invertXconj1d(std::vector &out, XconjugateMobiusFermionD &action){ + GridBase* UGrid = action.GaugeGrid(); + GridBase* FGrid = action.FermionGrid(); + GridBase* FrbGrid = action.FermionRedBlackGrid(); + + ConjugateGradient cg(1e-08,10000); + SchurRedBlackDiagMooeeSolve solver(cg); + + typedef FermionField1f::scalar_object Spinor; + + LatticeInteger tcoor4d(UGrid); + LatticeCoordinate(tcoor4d,3); + + FermionField1f tmp4d(UGrid); + FermionField1f src4d(UGrid); + FermionField1f src5d(FGrid), sol5d(FGrid); + FermionField1f src5d_e(FrbGrid), src5d_o(FrbGrid), sol5d_o(FrbGrid); + + SCmatrixField tmp_vplus_f0_4d(UGrid); + SCmatrixField tmp_vminus_f0_4d(UGrid); + + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + + for(int s=0;s<4;s++){ + for(int c=0;c<3;c++){ + Spinor vb = Zero(); + vb()(s)(c) = 1.; + + Spinor vplus = mulPplusLeft(vb); + Spinor vminus = mulPminusLeft(vb); + + Spinor* vpm[2] = {&vplus, &vminus}; + SCmatrixField* solpm[2] = {&tmp_vplus_f0_4d, &tmp_vminus_f0_4d}; + + for(int pm=0;pm<2;pm++){ + tmp4d = *vpm[pm]; + src4d = Zero(); + src4d = where( tcoor4d == Integer(0), tmp4d, src4d); //vector on every site on timeslice 0, zero elsewhere + + action.ImportPhysicalFermionSource(src4d, src5d); + solver.RedBlackSource(action, src5d, src5d_e, src5d_o); + solver.RedBlackSolve(action, src5d_o, sol5d_o); + solver.RedBlackSolution(action, sol5d_o, src5d_e, sol5d); + action.ExportPhysicalFermionSolution(sol5d, tmp4d); + pokeSpinColorColumn(*solpm[pm], tmp4d, s,c); + } + } + } + + SCmatrixField vplus_f0_Pplus = mulPplusRight(tmp_vplus_f0_4d); + SCmatrixField vplus_f0_Pminus = mulPminusRight(tmp_vplus_f0_4d); + + SCmatrixField vminus_f0_Pplus = mulPplusRight(tmp_vminus_f0_4d); + SCmatrixField vminus_f0_Pminus = mulPminusRight(tmp_vminus_f0_4d); + + SCmatrixField Minv11 = vplus_f0_Pplus + vminus_f0_Pminus; + SCmatrixField Minv12 = vplus_f0_Pminus*X + vminus_f0_Pplus*X; + SCmatrixField Minv22 = -(X*(conjugate(Minv11)*X)); + SCmatrixField Minv21 = X*(conjugate(Minv12)*X); + + out = std::vector(24, UGrid); + + SCmatrixField* Minv[2][2] = { {&Minv11,&Minv12}, {&Minv21,&Minv22} }; + + for(int fc=0;fc<2;fc++){ + for(int sc=0;sc<4;sc++){ + for(int cc=0;cc<3;cc++){ + int out_idx = cc+3*(sc + 4*fc); + + for(int fr=0;fr<2;fr++){ + peekSpinColorColumn(tmp4d, *Minv[fr][fc] , sc, cc); + + PokeIndex(out[out_idx], tmp4d, fr); + } + } + } + } +} + + +template +T mulURight(const T &in){ + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static GparityFlavour sigma3 = GparityFlavour(GparityFlavour::Algebra::SigmaZ); + + T out = 0.5*(in + in*X) + 0.5*(in*sigma3 - (in*X)*sigma3); + return out; +} + + + +void invertXconj1d_matrix(std::vector &out, XconjugateMobiusFermionD &action){ + GridBase* UGrid = action.GaugeGrid(); + GridBase* FGrid = action.FermionGrid(); + GridBase* FrbGrid = action.FermionRedBlackGrid(); + + ConjugateGradient cg(1e-08,10000); + SchurRedBlackDiagMooeeSolve solver(cg); + + typedef FermionField1f::scalar_object Spinor; + + LatticeInteger tcoor4d(UGrid); + LatticeCoordinate(tcoor4d,3); + + FermionField1f tmp4d(UGrid); + FermionField1f src4d(UGrid); + FermionField1f src5d(FGrid), sol5d(FGrid); + FermionField1f src5d_e(FrbGrid), src5d_o(FrbGrid), sol5d_o(FrbGrid); + + FermionField2f tmp2f(UGrid); + + SCFmatrixField V_4d(UGrid); + + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static GparityFlavour sigma1 = GparityFlavour(GparityFlavour::Algebra::SigmaX); + + + for(int s=0;s<4;s++){ + for(int c=0;c<3;c++){ + Spinor vb = Zero(); + vb()(s)(c) = 1.; + + Spinor vplus = mulPplusLeft(vb); + Spinor vminus = mulPminusLeft(vb); + + Spinor* vpm[2] = {&vplus, &vminus}; + + for(int pm=0;pm<2;pm++){ + tmp4d = *vpm[pm]; + src4d = Zero(); + src4d = where( tcoor4d == Integer(0), tmp4d, src4d); //vector on every site on timeslice 0, zero elsewhere + + action.ImportPhysicalFermionSource(src4d, src5d); + solver.RedBlackSource(action, src5d, src5d_e, src5d_o); + solver.RedBlackSolve(action, src5d_o, sol5d_o); + solver.RedBlackSolution(action, sol5d_o, src5d_e, sol5d); + action.ExportPhysicalFermionSolution(sol5d, tmp4d); + + //Generate 2f X-conjugate output + PokeIndex(tmp2f, tmp4d, 0); + tmp4d = -(X*conjugate(tmp4d)); + PokeIndex(tmp2f, tmp4d, 1); + + pokeSpinColorFlavorColumn(V_4d, tmp2f, s,c,pm); + } + } + } + + SCFmatrixField tmp = mulPplusRight(V_4d) + mulPminusRight(V_4d)*sigma1; + SCFmatrixField Minv = mulURight(tmp); + + out = std::vector(24, UGrid); + + for(int fc=0;fc<2;fc++){ + for(int sc=0;sc<4;sc++){ + for(int cc=0;cc<3;cc++){ + int out_idx = cc+3*(sc + 4*fc); + peekSpinColorFlavorColumn(out[out_idx], Minv, sc, cc, fc); + } + } + } +} + + +//Use the full 2x2 G-parity Dirac op inverted on a complex source with flavor structure +// |\phi 0 | +// | 0 \phi^* | +void invertGparityComplexSrc(std::vector &out, GparityMobiusFermionD &action, ComplexD phi){ + GridBase* UGrid = action.GaugeGrid(); + GridBase* FGrid = action.FermionGrid(); + GridBase* FrbGrid = action.FermionRedBlackGrid(); + + ConjugateGradient cg(1e-08,10000); + SchurRedBlackDiagMooeeSolve solver(cg); + + typedef FermionField2f::scalar_object Spinor; + + LatticeInteger tcoor4d(UGrid); + LatticeCoordinate(tcoor4d,3); + + FermionField2f tmp4d(UGrid); + FermionField2f src4d(UGrid); + FermionField2f src5d(FGrid), sol5d(FGrid); + FermionField2f src5d_e(FrbGrid), src5d_o(FrbGrid), sol5d_o(FrbGrid); + + out = std::vector(24, UGrid); + + for(int f=0;f<2;f++){ + for(int s=0;s<4;s++){ + for(int c=0;c<3;c++){ + Spinor v = Zero(); + v(f)(s)(c) = f==0 ? phi : conjugate(phi); + + tmp4d = v; + src4d = Zero(); + src4d = where( tcoor4d == Integer(0), tmp4d, src4d); //unit vector on every site on timeslice 0, zero elsewhere + + action.ImportPhysicalFermionSource(src4d, src5d); + solver.RedBlackSource(action, src5d, src5d_e, src5d_o); + solver.RedBlackSolve(action, src5d_o, sol5d_o); + solver.RedBlackSolution(action, sol5d_o, src5d_e, sol5d); + action.ExportPhysicalFermionSolution(sol5d, out[c+3*(s + 4*f)]); + } + } + } +} + +//Do the same as the above but with the X-conjugate Dirac op +void invertXconj1dComplexSrc(std::vector &out, XconjugateMobiusFermionD &action, ComplexD phi){ + GridBase* UGrid = action.GaugeGrid(); + GridBase* FGrid = action.FermionGrid(); + GridBase* FrbGrid = action.FermionRedBlackGrid(); + + ConjugateGradient cg(1e-08,10000); + SchurRedBlackDiagMooeeSolve solver(cg); + + typedef FermionField1f::scalar_object Spinor; + + LatticeInteger tcoor4d(UGrid); + LatticeCoordinate(tcoor4d,3); + + FermionField1f tmp4d(UGrid); + FermionField1f src4d(UGrid); + FermionField1f src5d(FGrid), sol5d(FGrid); + FermionField1f src5d_e(FrbGrid), src5d_o(FrbGrid), sol5d_o(FrbGrid); + + FermionField2f tmp2f(UGrid); + + SCFmatrixField V_4d(UGrid); + + static Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + static Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + static Gamma X = C*g5; + static GparityFlavour sigma1 = GparityFlavour(GparityFlavour::Algebra::SigmaX); + + + for(int s=0;s<4;s++){ + for(int c=0;c<3;c++){ + Spinor vb = Zero(); + vb()(s)(c) = phi; + + Spinor vplus = mulPplusLeft(vb); + Spinor vminus = mulPminusLeft(vb); + + Spinor* vpm[2] = {&vplus, &vminus}; + + for(int pm=0;pm<2;pm++){ + tmp4d = *vpm[pm]; + src4d = Zero(); + src4d = where( tcoor4d == Integer(0), tmp4d, src4d); //vector on every site on timeslice 0, zero elsewhere + + action.ImportPhysicalFermionSource(src4d, src5d); + solver.RedBlackSource(action, src5d, src5d_e, src5d_o); + solver.RedBlackSolve(action, src5d_o, sol5d_o); + solver.RedBlackSolution(action, sol5d_o, src5d_e, sol5d); + action.ExportPhysicalFermionSolution(sol5d, tmp4d); + + //Generate 2f X-conjugate output + PokeIndex(tmp2f, tmp4d, 0); + tmp4d = -(X*conjugate(tmp4d)); + PokeIndex(tmp2f, tmp4d, 1); + + pokeSpinColorFlavorColumn(V_4d, tmp2f, s,c,pm); + } + } + } + + SCFmatrixField tmp = mulPplusRight(V_4d) + mulPminusRight(V_4d)*sigma1; + SCFmatrixField Minv = mulURight(tmp); + + out = std::vector(24, UGrid); + + for(int fc=0;fc<2;fc++){ + for(int sc=0;sc<4;sc++){ + for(int cc=0;cc<3;cc++){ + int out_idx = cc+3*(sc + 4*fc); + peekSpinColorFlavorColumn(out[out_idx], Minv, sc, cc, fc); + } + } + } +} + + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + RealD tol = 1e-10; + for(int i=1;i> tol; + std::cout << "Set tolerance to " << tol << std::endl; + } + } + + const int Ls=4; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeGaugeField Umu(UGrid); + SU::HotConfiguration(RNG4, Umu); + + Gamma C = Gamma(Gamma::Algebra::MinusGammaY) * Gamma(Gamma::Algebra::GammaT); + Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + Gamma X = C*g5; + + //Set up a regular MDWF action instance as well as X-conj and Xbar-conj versions + RealD mass=0.04; + RealD M5=1.8; + RealD mob_b=1.5; + GparityMobiusFermionD ::ImplParams params; + std::vector twists({1,1,1,0}); + params.twists = twists; + + GparityMobiusFermionD reg_action(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); + + XconjugateMobiusFermionD::ImplParams xparams; + xparams.twists = twists; + xparams.boundary_phase = 1.0; + + XconjugateMobiusFermionD xconj_action(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,xparams); + + xparams.boundary_phase = -1.0; + XconjugateMobiusFermionD xbarconj_action(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,xparams); + + FermionField2f tmp2f(UGrid); + std::cout.precision(12); + + { + std::cout << "Tests with real source" << std::endl; + std::cout << "-------------------------------------------------" << std::endl; + + std::vector sol_gp; + invertGparity(sol_gp, reg_action); + + std::vector sol_Xconj2d; + invertXconj2d(sol_Xconj2d, reg_action); + + std::vector sol_Xconj1d; + invertXconj1d(sol_Xconj1d, xconj_action); + + std::vector sol_Xconj1d_matrix; + invertXconj1d_matrix(sol_Xconj1d_matrix, xconj_action); + + std::cout << "Comparisons for real source" << std::endl; + for(int fc=0;fc<2;fc++){ + for(int sc=0;sc<4;sc++){ + for(int cc=0;cc<3;cc++){ + int out_idx = cc+3*(sc + 4*fc); + tmp2f = sol_gp[out_idx] - sol_Xconj2d[out_idx]; + RealD nrm2d = norm2(tmp2f); + tmp2f = sol_gp[out_idx] - sol_Xconj1d[out_idx]; + RealD nrm1d = norm2(tmp2f); + tmp2f = sol_gp[out_idx] - sol_Xconj1d_matrix[out_idx]; + RealD nrm1dm = norm2(tmp2f); + + std::cout << fc << " " << sc << " " << cc << " " << nrm2d << " " << nrm1d << " " << nrm1dm << std::endl; + assert(nrm2d < tol); + assert(nrm1d < tol); + assert(nrm1dm < tol); + } + } + } + } + + { + std::cout << "Tests with complex source" << std::endl; + std::cout << "-------------------------------------------------" << std::endl; + + ComplexD phi(1.234, -6.444); //"random" phase + + std::vector sol_gp; + invertGparityComplexSrc(sol_gp, reg_action, phi); + + std::vector sol_Xconj1d_matrix; + invertXconj1dComplexSrc(sol_Xconj1d_matrix, xconj_action, phi); + + std::cout << "Comparisons for complex source" << std::endl; + for(int fc=0;fc<2;fc++){ + for(int sc=0;sc<4;sc++){ + for(int cc=0;cc<3;cc++){ + int out_idx = cc+3*(sc + 4*fc); + tmp2f = sol_gp[out_idx] - sol_Xconj1d_matrix[out_idx]; + RealD nrm1dm = norm2(tmp2f); + + std::cout << fc << " " << sc << " " << cc << " " << nrm1dm << std::endl; + assert(nrm1dm < tol); + } + } + } + } + + std::cout << "All tests passed" << std::endl; + + Grid_finalize(); +}