From b58fd80379f48fd06d466b9065793e7d02a8b48c Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 6 Oct 2023 13:43:46 -0400 Subject: [PATCH] I/O for coarse op and reorganise multigrid headers --- tests/debug/Test_general_coarse.cc | 1 - tests/debug/Test_general_coarse_hdcg.cc | 144 +++++++++++++++------- tests/debug/Test_general_coarse_pvdagm.cc | 1 - 3 files changed, 99 insertions(+), 47 deletions(-) diff --git a/tests/debug/Test_general_coarse.cc b/tests/debug/Test_general_coarse.cc index 240461b7..6f7f0d6d 100644 --- a/tests/debug/Test_general_coarse.cc +++ b/tests/debug/Test_general_coarse.cc @@ -28,7 +28,6 @@ Author: Peter Boyle #include #include #include -#include #include #include diff --git a/tests/debug/Test_general_coarse_hdcg.cc b/tests/debug/Test_general_coarse_hdcg.cc index 2fe0b90a..29b5fdbd 100644 --- a/tests/debug/Test_general_coarse_hdcg.cc +++ b/tests/debug/Test_general_coarse_hdcg.cc @@ -28,12 +28,71 @@ Author: Peter Boyle #include #include #include -#include +//#include #include using namespace std; using namespace Grid; +template +void SaveOperator(Coarsened &Operator,std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + ScidacWriter WR(Operator.Grid()->IsBoss()); + assert(Operator._A.size()==Operator.geom.npoint); + WR.open(file); + for(int p=0;p +void LoadOperator(Coarsened Operator,std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + Grid::ScidacReader RD ; + RD.open(file); + assert(Operator._A.size()==Operator.geom.npoint); + for(int p=0;p +void SaveBasis(aggregation &Agg,std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + ScidacWriter WR(Agg.FineGrid->IsBoss()); + WR.open(file); + for(int b=0;b +void LoadBasis(aggregation &Agg, std::string file) +{ +#ifdef HAVE_LIME + emptyUserRecord record; + ScidacReader RD ; + RD.open(file); + for(int b=0;b class TestSolver : public LinearFunction { public: TestSolver() {}; @@ -154,47 +213,35 @@ int main (int argc, char ** argv) // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp typedef Aggregation Subspace; Subspace Aggregates(Coarse5d,FrbGrid,cb); -#if 1 - Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis, - 95.0,0.1, - // 400,200,200 -- 48 iters - // 600,200,200 -- 38 iters, 162s - // 600,200,100 -- 38 iters, 169s - // 600,200,50 -- 88 iters. 370s - 600, - 200, - 100, - 0.0); -#else - Aggregates.CreateSubspaceChebyshev(RNG5, - HermOpEO, - nbasis, - // 100.0, - // 0.1, // Low pass is pretty high still -- 311 iters - // 250.0, - // 0.01, // subspace too low filter power wrong - // 250.0, - // 0.2, // slower - 95.0, - // 0.05, // nbasis 12 - 311 -- wrong coarse inv - // 0.05, // nbasis 12 - 154 -- right filt - // 0.1, // nbasis 12 - 169 oops - // 0.05, // nbasis 16 -- 127 iters - // 0.03, // nbasis 16 -- 13- - // 0.1, // nbasis 16 -- 142; sloppy solve - 0.1, // nbasis 24 - 300); -#endif + //////////////////////////////////////////////////////////// // Need to check about red-black grid coarsening //////////////////////////////////////////////////////////// LittleDiracOperator LittleDiracOp(geom,FrbGrid,Coarse5d); - LittleDiracOp.CoarsenOperator(FineHermOp,Aggregates); + bool load=false; + if ( load ) { + LoadBasis(Aggregates,"Subspace.scidac"); + LoadOperator(LittleDiracOp,"LittleDiracOp.scidac"); + } else { + Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis, + 95.0,0.1, + // 400,200,200 -- 48 iters + // 600,200,200 -- 38 iters, 162s + // 600,200,100 -- 38 iters, 169s + // 600,200,50 -- 88 iters. 370s + 600, + 200, + 100, + 0.0); + LittleDiracOp.CoarsenOperator(FineHermOp,Aggregates); + SaveBasis(Aggregates,"Subspace.scidac"); + SaveOperator(LittleDiracOp,"LittleDiracOp.scidac"); + } + // Try projecting to one hop only - std::cout << " projecting coarse matrix "< HermMatrix; HermMatrix CoarseOp (LittleDiracOp); @@ -216,6 +263,7 @@ int main (int argc, char ** argv) std::vector evec(Nm,Coarse5d); CoarseVector c_src(Coarse5d); c_src=1.0; CoarseVector c_res(Coarse5d); + CoarseVector c_ref(Coarse5d); PowerMethod cPM; cPM(CoarseOp,c_src); @@ -234,7 +282,7 @@ int main (int argc, char ** argv) // HPDSolver HPDSolve(CoarseOp,CG,CoarseZeroGuesser); HPDSolver HPDSolve(CoarseOp,CG,DeflCoarseGuesser); c_res=Zero(); - HPDSolve(c_src,c_res); + HPDSolve(c_src,c_res); c_ref = c_res; ////////////////////////////////////////////////////////////////////////// // Deflated (with real op EV's) solve for the projected coarse op @@ -243,12 +291,14 @@ int main (int argc, char ** argv) HPDSolver HPDSolveProj(CoarseOpProj,CG,DeflCoarseGuesser); c_res=Zero(); HPDSolveProj(c_src,c_res); + c_res = c_res - c_ref; + std::cout << "Projected solver error "< CoarseSmoother(2.0,40.,8,CoarseOpProj); // 311 - ChebyshevSmoother CoarseSmoother(2.0,36.,12,CoarseOpProj); // 311 + ChebyshevSmoother CoarseSmoother(4.0,45.,16,CoarseOpProj); // 311 + // ChebyshevSmoother CoarseSmoother(0.5,36.,10,CoarseOpProj); // 311 //////////////////////////////////////////////////////// // CG, Cheby mode spacing 200,200 @@ -275,10 +325,14 @@ int main (int argc, char ** argv) c_res=Zero(); cADEF1(c_src,c_res); + c_res = c_res - c_ref; + std::cout << "cADEF1 solver error "< - HDCGdefl(1.0e-8, 3000, - FineHermOp, - Smoother, - cADEF1, - HPDSolve, - Aggregates); + HDCGdefl(1.0e-8, 100, + FineHermOp, + Smoother, + cADEF1, + HPDSolve, + Aggregates); result=Zero(); HDCGdefl(src,result); diff --git a/tests/debug/Test_general_coarse_pvdagm.cc b/tests/debug/Test_general_coarse_pvdagm.cc index 269229a7..b0b9d4d8 100644 --- a/tests/debug/Test_general_coarse_pvdagm.cc +++ b/tests/debug/Test_general_coarse_pvdagm.cc @@ -28,7 +28,6 @@ Author: Peter Boyle #include #include #include -#include #include #include