CHROMA
unprec_zolo_nef_fermact_array_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Unpreconditioned NEF fermion action
3  */
4 
5 #include "chromabase.h"
9 
12 
14 
17 
18 namespace Chroma
19 {
20  //! Hooks to register the class with the fermact factory
21  namespace UnprecZoloNEFFermActArrayEnv
22  {
23  //! Callback function
24  WilsonTypeFermAct5D<LatticeFermion,
25  multi1d<LatticeColorMatrix>,
26  multi1d<LatticeColorMatrix> >* createFermAct5D(XMLReader& xml_in,
27  const std::string& path)
28  {
29  return new UnprecZoloNEFFermActArray(CreateFermStateEnv::reader(xml_in, path),
30  UnprecZoloNEFFermActArrayParams(xml_in, path));
31  }
32 
33  //! Callback function
34  /*! Differs in return type */
35  FermionAction<LatticeFermion,
36  multi1d<LatticeColorMatrix>,
37  multi1d<LatticeColorMatrix> >* createFermAct(XMLReader& xml_in,
38  const std::string& path)
39  {
40  return createFermAct5D(xml_in, path);
41  }
42 
43  //! Name to be used
44  const std::string name = "UNPRECONDITIONED_ZOLO_NEF";
45 
46  //! Local registration flag
47  static bool registered = false;
48 
49  //! Register all the factories
50  bool registerAll()
51  {
52  bool success = true;
53  if (! registered)
54  {
55  success &= Chroma::TheFermionActionFactory::Instance().registerObject(name, createFermAct);
57  registered = true;
58  }
59  return success;
60  }
61  }
62 
63 
64  //! Read parameters
66  const std::string& path)
67  {
68  XMLReader paramtop(xml, path);
69  try {
70  // Read the stuff for the action
71  read(paramtop, "OverMass", OverMass);
72  read(paramtop, "Mass", Mass);
73  read(paramtop, "b5", b5);
74  read(paramtop, "c5", c5);
75  read(paramtop, "N5", N5);
76 
77  if( paramtop.count("ApproximationType") == 1 )
78  {
79  read(paramtop, "ApproximationType", approximation_type);
80  }
81  else
82  {
83  // Default coeffs are unscaled tanh
85  }
86 
88  {
89  read(paramtop, "ApproxMin", ApproxMin);
90  read(paramtop, "ApproxMax", ApproxMax);
91  }
92  else
93  {
94  ApproxMin = ApproxMax = 0.0;
95  }
96  }
97  catch(const std::string& e) {
98  QDPIO::cerr << "Caught Exception : " << e << std::endl;
99  }
100  }
101 
102 
103  //! Read parameters
104  void read(XMLReader& xml, const std::string& path, UnprecZoloNEFFermActArrayParams& param)
105  {
107  param = tmp;
108  }
109 
110 
111 
112  void UnprecZoloNEFFermActArray::initCoeffs(multi1d<Real>& b5_arr,
113  multi1d<Real>& c5_arr) const
114  {
115  b5_arr.resize(params.N5);
116  c5_arr.resize(params.N5);
117 
118  Real approxMin = params.ApproxMin;
119  Real approxMax = params.ApproxMax;
120 
121  zolotarev_data *rdata;
122  Real epsilon;
123  Real scale_fac;
124 
125  switch(params.approximation_type)
126  {
128  epsilon = approxMin / approxMax;
129  QDPIO::cout << "Initing Linop with Zolotarev Coefficients: epsilon = " << epsilon << std::endl;
130  rdata = zolotarev(toFloat(epsilon), params.N5, 0);
131  scale_fac = approxMax;
132  break;
133 
135  epsilon = approxMin;
136  QDPIO::cout << "Initing Linop with Higham Rep tanh Coefficients" << std::endl;
137  rdata = higham(toFloat(epsilon), params.N5);
138  scale_fac = Real(1);
139  break;
140 
141  default:
142  // The std::map system should ensure that we never get here but
143  // just for style
144  QDPIO::cerr << "Unknown coefficient type: " << params.approximation_type
145  << std::endl;
146  QDP_abort(1);
147  }
148 
149  if( rdata->n != params.N5 ) {
150  QDPIO::cerr << "Error:rdata->n != N5" << std::endl;
151  QDP_abort(1);
152  }
153 
154  Real maxerr = rdata->Delta;
155 
156  multi1d<Real> gamma(params.N5);
157  for(int i=0; i < params.N5; i++) {
158  gamma[i] = Real(rdata->gamma[i])*scale_fac;
159  }
160 
161  zolotarev_free(rdata);
162 
163  QDPIO::cout << "Initing Zolo NEF Linop: N5=" << params.N5
164  << " b5+c5=" << params.b5+params.c5
165  << " b5-c5=" << params.b5-params.c5
166  << " epsilon=" << epsilon
167  << " scale=" << approxMax
168  << " Maxerr=" << maxerr << std::endl;
169 
170  for(int i=0; i < params.N5; i++) {
171  QDPIO::cout << "gamma[" << i << "] = " << gamma[i] << std::endl;
172  }
173 
174  for(int i = 0; i < params.N5; i++) {
175  Real omega = Real(1)/gamma[i];
176 
177  b5_arr[i] = Real(0.5)*( (omega + Real(1))*params.b5 + (omega - Real(1))*params.c5 );
178  c5_arr[i] = Real(0.5)*( (omega - Real(1))*params.b5 + (omega + Real(1))*params.c5 );
179 
180  QDPIO::cout << "i=" << i << " omega["<<i<<"]=" << omega
181  << " b5["<< i << "] ="<< b5_arr[i]
182  << " c5["<< i << "] ="<< c5_arr[i] << std::endl;
183  QDPIO::cout << "i=" << i << " b5["<<i<<"]+c5["<<i<<"]/omega["<<i<<"] ="
184  << (b5_arr[i]+c5_arr[i])/omega
185  << " b5["<<i<<"]-c5["<<i<<"]=" << b5_arr[i]-c5_arr[i] << std::endl;
186  }
187  }
188 
189 
190  //! Produce an unpreconditioned linear operator for this action with arbitrary quark mass
191  UnprecDWLikeLinOpBaseArray<LatticeFermion,
192  multi1d<LatticeColorMatrix>,
193  multi1d<LatticeColorMatrix> >*
195  const Real& m_q) const
196  {
197  multi1d<Real> b5_arr;
198  multi1d<Real> c5_arr;
199 
200  initCoeffs(b5_arr,c5_arr);
201 
202  return new UnprecNEFDWLinOpArray(state,params.OverMass,b5_arr,c5_arr,m_q,params.N5);
203  }
204 
205 
206 
207  // Given a complete propagator as a source, this does all the inversions needed
208  void
209  UnprecZoloNEFFermActArray::quarkProp(LatticePropagator& q_sol,
210  XMLWriter& xml_out,
211  const LatticePropagator& q_src,
212  int t_src, int j_decay,
214  const GroupXML_t& invParam,
215  QuarkSpinType quarkSpinType,
216  bool obsvP,
217  int& ncg_had) const
218  {
219  if (obsvP && (quarkSpinType == QUARK_SPIN_TYPE_FULL))
220  nef_quarkProp4(q_sol, xml_out, q_src, t_src, j_decay, *this, state, invParam, ncg_had);
221  else
222  {
224  quarkProp4(q_sol, xml_out, q_src, qprop, quarkSpinType, ncg_had);
225  }
226  }
227 
228 
229 }
Primary include file for CHROMA library code.
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
SystemSolver< LatticeFermion > * qprop(Handle< FermState< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > state, const GroupXML_t &invParam) const
Define quark propagator routine for 4D fermions.
Unpreconditioned domain-wall Dirac operator.
Unpreconditioned domain-wall Dirac operator.
Unpreconditioned NEF fermion action.
UnprecDWLikeLinOpBaseArray< T, P, Q > * unprecLinOp(Handle< FermState< T, P, Q > > state, const Real &m_q) const
Produce an unpreconditioned linear operator for this action with arbitrary quark mass.
void quarkProp(LatticePropagator &q_sol, XMLWriter &xml_out, const LatticePropagator &q_src, int t_src, int j_decay, Handle< FermState< T, P, Q > > state, const GroupXML_t &invParam, QuarkSpinType quarkSpinType, bool obsvP, int &ncg_had) const
Given a complete propagator as a source, this does all the inversions needed.
void initCoeffs(multi1d< Real > &b5_arr, multi1d< Real > &c5_arr) const
Hide =.
Wilson-like fermion actions.
Definition: fermact.orig.h:403
All ferm create-state method.
Fermion action factories.
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
QuarkSpinType
Quark spin type.
@ COEFF_TYPE_ZOLOTAREV
@ COEFF_TYPE_TANH_UNSCALED
void nef_quarkProp4(LatticePropagator &q_sol, XMLWriter &xml_out, const LatticePropagator &q_src, int t_src, int j_decay, const UnprecDWFermActBaseArray< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > &S_f, Handle< FermState< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > state, const GroupXML_t &invParam, int &ncg_had)
Given a complete propagator as a source, this does all the inversions needed.
void quarkProp4(LatticeStaggeredPropagator &q_sol, XMLWriter &xml_out, const LatticeStaggeredPropagator &q_src, const StaggeredTypeFermAct< LatticeStaggeredFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > &S_f, Handle< FermState< LatticeStaggeredFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > state, const GroupXML_t &invParam, QuarkSpinType quarkSpinType, int &ncg_had)
Given a complete propagator as a source, this does all the inversions needed.
int j_decay
Definition: meslate.cc:22
Handle< CreateFermState< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > reader(XMLReader &xml_in, const std::string &path)
Helper function for the CreateFermState readers.
int epsilon(int i, int j, int k)
FermionAction< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > * createFermAct(XMLReader &xml_in, const std::string &path)
Callback function.
WilsonTypeFermAct5D< LatticeFermion, multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > * createFermAct5D(XMLReader &xml_in, const std::string &path)
Callback function.
static bool registered
Local registration flag.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
int i
Definition: pbg5p_w.cc:55
Complex omega
Definition: invbicg.cc:97
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
::std::string string
Definition: gtest.h:1979
Full quark propagator solver for domain wall fermions.
Full quark propagator solver.
Hold group xml and type id.
CoeffType approximation_type
ZOLOTAREV | TANH | Other approximation coeffs.
Real ApproxMax
Approximate max eigenvalue of H_T.
Real OverMass
Mass of auxiliary Wilson action.
Real ApproxMin
Approximate min eigenvalue of H_T.
Unpreconditioned NEF domain-wall fermion linear operator.
Unpreconditioned Wilson fermion action.
Unpreconditioned NEF domain-wall fermion action.
void zolotarev_free(ZOLOTAREV_DATA *f)
ZOLOTAREV_DATA * higham(PRECISION epsilon, int n)
ZOLOTAREV_DATA * zolotarev(PRECISION epsilon, int n, int type)