From 24c07694bca49cd99b3bce460ae4167b4399e8ed Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sun, 14 Oct 2018 00:22:52 +0100 Subject: [PATCH] Mixed precision now supported in MADWF --- Grid/qcd/action/fermion/MADWF.h | 24 +++++++++++--- tests/debug/Test_cayley_cg.cc | 57 ++++++++++++++++++++++----------- 2 files changed, 57 insertions(+), 24 deletions(-) diff --git a/Grid/qcd/action/fermion/MADWF.h b/Grid/qcd/action/fermion/MADWF.h index f7ffea2f..064b13a8 100644 --- a/Grid/qcd/action/fermion/MADWF.h +++ b/Grid/qcd/action/fermion/MADWF.h @@ -30,6 +30,17 @@ Author: Peter Boyle namespace Grid { namespace QCD { +template X=0> +inline void convert(const Fieldi &from,Fieldo &to) +{ + precisionChange(to,from); +} +template X=0> +inline void convert(const Fieldi &from,Fieldo &to) +{ + to=from; +} + template class MADWF { @@ -70,8 +81,10 @@ class MADWF std::cout << GridLogMessage<< " MADWF-like algorithm " << std::endl; std::cout << GridLogMessage<< " ************************************************" << std::endl; - FermionFieldo c0(Mato.GaugeGrid()); // 4d - FermionFieldo y0(Mato.GaugeGrid()); // 4d + FermionFieldi c0i(Mati.GaugeGrid()); // 4d + FermionFieldi y0i(Mati.GaugeGrid()); // 4d + FermionFieldo c0 (Mato.GaugeGrid()); // 4d + FermionFieldo y0 (Mato.GaugeGrid()); // 4d FermionFieldo A(Mato.FermionGrid()); // Temporary outer FermionFieldo B(Mato.FermionGrid()); // Temporary outer @@ -110,15 +123,16 @@ class MADWF // Solve the inner system with surface term c0 //////////////////////////////////////////////// ci = zero; - InsertSlice(c0,ci,0, 0); + convert(c0,c0i); // Possible precison change + InsertSlice(c0i,ci,0, 0); // Dwm P y = Dwm x = D(1) P (c0,0,0,0)^T - Mati.P(ci,Ai); Mati.SetMass(1.0); Mati.M(Ai,srci); Mati.SetMass(m); SchurSolveri(Mati,srci,xi,Guesseri); Mati.Pdag(xi,yi); - ExtractSlice(y0, yi, 0 , 0); + ExtractSlice(y0i, yi, 0 , 0); + convert(y0i,y0); // Possible precision change ////////////////////////////////////// // Propagate solution back to outer system diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc index 361b07e9..5150268f 100644 --- a/tests/debug/Test_cayley_cg.cc +++ b/tests/debug/Test_cayley_cg.cc @@ -85,8 +85,9 @@ void TestReconstruct5D(What & Ddwf, GridParallelRNG *RNG4, GridParallelRNG *RNG5); -template +template void TestReconstruct5DFA(What & Ddwf, + WhatF & DdwfF, LatticeGaugeField &Umu, GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, @@ -102,19 +103,31 @@ int main (int argc, char ** argv) std::cout< seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); LatticeGaugeField Umu(UGrid); + LatticeGaugeFieldF UmuF(UGridF); SU3::HotConfiguration(RNG4,Umu); + precisionChange(UmuF,Umu); std::vector U(4,UGrid); RealD mass=0.1; @@ -123,8 +136,9 @@ int main (int argc, char ** argv) std::cout<(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Ddwf,DdwfF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); RealD b=1.5;// Scale factor b+c=2, b-c=1 RealD c=0.5; @@ -134,8 +148,9 @@ int main (int argc, char ** argv) std::cout<(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dmob,DmobF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dsham,DshamF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); - TestReconstruct5DFA(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + TestReconstruct5DFA(Dov,DovF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout< +template void TestReconstruct5DFA(What & Ddwf, - LatticeGaugeField & Umu, - GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, - GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, - RealD mass, RealD M5, - GridParallelRNG *RNG4, - GridParallelRNG *RNG5) + WhatF & DdwfF, + LatticeGaugeField & Umu, + 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; @@ -326,7 +344,7 @@ void TestReconstruct5DFA(What & Ddwf, double Resid = 1.0e-12; double Residi = 1.0e-5; ConjugateGradient CG(Resid,10000); - ConjugateGradient CGi(Residi,10000); + ConjugateGradient CGi(Residi,10000); Ddwf.ImportPhysicalFermionSource(src4,src); Ddwf.Mdag(src,src_NE); @@ -344,7 +362,8 @@ void TestReconstruct5DFA(What & Ddwf, // Fourier accel PV inverse //////////////////////////// typedef LatticeFermion Field; - typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; + typedef LatticeFermionF FieldF; + typedef SchurRedBlackDiagTwoSolve SchurSolverTypei; typedef PauliVillarsSolverFourierAccel PVinverter; PVinverter PVinverse(Umu,CG); @@ -362,9 +381,9 @@ void TestReconstruct5DFA(What & Ddwf, // Now try MADWF ////////////////////////////// SchurSolverTypei SchurSolver(CGi); - ZeroGuesser Guess; - MADWF > - madwf(Ddwf,Ddwf,PVinverse,SchurSolver,Guess,Resid,10); + ZeroGuesser Guess; + MADWF > + madwf(Ddwf,DdwfF,PVinverse,SchurSolver,Guess,Resid,10); madwf(src4,result_madwf); result_madwf = result_madwf - result;