CHROMA
one_flavor_rat_monomial5d_w.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 
3 /*! @file
4  * @brief One flavor monomials using RHMC
5  */
6 
7 #ifndef __one_flavor_rat_monomial5d_w_h__
8 #define __one_flavor_rat_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, "OneFlavorRatExactWilsonTypeFermMonomial5D");
47 
48  /**** Identical code for unprec and even-odd prec case *****/
49 
50  // S_f = chi^dag*N(M^dag*M)/D(M^dag*M)*chi
51  //
52  // Here, M is some 5D operator
53  // N and D makeup the rat. poly of the M 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  {
77  // Get multi-shift system solver
79 
80  // Get linear operator
82 
83  // Partial Fraction Expansion coeffs for force
84  const RemezCoeff_t& fpfe = getFPFE();
85 
86  multi1d< multi1d<Phi> > X;
87  multi1d<Phi> Y;
88  P F_1, F_2, F_tmp(Nd);
89  multi1d<int> n_m_count(getNPF());
90 
91  // Loop over the pseudoferms
92  for(int n=0; n < getNPF(); ++n)
93  {
94  // The multi-shift inversion
95  SystemSolverResults_t res = (*invMdagM)(X, fpfe.pole, getPhi()[n]);
96  n_m_count[n] = res.n_count;
97 
98  // Loop over solns and accumulate force contributions
99  F_tmp = zero;
100  for(int i=0; i < X.size(); ++i)
101  {
102  (*M)(Y, X[i], PLUS);
103 
104  // The d(M^dag)*M term
105  M->deriv(F_1, X[i], Y, MINUS);
106 
107  // The M^dag*d(M) term
108  M->deriv(F_2, Y, X[i], PLUS);
109  F_1 += F_2;
110 
111  // Reweight each contribution in partial fraction
112  for(int mu=0; mu < F.size(); mu++)
113  F_tmp[mu] -= fpfe.res[i] * F_1[mu];
114  }
115 
116  // We are not monitoring by pole anymore, so I will
117  // Just accumulate this and take one derivative at the end
118  F += F_tmp; // add on base force terms
119  }
120 
121  // Take the derivative with respect to possibly fat links
122  state->deriv(F);
123 
124  // Write out the n_count array and the Forces for the Base Operator
125  write(xml_out, "n_m_count", n_m_count);
126  monitorForces(xml_out, "ForcesOperator", F);
127  }
128 
129  monitorForces(xml_out, "Forces", F);
130 
131  pop(xml_out);
132 
133  END_CODE();
134  }
135 
136 
137  //! Refresh pseudofermions
138  /*!
139  * This routine calculates the pseudofermion field (chi) for the case
140  * of rational evolution
141  *
142  * chi = n(Q)*[d(Q)]^(-1) * eta
143  * Where: Q = M^dag*M
144  * d(Q) = (Q+q_1)*(Q+q_2)*...*(Q+q_m)
145  * n(Q) = (Q+p_1)*(Q+p_2)*... *(Q+p_m)
146  * m = SIRatDeg
147  *
148  * The rational function n(x)/d(x) is the optimal rational
149  * approximation to the inverse square root of N(x)/D(x) which in
150  * turn is the optimal rational approximation to x^(-alpha).
151  * Here, alpha = 1/2
152  *
153  * To solve {n(Q)*[d(Q)]^(-1) * eta} the partial fraction expansion is
154  * used in combination with a multishift solver.
155  */
157  {
158  START_CODE();
159 
160  // SelfIdentification/Encapsultaion Rule
161  XMLWriter& xml_out = TheXMLLogWriter::Instance();
162  push(xml_out, "OneFlavorRatExactWilsonTypeFermMonomial5DRefresh");
163 
164  // Heatbath all the fields
165 
166  // Get at the ferion action for piece i
168 
169  // Create a Connect State, apply fermionic boundaries
170  Handle< FermState<Phi,P,Q> > f_state(FA.createState(s.getQ()));
171 
172  // Force terms
173  const int N5 = FA.size();
174 
175  // Pseudofermions for M term
176  multi1d<int> n_m_count(getNPF());
177  getPhi().resize(getNPF()); // Will hold nth-root pseudoferms
178  {
179  // Get multi-shift system solver
181 
182  // Get linear operator
184 
185  // Partial Fraction Expansion coeffs for heat-bath
186  const RemezCoeff_t& sipfe = getSIPFE();
187 
188  multi1d<Phi> eta(N5);
189 
190  // Loop over the pseudoferms
191  for(int n=0; n < getNPF(); ++n)
192  {
193  // Fill the eta field with gaussian noise
194  eta = zero;
195  for(int i=0; i < N5; ++i)
196  gaussian(eta[i], M->subset());
197 
198  // Account for fermion BC by modifying the proposed field
199  FA.getFermBC().modifyF(eta);
200 
201  // Temporary: Move to correct normalisation
202  for(int i=0; i < N5; ++i)
203  eta[i][M->subset()] *= sqrt(0.5);
204 
205  // The multi-shift inversion
206  multi1d< multi1d<Phi> > X;
207  SystemSolverResults_t res = (*invMdagM)(X, sipfe.pole, eta);
208  n_m_count[n] = res.n_count;
209 
210  // Sanity checks
211  if (X.size() != sipfe.pole.size())
212  QDP_error_exit("%s : sanity failure, internal size not getSIPartFracRoot size", __func__);
213 
214  if (X[0].size() != N5)
215  QDP_error_exit("%s : sanity failure, internal size not N5", __func__);
216 
217  // Weight solns to make final PF field
218  getPhi()[n].resize(N5);
219  getPhi()[n] = zero;
220  // Loop over each 5d slice
221  for(int j=0; j < N5; ++j)
222  {
223  getPhi()[n][j][M->subset()] = sipfe.norm * eta[j];
224  for(int i=0; i < X.size(); ++i)
225  getPhi()[n][j][M->subset()] += sipfe.res[i] * X[i][j];
226  }
227  }
228  }
229 
230  write(xml_out, "n_m_count", n_m_count);
231  pop(xml_out);
232 
233  END_CODE();
234  }
235 
236 
237  //! Copy internal fields
238  virtual void setInternalFields(const Monomial<P,Q>& m)
239  {
240  START_CODE();
241 
242  try {
244 
245  // Do a resize here -- otherwise if the fields have not yet
246  // been refreshed there may be trouble
247  getPhi().resize(fm.getPhi().size());
248  for(int i=0 ; i < fm.getPhi().size(); i++) {
249  (getPhi())[i] = (fm.getPhi())[i];
250  }
251  }
252  catch(std::bad_cast) {
253  QDPIO::cerr << "Failed to cast input Monomial to OneFlavorRatExactWilsonTypeFermMonomial5D" << std::endl;
254  QDP_abort(1);
255  }
256 
257  END_CODE();
258  }
259 
260 
261  //! Compute action on the appropriate subset
262  /*!
263  * This may only be a Piece Of The Action
264  *
265  * This function measures the pseudofermion contribution to the Hamiltonian
266  * for the case of rational evolution (with polynomials n(x) and d(x),
267  * of degree SRatDeg
268  *
269  * S_f = chi_dag * (n(A)*d(A)^(-1))^2* chi
270  *
271  * where A is M^dag*M
272  *
273  * The rational function n(x)/d(x) is the optimal rational
274  * approximation to the square root of N(x)/D(x) which in
275  * turn is the optimal rational approximation to x^(-alpha).
276  * Here, alpha = 1/2
277  */
278  virtual Double S_subset(const AbsFieldState<P,Q>& s) const
279  {
280  START_CODE();
281 
282  XMLWriter& xml_out = TheXMLLogWriter::Instance();
283  push(xml_out, "S_subset");
284 
286 
287  // Create a Connect State, apply fermionic boundaries
288  Handle< FermState<Phi,P,Q> > bc_g_state(FA.createState(s.getQ()));
289 
290  // Force terms
291  const int N5 = FA.size();
292 
293  // Action for M term
294  Double action_m = zero;
295  multi1d<int> n_m_count(getNPF());
296  {
297  // Get multi-shift system solver
299 
300  // Get linear operator
302 
303  // Partial Fraction Expansion coeffs for action
304  const RemezCoeff_t& spfe = getSPFE();
305 
306  // Get X out here via multisolver
307  multi1d< multi1d<Phi> > X;
308  multi1d<Phi> tmp(N5);
309 
310  // Loop over the pseudoferms
311  for(int n=0; n < getNPF(); ++n)
312  {
313  // The multi-shift inversion
314  SystemSolverResults_t res = (*invMdagM)(X, spfe.pole, getPhi()[n]);
315  n_m_count[n] = res.n_count;
316 
317  // Sanity checks
318  if (X.size() != spfe.pole.size())
319  QDP_error_exit("%s : sanity failure, internal size not getSPartFracRoot size", __func__);
320 
321  if (X[0].size() != N5)
322  QDP_error_exit("%s : sanity failure, internal size not N5", __func__);
323 
324  // Weight solns
325  // Loop over each 5d slice
326  for(int j=0; j < N5; ++j)
327  {
328  tmp[j][M->subset()] = spfe.norm * getPhi()[n][j];
329  for(int i=0; i < X.size(); ++i)
330  tmp[j][M->subset()] += spfe.res[i] * X[i][j];
331  }
332 
333  // Action on the subset
334  action_m += norm2(tmp, M->subset());
335  }
336  }
337 
338  write(xml_out, "n_m_count", n_m_count);
339  write(xml_out, "S_m", action_m);
340  Double action = action_m;
341  write(xml_out, "S", action);
342  pop(xml_out);
343 
344  END_CODE();
345 
346  return action;
347  }
348 
349 
350  protected:
351  //! Get at fermion action
352  virtual const WilsonTypeFermAct5D<Phi,P,Q>& getFermAct(void) const = 0;
353 
354  //! Get inverter params
355  virtual const GroupXML_t& getActionInvParams(void) const = 0;
356 
357  //! Get inverter params
358  virtual const GroupXML_t& getForceInvParams(void) const = 0;
359 
360  //! Accessor for pseudofermion (read only)
361  virtual const multi1d< multi1d<Phi> >& getPhi(void) const = 0;
362 
363  //! mutator for pseudofermion
364  virtual multi1d< multi1d<Phi> >& getPhi(void) = 0;
365 
366  //! Return number of pseudofermions
367  virtual int getNPF() const = 0;
368 
369  //! Return the partial fraction expansion for the force calc
370  virtual const RemezCoeff_t& getFPFE() const = 0;
371 
372  //! Return the partial fraction expansion for the action calc
373  virtual const RemezCoeff_t& getSPFE() const = 0;
374 
375  //! Return the partial fraction expansion for the heat-bath
376  virtual const RemezCoeff_t& getSIPFE() const = 0;
377  };
378 
379 
380  //-------------------------------------------------------------------------------------------
381  //! Exact 1 flavor unpreconditioned fermact monomial living in extra dimensions
382  /*! @ingroup monomial
383  *
384  * Exact 1 flavor unpreconditioned fermact monomial.
385  */
386  template<typename P, typename Q, typename Phi>
388  {
389  public:
390  //! virtual destructor:
392 
393  //! Compute the total action
394  virtual Double S(const AbsFieldState<P,Q>& s)
395  {
396  START_CODE();
397 
398  XMLWriter& xml_out=TheXMLLogWriter::Instance();
399  push(xml_out, "OneFlavorRatExactUnprecWilsonTypeFermMonomial5D");
400 
401  Double action = this->S_subset(s);
402 
403  write(xml_out, "S", action);
404  pop(xml_out);
405 
406  END_CODE();
407 
408  return action;
409  }
410 
411 
412  protected:
413  //! Get at fermion action
414  virtual const UnprecWilsonTypeFermAct5D<Phi,P,Q>& getFermAct(void) const = 0;
415 
416  //! Get inverter params
417  virtual const GroupXML_t& getActionInvParams(void) const = 0;
418 
419  //! Get inverter params
420  virtual const GroupXML_t& getForceInvParams(void) const = 0;
421 
422  //! Return number of roots in used
423  virtual int getNPF() const = 0;
424 
425  //! Accessor for pseudofermion (read only)
426  virtual const multi1d< multi1d<Phi> >& getPhi(void) const = 0;
427 
428  //! mutator for pseudofermion
429  virtual multi1d< multi1d<Phi> >& getPhi(void) = 0;
430 
431  //! Return the partial fraction expansion for the force calc
432  virtual const RemezCoeff_t& getFPFE() const = 0;
433 
434  //! Return the partial fraction expansion for the action calc
435  virtual const RemezCoeff_t& getSPFE() const = 0;
436 
437  //! Return the partial fraction expansion for the heat-bath
438  virtual const RemezCoeff_t& getSIPFE() const = 0;
439 
440  };
441 
442 
443  //-------------------------------------------------------------------------------------------
444  //! Exact 1 flavor even-odd preconditioned fermact monomial living in extra dimensions
445  /*! @ingroup monomial
446  *
447  * Exact 1 flavor even-odd preconditioned fermact monomial.
448  * Can supply a default dsdq algorithm
449  */
450  template<typename P, typename Q, typename Phi>
452  {
453  public:
454  //! virtual destructor:
456 
457  //! Even even contribution (eg ln det Clover)
458  virtual Double S_even_even(const AbsFieldState<P,Q>& s) = 0;
459 
460  //! Compute the odd odd contribution (eg
462  {
463  return this->S_subset(s);
464  }
465 
466  //! Compute the total action
468  {
469  START_CODE();
470 
471  XMLWriter& xml_out=TheXMLLogWriter::Instance();
472  push(xml_out, "OneFlavorRatExactEvenOddPrecWilsonTypeFermMonomial5D");
473 
474  Double action_e = S_even_even(s);
475  Double action_o = S_odd_odd(s);
476  Double action = action_e + action_o;
477 
478  write(xml_out, "S_even_even", action_e);
479  write(xml_out, "S_odd_odd", action_o);
480  write(xml_out, "S", action);
481  pop(xml_out);
482 
483  END_CODE();
484 
485  return action;
486  }
487 
488  protected:
489  //! Get at fermion action
491 
492  //! Get inverter params
493  virtual const GroupXML_t& getActionInvParams(void) const = 0;
494 
495  //! Get inverter params
496  virtual const GroupXML_t& getForceInvParams(void) const = 0;
497 
498  //! Return number of roots in used
499  virtual int getNPF() const = 0;
500 
501  //! Accessor for pseudofermion (read only)
502  virtual const multi1d< multi1d<Phi> >& getPhi(void) const = 0;
503 
504  //! mutator for pseudofermion
505  virtual multi1d< multi1d<Phi> >& getPhi(void) = 0;
506 
507  //! Return the partial fraction expansion for the force calc
508  virtual const RemezCoeff_t& getFPFE() const = 0;
509 
510  //! Return the partial fraction expansion for the action calc
511  virtual const RemezCoeff_t& getSPFE() const = 0;
512 
513  //! Return the partial fraction expansion for the heat-bath
514  virtual const RemezCoeff_t& getSIPFE() const = 0;
515  };
516 
517  //-------------------------------------------------------------------------------------------
518  //! Exact 1 flavor even-odd preconditioned fermact monomial living in extra dimensions
519  /*! @ingroup monomial
520  *
521  * Exact 1 flavor even-odd preconditioned fermact monomial.
522  * Can supply a default dsdq algorithm
523  */
524  template<typename P, typename Q, typename Phi>
527  {
528  public:
529  //! virtual destructor:
531 
532  //! Even even contribution (eg ln det Clover)
534  return Double(0);
535  }
536 
537 
538  protected:
539  //! Get at fermion action
541 
542  };
543 
544 }
545 
546 
547 #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 > & getFermAct() const =0
Get at fermion action.
Exact 1 flavor even-odd preconditioned fermact monomial living in extra dimensions.
virtual const RemezCoeff_t & getSPFE() const =0
Return the partial fraction expansion for the action calc.
virtual const GroupXML_t & getActionInvParams(void) const =0
Get inverter params.
virtual const RemezCoeff_t & getFPFE() const =0
Return the partial fraction expansion for the force calc.
virtual const EvenOddPrecWilsonTypeFermAct5D< Phi, P, Q > & getFermAct() const =0
Get at fermion action.
virtual multi1d< multi1d< Phi > > & getPhi(void)=0
mutator for pseudofermion
virtual const multi1d< multi1d< Phi > > & getPhi(void) const =0
Accessor for pseudofermion (read only)
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 GroupXML_t & getForceInvParams(void) const =0
Get inverter params.
virtual int getNPF() const =0
Return number of roots in used.
virtual const RemezCoeff_t & getSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
Exact 1 flavor unpreconditioned fermact monomial living in extra dimensions.
virtual const multi1d< multi1d< Phi > > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual const RemezCoeff_t & getSPFE() const =0
Return the partial fraction expansion for the action calc.
virtual const UnprecWilsonTypeFermAct5D< Phi, P, Q > & getFermAct(void) const =0
Get at fermion action.
virtual int getNPF() const =0
Return number of roots in used.
virtual const RemezCoeff_t & getSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
virtual multi1d< multi1d< Phi > > & getPhi(void)=0
mutator for pseudofermion
virtual const GroupXML_t & getForceInvParams(void) const =0
Get inverter params.
virtual Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
virtual const GroupXML_t & getActionInvParams(void) const =0
Get inverter params.
virtual const RemezCoeff_t & getFPFE() const =0
Return the partial fraction expansion for the force calc.
Exact 1 flavor fermact monomial in extra dimensions.
virtual const RemezCoeff_t & getSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
virtual void dsdq(P &F, const AbsFieldState< P, Q > &s)
Compute dsdq for the system...
virtual const WilsonTypeFermAct5D< Phi, P, Q > & getFermAct(void) const =0
Get at fermion action.
virtual int getNPF() const =0
Return number of pseudofermions.
virtual const GroupXML_t & getForceInvParams(void) const =0
Get inverter params.
virtual multi1d< multi1d< Phi > > & getPhi(void)=0
mutator for pseudofermion
virtual const GroupXML_t & getActionInvParams(void) const =0
Get inverter params.
virtual const RemezCoeff_t & getSPFE() const =0
Return the partial fraction expansion for the action calc.
virtual const multi1d< multi1d< Phi > > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual void setInternalFields(const Monomial< P, Q > &m)
Copy internal fields.
virtual Double S(const AbsFieldState< P, Q > &s)=0
Compute the total action.
virtual const RemezCoeff_t & getFPFE() const =0
Return the partial fraction expansion for the force calc.
virtual Double S_subset(const AbsFieldState< P, Q > &s) const
Compute action on the appropriate subset.
virtual void refreshInternalFields(const AbsFieldState< P, Q > &s)
Refresh pseudofermions.
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