portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-14 09:45:36 +00:00
Code Releases Activity

4d 5d reconstruction code & test

Browse Source
This commit is contained in:
Peter Boyle 2018-10-09 18:37:20 +01:00
parent e069fd5ed8
commit be40aaf751
4 changed files with 91 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Grid/algorithms/Algorithms.h
View File

@ -52,6 +52,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <Grid/algorithms/CoarsenedMatrix.h> #include <Grid/algorithms/CoarsenedMatrix.h>
#include <Grid/algorithms/FFT.h> #include <Grid/algorithms/FFT.h>
// EigCg // EigCg
// Pcg // Pcg
// Hdcg // Hdcg

20
Grid/qcd/action/fermion/CayleyFermion5D.cc
View File

@ -68,6 +68,26 @@ void CayleyFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &so
ExtractSlice(exported4d, tmp, 0, 0); ExtractSlice(exported4d, tmp, 0, 0);
} }
template<class Impl> template<class Impl>
void CayleyFermion5D<Impl>::P(const FermionField &psi, FermionField &chi)
{
int Ls= this->Ls;
chi=zero;
for(int s=0;s<Ls;s++){
axpby_ssp_pminus(chi,1.0,chi,1.0,psi,s,s);
axpby_ssp_pplus (chi,1.0,chi,1.0,psi,s,(s+1)%Ls);
}
}
template<class Impl>
void CayleyFermion5D<Impl>::Pdag(const FermionField &psi, FermionField &chi)
{
int Ls= this->Ls;
chi=zero;
for(int s=0;s<Ls;s++){
axpby_ssp_pminus(chi,1.0,chi,1.0,psi,s,s);
axpby_ssp_pplus (chi,1.0,chi,1.0,psi,s,(s-1+Ls)%Ls);
}
}
template<class Impl>
void CayleyFermion5D<Impl>::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d) void CayleyFermion5D<Impl>::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d)
{ {
int Ls = this->Ls; int Ls = this->Ls;

8
Grid/qcd/action/fermion/CayleyFermion5D.h
View File

@ -93,6 +93,14 @@ namespace Grid {
virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d); virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);
virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d); virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d);
///////////////////////////////////////////////////////////////
// Support for MADWF tricks
///////////////////////////////////////////////////////////////
RealD Mass(void) { return mass; };
void SetMass(RealD _mass) { mass=_mass; } ;
void P(const FermionField &psi, FermionField &chi);
void Pdag(const FermionField &psi, FermionField &chi);
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// Instantiate different versions depending on Impl // Instantiate different versions depending on Impl
///////////////////////////////////////////////////// /////////////////////////////////////////////////////

62
tests/debug/Test_cayley_cg.cc
View File

@ -27,6 +27,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid/Grid.h> #include <Grid/Grid.h>
#include <Grid/algorithms/iterative/Reconstruct5Dprop.h>
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
@ -75,6 +76,14 @@ void TestCGprec(What & Ddwf,
GridParallelRNG *RNG4, GridParallelRNG *RNG4,
GridParallelRNG *RNG5); GridParallelRNG *RNG5);
template<class What>
void TestReconstruct5D(What & Ddwf,
GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid,
GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid,
RealD mass, RealD M5,
GridParallelRNG *RNG4,
GridParallelRNG *RNG5);
int main (int argc, char ** argv) int main (int argc, char ** argv)
{ {
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
@ -152,6 +161,9 @@ void TestCGinversions(What & Ddwf,
TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
std::cout<<GridLogMessage << "Testing red black Schur inverter"<<std::endl; std::cout<<GridLogMessage << "Testing red black Schur inverter"<<std::endl;
TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
std::cout<<GridLogMessage << "Testing 5D PV reconstruction"<<std::endl;
TestReconstruct5D<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
} }
template<class What> template<class What>
@ -189,6 +201,56 @@ void TestCGprec(What & Ddwf,
CG(HermOpEO,src_o,result_o); CG(HermOpEO,src_o,result_o);
} }
template<class What>
void TestReconstruct5D(What & Ddwf,
GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid,
GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid,
RealD mass, RealD M5,
GridParallelRNG *RNG4,
GridParallelRNG *RNG5)
{
LatticeFermion src4 (UGrid); random(*RNG4,src4);
LatticeFermion res4 (UGrid); res4 = zero;
LatticeFermion src (FGrid);
LatticeFermion src_NE(FGrid);
LatticeFermion result(FGrid);
LatticeFermion result_rec(FGrid);
MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-12,10000);
Ddwf.ImportPhysicalFermionSource(src4,src);
Ddwf.Mdag(src,src_NE);
CG(HermOp,src_NE,result);
Ddwf.ExportPhysicalFermionSolution(result, res4);
Ddwf.M(result,src_NE);
src_NE = src_NE - src;
std::cout <<GridLogMessage<< " True residual is " << norm2(src_NE)<<std::endl;
Reconstruct5DfromPhysical<LatticeFermion> reconstructor(CG);
std::cout <<GridLogMessage<< " True result " << norm2(result)<<std::endl;
std::cout <<GridLogMessage<< " 4d result " << norm2(res4)<<std::endl;
std::cout <<GridLogMessage<< " Reconstructing " <<std::endl;
result_rec = result;
reconstructor(Ddwf,res4,src4,result_rec);
std::cout <<GridLogMessage << "Result "<<norm2(result)<<std::endl;
std::cout <<GridLogMessage << "Result_rec "<<norm2(result_rec)<<std::endl;
result_rec = result_rec - result;
std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
// reconstructor.SliceDump(result_rec);
}
template<class What> template<class What>
void TestCGschur(What & Ddwf, void TestCGschur(What & Ddwf,