CHROMA
one_flavor_rat_monomial_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_monomial_w_h__
8 #define __one_flavor_rat_monomial_w_h__
9 
16 #include <typeinfo>
17 
18 namespace Chroma
19 {
20  //-------------------------------------------------------------------------------------------
21  //! Exact 1 flavor fermact monomial using rational polynomials
22  /*! @ingroup monomial
23  *
24  * Exact 1 flavor fermact monomial using Rational Polynomial.
25  * Preconditioning is not 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  /*! Actions 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  // Self Description/Encapsulation Rule
45  XMLWriter& xml_out = TheXMLLogWriter::Instance();
46  push(xml_out, "OneFlavorRatExactWilsonTypeFermMonomial");
47 
48  /**** Identical code for unprec and even-odd prec case *****/
49 
50  /* The pseudofermionic action is of the form S_pf = chi^dag N*D^(-1) chi
51  * where N(Q) and D(Q) are polynomials in the rescaled kernel:
52  *
53  * Q = M^dag*M
54  * M = is some operator
55  *
56  * The function x^(-alpha) is approximated by c*N(x)/D(x)
57  * with c=FRatNorm. This has to be taken into account when comparing
58  * the force with the standard HMC routine dsduf for alpha=1. To solve
59  * ( N(Q)/D(Q) )psi the rational function is reformed as a partial fraction
60  * expansion and a multishift solver used to find the solution.
61  *
62  * Need
63  * dS_f/dU = - \sum_i psi_i^dag * p_i * [d(M^dag)*M + M^dag*dM] * psi
64  * where psi_i = (M^dag*M + q_i)^(-1)*chi
65  *
66  * In Balint's notation, the result is
67  * \dot{S} = -\sum_i p_i [ X_i^dag*\dot{M}^\dag*Y_i - Y_i^dag\dot{M}*X_i]
68  * where
69  * X_i = (M^dag*M + q_i)^(-1)*chi Y_i = M*X_i
70  * In Robert's notation, X_i -> psi_i .
71  */
73 
74  // Create a state for linop
76 
77  // Need way to get gauge state from AbsFieldState<P,Q>
79 
80  // Get multi-shift system solver
82 
83  // Partial Fraction Expansion coeffs for force
84  const RemezCoeff_t& fpfe = getFPFE();
85 
86  multi1d<Phi> X;
87  Phi Y;
88 
89  P F_1;
90  F.resize(Nd);
91  F = zero;
92 
93  // Loop over all the pseudoferms
94  multi1d<int> n_count(getNPF());
95  QDPIO::cout << "num_pf = " << getNPF() << std::endl;
96 
97  for(int n=0; n < getNPF(); ++n)
98  {
99  // The multi-shift inversion
100  SystemSolverResults_t res = (*invMdagM)(X, fpfe.pole, getPhi()[n]);
101  n_count[n] = res.n_count;
102 
103  // Loop over solns and accumulate force contributions
104 
105 
106 #if 0
107 
108  P F_2;
109  P F_tmp(Nd);
110 
111  F_tmp = zero;
112  for(int i=0; i < X.size(); ++i)
113  {
114  (*lin)(Y, X[i], PLUS);
115 
116  // The d(M^dag)*M term
117  lin->deriv(F_1, X[i], Y, MINUS);
118 
119  // The M^dag*d(M) term
120  lin->deriv(F_2, Y, X[i], PLUS);
121  F_1 += F_2;
122 
123  // Reweight each contribution in partial fraction
124  for(int mu=0; mu < F.size(); mu++) {
125  F_tmp[mu] -= fpfe.res[i] * F_1[mu];
126  }
127  }
128  F += F_tmp;
129 #else
130  // New code with new force term
131  multi1d<Phi> Y(X.size());
132  for(int i=0; i < X.size(); i++) {
133  (*lin)(Y[i], X[i], PLUS);
134  Y[i]*= -fpfe.res[i];
135  }
136 
137  lin->derivMultipole(F_1, X, Y, MINUS);
138  F += F_1;
139  lin->derivMultipole(F_1, Y, X, PLUS);
140  F += F_1;
141 #endif
142  // Concious choice. Don't monitor forces by pole
143  // Just accumulate it
144 
145  }
146 
147  state->deriv(F);
148  write(xml_out, "n_count", n_count);
149  monitorForces(xml_out, "Forces", F);
150  pop(xml_out);
151 
152  END_CODE();
153  }
154 
155 
156  //! Refresh pseudofermions
157  /*!
158  * This routine calculates the pseudofermion field (chi) for the case
159  * of rational evolution
160  *
161  * chi = n(Q)*[d(Q)]^(-1) * eta
162  * Where: Q = M^dag*M
163  * d(Q) = (Q+q_1)*(Q+q_2)*...*(Q+q_m)
164  * n(Q) = (Q+p_1)*(Q+p_2)*... *(Q+p_m)
165  * m = HBRatDeg
166  *
167  * The rational function n(x)/d(x) is the optimal rational
168  * approximation to the inverse square root of N(x)/D(x) which in
169  * turn is the optimal rational approximation to x^(-alpha).
170  * Here, alpha = 1/2
171  *
172  * To solve {n(Q)*[d(Q)]^(-1) * eta} the partial fraction expansion is
173  * used in combination with a multishift solver.
174  */
176  {
177  START_CODE();
178 
179  // Heatbath all the fields
180 
181  // Self Description/Encapsulation Rule
182  XMLWriter& xml_out = TheXMLLogWriter::Instance();
183  push(xml_out, "OneFlavorRatExactWilsonTypeFermMonomial");
184 
185  // Get at the ferion action for piece i
187 
188  // Create a Connect State, apply fermionic boundaries
189  Handle< FermState<Phi,P,Q> > f_state(FA.createState(s.getQ()));
190 
191  // Create a linear operator
193 
194  // Get multi-shift system solver
195 #if 1
197 #else
199 #endif
200  // Partial Fraction Expansion coeffs for heat-bath
201  const RemezCoeff_t& sipfe = getSIPFE();
202 
203  // Loop over pseudoferms
204  getPhi().resize(getNPF());
205  multi1d<int> n_count(getNPF());
206  Phi eta;
207 
208  for(int n=0; n < getNPF(); ++n)
209  {
210  // Fill the eta field with gaussian noise
211  eta = zero;
212  gaussian(eta, M->subset());
213 
214  // Account for fermion BC by modifying the proposed field
215  FA.getFermBC().modifyF(eta);
216 
217  // Temporary: Move to correct normalisation
218  eta *= sqrt(0.5);
219 
220  // The multi-shift inversion
221 #if 1
222  multi1d<Phi> X;
223  SystemSolverResults_t res = (*invMdagM)(X, sipfe.pole, eta);
224 #else
225  SystemSolverResults_t res = (*invMdagM)(getPhi()[n], sipfe.norm, sipfe.res,sipfe.pole, eta);
226 #endif
227  n_count[n] = res.n_count;
228 
229  // Weight solns to make final PF field
230 #if 1
231  getPhi()[n][M->subset()] = sipfe.norm * eta;
232  for(int i=0; i < X.size(); ++i)
233  getPhi()[n][M->subset()] += sipfe.res[i] * X[i];
234 #endif
235  }
236 
237  write(xml_out, "n_count", n_count);
238  pop(xml_out);
239 
240  END_CODE();
241  }
242 
243  //! Copy pseudofermions if any
244  virtual void setInternalFields(const Monomial<P,Q>& m)
245  {
246  START_CODE();
247 
248  try {
250 
251  getPhi() = fm.getPhi();
252  }
253  catch(std::bad_cast) {
254  QDPIO::cerr << "Failed to cast input Monomial to OneFlavorRatExactWilsonTypeFermMonomial " << std::endl;
255  QDP_abort(1);
256  }
257 
258  END_CODE();
259  }
260 
261 
262  //! Compute the action on the appropriate subset
263  /*!
264  * This 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  Handle< FermState<Phi,P,Q> > bc_g_state = FA.createState(s.getQ());
287 
288  // Need way to get gauge state from AbsFieldState<P,Q>
289  Handle< DiffLinearOperator<Phi,P,Q> > lin(FA.linOp(bc_g_state));
290 
291  // Get multi-shift system solver
292 #if 0
294 #else
296 #endif
297 
298  // Partial Fraction Expansion coeffs for action
299  const RemezCoeff_t& spfe = getSPFE();
300 
301  // Compute energy
302  // Get X out here via multisolver
303 #if 1
304  multi1d<Phi> X;
305 #endif
306  // Loop over all the pseudoferms
307  multi1d<int> n_count(getNPF());
308  Double action = zero;
309  Phi psi;
310 
311  for(int n=0; n < getNPF(); ++n)
312  {
313 #if 0
314  // The multi-shift inversion
315  SystemSolverResults_t res = (*invMdagM)(psi, spfe.norm, spfe.res,spfe.pole, getPhi()[n]);
316 #else
317  // The multi-shift inversion
318  SystemSolverResults_t res = (*invMdagM)(X, spfe.pole, getPhi()[n]);
319 #endif
320  n_count[n] = res.n_count;
321  LatticeDouble site_S=zero;
322 
323  // Take a volume factor out - redefine zero point energy
324  // this constant should have no physical effect, but by
325  // making S fluctuate around 0, it should remove a volume
326  // factor from the energies
327  site_S[ lin->subset() ] = -Double(12);
328 #if 1
329  psi[lin->subset()] = spfe.norm * getPhi()[n];
330  for(int i=0; i < X.size(); ++i) {
331  psi[lin->subset()] += spfe.res[i] * X[i];
332  }
333 #endif
334  // Accumulate locally
335  site_S[ lin->subset()] += localNorm2(psi);
336 
337  // action += norm2(psi, lin->subset());
338 
339  // Sum
340  action += sum(site_S, lin->subset());
341  }
342 
343  write(xml_out, "n_count", n_count);
344  write(xml_out, "S", action);
345  pop(xml_out);
346 
347  END_CODE();
348 
349  return action;
350  }
351 
352 
353  protected:
354  //! Get at fermion action
355  virtual const WilsonTypeFermAct<Phi,P,Q>& getFermAct(void) const = 0;
356 
357  //! Get inverter params
358  virtual const GroupXML_t& getActionInvParams(void) const = 0;
359 
360  //! Get inverter params
361  virtual const GroupXML_t& getForceInvParams(void) const = 0;
362 
363  //! Return number of roots in used
364  virtual int getNPF() const = 0;
365 
366  //! Return the partial fraction expansion for the force calc
367  virtual const RemezCoeff_t& getFPFE() const = 0;
368 
369  //! Return the partial fraction expansion for the action calc
370  virtual const RemezCoeff_t& getSPFE() const = 0;
371 
372  //! Return the partial fraction expansion for the heat-bath
373  virtual const RemezCoeff_t& getSIPFE() const = 0;
374 
375  //! Accessor for pseudofermion (read only)
376  virtual const multi1d<Phi>& getPhi(void) const = 0;
377 
378  //! mutator for pseudofermion
379  virtual multi1d<Phi>& getPhi(void) = 0;
380 
381  };
382 
383 
384  //-------------------------------------------------------------------------------------------
385  //! Exact 1 flavor unpreconditioned fermact monomial
386  /*! @ingroup monomial
387  *
388  * Exact 1 flavor unpreconditioned fermact monomial.
389  */
390  template<typename P, typename Q, typename Phi>
392  {
393  public:
394  //! virtual destructor:
396 
397  //! Compute the total action
398  virtual Double S(const AbsFieldState<P,Q>& s)
399  {
400  START_CODE();
401 
402  // Self identification/encapsulation Rule
403  XMLWriter& xml_out = TheXMLLogWriter::Instance();
404  push(xml_out, "OneFlavorRatExactUnprecWilsonTypeFermMonomial");
405 
406  Double action = this->S_subset(s);
407 
408  write(xml_out, "S", action);
409  pop(xml_out);
410 
411  END_CODE();
412 
413  return action;
414  }
415 
416  protected:
417  //! Get at fermion action
418  virtual const WilsonTypeFermAct<Phi,P,Q>& getFermAct(void) const = 0;
419 
420  //! Get inverter params
421  virtual const GroupXML_t& getActionInvParams(void) const = 0;
422 
423  //! Get inverter params
424  virtual const GroupXML_t& getForceInvParams(void) const = 0;
425 
426  //! Return number of roots in used
427  virtual int getNPF() const = 0;
428 
429  //! Return the partial fraction expansion for the force calc
430  virtual const RemezCoeff_t& getFPFE() const = 0;
431 
432  //! Return the partial fraction expansion for the action calc
433  virtual const RemezCoeff_t& getSPFE() const = 0;
434 
435  //! Return the partial fraction expansion for the heat-bath
436  virtual const RemezCoeff_t& getSIPFE() const = 0;
437 
438  //! Accessor for pseudofermion (read only)
439  virtual const multi1d<Phi>& getPhi(void) const = 0;
440 
441  //! mutator for pseudofermion
442  virtual multi1d<Phi>& getPhi(void) = 0;
443  };
444 
445 
446  //-------------------------------------------------------------------------------------------
447  //! Exact 1 flavor even-odd preconditioned fermact monomial
448  /*! @ingroup monomial
449  *
450  * Exact 1 flavor even-odd preconditioned fermact monomial.
451  * Can supply a default dsdq algorithm
452  */
453  template<typename P, typename Q, typename Phi,
454  template<class, class,class> class EOFermActT>
456  {
457  public:
458  //! virtual destructor:
460 
461  //! Even even contribution (eg ln det Clover)
462  virtual Double S_even_even(const AbsFieldState<P,Q>& s) = 0;
463 
464  //! Compute the odd odd contribution (eg
466  {
467  return this->S_subset(s);
468  }
469 
470  //! Compute the total action
472  {
473  START_CODE();
474 
475  XMLWriter& xml_out = TheXMLLogWriter::Instance();
476  push(xml_out, "OneFlavorRatExactEvenOddPrecWilsonTypeFermMonomial");
477 
478  Double action_e = S_even_even(s);
479  Double action_o = S_odd_odd(s);
480  Double action = action_e + action_o;
481 
482  write(xml_out, "S_even_even", action_e);
483  write(xml_out, "S_odd_odd", action_o);
484  write(xml_out, "S", action);
485  pop(xml_out);
486 
487  END_CODE();
488 
489  return action;
490  }
491 
492  protected:
493  //! Get at fermion action
494  virtual const EOFermActT<Phi,P,Q>& getFermAct() const = 0;
495 
496  //! Get inverter params
497  virtual const GroupXML_t& getActionInvParams(void) const = 0;
498 
499  //! Get inverter params
500  virtual const GroupXML_t& getForceInvParams(void) const = 0;
501 
502  //! Return number of roots in used
503  virtual int getNPF() const = 0;
504 
505  //! Return the partial fraction expansion for the force calc
506  virtual const RemezCoeff_t& getFPFE() const = 0;
507 
508  //! Return the partial fraction expansion for the action calc
509  virtual const RemezCoeff_t& getSPFE() const = 0;
510 
511  //! Return the partial fraction expansion for the heat-bath
512  virtual const RemezCoeff_t& getSIPFE() const = 0;
513 
514  //! Accessor for pseudofermion (read only)
515  virtual const multi1d<Phi>& getPhi(void) const = 0;
516 
517  //! mutator for pseudofermion
518  virtual multi1d<Phi>& getPhi(void) = 0;
519  };
520 
521  template<typename P, typename Q, typename Phi>
524 
525  template<typename P, typename Q, typename Phi>
528 
529  //-------------------------------------------------------------------------------------------
530  //! Exact 1 flavor even-odd preconditioned fermact monomial constant determinant
531  // Can fill out the S_odd_odd piece
532 
533  /*! @ingroup monomial
534  *
535  * Exact 1 flavor even-odd preconditioned fermact monomial.
536  * Can supply a default dsdq algorithm
537  */
538  template<typename P, typename Q, typename Phi,
539  template<class,class,class> class EOFermActT>
541  public OneFlavorRatExactEOPrecWilsonTypeFermMonomialT<P,Q,Phi,EOFermActT>
542  {
543  public:
544  //! virtual destructor:
546 
547  //! Even even contribution (eg ln det Clover)
549  return Double(0);
550  }
551 
552  protected:
553  //! Replace thiw with PrecConstDet
554  virtual const EOFermActT<Phi,P,Q>& getFermAct() const = 0;
555  };
556 
557  template<typename P, typename Q, typename Phi>
560 
561  template<typename P, typename Q, typename Phi>
564 }
565 
566 
567 #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.
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 EOFermActT< Phi, P, Q > & getFermAct() const =0
Replace thiw with PrecConstDet.
virtual Double S_even_even(const AbsFieldState< P, Q > &s)
Even even contribution (eg ln det Clover)
Exact 1 flavor even-odd preconditioned fermact monomial.
virtual const multi1d< Phi > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual const GroupXML_t & getForceInvParams(void) const =0
Get inverter params.
virtual int getNPF() const =0
Return number of roots in used.
Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
virtual const EOFermActT< Phi, P, Q > & getFermAct() const =0
Get at fermion action.
virtual Double S_even_even(const AbsFieldState< P, Q > &s)=0
Even even contribution (eg ln det Clover)
virtual Double S_odd_odd(const AbsFieldState< P, Q > &s)
Compute the odd odd contribution (eg.
virtual const RemezCoeff_t & getSPFE() const =0
Return the partial fraction expansion for the action calc.
virtual const RemezCoeff_t & getFPFE() const =0
Return the partial fraction expansion for the force calc.
virtual const GroupXML_t & getActionInvParams(void) const =0
Get inverter params.
virtual multi1d< Phi > & getPhi(void)=0
mutator for pseudofermion
virtual const RemezCoeff_t & getSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
Exact 1 flavor unpreconditioned fermact monomial.
virtual const WilsonTypeFermAct< Phi, P, Q > & getFermAct(void) const =0
Get at fermion action.
virtual const multi1d< Phi > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual const GroupXML_t & getForceInvParams(void) const =0
Get inverter params.
virtual const RemezCoeff_t & getSPFE() const =0
Return the partial fraction expansion for the action calc.
virtual multi1d< Phi > & getPhi(void)=0
mutator for pseudofermion
virtual const RemezCoeff_t & getFPFE() const =0
Return the partial fraction expansion for the force calc.
virtual const GroupXML_t & getActionInvParams(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.
virtual Double S(const AbsFieldState< P, Q > &s)
Compute the total action.
Exact 1 flavor fermact monomial using rational polynomials.
virtual void setInternalFields(const Monomial< P, Q > &m)
Copy pseudofermions if any.
virtual int getNPF() const =0
Return number of roots in used.
virtual Double S_subset(const AbsFieldState< P, Q > &s) const
Compute the action on the appropriate subset.
virtual Double S(const AbsFieldState< P, Q > &s)=0
Compute the total action.
virtual const multi1d< Phi > & getPhi(void) const =0
Accessor for pseudofermion (read only)
virtual const WilsonTypeFermAct< Phi, P, Q > & getFermAct(void) const =0
Get at fermion action.
virtual const GroupXML_t & getForceInvParams(void) const =0
Get inverter params.
virtual const RemezCoeff_t & getSIPFE() const =0
Return the partial fraction expansion for the heat-bath.
virtual void refreshInternalFields(const AbsFieldState< P, Q > &s)
Refresh pseudofermions.
virtual void dsdq(P &F, const AbsFieldState< P, Q > &s)
Compute dsdq for the system...
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 multi1d< Phi > & getPhi(void)=0
mutator for pseudofermion
virtual const RemezCoeff_t & getFPFE() const =0
Return the partial fraction expansion for the force 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.
virtual MdagMMultiSystemSolverAccumulate< T > * mInvMdagMAcc(Handle< FermState< T, P, Q > > state, const GroupXML_t &invParam) const
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
multi1d< LatticeColorMatrix > P
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
int n_count
Definition: invbicg.cc:78
LinOpSysSolverMGProtoClover::Q Q
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
Double sum
Definition: qtopcor.cc:37
Remez algorithm coefficients.
Symmetric even-odd const determinant Wilson-like fermact.
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