CHROMA
syssolver_linop_clover_qphix_w.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /*! \file
3  * \brief Solve a M*psi=chi linear system by BiCGStab
4  */
5 
6 #ifndef __syssolver_linop_clover_qphix_h__
7 #define __syssolver_linop_clover_qphix_h__
8 
9 #include "chroma_config.h"
10 
11 #ifdef BUILD_QPHIX
12 #include "handle.h"
13 #include "state.h"
14 #include "syssolver.h"
15 #include "linearop.h"
16 #include "actions/ferm/fermbcs/simple_fermbc.h"
17 #include "actions/ferm/fermstates/periodic_fermstate.h"
18 #include "actions/ferm/invert/qphix/syssolver_qphix_clover_params.h"
19 #include "actions/ferm/linop/clover_term_w.h"
20 #include "actions/ferm/linop/eoprec_clover_linop_w.h"
21 #include "meas/gfix/temporal_gauge.h"
22 #include "io/aniso_io.h"
23 #include <string>
24 
25 #include "chroma_config.h"
26 #include "util/gauge/reunit.h"
27 #include "io/aniso_io.h"
28 
29 // Header files from the dslash package
30 #include "qphix/geometry.h"
31 #include "qphix/qdp_packer.h"
32 #include "qphix/clover.h"
33 #include "qphix/invcg.h"
34 
35 #include "qphix/invbicgstab.h"
36 
37 #include "actions/ferm/invert/qphix/qphix_vec_traits.h"
38 #include "qphix_singleton.h"
39 using namespace QDP;
40 
41 namespace Chroma
42 {
43 
44  //! Richardson system solver namespace
45  namespace LinOpSysSolverQPhiXCloverEnv
46  {
47  //! Register the syssolver
48  bool registerAll();
49 
50 
51  }
52 
53 
54 
55  //! Solve a Clover Fermion System using the QPhiX inverter
56  /*! \ingroup invert
57  *** WARNING THIS SOLVER WORKS FOR Clover FERMIONS ONLY ***
58  */
59  using namespace QPhiXVecTraits;
60 
61  template<typename T, typename U>
62  class LinOpSysSolverQPhiXClover : public LinOpSystemSolver<T>
63  {
64  public:
65  typedef multi1d<U> Q;
66  typedef typename WordType<T>::Type_t REALT;
67  typedef typename QPhiX::Geometry<REALT,VecTraits<REALT>::Vec,VecTraits<REALT>::Soa,VecTraits<REALT>::compress12>::FourSpinorBlock QPhiX_Spinor;
68  typedef typename QPhiX::Geometry<REALT,VecTraits<REALT>::Vec,VecTraits<REALT>::Soa,VecTraits<REALT>::compress12>::SU3MatrixBlock QPhiX_Gauge;
69  typedef typename QPhiX::Geometry<REALT,VecTraits<REALT>::Vec,VecTraits<REALT>::Soa,VecTraits<REALT>::compress12>::CloverBlock QPhiX_Clover;
70 
71 
72  //! Constructor
73  /*!
74  * \param M_ Linear operator ( Read )
75  * \param invParam inverter parameters ( Read )
76  */
77  LinOpSysSolverQPhiXClover(Handle< LinearOperator<T> > A_,
78  Handle< FermState<T,Q,Q> > state_,
79  const SysSolverQPhiXCloverParams& invParam_) :
80  A(A_), invParam(invParam_), clov(new CloverTermT<T, U>()), invclov(new CloverTermT<T, U>())
81  {
82  QDPIO::cout << "LinOpSysSolverQPhiXClover:" << std::endl;
83  QDPIO::cout << "AntiPeriodicT is: " << invParam.AntiPeriodicT << std::endl;
84 
85 
86  QDPIO::cout << "Veclen is " << VecTraits<REALT>::Vec << std::endl;
87  QDPIO::cout << "Soalen is " << VecTraits<REALT>::Soa << std::endl;
88  if ( VecTraits<REALT>::Soa > VecTraits<REALT>::Vec ) {
89  QDPIO::cerr << "PROBLEM: Soalen > Veclen. Please set soalen appropriately (<=VECLEN) at compile time" << std::endl;
90  QDP_abort(1);
91  }
92 
93  Q u(Nd);
94  for(int mu=0; mu < Nd; mu++) {
95  u[mu] = state_->getLinks()[mu];
96  }
97 
98  // Set up aniso coefficients
99  multi1d<Real> aniso_coeffs(Nd);
100  for(int mu=0; mu < Nd; mu++) aniso_coeffs[mu] = Real(1);
101 
102  bool anisotropy = invParam.CloverParams.anisoParam.anisoP;
103  if( anisotropy ) {
104  aniso_coeffs = makeFermCoeffs( invParam.CloverParams.anisoParam );
105  }
106 
107  double t_boundary=(double)(1);
108  // NB: In this case, the state will have boundaries applied.
109  // So we only need to apply our own boundaries if Compression is enabled
110  //
111  if (invParam.AntiPeriodicT) {
112  t_boundary=(double)(-1);
113  // Flip off the boundaries -- Dslash expects them off..
114  u[3] *= where(Layout::latticeCoordinate(3) == (Layout::lattSize()[3]-1),
115  Real(t_boundary), Real(1));
116  }
117 
118  QDPIO::cout << "About to grap a Dslash" << std::endl;
119  const QPhiX::QPhiXCLIArgs& QPhiXParams = TheQPhiXParams::Instance();
120  cbsize_in_blocks = rb[0].numSiteTable() / VecTraits<REALT>::Soa;
121 
122 #ifdef QDP_IS_QDPJIT
123  int pad_xy = 0;
124  int pad_xyz = 0;
125 #else
126  int pad_xy = QPhiXParams.getPxy();
127  int pad_xyz = QPhiXParams.getPxyz();
128 #endif
129 
130  geom = new QPhiX::Geometry<REALT, VecTraits<REALT>::Vec, VecTraits<REALT>::Soa,VecTraits<REALT>::compress12>(Layout::subgridLattSize().slice(),
131  QPhiXParams.getBy(),
132  QPhiXParams.getBz(),
133  QPhiXParams.getNCores(),
134  QPhiXParams.getSy(),
135  QPhiXParams.getSz(),
136  pad_xy,
137  pad_xyz,
138  QPhiXParams.getMinCt());
139 
140 
141 #ifndef QDP_IS_QDPJIT
142  QDPIO::cout << " Allocating p and c" << std::endl << std::flush ;
143  psi_qphix=(QPhiX_Spinor *)geom->allocCBFourSpinor();
144  chi_qphix=(QPhiX_Spinor *)geom->allocCBFourSpinor();
145 #else
146  psi_qphix = nullptr;
147  chi_qphix = nullptr;
148 #endif
149  QDPIO::cout << " Allocating Clover" << std::endl << std::flush ;
150  QPhiX_Clover* A_cb0=(QPhiX_Clover *)geom->allocCBClov();
151  QPhiX_Clover* A_cb1=(QPhiX_Clover *)geom->allocCBClov();
152  clov_packed[0] = A_cb0;
153  clov_packed[1] = A_cb1;
154 
155  QDPIO::cout << " Allocating CloverInv " << std::endl << std::flush ;
156  QPhiX_Clover* A_inv_cb0=(QPhiX_Clover *)geom->allocCBClov();
157  QPhiX_Clover* A_inv_cb1=(QPhiX_Clover *)geom->allocCBClov();
158  invclov_packed[0] = A_inv_cb0;
159  invclov_packed[1] = A_inv_cb1;
160 
161  // Pack the gauge field
162  QDPIO::cout << "Packing gauge field..." << std::endl << std::flush ;
163  QPhiX_Gauge* packed_gauge_cb0=(QPhiX_Gauge *)geom->allocCBGauge();
164  QPhiX_Gauge* packed_gauge_cb1=(QPhiX_Gauge *)geom->allocCBGauge();
165 
166  QPhiX::qdp_pack_gauge<>(u, packed_gauge_cb0,packed_gauge_cb1, *geom);
167  u_packed[0] = packed_gauge_cb0;
168  u_packed[1] = packed_gauge_cb1;
169 
170 
171  QDPIO::cout << "Creating Clover Term" << std::endl;
172  CloverTerm clov_qdp;
173  clov->create(state_, invParam.CloverParams);
174  QDPIO::cout << "Inverting Clover Term" << std::endl;
175  invclov->create(state_, invParam.CloverParams, (*clov));
176  for(int cb=0; cb < 2; cb++) {
177  invclov->choles(cb);
178  }
179  QDPIO::cout << "Done" << std::endl;
180  QDPIO::cout << "Packing Clover term..." << std::endl;
181 
182  for(int cb=0; cb < 2; cb++) {
183  QPhiX::qdp_pack_clover<>((*invclov), invclov_packed[cb], *geom, cb);
184  }
185 
186  for(int cb=0; cb < 2; cb++) {
187  QPhiX::qdp_pack_clover<>((*clov), clov_packed[cb], *geom, cb);
188  }
189  QDPIO::cout << "Done" << std::endl;
190 
191  QDPIO::cout << "Creating the Even Odd Operator" << std::endl;
192  M=new QPhiX::EvenOddCloverOperator<REALT,VecTraits<REALT>::Vec,VecTraits<REALT>::Soa,VecTraits<REALT>::compress12>(u_packed,
193  clov_packed[1],
194  invclov_packed[0],
195  geom,
196  t_boundary,
197  toDouble(aniso_coeffs[0]),
198  toDouble(aniso_coeffs[3]));
199 
200 
201  switch( invParam.SolverType ) {
202  case CG:
203  {
204  QDPIO::cout << "Creating the CG Solver" << std::endl;
205  cg_solver = new QPhiX::InvCG<REALT,VecTraits<REALT>::Vec, VecTraits<REALT>::Soa, VecTraits<REALT>::compress12>((*M), invParam.MaxIter);
206  }
207  break;
208 
209  case BICGSTAB:
210  {
211  QDPIO::cout << "Creating the BiCGStab Solver" << std::endl;
212  bicgstab_solver = new QPhiX::InvBiCGStab<REALT,VecTraits<REALT>::Vec, VecTraits<REALT>::Soa, VecTraits<REALT>::compress12>((*M), invParam.MaxIter);
213  }
214  break;
215  default:
216  QDPIO::cerr << "UNKNOWN Solver Type" << std::endl;
217  QDP_abort(1);
218  }
219  }
220 
221 
222 
223  //! Destructor is automatic
224  ~LinOpSysSolverQPhiXClover()
225  {
226 
227  // Need to unalloc all the memory...
228  QDPIO::cout << "Destructing" << std::endl;
229 
230 #ifndef QDP_IS_QDPJIX
231  geom->free(psi_qphix);
232  geom->free(chi_qphix);
233 #endif
234  psi_qphix = nullptr;
235  chi_qphix = nullptr;
236 
237  geom->free(invclov_packed[0]);
238  geom->free(invclov_packed[1]);
239  invclov_packed[0] = nullptr;
240  invclov_packed[1] = nullptr;
241 
242  geom->free(clov_packed[0]);
243  geom->free(clov_packed[1]);
244  clov_packed[0] = nullptr;
245  clov_packed[1] = nullptr;
246 
247  geom->free(u_packed[0]);
248  geom->free(u_packed[1]);
249  u_packed[0] = nullptr;
250  u_packed[1] = nullptr;
251 
252  delete geom;
253 
254  // Evertyhing else freed as handles freed.
255  }
256 
257  //! Return the subset on which the operator acts
258  const Subset& subset() const {return A->subset();}
259 
260 
261 
262  //! Solver the linear system
263  /*!
264  * \param psi solution ( Modify )
265  * \param chi source ( Read )
266  * \return syssolver results
267  */
268  SystemSolverResults_t operator()(T& psi, const T& chi) const
269  {
270  /* Factories here later? */
271  SystemSolverResults_t res;
272  switch( invParam.SolverType ) {
273  case CG:
274  {
275  res = cgnrSolve(psi,chi);
276  }
277  break;
278  case BICGSTAB:
279  {
280  res = biCGStabSolve(psi,chi);
281  }
282  break;
283  default:
284  QDPIO::cout << "Unknown Solver " << std::endl;
285  break;
286  }
287 
288  return res;
289  }
290 
291 
292 
293 
294 
295  private:
296  // Hide default constructor
297  LinOpSysSolverQPhiXClover() {}
298 
299 
300 
301  Handle< LinearOperator<T> > A;
302  const SysSolverQPhiXCloverParams invParam;
303  Handle< CloverTermT<T, U> > clov;
304  Handle< CloverTermT<T, U> > invclov;
305 
306  QPhiX::Geometry<REALT, VecTraits<REALT>::Vec, VecTraits<REALT>::Soa, VecTraits<REALT>::compress12>* geom;
307 
308  Handle< QPhiX::EvenOddCloverOperator<REALT, VecTraits<REALT>::Vec, VecTraits<REALT>::Soa, VecTraits<REALT>::compress12> > M;
309 
310 
311  Handle< QPhiX::InvCG<REALT,VecTraits<REALT>::Vec, VecTraits<REALT>::Soa, VecTraits<REALT>::compress12> > cg_solver;
312 
313  Handle< QPhiX::InvBiCGStab<REALT,VecTraits<REALT>::Vec, VecTraits<REALT>::Soa, VecTraits<REALT>::compress12> > bicgstab_solver;
314 
315  QPhiX_Clover* invclov_packed[2];
316  QPhiX_Clover* clov_packed[2];
317  QPhiX_Gauge* u_packed[2];
318 
319  mutable QPhiX_Spinor* psi_qphix;
320  mutable QPhiX_Spinor* chi_qphix;
321  size_t cbsize_in_blocks;
322 
323  SystemSolverResults_t cgnrSolve(T& psi, const T& chi) const
324  {
325 
326  SystemSolverResults_t res;
327  psi = zero;
328 
329  // This is LinOpSolver so CGNE
330  // Step 1 tranform RHS: chi -> M^\dagger chi
331  // Move this into library later.
332  T mdag_chi;
333  (*A)(mdag_chi, chi, MINUS);
334 
335 
336  // QDPIO::cout << "Allocating Spinor fields" << std::endl;
337  // Pack Spinors psi and chi
338 #ifndef QDP_IS_QDPJIT
339  QPhiX::qdp_pack_cb_spinor<>(psi, psi_qphix, *geom,1);
340  QPhiX::qdp_pack_cb_spinor<>(mdag_chi, chi_qphix, *geom,1);
341 #else
342  psi_qphix = (QPhiX_Spinor*)(psi.getFjit()) + cbsize_in_blocks;
343  chi_qphix = (QPhiX_Spinor*)(chi.getFjit()) + cbsize_in_blocks;
344 #endif
345 
346  double rsd_final;
347  unsigned long site_flops=0;
348  unsigned long mv_apps=0;
349  int my_isign=1;
350 
351  double start = omp_get_wtime();
352  (*cg_solver)(psi_qphix,chi_qphix, toDouble(invParam.RsdTarget), res.n_count, rsd_final, site_flops, mv_apps, my_isign,invParam.VerboseP);
353  double end = omp_get_wtime();
354 
355  QDPIO::cout << "QPHIX_CLOVER_CG_SOLVER: " << res.n_count << " iters, rsd_sq_final=" << rsd_final << std::endl;
356 #ifndef QDP_IS_QDPJIT
357  QPhiX::qdp_unpack_cb_spinor<>(psi_qphix, psi,*geom,1);
358 #endif
359 
360  // Chi Should now hold the result spinor
361  // Check it against chroma.
362  T r = chi;
363  T tmp;
364  (*A)(tmp, psi, PLUS);
365  r[ A->subset() ] -= tmp;
366 
367  Double r2 = norm2(r,A->subset());
368  Double b2 = norm2(chi, A->subset());
369  Double rel_resid = sqrt(r2/b2);
370 
371  QDPIO::cout << "QPHIX_CLOVER_CG_SOLVER: || r || / || b || = " << rel_resid << std::endl;
372 
373 
374  if ( !toBool ( rel_resid < invParam.RsdTarget*invParam.RsdToleranceFactor ) ) {
375  QDPIO::cout << "SOLVE FAILED" << std::endl;
376  QDP_abort(1);
377  }
378 
379 
380  int num_cb_sites = Layout::vol()/2;
381  unsigned long total_flops = (site_flops + (1320+504+1320+504+48)*mv_apps)*num_cb_sites;
382  double gflops = (double)(total_flops)/(1.0e9);
383 
384  double total_time = end - start;
385  QDPIO::cout << "QPHIX_CLOVER_CG_SOLVER: Solver Time="<< total_time <<" (sec) Performance=" << gflops / total_time << " GFLOPS" << std::endl;
386 
387  END_CODE();
388  return res;
389 
390  }
391 
392  SystemSolverResults_t biCGStabSolve(T& psi, const T& chi) const
393  {
394  SystemSolverResults_t res;
395  psi=zero;
396 
397  // Pack Spinors psi and chi
398 #ifndef QDP_IS_QDPJIT
399  QDPIO::cout << "Packing" << std::endl << std::flush ;
400  QPhiX::qdp_pack_cb_spinor<>(psi, psi_qphix, *geom,1);
401  QPhiX::qdp_pack_cb_spinor<>(chi, chi_qphix, *geom,1);
402  QDPIO::cout << "Done" << std::endl << std::flush;
403 #else
404  psi_qphix = (QPhiX_Spinor *)(psi.getFjit())+cbsize_in_blocks;
405  chi_qphix = (QPhiX_Spinor *)(chi.getFjit())+cbsize_in_blocks;
406 #endif
407 
408  double rsd_final;
409  unsigned long site_flops=0;
410  unsigned long mv_apps=0;
411 
412  QDPIO::cout << "Starting solve" << std::endl << std::flush ;
413  double start = omp_get_wtime();
414  (*bicgstab_solver)(psi_qphix,chi_qphix, toDouble(invParam.RsdTarget), res.n_count, rsd_final, site_flops, mv_apps, 1, invParam.VerboseP);
415  double end = omp_get_wtime();
416 
417  QDPIO::cout << "QPHIX_CLOVER_BICGSTAB_SOLVER: " << res.n_count << " iters, rsd_sq_final=" << rsd_final << std::endl;
418  // QPhiX::qdp_unpack_spinor<>(psi_s[0], psi_qphix, psi, *geom);
419 #ifndef QDP_IS_QDPJIT
420  QPhiX::qdp_unpack_cb_spinor<>( psi_qphix, psi, *geom,1);
421 #endif
422 
423  // Chi Should now hold the result spinor
424  // Check it against chroma.
425  T r = chi;
426  T tmp;
427  (*A)(tmp, psi, PLUS);
428  r[ A->subset() ] -= tmp;
429 
430  Double r2 = norm2(r,A->subset());
431  Double b2 = norm2(chi, A->subset());
432  Double rel_resid = sqrt(r2/b2);
433 
434  QDPIO::cout << "QPHIX_CLOVER_BICGSTAB_SOLVER: || r || / || b || = " << rel_resid << std::endl;
435 
436 #if 0
437  if ( !toBool ( rel_resid < invParam.RsdTarget*invParam.RsdToleranceFactor ) ) {
438  QDPIO::cout << "SOLVE FAILED" << std::endl;
439  QDP_abort(1);
440  }
441 #endif
442 
443  int num_cb_sites = Layout::vol()/2;
444  unsigned long total_flops = (site_flops + (1320+504+1320+504+48)*mv_apps)*num_cb_sites;
445  double gflops = (double)(total_flops)/(1.0e9);
446 
447  double total_time = end - start;
448  QDPIO::cout << "QPHIX_CLOVER_BICGSTAB_SOLVER: Solver Time="<< total_time <<" (sec) Performance=" << gflops / total_time << " GFLOPS" << std::endl;
449 
450  END_CODE();
451  return res;
452  }
453 
454  };
455 
456 
457 } // End namespace
458 
459 #endif // BUILD_QPHIX
460 #endif
461 
aniso_io.h
Anisotropy parameters.
chroma_config.h
END_CODE
#define END_CODE()
Definition: chromabase.h:65
Chroma::QDPCloverTermT::create
void create(Handle< FermState< T, multi1d< U >, multi1d< U > > > fs, const CloverFermActParams &param_)
Creation routine.
Definition: clover_term_qdp_w.h:245
clover_term_w.h
Include possibly optimized Clover terms.
mu
int mu
Definition: cool.cc:24
eoprec_clover_linop_w.h
Even-odd preconditioned Clover fermion linear operator.
Chroma::CG
@ CG
Definition: syssolver_qphix_clover_params.h:16
Chroma::BICGSTAB
@ BICGSTAB
Definition: syssolver_qphix_clover_params.h:17
handle.h
Class for counted reference semantics.
linearop.h
Linear Operators.
tmp
Double tmp
Definition: meslate.cc:60
Nd
Nd
Definition: meslate.cc:74
Chroma::AsqtadFermActEnv::registerAll
bool registerAll()
Register all the factories.
Definition: asqtad_fermact_s.cc:52
Chroma::SimpleBaryonSeqSourceEnv::operator()
multi1d< Hadron2PtContraction_t > operator()(const multi1d< LatticeColorMatrix > &u)
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::Q
LinOpSysSolverMGProtoClover::Q Q
Definition: syssolver_linop_clover_mg_proto.cc:64
Chroma::CloverTerm
QDPCloverTerm CloverTerm
Definition: clover_term_w.h:92
Chroma::MINUS
@ MINUS
Definition: chromabase.h:45
Chroma::PLUS
@ PLUS
Definition: chromabase.h:45
Chroma::CloverTermT
QDPCloverTermT< T, U > CloverTermT
Definition: clover_term_w.h:97
Chroma::A
A(A, psi, r, Ncb, PLUS)
Chroma::cb
int cb
Definition: invbicg.cc:120
Chroma::zero
Double zero
Definition: invbicg.cc:106
Chroma::makeFermCoeffs
multi1d< Real > makeFermCoeffs(const AnisoParam_t &aniso)
Make fermion coefficients.
Definition: aniso_io.cc:63
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
r
multi1d< LatticeFermion > r(Ncb)
chi
chi
Definition: pade_trln_w.cc:24
psi
psi
Definition: pade_trln_w.cc:191
periodic_fermstate.h
Periodic ferm state and a creator.
qphix_singleton.h
Fermion action factories.
qphix_vec_traits.h
reunit.h
Reunitarize in place a color matrix to SU(N)
simple_fermbc.h
Simple fermionic BC.
state.h
Support class for fermion actions and linear operators.
syssolver.h
Linear system solvers.
syssolver_qphix_clover_params.h
U
multi1d< LatticeColorMatrix > U
Definition: t_aniso_gaugeact.cc:11
T
LatticeFermion T
Definition: t_clover.cc:11
Q
multi1d< LatticeColorMatrix > Q
Definition: t_clover.cc:12
temporal_gauge.h
Axial gauge fixing.
Generated by doxygen 1.9.1