CHROMA
one_flavor_ratio_rat_conv_monomial5d_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_monomial5d_w_h__
8 #define __one_flavor_ratio_rat_conv_monomial5d_w_h__
9 
15 
16 #include <typeinfo>
17 
18 namespace Chroma
19 {
20  //-------------------------------------------------------------------------------------------
21  //! Exact 1 flavor fermact monomial in extra dimensions
22  /*! @ingroup monomial
23  *
24  * Exact 1 flavor fermact monomial. Preconditioning is not
25  * specified yet.
26  * Can supply a default dsdq and pseudoferm refresh algorithm
27  */
28  template<typename P, typename Q, typename Phi>
30  {
31  public:
32  //! virtual destructor:
34 
35  //! Compute the total action
36  virtual Double S(const AbsFieldState<P,Q>& s) = 0;
37 
38  //! Compute dsdq for the system...
39  /*! Monomial of the form chi^dag*(M^dag*M)*chi */
40  virtual void dsdq(P& F, const AbsFieldState<P,Q>& s)
41  {
42  START_CODE();
43 
44  // SelfIdentification/Encapsultaion Rule
45  XMLWriter& xml_out = TheXMLLogWriter::Instance();
46  push(xml_out, "OneFlavorRatioRatConvExactWilsonTypeFermMonomial5D");
47 
48  /**** Identical code for unprec and even-odd prec case *****/
49 
50  // S_f = chi^dag* P(V^dag*V)/Q(V^dag*V)* N(M^dag*M)/D(M^dag*M)* P(V^dag*V)/Q(V^dag*V)* chi
51  //
52  // Here, M is some 5D operator and V is the Pauli-Villars field
53  // N and D makeup the rat. poly of the M term and P and & makeup the rat.poly of the den term
54  //
55  // Need
56  // dS_f/dU = - \sum_i psi_i^dag * n_i * [d(M^dag)*M + M^dag*dM] * psi
57  //
58  // where psi_i = (M^dag*M + n_i)^(-1)*chi
59  //
60  // In Balint's notation, the result is
61  // \dot{S} = -\sum_i p_i [ X_i^dag*\dot{M}^\dag*Y_i - Y_i^dag\dot{M}*X_i]
62  // where
63  // X_i = (M^dag*M + q_i)^(-1)*chi Y_i = M*X_i
64  // In Robert's notation, X_i -> psi_i .
65  //
67 
68  // Create a state for linop
70 
71  // Start the force
72  F.resize(Nd);
73  F = zero;
74 
75  // Force term for the linop
76  // Get multi-shift system solver
78 
79  // Get linear operator
81 
82  // Partial Fraction Expansion coeffs for force
83  const RemezCoeff_t& fpfe = getNumerFPFE();
84 
85  multi1d< multi1d<Phi> > X;
86  multi1d<Phi> Y;
87  P F_1, F_2, F_tmp(Nd);
88  multi1d<int> n_m_count(getNPF());
89 
90  // Loop over the pseudoferms
91  for(int n=0; n < getNPF(); ++n)
92  {
93  // The multi-shift inversion
94  SystemSolverResults_t res = (*invMdagM)(X, fpfe.pole, getPhi()[n]);
95  n_m_count[n] = res.n_count;
96 
97  // Loop over solns and accumulate force contributions
98  F_tmp = zero;
99  for(int i=0; i < X.size(); ++i)
100  {
101  (*M)(Y, X[i], PLUS);
102 
103  // The d(M^dag)*M term
104  M->deriv(F_1, X[i], Y, MINUS);
105 
106  // The M^dag*d(M) term
107  M->deriv(F_2, Y, X[i], PLUS);
108  F_1 += F_2;
109 
110  // Reweight each contribution in partial fraction
111  for(int mu=0; mu < F.size(); mu++)
112  F_tmp[mu] -= fpfe.res[i] * F_1[mu];
113  }
114 
115  // We are not monitoring by pole anymore, so I will
116  // Just accumulate this and take one derivative at the end
117  F += F_tmp; // add on base force terms
118  }
119 
120  // Take the derivative with respect to possibly fat links
121  state->deriv(F);
122 
123  // Write out the n_count array and the Forces for the Base Operator
124  write(xml_out, "n_m_count", n_m_count);
125  monitorForces(xml_out, "ForcesOperator", F);
126 
127  monitorForces(xml_out, "Forces", F);
128 
129  pop(xml_out);
130 
131  END_CODE();
132  }
133 
134 
135  //! Refresh pseudofermions
136  /*!
137  * This routine calculates the pseudofermion field (chi) for the case
138  * of rational evolution
139  *
140  * chi = n(Q)*[d(Q)]^(-1) * eta
141  * Where: Q = M^dag*M
142  * d(Q) = (Q+q_1)*(Q+q_2)*...*(Q+q_m)
143  * n(Q) = (Q+p_1)*(Q+p_2)*... *(Q+p_m)
144  * m = SIRatDeg
145  *
146  * The rational function n(x)/d(x) is the optimal rational
147  * approximation to the inverse square root of N(x)/D(x) which in
148  * turn is the optimal rational approximation to x^(-alpha).
149  * Here, alpha = 1/2
150  *
151  * To solve {n(Q)*[d(Q)]^(-1) * eta} the partial fraction expansion is
152  * used in combination with a multishift solver.
153  */
155  {
156  START_CODE();
157 
158  // SelfIdentification/Encapsultaion Rule
159  XMLWriter& xml_out = TheXMLLogWriter::Instance();
160  push(xml_out, "OneFlavorRatioRatConvExactWilsonTypeFermMonomial5DRefresh");
161 
162  // Heatbath all the fields
163 
164  // Get at the ferion action for piece i
166 
167  // Create a Connect State, apply fermionic boundaries
168  Handle< FermState<Phi,P,Q> > f_state(FA.createState(s.getQ()));
169 
170  // Force terms
171  const int N5 = FA.size();
172 
173  // Pseudofermions for M term
174  multi1d<int> n_m_count(getNPF());
175  getPhi().resize(getNPF()); // Will hold nth-root pseudoferms
176 
177  // Get multi-shift system solver
179 
180  // Get linear operator
182 
183  // Partial Fraction Expansion coeffs for heat-bath
184  const RemezCoeff_t& sipfe = getNumerSIPFE();
185 
186  multi1d<Phi> eta(N5);
187 
188  // Loop over the pseudoferms
189  for(int n=0; n < getNPF(); ++n)
190  {
191  // Fill the eta field with gaussian noise
192  eta = zero;
193  for(int i=0; i < N5; ++i)
194  gaussian(eta[i], M->subset());
195 
196  // Account for fermion BC by modifying the proposed field
197  FA.getFermBC().modifyF(eta);
198 
199  // Temporary: Move to correct normalisation
200  for(int i=0; i < N5; ++i)
201  eta[i][M->subset()] *= sqrt(0.5);
202 
203  // The multi-shift inversion
204  multi1d< multi1d<Phi> > X;
205  SystemSolverResults_t res = (*invMdagM)(X, sipfe.pole, eta);
206  n_m_count[n] = res.n_count;
207 
208  // Sanity checks
209  if (X.size() != sipfe.pole.size())
210  QDP_error_exit("%s : sanity failure, internal size not getSIPartFracRoot size", __func__);
211 
212  if (X[0].size() != N5)
213  QDP_error_exit("%s : sanity failure, internal size not N5", __func__);
214 
215  // Weight solns to make final PF field
216  getPhi()[n].resize(N5);
217  getPhi()[n] = zero;
218  // Loop over each 5d slice
219  for(int j=0; j < N5; ++j)
220  {
221  getPhi()[n][j][M->subset()] = sipfe.norm * eta[j];
222  for(int i=0; i < X.size(); ++i)
223  getPhi()[n][j][M->subset()] += sipfe.res[i] * X[i][j];
224  }
225  }
226 
227 
228  write(xml_out, "n_m_count", n_m_count);
229  pop(xml_out);
230 
231  END_CODE();
232  }
233 
234 
235  //! Copy internal fields
236  virtual void setInternalFields(const Monomial<P,Q>& m)
237  {
238  START_CODE();
239 
240  try
241  {
243 
244  // Do a resize here -- otherwise if the fields have not yet
245  // been refreshed there may be trouble
246  getPhi().resize(fm.getPhi().size());
247  for(int i=0 ; i < fm.getPhi().size(); i++) {
248  (getPhi())[i] = (fm.getPhi())[i];
249  }
250  }
251  catch(std::bad_cast) {
252  QDPIO::cerr << "Failed to cast input Monomial to OneFlavorRatioRatConvExactWilsonTypeFermMonomial5D" << std::endl;
253  QDP_abort(1);
254  }
255 
256  END_CODE();
257  }
258 
259 
260  //! Compute action on the appropriate subset
261  /*!
262  * This may only be a Piece Of The Action
263  *
264  * This function measures the pseudofermion contribution to the Hamiltonian
265  * for the case of rational evolution (with polynomials n(x) and d(x),
266  * of degree SRatDeg
267  *
268  * S_f = chi_dag * (n(A)*d(A)^(-1))^2* chi
269  *
270  * where A is M^dag*M
271  *
272  * The rational function n(x)/d(x) is the optimal rational
273  * approximation to the square root of N(x)/D(x) which in
274  * turn is the optimal rational approximation to x^(-alpha).
275  * Here, alpha = 1/2
276  */
277  virtual Double S_subset(const AbsFieldState<P,Q>& s) const
278  {
279  START_CODE();
280 
281  XMLWriter& xml_out = TheXMLLogWriter::Instance();
282  push(xml_out, "S_subset");
283 
285 
286  // Create a Connect State, apply fermionic boundaries
287  Handle< FermState<Phi,P,Q> > bc_g_state(FA.createState(s.getQ()));
288 
289  // Force terms
290  const int N5 = FA.size();
291 
292  // Action for M term
293  Double action_m = zero;
294  multi1d<int> n_m_count(getNPF());
295 
296  // Get multi-shift system solver
298 
299  // Get linear operator
301 
302  // Partial Fraction Expansion coeffs for action
303  const RemezCoeff_t& spfe = getNumerSPFE();
304 
305  // Get X out here via multisolver
306  multi1d< multi1d<Phi> > X;
307  multi1d<Phi> tmp(N5);
308 
309  // Loop over the pseudoferms
310  for(int n=0; n < getNPF(); ++n)
311  {
312  // The multi-shift inversion
313  SystemSolverResults_t res = (*invMdagM)(X, spfe.pole, getPhi()[n]);
314  n_m_count[n] = res.n_count;
315 
316  // Sanity checks
317  if (X.size() != spfe.pole.size())
318  QDP_error_exit("%s : sanity failure, internal size not getSPartFracRoot size", __func__);
319 
320  if (X[0].size() != N5)
321  QDP_error_exit("%s : sanity failure, internal size not N5", __func__);
322 
323  // Weight solns
324  // Loop over each 5d slice
325  for(int j=0; j < N5; ++j)
326  {
327  tmp[j][M->subset()] = spfe.norm * getPhi()[n][j];
328  for(int i=0; i < X.size(); ++i)
329  tmp[j][M->subset()] += spfe.res[i] * X[i][j];
330  }
331 
332  // Action on the subset
333  action_m += norm2(tmp, M->subset());
334  }
335 
336 
337  write(xml_out, "n_m_count", n_m_count);
338  write(xml_out, "S_m", action_m);
339  Double action = action_m;
340  write(xml_out, "S", action);
341  pop(xml_out);
342 
343  END_CODE();
344 
345  return action;
346  }
347 
348 
349  protected:
350  //! Get at fermion action
352  {return getNumerFermAct();}
353 
354  //! Get at fermion action
355  virtual const WilsonTypeFermAct5D<Phi,P,Q>& getNumerFermAct() const = 0;
356 
357  //! Get at fermion action
358  virtual const WilsonTypeFermAct5D<Phi,P,Q>& getDenomFermAct() const = 0;
359 
360  //! Get inverter params
361  virtual const GroupXML_t& getNumerActionInvParams(void) const = 0;
362 
363  //! Get inverter params
364  virtual const GroupXML_t& getNumerForceInvParams(void) const = 0;
365 
366  //! Return the partial fraction expansion for the force calc
367  virtual const RemezCoeff_t& getNumerFPFE() const = 0;
368 
369  //! Return the partial fraction expansion for the action calc
370  virtual const RemezCoeff_t& getNumerSPFE() const = 0;
371 
372  //! Return the partial fraction expansion for the heat-bath
373  virtual const RemezCoeff_t& getNumerSIPFE() const = 0;
374 
375  //! Accessor for pseudofermion (read only)
376  virtual const multi1d< multi1d<Phi> >& getPhi(void) const = 0;
377 
378  //! mutator for pseudofermion
379  virtual multi1d< multi1d<Phi> >& getPhi(void) = 0;
380 
381  //! Return number of pseudofermions
382  virtual int getNPF() const = 0;
383  };
384 
385 
386  //-------------------------------------------------------------------------------------------
387  //! Exact 1 flavor unpreconditioned fermact monomial living in extra dimensions
388  /*! @ingroup monomial
389  *
390  * Exact 1 flavor unpreconditioned fermact monomial.
391  */
392  template<typename P, typename Q, typename Phi>
394  {
395  public:
396  //! virtual destructor:
398 
399  //! Compute the total action
400  virtual Double S(const AbsFieldState<P,Q>& s)
401  {
402  START_CODE();
403 
404  XMLWriter& xml_out=TheXMLLogWriter::Instance();
405  push(xml_out, "OneFlavorRatioRatConvExactUnprecWilsonTypeFermMonomial5D");
406 
407  Double action = this->S_subset(s);
408 
409  write(xml_out, "S", action);
410  pop(xml_out);
411 
412  END_CODE();
413 
414  return action;
415  }
416 
417 
418  protected:
419  //! Get at fermion action
421 
422  //! Get at fermion action
424  };
425 
426 
427  //-------------------------------------------------------------------------------------------
428  //! Exact 1 flavor even-odd preconditioned fermact monomial living in extra dimensions
429  /*! @ingroup monomial
430  *
431  * Exact 1 flavor even-odd preconditioned fermact monomial.
432  * Can supply a default dsdq algorithm
433  */
434  template<typename P, typename Q, typename Phi>
436  {
437  public:
438  //! virtual destructor:
440 
441  //! Even even contribution (eg ln det Clover)
442  virtual Double S_even_even(const AbsFieldState<P,Q>& s) = 0;
443 
444  //! Compute the odd odd contribution (eg
446  {
447  return this->S_subset(s);
448  }
449 
450  //! Compute the total action
452  {
453  START_CODE();
454 
455  XMLWriter& xml_out=TheXMLLogWriter::Instance();
456  push(xml_out, "OneFlavorRatioRatConvExactEvenOddPrecWilsonTypeFermMonomial5D");
457 
458  Double action_e = S_even_even(s);
459  Double action_o = S_odd_odd(s);
460  Double action = action_e + action_o;
461 
462  write(xml_out, "S_even_even", action_e);
463  write(xml_out, "S_odd_odd", action_o);
464  write(xml_out, "S", action);
465  pop(xml_out);
466 
467  END_CODE();
468 
469  return action;
470  }
471 
472  protected:
473  //! Get at fermion action
475 
476  //! Get at fermion action
478  };
479 
480  //-------------------------------------------------------------------------------------------
481  //! Exact 1 flavor even-odd preconditioned fermact monomial living in extra dimensions
482  /*! @ingroup monomial
483  *
484  * Exact 1 flavor even-odd preconditioned fermact monomial.
485  * Can supply a default dsdq algorithm
486  */
487  template<typename P, typename Q, typename Phi>
490  {
491  public:
492  //! virtual destructor:
494 
495  //! Even even contribution (eg ln det Clover)
497  return Double(0);
498  }
499 
500 
501  protected:
502  //! Get at fermion action
504 
505  //! Get at fermion action
507  };
508 
509 }
510 
511 
512 #endif
Monomials - gauge action or fermion binlinear contributions for HMC.
Abstract field state.
Definition: field_state.h:27
virtual DiffLinearOperatorArray< T, P, Q > * 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:280
virtual int size() const =0
Expected length of array index.
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 living in extra dimensions.
virtual Double S_even_even(const AbsFieldState< P, Q > &s)
Even even contribution (eg ln det Clover)
virtual const EvenOddPrecWilsonTypeFermAct5D< Phi, P, Q > & getNumerFermAct() const =0
Get at fermion action.
virtual const EvenOddPrecWilsonTypeFermAct5D< Phi, P, Q > & getDenomFermAct() const =0
Get at fermion action.
Exact 1 flavor even-odd preconditioned fermact monomial living in extra dimensions.
virtual Double S_odd_odd(const AbsFieldState< P, Q > &s)
Compute the odd odd contribution (eg.
virtual const EvenOddPrecWilsonTypeFermAct5D< Phi, P, Q > & getNumerFermAct() const =0
Get at fermion action.
virtual const EvenOddPrecWilsonTypeFermAct5D< Phi, P, Q > & getDenomFermAct() const =0
Get at fermion action.
virtual Double S_even_even(const AbsFieldState< P, Q > &s)=0
Even even contribution (eg ln det Clover)
Exact 1 flavor unpreconditioned fermact monomial living in extra dimensions.
virtual Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
virtual const UnprecWilsonTypeFermAct5D< Phi, P, Q > & getNumerFermAct() const =0
Get at fermion action.
virtual const UnprecWilsonTypeFermAct5D< Phi, P, Q > & getDenomFermAct() const =0
Get at fermion action.
virtual const RemezCoeff_t & getNumerSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
virtual Double S_subset(const AbsFieldState< P, Q > &s) const
Compute action on the appropriate subset.
virtual const WilsonTypeFermAct5D< Phi, P, Q > & getDenomFermAct() const =0
Get at fermion action.
virtual void setInternalFields(const Monomial< P, Q > &m)
Copy internal fields.
virtual void refreshInternalFields(const AbsFieldState< P, Q > &s)
Refresh pseudofermions.
virtual const multi1d< multi1d< Phi > > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual int getNPF() const =0
Return number of pseudofermions.
virtual const WilsonTypeFermAct5D< Phi, P, Q > & getNumerFermAct() const =0
Get at fermion action.
virtual const GroupXML_t & getNumerForceInvParams(void) const =0
Get inverter params.
virtual const RemezCoeff_t & getNumerSPFE() const =0
Return the partial fraction expansion for the action calc.
virtual const RemezCoeff_t & getNumerFPFE() const =0
Return the partial fraction expansion for the force calc.
virtual void dsdq(P &F, const AbsFieldState< P, Q > &s)
Compute dsdq for the system...
virtual const GroupXML_t & getNumerActionInvParams(void) const =0
Get inverter params.
virtual const WilsonTypeFermAct5D< Phi, P, Q > & getFermAct() const
Get at fermion action.
virtual multi1d< multi1d< Phi > > & getPhi(void)=0
mutator for pseudofermion
virtual Double S(const AbsFieldState< P, Q > &s)=0
Compute the total action.
static T & Instance()
Definition: singleton.h:432
Unpreconditioned Wilson-like fermion actions in extra dims with derivatives.
Definition: fermact.orig.h:571
Wilson-like fermion actions.
Definition: fermact.orig.h:403
virtual MdagMMultiSystemSolverArray< 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 j
Definition: ldumul_w.cc:35
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)
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
push(xml_out,"Condensates")
int i
Definition: pbg5p_w.cc:55
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
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