mirror of
https://github.com/paboyle/Grid.git
synced 2024-09-20 09:15:38 +01:00
Namespace, indent
This commit is contained in:
paboyle
parent
1619e42d90
commit
83c5f05094
@ -28,411 +28,410 @@ with this program; if not, write to the Free Software Foundation, Inc.,
|
|||||||
|
|
||||||
See the full license in the file "LICENSE" in the top level distribution directory
|
See the full license in the file "LICENSE" in the top level distribution directory
|
||||||
*************************************************************************************/
|
*************************************************************************************/
|
||||||
/* END LEGAL */
|
/* END LEGAL */
|
||||||
|
|
||||||
#include <Grid/Grid_Eigen_Dense.h>
|
#include <Grid/Grid_Eigen_Dense.h>
|
||||||
#include <Grid/qcd/action/fermion/FermionCore.h>
|
#include <Grid/qcd/action/fermion/FermionCore.h>
|
||||||
#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
|
#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
|
||||||
|
|
||||||
namespace Grid {
|
NAMESPACE_BEGIN(Grid);
|
||||||
namespace QCD {
|
|
||||||
|
|
||||||
template<class Impl>
|
template<class Impl>
|
||||||
DomainWallEOFAFermion<Impl>::DomainWallEOFAFermion(
|
DomainWallEOFAFermion<Impl>::DomainWallEOFAFermion(
|
||||||
GaugeField &_Umu,
|
GaugeField &_Umu,
|
||||||
GridCartesian &FiveDimGrid,
|
GridCartesian &FiveDimGrid,
|
||||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||||
GridCartesian &FourDimGrid,
|
GridCartesian &FourDimGrid,
|
||||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||||
RealD _mq1, RealD _mq2, RealD _mq3,
|
RealD _mq1, RealD _mq2, RealD _mq3,
|
||||||
RealD _shift, int _pm, RealD _M5, const ImplParams &p) :
|
RealD _shift, int _pm, RealD _M5, const ImplParams &p) :
|
||||||
AbstractEOFAFermion<Impl>(_Umu, FiveDimGrid, FiveDimRedBlackGrid,
|
AbstractEOFAFermion<Impl>(_Umu, FiveDimGrid, FiveDimRedBlackGrid,
|
||||||
FourDimGrid, FourDimRedBlackGrid, _mq1, _mq2, _mq3,
|
FourDimGrid, FourDimRedBlackGrid, _mq1, _mq2, _mq3,
|
||||||
_shift, _pm, _M5, 1.0, 0.0, p)
|
_shift, _pm, _M5, 1.0, 0.0, p)
|
||||||
{
|
{
|
||||||
RealD eps = 1.0;
|
RealD eps = 1.0;
|
||||||
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);
|
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);
|
||||||
assert(zdata->n == this->Ls);
|
assert(zdata->n == this->Ls);
|
||||||
|
|
||||||
std::cout << GridLogMessage << "DomainWallEOFAFermion with Ls=" << this->Ls << std::endl;
|
std::cout << GridLogMessage << "DomainWallEOFAFermion with Ls=" << this->Ls << std::endl;
|
||||||
this->SetCoefficientsTanh(zdata, 1.0, 0.0);
|
this->SetCoefficientsTanh(zdata, 1.0, 0.0);
|
||||||
|
|
||||||
|
Approx::zolotarev_free(zdata);
|
||||||
|
}
|
||||||
|
|
||||||
|
/***************************************************************
|
||||||
|
* Additional EOFA operators only called outside the inverter.
|
||||||
|
* Since speed is not essential, simple axpby-style
|
||||||
|
* implementations should be fine.
|
||||||
|
***************************************************************/
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& Din, int sign, int dag)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
|
||||||
|
Din = zero;
|
||||||
|
if((sign == 1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, Ls-1, 0); }
|
||||||
|
else if((sign == -1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); }
|
||||||
|
else if((sign == 1 ) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, Ls-1); }
|
||||||
|
else if((sign == -1) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); }
|
||||||
|
}
|
||||||
|
|
||||||
|
// This is just the identity for DWF
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::Dtilde(const FermionField& psi, FermionField& chi){ chi = psi; }
|
||||||
|
|
||||||
|
// This is just the identity for DWF
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionField& chi){ chi = psi; }
|
||||||
|
|
||||||
|
/*****************************************************************************************************/
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
|
||||||
|
FermionField Din(psi._grid);
|
||||||
|
|
||||||
|
this->Meooe5D(psi, Din);
|
||||||
|
this->DW(Din, chi, DaggerNo);
|
||||||
|
axpby(chi, 1.0, 1.0, chi, psi);
|
||||||
|
this->M5D(psi, chi);
|
||||||
|
return(norm2(chi));
|
||||||
|
}
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
|
||||||
|
FermionField Din(psi._grid);
|
||||||
|
|
||||||
|
this->DW(psi, Din, DaggerYes);
|
||||||
|
this->MeooeDag5D(Din, chi);
|
||||||
|
this->M5Ddag(psi, chi);
|
||||||
|
axpby(chi, 1.0, 1.0, chi, psi);
|
||||||
|
return(norm2(chi));
|
||||||
|
}
|
||||||
|
|
||||||
|
/********************************************************************
|
||||||
|
* Performance critical fermion operators called inside the inverter
|
||||||
|
********************************************************************/
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
int pm = this->pm;
|
||||||
|
RealD shift = this->shift;
|
||||||
|
RealD mq1 = this->mq1;
|
||||||
|
RealD mq2 = this->mq2;
|
||||||
|
RealD mq3 = this->mq3;
|
||||||
|
|
||||||
|
// coefficients for shift operator ( = shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} )
|
||||||
|
Coeff_t shiftp(0.0), shiftm(0.0);
|
||||||
|
if(shift != 0.0){
|
||||||
|
if(pm == 1){ shiftp = shift*(mq3-mq2); }
|
||||||
|
else{ shiftm = -shift*(mq3-mq2); }
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<Coeff_t> diag(Ls,1.0);
|
||||||
|
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm;
|
||||||
|
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftp;
|
||||||
|
|
||||||
|
#if(0)
|
||||||
|
std::cout << GridLogMessage << "DomainWallEOFAFermion::M5D(FF&,FF&):" << std::endl;
|
||||||
|
for(int i=0; i<diag.size(); ++i){
|
||||||
|
std::cout << GridLogMessage << "diag[" << i << "] =" << diag[i] << std::endl;
|
||||||
|
}
|
||||||
|
for(int i=0; i<upper.size(); ++i){
|
||||||
|
std::cout << GridLogMessage << "upper[" << i << "] =" << upper[i] << std::endl;
|
||||||
|
}
|
||||||
|
for(int i=0; i<lower.size(); ++i){
|
||||||
|
std::cout << GridLogMessage << "lower[" << i << "] =" << lower[i] << std::endl;
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
|
this->M5D(psi, chi, chi, lower, diag, upper);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField& chi)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
int pm = this->pm;
|
||||||
|
RealD shift = this->shift;
|
||||||
|
RealD mq1 = this->mq1;
|
||||||
|
RealD mq2 = this->mq2;
|
||||||
|
RealD mq3 = this->mq3;
|
||||||
|
|
||||||
|
// coefficients for shift operator ( = shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} )
|
||||||
|
Coeff_t shiftp(0.0), shiftm(0.0);
|
||||||
|
if(shift != 0.0){
|
||||||
|
if(pm == 1){ shiftp = shift*(mq3-mq2); }
|
||||||
|
else{ shiftm = -shift*(mq3-mq2); }
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<Coeff_t> diag(Ls,1.0);
|
||||||
|
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
|
||||||
|
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftm;
|
||||||
|
|
||||||
|
#if(0)
|
||||||
|
std::cout << GridLogMessage << "DomainWallEOFAFermion::M5Ddag(FF&,FF&):" << std::endl;
|
||||||
|
for(int i=0; i<diag.size(); ++i){
|
||||||
|
std::cout << GridLogMessage << "diag[" << i << "] =" << diag[i] << std::endl;
|
||||||
|
}
|
||||||
|
for(int i=0; i<upper.size(); ++i){
|
||||||
|
std::cout << GridLogMessage << "upper[" << i << "] =" << upper[i] << std::endl;
|
||||||
|
}
|
||||||
|
for(int i=0; i<lower.size(); ++i){
|
||||||
|
std::cout << GridLogMessage << "lower[" << i << "] =" << lower[i] << std::endl;
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
|
this->M5Ddag(psi, chi, chi, lower, diag, upper);
|
||||||
|
}
|
||||||
|
|
||||||
|
// half checkerboard operations
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& chi)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
|
||||||
|
std::vector<Coeff_t> diag = this->bee;
|
||||||
|
std::vector<Coeff_t> upper(Ls);
|
||||||
|
std::vector<Coeff_t> lower(Ls);
|
||||||
|
|
||||||
|
for(int s=0; s<Ls; s++){
|
||||||
|
upper[s] = -this->cee[s];
|
||||||
|
lower[s] = -this->cee[s];
|
||||||
|
}
|
||||||
|
upper[Ls-1] = this->dm;
|
||||||
|
lower[0] = this->dp;
|
||||||
|
|
||||||
|
this->M5D(psi, psi, chi, lower, diag, upper);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField& chi)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
|
||||||
|
std::vector<Coeff_t> diag = this->bee;
|
||||||
|
std::vector<Coeff_t> upper(Ls);
|
||||||
|
std::vector<Coeff_t> lower(Ls);
|
||||||
|
|
||||||
|
for(int s=0; s<Ls; s++){
|
||||||
|
upper[s] = -this->cee[s];
|
||||||
|
lower[s] = -this->cee[s];
|
||||||
|
}
|
||||||
|
upper[Ls-1] = this->dp;
|
||||||
|
lower[0] = this->dm;
|
||||||
|
|
||||||
|
this->M5Ddag(psi, psi, chi, lower, diag, upper);
|
||||||
|
}
|
||||||
|
|
||||||
|
/****************************************************************************************/
|
||||||
|
|
||||||
|
//Zolo
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
int pm = this->pm;
|
||||||
|
RealD mq1 = this->mq1;
|
||||||
|
RealD mq2 = this->mq2;
|
||||||
|
RealD mq3 = this->mq3;
|
||||||
|
RealD shift = this->shift;
|
||||||
|
|
||||||
|
////////////////////////////////////////////////////////
|
||||||
|
// Constants for the preconditioned matrix Cayley form
|
||||||
|
////////////////////////////////////////////////////////
|
||||||
|
this->bs.resize(Ls);
|
||||||
|
this->cs.resize(Ls);
|
||||||
|
this->aee.resize(Ls);
|
||||||
|
this->aeo.resize(Ls);
|
||||||
|
this->bee.resize(Ls);
|
||||||
|
this->beo.resize(Ls);
|
||||||
|
this->cee.resize(Ls);
|
||||||
|
this->ceo.resize(Ls);
|
||||||
|
|
||||||
|
for(int i=0; i<Ls; ++i){
|
||||||
|
this->bee[i] = 4.0 - this->M5 + 1.0;
|
||||||
|
this->cee[i] = 1.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int i=0; i<Ls; ++i){
|
||||||
|
this->aee[i] = this->cee[i];
|
||||||
|
this->bs[i] = this->beo[i] = 1.0;
|
||||||
|
this->cs[i] = this->ceo[i] = 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
//////////////////////////////////////////
|
||||||
|
// EOFA shift terms
|
||||||
|
//////////////////////////////////////////
|
||||||
|
if(pm == 1){
|
||||||
|
this->dp = mq1*this->cee[0] + shift*(mq3-mq2);
|
||||||
|
this->dm = mq1*this->cee[Ls-1];
|
||||||
|
} else if(this->pm == -1) {
|
||||||
|
this->dp = mq1*this->cee[0];
|
||||||
|
this->dm = mq1*this->cee[Ls-1] - shift*(mq3-mq2);
|
||||||
|
} else {
|
||||||
|
this->dp = mq1*this->cee[0];
|
||||||
|
this->dm = mq1*this->cee[Ls-1];
|
||||||
|
}
|
||||||
|
|
||||||
|
//////////////////////////////////////////
|
||||||
|
// LDU decomposition of eeoo
|
||||||
|
//////////////////////////////////////////
|
||||||
|
this->dee.resize(Ls+1);
|
||||||
|
this->lee.resize(Ls);
|
||||||
|
this->leem.resize(Ls);
|
||||||
|
this->uee.resize(Ls);
|
||||||
|
this->ueem.resize(Ls);
|
||||||
|
|
||||||
|
for(int i=0; i<Ls; ++i){
|
||||||
|
|
||||||
|
if(i < Ls-1){
|
||||||
|
|
||||||
|
this->lee[i] = -this->cee[i+1]/this->bee[i]; // sub-diag entry on the ith column
|
||||||
|
|
||||||
|
this->leem[i] = this->dm/this->bee[i];
|
||||||
|
for(int j=0; j<i; j++){ this->leem[i] *= this->aee[j]/this->bee[j]; }
|
||||||
|
|
||||||
|
this->dee[i] = this->bee[i];
|
||||||
|
|
||||||
|
this->uee[i] = -this->aee[i]/this->bee[i]; // up-diag entry on the ith row
|
||||||
|
|
||||||
|
this->ueem[i] = this->dp / this->bee[0];
|
||||||
|
for(int j=1; j<=i; j++){ this->ueem[i] *= this->cee[j]/this->bee[j]; }
|
||||||
|
|
||||||
|
} else {
|
||||||
|
|
||||||
|
this->lee[i] = 0.0;
|
||||||
|
this->leem[i] = 0.0;
|
||||||
|
this->uee[i] = 0.0;
|
||||||
|
this->ueem[i] = 0.0;
|
||||||
|
|
||||||
Approx::zolotarev_free(zdata);
|
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
/***************************************************************
|
{
|
||||||
* Additional EOFA operators only called outside the inverter.
|
Coeff_t delta_d = 1.0 / this->bee[0];
|
||||||
* Since speed is not essential, simple axpby-style
|
for(int j=1; j<Ls-1; j++){ delta_d *= this->cee[j] / this->bee[j]; }
|
||||||
* implementations should be fine.
|
this->dee[Ls-1] = this->bee[Ls-1] + this->cee[0] * this->dm * delta_d;
|
||||||
***************************************************************/
|
this->dee[Ls] = this->bee[Ls-1] + this->cee[Ls-1] * this->dp * delta_d;
|
||||||
template<class Impl>
|
}
|
||||||
void DomainWallEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& Din, int sign, int dag)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
|
|
||||||
Din = zero;
|
int inv = 1;
|
||||||
if((sign == 1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, Ls-1, 0); }
|
this->MooeeInternalCompute(0, inv, this->MatpInv, this->MatmInv);
|
||||||
else if((sign == -1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); }
|
this->MooeeInternalCompute(1, inv, this->MatpInvDag, this->MatmInvDag);
|
||||||
else if((sign == 1 ) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, Ls-1); }
|
}
|
||||||
else if((sign == -1) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); }
|
|
||||||
|
// Recompute Cayley-form coefficients for different shift
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::RefreshShiftCoefficients(RealD new_shift)
|
||||||
|
{
|
||||||
|
this->shift = new_shift;
|
||||||
|
Approx::zolotarev_data *zdata = Approx::higham(1.0, this->Ls);
|
||||||
|
this->SetCoefficientsTanh(zdata, 1.0, 0.0);
|
||||||
|
}
|
||||||
|
|
||||||
|
template<class Impl>
|
||||||
|
void DomainWallEOFAFermion<Impl>::MooeeInternalCompute(int dag, int inv,
|
||||||
|
Vector<iSinglet<Simd> >& Matp, Vector<iSinglet<Simd> >& Matm)
|
||||||
|
{
|
||||||
|
int Ls = this->Ls;
|
||||||
|
|
||||||
|
GridBase* grid = this->FermionRedBlackGrid();
|
||||||
|
int LLs = grid->_rdimensions[0];
|
||||||
|
|
||||||
|
if(LLs == Ls){ return; } // Not vectorised in 5th direction
|
||||||
|
|
||||||
|
Eigen::MatrixXcd Pplus = Eigen::MatrixXcd::Zero(Ls,Ls);
|
||||||
|
Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls);
|
||||||
|
|
||||||
|
for(int s=0; s<Ls; s++){
|
||||||
|
Pplus(s,s) = this->bee[s];
|
||||||
|
Pminus(s,s) = this->bee[s];
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int s=0; s<Ls-1; s++){
|
||||||
|
Pminus(s,s+1) = -this->cee[s];
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int s=0; s<Ls-1; s++){
|
||||||
|
Pplus(s+1,s) = -this->cee[s+1];
|
||||||
|
}
|
||||||
|
|
||||||
|
Pplus (0,Ls-1) = this->dp;
|
||||||
|
Pminus(Ls-1,0) = this->dm;
|
||||||
|
|
||||||
|
Eigen::MatrixXcd PplusMat ;
|
||||||
|
Eigen::MatrixXcd PminusMat;
|
||||||
|
|
||||||
|
#if(0)
|
||||||
|
std::cout << GridLogMessage << "Pplus:" << std::endl;
|
||||||
|
for(int s=0; s<Ls; ++s){
|
||||||
|
for(int ss=0; ss<Ls; ++ss){
|
||||||
|
std::cout << Pplus(s,ss) << "\t";
|
||||||
}
|
}
|
||||||
|
std::cout << std::endl;
|
||||||
// This is just the identity for DWF
|
}
|
||||||
template<class Impl>
|
std::cout << GridLogMessage << "Pminus:" << std::endl;
|
||||||
void DomainWallEOFAFermion<Impl>::Dtilde(const FermionField& psi, FermionField& chi){ chi = psi; }
|
for(int s=0; s<Ls; ++s){
|
||||||
|
for(int ss=0; ss<Ls; ++ss){
|
||||||
// This is just the identity for DWF
|
std::cout << Pminus(s,ss) << "\t";
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionField& chi){ chi = psi; }
|
|
||||||
|
|
||||||
/*****************************************************************************************************/
|
|
||||||
|
|
||||||
template<class Impl>
|
|
||||||
RealD DomainWallEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
|
|
||||||
FermionField Din(psi._grid);
|
|
||||||
|
|
||||||
this->Meooe5D(psi, Din);
|
|
||||||
this->DW(Din, chi, DaggerNo);
|
|
||||||
axpby(chi, 1.0, 1.0, chi, psi);
|
|
||||||
this->M5D(psi, chi);
|
|
||||||
return(norm2(chi));
|
|
||||||
}
|
}
|
||||||
|
std::cout << std::endl;
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
template<class Impl>
|
if(inv) {
|
||||||
RealD DomainWallEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
|
PplusMat = Pplus.inverse();
|
||||||
{
|
PminusMat = Pminus.inverse();
|
||||||
int Ls = this->Ls;
|
} else {
|
||||||
|
PplusMat = Pplus;
|
||||||
|
PminusMat = Pminus;
|
||||||
|
}
|
||||||
|
|
||||||
FermionField Din(psi._grid);
|
if(dag){
|
||||||
|
PplusMat.adjointInPlace();
|
||||||
|
PminusMat.adjointInPlace();
|
||||||
|
}
|
||||||
|
|
||||||
this->DW(psi, Din, DaggerYes);
|
typedef typename SiteHalfSpinor::scalar_type scalar_type;
|
||||||
this->MeooeDag5D(Din, chi);
|
const int Nsimd = Simd::Nsimd();
|
||||||
this->M5Ddag(psi, chi);
|
Matp.resize(Ls*LLs);
|
||||||
axpby(chi, 1.0, 1.0, chi, psi);
|
Matm.resize(Ls*LLs);
|
||||||
return(norm2(chi));
|
|
||||||
}
|
|
||||||
|
|
||||||
/********************************************************************
|
for(int s2=0; s2<Ls; s2++){
|
||||||
* Performance critical fermion operators called inside the inverter
|
for(int s1=0; s1<LLs; s1++){
|
||||||
********************************************************************/
|
int istride = LLs;
|
||||||
|
int ostride = 1;
|
||||||
|
Simd Vp;
|
||||||
|
Simd Vm;
|
||||||
|
scalar_type *sp = (scalar_type*) &Vp;
|
||||||
|
scalar_type *sm = (scalar_type*) &Vm;
|
||||||
|
for(int l=0; l<Nsimd; l++){
|
||||||
|
if(switcheroo<Coeff_t>::iscomplex()) {
|
||||||
|
sp[l] = PplusMat (l*istride+s1*ostride,s2);
|
||||||
|
sm[l] = PminusMat(l*istride+s1*ostride,s2);
|
||||||
|
} else {
|
||||||
|
// if real
|
||||||
|
scalar_type tmp;
|
||||||
|
tmp = PplusMat (l*istride+s1*ostride,s2);
|
||||||
|
sp[l] = scalar_type(tmp.real(),tmp.real());
|
||||||
|
tmp = PminusMat(l*istride+s1*ostride,s2);
|
||||||
|
sm[l] = scalar_type(tmp.real(),tmp.real());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
Matp[LLs*s2+s1] = Vp;
|
||||||
|
Matm[LLs*s2+s1] = Vm;
|
||||||
|
}}
|
||||||
|
}
|
||||||
|
|
||||||
template<class Impl>
|
FermOpTemplateInstantiate(DomainWallEOFAFermion);
|
||||||
void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)
|
GparityFermOpTemplateInstantiate(DomainWallEOFAFermion);
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
int pm = this->pm;
|
|
||||||
RealD shift = this->shift;
|
|
||||||
RealD mq1 = this->mq1;
|
|
||||||
RealD mq2 = this->mq2;
|
|
||||||
RealD mq3 = this->mq3;
|
|
||||||
|
|
||||||
// coefficients for shift operator ( = shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} )
|
NAMESPACE_END(Grid);
|
||||||
Coeff_t shiftp(0.0), shiftm(0.0);
|
|
||||||
if(shift != 0.0){
|
|
||||||
if(pm == 1){ shiftp = shift*(mq3-mq2); }
|
|
||||||
else{ shiftm = -shift*(mq3-mq2); }
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<Coeff_t> diag(Ls,1.0);
|
|
||||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftm;
|
|
||||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftp;
|
|
||||||
|
|
||||||
#if(0)
|
|
||||||
std::cout << GridLogMessage << "DomainWallEOFAFermion::M5D(FF&,FF&):" << std::endl;
|
|
||||||
for(int i=0; i<diag.size(); ++i){
|
|
||||||
std::cout << GridLogMessage << "diag[" << i << "] =" << diag[i] << std::endl;
|
|
||||||
}
|
|
||||||
for(int i=0; i<upper.size(); ++i){
|
|
||||||
std::cout << GridLogMessage << "upper[" << i << "] =" << upper[i] << std::endl;
|
|
||||||
}
|
|
||||||
for(int i=0; i<lower.size(); ++i){
|
|
||||||
std::cout << GridLogMessage << "lower[" << i << "] =" << lower[i] << std::endl;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
this->M5D(psi, chi, chi, lower, diag, upper);
|
|
||||||
}
|
|
||||||
|
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField& chi)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
int pm = this->pm;
|
|
||||||
RealD shift = this->shift;
|
|
||||||
RealD mq1 = this->mq1;
|
|
||||||
RealD mq2 = this->mq2;
|
|
||||||
RealD mq3 = this->mq3;
|
|
||||||
|
|
||||||
// coefficients for shift operator ( = shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} )
|
|
||||||
Coeff_t shiftp(0.0), shiftm(0.0);
|
|
||||||
if(shift != 0.0){
|
|
||||||
if(pm == 1){ shiftp = shift*(mq3-mq2); }
|
|
||||||
else{ shiftm = -shift*(mq3-mq2); }
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<Coeff_t> diag(Ls,1.0);
|
|
||||||
std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
|
|
||||||
std::vector<Coeff_t> lower(Ls,-1.0); lower[0] = mq1 + shiftm;
|
|
||||||
|
|
||||||
#if(0)
|
|
||||||
std::cout << GridLogMessage << "DomainWallEOFAFermion::M5Ddag(FF&,FF&):" << std::endl;
|
|
||||||
for(int i=0; i<diag.size(); ++i){
|
|
||||||
std::cout << GridLogMessage << "diag[" << i << "] =" << diag[i] << std::endl;
|
|
||||||
}
|
|
||||||
for(int i=0; i<upper.size(); ++i){
|
|
||||||
std::cout << GridLogMessage << "upper[" << i << "] =" << upper[i] << std::endl;
|
|
||||||
}
|
|
||||||
for(int i=0; i<lower.size(); ++i){
|
|
||||||
std::cout << GridLogMessage << "lower[" << i << "] =" << lower[i] << std::endl;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
this->M5Ddag(psi, chi, chi, lower, diag, upper);
|
|
||||||
}
|
|
||||||
|
|
||||||
// half checkerboard operations
|
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& chi)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
|
|
||||||
std::vector<Coeff_t> diag = this->bee;
|
|
||||||
std::vector<Coeff_t> upper(Ls);
|
|
||||||
std::vector<Coeff_t> lower(Ls);
|
|
||||||
|
|
||||||
for(int s=0; s<Ls; s++){
|
|
||||||
upper[s] = -this->cee[s];
|
|
||||||
lower[s] = -this->cee[s];
|
|
||||||
}
|
|
||||||
upper[Ls-1] = this->dm;
|
|
||||||
lower[0] = this->dp;
|
|
||||||
|
|
||||||
this->M5D(psi, psi, chi, lower, diag, upper);
|
|
||||||
}
|
|
||||||
|
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField& chi)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
|
|
||||||
std::vector<Coeff_t> diag = this->bee;
|
|
||||||
std::vector<Coeff_t> upper(Ls);
|
|
||||||
std::vector<Coeff_t> lower(Ls);
|
|
||||||
|
|
||||||
for(int s=0; s<Ls; s++){
|
|
||||||
upper[s] = -this->cee[s];
|
|
||||||
lower[s] = -this->cee[s];
|
|
||||||
}
|
|
||||||
upper[Ls-1] = this->dp;
|
|
||||||
lower[0] = this->dm;
|
|
||||||
|
|
||||||
this->M5Ddag(psi, psi, chi, lower, diag, upper);
|
|
||||||
}
|
|
||||||
|
|
||||||
/****************************************************************************************/
|
|
||||||
|
|
||||||
//Zolo
|
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
int pm = this->pm;
|
|
||||||
RealD mq1 = this->mq1;
|
|
||||||
RealD mq2 = this->mq2;
|
|
||||||
RealD mq3 = this->mq3;
|
|
||||||
RealD shift = this->shift;
|
|
||||||
|
|
||||||
////////////////////////////////////////////////////////
|
|
||||||
// Constants for the preconditioned matrix Cayley form
|
|
||||||
////////////////////////////////////////////////////////
|
|
||||||
this->bs.resize(Ls);
|
|
||||||
this->cs.resize(Ls);
|
|
||||||
this->aee.resize(Ls);
|
|
||||||
this->aeo.resize(Ls);
|
|
||||||
this->bee.resize(Ls);
|
|
||||||
this->beo.resize(Ls);
|
|
||||||
this->cee.resize(Ls);
|
|
||||||
this->ceo.resize(Ls);
|
|
||||||
|
|
||||||
for(int i=0; i<Ls; ++i){
|
|
||||||
this->bee[i] = 4.0 - this->M5 + 1.0;
|
|
||||||
this->cee[i] = 1.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
for(int i=0; i<Ls; ++i){
|
|
||||||
this->aee[i] = this->cee[i];
|
|
||||||
this->bs[i] = this->beo[i] = 1.0;
|
|
||||||
this->cs[i] = this->ceo[i] = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
//////////////////////////////////////////
|
|
||||||
// EOFA shift terms
|
|
||||||
//////////////////////////////////////////
|
|
||||||
if(pm == 1){
|
|
||||||
this->dp = mq1*this->cee[0] + shift*(mq3-mq2);
|
|
||||||
this->dm = mq1*this->cee[Ls-1];
|
|
||||||
} else if(this->pm == -1) {
|
|
||||||
this->dp = mq1*this->cee[0];
|
|
||||||
this->dm = mq1*this->cee[Ls-1] - shift*(mq3-mq2);
|
|
||||||
} else {
|
|
||||||
this->dp = mq1*this->cee[0];
|
|
||||||
this->dm = mq1*this->cee[Ls-1];
|
|
||||||
}
|
|
||||||
|
|
||||||
//////////////////////////////////////////
|
|
||||||
// LDU decomposition of eeoo
|
|
||||||
//////////////////////////////////////////
|
|
||||||
this->dee.resize(Ls+1);
|
|
||||||
this->lee.resize(Ls);
|
|
||||||
this->leem.resize(Ls);
|
|
||||||
this->uee.resize(Ls);
|
|
||||||
this->ueem.resize(Ls);
|
|
||||||
|
|
||||||
for(int i=0; i<Ls; ++i){
|
|
||||||
|
|
||||||
if(i < Ls-1){
|
|
||||||
|
|
||||||
this->lee[i] = -this->cee[i+1]/this->bee[i]; // sub-diag entry on the ith column
|
|
||||||
|
|
||||||
this->leem[i] = this->dm/this->bee[i];
|
|
||||||
for(int j=0; j<i; j++){ this->leem[i] *= this->aee[j]/this->bee[j]; }
|
|
||||||
|
|
||||||
this->dee[i] = this->bee[i];
|
|
||||||
|
|
||||||
this->uee[i] = -this->aee[i]/this->bee[i]; // up-diag entry on the ith row
|
|
||||||
|
|
||||||
this->ueem[i] = this->dp / this->bee[0];
|
|
||||||
for(int j=1; j<=i; j++){ this->ueem[i] *= this->cee[j]/this->bee[j]; }
|
|
||||||
|
|
||||||
} else {
|
|
||||||
|
|
||||||
this->lee[i] = 0.0;
|
|
||||||
this->leem[i] = 0.0;
|
|
||||||
this->uee[i] = 0.0;
|
|
||||||
this->ueem[i] = 0.0;
|
|
||||||
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
{
|
|
||||||
Coeff_t delta_d = 1.0 / this->bee[0];
|
|
||||||
for(int j=1; j<Ls-1; j++){ delta_d *= this->cee[j] / this->bee[j]; }
|
|
||||||
this->dee[Ls-1] = this->bee[Ls-1] + this->cee[0] * this->dm * delta_d;
|
|
||||||
this->dee[Ls] = this->bee[Ls-1] + this->cee[Ls-1] * this->dp * delta_d;
|
|
||||||
}
|
|
||||||
|
|
||||||
int inv = 1;
|
|
||||||
this->MooeeInternalCompute(0, inv, this->MatpInv, this->MatmInv);
|
|
||||||
this->MooeeInternalCompute(1, inv, this->MatpInvDag, this->MatmInvDag);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Recompute Cayley-form coefficients for different shift
|
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::RefreshShiftCoefficients(RealD new_shift)
|
|
||||||
{
|
|
||||||
this->shift = new_shift;
|
|
||||||
Approx::zolotarev_data *zdata = Approx::higham(1.0, this->Ls);
|
|
||||||
this->SetCoefficientsTanh(zdata, 1.0, 0.0);
|
|
||||||
}
|
|
||||||
|
|
||||||
template<class Impl>
|
|
||||||
void DomainWallEOFAFermion<Impl>::MooeeInternalCompute(int dag, int inv,
|
|
||||||
Vector<iSinglet<Simd> >& Matp, Vector<iSinglet<Simd> >& Matm)
|
|
||||||
{
|
|
||||||
int Ls = this->Ls;
|
|
||||||
|
|
||||||
GridBase* grid = this->FermionRedBlackGrid();
|
|
||||||
int LLs = grid->_rdimensions[0];
|
|
||||||
|
|
||||||
if(LLs == Ls){ return; } // Not vectorised in 5th direction
|
|
||||||
|
|
||||||
Eigen::MatrixXcd Pplus = Eigen::MatrixXcd::Zero(Ls,Ls);
|
|
||||||
Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls);
|
|
||||||
|
|
||||||
for(int s=0; s<Ls; s++){
|
|
||||||
Pplus(s,s) = this->bee[s];
|
|
||||||
Pminus(s,s) = this->bee[s];
|
|
||||||
}
|
|
||||||
|
|
||||||
for(int s=0; s<Ls-1; s++){
|
|
||||||
Pminus(s,s+1) = -this->cee[s];
|
|
||||||
}
|
|
||||||
|
|
||||||
for(int s=0; s<Ls-1; s++){
|
|
||||||
Pplus(s+1,s) = -this->cee[s+1];
|
|
||||||
}
|
|
||||||
|
|
||||||
Pplus (0,Ls-1) = this->dp;
|
|
||||||
Pminus(Ls-1,0) = this->dm;
|
|
||||||
|
|
||||||
Eigen::MatrixXcd PplusMat ;
|
|
||||||
Eigen::MatrixXcd PminusMat;
|
|
||||||
|
|
||||||
#if(0)
|
|
||||||
std::cout << GridLogMessage << "Pplus:" << std::endl;
|
|
||||||
for(int s=0; s<Ls; ++s){
|
|
||||||
for(int ss=0; ss<Ls; ++ss){
|
|
||||||
std::cout << Pplus(s,ss) << "\t";
|
|
||||||
}
|
|
||||||
std::cout << std::endl;
|
|
||||||
}
|
|
||||||
std::cout << GridLogMessage << "Pminus:" << std::endl;
|
|
||||||
for(int s=0; s<Ls; ++s){
|
|
||||||
for(int ss=0; ss<Ls; ++ss){
|
|
||||||
std::cout << Pminus(s,ss) << "\t";
|
|
||||||
}
|
|
||||||
std::cout << std::endl;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
|
|
||||||
if(inv) {
|
|
||||||
PplusMat = Pplus.inverse();
|
|
||||||
PminusMat = Pminus.inverse();
|
|
||||||
} else {
|
|
||||||
PplusMat = Pplus;
|
|
||||||
PminusMat = Pminus;
|
|
||||||
}
|
|
||||||
|
|
||||||
if(dag){
|
|
||||||
PplusMat.adjointInPlace();
|
|
||||||
PminusMat.adjointInPlace();
|
|
||||||
}
|
|
||||||
|
|
||||||
typedef typename SiteHalfSpinor::scalar_type scalar_type;
|
|
||||||
const int Nsimd = Simd::Nsimd();
|
|
||||||
Matp.resize(Ls*LLs);
|
|
||||||
Matm.resize(Ls*LLs);
|
|
||||||
|
|
||||||
for(int s2=0; s2<Ls; s2++){
|
|
||||||
for(int s1=0; s1<LLs; s1++){
|
|
||||||
int istride = LLs;
|
|
||||||
int ostride = 1;
|
|
||||||
Simd Vp;
|
|
||||||
Simd Vm;
|
|
||||||
scalar_type *sp = (scalar_type*) &Vp;
|
|
||||||
scalar_type *sm = (scalar_type*) &Vm;
|
|
||||||
for(int l=0; l<Nsimd; l++){
|
|
||||||
if(switcheroo<Coeff_t>::iscomplex()) {
|
|
||||||
sp[l] = PplusMat (l*istride+s1*ostride,s2);
|
|
||||||
sm[l] = PminusMat(l*istride+s1*ostride,s2);
|
|
||||||
} else {
|
|
||||||
// if real
|
|
||||||
scalar_type tmp;
|
|
||||||
tmp = PplusMat (l*istride+s1*ostride,s2);
|
|
||||||
sp[l] = scalar_type(tmp.real(),tmp.real());
|
|
||||||
tmp = PminusMat(l*istride+s1*ostride,s2);
|
|
||||||
sm[l] = scalar_type(tmp.real(),tmp.real());
|
|
||||||
}
|
|
||||||
}
|
|
||||||
Matp[LLs*s2+s1] = Vp;
|
|
||||||
Matm[LLs*s2+s1] = Vm;
|
|
||||||
}}
|
|
||||||
}
|
|
||||||
|
|
||||||
FermOpTemplateInstantiate(DomainWallEOFAFermion);
|
|
||||||
GparityFermOpTemplateInstantiate(DomainWallEOFAFermion);
|
|
||||||
|
|
||||||
}}
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user