CHROMA
one_flavor_ratio_rat_conv_monomial_w.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 
3 /*! @file
4  * @brief One flavor ratio of rational monomials using RHMC
5  */
6 
7 #ifndef __one_flavor_ratio_rat_conv_monomial_w_h__
8 #define __one_flavor_ratio_rat_conv_monomial_w_h__
9 
15 #include <typeinfo>
16 
17 namespace Chroma
18 {
19  //-------------------------------------------------------------------------------------------
20  //! Exact 1 flavor fermact monomial using rational polynomials
21  /*! @ingroup monomial
22  *
23  * Exact 1 flavor fermact monomial using Rational Polynomial.
24  * Preconditioning is not specified yet.
25  * Can supply a default dsdq and pseudoferm refresh algorithm
26  */
27  template<typename P, typename Q, typename Phi>
29  {
30  public:
31  //! virtual destructor:
33 
34  //! Compute the total action
35  virtual Double S(const AbsFieldState<P,Q>& s) = 0;
36 
37  //! Compute dsdq for the system...
38  /*! Actions of the form chi^dag*f_2(A_2)*f_1(A_1)*f_2(A_2)*chi */
39  virtual void dsdq(P& F, const AbsFieldState<P,Q>& s)
40  {
41  START_CODE();
42 
43  // Self Description/Encapsulation Rule
44  XMLWriter& xml_out = TheXMLLogWriter::Instance();
45  push(xml_out, "OneFlavorRatioRatConvExactWilsonTypeFermMonomial");
46 
47  /**** Identical code for unprec and even-odd prec case *****/
48 
49  // S_f = chi^dag*M_2*f_1(A_1)*M_2^dag*chi
50  // Here, A_1=M_1^dag*M_1 is some 4D operator and M_2
51  // is the preconditioning linop. The function
52  // f(A) = norm + \sum_i p_i/(A + q_i)
53  //
54  // We apply f(A)*chi and define the multi-shift solutions
55  // (A+q_i)\hat{chi}_i = chi
56  // (A+q_i)\hat{phi}_i = phi
57  //
58  // A core piece of math is that
59  // S_sample = chi^dag*f(A)*phi
60  //
61  // dS_sample/dU
62  // = chi^dag*df(A)*phi
63  // = - \sum_i p_i * chi^dag* (A+q_i)^{-1} * dA * (A+q_i)^{-1}*phi
64  // = - \sum_i p_i * \hat{chi}_i^dag * dA * \hat{phi}_i
65  // so, two solutions are needed
66  // Also, dA = dM^dag*M + M^dag*dM
67  //
68  // Need
69  // dS_f/dU = chi^dag*df_2(A_2)*f_1(A_1)*f_2(A_2)*chi
70  // + chi^dag*f_2(A_2)*f_1(A_1)*df_2(A_2)*chi
71  // + chi^dag*f_2(A_2)*df_1(A_1)*f_2(A_2)*chi
72  //
73  // Fermion action
75 
76  // Create a state for linop
78 
79  // Need way to get gauge state from AbsFieldState<P,Q>
81 
82  // Get multi-shift system solver
84 
85  // Fermion action
87 
88  // Need way to get gauge state from AbsFieldState<P,Q>
90 
91  // Partial Fraction Expansion coeffs for force
92  const RemezCoeff_t& fpfe_num = getNumerFPFE();
93 
94  multi1d<Phi> X;
95  Phi Y;
96 
97  P F_1, F_2, F_tmp(Nd);
98  F.resize(Nd);
99  F = zero;
100 
101  // Loop over all the pseudoferms
102  multi1d<int> n_count(getNPF());
103 
104  for(int n=0; n < getNPF(); ++n)
105  {
106  // M_dag_den phi = M^{dag}_den \phi - the RHS
107  Phi M_dag_den_phi;
108  (*M_den)(M_dag_den_phi, getPhi()[n], MINUS);
109 
110  // The multi-shift inversion
111  SystemSolverResults_t res = (*invMdagM_num)(X, fpfe_num.pole, M_dag_den_phi);
112  n_count[n] = res.n_count;
113 
114  // Loop over solns and accumulate force contributions
115  F_tmp = zero;
116  for(int i=0; i < X.size(); ++i)
117  {
118  (*M_num)(Y, X[i], PLUS);
119 
120  // The d(M^dag)*M term
121  M_num->deriv(F_1, X[i], Y, MINUS);
122 
123  // The M^dag*d(M) term
124  M_num->deriv(F_2, Y, X[i], PLUS);
125  F_1 += F_2;
126 
127  // Reweight each contribution in partial fraction
128  for(int mu=0; mu < F.size(); mu++)
129  F_tmp[mu] -= fpfe_num.res[i] * F_1[mu];
130  }
131 
132  // Concious choice. Don't monitor forces by pole
133  // Just accumulate it
134  F += F_tmp;
135  }
136 
137  state->deriv(F);
138  write(xml_out, "n_count", n_count);
139  monitorForces(xml_out, "Forces", F);
140  pop(xml_out);
141 
142  END_CODE();
143  }
144 
145 
146  //! Refresh pseudofermions
147  /*!
148  * This routine calculates the pseudofermion field (chi) for the case
149  * of rational evolution
150  *
151  * chi = n(Q)*[d(Q)]^(-1) * eta
152  * Where: Q = M^dag*M
153  * d(Q) = (Q+q_1)*(Q+q_2)*...*(Q+q_m)
154  * n(Q) = (Q+p_1)*(Q+p_2)*... *(Q+p_m)
155  * m = HBRatDeg
156  *
157  * The rational function n(x)/d(x) is the optimal rational
158  * approximation to the inverse square root of N(x)/D(x) which in
159  * turn is the optimal rational approximation to x^(-alpha).
160  * Here, alpha = 1/2
161  *
162  * To solve {n(Q)*[d(Q)]^(-1) * eta} the partial fraction expansion is
163  * used in combination with a multishift solver.
164  */
166  {
167  START_CODE();
168 
169  // Heatbath all the fields
170 
171  // Self Description/Encapsulation Rule
172  XMLWriter& xml_out = TheXMLLogWriter::Instance();
173  push(xml_out, "OneFlavorRatioRatConvExactWilsonTypeFermMonomial");
174 
175  // Get at the fermion action for piece i
177 
178  // Create a Connect State, apply fermionic boundaries
179  Handle< FermState<Phi,P,Q> > state(FA_num.createState(s.getQ()));
180 
181  // Create a linear operator
183 
184  // Get multi-shift system solver
186 
187  // Get at the fermion action for piece i
189 
190  // Create a linear operator
192 
193  // Partial Fraction Expansion coeffs for heat-bath
194  const RemezCoeff_t& sipfe = getNumerSIPFE();
195 
196  // Loop over pseudoferms
197  getPhi().resize(getNPF());
198  multi1d<int> n_count(getNPF());
199  Phi eta;
200 
201  for(int n=0; n < getNPF(); ++n)
202  {
203  // Fill the eta field with gaussian noise
204  eta = zero;
205  gaussian(eta, M_num->subset());
206 
207  // Account for fermion BC by modifying the proposed field
208  FA_num.getFermBC().modifyF(eta);
209 
210  // Temporary: Move to correct normalisation
211  eta *= sqrt(0.5);
212 
213  // The multi-shift inversion
214  multi1d<Phi> X;
215  SystemSolverResults_t res = (*invMdagM_num)(X, sipfe.pole, eta);
216  n_count[n] = res.n_count;
217 
218  // Weight solns to make final PF field
219  getPhi()[n][M_num->subset()] = sipfe.norm * eta;
220  for(int i=0; i < X.size(); ++i)
221  getPhi()[n][M_num->subset()] += sipfe.res[i] * X[i];
222  }
223 
224  write(xml_out, "n_count", n_count);
225  pop(xml_out);
226 
227  END_CODE();
228  }
229 
230  //! Copy pseudofermions if any
231  virtual void setInternalFields(const Monomial<P,Q>& m)
232  {
233  START_CODE();
234 
235  try {
237 
238  getPhi() = fm.getPhi();
239  }
240  catch(std::bad_cast) {
241  QDPIO::cerr << "Failed to cast input Monomial to OneFlavorRatioRatConvExactWilsonTypeFermMonomial " << std::endl;
242  QDP_abort(1);
243  }
244 
245  END_CODE();
246  }
247 
248 
249  //! Compute the action on the appropriate subset
250  /*!
251  * This measures the pseudofermion contribution to the Hamiltonian
252  * for the case of rational evolution (with polynomials n(x) and d(x),
253  * of degree SRatDeg
254  *
255  * S_f = chi_dag * (n(A)*d(A)^(-1))^2* chi
256  *
257  * where A is M^dag*M
258  *
259  * The rational function n(x)/d(x) is the optimal rational
260  * approximation to the square root of N(x)/D(x) which in
261  * turn is the optimal rational approximation to x^(-alpha).
262  * Here, alpha = 1/2
263  */
264  virtual Double S_subset(const AbsFieldState<P,Q>& s) const
265  {
266  START_CODE();
267 
268  XMLWriter& xml_out = TheXMLLogWriter::Instance();
269  push(xml_out, "S_subset");
270 
272 
273  Handle< FermState<Phi,P,Q> > state = FA_num.createState(s.getQ());
274 
275  // Need way to get gauge state from AbsFieldState<P,Q>
277 
278  // Get multi-shift system solver
280 
281  // Fermion action
283 
284  // Need way to get gauge state from AbsFieldState<P,Q>
286 
287  // Partial Fraction Expansion coeffs for action
288  const RemezCoeff_t& spfe = getNumerSPFE();
289 
290  // Compute energy
291  // Get X out here via multisolver
292  multi1d<Phi> X;
293 
294  // Loop over all the pseudoferms
295  multi1d<int> n_count(getNPF());
296  Double action = zero;
297  Phi psi;
298 
299  for(int n=0; n < getNPF(); ++n)
300  {
301  // The multi-shift inversion
302  SystemSolverResults_t res = (*invMdagM_num)(X, spfe.pole, getPhi()[n]);
303  n_count[n] = res.n_count;
304 
305  // Weight solns to make final PF field
306  psi[M_num->subset()] = spfe.norm * getPhi()[n];
307  for(int i=0; i < X.size(); ++i)
308  psi[M_num->subset()] += spfe.res[i] * X[i];
309 
310  action += norm2(psi, M_num->subset());
311  }
312 
313  write(xml_out, "n_count", n_count);
314  write(xml_out, "S", action);
315  pop(xml_out);
316 
317  END_CODE();
318 
319  return action;
320  }
321 
322 
323  protected:
324  //! Get at fermion action
325  virtual const WilsonTypeFermAct<Phi,P,Q>& getFermAct(void) const
326  {return getNumerFermAct();}
327 
328  //! Get at fermion action
329  virtual const WilsonTypeFermAct<Phi,P,Q>& getNumerFermAct(void) const = 0;
330 
331  //! Get at fermion action
332  virtual const WilsonTypeFermAct<Phi,P,Q>& getDenomFermAct(void) const = 0;
333 
334  //! Get inverter params
335  virtual const GroupXML_t& getNumerActionInvParams(void) const = 0;
336 
337  //! Get inverter params
338  virtual const GroupXML_t& getNumerForceInvParams(void) const = 0;
339 
340  //! Return the partial fraction expansion for the force calc
341  virtual const RemezCoeff_t& getNumerFPFE() const = 0;
342 
343  //! Return the partial fraction expansion for the action calc
344  virtual const RemezCoeff_t& getNumerSPFE() const = 0;
345 
346  //! Return the partial fraction expansion for the heat-bath
347  virtual const RemezCoeff_t& getNumerSIPFE() const = 0;
348 
349  //! Return number of roots in used
350  virtual int getNPF() const = 0;
351 
352  //! Accessor for pseudofermion (read only)
353  virtual const multi1d<Phi>& getPhi(void) const = 0;
354 
355  //! mutator for pseudofermion
356  virtual multi1d<Phi>& getPhi(void) = 0;
357 
358  };
359 
360 
361  //-------------------------------------------------------------------------------------------
362  //! Exact 1 flavor unpreconditioned fermact monomial
363  /*! @ingroup monomial
364  *
365  * Exact 1 flavor unpreconditioned fermact monomial.
366  */
367  template<typename P, typename Q, typename Phi>
369  {
370  public:
371  //! virtual destructor:
373 
374  //! Compute the total action
375  virtual Double S(const AbsFieldState<P,Q>& s)
376  {
377  START_CODE();
378 
379  // Self identification/encapsulation Rule
380  XMLWriter& xml_out = TheXMLLogWriter::Instance();
381  push(xml_out, "OneFlavorRatioRatConvExactUnprecWilsonTypeFermMonomial");
382 
383  Double action = this->S_subset(s);
384 
385  write(xml_out, "S", action);
386  pop(xml_out);
387 
388  END_CODE();
389 
390  return action;
391  }
392 
393  protected:
394  //! Get at fermion action
395  virtual const WilsonTypeFermAct<Phi,P,Q>& getNumerFermAct(void) const = 0;
396 
397  //! Get at fermion action
398  virtual const WilsonTypeFermAct<Phi,P,Q>& getDenomFermAct(void) const = 0;
399  };
400 
401 
402  //-------------------------------------------------------------------------------------------
403  //! Exact 1 flavor even-odd preconditioned fermact monomial
404  /*! @ingroup monomial
405  *
406  * Exact 1 flavor even-odd preconditioned fermact monomial.
407  * Can supply a default dsdq algorithm
408  */
409  template<typename P, typename Q, typename Phi>
411  {
412  public:
413  //! virtual destructor:
415 
416  //! Even even contribution (eg ln det Clover)
417  virtual Double S_even_even(const AbsFieldState<P,Q>& s) = 0;
418 
419  //! Compute the odd odd contribution (eg
421  {
422  return this->S_subset(s);
423  }
424 
425  //! Compute the total action
427  {
428  START_CODE();
429 
430  XMLWriter& xml_out = TheXMLLogWriter::Instance();
431  push(xml_out, "OneFlavorRatioRatConvExactEvenOddPrecWilsonTypeFermMonomial");
432 
433  Double action_e = S_even_even(s);
434  Double action_o = S_odd_odd(s);
435  Double action = action_e + action_o;
436 
437  write(xml_out, "S_even_even", action_e);
438  write(xml_out, "S_odd_odd", action_o);
439  write(xml_out, "S", action);
440  pop(xml_out);
441 
442  END_CODE();
443 
444  return action;
445  }
446 
447  protected:
448  //! Get at fermion action
450 
451  //! Get at fermion action
453  };
454 
455  //-------------------------------------------------------------------------------------------
456  //! Exact 1 flavor even-odd preconditioned fermact monomial constant determinant
457  // Can fill out the S_odd_odd piece
458 
459  /*! @ingroup monomial
460  *
461  * Exact 1 flavor even-odd preconditioned fermact monomial.
462  * Can supply a default dsdq algorithm
463  */
464  template<typename P, typename Q, typename Phi>
467  {
468  public:
469  //! virtual destructor:
471 
472  //! Even even contribution (eg ln det Clover)
474  return Double(0);
475  }
476 
477  protected:
478  //! Get at fermion action
480 
481  //! Get at fermion action
483  };
484 
485 }
486 
487 
488 #endif
Monomials - gauge action or fermion binlinear contributions for HMC.
Abstract field state.
Definition: field_state.h:27
virtual DiffLinearOperator< T, Q, P > * linOp(Handle< FermState< T, P, Q > > state) const =0
Produce a linear operator for this action.
Even-odd preconditioned Wilson-like fermion actions including derivatives.
Fermionic monomials (binlinears in fermion fields)
Definition: abs_monomial.h:243
virtual const FermBC< T, P, Q > & getFermBC() const
Return the fermion BC object for this action.
Definition: fermact.h:79
virtual FermState< T, P, Q > * createState(const Q &q) const
Given links (coordinates Q) create the state needed for the linear operators.
Definition: fermact.h:59
Class for counted reference semantics.
Definition: handle.h:33
An abstract monomial class, for inexact algorithms.
Definition: abs_monomial.h:43
Exact 1 flavor even-odd preconditioned fermact monomial constant determinant.
virtual const EvenOddPrecWilsonTypeFermAct< Phi, P, Q > & getNumerFermAct() const =0
Get at fermion action.
virtual Double S_even_even(const AbsFieldState< P, Q > &s)
Even even contribution (eg ln det Clover)
virtual const EvenOddPrecWilsonTypeFermAct< Phi, P, Q > & getDenomFermAct() const =0
Get at fermion action.
virtual const EvenOddPrecWilsonTypeFermAct< Phi, P, Q > & getDenomFermAct() const =0
Get at fermion action.
virtual const EvenOddPrecWilsonTypeFermAct< Phi, P, Q > & getNumerFermAct() const =0
Get at fermion action.
Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
virtual Double S_odd_odd(const AbsFieldState< P, Q > &s)
Compute the odd odd contribution (eg.
virtual Double S_even_even(const AbsFieldState< P, Q > &s)=0
Even even contribution (eg ln det Clover)
virtual const WilsonTypeFermAct< Phi, P, Q > & getDenomFermAct(void) const =0
Get at fermion action.
virtual Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
virtual const WilsonTypeFermAct< Phi, P, Q > & getNumerFermAct(void) const =0
Get at fermion action.
Exact 1 flavor fermact monomial using rational polynomials.
virtual Double S_subset(const AbsFieldState< P, Q > &s) const
Compute the action on the appropriate subset.
virtual int getNPF() const =0
Return number of roots in used.
virtual const WilsonTypeFermAct< Phi, P, Q > & getDenomFermAct(void) const =0
Get at fermion action.
virtual Double S(const AbsFieldState< P, Q > &s)=0
Compute the total action.
virtual void dsdq(P &F, const AbsFieldState< P, Q > &s)
Compute dsdq for the system...
virtual const GroupXML_t & getNumerForceInvParams(void) const =0
Get inverter params.
virtual multi1d< Phi > & getPhi(void)=0
mutator for pseudofermion
virtual const RemezCoeff_t & getNumerSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
virtual void setInternalFields(const Monomial< P, Q > &m)
Copy pseudofermions if any.
virtual const GroupXML_t & getNumerActionInvParams(void) const =0
Get inverter params.
virtual void refreshInternalFields(const AbsFieldState< P, Q > &s)
Refresh pseudofermions.
virtual const WilsonTypeFermAct< Phi, P, Q > & getNumerFermAct(void) const =0
Get at fermion action.
virtual const multi1d< Phi > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual const WilsonTypeFermAct< Phi, P, Q > & getFermAct(void) const
Get at fermion action.
virtual const RemezCoeff_t & getNumerFPFE() const =0
Return the partial fraction expansion for the force calc.
virtual const RemezCoeff_t & getNumerSPFE() const =0
Return the partial fraction expansion for the action calc.
static T & Instance()
Definition: singleton.h:432
Wilson-like fermion actions.
Definition: fermact.orig.h:344
virtual MdagMMultiSystemSolver< T > * mInvMdagM(Handle< FermState< T, P, Q > > state, const GroupXML_t &invParam) const
Return a multi-shift linear operator solver for this action to solve (MdagM+shift)*psi=chi.
int mu
Definition: cool.cc:24
Even-odd const determinant Wilson-like fermact.
Helper function for calculating forces.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void monitorForces(XMLWriter &xml_out, const std::string &path, const multi1d< LatticeColorMatrix > &F)
Calculate and write out forces.
unsigned n
Definition: ldumul_w.cc:36
static int m[4]
Definition: make_seeds.cc:16
Nd
Definition: meslate.cc:74
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
int n_count
Definition: invbicg.cc:78
push(xml_out,"Condensates")
int i
Definition: pbg5p_w.cc:55
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
LatticeFermion eta
Definition: mespbg5p_w.cc:37
START_CODE()
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
Double zero
Definition: invbicg.cc:106
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
Definition: gtest.h:7351
Remez algorithm coefficients.
Hold group xml and type id.
Convenient structure to package Remez coeffs.
Definition: remez_coeff.h:19
multi1d< Real > pole
Definition: remez_coeff.h:21
multi1d< Real > res
Definition: remez_coeff.h:20
Holds return info from SystemSolver call.
Definition: syssolver.h:17
multi1d< LatticeColorMatrix > P
Definition: t_clover.cc:13
Wilson-like fermion actions.
static INTERNAL_PRECISION F