portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00
Code Releases Activity

Mobius Caley form, Mobius Zolotarev operators. Pass Even Odd vs unprec test and hermiticity checks

Browse Source 
in tests/Grid_any_evenodd.cc; will work on inversion tests shortly.
This commit is contained in:
Peter Boyle
2015-06-03 09:36:26 +01:00
parent 494d2b8b61
commit 68e26140ee
10 changed files with 507 additions and 95 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
lib/qcd/action/Actions.h
View File

@ -38,6 +38,9 @@
//#include <qcd/action/fermion/PartialFraction.h>
#include <qcd/action/fermion/DomainWallFermion.h>
#include <qcd/action/fermion/DomainWallFermion.h>
#include <qcd/action/fermion/MobiusFermion.h>
#include <qcd/action/fermion/MobiusZolotarevFermion.h>
//#include <qcd/action/fermion/ScaledShamirCayleyTanh.h>
@ -70,11 +73,6 @@
class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
typedef multi1d<LatticeColorMatrix> P;
typedef multi1d<LatticeColorMatrix> Q;
virtual void staple(LatticeColorMatrix& result,
const Handle< GaugeState<P,Q> >& state,
int mu, int cb) const = 0;
*/
#endif

110
lib/qcd/action/fermion/CayleyFermion5D.cc
View File

@ -15,7 +15,6 @@ namespace QCD {
FourDimRedBlackGrid,_M5),
mass(_mass)
{
std::cout << "Constructing a CayleyFermion5D"<<std::endl;
}
// override multiply
@ -229,7 +228,112 @@ namespace QCD {
axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1); // chi[Ls]
}
}
void CayleyFermion5D::SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c)
{
///////////////////////////////////////////////////////////
// The Cayley coeffs (unprec)
///////////////////////////////////////////////////////////
omega.resize(Ls);
bs.resize(Ls);
cs.resize(Ls);
as.resize(Ls);
//
// Ts = ( [bs+cs]Dw )^-1 ( (bs+cs) Dw )
// -(g5 ------- -1 ) ( g5 --------- + 1 )
// ( {2+(bs-cs)Dw} ) ( 2+(bs-cs) Dw )
//
// bs = 1/2( (1/omega_s + 1)*b + (1/omega - 1)*c ) = 1/2( 1/omega(b+c) + (b-c) )
// cs = 1/2( (1/omega_s - 1)*b + (1/omega + 1)*c ) = 1/2( 1/omega(b+c) - (b-c) )
//
// bs+cs = 0.5*( 1/omega(b+c) + (b-c) + 1/omega(b+c) - (b-c) ) = 1/omega(b+c)
// bs-cs = 0.5*( 1/omega(b+c) + (b-c) - 1/omega(b+c) + (b-c) ) = b-c
//
// So
//
// Ts = ( [b+c]Dw/omega_s )^-1 ( (b+c) Dw /omega_s )
// -(g5 ------- -1 ) ( g5 --------- + 1 )
// ( {2+(b-c)Dw} ) ( 2+(b-c) Dw )
//
// Ts = ( [b+c]Dw )^-1 ( (b+c) Dw )
// -(g5 ------- -omega_s) ( g5 --------- + omega_s )
// ( {2+(b-c)Dw} ) ( 2+(b-c) Dw )
//
double bpc = b+c;
double bmc = b-c;
for(int i=0; i < Ls; i++){
as[i] = 1.0;
omega[i] = ((double)zdata->gamma[i]); //NB reciprocal relative to Chroma NEF code
bs[i] = 0.5*(bpc/omega[i] + bmc);
cs[i] = 0.5*(bpc/omega[i] - bmc);
}
////////////////////////////////////////////////////////
// Constants for the preconditioned matrix Cayley form
////////////////////////////////////////////////////////
bee.resize(Ls);
cee.resize(Ls);
beo.resize(Ls);
ceo.resize(Ls);
for(int i=0;i<Ls;i++){
bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0);
cee[i]=as[i]*(1.0-cs[i]*(4.0-M5));
beo[i]=as[i]*bs[i];
ceo[i]=-as[i]*cs[i];
}
aee.resize(Ls);
aeo.resize(Ls);
for(int i=0;i<Ls;i++){
aee[i]=cee[i];
aeo[i]=ceo[i];
}
//////////////////////////////////////////
// LDU decomposition of eeoo
//////////////////////////////////////////
dee.resize(Ls);
lee.resize(Ls);
leem.resize(Ls);
uee.resize(Ls);
ueem.resize(Ls);
for(int i=0;i<Ls;i++){
dee[i] = bee[i];
if ( i < Ls-1 ) {
lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
leem[i]=mass*cee[Ls-1]/bee[0];
for(int j=0;j<i;j++) leem[i]*= aee[j]/bee[j+1];
uee[i] =-aee[i]/bee[i]; // up-diag entry on the ith row
ueem[i]=mass;
for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
ueem[i]*= aee[0]/bee[0];
} else {
lee[i] =0.0;
leem[i]=0.0;
uee[i] =0.0;
ueem[i]=0.0;
}
}
{
double delta_d=mass*cee[Ls-1];
for(int j=0;j<Ls-1;j++) delta_d *= cee[j]/bee[j];
dee[Ls-1] += delta_d;
}
}
}}
}
}

6
lib/qcd/action/fermion/CayleyFermion5D.h
View File

@ -22,10 +22,8 @@ namespace Grid {
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
// protected:
Approx::zolotarev_data *zdata;
RealD mass;
// Cayley form Moebius (tanh and zolotarev)
std::vector<RealD> omega;
std::vector<RealD> bs; // S dependent coeffs
@ -53,6 +51,8 @@ namespace Grid {
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5);
protected:
void SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c);
};
}

83
lib/qcd/action/fermion/DomainWallFermion.h
View File

@ -28,86 +28,13 @@ namespace Grid {
{
RealD eps = 1.0;
zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficients(1.0,zdata,1.0,0.0);
///////////////////////////////////////////////////////////
// The Cayley coeffs (unprec)
///////////////////////////////////////////////////////////
this->omega.resize(this->Ls);
this->bs.resize(this->Ls);
this->cs.resize(this->Ls);
this->as.resize(this->Ls);
for(int i=0; i < this->Ls; i++){
this->as[i] = 1.0;
this->omega[i] = ((double)zdata -> gamma[i]);
double bb=1.0;
this->bs[i] = 0.5*(bb/(this->omega[i]) + 1.0);
this->cs[i] = 0.5*(bb/(this->omega[i]) - 1.0);
}
////////////////////////////////////////////////////////
// Constants for the preconditioned matrix Cayley form
////////////////////////////////////////////////////////
this->bee.resize(this->Ls);
this->cee.resize(this->Ls);
this->beo.resize(this->Ls);
this->ceo.resize(this->Ls);
for(int i=0;i<this->Ls;i++){
this->bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0);
this->cee[i]=as[i]*(1.0-cs[i]*(4.0-M5));
this->beo[i]=as[i]*bs[i];
this->ceo[i]=-as[i]*cs[i];
}
aee.resize(this->Ls);
aeo.resize(this->Ls);
for(int i=0;i<this->Ls;i++){
aee[i]=cee[i];
aeo[i]=ceo[i];
}
//////////////////////////////////////////
// LDU decomposition of eeoo
//////////////////////////////////////////
dee.resize(this->Ls);
lee.resize(this->Ls);
leem.resize(this->Ls);
uee.resize(this->Ls);
ueem.resize(this->Ls);
for(int i=0;i<this->Ls;i++){
dee[i] = bee[i];
if ( i < this->Ls-1 ) {
lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
leem[i]=this->mass*cee[this->Ls-1]/bee[0];
for(int j=0;j<i;j++) leem[i]*= aee[j]/bee[j+1];
uee[i] =-aee[i]/bee[i]; // up-diag entry on the ith row
ueem[i]=this->mass;
for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
ueem[i]*= aee[0]/bee[0];
} else {
lee[i] =0.0;
leem[i]=0.0;
uee[i] =0.0;
ueem[i]=0.0;
}
}
{
double delta_d=mass*cee[this->Ls-1];
for(int j=0;j<this->Ls-1;j++) delta_d *= cee[j]/bee[j];
dee[this->Ls-1] += delta_d;
}
}
};

46
lib/qcd/action/fermion/MobiusFermion.h Normal file
View File

@ -0,0 +1,46 @@
#ifndef GRID_QCD_MOBIUS_FERMION_H
#define GRID_QCD_MOBIUS_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class MobiusFermion : public CayleyFermion5D
{
public:
// Constructors
MobiusFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD b, RealD c) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
RealD eps = 1.0;
std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
// Call base setter
this->CayleyFermion5D::SetCoefficients(1.0,zdata,b,c);
}
};
}
}
#endif

48
lib/qcd/action/fermion/MobiusZolotarevFermion.h Normal file
View File

@ -0,0 +1,48 @@
#ifndef GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
#define GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class MobiusZolotarevFermion : public CayleyFermion5D
{
public:
// Constructors
MobiusZolotarevFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5,
RealD b, RealD c,
RealD lo, RealD hi) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
RealD eps = lo/hi;
Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,this->Ls,0);// eps is ignored for higham
assert(zdata->n==this->Ls);
std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficients(1.0,zdata,b,c);
}
};
}
}
#endif

51
lib/qcd/action/fermion/ScaledShamir.h Normal file
View File

@ -0,0 +1,51 @@
#ifndef GRID_QCD_DOMAIN_WALL_FERMION_H
#define GRID_QCD_DOMAIN_WALL_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class ScaledShamirFermion : public CayleyFermion5D
{
public:
// Constructors
ScaledShamirFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5, RealD scale) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5,
RealD b,
RealD c)
{
RealD eps = 1.0;
Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
//b+c = scale;
//b-c = 1
//b = 0.5(scale+1);
//c = 0.5(scale-1);
// Call base setter
this->CayleyFermion5D::SetCoefficients(1.0,zdata,0.5*(scale+1.0),0.5*(scale-1.0));
}
};
}
}
#endif

12
lib/qcd/action/fermion/WilsonFermion5D.h
View File

@ -14,6 +14,18 @@ namespace Grid {
// i.e. even even contains fifth dim hopping term.
//
// [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ]
////////////////////////////
//ContFrac:
// Ls always odd. Rational poly deg is either Ls or Ls-1
//PartFrac
// Ls always odd. Rational poly deg is either Ls or Ls-1
//
//Cayley: Ls always even, Rational poly deg is Ls
//
// Just set nrational as Ls. Forget about Ls-1 cases.
//
// Require odd Ls for cont and part frac
////////////////////////////
////////////////////////////////////////////////////////////////////////////////
class WilsonFermion5D : public FermionOperator<LatticeFermion,LatticeGaugeField>
{