From 31cb8db506e0a67bcd2a7d79fa16b8e581aa0809 Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Thu, 14 Mar 2019 16:06:09 -0400 Subject: [PATCH 1/7] fixed a bug of XML parsing --- lib/meas/inline/io/inline_qio_read_obj.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/meas/inline/io/inline_qio_read_obj.cc b/lib/meas/inline/io/inline_qio_read_obj.cc index c6546b394..f2142a90d 100644 --- a/lib/meas/inline/io/inline_qio_read_obj.cc +++ b/lib/meas/inline/io/inline_qio_read_obj.cc @@ -364,7 +364,7 @@ namespace Chroma // BJOO: Yes that's all very well. but RecordXML has to hold something XMLBufferWriter dummy; push(dummy,"DummyRecordXML"); - write(dummy, "std::mapSize", num_vecs); + write(dummy, "std_mapSize", num_vecs); pop(dummy); record_xml.open(dummy); From 160f7b9ad443004487b9fc5d86b2783d8c63b840 Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Sun, 29 Sep 2019 21:52:06 -0500 Subject: [PATCH 2/7] add multigrid MdagM solver based on mg_proto --- lib/Makefile.am | 6 +- .../invert/mg_proto/mg_proto_qphix_helpers.h | 5 + ...yssolver_mdagm_clover_mg_proto_qphix_eo.cc | 291 ++++++++++++++++++ ...syssolver_mdagm_clover_mg_proto_qphix_eo.h | 82 +++++ .../qphix/syssolver_mdagm_clover_qphix_w.h | 142 ++++----- .../ferm/invert/syssolver_mdagm_aggregate.cc | 13 + other_libs/cpp_wilson_dslash | 2 +- 7 files changed, 467 insertions(+), 74 deletions(-) create mode 100644 lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc create mode 100644 lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.h diff --git a/lib/Makefile.am b/lib/Makefile.am index 1f8dc30af..9cd7c5edd 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -2121,13 +2121,15 @@ nobase_include_HEADERS += \ actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.h \ meas/inline/io/inline_erase_mg_proto_qphix_space.h \ actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix.h \ - actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.h + actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.h \ + actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.h libchroma_a_SOURCES += \ actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.cc \ meas/inline/io/inline_erase_mg_proto_qphix_space.cc \ actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix.cc \ - actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.cc + actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.cc \ + actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc endif endif diff --git a/lib/actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.h b/lib/actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.h index 07dc7d586..9ca854ad4 100644 --- a/lib/actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.h +++ b/lib/actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.h @@ -40,6 +40,11 @@ using MGPreconditionerEO = MGPreconditionerT +std::shared_ptr +createFineLinOpT( const MGProtoSolverParams& params, const multi1d& u, + const MG::LatticeInfo& info); + std::shared_ptr createFineLinOp( const MGProtoSolverParams& params, const multi1d& u, const MG::LatticeInfo& info); diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc new file mode 100644 index 000000000..dc17c8a07 --- /dev/null +++ b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc @@ -0,0 +1,291 @@ +/* + * solve a MdagM*psi=chi system by mg_proto + */ + +#include "chromabase.h" +#include "handle.h" +#include "state.h" +#include "actions/ferm/invert/syssolver_mdagm_factory.h" +#include "actions/ferm/invert/syssolver_mdagm_aggregate.h" +#include "actions/ferm/invert/syssolver_linop_factory.h" +#include "actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.h" +#include "lattice/qphix/qphix_blas_wrappers.h" + +using namespace QDP; + +namespace Chroma +{ + namespace MdagMSysSolverMGProtoQPhiXEOCloverEnv + { + + //! Anonymous namespace + namespace + { + //! Name to be used + const std::string name("MG_PROTO_QPHIX_EO_CLOVER_INVERTER"); + + //! Local registration flag + bool registered = false; + } + + + + // Double precision + MdagMSystemSolver* createFerm(XMLReader& xml_in, + const std::string& path, + Handle< FermState< LatticeFermion, multi1d, multi1d > > state, + Handle< LinearOperator > A) + { + return new MdagMSysSolverMGProtoQPhiXEOClover(A,state,MGProtoSolverParams(xml_in, path)); + } + + //! Register all the factories + bool registerAll() + { + bool success = true; + if (! registered) + { + + success &= Chroma::TheMdagMFermSystemSolverFactory::Instance().registerObject(name, createFerm); + registered = true; + } + return success; + } + }; + + // Save me typing, by exposing this file level from here + using T = MdagMSysSolverMGProtoQPhiXEOClover::T; + using Q = MdagMSysSolverMGProtoQPhiXEOClover::Q; + + + // Constructor + MdagMSysSolverMGProtoQPhiXEOClover::MdagMSysSolverMGProtoQPhiXEOClover(Handle< LinearOperator > A_, + Handle< FermState > state_, + const MGProtoSolverParams& invParam_) : + A(A_), state(state_), invParam(invParam_), subspaceId(invParam_.SubspaceId){ + using namespace MGProtoHelpersQPhiX; + + // Create the new Dirac operator with latest gauge field, which is necessary for + // every MD step in HMC, while keep the same multigrid near null subspace + // until the solver failed, then recalculate the subspace with current gauge field. + // This may save time for the multigrid setup stage. + auto mg_levels = std::make_shared(); + + IndexArray latdims = {{QDP::Layout::subgridLattSize()[0], + QDP::Layout::subgridLattSize()[1], + QDP::Layout::subgridLattSize()[2], + QDP::Layout::subgridLattSize()[3]}}; + + auto info = std::make_shared(latdims, 4, 3, NodeInfo()); + M_ptr = createFineLinOpT(invParam, state->getLinks(), *info); + + mg_pointer = MGProtoHelpersQPhiX::getMGPreconditionerEO(subspaceId); + if ( ! mg_pointer ) { + QDPIO::cout << "EO MG Preconditioner not found in Named Obj. Creating" << std::endl; + + // Check on the links -- they are ferm state and may already have BC's applied? need to figure that out. + MGProtoHelpersQPhiX::createMGPreconditionerEO(invParam, state->getLinks()); + + // Now get the setup + mg_pointer = MGProtoHelpersQPhiX::getMGPreconditionerEO(subspaceId); + } + + // Next step is to create a solver instance: + MG::FGMRESParams fine_solve_params; + fine_solve_params.MaxIter=invParam.OuterSolverMaxIters; + fine_solve_params.RsdTarget=toDouble(invParam.OuterSolverRsdTarget); + fine_solve_params.VerboseP =invParam.OuterSolverVerboseP; + fine_solve_params.NKrylov = invParam.OuterSolverNKrylov; + + // Internal one with EO preconditioning + using EoFGMRES = const MG::FGMRESSolverQPhiX; + + eo_solver = std::make_shared(*M_ptr, fine_solve_params, (mg_pointer->v_cycle).get()); + wrapped= std::make_shared(eo_solver, M_ptr); + + } + + // Destructor + MdagMSysSolverMGProtoQPhiXEOClover::~MdagMSysSolverMGProtoQPhiXEOClover(){} + + //! Return the subset on which the operator acts + const Subset& + MdagMSysSolverMGProtoQPhiXEOClover::subset(void) const + { + return A->subset(); + } + + SystemSolverResults_t + MdagMSysSolverMGProtoQPhiXEOClover::operator()(T& psi, const T& chi) const + { + START_CODE(); + + QDPIO::cout << "Jolly Greetings from Even-Odd Multigridland" << std::endl; + const Subset& s = A->subset(); + StopWatch swatch; + + swatch.reset(); + swatch.start(); + + QDPIO::cout << "DEBUG: Norm2 Chi Before=" << norm2(chi,s) << std::endl; + const LatticeInfo& info = M_ptr->GetInfo(); + QPhiXSpinor qphix_in(info); + QPhiXSpinor qphix_out(info); + + + // Solve the MdagM psi = chi by two-step + // For asymmetric conditioned operator M_a = A_oo - D_oe A^-1_ee D_eo + // Use M_a^\dagger = \gamma_5 M_a \gamma_5 + // + // First solve M_a Y = \gamma_5 chi', with Y = \gamma_5 M_a psi and chi'=\gamma_5 chi + // Then Solve M_a psi = Y', with Y' = \gamma_5 Y + // psi is what we need + QDPIO::cout<<"MG_PROTO Two Step Solver: "<getLinks()); + mg_pointer = MGProtoHelpersQPhiX::getMGPreconditionerEO(subspaceId); + + // Next step is to create a solver instance: + MG::FGMRESParams fine_solve_params; + fine_solve_params.MaxIter=invParam.OuterSolverMaxIters; + fine_solve_params.RsdTarget=toDouble(invParam.OuterSolverRsdTarget); + fine_solve_params.VerboseP =invParam.OuterSolverVerboseP; + fine_solve_params.NKrylov = invParam.OuterSolverNKrylov; + + // Internal one with EO preconditioning + using EoFGMRES = const MG::FGMRESSolverQPhiX; + M_ptr = mg_pointer->M; + eo_solver = std::make_shared(*M_ptr, fine_solve_params, (mg_pointer->v_cycle).get()); + wrapped= std::make_shared(eo_solver, M_ptr); + + QDPIO::cout<<"MG_PROTO Two Step Solver: "<& predictor) const { + START_CODE(); + + SystemSolverResults_t res; + res = (*this)(psi, chi); + + END_CODE(); + return res; + } + +}; + diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.h b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.h new file mode 100644 index 000000000..d4f5e7295 --- /dev/null +++ b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.h @@ -0,0 +1,82 @@ +/* + * Solve MdagM*psi=chi system by mg_proto + */ + +#ifndef LIB_ACTIONS_FERM_INVERT_MG_PROTO_SYSSOLVER_MDAGM_CLOVER_MG_PROTO_QPHIX_EO_H_ +#define LIB_ACTIONS_FERM_INVERT_MG_PROTO_SYSSOLVER_MDAGM_CLOVER_MG_PROTO_QPHIX_EO_H_ + +#include "chromabase.h" +#include "handle.h" +#include "state.h" +#include "syssolver.h" +#include "linearop.h" +#include "actions/ferm/invert/mg_proto/mgproto_solver_params.h" +#include "actions/ferm/invert/syssolver_linop.h" +#include "actions/ferm/invert/mg_proto/mg_proto_qphix_helpers.h" +#include "lattice/solver.h" +#include "lattice/fgmres_common.h" +#include "lattice/qphix/invfgmres_qphix.h" +#include "lattice/qphix/qphix_qdp_utils.h" +#include "lattice/qphix/qphix_eo_clover_linear_operator.h" +#include "actions/ferm/invert/mg_solver_exception.h" +#include + +using namespace QDP; + +namespace Chroma { + +//! Registration and other yuckies + namespace MdagMSysSolverMGProtoQPhiXEOCloverEnv + { + //! Register the syssolver + bool registerAll(); + + + } + + using EoFGMRES = const MG::FGMRESSolverQPhiX; + + class MdagMSysSolverMGProtoQPhiXEOClover : public MdagMSystemSolver + { + public: + using T = LatticeFermion; + using Q = multi1d; + + MdagMSysSolverMGProtoQPhiXEOClover(Handle< LinearOperator > A_, + Handle< FermState > state_, + const MGProtoSolverParams& invParam_); + + ~MdagMSysSolverMGProtoQPhiXEOClover(); + + + //! Return the subset on which the operator acts + const Subset& subset() const; + + //! Solve It! + SystemSolverResults_t operator()(T& psi, const T& chi) const override; + + SystemSolverResults_t operator()(T& psi, const T& chi, + AbsChronologicalPredictor4D& predictor) const override; + private: + Handle< LinearOperator< T > > A; + Handle< FermState > state; + MGProtoSolverParams invParam; + std::string subspaceId; + mutable std::shared_ptr mg_pointer; + mutable std::shared_ptr M_ptr; + + // Shorthand for the UnprecWrapper + using UnprecFGMRES = MG::UnprecFGMRESSolverQPhiXWrapper; + + + mutable std::shared_ptr wrapped; + mutable std::shared_ptr eo_solver; + + }; + +}; + + + + +#endif /* LIB_ACTIONS_FERM_INVERT_MG_PROTO_SYSSOLVER_MDAGM_CLOVER_MG_PROTO_H_ */ diff --git a/lib/actions/ferm/invert/qphix/syssolver_mdagm_clover_qphix_w.h b/lib/actions/ferm/invert/qphix/syssolver_mdagm_clover_qphix_w.h index 643aa9b04..a91bb200d 100644 --- a/lib/actions/ferm/invert/qphix/syssolver_mdagm_clover_qphix_w.h +++ b/lib/actions/ferm/invert/qphix/syssolver_mdagm_clover_qphix_w.h @@ -58,7 +58,7 @@ namespace Chroma */ using namespace QPhiXVecTraits; - template + template class MdagMSysSolverQPhiXClover : public MdagMSystemSolver { public: @@ -68,7 +68,7 @@ namespace Chroma typedef typename QPhiX::Geometry::Vec,VecTraits::Soa,VecTraits::compress12>::SU3MatrixBlock QPhiX_Gauge; typedef typename QPhiX::Geometry::Vec,VecTraits::Soa,VecTraits::compress12>::CloverBlock QPhiX_Clover; - + //! Constructor /*! * \param M_ Linear operator ( Read ) @@ -76,20 +76,20 @@ namespace Chroma */ MdagMSysSolverQPhiXClover(Handle< LinearOperator > A_, Handle< FermState > state_, - const SysSolverQPhiXCloverParams& invParam_) : + const SysSolverQPhiXCloverParams& invParam_) : A(A_), invParam(invParam_), clov(new CloverTermT()), invclov(new CloverTermT()) { QDPIO::cout << "MdagMSysSolverQPhiXClover:" << std::endl; QDPIO::cout << "AntiPeriodicT is: " << invParam.AntiPeriodicT << std::endl; - + QDPIO::cout << "Veclen is " << VecTraits::Vec << std::endl; QDPIO::cout << "Soalen is " << VecTraits::Soa << std::endl; - if ( VecTraits::Soa > VecTraits::Vec ) { + if ( VecTraits::Soa > VecTraits::Vec ) { QDPIO::cerr << "PROBLEM: Soalen > Veclen. Please set soalen appropriately (<=VECLEN) at compile time" << std::endl; QDP_abort(1); } - + Q u(Nd); for(int mu=0; mu < Nd; mu++) { u[mu] = state_->getLinks()[mu]; @@ -98,12 +98,12 @@ namespace Chroma // Set up aniso coefficients multi1d aniso_coeffs(Nd); for(int mu=0; mu < Nd; mu++) aniso_coeffs[mu] = Real(1); - + bool anisotropy = invParam.CloverParams.anisoParam.anisoP; - if( anisotropy ) { + if( anisotropy ) { aniso_coeffs = makeFermCoeffs( invParam.CloverParams.anisoParam ); } - + double t_boundary=(double)(1); // NB: In this case, the state will have boundaries applied. // So we only need to apply our own boundaries if Compression is enabled @@ -114,7 +114,7 @@ namespace Chroma u[3] *= where(Layout::latticeCoordinate(3) == (Layout::lattSize()[3]-1), Real(t_boundary), Real(1)); } - + cbsize_in_blocks = rb[0].numSiteTable()/VecTraits::Soa; const QPhiX::QPhiXCLIArgs& QPhiXParams = TheQPhiXParams::Instance(); @@ -138,7 +138,7 @@ namespace Chroma pad_xy, pad_xyz, QPhiXParams.getMinCt()); - + QDPIO::cout << " Allocating p and c" << std::endl << std::flush ; #ifndef QPD_IS_QDPJIT @@ -155,61 +155,61 @@ namespace Chroma tmp_qphix = nullptr; #endif - - + + QDPIO::cout << " Allocating Clover" << std::endl << std::flush ; QPhiX_Clover* A_cb0=(QPhiX_Clover *)geom->allocCBClov(); QPhiX_Clover* A_cb1=(QPhiX_Clover *)geom->allocCBClov(); clov_packed[0] = A_cb0; clov_packed[1] = A_cb1; - + QDPIO::cout << " Allocating CloverInv " << std::endl << std::flush ; QPhiX_Clover* A_inv_cb0=(QPhiX_Clover *)geom->allocCBClov(); QPhiX_Clover* A_inv_cb1=(QPhiX_Clover *)geom->allocCBClov(); invclov_packed[0] = A_inv_cb0; invclov_packed[1] = A_inv_cb1; - + // Pack the gauge field QDPIO::cout << "Packing gauge field..." << std::endl << std::flush ; QPhiX_Gauge* packed_gauge_cb0=(QPhiX_Gauge *)geom->allocCBGauge(); QPhiX_Gauge* packed_gauge_cb1=(QPhiX_Gauge *)geom->allocCBGauge(); - + QPhiX::qdp_pack_gauge<>(u, packed_gauge_cb0,packed_gauge_cb1, *geom); u_packed[0] = packed_gauge_cb0; u_packed[1] = packed_gauge_cb1; - - + + QDPIO::cout << "Creating Clover Term" << std::endl; CloverTerm clov_qdp; clov->create(state_, invParam.CloverParams); QDPIO::cout << "Inverting Clover Term" << std::endl; invclov->create(state_, invParam.CloverParams, (*clov)); - for(int cb=0; cb < 2; cb++) { + for(int cb=0; cb < 2; cb++) { invclov->choles(cb); } QDPIO::cout << "Done" << std::endl; QDPIO::cout << "Packing Clover term..." << std::endl; - - for(int cb=0; cb < 2; cb++) { + + for(int cb=0; cb < 2; cb++) { QPhiX::qdp_pack_clover<>((*invclov), invclov_packed[cb], *geom, cb); } - - for(int cb=0; cb < 2; cb++) { + + for(int cb=0; cb < 2; cb++) { QPhiX::qdp_pack_clover<>((*clov), clov_packed[cb], *geom, cb); } QDPIO::cout << "Done" << std::endl; - + QDPIO::cout << "Creating the Even Odd Operator" << std::endl; - M=new QPhiX::EvenOddCloverOperator::Vec,VecTraits::Soa,VecTraits::compress12>(u_packed, - clov_packed[1], - invclov_packed[0], + M=new QPhiX::EvenOddCloverOperator::Vec,VecTraits::Soa,VecTraits::compress12>(u_packed, + clov_packed[1], + invclov_packed[0], geom, t_boundary, toDouble(aniso_coeffs[0]), toDouble(aniso_coeffs[3])); - - - switch( invParam.SolverType ) { + + + switch( invParam.SolverType ) { case CG: { QDPIO::cout << "Creating the CG Solver" << std::endl; @@ -227,12 +227,12 @@ namespace Chroma } } - + //! Destructor is automatic - ~MdagMSysSolverQPhiXClover() + ~MdagMSysSolverQPhiXClover() { - + // Need to unalloc all the memory... QDPIO::cout << "Destructing" << std::endl; @@ -247,17 +247,17 @@ namespace Chroma psi_qphix = nullptr; chi_qphix = nullptr; tmp_qphix = nullptr; - + geom->free(invclov_packed[0]); geom->free(invclov_packed[1]); invclov_packed[0] = nullptr; invclov_packed[1] = nullptr; - + geom->free(clov_packed[0]); geom->free(clov_packed[1]); clov_packed[0] = nullptr; clov_packed[1] = nullptr; - + geom->free(u_packed[0]); geom->free(u_packed[1]); u_packed[0] = nullptr; @@ -270,18 +270,18 @@ namespace Chroma //! Return the subset on which the operator acts const Subset& subset() const {return A->subset();} - - SystemSolverResults_t operator()(T& psi, const T& chi, - Chroma::AbsChronologicalPredictor4D& predictor) const + + SystemSolverResults_t operator()(T& psi, const T& chi, + Chroma::AbsChronologicalPredictor4D& predictor) const { - + START_CODE(); StopWatch swatch; swatch.reset(); swatch.start(); /* Factories here later? */ SystemSolverResults_t res; - switch( invParam.SolverType ) { + switch( invParam.SolverType ) { case CG: { res = cgSolve(psi,chi,predictor); @@ -303,7 +303,7 @@ namespace Chroma return res; } - + //! Solver the linear system /*! * \param psi solution ( Modify ) @@ -311,14 +311,14 @@ namespace Chroma * \return syssolver results */ - // NEED THE CHRONO Predictor argument here... + // NEED THE CHRONO Predictor argument here... SystemSolverResults_t operator()(T& psi, const T& chi) const { START_CODE(); SystemSolverResults_t res; Null4DChronoPredictor not_predicting; res=(*this)(psi,chi, not_predicting); - + END_CODE(); return res; } @@ -330,7 +330,7 @@ namespace Chroma private: // Hide default constructor MdagMSysSolverQPhiXClover() {} - + Handle< LinearOperator > A; @@ -339,17 +339,17 @@ namespace Chroma Handle< CloverTermT > invclov; QPhiX::Geometry::Vec, VecTraits::Soa, VecTraits::compress12>* geom; - + Handle< QPhiX::EvenOddCloverOperator::Vec, VecTraits::Soa, VecTraits::compress12> > M; Handle< QPhiX::InvCG::Vec, VecTraits::Soa, VecTraits::compress12> > cg_solver; Handle< QPhiX::InvBiCGStab::Vec, VecTraits::Soa, VecTraits::compress12> > bicgstab_solver; - + QPhiX_Clover* invclov_packed[2]; QPhiX_Clover* clov_packed[2]; QPhiX_Gauge* u_packed[2]; - + mutable QPhiX_Spinor* psi_qphix; mutable QPhiX_Spinor* chi_qphix; mutable QPhiX_Spinor* tmp_qphix; @@ -374,7 +374,7 @@ namespace Chroma double rsd_final; unsigned long site_flops=0; unsigned long mv_apps=0; - + double start = omp_get_wtime(); int my_isign=1; (*cg_solver)(psi_qphix, chi_qphix, toDouble(invParam.RsdTarget), res.n_count, rsd_final, site_flops, mv_apps, my_isign, invParam.VerboseP); @@ -386,7 +386,7 @@ namespace Chroma predictor.newVector(psi); - // Chi Should now hold the result spinor + // Chi Should now hold the result spinor // Check it against chroma. { T r = chi; @@ -394,22 +394,22 @@ namespace Chroma (*A)(tmp, psi, PLUS); (*A)(tmp2, tmp, MINUS); r[ A->subset() ] -= tmp2; - + Double r2 = norm2(r,A->subset()); Double b2 = norm2(chi, A->subset()); Double rel_resid = sqrt(r2/b2); res.resid = rel_resid; - QDPIO::cout << "QPHIX_CLOVER_CG_SOLVER: " << res.n_count << " iters, rsd_sq_final=" << rel_resid << std::endl; + QDPIO::cout << "QPHIX_CLOVER_CG_SOLVER: " << res.n_count << " iters, rsd_sq_final=" << rel_resid << std::endl; QDPIO::cout << "QPHIX_CLOVER_CG_SOLVER: || r || / || b || = " << rel_resid << std::endl; - -#if 0 + +#if 1 if ( !toBool ( rel_resid < invParam.RsdTarget*invParam.RsdToleranceFactor ) ) { QDPIO::cout << "SOLVE FAILED" << std::endl; QDP_abort(1); } #endif - } + } int num_cb_sites = Layout::vol()/2; unsigned long total_flops = (site_flops + (1320+504+1320+504+48)*mv_apps)*num_cb_sites; @@ -432,23 +432,23 @@ namespace Chroma SystemSolverResults_t res; Handle< LinearOperator > MdagM( new MdagMLinOp(A) ); - + T Y; Y[ A->subset() ] = psi; // Y is initial guess - + try { - // Try to cast the predictor to a two step predictor - AbsTwoStepChronologicalPredictor4D& two_step_predictor = + // Try to cast the predictor to a two step predictor + AbsTwoStepChronologicalPredictor4D& two_step_predictor = dynamic_cast& >(predictor); - - + + // Predict Y and X separately two_step_predictor.predictY(Y,*A,chi); two_step_predictor.predictX(psi,*MdagM, chi); } catch( std::bad_cast) { - // Not a 2 step predictor. Predict X + // Not a 2 step predictor. Predict X // Then MX = Y is a good guess. predictor(psi,*MdagM, chi); (*A)(Y,psi,PLUS); @@ -516,28 +516,28 @@ namespace Chroma #endif try { - // Try to cast the predictor to a two step predictor - AbsTwoStepChronologicalPredictor4D& two_step_predictor = + // Try to cast the predictor to a two step predictor + AbsTwoStepChronologicalPredictor4D& two_step_predictor = dynamic_cast& >(predictor); - two_step_predictor.newYVector(Y); + two_step_predictor.newYVector(Y); two_step_predictor.newXVector(psi); } catch( std::bad_cast) { - // Not a 2 step predictor. Predict X + // Not a 2 step predictor. Predict X // Then MX = Y is a good guess. predictor.newVector(psi); } - // Chi Should now hold the result spinor + // Chi Should now hold the result spinor // Check it against chroma. -- reuse Y as the residuum Y[ A->subset() ] = chi; { T tmp,tmp2; (*A)(tmp, psi, PLUS); (*A)(tmp2, tmp, MINUS); - + Y[ A->subset() ] -= tmp2; } @@ -549,13 +549,13 @@ namespace Chroma QDPIO::cout << "QPHIX_CLOVER_BICGSTAB_SOLVER: total_iters="< Date: Mon, 30 Sep 2019 16:26:45 -0500 Subject: [PATCH 3/7] add RsdToleranceFactor parameter and a test input --- .../invert/mg_proto/mgproto_solver_params.cc | 5 + .../invert/mg_proto/mgproto_solver_params.h | 1 + ...yssolver_mdagm_clover_mg_proto_qphix_eo.cc | 4 +- tests/hmc/hmc.eoprec_clover.mgproto.ini.xml | 229 ++++++++++++++++++ 4 files changed, 237 insertions(+), 2 deletions(-) create mode 100644 tests/hmc/hmc.eoprec_clover.mgproto.ini.xml diff --git a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc index 2522c5c15..ffd288117 100644 --- a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc +++ b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc @@ -67,6 +67,11 @@ MGProtoSolverParams::MGProtoSolverParams(XMLReader& xml, const std::string& path read( paramtop, "OuterSolverRsdTarget", OuterSolverRsdTarget); read( paramtop, "OuterSolverMaxIters", OuterSolverMaxIters); read( paramtop, "OuterSolverVerboseP", OuterSolverVerboseP); + if( paramtop.count("RsdToleranceFactor") > 0 ) { + read(paramtop, "RsdToleranceFactor", RsdToleranceFactor); + } else { + RsdToleranceFactor = 10; + } read( paramtop, "VCyclePreSmootherMaxIters", VCyclePreSmootherMaxIters, MGLevels-1); read( paramtop, "VCyclePreSmootherRsdTarget", VCyclePreSmootherRsdTarget, MGLevels-1); diff --git a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h index dd80b2900..541e64d1f 100644 --- a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h +++ b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h @@ -39,6 +39,7 @@ struct MGProtoSolverParams { // Details of Outer Solver int OuterSolverNKrylov; Double OuterSolverRsdTarget; + Double RsdToleranceFactor; int OuterSolverMaxIters; bool OuterSolverVerboseP; diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc index dc17c8a07..c4cc422bd 100644 --- a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc +++ b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc @@ -187,7 +187,7 @@ namespace Chroma Double n2rel = n2 / norm2(chi, s); QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: iters = " << res1.n_count + res2.n_count << " rel resid = " << sqrt(n2rel) << std::endl; - if( toBool( sqrt(n2rel) > invParam.OuterSolverRsdTarget ) ) { + if( toBool( sqrt(n2rel) > invParam.OuterSolverRsdTarget * invParam.RsdToleranceFactor ) ) { QDPIO::cout<<"Error in MG_PROTO convergence, retrying..."< + + + + DISORDERED + DUMMY + false + + + + 2846 + 1259 + 980 + 108 + + + 0 + 0 + 5 + 5 + 5 + dummy_run_ + SINGLEFILE + false + false + 20 + false + 1 + true + + + PLAQUETTE + 1 + + 2 + + SIMPLE_GAUGE_STATE + + PERIODIC_GAUGEBC + + + + + default_gauge_field + + + + + + 8 8 8 8 + + + N_FLAVOR_LOGDET_EVEN_EVEN_FERM_MONOMIAL + + CLOVER + -0.0640 + 1.58932722549812 + 0.902783591772098 + + true + 3 + 4.3 + 1.265 + + + STOUT_FERM_STATE + 0.14 + 2 + 3 + + SIMPLE_FERMBC + 1 1 1 -1 + + + + 2 + + HMC::logdet_2flav + + + + TWO_FLAVOR_EOPREC_CONSTDET_FERM_MONOMIAL + + CLOVER + -0.0640 + 1.58932722549812 + 0.902783591772098 + + true + 3 + 4.3 + 1.265 + + + STOUT_FERM_STATE + 0.14 + 2 + 3 + + SIMPLE_FERMBC + 1 1 1 -1 + + + + + MG_PROTO_QPHIX_EO_CLOVER_INVERTER + + -0.0640 + 1.58932722549812 + 0.902783591772098 + + true + 3 + 4.3 + 1.265 + + + true + + 3 + + 2 2 2 2 + 2 2 2 2 + + + 24 16 + 250 250 + 5e-6 5e-6 + 0 0 + + 25 + 1.0e-7 + 100 + true + 5.0 + + 0 0 + 0.1 0.1 + 1.1 1.1 + 0 0 + + 8 8 + 0.1 0.1 + 1.1 1.1 + 0 0 + + 25 25 + 0.1 0.1 + 10 10 + 0 0 + + 5 5 + 0.1 0.1 + 0 0 + + foo + + + HMC::conv_2flav + + + + GAUGE_MONOMIAL + + ANISO_SYM_SPATIAL_GAUGEACT + 1.5 + 0.73356648935927 + 0.0 + + true + 3 + 4.3 + + + PERIODIC_GAUGEBC + + + + HMC::gauge_s + + + + GAUGE_MONOMIAL + + ANISO_SYM_TEMPORAL_GAUGEACT + 1.5 + 0.73356648935927 + 1.0 + 0.0 + + true + 3 + 4.3 + + + PERIODIC_GAUGEBC + + + + HMC::gauge_t + + + + + + HMC::logdet_2flav + HMC::conv_2flav + HMC::gauge_s + HMC::gauge_t + + + + 0.5 + true + 3 + 3.5 + + LCM_STS_LEAPFROG + 10 + + HMC::logdet_2flav + HMC::conv_2flav + HMC::gauge_s + HMC::gauge_t + + + + + From 30e59f26bb57e9e2630f205f9397c9b3782ab667 Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Thu, 3 Oct 2019 15:16:56 -0500 Subject: [PATCH 4/7] recreate multigrid subspace if number of iterations exceed some threshold count --- .../invert/mg_proto/mgproto_solver_params.cc | 6 ++++ .../invert/mg_proto/mgproto_solver_params.h | 1 + ...yssolver_mdagm_clover_mg_proto_qphix_eo.cc | 28 +++++++++++++++---- 3 files changed, 29 insertions(+), 6 deletions(-) diff --git a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc index ffd288117..9e9b431d5 100644 --- a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc +++ b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.cc @@ -73,6 +73,12 @@ MGProtoSolverParams::MGProtoSolverParams(XMLReader& xml, const std::string& path RsdToleranceFactor = 10; } + if( paramtop.count("ThresholdCount") > 0 ) { + read(paramtop, "ThresholdCount", ThresholdCount); + } else { + ThresholdCount = 50; // current setup for test + } + read( paramtop, "VCyclePreSmootherMaxIters", VCyclePreSmootherMaxIters, MGLevels-1); read( paramtop, "VCyclePreSmootherRsdTarget", VCyclePreSmootherRsdTarget, MGLevels-1); read( paramtop, "VCyclePreSmootherRelaxOmega", VCyclePreSmootherRelaxOmega, MGLevels-1); diff --git a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h index 541e64d1f..40c80706b 100644 --- a/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h +++ b/lib/actions/ferm/invert/mg_proto/mgproto_solver_params.h @@ -42,6 +42,7 @@ struct MGProtoSolverParams { Double RsdToleranceFactor; int OuterSolverMaxIters; bool OuterSolverVerboseP; + int ThresholdCount; // VCycle Details (MGLevels - 1 of these) // Presmoother diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc index c4cc422bd..f5f4317ec 100644 --- a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc +++ b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc @@ -173,6 +173,7 @@ namespace Chroma QDPIO::cout << "DEBUG: Norm2 qphix_out_cb_0 = " << qphix_out_norm_cb0 << " Norm psi[0]="<< psi_norm_cb0 << std::endl; QDPIO::cout << "DEBUG: Norm2 qphix_out_cb_1 = " << qphix_out_norm_cb1 << " Norm psi[1]="<< psi_norm_cb1 << std::endl; + res.n_count = res1.n_count + res2.n_count; bool solution_good = true; { // Chroma level check (may be slow) @@ -186,14 +187,30 @@ namespace Chroma Double n2 = norm2(tmp1, s); Double n2rel = n2 / norm2(chi, s); QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: iters = " - << res1.n_count + res2.n_count << " rel resid = " << sqrt(n2rel) << std::endl; + << res.n_count << " rel resid = " << sqrt(n2rel) << std::endl; if( toBool( sqrt(n2rel) > invParam.OuterSolverRsdTarget * invParam.RsdToleranceFactor ) ) { QDPIO::cout<<"Error in MG_PROTO convergence, retrying..."<= invParam.ThresholdCount){ + QDPIO::cout<<"Iteration Threshold Exceeded! iters = "<getLinks()); + refresh.stop(); + + QDPIO::cout<<"Subspace Refreshing Time = "<getLinks()); @@ -253,7 +270,7 @@ namespace Chroma tmp1[s] -= chi; Double n2 = norm2(tmp1, s); Double n2rel = n2 / norm2(chi, s); - QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: iters = "<< res1.n_count + res2.n_count<< " rel resid = " << sqrt(n2rel) << std::endl; + QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: retry iters = "<< res1.n_count + res2.n_count<< " rel resid = " << sqrt(n2rel) << std::endl; if( toBool( sqrt(n2rel) > invParam.OuterSolverRsdTarget * invParam.RsdToleranceFactor ) ) { QDPIO::cout<<"Error in MG_PROTO convergence, exiting"< Date: Fri, 4 Oct 2019 16:38:58 -0500 Subject: [PATCH 5/7] check convergence after each step of two-step solver --- ...yssolver_mdagm_clover_mg_proto_qphix_eo.cc | 171 ++++++++++++------ ...syssolver_mdagm_clover_mg_proto_qphix_eo.h | 6 +- tests/hmc/hmc.eoprec_clover.mgproto.ini.xml | 1 + 3 files changed, 123 insertions(+), 55 deletions(-) diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc index f5f4317ec..070db6290 100644 --- a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc +++ b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc @@ -91,7 +91,6 @@ namespace Chroma } // Next step is to create a solver instance: - MG::FGMRESParams fine_solve_params; fine_solve_params.MaxIter=invParam.OuterSolverMaxIters; fine_solve_params.RsdTarget=toDouble(invParam.OuterSolverRsdTarget); fine_solve_params.VerboseP =invParam.OuterSolverVerboseP; @@ -101,8 +100,6 @@ namespace Chroma using EoFGMRES = const MG::FGMRESSolverQPhiX; eo_solver = std::make_shared(*M_ptr, fine_solve_params, (mg_pointer->v_cycle).get()); - wrapped= std::make_shared(eo_solver, M_ptr); - } // Destructor @@ -137,13 +134,14 @@ namespace Chroma // For asymmetric conditioned operator M_a = A_oo - D_oe A^-1_ee D_eo // Use M_a^\dagger = \gamma_5 M_a \gamma_5 // - // First solve M_a Y = \gamma_5 chi', with Y = \gamma_5 M_a psi and chi'=\gamma_5 chi + // First solve M_a Y = chi', with Y = \gamma_5 M_a psi and chi'=\gamma_5 chi // Then Solve M_a psi = Y', with Y' = \gamma_5 Y // psi is what we need QDPIO::cout<<"MG_PROTO Two Step Solver: "<= invParam.ThresholdCount){ - QDPIO::cout<<"Iteration Threshold Exceeded! iters = "<getLinks()); + mg_pointer = MGProtoHelpersQPhiX::getMGPreconditionerEO(subspaceId); + using EoFGMRES = const MG::FGMRESSolverQPhiX; + M_ptr = mg_pointer->M; + eo_solver = std::make_shared(*M_ptr, fine_solve_params, (mg_pointer->v_cycle).get()); + refresh.stop(); QDPIO::cout<<"Subspace Refreshing Time = "<getLinks()); mg_pointer = MGProtoHelpersQPhiX::getMGPreconditionerEO(subspaceId); - // Next step is to create a solver instance: - MG::FGMRESParams fine_solve_params; - fine_solve_params.MaxIter=invParam.OuterSolverMaxIters; - fine_solve_params.RsdTarget=toDouble(invParam.OuterSolverRsdTarget); - fine_solve_params.VerboseP =invParam.OuterSolverVerboseP; - fine_solve_params.NKrylov = invParam.OuterSolverNKrylov; - - // Internal one with EO preconditioning using EoFGMRES = const MG::FGMRESSolverQPhiX; M_ptr = mg_pointer->M; eo_solver = std::make_shared(*M_ptr, fine_solve_params, (mg_pointer->v_cycle).get()); - wrapped= std::make_shared(eo_solver, M_ptr); - - QDPIO::cout<<"MG_PROTO Two Step Solver: "<= invParam.ThresholdCount){ + QDPIO::cout<<"Solver-X Iteration Threshold Exceeded! iters = "<getLinks()); + mg_pointer = MGProtoHelpersQPhiX::getMGPreconditionerEO(subspaceId); + + using EoFGMRES = const MG::FGMRESSolverQPhiX; + M_ptr = mg_pointer->M; + eo_solver = std::make_shared(*M_ptr, fine_solve_params, (mg_pointer->v_cycle).get()); ZeroVec(qphix_out, SUBSET_ALL); res2 = (*eo_solver)(qphix_out, qphix_in, RELATIVE); @@ -255,7 +323,7 @@ namespace Chroma Double psi_norm_cb1 = norm2(psi,rb[1]); Double chi_norm_after = norm2(chi,s); - QDPIO::cout << "DEBUG: After solve Norm2 chi = " << chi_norm_after << std::endl; + QDPIO::cout << "DEBUG: After Resolve-X Norm2 chi = " << chi_norm_after << std::endl; QDPIO::cout << "DEBUG: Norm2 qphix_out_cb_0 = " << qphix_out_norm_cb0 << " Norm psi[0]="<< psi_norm_cb0 << std::endl; QDPIO::cout << "DEBUG: Norm2 qphix_out_cb_1 = " << qphix_out_norm_cb1 << " Norm psi[1]="<< psi_norm_cb1 << std::endl; @@ -270,23 +338,24 @@ namespace Chroma tmp1[s] -= chi; Double n2 = norm2(tmp1, s); Double n2rel = n2 / norm2(chi, s); - QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: retry iters = "<< res1.n_count + res2.n_count<< " rel resid = " << sqrt(n2rel) << std::endl; + QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: Resolve-X iters = "< > A; Handle< FermState > state; MGProtoSolverParams invParam; - std::string subspaceId; + MG::FGMRESParams fine_solve_params; + std::string subspaceId; mutable std::shared_ptr mg_pointer; mutable std::shared_ptr M_ptr; // Shorthand for the UnprecWrapper using UnprecFGMRES = MG::UnprecFGMRESSolverQPhiXWrapper; - - - mutable std::shared_ptr wrapped; mutable std::shared_ptr eo_solver; }; diff --git a/tests/hmc/hmc.eoprec_clover.mgproto.ini.xml b/tests/hmc/hmc.eoprec_clover.mgproto.ini.xml index 1cecd975f..a50fa7490 100644 --- a/tests/hmc/hmc.eoprec_clover.mgproto.ini.xml +++ b/tests/hmc/hmc.eoprec_clover.mgproto.ini.xml @@ -133,6 +133,7 @@ 100 true 5.0 + 10 0 0 0.1 0.1 From f9600e6c6f84df7bdd5da75b7536d65b08badcc4 Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Thu, 14 Nov 2019 16:12:22 -0600 Subject: [PATCH 6/7] fix a bug to avoid recreate the subspace at X solve where the same subspace already generated at first Y solve --- .../invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc index 070db6290..9c56511cc 100644 --- a/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc +++ b/lib/actions/ferm/invert/mg_proto/syssolver_mdagm_clover_mg_proto_qphix_eo.cc @@ -263,6 +263,7 @@ namespace Chroma psi = zero; QPhiXSpinorToQDPSpinor(qphix_out, psi); + solution_good = true; res.n_count += res2.n_count; res.resid = res2.resid; { From cfe6572e6f8e9326c589437114a2339119eac76f Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Mon, 25 Nov 2019 09:34:20 -0600 Subject: [PATCH 7/7] add RsdToleranceFactor for linop --- .../syssolver_linop_clover_mg_proto_qphix_eo.cc | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/lib/actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.cc b/lib/actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.cc index dec3de74f..69436921c 100644 --- a/lib/actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.cc +++ b/lib/actions/ferm/invert/mg_proto/syssolver_linop_clover_mg_proto_qphix_eo.cc @@ -112,7 +112,7 @@ namespace Chroma LinOpSysSolverMGProtoQPhiXEOClover::operator()(T& psi, const T& chi) const { QDPIO::cout << "Jolly Greetings from Even-Odd Multigridland" << std::endl; - const Subset& s = A->subset(); + const Subset& s = A->subset(); StopWatch swatch; StopWatch swatch2; @@ -146,10 +146,10 @@ namespace Chroma MG::LinearSolverResults res=(*eo_solver)(qphix_out,qphix_in, RELATIVE); swatch2.stop(); - + double qphix_out_norm_cb0 = MG::Norm2Vec(qphix_out, SUBSET_EVEN); double qphix_out_norm_cb1 = MG::Norm2Vec(qphix_out, SUBSET_ODD); - + psi = zero; QPhiXSpinorToQDPSpinor(qphix_out,psi); Double psi_norm_cb0 = norm2(psi,rb[0]); @@ -170,8 +170,9 @@ namespace Chroma Double n2 = norm2(tmp,s); Double n2rel = n2 / norm2(chi,s); QDPIO::cout << "MG_PROTO_QPHIX_EO_CLOVER_INVERTER: iters = "<< res.n_count << " rel resid = " << sqrt(n2rel) << std::endl; - if( toBool( sqrt(n2rel) > invParam.OuterSolverRsdTarget ) ) { - MGSolverException convergence_fail(invParam.CloverParams.Mass, + if( toBool( sqrt(n2rel) > invParam.OuterSolverRsdTarget * invParam.RsdToleranceFactor) ) { + QDPIO::cout<<"MG_PROTO Solver Not Converged!"<