mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-13 04:37:05 +01:00
Mobius Caley form, Mobius Zolotarev operators. Pass Even Odd vs unprec test and hermiticity checks
in tests/Grid_any_evenodd.cc; will work on inversion tests shortly.
This commit is contained in:
Peter Boyle
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
}}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -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);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -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
46
lib/qcd/action/fermion/MobiusFermion.h
Normal 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
48
lib/qcd/action/fermion/MobiusZolotarevFermion.h
Normal 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
51
lib/qcd/action/fermion/ScaledShamir.h
Normal 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
|
||||
@ -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>
|
||||
{
|
||||
|
||||
Reference in New Issue
Block a user