CHROMA
eoprec_ovlap_contfrac5d_fermact_array_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Unpreconditioned extended-Overlap (5D) (Naryanan&Neuberger) action
3  */
4 
5 #include "chromabase.h"
6 
9 //#include "actions/ferm/linop/eoprec_ovlap_contfrac5d_linop_array_w.h"
14 
17 
18 #include "io/enum_io/enum_io.h"
19 #include "io/overlap_state_info.h"
20 
22 
23 namespace Chroma
24 {
25 
26  //! Hooks to register the class with the fermact factory
27  namespace EvenOddPrecOvlapContFrac5DFermActArrayEnv
28  {
29  //! Callback function
30  WilsonTypeFermAct5D<LatticeFermion,
31  multi1d<LatticeColorMatrix>,
32  multi1d<LatticeColorMatrix> >* createFermAct5D(XMLReader& xml_in,
33  const std::string& path)
34  {
37  }
38 
39  //! Callback function
40  /*! Differs in return type */
41  FermionAction<LatticeFermion,
42  multi1d<LatticeColorMatrix>,
43  multi1d<LatticeColorMatrix> >* createFermAct(XMLReader& xml_in,
44  const std::string& path)
45  {
46  return createFermAct5D(xml_in, path);
47  }
48 
49  //! Name to be used
50  const std::string name = "OVERLAP_CONTINUED_FRACTION_5D";
51 
52  //! Local registration flag
53  static bool registered = false;
54 
55  //! Register all the factories
56  bool registerAll()
57  {
58  bool success = true;
59  if (! registered)
60  {
61  success &= Chroma::TheFermionActionFactory::Instance().registerObject(name, createFermAct);
63  registered = true;
64  }
65  return success;
66  }
67  } // End Namespace EvenOddPrecOvlapContFrac5DFermActArrayEnv
68 
69 
70  //! Read XML
72  const std::string& path)
73  {
74  XMLReader in(xml, path);
75 
76  try
77  {
78  read(in, "Mass", Mass);
79  read(in, "RatPolyDeg", RatPolyDeg);
80  read(in, "OverMass", OverMass);
81 
82  if( in.count("ApproximationType") == 1 )
83  {
84  read(in, "ApproximationType", approximation_type);
85  }
86  else
87  {
88  // Default coeffs are unscaled tanh
90  }
91 
93  {
94  read(in, "ApproxMin", ApproxMin);
95  read(in, "ApproxMax", ApproxMax);
96  }
97  else
98  {
99  ApproxMin = ApproxMax = 0.0;
100  }
101  }
102  catch( const std::string &e ) {
103  QDPIO::cerr << "Caught Exception reading prec ContFrac Fermact params: " << e << std::endl;
104  QDP_abort(1);
105  }
106  }
107 
108  void read(XMLReader& xml_in, const std::string& path,
110  {
112  param = tmp;
113  }
114 
115  void write(XMLWriter& xml_out, const std::string& path, const EvenOddPrecOvlapContFrac5DFermActParams& p)
116  {
117  if ( path != "." ) {
118  push( xml_out, path);
119  }
120 
121  write(xml_out, "Mass", p.Mass);
122  write(xml_out, "OverMass", p.OverMass);
123  write(xml_out, "RatPolyDeg", p.RatPolyDeg);
124  write(xml_out, "ApproximationType", p.approximation_type);
125  write(xml_out, "ApproxMin", p.ApproxMin);
126  write(xml_out, "ApproxMax", p.ApproxMax);
127 
128  pop(xml_out);
129 
130  if( path != "." ) {
131  pop(xml_out);
132  }
133  }
134 
135 
136 
137  // Construct the action out of a parameter structure
139  Handle< CreateFermState<T,P,Q> > cfs_a_,
141  cfs(cfs_a_), params(params_)
142  {
143 
144  // WHAT IS BELOW ONLY WORKS FOR TYPE=0 approximations
145  // which is what we use. Forget TYPE=1
146  // the Tanh approximation (Higham) is of type TYPE=0
147  // We have two cases.
148  bool isEvenRatPolyDeg = ( params.RatPolyDeg % 2 == 0);
149 
150  if( isEvenRatPolyDeg ) {
151  N5 = params.RatPolyDeg+1;
152  isLastZeroP = true;
153  }
154  else {
155  N5 = params.RatPolyDeg;
156  isLastZeroP = false;
157  }
158  }
159 
160 
161  void
163  multi1d<Real>& alpha,
164  multi1d<Real>& beta) const
165  {
166  int type = 0;
167  zolotarev_data *rdata;
168  Real epsilon;
169 
170  Real approxMin = params.ApproxMin;
171  Real approxMax = params.ApproxMax;
172 
173  switch(params.approximation_type)
174  {
176  epsilon = approxMin / approxMax;
177  QDPIO::cout << "Initing Linop with Zolotarev Coefficients: epsilon = " << epsilon << std::endl;
178  rdata = zolotarev(toFloat(epsilon), params.RatPolyDeg, type);
179  scale_fac = Real(1) / approxMax;
180  break;
181 
183  epsilon = approxMin;
184  QDPIO::cout << "Initing Linop with Higham Rep tanh Coefficients" << std::endl;
185  rdata = higham(toFloat(epsilon), params.RatPolyDeg);
186  scale_fac = Real(1);
187  break;
188 
189  default:
190  // The std::map system should ensure that we never get here but
191  // just for style
192  QDPIO::cerr << "Unknown coefficient type: " << params.approximation_type
193  << std::endl;
194  QDP_abort(1);
195  }
196 
197  Real maxerr = (Real)(rdata->Delta);
198 
199  // Check N5 is good:
200  if( N5 != rdata->db ) {
201  QDPIO::cerr << "Mismatch between N5 and N5 from Coefficient Code" << std::endl;
202  QDPIO::cerr << "N5 = " << N5 << " rdata->db=" << rdata->db << std::endl;
203  QDP_abort(1);
204  }
205 
206 
207  // The coefficients from the continued fraction
208  beta.resize(N5);
209  for(int i = 0; i < N5; i++) {
210  beta[i] = rdata->beta[i];
211  }
212 
213  alpha.resize(N5-1);
214  for(int i=0; i < N5-1; i++) {
215  alpha[i] = Real(1);
216  }
217 
218  // The gamma's are the equivalence transforms
219  // There are N5-1 of them and they appear in the
220  // diagonal terms as gamma^2
221  // and in the off diagonal terms as gamma_i gamma_i+1
222  // except in the last one which is just gamma
223  //
224  // For the moment choose gamma_i = 1/sqrt(beta_i) */
225 
226  multi1d<Real> gamma(N5-1);
227  for(int i=0; i < N5-1; i++) {
228  gamma[i] = Real(1)/ sqrt( beta[i] );
229  }
230 
231  // Now perform the equivalence transformation
232  //
233  // On the diagonal coefficients
234  // Note that beta[N5-1] is NOT changed
235  for(int i=0; i < N5-1; i++) {
236  beta[i] = beta[i]*gamma[i]*gamma[i];
237  }
238 
239  // and the off diagonal ones
240  // from 0..N5-3 we have gamma_i gamma_i+1
241  // and on N5-2 we have gamma_i
242  for(int i=0; i < N5-2; i++) {
243  alpha[i] *= gamma[i]*gamma[i+1];
244  }
245  alpha[N5-2] *= gamma[N5-2];
246 
247  QDPIO::cout << "EvenOddPrecOvlapContfrac5d: " << std::endl
248  << " Degree="<< params.RatPolyDeg
249  << " N5=" << N5 << " scale=" << scale_fac
250  << " Mass=" << params.Mass << std::endl
251  << " OverMass=" << params.OverMass
252  << " IsLastZeroP=" << isLastZeroP << std::endl;
253 
254  QDPIO::cout << "Approximation on [-1,eps] U [eps,1] with eps = " << epsilon <<std::endl;
255 
256  QDPIO::cout << "Maximum error | R(x) - sgn(x) | <= Delta = " << maxerr << std::endl;
257  /*
258  for(int i=0; i < N5; i++) {
259  QDPIO::cout << "beta["<<i<<"] = " << beta[i] << std::endl;
260  }
261  for(int i=0; i < N5; i++) {
262  QDPIO::cout << "alpha["<<i<<"] = " << alpha[i] << std::endl;
263  }
264  */
265 
266  switch( params.approximation_type) {
268  QDPIO::cout << "Coefficients from Zolotarev" << std::endl;
269 
270  if(type == 0) {
271  QDPIO::cout << "Approximation type " << type << " with R(0) = 0"
272  << std::endl;
273  }
274  else {
275  QDPIO::cout << "Approximation type " << type << " with R(0) = infinity" << std::endl;
276  }
277 
278  break;
279  case COEFF_TYPE_TANH:
280  QDPIO::cout << "Coefficients from Higham Tanh representation" << std::endl;
281  break;
283  QDPIO::cout << "Coefficients from Unscaled Higham Tanh representation" << std::endl;
284  break;
285  default:
286  QDPIO::cerr << "Unknown coefficient type " << params.approximation_type
287  << std::endl;
288  QDP_abort(1);
289  break;
290  }
291 
292  // free the arrays allocated by Tony's Zolo
293  zolotarev_free(rdata);
294  }
295 
296  //! Produce a linear operator for this action
297  /*!
298  * \param state gauge field (Read)
299  */
301  multi1d<LatticeColorMatrix>,
302  multi1d<LatticeColorMatrix> >*
304  {
305  START_CODE();
306 
307  multi1d<Real> alpha;
308  multi1d<Real> beta;
309  Real scale_factor;
310 
311  init(scale_factor, alpha, beta);
312 
313  return new EvenOddPrecOvlapContFrac5DLinOpArray(state_,
314  params.Mass,
316  N5,
317  scale_factor,
318  alpha,
319  beta,
320  isLastZeroP);
321  }
322 
323 
324  //! Produce a Pauli-Villars linear operator for this action
325  /*!
326  * \param state gauge field (Read)
327  */
329  multi1d<LatticeColorMatrix>,
330  multi1d<LatticeColorMatrix> >*
332  {
333  multi1d<Real> alpha;
334  multi1d<Real> beta;
335  Real scale_factor;
336 
337  init(scale_factor, alpha, beta);
338 
339  // Hmm, not sure about what all the rescaling does to the PV....
340  return new EvenOddPrecOvlapContFrac5DPVLinOpArray(state_,
341  params.Mass,
343  N5,
344  scale_factor,
345  alpha,
346  beta,
347  isLastZeroP);
348  }
349 
350 
351 
352  //! Propagator of an un-preconditioned Extended-Overlap linear operator
353  /*! \ingroup qprop
354  *
355  * Propagator solver for DWF-like fermions
356  */
357  template<typename T, typename P, typename Q>
358  class ContFrac5DQprop : public SystemSolver<T>
359  {
360  public:
361  //! Constructor
362  /*!
363  * \param A_ Linear operator ( Read )
364  * \param Mass_ quark mass ( Read )
365  */
367  const Real& Mass_,
368  const SysSolverCGParams& invParam_) :
369  A(A_), Mass(Mass_), invParam(invParam_) {}
370 
371  //! Destructor is automatic
373 
374  //! Return the subset on which the operator acts
375  const Subset& subset() const {return all;}
376 
377  //! Solver the linear system
378  /*!
379  * \param psi quark propagator ( Modify )
380  * \param chi source ( Read )
381  * \return number of CG iterations
382  */
384  {
385  START_CODE();
386 
388  const int N5 = A->size();
389 
390  int G5 = Ns*Ns - 1;
391 
392  // Initialize the 5D fields
393  multi1d<T> chi5(N5);
394  multi1d<T> psi5(N5);
395 
396 // if( invType == "CG_INVERTER")
397  {
398  multi1d<T> tmp5_1(N5);
399  {
400  multi1d<T> tmp5_2(N5);
401  multi1d<T> tmp5_3(N5);
402 
403  chi5 = zero;
404  psi5 = zero;
405  tmp5_1 = zero;
406 
407  // We need to solve D_5 psi = (0,0,0,...,gamma_5 chi)^T
408  // Use both subsets
409  tmp5_1[N5-1] = Gamma(G5)*chi;
410 
411 
412  // tmp5_3_odd = Qoe Qee^{-1} S_e
413  A->evenEvenInvLinOp(tmp5_2, tmp5_1, PLUS);
414  A->oddEvenLinOp(tmp5_3, tmp5_2, PLUS);
415 
416 
417  // chi5_e = S_e
418  // chi5_o = S_o - QoeQee^{-1} S_e
419  for(int i=0; i < N5; i++) {
420  chi5[i][rb[0]] = tmp5_1[i];
421  chi5[i][rb[1]] = tmp5_1[i] - tmp5_3[i];
422  }
423  } // tmp5_2 and tmp5_3 go away
424 
425  psi5[N5-1][rb[1]] = psi;
426  (*A)(tmp5_1, chi5, MINUS);
427 
428  // Solve M^{+}M psi = M^{+} chi
429  res = InvCG2(*A, tmp5_1, psi5, invParam.RsdCG, invParam.MaxCG);
430 
431 
432  // psi[N5-1]_odd now holds the desired piece.
433 
434  // Reconstruct psi[N5-1]_e = Q_ee^{-1} S_e - Q_ee^{-1}Q_eo psi[N5-1]_o
435  // = Q_ee^{-1} ( S_e - Q_eo psi_o )
436  {
437 
438  // Dont need to initialise as the parts I use
439  // will be over written the other parts I ignore
440  multi1d<T> tmp5_2(N5);
441  multi1d<T> tmp5_3(N5);
442 
443  // tmp5_2_e = Qeo psi_o
444  A->evenOddLinOp(tmp5_2, psi5, PLUS);
445  for(int i=0; i < N5; i++) {
446 
447  // tmp5_3_e = S_e - Qeo psi_o
448  // = - tmp5_2_e
449  tmp5_3[i][rb[0]] = chi5[i] - tmp5_2[i];
450 
451  }
452 
453  // tmp5_1_e = Qee^{-1} ( S_e - Qeo psi_o )
454  A->evenEvenInvLinOp(tmp5_1, tmp5_3, PLUS);
455 
456  // psi5_e = tmp5_1
457  for(int i=0; i < N5; i++) {
458  psi5[i][rb[0]] = tmp5_1[i];
459  }
460  } // tmp5_2, tmp5_3 disappear
461  } // tmp5_1 disappears
462 // else
463 // {
464 // QDPIO::cerr << EvenOddPrecOvlapContFrac5DFermActArrayEnv::name
465 // << "Unsupported inverter type =" << invType << std::endl;
466 // QDP_abort(1);
467 // }
468 
469  if ( res.n_count == invParam.MaxCG )
470  QDP_error_exit("no convergence in the inverter", res.n_count);
471 
472  // Compute residual
473  {
474  multi1d<T> r(N5);
475  A->unprecLinOp(r, psi5, PLUS);
476  chi5 = zero;
477  chi5[N5-1] = Gamma(G5)*chi;
478  r -= chi5;
479  res.resid = sqrt(norm2(r));
480  }
481 
482  // Solution now lives in chi5
483 
484  // Multiply back in factor 2/(1-m) to return to
485  // (1/2)( 1 + m + (1-m) gamma_5 epsilon )
486  // normalisation
487  psi5[N5-1] *= Real(2)/(Real(1)-Mass);
488 
489  // Remove contact term
490  psi = psi5[N5-1] - chi;
491 
492  // Overall normalization
493  Real ftmp1 = Real(1) / (Real(1) - Mass);
494  psi *= ftmp1;
495 
496  END_CODE();
497 
498  return res;
499  }
500 
501  private:
502  // Hide default constructor
504 
506  Real Mass;
508  };
509 
510 
511  //! Propagator of an un-preconditioned Extended-Overlap linear operator
514  const GroupXML_t& invParam) const
515  {
516  std::istringstream is(invParam.xml);
517  XMLReader paramtop(is);
518 
519  return new ContFrac5DQprop<T,P,Q>(
521  getQuarkMass(),
522  SysSolverCGParams(paramtop,invParam.path));
523  }
524 
525 } // End namespace Chroma
Primary include file for CHROMA library code.
Propagator of an un-preconditioned Extended-Overlap linear operator.
ContFrac5DQprop(Handle< EvenOddPrecConstDetLinearOperatorArray< T, P, Q > > A_, const Real &Mass_, const SysSolverCGParams &invParam_)
Constructor.
Handle< EvenOddPrecConstDetLinearOperatorArray< T, P, Q > > A
const Subset & subset() const
Return the subset on which the operator acts.
SystemSolverResults_t operator()(T &psi, const T &chi) const
Solver the linear system.
Create a fermion connection state.
Definition: create_state.h:69
Even-odd preconditioned linear operator including derivatives for arrays.
5D continued fraction overlap action (Borici,Wenger, Edwards)
void init(Real &scale_fac, multi1d< Real > &alpha, multi1d< Real > &beta) const
Helper in construction.
EvenOddPrecConstDetLinearOperatorArray< T, P, Q > * linOp(Handle< FermState< T, P, Q > > state) const
Produce a linear operator for this action.
SystemSolver< LatticeFermion > * qprop(Handle< FermState< T, P, Q > > state, const GroupXML_t &invParam) const
Compute quark propagator over base type.
EvenOddPrecConstDetLinearOperatorArray< T, P, Q > * linOpPV(Handle< FermState< T, P, Q > > state) const
Produce a Pauli-Villars linear operator for this action.
Even-odd preconditioned Pauli-Villars Continued Fraction 5D.
Support class for fermion actions and linear operators.
Definition: state.h:94
Base class for quadratic matter actions (e.g., fermions)
Definition: fermact.h:53
Class for counted reference semantics.
Definition: handle.h:33
static T & Instance()
Definition: singleton.h:432
Linear system solvers.
Definition: syssolver.h:34
Wilson-like fermion actions.
Definition: fermact.orig.h:403
Enums.
Even-odd preconditioned Continued Fraction 5D.
Even-odd preconditioned Pauli-Villars Continued Fraction 5D.
All ferm create-state method.
Fermion action factories.
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
SystemSolverResults_t InvCG2(const LinearOperator< LatticeFermionF > &M, const LatticeFermionF &chi, LatticeFermionF &psi, const Real &RsdCG, int MaxCG)
Conjugate-Gradient (CGNE) algorithm for a generic Linear Operator.
Definition: invcg2.cc:240
@ COEFF_TYPE_ZOLOTAREV
@ COEFF_TYPE_TANH
@ COEFF_TYPE_TANH_UNSCALED
Params params
Conjugate-Gradient algorithm for a generic Linear Operator.
Handle< CreateFermState< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > reader(XMLReader &xml_in, const std::string &path)
Helper function for the CreateFermState readers.
WilsonTypeFermAct5D< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > * createFermAct5D(XMLReader &xml_in, const std::string &path)
Callback function.
FermionAction< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > * createFermAct(XMLReader &xml_in, const std::string &path)
Callback function.
int epsilon(int i, int j, int k)
static const LatticeInteger & beta(const int dim)
Definition: stag_phases_s.h:47
static const LatticeInteger & alpha(const int dim)
Definition: stag_phases_s.h:43
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
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)
int G5
Definition: pbg5p_w.cc:57
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
push(xml_out,"Condensates")
LinOpSysSolverMGProtoClover::T T
int i
Definition: pbg5p_w.cc:55
QDPEvenOddPrecOvlapContFrac5DLinOpArray EvenOddPrecOvlapContFrac5DLinOpArray
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
START_CODE()
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
Double zero
Definition: invbicg.cc:106
static QDP_ColorVector * in
::std::string string
Definition: gtest.h:1979
Include possibly optimized partfrac5d.
CoeffType approximation_type
ZOLOTAREV | TANH | Other approximation coeffs.
Hold group xml and type id.
Params for CG inverter.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Solve a CG1 system.
Double ftmp1
Definition: topol.cc:29
Unpreconditioned Wilson fermion action.
void zolotarev_free(ZOLOTAREV_DATA *f)
ZOLOTAREV_DATA * higham(PRECISION epsilon, int n)
ZOLOTAREV_DATA * zolotarev(PRECISION epsilon, int n, int type)