CHROMA
syssolver_linop_clover_quda_w.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /*! \file
3  * \brief Solve a MdagM*psi=chi linear system by BiCGStab
4  */
5 
6 #ifndef __syssolver_linop_quda_clover_h__
7 #define __syssolver_linop_quda_clover_h__
8 
9 #include "chroma_config.h"
10 
11 #ifdef BUILD_QUDA
12 #include <quda.h>
13 
14 #include "handle.h"
15 #include "state.h"
16 #include "syssolver.h"
17 #include "linearop.h"
18 #include "actions/ferm/fermbcs/simple_fermbc.h"
19 #include "actions/ferm/fermstates/periodic_fermstate.h"
20 #include "actions/ferm/invert/quda_solvers/syssolver_quda_clover_params.h"
21 #include "actions/ferm/linop/clover_term_w.h"
22 #include "meas/gfix/temporal_gauge.h"
23 #include "io/aniso_io.h"
24 #include <string>
25 
26 #include "util/gauge/reunit.h"
27 #ifdef QDP_IS_QDPJIT
28 #include "actions/ferm/invert/quda_solvers/qdpjit_memory_wrapper.h"
29 #endif
30 
31 //#include <util_quda.h>
32 
33 namespace Chroma
34 {
35 
36 //! Richardson system solver namespace
37 namespace LinOpSysSolverQUDACloverEnv
38 {
39 //! Register the syssolver
40 bool registerAll();
41 }
42 
43 
44 
45 //! Solve a Clover Fermion System using the QUDA inverter
46 /*! \ingroup invert
47  *** WARNING THIS SOLVER WORKS FOR Clover FERMIONS ONLY ***
48  */
49 
50 class LinOpSysSolverQUDAClover : public LinOpSystemSolver<LatticeFermion>
51 {
52 public:
53  typedef LatticeFermion T;
54  typedef LatticeColorMatrix U;
55  typedef multi1d<LatticeColorMatrix> Q;
56 
57  typedef LatticeFermionF TF;
58  typedef LatticeColorMatrixF UF;
59  typedef multi1d<LatticeColorMatrixF> QF;
60 
61  typedef LatticeFermionF TD;
62  typedef LatticeColorMatrixF UD;
63  typedef multi1d<LatticeColorMatrixF> QD;
64 
65  typedef WordType<T>::Type_t REALT;
66  //! Constructor
67  /*!
68  * \param M_ Linear operator ( Read )
69  * \param invParam inverter parameters ( Read )
70  */
71  LinOpSysSolverQUDAClover(Handle< LinearOperator<T> > A_,
72  Handle< FermState<T,Q,Q> > state_,
73  const SysSolverQUDACloverParams& invParam_) :
74  A(A_), invParam(invParam_), clov(new CloverTermT<T, U>() ), invclov(new CloverTermT<T, U>())
75  {
76  QDPIO::cout << "LinOpSysSolverQUDAClover:" << std::endl;
77 
78  // FOLLOWING INITIALIZATION in test QUDA program
79 
80  // 1) work out cpu_prec, cuda_prec, cuda_prec_sloppy
81  int s = sizeof( WordType<T>::Type_t );
82  if (s == 4) {
83  cpu_prec = QUDA_SINGLE_PRECISION;
84  }
85  else {
86  cpu_prec = QUDA_DOUBLE_PRECISION;
87  }
88 
89 
90  // Work out GPU precision
91  switch( invParam.cudaPrecision ) {
92  case HALF:
93  gpu_prec = QUDA_HALF_PRECISION;
94  break;
95  case SINGLE:
96  gpu_prec = QUDA_SINGLE_PRECISION;
97  break;
98  case DOUBLE:
99  gpu_prec = QUDA_DOUBLE_PRECISION;
100  break;
101  default:
102  gpu_prec = cpu_prec;
103  break;
104  }
105 
106  // Work out GPU Sloppy precision
107  // Default: No Sloppy
108  switch( invParam.cudaSloppyPrecision ) {
109  case HALF:
110  gpu_half_prec = QUDA_HALF_PRECISION;
111  break;
112  case SINGLE:
113  gpu_half_prec = QUDA_SINGLE_PRECISION;
114  break;
115  case DOUBLE:
116  gpu_half_prec = QUDA_DOUBLE_PRECISION;
117  break;
118  default:
119  gpu_half_prec = gpu_prec;
120  break;
121  }
122 
123  // 2) pull 'new; GAUGE and Invert params
124  q_gauge_param = newQudaGaugeParam();
125  quda_inv_param = newQudaInvertParam();
126 
127  // 3) set lattice size
128  const multi1d<int>& latdims = Layout::subgridLattSize();
129 
130  q_gauge_param.X[0] = latdims[0];
131  q_gauge_param.X[1] = latdims[1];
132  q_gauge_param.X[2] = latdims[2];
133  q_gauge_param.X[3] = latdims[3];
134 
135  // 4) - deferred (anisotropy)
136 
137  // 5) - set QUDA_WILSON_LINKS, QUDA_GAUGE_ORDER
138  q_gauge_param.type = QUDA_WILSON_LINKS;
139 #ifndef BUILD_QUDA_DEVIFACE_GAUGE
140  q_gauge_param.gauge_order = QUDA_QDP_GAUGE_ORDER; // gauge[mu], p
141 #else
142  QDPIO::cout << "MDAGM Using QDP-JIT gauge order" << std::endl;
143  q_gauge_param.location = QUDA_CUDA_FIELD_LOCATION;
144  q_gauge_param.gauge_order = QUDA_QDPJIT_GAUGE_ORDER;
145 #endif
146 
147  // 6) - set t_boundary
148  // Convention: BC has to be applied already
149  // This flag just tells QUDA that this is so,
150  // so that QUDA can take care in the reconstruct
151  if( invParam.AntiPeriodicT ) {
152  q_gauge_param.t_boundary = QUDA_ANTI_PERIODIC_T;
153  }
154  else {
155  q_gauge_param.t_boundary = QUDA_PERIODIC_T;
156  }
157 
158  // Set cpu_prec, cuda_prec, reconstruct and sloppy versions
159  q_gauge_param.cpu_prec = cpu_prec;
160  q_gauge_param.cuda_prec = gpu_prec;
161 
162 
163  switch( invParam.cudaReconstruct ) {
164  case RECONS_NONE:
165  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_NO;
166  break;
167  case RECONS_8:
168  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_8;
169  break;
170  case RECONS_12:
171  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_12;
172  break;
173  default:
174  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_12;
175  break;
176  };
177 
178  q_gauge_param.cuda_prec_sloppy = gpu_half_prec;
179 
180  // Default for no preconditioner -- may be overwritten based
181  // on innerParams
182  q_gauge_param.cuda_prec_precondition = gpu_half_prec;
183 
184  switch( invParam.cudaSloppyReconstruct ) {
185  case RECONS_NONE:
186  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_NO;
187  break;
188  case RECONS_8:
189  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_8;
190  break;
191  case RECONS_12:
192  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_12;
193  break;
194  default:
195  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_12;
196  break;
197  };
198 
199  // Default. This may be overrridden later.
200  q_gauge_param.reconstruct_precondition = q_gauge_param.reconstruct_sloppy;
201  // Gauge fixing:
202 
203  // These are the links
204  // They may be smeared and the BC's may be applied
205  Q links_single(Nd);
206 
207  // Now downcast to single prec fields.
208  for(int mu=0; mu < Nd; mu++) {
209  links_single[mu] = (state_->getLinks())[mu];
210  }
211 
212  // GaugeFix
213  if( invParam.axialGaugeP ) {
214  QDPIO::cout << "Fixing Temporal Gauge" << std::endl;
215  temporalGauge(links_single, GFixMat, Nd-1);
216  for(int mu=0; mu < Nd; mu++){
217  links_single[mu] = GFixMat*(state_->getLinks())[mu]*adj(shift(GFixMat, FORWARD, mu));
218  }
219  q_gauge_param.gauge_fix = QUDA_GAUGE_FIXED_YES;
220  }
221  else {
222  // No GaugeFix
223  q_gauge_param.gauge_fix = QUDA_GAUGE_FIXED_NO; // No Gfix yet
224  }
225 
226  // deferred 4) Gauge Anisotropy
227  const AnisoParam_t& aniso = invParam.CloverParams.anisoParam;
228  if( aniso.anisoP ) { // Anisotropic case
229  Real gamma_f = aniso.xi_0 / aniso.nu;
230  q_gauge_param.anisotropy = toDouble(gamma_f);
231  }
232  else {
233  q_gauge_param.anisotropy = 1.0;
234  }
235 
236  // MAKE FSTATE BEFORE RESCALING links_single
237  // Because the clover term expects the unrescaled links...
238  Handle<FermState<T,Q,Q> > fstate( new PeriodicFermState<T,Q,Q>(links_single));
239 
240  if( aniso.anisoP ) { // Anisotropic case
241  multi1d<Real> cf=makeFermCoeffs(aniso);
242  for(int mu=0; mu < Nd; mu++) {
243  links_single[mu] *= cf[mu];
244  }
245  }
246 
247  // Now onto the inv param:
248  // Dslash type
249  quda_inv_param.dslash_type = QUDA_CLOVER_WILSON_DSLASH;
250 
251  // Invert type:
252  switch( invParam.solverType ) {
253  case CG:
254  quda_inv_param.inv_type = QUDA_CG_INVERTER;
255  solver_string = "CG";
256  break;
257  case BICGSTAB:
258  quda_inv_param.inv_type = QUDA_BICGSTAB_INVERTER;
259  solver_string = "BICGSTAB";
260  break;
261  case GCR:
262  quda_inv_param.inv_type = QUDA_GCR_INVERTER;
263  solver_string = "GCR";
264  break;
265  default:
266  QDPIO::cerr << "Unknown SOlver type" << std::endl;
267  QDP_abort(1);
268  break;
269  }
270 
271  // Mass
272 
273  // Fiendish idea from Ron. Set the kappa=1/2 and use
274  // unmodified clover term, and ask for Kappa normalization
275  // This should give us A - (1/2)D as the unpreconditioned operator
276  // and probabl 1 - {1/4} A^{-1} D A^{-1} D as the preconditioned
277  // op. Apart from the A_oo stuff on the antisymmetric we have
278  // nothing to do...
279  quda_inv_param.kappa = 0.5;
280 
281  // FIXME: If we want QUDA to compute the clover coeff, we need to be able to deal
282  // with awfuless of anisotropy
283  // The value below is a dummy one.
284  quda_inv_param.clover_coeff = 1.0; // Dummy tree level value. Not used
285  quda_inv_param.tol = toDouble(invParam.RsdTarget);
286  quda_inv_param.maxiter = invParam.MaxIter;
287  quda_inv_param.reliable_delta = toDouble(invParam.Delta);
288  quda_inv_param.pipeline = invParam.Pipeline;
289 
290  // Solution type
291  quda_inv_param.solution_type = QUDA_MATPC_SOLUTION;
292 
293  // Solve type
294  switch( invParam.solverType ) {
295  case CG:
296  quda_inv_param.solve_type = QUDA_NORMOP_PC_SOLVE;
297  break;
298  case BICGSTAB:
299  quda_inv_param.solve_type = QUDA_DIRECT_PC_SOLVE;
300  break;
301  case GCR:
302  quda_inv_param.solve_type = QUDA_DIRECT_PC_SOLVE;
303  break;
304  case CA_GCR:
305  quda_inv_param.solve_type = QUDA_DIRECT_PC_SOLVE;
306  break;
307  case MR:
308  quda_inv_param.solve_type = QUDA_DIRECT_PC_SOLVE;
309  break;
310 
311  default:
312  quda_inv_param.solve_type = QUDA_NORMOP_PC_SOLVE;
313 
314  break;
315  }
316 
317  if( invParam.asymmetricP ) {
318  QDPIO::cout << "Using Asymmetric Linop: A_oo - D A^{-1}_ee D" << std::endl;
319  quda_inv_param.matpc_type = QUDA_MATPC_ODD_ODD_ASYMMETRIC;
320  }
321  else {
322  QDPIO::cout << "Using Symmetric Linop: 1 - A^{-1}_oo D A^{-1}_ee D" << std::endl;
323  quda_inv_param.matpc_type = QUDA_MATPC_ODD_ODD;
324  }
325 
326  quda_inv_param.dagger = QUDA_DAG_NO;
327  quda_inv_param.mass_normalization = QUDA_KAPPA_NORMALIZATION;
328 
329  quda_inv_param.cpu_prec = cpu_prec;
330  quda_inv_param.cuda_prec = gpu_prec;
331  quda_inv_param.cuda_prec_sloppy = gpu_half_prec;
332 
333  // Default. May be overridden by inner params
334  quda_inv_param.cuda_prec_precondition = gpu_half_prec;
335 
336 
337  quda_inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
338  quda_inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
339 
340 #ifndef BUILD_QUDA_DEVIFACE_SPINOR
341  quda_inv_param.dirac_order = QUDA_DIRAC_ORDER;
342 #else
343  QDPIO::cout << "MDAGM Using QDP-JIT spinor order" << std::endl;
344  quda_inv_param.dirac_order = QUDA_QDPJIT_DIRAC_ORDER;
345  quda_inv_param.input_location = QUDA_CUDA_FIELD_LOCATION;
346  quda_inv_param.output_location = QUDA_CUDA_FIELD_LOCATION;
347 #endif
348 
349 
350  // Clover precision and order
351  quda_inv_param.clover_cpu_prec = cpu_prec;
352  quda_inv_param.clover_cuda_prec = gpu_prec;
353  quda_inv_param.clover_cuda_prec_sloppy = gpu_half_prec;
354 
355  // Default. may be overrridden by inner params
356  quda_inv_param.clover_cuda_prec_precondition = gpu_half_prec;
357 
358 #ifndef BUILD_QUDA_DEVIFACE_CLOVER
359 #warning "NOT USING QUDA DEVICE IFACE"
360  quda_inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
361 #else
362 #warning "USING QUDA DEVICE IFACE"
363  QDPIO::cout << "MDAGM clover CUDA location\n";
364  quda_inv_param.clover_location = QUDA_CUDA_FIELD_LOCATION;
365  quda_inv_param.clover_order = QUDA_QDPJIT_CLOVER_ORDER;
366 #endif
367 
368  // Autotuning
369  if( invParam.tuneDslashP ) {
370  QDPIO::cout << "Enabling Dslash Autotuning" << std::endl;
371 
372  quda_inv_param.tune = QUDA_TUNE_YES;
373  }
374  else {
375  QDPIO::cout << "Disabling Dslash Autotuning" << std::endl;
376 
377  quda_inv_param.tune = QUDA_TUNE_NO;
378  }
379 
380 
381  // Setup padding
382  multi1d<int> face_size(4);
383  face_size[0] = latdims[1]*latdims[2]*latdims[3]/2;
384  face_size[1] = latdims[0]*latdims[2]*latdims[3]/2;
385  face_size[2] = latdims[0]*latdims[1]*latdims[3]/2;
386  face_size[3] = latdims[0]*latdims[1]*latdims[2]/2;
387 
388  int max_face = face_size[0];
389  for(int i=1; i <=3; i++) {
390  if ( face_size[i] > max_face ) {
391  max_face = face_size[i];
392  }
393  }
394 
395 
396  q_gauge_param.ga_pad = max_face;
397  // PADDING
398  quda_inv_param.sp_pad = 0;
399  quda_inv_param.cl_pad = 0;
400 
401  if( invParam.innerParamsP ) {
402  QDPIO::cout << "Setting inner solver params" << std::endl;
403  // Dereference handle
404  const GCRInnerSolverParams& ip = *(invParam.innerParams);
405 
406  // Set preconditioner precision
407  switch( ip.precPrecondition ) {
408  case HALF:
409  quda_inv_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
410  quda_inv_param.clover_cuda_prec_precondition = QUDA_HALF_PRECISION;
411  q_gauge_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
412  break;
413 
414  case SINGLE:
415  quda_inv_param.cuda_prec_precondition = QUDA_SINGLE_PRECISION;
416  quda_inv_param.clover_cuda_prec_precondition = QUDA_SINGLE_PRECISION;
417  q_gauge_param.cuda_prec_precondition = QUDA_SINGLE_PRECISION;
418  break;
419 
420  case DOUBLE:
421  quda_inv_param.cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
422  quda_inv_param.clover_cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
423  q_gauge_param.cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
424  break;
425  default:
426  quda_inv_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
427  quda_inv_param.clover_cuda_prec_precondition = QUDA_HALF_PRECISION;
428  q_gauge_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
429  break;
430  }
431 
432  switch( ip.reconstructPrecondition ) {
433  case RECONS_NONE:
434  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_NO;
435  break;
436  case RECONS_8:
437  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_8;
438  break;
439  case RECONS_12:
440  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_12;
441  break;
442  default:
443  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_12;
444  break;
445  };
446 
447  quda_inv_param.tol_precondition = toDouble(ip.tolPrecondition);
448  quda_inv_param.maxiter_precondition = ip.maxIterPrecondition;
449  quda_inv_param.gcrNkrylov = ip.gcrNkrylov;
450  switch( ip.schwarzType ) {
451  case ADDITIVE_SCHWARZ :
452  quda_inv_param.schwarz_type = QUDA_ADDITIVE_SCHWARZ;
453  break;
454  case MULTIPLICATIVE_SCHWARZ :
455  quda_inv_param.schwarz_type = QUDA_MULTIPLICATIVE_SCHWARZ;
456  break;
457  default:
458  quda_inv_param.schwarz_type = QUDA_ADDITIVE_SCHWARZ;
459  break;
460  }
461  quda_inv_param.precondition_cycle = ip.preconditionCycle;
462 
463  if( ip.verboseInner ) {
464  quda_inv_param.verbosity_precondition = QUDA_VERBOSE;
465  }
466  else {
467  quda_inv_param.verbosity_precondition = QUDA_SILENT;
468  }
469 
470  switch( ip.invTypePrecondition ) {
471  case CG:
472  quda_inv_param.inv_type_precondition = QUDA_CG_INVERTER;
473  break;
474  case BICGSTAB:
475  quda_inv_param.inv_type_precondition = QUDA_BICGSTAB_INVERTER;
476 
477  break;
478 
479  case CA_GCR:
480  quda_inv_param.inv_type_precondition= QUDA_CA_GCR_INVERTER;
481 
482  case MR:
483  quda_inv_param.inv_type_precondition= QUDA_MR_INVERTER;
484  break;
485 
486  default:
487  quda_inv_param.inv_type_precondition = QUDA_MR_INVERTER;
488  break;
489  }
490  }
491  else {
492  QDPIO::cout << "Setting Precondition stuff to defaults for not using" << std::endl;
493  quda_inv_param.inv_type_precondition= QUDA_INVALID_INVERTER;
494  quda_inv_param.tol_precondition = 1.0e-1;
495  quda_inv_param.maxiter_precondition = 1000;
496  quda_inv_param.verbosity_precondition = QUDA_SILENT;
497  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_NO;
498  quda_inv_param.gcrNkrylov = 1;
499  }
500 
501 
502  if( invParam.verboseP ) {
503  quda_inv_param.verbosity = QUDA_VERBOSE;
504  }
505  else {
506  quda_inv_param.verbosity = QUDA_SUMMARIZE;
507  }
508 
509  // Set up the links
510  void* gauge[4];
511 
512 #ifndef BUILD_QUDA_DEVIFACE_GAUGE
513  for(int mu=0; mu < Nd; mu++) {
514  gauge[mu] = (void *)&(links_single[mu].elem(all.start()).elem().elem(0,0).real());
515  }
516 #else
517  GetMemoryPtrGauge(gauge,links_single);
518  // gauge[mu] = GetMemoryPtr( links_single[mu].getId() );
519  // QDPIO::cout << "MDAGM CUDA gauge[" << mu << "] in = " << gauge[mu] << "\n";
520 #endif
521 
522  loadGaugeQuda((void *)gauge, &q_gauge_param);
523 
524  // Setup the clover term...
525  QDPIO::cout << "Creating CloverTerm" << std::endl;
526  clov->create(fstate, invParam_.CloverParams);
527  // Don't recompute, just copy
528  invclov->create(fstate, invParam_.CloverParams);
529 
530  QDPIO::cout << "Inverting CloverTerm" << std::endl;
531  invclov->choles(0);
532  invclov->choles(1);
533 
534 
535 #ifndef BUILD_QUDA_DEVIFACE_CLOVER
536  multi1d<QUDAPackedClovSite<REALT> > packed_clov;
537 
538 
539  packed_clov.resize(all.siteTable().size());
540 
541  clov->packForQUDA(packed_clov, 0);
542  clov->packForQUDA(packed_clov, 1);
543 
544 
545  // Always need inverse
546  multi1d<QUDAPackedClovSite<REALT> > packed_invclov(all.siteTable().size());
547  invclov->packForQUDA(packed_invclov, 0);
548  invclov->packForQUDA(packed_invclov, 1);
549 
550 
551  loadCloverQuda(&(packed_clov[0]), &(packed_invclov[0]),&quda_inv_param);
552 
553 #else
554  void *clover[2];
555  void *cloverInv[2];
556 
557  GetMemoryPtrClover(clov->getOffId(),clov->getDiaId(),invclov->getOffId(),invclov->getDiaId());
558 
559  loadCloverQuda( (void*)(clover) , (void*)(cloverInv) ,&quda_inv_param);
560 
561 #endif
562 
563 
564 
565  }
566 
567 
568  //! Destructor is automatic
569  ~LinOpSysSolverQUDAClover()
570  {
571  QDPIO::cout << "Destructing" << std::endl;
572  freeGaugeQuda();
573  freeCloverQuda();
574  }
575 
576  //! Return the subset on which the operator acts
577  const Subset& subset() const {return A->subset();}
578 
579  //! Solver the linear system
580  /*!
581  * \param psi solution ( Modify )
582  * \param chi source ( Read )
583  * \return syssolver results
584  */
585  SystemSolverResults_t operator() (T& psi, const T& chi) const
586  {
587  SystemSolverResults_t res;
588 
589  START_CODE();
590  StopWatch swatch;
591  swatch.start();
592 
593  // T MdagChi;
594 
595  // This is a CGNE. So create new RHS
596  // (*A)(MdagChi, chi, MINUS);
597  // Handle< LinearOperator<T> > MM(new MdagMLinOp<T>(A));
598  if ( invParam.axialGaugeP ) {
599  T g_chi,g_psi;
600 
601  // Gauge Fix source and initial guess
602  QDPIO::cout << "Gauge Fixing source and initial guess" << std::endl;
603  g_chi[ rb[1] ] = GFixMat * chi;
604  g_psi[ rb[1] ] = GFixMat * psi;
605  QDPIO::cout << "Solving" << std::endl;
606  res = qudaInvert(*clov,
607  *invclov,
608  g_chi,
609  g_psi);
610  QDPIO::cout << "Untransforming solution." << std::endl;
611  psi[ rb[1]] = adj(GFixMat)*g_psi;
612 
613  }
614  else {
615  res = qudaInvert(*clov,
616  *invclov,
617  chi,
618  psi);
619  }
620 
621  swatch.stop();
622 
623 
624  {
625  T r;
626  r[A->subset()]=chi;
627  T tmp;
628  (*A)(tmp, psi, PLUS);
629  r[A->subset()] -= tmp;
630  res.resid = sqrt(norm2(r, A->subset()));
631  }
632 
633  Double rel_resid = res.resid/sqrt(norm2(chi,A->subset()));
634 
635  QDPIO::cout << "QUDA_"<< solver_string <<"_CLOVER_SOLVER: " << res.n_count << " iterations. Rsd = " << res.resid << " Relative Rsd = " << rel_resid << std::endl;
636 
637  // Convergence Check/Blow Up
638  if ( ! invParam.SilentFailP ) {
639  if ( toBool( rel_resid > invParam.RsdToleranceFactor*invParam.RsdTarget) ) {
640  QDPIO::cerr << "ERROR: QUDA Solver residuum is outside tolerance: QUDA resid="<< rel_resid << " Desired =" << invParam.RsdTarget << " Max Tolerated = " << invParam.RsdToleranceFactor*invParam.RsdTarget << std::endl;
641  QDP_abort(1);
642  }
643  }
644 
645  END_CODE();
646  return res;
647  }
648 
649 
650 private:
651  // Hide default constructor
652  LinOpSysSolverQUDAClover() {}
653 
654 #if 1
655  Q links_orig;
656 #endif
657 
658  U GFixMat;
659  QudaPrecision_s cpu_prec;
660  QudaPrecision_s gpu_prec;
661  QudaPrecision_s gpu_half_prec;
662 
663  Handle< LinearOperator<T> > A;
664  const SysSolverQUDACloverParams invParam;
665  QudaGaugeParam q_gauge_param;
666  QudaInvertParam quda_inv_param;
667 
668  Handle< CloverTermT<T, U> > clov;
669  Handle< CloverTermT<T, U> > invclov;
670 
671  SystemSolverResults_t qudaInvert(const CloverTermT<T, U>& clover,
672  const CloverTermT<T, U>& inv_clov,
673  const T& chi_s,
674  T& psi_s
675  )const ;
676 
677  std::string solver_string;
678 };
679 
680 
681 } // End namespace
682 
683 #endif // BUILD_QUDA
684 #endif
685 
aniso_io.h
Anisotropy parameters.
chroma_config.h
clover_term_w.h
Include possibly optimized Clover terms.
mu
int mu
Definition: cool.cc:24
Chroma::temporalGauge
void temporalGauge(multi1d< LatticeColorMatrix > &ug, LatticeColorMatrix &g, int decay_dir)
Temporal gauge fixing.
Definition: temporal_gauge.cc:26
Chroma::MR
@ MR
Definition: enum_quda_io.h:26
Chroma::GCR
@ GCR
Definition: enum_quda_io.h:24
Chroma::CG
@ CG
Definition: syssolver_qphix_clover_params.h:16
Chroma::CA_GCR
@ CA_GCR
Definition: enum_quda_io.h:25
Chroma::BICGSTAB
@ BICGSTAB
Definition: syssolver_qphix_clover_params.h:17
handle.h
Class for counted reference semantics.
linearop.h
Linear Operators.
Nd
Nd
Definition: meslate.cc:74
Chroma::LinOpSysSolverQUDACloverEnv::registerAll
bool registerAll()
Register all the factories.
Definition: syssolver_linop_clover_quda_w.cc:47
Chroma::SimpleBaryonSeqSourceEnv::operator()
multi1d< Hadron2PtContraction_t > operator()(const multi1d< LatticeColorMatrix > &u)
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::Q
LinOpSysSolverMGProtoClover::Q Q
Definition: syssolver_linop_clover_mg_proto.cc:64
Chroma::tmp
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
Chroma::T
LinOpSysSolverMGProtoClover::T T
Definition: syssolver_linop_clover_mg_proto.cc:63
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::PLUS
@ PLUS
Definition: chromabase.h:45
Chroma::chi
multi1d< LatticeFermion > chi(Ncb)
Chroma::psi
LatticeFermion psi
Definition: mespbg5p_w.cc:35
Chroma::CloverTermT
QDPCloverTermT< T, U > CloverTermT
Definition: clover_term_w.h:97
Chroma::START_CODE
START_CODE()
Chroma::A
A(A, psi, r, Ncb, PLUS)
Chroma::makeFermCoeffs
multi1d< Real > makeFermCoeffs(const AnisoParam_t &aniso)
Make fermion coefficients.
Definition: aniso_io.cc:63
Chroma::s
multi1d< LatticeFermion > s(Ncb)
Chroma::ADDITIVE_SCHWARZ
@ ADDITIVE_SCHWARZ
Definition: enum_quda_io.h:103
Chroma::MULTIPLICATIVE_SCHWARZ
@ MULTIPLICATIVE_SCHWARZ
Definition: enum_quda_io.h:104
Chroma::RECONS_12
@ RECONS_12
Definition: enum_quda_io.h:80
Chroma::RECONS_NONE
@ RECONS_NONE
Definition: enum_quda_io.h:78
Chroma::RECONS_8
@ RECONS_8
Definition: enum_quda_io.h:79
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
testing::internal::string
::std::string string
Definition: gtest.h:1979
periodic_fermstate.h
Periodic ferm state and a creator.
FORWARD
#define FORWARD
Definition: primitives.h:82
qdpjit_memory_wrapper.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_quda_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
QF
multi1d< LatticeColorMatrixF > QF
Definition: t_quda_tprec.cc:19
UF
LatticeColorMatrixF UF
Definition: t_quda_tprec.cc:18
TD
LatticeFermionD TD
Definition: t_quda_tprec.cc:22
UD
LatticeColorMatrixD UD
Definition: t_quda_tprec.cc:23
TF
LatticeFermionF TF
Definition: t_quda_tprec.cc:17
QD
multi1d< LatticeColorMatrixD > QD
Definition: t_quda_tprec.cc:24
temporal_gauge.h
Axial gauge fixing.
Generated by doxygen 1.9.1