diff --git a/TODO b/TODO index 064e99b4..46c12838 100644 --- a/TODO +++ b/TODO @@ -1,5 +1,6 @@ * - BinaryWriter, TextWriter etc... - - protocol buffers? replace xmlReader/Writer ec.. + - use protocol buffers? replace xmlReader/Writer ec.. + - Binary use htonll, htonl * Stencil operator support -----Initial thoughts, trial implementation DONE. -----some simple tests that Stencil matches Cshift. @@ -15,10 +16,9 @@ boost::multi_array A()... to replace multi1d, multi2d etc.. -* norm2l is a hack. figure out syntax error and make this norm2. c.f. tests/Grid_gamma.cc +* norm2l is a hack. figure out syntax error and make this norm2 c.f. tests/Grid_gamma.cc - -* std::vector replacement; +**** std::vector replacement; Had to change "reserve" back to "resize" on std::vector in Lattice class. @@ -59,7 +59,6 @@ * rb4d support for 5th dimension in Mobius. - * Check for missing functionality - partially audited against QDP++ layout Algorithms @@ -70,6 +69,7 @@ Algorithms * Pcg * fPcg * MCR +* HDCG * etc.. AUDITS: diff --git a/lib/qcd/Grid_qcd_wilson_dop.cc b/lib/qcd/Grid_qcd_wilson_dop.cc new file mode 100644 index 00000000..73b19c12 --- /dev/null +++ b/lib/qcd/Grid_qcd_wilson_dop.cc @@ -0,0 +1,157 @@ +#ifnfdef GRID_QCD_WILSON_DOP_H +#define GRID_QCD_WILSON_DOP_H + +#include + +namespace Grid { +namespace QCD { + + +const std::vector WilsonMatrix::directions ({0,1,2,3, 0, 1, 2, 3,0}); +const std::vector WilsonMatrix::displacements({1,1,1,1,-1,-1,-1,-1,0}); + + // Should be in header? +static const int WilsonMatrix::Xp = 0; +static const int WilsonMatrix::Yp = 1; +static const int WilsonMatrix::Zp = 2; +static const int WilsonMatrix::Tp = 3; +static const int WilsonMatrix::Xm = 4; +static const int WilsonMatrix::Ym = 5; +static const int WilsonMatrix::Zm = 6; +static const int WilsonMatrix::Tm = 7; +static const int WilsonMatrix::X0 = 8; +static const int WilsonMatrix::npoint=9; + + +WilsonMatrix::WilsonMatrix(LatticeGaugeField &_Umu,int _mass) + : Stencil((&Umu._grid,npoint,0,directions,displacements), + mass(_mass), + Umu(_Umu) +{ + // Allocate the required comms buffer + grid = _Umu._grid; + comm_buf.resize(Stencil._unified_buffer_size); +} +void WilsonMatrix::multiply(const LatticeFermion &in, LatticeFermion &out) +{ + +} +void WilsonMatrix::Dhop(const LatticeFermion &in, LatticeFermion &out) +{ + Stencil.HaloExchange(in,comm_buf); + + for(int ss=0;ss<_grid->oSites();ss++){ + + int offset,local; + + vSpinColourVector result; + vHalfSpinColourVector UChi; + + // Xp + offset = Stencil._offsets [Xp][ss]; + local = Stencil._is_local[Xp][ss]; + if ( local ) { + Uchi = U[]*spProjXp(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result = ReconXp(Uchi); + + // Yp + offset = Stencil._offsets [Yp][ss]; + local = Stencil._is_local[Yp][ss]; + if ( local ) { + Uchi = U[]*spProjYp(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconYp(Uchi); + + // Zp + offset = Stencil._offsets [Zp][ss]; + local = Stencil._is_local[Zp][ss]; + if ( local ) { + Uchi = U[]*spProjZp(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconZp(Uchi); + + // Tp + offset = Stencil._offsets [Tp][ss]; + local = Stencil._is_local[Tp][ss]; + if ( local ) { + Uchi = U[]*spProjTp(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconTp(Uchi); + + // Xm + offset = Stencil._offsets [Xm][ss]; + local = Stencil._is_local[Xm][ss]; + if ( local ) { + Uchi = U[]*spProjXm(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconXm(Uchi); + + // Ym + offset = Stencil._offsets [Ym][ss]; + local = Stencil._is_local[Ym][ss]; + if ( local ) { + Uchi = U[]*spProjYm(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconYm(Uchi); + + // Zm + offset = Stencil._offsets [Zm][ss]; + local = Stencil._is_local[Zm][ss]; + if ( local ) { + Uchi = U[]*spProjZm(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconZm(Uchi); + + // Tm + offset = Stencil._offsets [Tm][ss]; + local = Stencil._is_local[Tm][ss]; + if ( local ) { + Uchi = U[]*spProjTm(in._odata[offset]); + } else { + Uchi = U[]*comm_buf._odata[offset] + } + result+= ReconTm(Uchi); + + out._odata[ss] = result; + } + +} +void WilsonMatrix::Dw(const LatticeFermion &in, LatticeFermion &out) +{ + +} +void WilsonMatrix::MpcDag (const LatticeFermion &in, LatticeFermion &out) +{ + +} +void WilsonMatrix::Mpc (const LatticeFermion &in, LatticeFermion &out) +{ + +} +void WilsonMatrix::MpcDagMpc(const LatticeFermion &in, LatticeFermion &out) +{ + +} +void WilsonMatrix::MDagM (const LatticeFermion &in, LatticeFermion &out) +{ + +} + + +}} +#endif diff --git a/lib/qcd/Grid_qcd_wilson_dop.h b/lib/qcd/Grid_qcd_wilson_dop.h new file mode 100644 index 00000000..f6775842 --- /dev/null +++ b/lib/qcd/Grid_qcd_wilson_dop.h @@ -0,0 +1,62 @@ +#ifnfdef GRID_QCD_WILSON_DOP_H +#define GRID_QCD_WILSON_DOP_H + +#include + +namespace Grid { + + namespace QCD { + + + template class LinearOperatorBase { + public: + void multiply(const Lattice &in, Lattice &out){ assert(0);} + }; + + class WilsonMatrix : public LinearOperatorBase + { + //NB r=1; + public: + double mass; + GridBase *grid; + + // Copy of the gauge field + LatticeGaugeField Umu; + + //Defines the stencil + CartesianStencil Stencil; + static const int npoint=9; + static const std::vector directions ; + static const std::vector displacements; + + static const int Xp,Xm,Yp,Ym,Zp,Zm,Tp,Tm; + + // Comms buffer + std::vector > comm_buf; + + // Constructor + WilsonMatrix(LatticeGaugeField &Umu,int mass); + + // override multiply + void multiply(const LatticeFermion &in, LatticeFermion &out); + + // non-hermitian hopping term; half cb or both + void Dhop(const LatticeFermion &in, LatticeFermion &out); + + // m+4r -1/2 Dhop; both cb's + void Dw(const LatticeFermion &in, LatticeFermion &out); + + // half checkerboard operaions + void MpcDag (const LatticeFermion &in, LatticeFermion &out); + void Mpc (const LatticeFermion &in, LatticeFermion &out); + void MpcDagMpc(const LatticeFermion &in, LatticeFermion &out); + + // full checkerboard hermitian + void MDagM (const LatticeFermion &in, LatticeFermion &out); + + + }; + + } +} +#endif