diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h index 99d97895..ad2f383d 100644 --- a/lib/qcd/action/fermion/Fermion.h +++ b/lib/qcd/action/fermion/Fermion.h @@ -265,9 +265,9 @@ typedef MobiusFermion GparityMobiusFermionRL; typedef MobiusFermion GparityMobiusFermionFH; typedef MobiusFermion GparityMobiusFermionDF; -typedef MobiusEOFAFermion GparityMobiusFermionR; -typedef MobiusEOFAFermion GparityMobiusFermionF; -typedef MobiusEOFAFermion GparityMobiusFermionD; +typedef MobiusEOFAFermion GparityMobiusEOFAFermionR; +typedef MobiusEOFAFermion GparityMobiusEOFAFermionF; +typedef MobiusEOFAFermion GparityMobiusEOFAFermionD; typedef MobiusEOFAFermion GparityMobiusEOFAFermionRL; typedef MobiusEOFAFermion GparityMobiusEOFAFermionFH; diff --git a/lib/qcd/action/fermion/MobiusEOFAFermion.cc b/lib/qcd/action/fermion/MobiusEOFAFermion.cc index c0837d78..085fa988 100644 --- a/lib/qcd/action/fermion/MobiusEOFAFermion.cc +++ b/lib/qcd/action/fermion/MobiusEOFAFermion.cc @@ -142,7 +142,7 @@ namespace QCD { RealD DtInv_p(0.0), DtInv_m(0.0); RealD N = std::pow(c+d,Ls) + m*std::pow(c-d,Ls); - FermionField tmp = zero; + FermionField tmp(this->FermionGrid()); for(int s=0; s sp) ? 0.0 : std::pow(-1.0,sp-s) * std::pow(c-d,sp-s) / std::pow(c+d,sp-s+1); - if(dag){ - RealD tmp(DtInv_p); - DtInv_p = DtInv_m; - DtInv_m = tmp; + if(sp == 0){ + axpby_ssp_pplus (tmp, 0.0, tmp, DtInv_p, psi, s, sp); + axpby_ssp_pminus(tmp, 0.0, tmp, DtInv_m, psi, s, sp); + } else { + axpby_ssp_pplus (tmp, 1.0, tmp, DtInv_p, psi, s, sp); + axpby_ssp_pminus(tmp, 1.0, tmp, DtInv_m, psi, s, sp); } - axpby_ssp_pplus (tmp, 1.0, tmp, DtInv_p, psi, s, sp); - axpby_ssp_pminus(tmp, 1.0, tmp, DtInv_m, psi, s, sp); - }} } @@ -217,11 +216,11 @@ namespace QCD { template void MobiusEOFAFermion::M5Ddag(const FermionField& psi, FermionField& chi) { - int Ls = this->Ls; + int Ls = this->Ls; std::vector diag(Ls,1.0); - std::vector upper(Ls,-1.0); upper[Ls-1] = this->mq1 + shiftp; - std::vector lower(Ls,-1.0); lower[0] = this->mq1 + shiftm; + std::vector upper(Ls,-1.0); upper[Ls-1] = this->mq1; + std::vector lower(Ls,-1.0); lower[0] = this->mq1; // no shift term if(this->shift == 0.0){ this->M5Ddag(psi, chi, chi, lower, diag, upper); } diff --git a/lib/qcd/action/fermion/MobiusEOFAFermion.h b/lib/qcd/action/fermion/MobiusEOFAFermion.h index 0a8d1788..519b49e7 100644 --- a/lib/qcd/action/fermion/MobiusEOFAFermion.h +++ b/lib/qcd/action/fermion/MobiusEOFAFermion.h @@ -111,7 +111,7 @@ namespace QCD { }; }} -#define INSTANTIATE_DPERP_DWF_EOFA(A)\ +#define INSTANTIATE_DPERP_MOBIUS_EOFA(A)\ template void MobiusEOFAFermion::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi, \ std::vector& lower, std::vector& diag, std::vector& upper); \ template void MobiusEOFAFermion::M5D_shift(const FermionField& psi, const FermionField& phi, FermionField& chi, \ @@ -122,7 +122,7 @@ template void MobiusEOFAFermion::M5Ddag_shift(const FermionField& psi, const std::vector& lower, std::vector& diag, std::vector& upper, std::vector& shift_coeffs); \ template void MobiusEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi); \ template void MobiusEOFAFermion::MooeeInv_shift(const FermionField& psi, FermionField& chi); \ -template void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi); +template void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi); \ template void MobiusEOFAFermion::MooeeInvDag_shift(const FermionField& psi, FermionField& chi); #undef MOBIUS_EOFA_DPERP_DENSE diff --git a/lib/qcd/action/fermion/MobiusEOFAFermioncache.cc b/lib/qcd/action/fermion/MobiusEOFAFermioncache.cc index d184ebe2..420f6390 100644 --- a/lib/qcd/action/fermion/MobiusEOFAFermioncache.cc +++ b/lib/qcd/action/fermion/MobiusEOFAFermioncache.cc @@ -410,19 +410,19 @@ namespace QCD { #ifdef MOBIUS_EOFA_DPERP_CACHE - INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplF); - INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplD); - INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplF); - INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplD); - INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplF); - INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplD); - INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplFH); - INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplDF); - INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplFH); - INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplDF); - INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplFH); - INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplDF); #endif diff --git a/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc b/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc new file mode 100644 index 00000000..d66b8cd9 --- /dev/null +++ b/lib/qcd/action/fermion/MobiusEOFAFermiondense.cc @@ -0,0 +1,184 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/MobiusEOFAFermiondense.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include +#include +#include + +namespace Grid { +namespace QCD { + + /* + * Dense matrix versions of routines + */ + template + void MobiusEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInv_shift(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInvDag_shift(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv) + { + int Ls = this->Ls; + int LLs = psi._grid->_rdimensions[0]; + int vol = psi._grid->oSites()/LLs; + + int pm = this->pm; + RealD shift = this->shift; + RealD alpha = this->alpha; + RealD k = this->k; + RealD mq1 = this->mq1; + + chi.checkerboard = psi.checkerboard; + + assert(Ls==LLs); + + Eigen::MatrixXd Pplus = Eigen::MatrixXd::Zero(Ls,Ls); + Eigen::MatrixXd Pminus = Eigen::MatrixXd::Zero(Ls,Ls); + + for(int s=0;sbee[s]; + Pminus(s,s) = this->bee[s]; + } + + for(int s=0; scee[s]; + } + + for(int s=0; scee[s+1]; + } + Pplus (0,Ls-1) = mq1*this->cee[0]; + Pminus(Ls-1,0) = mq1*this->cee[Ls-1]; + + if(shift != 0.0){ + Coeff_t N = 2.0 * ( std::pow(alpha+1.0,Ls) + mq1*std::pow(alpha-1.0,Ls) ); + for(int s=0; s::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplDF); + + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + + #endif + +}} diff --git a/lib/qcd/action/fermion/MobiusEOFAFermionssp.cc b/lib/qcd/action/fermion/MobiusEOFAFermionssp.cc new file mode 100644 index 00000000..c86bb995 --- /dev/null +++ b/lib/qcd/action/fermion/MobiusEOFAFermionssp.cc @@ -0,0 +1,290 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/MobiusEOFAFermionssp.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include +#include + +namespace Grid { +namespace QCD { + + // FIXME -- make a version of these routines with site loop outermost for cache reuse. + // Pminus fowards + // Pplus backwards + template + void MobiusEOFAFermion::M5D(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; s + void MobiusEOFAFermion::M5D_shift(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper, + std::vector& shift_coeffs) + { + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; spm == 1){ axpby_ssp_pplus(chi, one, chi, shift_coeffs[s], psi, s, Ls-1); } + else{ axpby_ssp_pminus(chi, one, chi, shift_coeffs[s], psi, s, 0); } + } + } + + template + void MobiusEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; s + void MobiusEOFAFermion::M5Ddag_shift(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper, + std::vector& shift_coeffs) + { + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; spm == 1){ axpby_ssp_pplus(chi, one, chi, shift_coeffs[s], psi, Ls-1, s); } + else{ axpby_ssp_pminus(chi, one, chi, shift_coeffs[s], psi, 0, s); } + } + } + + template + void MobiusEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) + { + if(this->shift != 0.0){ MooeeInv_shift(psi,chi); return; } + + Coeff_t one(1.0); + Coeff_t czero(0.0); + chi.checkerboard = psi.checkerboard; + int Ls = this->Ls; + + // Apply (L^{\prime})^{-1} + axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] + for(int s=1; slee[s-1], chi, s, s-1);// recursion Psi[s] -lee P_+ chi[s-1] + } + + // L_m^{-1} + for(int s=0; sleem[s], chi, Ls-1, s); + } + + // U_m^{-1} D^{-1} + for(int s=0; sdee[s], chi, -this->ueem[s]/this->dee[Ls-1], chi, s, Ls-1); + } + axpby_ssp(chi, one/this->dee[Ls-1], chi, czero, chi, Ls-1, Ls-1); + + // Apply U^{-1} + for(int s=Ls-2; s>=0; s--){ + axpby_ssp_pminus(chi, one, chi, -this->uee[s], chi, s, s+1); // chi[Ls] + } + } + + template + void MobiusEOFAFermion::MooeeInv_shift(const FermionField& psi, FermionField& chi) + { + Coeff_t one(1.0); + Coeff_t czero(0.0); + chi.checkerboard = psi.checkerboard; + int Ls = this->Ls; + + FermionField tmp(psi._grid); + + // Apply (L^{\prime})^{-1} + axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] + axpby_ssp(tmp, czero, tmp, this->MooeeInv_shift_lc[0], psi, 0, 0); + for(int s=1; slee[s-1], chi, s, s-1);// recursion Psi[s] -lee P_+ chi[s-1] + axpby_ssp(tmp, one, tmp, this->MooeeInv_shift_lc[s], psi, 0, s); + } + + // L_m^{-1} + for(int s=0; sleem[s], chi, Ls-1, s); + } + + // U_m^{-1} D^{-1} + for(int s=0; sdee[s], chi, -this->ueem[s]/this->dee[Ls-1], chi, s, Ls-1); + } + axpby_ssp(chi, one/this->dee[Ls-1], chi, czero, chi, Ls-1, Ls-1); + + // Apply U^{-1} and add shift term + if(this->pm == 1){ axpby_ssp_pplus(chi, one, chi, this->MooeeInv_shift_norm[Ls-1], tmp, Ls-1, 0); } + else{ axpby_ssp_pminus(chi, one, chi, this->MooeeInv_shift_norm[Ls-1], tmp, Ls-1, 0); } + for(int s=Ls-2; s>=0; s--){ + axpby_ssp_pminus(chi, one, chi, -this->uee[s], chi, s, s+1); // chi[Ls] + if(this->pm == 1){ axpby_ssp_pplus(chi, one, chi, this->MooeeInv_shift_norm[s], tmp, s, 0); } + else{ axpby_ssp_pminus(chi, one, chi, this->MooeeInv_shift_norm[s], tmp, s, 0); } + } + } + + template + void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) + { + if(this->shift != 0.0){ MooeeInvDag_shift(psi,chi); return; } + + Coeff_t one(1.0); + Coeff_t czero(0.0); + chi.checkerboard = psi.checkerboard; + int Ls = this->Ls; + + // Apply (U^{\prime})^{-dagger} + axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] + for(int s=1; suee[s-1]), chi, s, s-1); + } + + // U_m^{-\dagger} + for(int s=0; sueem[s]), chi, Ls-1, s); + } + + // L_m^{-\dagger} D^{-dagger} + for(int s=0; sdee[s]), chi, -conjugate(this->leem[s]/this->dee[Ls-1]), chi, s, Ls-1); + } + axpby_ssp(chi, one/conjugate(this->dee[Ls-1]), chi, czero, chi, Ls-1, Ls-1); + + // Apply L^{-dagger} + for(int s=Ls-2; s>=0; s--){ + axpby_ssp_pplus(chi, one, chi, -conjugate(this->lee[s]), chi, s, s+1); // chi[Ls] + } + } + + template + void MobiusEOFAFermion::MooeeInvDag_shift(const FermionField& psi, FermionField& chi) + { + Coeff_t one(1.0); + Coeff_t czero(0.0); + chi.checkerboard = psi.checkerboard; + int Ls = this->Ls; + + FermionField tmp(psi._grid); + + // Apply (U^{\prime})^{-dagger} and accumulate (MooeeInvDag_shift_lc)_{j} \psi_{j} in tmp[0] + axpby_ssp(chi, one, psi, czero, psi, 0, 0); // chi[0]=psi[0] + axpby_ssp(tmp, czero, tmp, this->MooeeInvDag_shift_lc[0], psi, 0, 0); + for(int s=1; suee[s-1]), chi, s, s-1); + axpby_ssp(tmp, one, tmp, this->MooeeInvDag_shift_lc[s], psi, 0, s); + } + + // U_m^{-\dagger} + for(int s=0; sueem[s]), chi, Ls-1, s); + } + + // L_m^{-\dagger} D^{-dagger} + for(int s=0; sdee[s]), chi, -conjugate(this->leem[s]/this->dee[Ls-1]), chi, s, Ls-1); + } + axpby_ssp(chi, one/conjugate(this->dee[Ls-1]), chi, czero, chi, Ls-1, Ls-1); + + // Apply L^{-dagger} and add shift + if(this->pm == 1){ axpby_ssp_pplus(chi, one, chi, this->MooeeInvDag_shift_norm[Ls-1], tmp, Ls-1, 0); } + else{ axpby_ssp_pminus(chi, one, chi, this->MooeeInvDag_shift_norm[Ls-1], tmp, Ls-1, 0); } + for(int s=Ls-2; s>=0; s--){ + axpby_ssp_pplus(chi, one, chi, -conjugate(this->lee[s]), chi, s, s+1); // chi[Ls] + if(this->pm == 1){ axpby_ssp_pplus(chi, one, chi, this->MooeeInvDag_shift_norm[s], tmp, s, 0); } + else{ axpby_ssp_pminus(chi, one, chi, this->MooeeInvDag_shift_norm[s], tmp, s, 0); } + } + } + + #ifdef MOBIUS_EOFA_DPERP_LINALG + + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplD); + + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(WilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(GparityWilsonImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZWilsonImplDF); + + #endif + +}} diff --git a/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc b/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc new file mode 100644 index 00000000..59544e5a --- /dev/null +++ b/lib/qcd/action/fermion/MobiusEOFAFermionvec.cc @@ -0,0 +1,654 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/action/fermion/MobiusEOFAFermionvec.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include +#include + +namespace Grid { +namespace QCD { + + /* + * Dense matrix versions of routines + */ + template + void MobiusEOFAFermion::MooeeInv(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInv_shift(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerNo, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInvDag(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); + } + + template + void MobiusEOFAFermion::MooeeInvDag_shift(const FermionField& psi, FermionField& chi) + { + this->MooeeInternal(psi, chi, DaggerYes, InverseYes); + } + + template + void MobiusEOFAFermion::M5D(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + GridBase* grid = psi._grid; + int Ls = this->Ls; + int LLs = grid->_rdimensions[0]; + const int nsimd = Simd::Nsimd(); + + Vector> u(LLs); + Vector> l(LLs); + Vector> d(LLs); + + assert(Ls/LLs == nsimd); + assert(phi.checkerboard == psi.checkerboard); + + chi.checkerboard = psi.checkerboard; + + // just directly address via type pun + typedef typename Simd::scalar_type scalar_type; + scalar_type* u_p = (scalar_type*) &u[0]; + scalar_type* l_p = (scalar_type*) &l[0]; + scalar_type* d_p = (scalar_type*) &d[0]; + + for(int o=0; oM5Dcalls++; + this->M5Dtime -= usecond(); + + assert(Nc == 3); + + parallel_for(int ss=0; ssoSites(); ss+=LLs){ // adds LLs + + #if 0 + + alignas(64) SiteHalfSpinor hp; + alignas(64) SiteHalfSpinor hm; + alignas(64) SiteSpinor fp; + alignas(64) SiteSpinor fm; + + for(int v=0; v= v){ rotate(hm, hm, nsimd-1); } + + hp = 0.5*hp; + hm = 0.5*hm; + + spRecon5m(fp, hp); + spRecon5p(fm, hm); + + chi[ss+v] = d[v]*phi[ss+v]; + chi[ss+v] = chi[ss+v] + u[v]*fp; + chi[ss+v] = chi[ss+v] + l[v]*fm; + + } + + #else + + for(int v=0; v(hp_00.v); + hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v); + hp_02.v = Optimization::Rotate::tRotate<2>(hp_02.v); + hp_10.v = Optimization::Rotate::tRotate<2>(hp_10.v); + hp_11.v = Optimization::Rotate::tRotate<2>(hp_11.v); + hp_12.v = Optimization::Rotate::tRotate<2>(hp_12.v); + } + + if(vm >= v){ + hm_00.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_00.v); + hm_01.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_01.v); + hm_02.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_02.v); + hm_10.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_10.v); + hm_11.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_11.v); + hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v); + } + + // Can force these to real arithmetic and save 2x. + Simd p_00 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(0)) + switcheroo::mult(l[v]()()(), hm_00); + Simd p_01 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(1)) + switcheroo::mult(l[v]()()(), hm_01); + Simd p_02 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(2)) + switcheroo::mult(l[v]()()(), hm_02); + Simd p_10 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(0)) + switcheroo::mult(l[v]()()(), hm_10); + Simd p_11 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(1)) + switcheroo::mult(l[v]()()(), hm_11); + Simd p_12 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(2)) + switcheroo::mult(l[v]()()(), hm_12); + Simd p_20 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(0)) + switcheroo::mult(u[v]()()(), hp_00); + Simd p_21 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(1)) + switcheroo::mult(u[v]()()(), hp_01); + Simd p_22 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(2)) + switcheroo::mult(u[v]()()(), hp_02); + Simd p_30 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(0)) + switcheroo::mult(u[v]()()(), hp_10); + Simd p_31 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(1)) + switcheroo::mult(u[v]()()(), hp_11); + Simd p_32 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(2)) + switcheroo::mult(u[v]()()(), hp_12); + + vstream(chi[ss+v]()(0)(0), p_00); + vstream(chi[ss+v]()(0)(1), p_01); + vstream(chi[ss+v]()(0)(2), p_02); + vstream(chi[ss+v]()(1)(0), p_10); + vstream(chi[ss+v]()(1)(1), p_11); + vstream(chi[ss+v]()(1)(2), p_12); + vstream(chi[ss+v]()(2)(0), p_20); + vstream(chi[ss+v]()(2)(1), p_21); + vstream(chi[ss+v]()(2)(2), p_22); + vstream(chi[ss+v]()(3)(0), p_30); + vstream(chi[ss+v]()(3)(1), p_31); + vstream(chi[ss+v]()(3)(2), p_32); + } + + #endif + } + + this->M5Dtime += usecond(); + } + + template + void MobiusEOFAFermion::M5D_shift(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper, + std::vector& shift_coeffs) + { + this->M5D(psi, phi, chi, lower, diag, upper); + + // FIXME: possible gain from vectorizing shift operation as well? + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; spm == 1){ axpby_ssp_pplus(chi, one, chi, shift_coeffs[s], psi, s, Ls-1); } + else{ axpby_ssp_pminus(chi, one, chi, shift_coeffs[s], psi, s, 0); } + } + } + + template + void MobiusEOFAFermion::M5Ddag(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper) + { + GridBase* grid = psi._grid; + int Ls = this->Ls; + int LLs = grid->_rdimensions[0]; + int nsimd = Simd::Nsimd(); + + Vector > u(LLs); + Vector > l(LLs); + Vector > d(LLs); + + assert(Ls/LLs == nsimd); + assert(phi.checkerboard == psi.checkerboard); + + chi.checkerboard = psi.checkerboard; + + // just directly address via type pun + typedef typename Simd::scalar_type scalar_type; + scalar_type* u_p = (scalar_type*) &u[0]; + scalar_type* l_p = (scalar_type*) &l[0]; + scalar_type* d_p = (scalar_type*) &d[0]; + + for(int o=0; oM5Dcalls++; + this->M5Dtime -= usecond(); + + parallel_for(int ss=0; ssoSites(); ss+=LLs){ // adds LLs + + #if 0 + + alignas(64) SiteHalfSpinor hp; + alignas(64) SiteHalfSpinor hm; + alignas(64) SiteSpinor fp; + alignas(64) SiteSpinor fm; + + for(int v=0; v= v){ rotate(hm, hm, nsimd-1); } + + hp = hp*0.5; + hm = hm*0.5; + spRecon5p(fp, hp); + spRecon5m(fm, hm); + + chi[ss+v] = d[v]*phi[ss+v]+u[v]*fp; + chi[ss+v] = chi[ss+v] +l[v]*fm; + + } + + #else + + for(int v=0; v(hp_00.v); + hp_01.v = Optimization::Rotate::tRotate<2>(hp_01.v); + hp_02.v = Optimization::Rotate::tRotate<2>(hp_02.v); + hp_10.v = Optimization::Rotate::tRotate<2>(hp_10.v); + hp_11.v = Optimization::Rotate::tRotate<2>(hp_11.v); + hp_12.v = Optimization::Rotate::tRotate<2>(hp_12.v); + } + + if(vm >= v){ + hm_00.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_00.v); + hm_01.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_01.v); + hm_02.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_02.v); + hm_10.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_10.v); + hm_11.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_11.v); + hm_12.v = Optimization::Rotate::tRotate<2*Simd::Nsimd()-2>(hm_12.v); + } + + Simd p_00 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(0)) + switcheroo::mult(u[v]()()(), hp_00); + Simd p_01 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(1)) + switcheroo::mult(u[v]()()(), hp_01); + Simd p_02 = switcheroo::mult(d[v]()()(), phi[ss+v]()(0)(2)) + switcheroo::mult(u[v]()()(), hp_02); + Simd p_10 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(0)) + switcheroo::mult(u[v]()()(), hp_10); + Simd p_11 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(1)) + switcheroo::mult(u[v]()()(), hp_11); + Simd p_12 = switcheroo::mult(d[v]()()(), phi[ss+v]()(1)(2)) + switcheroo::mult(u[v]()()(), hp_12); + Simd p_20 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(0)) + switcheroo::mult(l[v]()()(), hm_00); + Simd p_21 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(1)) + switcheroo::mult(l[v]()()(), hm_01); + Simd p_22 = switcheroo::mult(d[v]()()(), phi[ss+v]()(2)(2)) + switcheroo::mult(l[v]()()(), hm_02); + Simd p_30 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(0)) + switcheroo::mult(l[v]()()(), hm_10); + Simd p_31 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(1)) + switcheroo::mult(l[v]()()(), hm_11); + Simd p_32 = switcheroo::mult(d[v]()()(), phi[ss+v]()(3)(2)) + switcheroo::mult(l[v]()()(), hm_12); + + vstream(chi[ss+v]()(0)(0), p_00); + vstream(chi[ss+v]()(0)(1), p_01); + vstream(chi[ss+v]()(0)(2), p_02); + vstream(chi[ss+v]()(1)(0), p_10); + vstream(chi[ss+v]()(1)(1), p_11); + vstream(chi[ss+v]()(1)(2), p_12); + vstream(chi[ss+v]()(2)(0), p_20); + vstream(chi[ss+v]()(2)(1), p_21); + vstream(chi[ss+v]()(2)(2), p_22); + vstream(chi[ss+v]()(3)(0), p_30); + vstream(chi[ss+v]()(3)(1), p_31); + vstream(chi[ss+v]()(3)(2), p_32); + + } + + #endif + + } + + this->M5Dtime += usecond(); + } + + template + void MobiusEOFAFermion::M5Ddag_shift(const FermionField& psi, const FermionField& phi, + FermionField& chi, std::vector& lower, std::vector& diag, std::vector& upper, + std::vector& shift_coeffs) + { + this->M5Ddag(psi, phi, chi, lower, diag, upper); + + // FIXME: possible gain from vectorizing shift operation as well? + Coeff_t one(1.0); + int Ls = this->Ls; + for(int s=0; spm == 1){ axpby_ssp_pplus(chi, one, chi, shift_coeffs[s], psi, Ls-1, s); } + else{ axpby_ssp_pminus(chi, one, chi, shift_coeffs[s], psi, 0, s); } + } + } + + #ifdef AVX512 + #include + #include + #include + #endif + + template + void MobiusEOFAFermion::MooeeInternalAsm(const FermionField& psi, FermionField& chi, + int LLs, int site, Vector >& Matp, Vector >& Matm) + { + #ifndef AVX512 + { + SiteHalfSpinor BcastP; + SiteHalfSpinor BcastM; + SiteHalfSpinor SiteChiP; + SiteHalfSpinor SiteChiM; + + // Ls*Ls * 2 * 12 * vol flops + for(int s1=0; s1); + + for(int s1=0; s1 + void MobiusEOFAFermion::MooeeInternalZAsm(const FermionField& psi, FermionField& chi, + int LLs, int site, Vector >& Matp, Vector >& Matm) + { + std::cout << "Error: zMobius not implemented for EOFA" << std::endl; + exit(-1); + }; + + template + void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv) + { + int Ls = this->Ls; + int LLs = psi._grid->_rdimensions[0]; + int vol = psi._grid->oSites()/LLs; + + chi.checkerboard = psi.checkerboard; + + Vector> Matp; + Vector> Matm; + Vector>* _Matp; + Vector>* _Matm; + + // MooeeInternalCompute(dag,inv,Matp,Matm); + if(inv && dag){ + _Matp = &this->MatpInvDag; + _Matm = &this->MatmInvDag; + } + + if(inv && (!dag)){ + _Matp = &this->MatpInv; + _Matm = &this->MatmInv; + } + + if(!inv){ + MooeeInternalCompute(dag, inv, Matp, Matm); + _Matp = &Matp; + _Matm = &Matm; + } + + assert(_Matp->size() == Ls*LLs); + + this->MooeeInvCalls++; + this->MooeeInvTime -= usecond(); + + if(switcheroo::iscomplex()){ + parallel_for(auto site=0; siteMooeeInvTime += usecond(); + } + + #ifdef MOBIUS_EOFA_DPERP_VEC + + INSTANTIATE_DPERP_MOBIUS_EOFA(DomainWallVec5dImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(DomainWallVec5dImplF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZDomainWallVec5dImplD); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZDomainWallVec5dImplF); + + INSTANTIATE_DPERP_MOBIUS_EOFA(DomainWallVec5dImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(DomainWallVec5dImplFH); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZDomainWallVec5dImplDF); + INSTANTIATE_DPERP_MOBIUS_EOFA(ZDomainWallVec5dImplFH); + + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + template void MobiusEOFAFermion::MooeeInternal(const FermionField& psi, FermionField& chi, int dag, int inv); + + #endif + +}} diff --git a/tests/core/Test_mobius_eofa_even_odd.cc b/tests/core/Test_mobius_eofa_even_odd.cc new file mode 100644 index 00000000..fceb58f6 --- /dev/null +++ b/tests/core/Test_mobius_eofa_even_odd.cc @@ -0,0 +1,241 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/core/Test_dwf_eofa_even_odd.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +template +struct scal { + d internal; +}; + +Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + const int Ls = 8; + // GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeFermion src (FGrid); random(RNG5, src); + LatticeFermion phi (FGrid); random(RNG5, phi); + LatticeFermion chi (FGrid); random(RNG5, chi); + LatticeFermion result(FGrid); result = zero; + LatticeFermion ref (FGrid); ref = zero; + LatticeFermion tmp (FGrid); tmp = zero; + LatticeFermion err (FGrid); err = zero; + LatticeGaugeField Umu (UGrid); SU3::HotConfiguration(RNG4, Umu); + std::vector U(4,UGrid); + + // Only one non-zero (y) + Umu = zero; + for(int nn=0; nn0){ U[nn] = zero; } + PokeIndex(Umu, U[nn], nn); + } + + RealD b = 2.5; + RealD c = 1.5; + RealD mq1 = 0.1; + RealD mq2 = 0.5; + RealD mq3 = 1.0; + RealD shift = 0.1234; + RealD M5 = 1.8; + int pm = 1; + MobiusEOFAFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mq1, mq2, mq3, shift, pm, M5, b, c); + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing that Meo + Moe + Moo + Mee = Munprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd, src_o, src); + + Ddwf.Meooe(src_e, r_o); std::cout << GridLogMessage << "Applied Meo" << std::endl; + Ddwf.Meooe(src_o, r_e); std::cout << GridLogMessage << "Applied Moe" << std::endl; + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + Ddwf.Mooee(src_e, r_e); std::cout << GridLogMessage << "Applied Mee" << std::endl; + Ddwf.Mooee(src_o, r_o); std::cout << GridLogMessage << "Applied Moo" << std::endl; + setCheckerboard(r_eeoo, r_e); + setCheckerboard(r_eeoo, r_o); + + r_eo = r_eo + r_eeoo; + Ddwf.M(src, ref); + + // std::cout << GridLogMessage << r_eo << std::endl; + // std::cout << GridLogMessage << ref << std::endl; + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl; + + LatticeComplex cerr(FGrid); + cerr = localInnerProduct(err,err); + // std::cout << GridLogMessage << cerr << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl; + std::cout << GridLogMessage << "= < phi | Deo | chi > * = < chi | Deo^dag| phi> " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + LatticeFermion chi_e (FrbGrid); + LatticeFermion chi_o (FrbGrid); + + LatticeFermion dchi_e(FrbGrid); + LatticeFermion dchi_o(FrbGrid); + + LatticeFermion phi_e (FrbGrid); + LatticeFermion phi_o (FrbGrid); + + LatticeFermion dphi_e(FrbGrid); + LatticeFermion dphi_o(FrbGrid); + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd , phi_o, phi); + + Ddwf.Meooe (chi_e, dchi_o); + Ddwf.Meooe (chi_o, dchi_e); + Ddwf.MeooeDag(phi_e, dphi_o); + Ddwf.MeooeDag(phi_o, dphi_e); + + ComplexD pDce = innerProduct(phi_e, dchi_e); + ComplexD pDco = innerProduct(phi_o, dchi_o); + ComplexD cDpe = innerProduct(chi_e, dphi_e); + ComplexD cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce-conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco-conj(cDpe) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv Mee = 1 " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + + Ddwf.Mooee (chi_e, src_e); + Ddwf.MooeeInv(src_e, phi_e); + + Ddwf.Mooee (chi_o, src_o); + Ddwf.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInvDag MeeDag = 1 " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + + Ddwf.MooeeDag (chi_e, src_e); + Ddwf.MooeeInvDag(src_e, phi_e); + + Ddwf.MooeeDag (chi_o, src_o); + Ddwf.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MpcDagMpc is Hermitian " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + random(RNG5, phi); + random(RNG5, chi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd , chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd , phi_o, phi); + RealD t1,t2; + + SchurDiagMooeeOperator HermOpEO(Ddwf); + HermOpEO.MpcDagMpc(chi_e, dchi_e, t1, t2); + HermOpEO.MpcDagMpc(chi_o, dchi_o, t1, t2); + + HermOpEO.MpcDagMpc(phi_e, dphi_e, t1, t2); + HermOpEO.MpcDagMpc(phi_o, dphi_o, t1, t2); + + pDce = innerProduct(phi_e, dchi_e); + pDco = innerProduct(phi_o, dchi_o); + cDpe = innerProduct(chi_e, dphi_e); + cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDco-conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDce-conj(cDpe) << std::endl; + + Grid_finalize(); +} diff --git a/tests/debug/Test_heatbath_mobius_eofa.cc b/tests/debug/Test_heatbath_mobius_eofa.cc new file mode 100644 index 00000000..a952873d --- /dev/null +++ b/tests/debug/Test_heatbath_mobius_eofa.cc @@ -0,0 +1,104 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/debug/Test_heatbath_dwf_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////////////////// +// This program sets up the initial pseudofermion field |Phi> = Meofa^{-1/2}*|eta>, and +// then uses this Phi to compute the action . +// If all is working, one should find that = . +////////////////////////////////////////////////////////////////////////////////////////// + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +// Parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Npoles = 12; +const RealD b = 2.5; +const RealD c = 1.5; +const RealD mf = 0.01; +const RealD mpv = 1.0; +const RealD M5 = 1.8; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is set up to use " << threads << " threads" << std::endl; + + // Initialize spacetime grid + std::cout << GridLogMessage << "Lattice dimensions: " << grid_dim << " Ls: " << Ls << std::endl; + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(grid_dim, + GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Set up RNGs + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + MobiusEOFAFermionR Lop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mpv, 0.0, -1, M5, b, c); + MobiusEOFAFermionR Rop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mpv, mf, mpv, -1.0, 1, M5, b, c); + + // Construct the action and test the heatbath (zero initial guess) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + // Construct the action and test the heatbath (forecasted initial guesses) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + return 0; +} diff --git a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc new file mode 100644 index 00000000..08c6d566 --- /dev/null +++ b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc @@ -0,0 +1,109 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/debug/Test_heatbath_dwf_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +////////////////////////////////////////////////////////////////////////////////////////// +// This program sets up the initial pseudofermion field |Phi> = Meofa^{-1/2}*|eta>, and +// then uses this Phi to compute the action . +// If all is working, one should find that = . +////////////////////////////////////////////////////////////////////////////////////////// + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +typedef GparityWilsonImplR FermionImplPolicy; +typedef GparityMobiusEOFAFermionR FermionAction; +typedef typename FermionAction::FermionField FermionField; + +// Parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Npoles = 12; +const RealD b = 2.5; +const RealD c = 1.5; +const RealD mf = 0.01; +const RealD mpv = 1.0; +const RealD M5 = 1.8; + +int main(int argc, char** argv) +{ + Grid_init(&argc, &argv); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is set up to use " << threads << " threads" << std::endl; + + // Initialize spacetime grid + std::cout << GridLogMessage << "Lattice dimensions: " << grid_dim << " Ls: " << Ls << std::endl; + GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(grid_dim, + GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian* FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Set up RNGs + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + FermionAction::ImplParams params; + FermionAction Lop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mpv, 0.0, -1, M5, b, c, params); + FermionAction Rop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mpv, mf, mpv, -1.0, 1, M5, b, c, params); + + // Construct the action and test the heatbath (zero initial guess) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + // Construct the action and test the heatbath (forecasted initial guesses) + { + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, Npoles); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); + + Meofa.refresh(Umu, RNG5); + printf(" = %1.15e\n", Meofa.S(Umu)); + } + + return 0; +} diff --git a/tests/debug/Test_reweight_dwf_eofa_gparity.cc b/tests/debug/Test_reweight_dwf_eofa_gparity.cc index b77ec33e..bb0fd98e 100644 --- a/tests/debug/Test_reweight_dwf_eofa_gparity.cc +++ b/tests/debug/Test_reweight_dwf_eofa_gparity.cc @@ -2,7 +2,7 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: ./tests/debug/Test_reweight_dwf_eofa.cc +Source file: ./tests/debug/Test_reweight_dwf_eofa_gparity.cc Copyright (C) 2017 diff --git a/tests/debug/Test_reweight_mobius_eofa.cc b/tests/debug/Test_reweight_mobius_eofa.cc new file mode 100644 index 00000000..c4fa78d0 --- /dev/null +++ b/tests/debug/Test_reweight_mobius_eofa.cc @@ -0,0 +1,215 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/debug/Test_reweight_dwf_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +// parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Nhits = 10; +const int max_iter = 5000; +const RealD b = 2.5; +const RealD c = 1.5; +const RealD mf = 0.1; +const RealD mb = 0.11; +const RealD M5 = 1.8; +const RealD stop_tol = 1.0e-12; + +RealD mean(const std::vector& data) +{ + int N = data.size(); + RealD mean(0.0); + for(int i=0; i& data, int sample) +{ + int N = data.size(); + RealD mean(0.0); + for(int i=0; i& jacks, RealD mean) +{ + int N = jacks.size(); + RealD std(0.0); + for(int i=0; i jack_stats(const std::vector& data) +{ + int N = data.size(); + std::vector jack_samples(N); + std::vector jack_stats(2); + + jack_stats[0] = mean(data); + for(int i=0; i seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + // Initialize RHMC fermion operators + MobiusFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5, b, c); + MobiusFermionR Ddwf_b(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, M5, b, c); + SchurDiagMooeeOperator MdagM(Ddwf_f); + SchurDiagMooeeOperator VdagV(Ddwf_b); + + // Degree 12 rational approximations to x^(1/4) and x^(-1/4) + double lo = 0.0001; + double hi = 95.0; + int precision = 64; + int degree = 12; + AlgRemez remez(lo, hi, precision); + std::cout << GridLogMessage << "Generating degree " << degree << " for x^(1/4)" << std::endl; + remez.generateApprox(degree, 1, 4); + MultiShiftFunction PowerQuarter(remez, stop_tol, false); + MultiShiftFunction PowerNegQuarter(remez, stop_tol, true); + + // Stochastically estimate reweighting factor via RHMC + RealD scale = std::sqrt(0.5); + std::vector rw_rhmc(Nhits); + ConjugateGradientMultiShift msCG_V(max_iter, PowerQuarter); + ConjugateGradientMultiShift msCG_M(max_iter, PowerNegQuarter); + std::cout.precision(12); + + for(int hit=0; hit tmp(2, Ddwf_f.FermionRedBlackGrid()); + gaussian(RNG5, Phi); + Phi = Phi*scale; + + pickCheckerboard(Odd, PhiOdd, Phi); + + // evaluate -log(rw) + msCG_V(VdagV, PhiOdd, tmp[0]); + msCG_M(MdagM, tmp[0], tmp[1]); + rw_rhmc[hit] = norm2(tmp[1]) - norm2(PhiOdd); + std::cout << std::endl << "==================================================" << std::endl; + std::cout << " --- RHMC: Hit " << hit << ": rw = " << rw_rhmc[hit]; + std::cout << std::endl << "==================================================" << std::endl << std::endl; + + } + + // Initialize EOFA fermion operators + RealD shift_L = 0.0; + RealD shift_R = -1.0; + int pm = 1; + MobiusEOFAFermionR Deofa_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5, b, c); + MobiusEOFAFermionR Deofa_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5, b, c); + MdagMLinearOperator LdagL(Deofa_L); + MdagMLinearOperator RdagR(Deofa_R); + + // Stochastically estimate reweighting factor via EOFA + RealD k = Deofa_L.k; + std::vector rw_eofa(Nhits); + ConjugateGradient CG(stop_tol, max_iter); + SchurRedBlackDiagMooeeSolve SchurSolver(CG); + + // Compute -log(Z), where: ( RHMC det ratio ) = Z * ( EOFA det ratio ) + RealD Z = std::pow(b+c+1.0,Ls) + mf*std::pow(b+c-1.0,Ls); + Z /= std::pow(b+c+1.0,Ls) + mb*std::pow(b+c-1.0,Ls); + Z = -12.0*grid_dim[0]*grid_dim[1]*grid_dim[2]*grid_dim[3]*std::log(Z); + + for(int hit=0; hit tmp(2, Deofa_L.FermionGrid()); + gaussian(RNG5, Phi); + Phi = Phi*scale; + + // evaluate -log(rw) + // LH term + for(int s=0; s rhmc_result = jack_stats(rw_rhmc); + std::vector eofa_result = jack_stats(rw_eofa); + std::cout << std::endl << "RHMC: rw = " << rhmc_result[0] << " +/- " << rhmc_result[1] << std::endl; + std::cout << std::endl << "EOFA: rw = " << eofa_result[0] << " +/- " << eofa_result[1] << std::endl; + + Grid_finalize(); +} diff --git a/tests/debug/Test_reweight_mobius_eofa_gparity.cc b/tests/debug/Test_reweight_mobius_eofa_gparity.cc new file mode 100644 index 00000000..11a242d2 --- /dev/null +++ b/tests/debug/Test_reweight_mobius_eofa_gparity.cc @@ -0,0 +1,218 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/debug/Test_reweight_dwf_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: paboyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +typedef typename GparityDomainWallFermionR::FermionField FermionField; + +// parameters for test +const std::vector grid_dim = { 8, 8, 8, 8 }; +const int Ls = 8; +const int Nhits = 10; +const int max_iter = 5000; +const RealD b = 2.5; +const RealD c = 1.5; +const RealD mf = 0.1; +const RealD mb = 0.11; +const RealD M5 = 1.8; +const RealD stop_tol = 1.0e-12; + +RealD mean(const std::vector& data) +{ + int N = data.size(); + RealD mean(0.0); + for(int i=0; i& data, int sample) +{ + int N = data.size(); + RealD mean(0.0); + for(int i=0; i& jacks, RealD mean) +{ + int N = jacks.size(); + RealD std(0.0); + for(int i=0; i jack_stats(const std::vector& data) +{ + int N = data.size(); + std::vector jack_samples(N); + std::vector jack_stats(2); + + jack_stats[0] = mean(data); + for(int i=0; i seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + + // Random gauge field + LatticeGaugeField Umu(UGrid); + SU3::HotConfiguration(RNG4, Umu); + + // Initialize RHMC fermion operators + GparityDomainWallFermionR::ImplParams params; + GparityMobiusFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5, b, c, params); + GparityMobiusFermionR Ddwf_b(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, M5, b, c, params); + SchurDiagMooeeOperator MdagM(Ddwf_f); + SchurDiagMooeeOperator VdagV(Ddwf_b); + + // Degree 12 rational approximations to x^(1/4) and x^(-1/4) + double lo = 0.0001; + double hi = 95.0; + int precision = 64; + int degree = 12; + AlgRemez remez(lo, hi, precision); + std::cout << GridLogMessage << "Generating degree " << degree << " for x^(1/4)" << std::endl; + remez.generateApprox(degree, 1, 4); + MultiShiftFunction PowerQuarter(remez, stop_tol, false); + MultiShiftFunction PowerNegQuarter(remez, stop_tol, true); + + // Stochastically estimate reweighting factor via RHMC + RealD scale = std::sqrt(0.5); + std::vector rw_rhmc(Nhits); + ConjugateGradientMultiShift msCG_V(max_iter, PowerQuarter); + ConjugateGradientMultiShift msCG_M(max_iter, PowerNegQuarter); + std::cout.precision(12); + + for(int hit=0; hit tmp(2, Ddwf_f.FermionRedBlackGrid()); + gaussian(RNG5, Phi); + Phi = Phi*scale; + + pickCheckerboard(Odd, PhiOdd, Phi); + + // evaluate -log(rw) + msCG_V(VdagV, PhiOdd, tmp[0]); + msCG_M(MdagM, tmp[0], tmp[1]); + rw_rhmc[hit] = norm2(tmp[1]) - norm2(PhiOdd); + std::cout << std::endl << "==================================================" << std::endl; + std::cout << " --- RHMC: Hit " << hit << ": rw = " << rw_rhmc[hit]; + std::cout << std::endl << "==================================================" << std::endl << std::endl; + + } + + // Initialize EOFA fermion operators + RealD shift_L = 0.0; + RealD shift_R = -1.0; + int pm = 1; + GparityMobiusEOFAFermionR Deofa_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5, b, c, params); + GparityMobiusEOFAFermionR Deofa_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5, b, c, params); + MdagMLinearOperator LdagL(Deofa_L); + MdagMLinearOperator RdagR(Deofa_R); + + // Stochastically estimate reweighting factor via EOFA + RealD k = Deofa_L.k; + std::vector rw_eofa(Nhits); + ConjugateGradient CG(stop_tol, max_iter); + SchurRedBlackDiagMooeeSolve SchurSolver(CG); + + // Compute -log(Z), where: ( RHMC det ratio ) = Z * ( EOFA det ratio ) + RealD Z = std::pow(b+c+1.0,Ls) + mf*std::pow(b+c-1.0,Ls); + Z /= std::pow(b+c+1.0,Ls) + mb*std::pow(b+c-1.0,Ls); + Z = -12.0*grid_dim[0]*grid_dim[1]*grid_dim[2]*grid_dim[3]*std::log(Z); + + for(int hit=0; hit tmp(2, Deofa_L.FermionGrid()); + gaussian(RNG5, Phi); + Phi = Phi*scale; + + // evaluate -log(rw) + // LH term + for(int s=0; s rhmc_result = jack_stats(rw_rhmc); + std::vector eofa_result = jack_stats(rw_eofa); + std::cout << std::endl << "RHMC: rw = " << rhmc_result[0] << " +/- " << rhmc_result[1] << std::endl; + std::cout << std::endl << "EOFA: rw = " << eofa_result[0] << " +/- " << eofa_result[1] << std::endl; + + Grid_finalize(); +} diff --git a/tests/forces/Test_dwf_gpforce_eofa.cc b/tests/forces/Test_dwf_gpforce_eofa.cc index 50789c89..3afeaa43 100644 --- a/tests/forces/Test_dwf_gpforce_eofa.cc +++ b/tests/forces/Test_dwf_gpforce_eofa.cc @@ -71,9 +71,9 @@ int main (int argc, char** argv) int threads = GridThread::GetThreads(); std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; - LatticeFermion phi (FGrid); gaussian(RNG5, phi); - LatticeFermion Mphi (FGrid); - LatticeFermion MphiPrime (FGrid); + FermionField phi (FGrid); gaussian(RNG5, phi); + FermionField Mphi (FGrid); + FermionField MphiPrime (FGrid); LatticeGaugeField U(UGrid); SU3::HotConfiguration(RNG4,U); diff --git a/tests/forces/Test_mobius_force_eofa.cc b/tests/forces/Test_mobius_force_eofa.cc new file mode 100644 index 00000000..2a5a7d04 --- /dev/null +++ b/tests/forces/Test_mobius_force_eofa.cc @@ -0,0 +1,166 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/forces/Test_dwf_force_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +int main (int argc, char** argv) +{ + Grid_init(&argc, &argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + const int Ls = 8; + + GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Want a different conf at every run + // First create an instance of an engine. + std::random_device rnd_device; + // Specify the engine and distribution. + std::mt19937 mersenne_engine(rnd_device()); + std::uniform_int_distribution dist(1, 100); + + auto gen = std::bind(dist, mersenne_engine); + std::vector seeds4(4); + generate(begin(seeds4), end(seeds4), gen); + + //std::vector seeds4({1,2,3,5}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + LatticeFermion phi (FGrid); gaussian(RNG5, phi); + LatticeFermion Mphi (FGrid); + LatticeFermion MphiPrime (FGrid); + + LatticeGaugeField U(UGrid); + SU3::HotConfiguration(RNG4,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD b = 2.5; + RealD c = 1.5; + RealD mf = 0.01; + RealD mb = 1.0; + RealD M5 = 1.8; + MobiusEOFAFermionR Lop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, 0.0, -1, M5, b, c); + MobiusEOFAFermionR Rop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, -1.0, 1, M5, b, c); + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, 12); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); + + Meofa.refresh(U, RNG5); + RealD S = Meofa.S(U); // pdag M p + + // get the deriv of phidag M phi with respect to "U" + LatticeGaugeField UdSdU(UGrid); + Meofa.deriv(U, UdSdU); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.0001; + + LatticeColourMatrix mommu(UGrid); + LatticeColourMatrix forcemu(UGrid); + LatticeGaugeField mom(UGrid); + LatticeGaugeField Uprime(UGrid); + + for(int mu=0; mu(mom, mommu, mu); + + // fourth order exponential approx + parallel_for(auto i=mom.begin(); i(UdSdU, mu); + mommu = Ta(mommu)*2.0; + PokeIndex(UdSdU, mommu, mu); + } + + for(int mu=0; mu(UdSdU, mu); + mommu = PeekIndex(mom, mu); + + // Update PF action density + dS = dS + trace(mommu*forcemu)*dt; + } + + ComplexD dSpred = sum(dS); + + /*std::cout << GridLogMessage << " S " << S << std::endl; + std::cout << GridLogMessage << " Sprime " << Sprime << std::endl; + std::cout << GridLogMessage << "dS " << Sprime-S << std::endl; + std::cout << GridLogMessage << "predict dS " << dSpred << std::endl;*/ + printf("\nS = %1.15e\n", S); + printf("Sprime = %1.15e\n", Sprime); + printf("dS = %1.15e\n", Sprime - S); + printf("real(dS_predict) = %1.15e\n", dSpred.real()); + printf("imag(dS_predict) = %1.15e\n\n", dSpred.imag()); + + assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ; + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); +} diff --git a/tests/forces/Test_mobius_gpforce_eofa.cc b/tests/forces/Test_mobius_gpforce_eofa.cc new file mode 100644 index 00000000..72f1dee2 --- /dev/null +++ b/tests/forces/Test_mobius_gpforce_eofa.cc @@ -0,0 +1,171 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/forces/Test_dwf_force_eofa.cc + +Copyright (C) 2017 + +Author: Peter Boyle +Author: David Murphy + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +typedef GparityWilsonImplR FermionImplPolicy; +typedef GparityMobiusEOFAFermionR FermionAction; +typedef typename FermionAction::FermionField FermionField; + +int main (int argc, char** argv) +{ + Grid_init(&argc, &argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + const int Ls = 8; + + GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid); + GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid); + + // Want a different conf at every run + // First create an instance of an engine. + std::random_device rnd_device; + // Specify the engine and distribution. + std::mt19937 mersenne_engine(rnd_device()); + std::uniform_int_distribution dist(1, 100); + + auto gen = std::bind(dist, mersenne_engine); + std::vector seeds4(4); + generate(begin(seeds4), end(seeds4), gen); + + //std::vector seeds4({1,2,3,5}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + FermionField phi (FGrid); gaussian(RNG5, phi); + FermionField Mphi (FGrid); + FermionField MphiPrime (FGrid); + + LatticeGaugeField U(UGrid); + SU3::HotConfiguration(RNG4,U); + + //////////////////////////////////// + // Unmodified matrix element + //////////////////////////////////// + RealD b = 2.5; + RealD c = 1.5; + RealD mf = 0.01; + RealD mb = 1.0; + RealD M5 = 1.8; + FermionAction::ImplParams params; + FermionAction Lop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, 0.0, -1, M5, b, c, params); + FermionAction Rop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, -1.0, 1, M5, b, c, params); + OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, 12); + ConjugateGradient CG(1.0e-12, 5000); + ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); + + Meofa.refresh(U, RNG5); + RealD S = Meofa.S(U); // pdag M p + + // get the deriv of phidag M phi with respect to "U" + LatticeGaugeField UdSdU(UGrid); + Meofa.deriv(U, UdSdU); + + //////////////////////////////////// + // Modify the gauge field a little + //////////////////////////////////// + RealD dt = 0.0001; + + LatticeColourMatrix mommu(UGrid); + LatticeColourMatrix forcemu(UGrid); + LatticeGaugeField mom(UGrid); + LatticeGaugeField Uprime(UGrid); + + for(int mu=0; mu(mom, mommu, mu); + + // fourth order exponential approx + parallel_for(auto i=mom.begin(); i(UdSdU, mu); + mommu = Ta(mommu)*2.0; + PokeIndex(UdSdU, mommu, mu); + } + + for(int mu=0; mu(UdSdU, mu); + mommu = PeekIndex(mom, mu); + + // Update PF action density + dS = dS + trace(mommu*forcemu)*dt; + } + + ComplexD dSpred = sum(dS); + + /*std::cout << GridLogMessage << " S " << S << std::endl; + std::cout << GridLogMessage << " Sprime " << Sprime << std::endl; + std::cout << GridLogMessage << "dS " << Sprime-S << std::endl; + std::cout << GridLogMessage << "predict dS " << dSpred << std::endl;*/ + printf("\nS = %1.15e\n", S); + printf("Sprime = %1.15e\n", Sprime); + printf("dS = %1.15e\n", Sprime - S); + printf("real(dS_predict) = %1.15e\n", dSpred.real()); + printf("imag(dS_predict) = %1.15e\n\n", dSpred.imag()); + + assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ; + + std::cout << GridLogMessage << "Done" << std::endl; + Grid_finalize(); +}