portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-18 07:47: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 35fdba81dd
commit 5916386242
10 changed files with 507 additions and 95 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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;
}
}
};