1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00

Starting the implementation of wilson; incomplete and committing non-functional code which

is not yet included from elsewhere or linked to the build system.
This commit is contained in:
Peter Boyle 2015-04-25 14:33:02 +01:00
parent 9dacdc947d
commit 51f0da7b93
3 changed files with 224 additions and 5 deletions

10
TODO
View File

@ -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<type, 3> 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:

View File

@ -0,0 +1,157 @@
#ifnfdef GRID_QCD_WILSON_DOP_H
#define GRID_QCD_WILSON_DOP_H
#include <Grid.h>
namespace Grid {
namespace QCD {
const std::vector<int> WilsonMatrix::directions ({0,1,2,3, 0, 1, 2, 3,0});
const std::vector<int> 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

View File

@ -0,0 +1,62 @@
#ifnfdef GRID_QCD_WILSON_DOP_H
#define GRID_QCD_WILSON_DOP_H
#include <Grid.h>
namespace Grid {
namespace QCD {
template<class vtype> class LinearOperatorBase {
public:
void multiply(const Lattice<vtype> &in, Lattice<vtype> &out){ assert(0);}
};
class WilsonMatrix : public LinearOperatorBase<LatticeFermion>
{
//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<int> directions ;
static const std::vector<int> displacements;
static const int Xp,Xm,Yp,Ym,Zp,Zm,Tp,Tm;
// Comms buffer
std::vector<vSpinColourVector,alignedAllocator<vSpinColourVector> > 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