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