From 784ef33a1aa651880a0c99207ad44df1fd1d2dec Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Fri, 19 Apr 2019 18:09:34 -0400 Subject: [PATCH 1/6] refactor the multi-Hasenbusch monomial by using MRHS proxy solver --- ...hmc.seoprec_clover_multihasenbusch.ini.xml | 40 +-- ..._clover_multihasenbusch.qudaclover.ini.xml | 259 ++++++++++++++++++ 2 files changed, 280 insertions(+), 19 deletions(-) create mode 100644 tests/hmc/hmc.seoprec_clover_multihasenbusch.qudaclover.ini.xml diff --git a/tests/hmc/hmc.seoprec_clover_multihasenbusch.ini.xml b/tests/hmc/hmc.seoprec_clover_multihasenbusch.ini.xml index 8ed3de5f9..9e121080d 100644 --- a/tests/hmc/hmc.seoprec_clover_multihasenbusch.ini.xml +++ b/tests/hmc/hmc.seoprec_clover_multihasenbusch.ini.xml @@ -80,10 +80,8 @@ TWO_FLAVOR_SEOPREC_CONSTDET_RATIO_CONV_CONV_MULTIHASEN_FERM_MONOMIAL - 0 0.001 0.002 0.003 - 3 - - + 0.0 + SEOPREC_CLOVER -0.0640 1.58932722549812 @@ -106,22 +104,26 @@ - CG_INVERTER - - -0.0640 - 1.58932722549812 - 0.902783591772098 - - true - 3 - 4.3 - 1.265 - - - 1.0e-7 - 1000 + MULTI_RHS_SEOPREC_TWISTED_PROXY_INVERTER + 3 + 0.001 0.002 0.003 + + CG_INVERTER + + -0.0640 + 1.58932722549812 + 0.902783591772098 + + true + 3 + 4.3 + 1.265 + + + 1.0e-7 + 1000 + - HMC::has diff --git a/tests/hmc/hmc.seoprec_clover_multihasenbusch.qudaclover.ini.xml b/tests/hmc/hmc.seoprec_clover_multihasenbusch.qudaclover.ini.xml new file mode 100644 index 000000000..020c0cda6 --- /dev/null +++ b/tests/hmc/hmc.seoprec_clover_multihasenbusch.qudaclover.ini.xml @@ -0,0 +1,259 @@ + + + + + DISORDERED + DUMMY + false + + + + 2846 + 1259 + 980 + 108 + + + 0 + 0 + 5 + 5 + 5 + dummy_run_multi + SINGLEFILE + false + false + 20 + true + 1 + true + + + PLAQUETTE + 1 + + 2 + + SIMPLE_GAUGE_STATE + + PERIODIC_GAUGEBC + + + + + default_gauge_field + + + + + + 4 4 4 4 + + + N_FLAVOR_LOGDET_DIAG_FERM_MONOMIAL + + SEOPREC_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_SEOPREC_CONSTDET_RATIO_CONV_CONV_MULTIHASEN_FERM_MONOMIAL + 0.0 + + SEOPREC_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 + + + + + MULTI_RHS_SEOPREC_TWISTED_PROXY_INVERTER + 3 + 0.001 0.002 0.003 + + QUDA_CLOVER_INVERTER + + -0.0640 + 1.58932722549812 + 0.902783591772098 + 666 + + true + 3 + 4.3 + 1.265 + + + 1.0e-7 + 1.0e-1 + 0 + 1000 + CG + 10.0 + true + false + false + RECONS_12 + SINGLE + RECONS_12 + false + true + true + + + + HMC::has + + + + TWO_FLAVOR_SEOPREC_CONSTDET_MULTIHASEN_CANCEL_FERM_MONOMIAL + 0.003 + + SEOPREC_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 + + + + + CG_INVERTER + + -0.0640 + 1.58932722549812 + 0.902783591772098 + + true + 3 + 4.3 + 1.265 + + + 1.0e-7 + 1000 + + + HMC::has_cancel_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::has + HMC::has_cancel_2flav + HMC::gauge_s + HMC::gauge_t + + + + 0.5 + true + 3 + 3.5 + + LCM_STS_FORCE_GRAD + 10 + + HMC::logdet_2flav + HMC::has + HMC::gauge_s + HMC::has_cancel_2flav + HMC::gauge_t + + + + + From 5dba4ce323f1d05683012d6589cd102041d1cf3a Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Fri, 19 Apr 2019 18:15:01 -0400 Subject: [PATCH 2/6] refactor the multi-Hasenbusch monomial by using MRHS proxy solver --- .../invert/syssolver_mrhs_twisted_params.h | 1 + .../invert/syssolver_mrhs_twisted_proxy.h | 2 +- .../two_flavor_multihasen_cancel_monomial_w.h | 2 +- ..._conv_conv_multihasen_monomial_params_w.cc | 12 +- ...o_conv_conv_multihasen_monomial_params_w.h | 8 +- ...or_ratio_conv_conv_multihasen_monomial_w.h | 321 ++++++++++-------- other_libs/cpp_wilson_dslash | 2 +- 7 files changed, 186 insertions(+), 162 deletions(-) diff --git a/lib/actions/ferm/invert/syssolver_mrhs_twisted_params.h b/lib/actions/ferm/invert/syssolver_mrhs_twisted_params.h index 639b05ee7..0d49d867b 100644 --- a/lib/actions/ferm/invert/syssolver_mrhs_twisted_params.h +++ b/lib/actions/ferm/invert/syssolver_mrhs_twisted_params.h @@ -20,6 +20,7 @@ namespace Chroma { int BlockSize; // Number of right hand sides multi1d Twists; GroupXML_t SubInverterXML; // XML to create the sub-integrator from + SysSolverMRHSTwistedParams() = default; SysSolverMRHSTwistedParams(XMLReader& xml_in, const std::string& path); }; diff --git a/lib/actions/ferm/invert/syssolver_mrhs_twisted_proxy.h b/lib/actions/ferm/invert/syssolver_mrhs_twisted_proxy.h index 46502b951..608fd09b7 100644 --- a/lib/actions/ferm/invert/syssolver_mrhs_twisted_proxy.h +++ b/lib/actions/ferm/invert/syssolver_mrhs_twisted_proxy.h @@ -23,7 +23,7 @@ namespace Chroma { namespace MRHSUtils { - void resetTwist(XMLReader& xml, const Real& Twist) + static void resetTwist(XMLReader& xml, const Real& Twist) { std::string sub_solver_type; read(xml, "invType", sub_solver_type); diff --git a/lib/update/molecdyn/monomial/two_flavor_multihasen_cancel_monomial_w.h b/lib/update/molecdyn/monomial/two_flavor_multihasen_cancel_monomial_w.h index 85f278e9f..94874c8ba 100644 --- a/lib/update/molecdyn/monomial/two_flavor_multihasen_cancel_monomial_w.h +++ b/lib/update/molecdyn/monomial/two_flavor_multihasen_cancel_monomial_w.h @@ -60,7 +60,7 @@ namespace Chroma START_CODE(); XMLWriter& xml_out = TheXMLLogWriter::Instance(); - push(xml_out, "S"); + push(xml_out, "PrecConstDetTwoFlavorMultihasenCancelMonomial"); const FAType& FA = getFermAct(); Handle > state = FA.createState(s.getQ()); diff --git a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.cc b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.cc index ec97e737a..680c77e6c 100644 --- a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.cc +++ b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.cc @@ -13,9 +13,9 @@ namespace Chroma // Get the top of the parameter XML tree XMLReader paramtop(xml_in, path); try{ - read(paramtop, "Action", fermactInv); - read(paramtop, "ShiftedMass", mu); - read(paramtop, "NumofHasenTerms", numHasenTerms); + fermact = readXMLGroup(paramtop, "FermionAction", "FermAct"); + read(paramtop, "InvertParam", inverter); + read(paramtop, "BaseTwist", base_twist); if(paramtop.count("./ChronologicalPredictor") == 0){ predictor.xml = ""; }else{ @@ -37,9 +37,9 @@ namespace Chroma //! Write Parameters void write(XMLWriter& xml, const std::string& path, const TwoFlavorRatioConvConvMultihasenWilsonTypeFermMonomialParams& params){ - write(xml, "Action", params.fermactInv); - write(xml, "ShiftedMass", params.mu); - write(xml, "NumofHasenTerms", params.numHasenTerms); + write(xml, "FermionAction", params.fermact.id); + write(xml, "BaseTwist", params.base_twist); + write(xml, "NumofHasenTerms", params.inverter.BlockSize); } }// end namespace Chroma diff --git a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.h b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.h index a9dde92c7..9b424e9ad 100644 --- a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.h +++ b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_params_w.h @@ -6,6 +6,7 @@ #ifndef __TWO_FLAVOR_RATIO_CONV_CONV_MULTIHASEN_MONOMIAL_PARAMS_W_H__ #define __TWO_FLAVOR_RATIO_CONV_CONV_MULTIHASEN_MONOMIAL_PARAMS_W_H__ #include "update/molecdyn/monomial/comp_approx.h" +#include "actions/ferm/invert/syssolver_mrhs_twisted_params.h" namespace Chroma { @@ -19,9 +20,10 @@ namespace Chroma // Read monomial from some root path TwoFlavorRatioConvConvMultihasenWilsonTypeFermMonomialParams(XMLReader& in, const std::string& path); - CompActionInv_t fermactInv; /*!< Fermion action and invert params */ - GroupXML_t predictor; /*!< The Chrono Predictor XML */ - multi1d mu; /*!< Shifted mass term */ + SysSolverMRHSTwistedParams inverter; /*!< Multi-right hand solver parameter */ + GroupXML_t fermact; /*!< Fermion action */ + GroupXML_t predictor; /*!< The Chrono Predictor XML */ + Real base_twist; /*!< Shifted mass of base operator */ int numHasenTerms; /*!< Number of Hasenbusch terms */ }; diff --git a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h index 19e45d270..72e7ae732 100644 --- a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h +++ b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h @@ -23,6 +23,11 @@ #include "eoprec_logdet_wilstype_fermact_w.h" #include "update/molecdyn/predictor/chrono_predictor_factory.h" #include "update/molecdyn/predictor/zero_guess_predictor.h" +#include "actions/ferm/invert/syssolver_mrhs_twisted_params.h" +#include "actions/ferm/invert/syssolver_mrhs_twisted_proxy.h" +#include "actions/ferm/invert/syssolver_linop_mrhs_factory.h" +#include "actions/ferm/invert/syssolver_mdagm_mrhs_factory.h" +#include "actions/ferm/linop/multi_twist_linop_w.h" namespace Chroma { @@ -59,46 +64,57 @@ namespace Chroma // Fermion action const FAType& FA = getFermAct(); + const int N = numHasenTerms; Double S=0; // Get the X fields - T X; + multi1d X(N); Handle > state(FA.createState(s.getQ())); Handle > base_op(FA.linOp(state)); - - for(int i=0; i > - M(new TwistedShiftedLinOp(*base_op, mu[i])); - Handle > - M_prec(new TwistedShiftedLinOp(*base_op, mu[i+1])); - - // Action calc doesnt use chrono predictor use zero guess - X[M->subset()] = zero; - // Reset chrono predictor - QDPIO::cout << - "TwoFlavorRatioConvConvMultihasenWilson4DMonomial: resetting Predictor before energy calc solve" << std::endl; - (getMDSolutionPredictor()).reset(); - // Get system solver - const GroupXML_t& invParam = getInvParams(); - std::istringstream xml(invParam.xml); - XMLReader paramtop(xml); - Handle > - invMdagM(TheMdagMFermSystemSolverFactory::Instance().createObject( - invParam.id, paramtop, invParam.path, state, M)); - // M_dag_prec phi = M^{dag}_prec \phi - the RHS - T M_dag_prec_phi; - (*M_prec)(M_dag_prec_phi, getPhi(i), MINUS); - - // Solve MdagM X = eta - SystemSolverResults_t res = (*invMdagM)(X, M_dag_prec_phi); - - T phi_tmp = zero; - (*M_prec)(phi_tmp, X, PLUS); - Double action = innerProductReal(getPhi(i), phi_tmp, M->subset()); + + // Construct Linop with multi-twist + // Denominator of multi-Hasenbusch monomial + Handle > + M_den(new MultiTwistLinOp(base_op, mu_den)); + // Numerator of multi-Hasenbusch monomial, not needed here, + // use fermion action directly. + //Handle > + // M_num(new MultiTwistLinOp(base_op, mu_num)); + + // Action calc doesnt use chrono predictor use zero guess + for(int i=0; isubset()] = zero; + } + // Reset chrono predictor + QDPIO::cout << + "TwoFlavorRatioConvConvMultihasenWilson4DMonomial: resetting Predictor before energy calc solve" << std::endl; + (getMDSolutionPredictor()).reset(); + + // A bit tedious work to set up the Twists for numerator + // MdagM solver + SysSolverMRHSTwistedParams paramnum = invParam; + paramnum.Twists = mu_num; + // Set up the MRHS solver for MdagM, where M is the + // twisted linop on numerator. + MdagMMRHSSysSolverTwistedProxy invMdagM(paramnum, FA, state); + + // M_dag_den_phi = M^{dag}_den \phi - the RHS + multi1d M_dag_den_phi(N); + (*M_den)(M_dag_den_phi, phi, MINUS); + + // Solve MdagM X = eta + SystemSolverResultsMRHS_t res = invMdagM(X, M_dag_den_phi); + multi1d phi_tmp(N); + phi_tmp = zero; + (*M_den)(phi_tmp, X, PLUS); + + for(int i=0; isubset()); // Write out inversion number and action for every hasenbusch term std::string n_count = "n_count_hasenterm_" + std::to_string(i+1); std::string s_action = "S_hasenterm_" + std::to_string(i+1); - write(xml_out, n_count, res.n_count); + write(xml_out, n_count, res.n_count[i]); write(xml_out, s_action, action); S += action; } @@ -110,87 +126,82 @@ namespace Chroma } // Total force of multi-hasen terms - virtual void dsdq(P& F, const AbsFieldState& s) + virtual void dsdq(P& F_t, const AbsFieldState& s) { START_CODE(); XMLWriter& xml_out = TheXMLLogWriter::Instance(); push(xml_out, "TwoFlavorRatioConvConvMultihasenWilsonTypeFermMonomial"); - P F_t; P F_tmp; F_t.resize(Nd); F_tmp.resize(Nd); - F.resize(Nd); for(int i=0; i& FA = getFermAct(); - // X = (M^\dag M)^(-1) M_prec^\dag \phi + // X = (M^\dag M)^(-1) M_den^\dag \phi // Y = M X - T X = zero; - T Y = zero; - // M_dag_prec_phi = M^\dag_prec \phi - T M_dag_prec_phi; + multi1d X(N), Y(N); + for(int i=0; i M_dag_den_phi(N); Handle > state(FA.createState(s.getQ())); Handle > base_op(FA.linOp(state)); - for(int i=0; i > - M(new TwistedShiftedLinOp(*base_op, mu[i])); - Handle > - M_prec(new TwistedShiftedLinOp(*base_op, mu[i+1])); - // Get system solver - const GroupXML_t& invParam = getInvParams(); - std::istringstream xml(invParam.xml); - XMLReader paramtop(xml); - Handle > - invMdagM(TheMdagMFermSystemSolverFactory::Instance().createObject( - invParam.id,paramtop, invParam.path, state, M)); - - (*M_prec)(M_dag_prec_phi, getPhi(i), MINUS); - SystemSolverResults_t res = (*invMdagM)(X, M_dag_prec_phi, getMDSolutionPredictor()); - (*M)(Y, X, PLUS); - - // cast linearop to difflinearop - const DiffLinearOperator& diffM = dynamic_cast&>(*M); - const DiffLinearOperator& diffM_prec = dynamic_cast&>(*M_prec); - - // deriv part 1: \phi^\dag \dot(M_prec) X - diffM_prec.deriv(F_t, getPhi(i), X, PLUS); - - // deriv part 2: - X^\dag \dot(M^\dag) Y - diffM.deriv(F_tmp, X, Y, MINUS); - F_t -= F_tmp; - - // deriv part 3: - Y^\dag \dot(M) X - diffM.deriv(F_tmp, Y, X, PLUS); - F_t -= F_tmp; - - // deriv part 4: X^\dag \dot(M_prec)^\dag \phi - diffM_prec.deriv(F_tmp, X, getPhi(i), MINUS); - F_t += F_tmp; - - // total force from all Hasenbusch terms - F += F_t; - - // F now holds derivative with respect to possibly fat links - // now derive it with respect to the thin links if needs be - state->deriv(F); - - // Write out inversion number and action for every hasenbusch term + // Construct Linop with multi-twist + // Denominator of multi-Hasenbusch monomial + Handle > + M_den(new MultiTwistLinOp(base_op, mu_den)); + // Numerator of multi-Hasenbusch monomial + Handle > + M_num(new MultiTwistLinOp(base_op, mu_num)); + + // A bit tedious work to set up the Twists for numerator + // MdagM solver + SysSolverMRHSTwistedParams paramnum = invParam; + paramnum.Twists = mu_num; + // Set up the MRHS solver for MdagM, where M is the + // twisted linop on numerator. + MdagMMRHSSysSolverTwistedProxy invMdagM(paramnum, FA, state); + + (*M_den)(M_dag_den_phi, phi, MINUS); + SystemSolverResultsMRHS_t res = invMdagM(X, M_dag_den_phi); + (*M_num)(Y, X, PLUS); + + // deriv part 1: \phi^\dag \dot(M_den) X + M_den->deriv(F_t, phi, X, PLUS); + + // deriv part 2: - X^\dag \dot(M_num^\dag) Y + M_num->deriv(F_tmp, X, Y, MINUS); + F_t -= F_tmp; + + // deriv part 3: - Y^\dag \dot(M_num) X + M_num->deriv(F_tmp, Y, X, PLUS); + F_t -= F_tmp; + + // deriv part 4: X^\dag \dot(M_prec)^\dag \phi + M_den->deriv(F_tmp, X, phi, MINUS); + F_t += F_tmp; + + // F now holds derivative with respect to possibly fat links + // now derive it with respect to the thin links if needs be + state->deriv(F_t); + + for(int i=0; i& FA = getFermAct(); - + const int N = numHasenTerms; + Handle > state(FA.createState(s.getQ())); Handle > base_op(FA.linOp(state)); - for(int i=0; i > - M(new TwistedShiftedLinOp(*base_op, mu[i])); - Handle > - M_prec(new TwistedShiftedLinOp(*base_op, mu[i+1])); - T eta = zero; - gaussian(eta, M->subset()); - - // Account for fermion BC by modifying the proposed field - FA.getFermBC().modifyF(eta); - - eta *= sqrt(0.5); - // Now we have to get the refreshment right: - // We have \eta^{\dag} \eta = \phi M_prec (M^dag M )^-1 M^dag_prec \phi - // so that \eta = (M^\dag)^{-1} M^{dag}_prec \phi - // So need to solve M^{dag}_prec \phi = M^{\dag} \eta - // Which we can solve as - // M^{dag}_prec M_prec ( M_prec^{-1} ) \phi = M^{\dag} \eta - // - // Zero out everything - to make sure there is no junk - // in uninitialised places - T eta_tmp = zero; - T phi_tmp = zero; + // Construct Linop with multi-twist + // Denominator of multi-Hasenbusch monomial + Handle > + M_den(new MultiTwistLinOp(base_op, mu_den)); + // Numerator of multi-Hasenbusch monomial + Handle > + M_num(new MultiTwistLinOp(base_op, mu_num)); + + multi1d eta(N); + for(int i=0; isubset()); + FA.getFermBC().modifyF(eta[i]); + eta[i] *= sqrt(0.5); + } + // Now we have to get the refreshment right: + // We have \eta^{\dag} \eta = \phi^{\dag} M_den (M_num^dag M_num )^-1 M^dag_den \phi + // so that \eta = (M_num^\dag)^{-1} M^{dag}_den \phi + // So need to solve M^{dag}_den \phi = M_num^{\dag} \eta + // Which we can solve as + // M^{dag}_den M_den ( M_den^{-1} ) \phi = M_num^{\dag} \eta + // + // Zero out everything - to make sure there is no junk + // in uninitialised places + multi1d eta_tmp(N), phi_tmp(N); + for(int i=0; i > - invMdagM(TheMdagMFermSystemSolverFactory::Instance().createObject( - invParam.id,paramtop, invParam.path, state, M_prec)); - - // Solve MdagM_prec X = eta - SystemSolverResults_t res = (*invMdagM)(phi_tmp, eta_tmp); - (*M_prec)(getPhi(i), phi_tmp, PLUS); // (Now get phi = M_prec (M_prec^{-1}\phi) - QDPIO::cout<<"TwoFlavRatioConvConvMultihasenWilson4DMonomial: resetting Predictor at end of field refresh"< invMdagM(invParam, FA, state); + // Solve MdagM_den X = eta + SystemSolverResultsMRHS_t res = invMdagM(phi_tmp, eta_tmp); + (*M_den)(phi, phi_tmp, PLUS); // (Now get phi = M_den (M_den^{-1}\phi) + QDPIO::cout<<"TwoFlavRatioConvConvMultihasenWilson4DMonomial: resetting Predictor at end of field refresh"< > fermact; - // Shifted mass mu - multi1d mu; - - // Number of Hasenbusch terms + // Shifted mass from denominator of multi-Hasenbusch term + multi1d mu_den; + // Shifted mass from numerator of multi-Hasenbusch term + multi1d mu_num; + // No. of Hasenbusch terms int numHasenTerms; - // The parameters for the inversion - GroupXML_t invParam; + SysSolverMRHSTwistedParams invParam; Handle > chrono_predictor; }; @@ -333,23 +349,23 @@ namespace Chroma START_CODE(); QDPIO::cout << "Constructor: " << __func__ << std::endl; - if( param.fermactInv.invParam.id == "NULL" ) { + if( param.fermact.id == "NULL" ) { QDPIO::cerr << "Fermact inverter params are NULL" << std::endl; QDP_abort(1); } - invParam = param.fermactInv.invParam; + invParam = param.inverter; // Fermion action (with original seoprec action) { - std::istringstream is(param.fermactInv.fermact.xml); + std::istringstream is(param.fermact.xml); XMLReader fermact_reader(is); - QDPIO::cout << "Construct fermion action= " << param.fermactInv.fermact.id << std::endl; + QDPIO::cout << "Construct fermion action= " << param.fermact.id << std::endl; WilsonTypeFermAct* tmp_act = - TheWilsonTypeFermActFactory::Instance().createObject(param.fermactInv.fermact.id, + TheWilsonTypeFermActFactory::Instance().createObject(param.fermact.id, fermact_reader, - param.fermactInv.fermact.path); + param.fermact.path); FAType* downcast=dynamic_cast* >(tmp_act); @@ -363,12 +379,17 @@ namespace Chroma } // Number of Hasenbusch term - numHasenTerms = param.numHasenTerms; - - // Shifted mass - mu.resize(numHasenTerms+1); - for(int i=0; i< numHasenTerms+1; ++i){ - mu[i] = param.mu[i]; + numHasenTerms = param.inverter.BlockSize; + // Shifted mass of denominator + mu_den.resize(numHasenTerms); + for(int i=0; i< numHasenTerms; ++i){ + mu_den[i] = param.inverter.Twists[i]; + } + // Shifted mass of numerator + mu_num.resize(numHasenTerms); + mu_num[0] = param.base_twist; + for(int i=1; i< numHasenTerms; ++i){ + mu_num[i] = param.inverter.Twists[i-1]; } // Pseudofermion field phi diff --git a/other_libs/cpp_wilson_dslash b/other_libs/cpp_wilson_dslash index 172632cef..56a4abf64 160000 --- a/other_libs/cpp_wilson_dslash +++ b/other_libs/cpp_wilson_dslash @@ -1 +1 @@ -Subproject commit 172632cef35898d6958782ab8532f832e054577b +Subproject commit 56a4abf64c4d586a73bc2fd7747067f6e92329c9 From 5b4efb3b640dd58ad6a3255fb3c72f57f1b57ff7 Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Fri, 23 Aug 2019 14:04:43 -0400 Subject: [PATCH 3/6] match the multi level definition of some variables in Quda multigrid --- .../ferm/invert/quda_solvers/quda_mg_utils.h | 9 +++++++-- .../syssolver_linop_wilson_quda_multigrid_w.h | 5 +++-- .../syssolver_mdagm_clover_quda_multigrid_w.cc | 16 ++++++++-------- .../syssolver_mdagm_clover_quda_multigrid_w.h | 11 ++++++++++- 4 files changed, 28 insertions(+), 13 deletions(-) diff --git a/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h b/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h index f7968635b..324ba5986 100644 --- a/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h +++ b/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h @@ -219,9 +219,13 @@ namespace Chroma { if (ip.check_multigrid_setup == true ) { mg_param.run_verify = QUDA_BOOLEAN_YES; + mg_param.run_low_mode_check = QUDA_BOOLEAN_NO; + mg_param.run_oblique_proj_check = QUDA_BOOLEAN_NO; } else { mg_param.run_verify = QUDA_BOOLEAN_NO; + mg_param.run_low_mode_check = QUDA_BOOLEAN_NO; + mg_param.run_oblique_proj_check = QUDA_BOOLEAN_NO; } for (int i=0; img_param).vec_outfile, subspace_prefix.c_str(), 256); + std::strncpy((subspace_pointers->mg_param).vec_outfile[0], subspace_prefix.c_str(), 256); // if the source string is too long result will not be null terminated, so null terminate in that case - if( subspace_prefix.size() > 255 ) { (subspace_pointers->mg_param).vec_outfile[255] = '\0'; } - QDPInternal::broadcast( (void*)(subspace_pointers->mg_param).vec_outfile, 256); + if( subspace_prefix.size() > 255 ) { (subspace_pointers->mg_param).vec_outfile[0][255] = '\0'; } + QDPInternal::broadcast( (void*)(subspace_pointers->mg_param).vec_outfile[0], 256); #ifdef QUDA_MG_DUMP_ENABLED dumpMultigridQuda(subspace_pointers->preconditioner, &(subspace_pointers->mg_param)); #endif - (subspace_pointers->mg_param).vec_outfile[0]='\0'; + (subspace_pointers->mg_param).vec_outfile[0][0]='\0'; } } @@ -324,15 +324,15 @@ namespace Chroma std::string subspace_prefix = file_prefix + "_subspace"; // Up to the length of the buffer (256) padded with zeros - std::strncpy((subspace_pointers->mg_param).vec_outfile, subspace_prefix.c_str(), 256); + std::strncpy((subspace_pointers->mg_param).vec_outfile[0], subspace_prefix.c_str(), 256); // If source string is too long it will be truncated and not null terminated, so null terminate - if( subspace_prefix.size() > 255 ) { (subspace_pointers->mg_param).vec_outfile[255] = '\0'; } - QDPInternal::broadcast( (void*)(subspace_pointers->mg_param).vec_outfile, 256); + if( subspace_prefix.size() > 255 ) { (subspace_pointers->mg_param).vec_outfile[0][255] = '\0'; } + QDPInternal::broadcast( (void*)(subspace_pointers->mg_param).vec_outfile[0], 256); #ifdef QUDA_MG_DUMP_ENABLED dumpMultigridQuda(subspace_pointers->preconditioner, &(subspace_pointers->mg_param)); #endif - (subspace_pointers->mg_param).vec_outfile[0] ='\0'; + (subspace_pointers->mg_param).vec_outfile[0][0] ='\0'; } } diff --git a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h index a877bcac5..ad3ec0b1a 100644 --- a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h +++ b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h @@ -481,7 +481,16 @@ class MdagMSysSolverQUDAMULTIGRIDClover : public MdagMSystemSolver Date: Mon, 9 Sep 2019 18:02:29 -0400 Subject: [PATCH 4/6] add Quda multigrid support for symmetric preconditioned twisted-Hasenbusch operator, mainly used in HMC to solve the MdagM psi = chi system, not support asymmetric preconditioning currently --- ...syssolver_mdagm_clover_quda_multigrid_w.cc | 2 +- .../syssolver_mdagm_clover_quda_multigrid_w.h | 45 ++- ...over_multihasenbusch.qudamultigrid.ini.xml | 365 ++++++++++++++++++ 3 files changed, 407 insertions(+), 5 deletions(-) create mode 100644 tests/hmc/hmc.seoprec_clover_multihasenbusch.qudamultigrid.ini.xml diff --git a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.cc b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.cc index 0340f275a..b81272d07 100644 --- a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.cc +++ b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.cc @@ -87,7 +87,7 @@ namespace Chroma // // So M x = b => A_oo (M_symm) x = b // => M_symm x = A^{-1}_oo b = chi_mod - invclov.apply(mod_chi, chi_s, PLUS, 1); + invclov.apply(mod_chi, chi_s, PLUS, 1); } else { // If we work with symmetric preconditioning nothing else needs done diff --git a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h index ad3ec0b1a..879592cce 100644 --- a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h +++ b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h @@ -654,17 +654,42 @@ QDPIO::cout << solver_string << " init_time = " // Y solve at 0.5 * Target Residuum -- Evan's bound quda_inv_param.tol = toDouble(Real(0.5)*invParam.RsdTarget); + + // For twisted-Hasenbusch operator, we have + // asymmetric: M_a(mu) = A_oo - kappa^2 D_oe A^{-1}_ee D_eo - i*mu*gamma5 + // symmetric: M_s(mu) = 1 - kappa^2 A^{-1}_oo D_oe A^{-1}_ee D_eo - i*mu*gamma5*A_oo + // + // NOTE: both M_a(mu) and M_s(mu) are not gamma5 hermitian, + // and current Quda just support the multigrid solver for symmetric preconditioned + // linear operator M, not Mdag + // + // To solve MdagM X = chi, use two step solver + // 1). solve Mdag Y = chi with M_s^\dagger(mu) = gamma5 A_oo M_s(-mu) A^{-1}_oo gamma5, + // therefore we need to solve the M_s(-mu) Y' = chi' system just by flip the sign of twisted term + // in Quda's multigrid parameter + // 2). Y = A_oo gamma5 Y', then flip back to solve M_s(mu) X = Y + // + if( invParam.asymmetricP == true ) { - res1 = qudaInvert(*clov, + if(quda_inv_param.dslash_type = QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH){ + QDPIO::cout<<"Multigrid solver for asymmetric preconditioned twisted-Hasenbusch operator not yet implemented!\n"; + QDP_abort(1); + }else{ + res1 = qudaInvert(*clov, *invclov, g5chi, Y_prime); - Y[rb[1]] = Gamma(Nd*Nd -1)*Y_prime; + Y[rb[1]] = Gamma(Nd*Nd -1)*Y_prime; + } } else { T tmp = zero; invclov->apply(tmp,g5chi,MINUS,1); + if(quda_inv_param.dslash_type = QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH) + // flip the sign of twisted term to solve the Mdag Y = chi' system + quda_inv_param.mu = 2.0 * toDouble(invParam.CloverParams.twisted_m); + res1 = qudaInvert(*clov, *invclov, tmp, @@ -839,6 +864,9 @@ QDPIO::cout << solver_string << " init_time = " quda_inv_param.tol = toDouble(invParam.RsdTarget); X_solve_timer.start(); // Solve for psi + if(quda_inv_param.dslash_type = QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH) + // flip back the twisted term to solve M X = Y system + quda_inv_param.mu = -2.0 * toDouble(invParam.CloverParams.twisted_m); res2 = qudaInvert(*clov, *invclov, Y, @@ -1078,16 +1106,23 @@ QDPIO::cout << solver_string << " init_time = " T g5chi = zero; g5chi[rb[1]]= Gamma(Nd*Nd-1)*chi; if( invParam.asymmetricP == true ) { - res1 = qudaInvert(*clov, + if(quda_inv_param.dslash_type = QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH){ + QDPIO::cout<<"Multigrid solver for asymmetric preconditioned twisted-Hasenbusch operator not yet implemented!\n"; + QDP_abort(1); + }else{ + res1 = qudaInvert(*clov, *invclov, g5chi, Y_prime); - Y[rb[1]] = Gamma(Nd*Nd -1)*Y_prime; + Y[rb[1]] = Gamma(Nd*Nd -1)*Y_prime; + } } else { T tmp = zero; invclov->apply(tmp,g5chi,MINUS,1); + if(quda_inv_param.dslash_type = QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH) + quda_inv_param.mu = 2.0 * toDouble(invParam.CloverParams.twisted_m); res1 = qudaInvert(*clov, *invclov, tmp, @@ -1274,6 +1309,8 @@ QDPIO::cout << solver_string << " init_time = " //psi[A->subset()]=zero; psi = zero; // Solve for psi + if(quda_inv_param.dslash_type = QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH) + quda_inv_param.mu = -2.0 * toDouble(invParam.CloverParams.twisted_m); res2 = qudaInvert(*clov, *invclov, Y, diff --git a/tests/hmc/hmc.seoprec_clover_multihasenbusch.qudamultigrid.ini.xml b/tests/hmc/hmc.seoprec_clover_multihasenbusch.qudamultigrid.ini.xml new file mode 100644 index 000000000..e73c3a52e --- /dev/null +++ b/tests/hmc/hmc.seoprec_clover_multihasenbusch.qudamultigrid.ini.xml @@ -0,0 +1,365 @@ + + + + + DISORDERED + DUMMY + false + + + + 2846 + 1259 + 980 + 108 + + + 0 + 0 + 5 + 1 + 5 + dummy_run_multi + SINGLEFILE + false + false + 20 + true + 1 + true + + + PLAQUETTE + 1 + + 2 + + SIMPLE_GAUGE_STATE + + PERIODIC_GAUGEBC + + + + + default_gauge_field + + + + ERASE_QUDA_MULTIGRID_SUBSPACE + mgsub + + + + + 8 8 8 8 + + + N_FLAVOR_LOGDET_DIAG_FERM_MONOMIAL + + SEOPREC_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_SEOPREC_CONSTDET_RATIO_CONV_CONV_MULTIHASEN_FERM_MONOMIAL + 0.0 + + SEOPREC_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 + + + + + MULTI_RHS_SEOPREC_TWISTED_PROXY_INVERTER + 3 + 0.01 0.02 0.03 + + QUDA_MULTIGRID_CLOVER_INVERTER + + 1.0 + true + 1.0e-1 + 10 + true + HALF + RECONS_12 + 24 32 + true + true + MG_RECURSIVE + 0 0 + 8 8 + ADDITIVE_SCHWARZ + + 2 2 2 2 + 2 2 2 2 + + + GCR + CA_GCR + + 0.1 0.1 0.1 + 12 12 8 + 1.0 1.0 1.0 + + CA_GCR + CA_GCR + CA_GCR + + 0.25 0.25 0.25 + + CG + CG + + 5e-06 5e-06 + 500 500 + 500 500 + 20 + 10 + + mgsub + 500 + 100 + + -0.0640 + 666 + 1.58932722549812 + 0.902783591772098 + + 3 + 4.3 + 1.265 + true + + + 1e-8 + 0.1 + 10 + false + true + GCR + true + false + RECONS_12 + SINGLE + RECONS_12 + false + true + true + + + + HMC::has + + + + TWO_FLAVOR_SEOPREC_CONSTDET_FERM_MONOMIAL + + SEOPREC_CLOVER + -0.0640 + 0.03 + 1.58932722549812 + 0.902783591772098 + + true + 3 + 4.3 + 1.265 + + + STOUT_FERM_STATE + 0.14 + 2 + 3 + + SIMPLE_FERMBC + 1 1 1 -1 + + + + + QUDA_MULTIGRID_CLOVER_INVERTER + + 1.0 + true + 1.0e-1 + 10 + true + HALF + RECONS_12 + 24 32 + true + true + MG_RECURSIVE + 0 0 + 8 8 + ADDITIVE_SCHWARZ + + 2 2 2 2 + 2 2 2 2 + + + GCR + CA_GCR + + 0.1 0.1 0.1 + 12 12 8 + 1.0 1.0 1.0 + + CA_GCR + CA_GCR + CA_GCR + + 0.25 0.25 0.25 + + CG + CG + + 5e-06 5e-06 + 500 500 + 500 500 + 20 + 10 + + mgsub + 500 + 100 + + -0.0640 + 0.03 + 1.58932722549812 + 0.902783591772098 + + 3 + 4.3 + 1.265 + true + + + 1e-8 + 0.1 + 10 + false + true + GCR + true + false + RECONS_12 + SINGLE + RECONS_12 + false + true + true + + + HMC::has_cancel_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::has + HMC::has_cancel_2flav + HMC::gauge_s + HMC::gauge_t + + + + 0.5 + true + 3 + 3.5 + + LCM_STS_FORCE_GRAD + 4 + + HMC::logdet_2flav + HMC::has + HMC::gauge_s + HMC::has_cancel_2flav + HMC::gauge_t + + + + + From 8e1c1a4b670f287d5efccba8a3d8c458b6cb6eb3 Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Sun, 24 Nov 2019 13:47:16 -0500 Subject: [PATCH 5/6] fix MG subspace creation bug for twisted operator, print out force from different twisted ratio term in the log file --- .../ferm/invert/quda_solvers/quda_mg_utils.h | 15 ++- .../quda_solvers/quda_multigrid_params.cc | 3 + .../quda_solvers/quda_multigrid_params.h | 7 +- .../syssolver_mdagm_clover_quda_multigrid_w.h | 112 +++++++++++------- .../syssolver_quda_multigrid_clover_params.h | 5 +- lib/actions/ferm/linop/multi_twist_linop_w.h | 13 ++ ...or_ratio_conv_conv_multihasen_monomial_w.h | 39 +++--- 7 files changed, 129 insertions(+), 65 deletions(-) diff --git a/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h b/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h index 324ba5986..034088049 100644 --- a/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h +++ b/lib/actions/ferm/invert/quda_solvers/quda_mg_utils.h @@ -116,8 +116,12 @@ namespace Chroma { QDPIO::cout<<"Creating MG subspace."<(paramtop, "MaxCoarseIterations", maxIterations, 12 ); + mu_factor.resize(mg_levels); + readArray(paramtop, "MuFactor", mu_factor, 1.0); + smootherType.resize(mg_levels); readArray(paramtop, "SmootherType", smootherType, MR); diff --git a/lib/actions/ferm/invert/quda_solvers/quda_multigrid_params.h b/lib/actions/ferm/invert/quda_solvers/quda_multigrid_params.h index e95a0aeed..361ee401a 100644 --- a/lib/actions/ferm/invert/quda_solvers/quda_multigrid_params.h +++ b/lib/actions/ferm/invert/quda_solvers/quda_multigrid_params.h @@ -43,7 +43,8 @@ namespace Chroma multi1d maxIterSubspaceCreate; multi1d rsdTargetSubspaceCreate; multi1d maxIterSubspaceRefresh; - + multi1d mu_factor; + MULTIGRIDSolverParams(XMLReader& xml, const std::string& path); MULTIGRIDSolverParams() { @@ -73,12 +74,13 @@ namespace Chroma tol.resize(mg_levels); maxIterations.resize(mg_levels); coarseSolverType.resize(mg_levels); - + smootherType.resize(mg_levels); smootherTol.resize(mg_levels); relaxationOmegaMG.resize(mg_levels); smootherSchwarzType.resize(mg_levels); smootherSchwarzCycle.resize(mg_levels); + mu_factor.resize(mg_levels); generate_nullspace = true; for(int l = 0; l < mg_levels - 1; l++) { @@ -103,6 +105,7 @@ namespace Chroma relaxationOmegaMG[l] = 0.85; smootherSchwarzType[l] = INVALID_SCHWARZ; smootherSchwarzCycle[l] = 1; + mu_factor[l] = 1.0; } outer_gcr_nkrylov = 12; precond_gcr_nkrylov = 12; diff --git a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h index 879592cce..f7fcf3e10 100644 --- a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h +++ b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h @@ -84,7 +84,8 @@ class MdagMSysSolverQUDAMULTIGRIDClover : public MdagMSystemSolver { ds_u += ds_tmp; } } + // return different component of force term + void deriv(multi1d

& ds_u, const multi1d& chi, const multi1d& psi, + enum PlusMinus isign) + { + assertArraySizes(chi, psi, _N); + for(int i=0; i<_N; ++i){ + ds_u[i].resize(Nd); + ds_u[i]= zero; + + TwistedShiftedLinOp theShiftedOp((*_theLinOp), _Twists[i]); + theShiftedOp.deriv(ds_u[i], chi[i], psi[i], isign); + } + } int size() const { diff --git a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h index 72e7ae732..0a2bf65cb 100644 --- a/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h +++ b/lib/update/molecdyn/monomial/two_flavor_ratio_conv_conv_multihasen_monomial_w.h @@ -132,16 +132,22 @@ namespace Chroma XMLWriter& xml_out = TheXMLLogWriter::Instance(); push(xml_out, "TwoFlavorRatioConvConvMultihasenWilsonTypeFermMonomial"); - P F_tmp; F_t.resize(Nd); - F_tmp.resize(Nd); for(int i=0; i F_tmp1(N), F_tmp2(N); + for(int i=0; i& FA = getFermAct(); // X = (M^\dag M)^(-1) M_den^\dag \phi // Y = M X @@ -150,7 +156,7 @@ namespace Chroma X[i] = zero; Y[i] = zero; } - // M_dag_den_phi = M^\dag_den \phi + // M_dag_den_phi = M^\dag_den \phi multi1d M_dag_den_phi(N); Handle > state(FA.createState(s.getQ())); @@ -177,20 +183,24 @@ namespace Chroma (*M_num)(Y, X, PLUS); // deriv part 1: \phi^\dag \dot(M_den) X - M_den->deriv(F_t, phi, X, PLUS); + M_den->deriv(F_tmp1, phi, X, PLUS); // deriv part 2: - X^\dag \dot(M_num^\dag) Y - M_num->deriv(F_tmp, X, Y, MINUS); - F_t -= F_tmp; + M_num->deriv(F_tmp2, X, Y, MINUS); + F_tmp1 -= F_tmp2; // deriv part 3: - Y^\dag \dot(M_num) X - M_num->deriv(F_tmp, Y, X, PLUS); - F_t -= F_tmp; + M_num->deriv(F_tmp2, Y, X, PLUS); + F_tmp1 -= F_tmp2; // deriv part 4: X^\dag \dot(M_prec)^\dag \phi - M_den->deriv(F_tmp, X, phi, MINUS); - F_t += F_tmp; - + M_den->deriv(F_tmp2, X, phi, MINUS); + F_tmp1 += F_tmp2; + + // Total force + for(int i=0; ideriv(F_t); @@ -199,9 +209,10 @@ namespace Chroma // Write out inversion number for every hasenbusch term std::string n_count = "n_count_hasenterm_" + std::to_string(i+1); write(xml_out, n_count, res.n_count[i]); + monitorForces(xml_out, "Forces_hasenterm_" + std::to_string(i+1), F_tmp1[i]); } // Total force from all Hasenbusch terms - monitorForces(xml_out, "Forces", F_t); + monitorForces(xml_out, "ForcesTotal", F_t); pop(xml_out); END_CODE(); } From 976452932ed900a02af40f5037bbdac01a5469aa Mon Sep 17 00:00:00 2001 From: Wei Sun Date: Mon, 9 Dec 2019 10:15:46 -0500 Subject: [PATCH 6/6] update MG subspace of X solver after every Y solver for twisted operator --- .../syssolver_mdagm_clover_quda_multigrid_w.h | 54 +++++++++++++++---- 1 file changed, 44 insertions(+), 10 deletions(-) diff --git a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h index f7fcf3e10..5fcd98036 100644 --- a/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h +++ b/lib/actions/ferm/invert/quda_solvers/syssolver_mdagm_clover_quda_multigrid_w.h @@ -504,7 +504,11 @@ if(TheNamedObjMap::Instance().check(invParam.SaveSubspaceID)) // Subspace ID exists add it to mg_state QDPIO::cout<< solver_string <<"Recovering subspace..."<(invParam.SaveSubspaceID); - for(int j=0; j < ip.mg_levels-1;++j) { + // Reset the twist parameter to current value + if(quda_inv_param.dslash_type == QUDA_CLOVER_HASENBUSCH_TWIST_DSLASH){ + (subspace_pointers->mg_inv_param).mu = toDouble(invParam.twist); + } + for(int j=0; j < ip.mg_levels-1;++j) { (subspace_pointers->mg_param).setup_maxiter_refresh[j] = 0; } updateMultigridQuda(subspace_pointers->preconditioner, &(subspace_pointers->mg_param)); @@ -865,7 +869,7 @@ QDPIO::cout << solver_string << " init_time = " Y_predictor_add_timer.start(); predictor.newYVector(Y); Y_predictor_add_timer.stop(); - + X_prediction_timer.start(); // Can predict psi in the usual way without reference to Y predictor.predictX(psi, (*MdagM), chi); @@ -880,8 +884,24 @@ QDPIO::cout << solver_string << " init_time = " { quda_inv_param.mu = -2.0 * toDouble(invParam.CloverParams.twisted_m); invParam.twist = quda_inv_param.mu; + + StopWatch subspaceX; + subspaceX.reset(); + subspaceX.start(); + + for(int i=0; img_param).setup_maxiter_refresh[i] = 0; + } + (subspace_pointers->mg_inv_param).mu = toDouble(invParam.twist); + updateMultigridQuda(subspace_pointers->preconditioner, &(subspace_pointers->mg_param)); + quda_inv_param.preconditioner = subspace_pointers->preconditioner; + + subspaceX.stop(); + QDPIO::cout<