CHROMA
syssolver_linop_wilson_quda_multigrid_w.h
Go to the documentation of this file.
1 
2 // -*- C++ -*-
3 /*! \file
4  * \QUDA MULTIGRID Wilson solver.
5  */
6 
7 #ifndef __syssolver_linop_quda_multigrid_wilson_h__
8 #define __syssolver_linop_quda_multigrid_wilson_h__
9 
10 #include "chroma_config.h"
11 
12 #ifdef BUILD_QUDA
13 #include <quda.h>
14 
15 #include "handle.h"
16 #include "state.h"
17 #include "syssolver.h"
18 #include "linearop.h"
19 #include "actions/ferm/fermbcs/simple_fermbc.h"
20 #include "actions/ferm/fermstates/periodic_fermstate.h"
21 #include "actions/ferm/invert/quda_solvers/syssolver_quda_multigrid_wilson_params.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 LinOpSysSolverQUDAMULTIGRIDWilsonEnv
38 {
39 //! Register the syssolver
40 bool registerAll();
41 }
42 
43 
44 
45 //! Solve a Wilson Fermion System using the QUDA inverter
46 /*! \ingroup invert
47  *** WARNING THIS SOLVER WORKS FOR Wilson FERMIONS ONLY ***
48  */
49 
50 class LinOpSysSolverQUDAMULTIGRIDWilson : 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  LinOpSysSolverQUDAMULTIGRIDWilson(Handle< LinearOperator<T> > A_,
72  Handle< FermState<T,Q,Q> > state_,
73  const SysSolverQUDAMULTIGRIDWilsonParams& invParam_) :
74  A(A_), invParam(invParam_)
75  {
76  QDPIO::cout << "LinOpSysSolverQUDAMULTIGRIDWilson:" << 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  mg_inv_param = newQudaInvertParam();
127  mg_param = newQudaMultigridParam();
128 
129 
130  // 3) set lattice size
131  const multi1d<int>& latdims = Layout::subgridLattSize();
132 
133  q_gauge_param.X[0] = latdims[0];
134  q_gauge_param.X[1] = latdims[1];
135  q_gauge_param.X[2] = latdims[2];
136  q_gauge_param.X[3] = latdims[3];
137 
138  // 4) - deferred (anisotropy)
139 
140  // 5) - set QUDA_WILSON_LINKS, QUDA_GAUGE_ORDER
141  q_gauge_param.type = QUDA_WILSON_LINKS;
142 #ifndef BUILD_QUDA_DEVIFACE_GAUGE
143  q_gauge_param.gauge_order = QUDA_QDP_GAUGE_ORDER; // gauge[mu], p
144 #else
145  QDPIO::cout << "MDAGM Using QDP-JIT gauge order" << std::endl;
146  q_gauge_param.location = QUDA_CUDA_FIELD_LOCATION;
147  q_gauge_param.gauge_order = QUDA_QDPJIT_GAUGE_ORDER;
148 #endif
149 
150  // 6) - set t_boundary
151  // Convention: BC has to be applied already
152  // This flag just tells QUDA that this is so,
153  // so that QUDA can take care in the reconstruct
154  if( invParam.AntiPeriodicT ) {
155  q_gauge_param.t_boundary = QUDA_ANTI_PERIODIC_T;
156  }
157  else {
158  q_gauge_param.t_boundary = QUDA_PERIODIC_T;
159  }
160 
161  // Set cpu_prec, cuda_prec, reconstruct and sloppy versions
162  q_gauge_param.cpu_prec = cpu_prec;
163  q_gauge_param.cuda_prec = gpu_prec;
164 
165 
166  switch( invParam.cudaReconstruct ) {
167  case RECONS_NONE:
168  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_NO;
169  break;
170  case RECONS_8:
171  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_8;
172  break;
173  case RECONS_12:
174  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_12;
175  break;
176  default:
177  q_gauge_param.reconstruct = QUDA_RECONSTRUCT_12;
178  break;
179  };
180 
181  q_gauge_param.cuda_prec_sloppy = gpu_half_prec;
182 
183  switch( invParam.cudaSloppyReconstruct ) {
184  case RECONS_NONE:
185  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_NO;
186  break;
187  case RECONS_8:
188  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_8;
189  break;
190  case RECONS_12:
191  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_12;
192  break;
193  default:
194  q_gauge_param.reconstruct_sloppy = QUDA_RECONSTRUCT_12;
195  break;
196  };
197 
198  // Gauge fixing:
199 
200  // These are the links
201  // They may be smeared and the BC's may be applied
202  Q links_single(Nd);
203 
204  // Now downcast to single prec fields.
205  for(int mu=0; mu < Nd; mu++) {
206  links_single[mu] = (state_->getLinks())[mu];
207  }
208 
209  // GaugeFix
210  if( invParam.axialGaugeP ) {
211  QDPIO::cout << "Fixing Temporal Gauge" << std::endl;
212  temporalGauge(links_single, GFixMat, Nd-1);
213  for(int mu=0; mu < Nd; mu++){
214  links_single[mu] = GFixMat*(state_->getLinks())[mu]*adj(shift(GFixMat, FORWARD, mu));
215  }
216  q_gauge_param.gauge_fix = QUDA_GAUGE_FIXED_YES;
217  }
218  else {
219  // No GaugeFix
220  q_gauge_param.gauge_fix = QUDA_GAUGE_FIXED_NO; // No Gfix yet
221  }
222 
223  // deferred 4) Gauge Anisotropy
224  const AnisoParam_t& aniso = invParam.WilsonParams.anisoParam;
225  if( aniso.anisoP ) { // Anisotropic case
226  Real gamma_f = aniso.xi_0 / aniso.nu;
227  q_gauge_param.anisotropy = toDouble(gamma_f);
228  }
229  else {
230  q_gauge_param.anisotropy = 1.0;
231  }
232 
233  // MAKE FSTATE BEFORE RESCALING links_single
234  // Because the clover term expects the unrescaled links...
235  Handle<FermState<T,Q,Q> > fstate( new PeriodicFermState<T,Q,Q>(links_single));
236 
237  if( aniso.anisoP ) { // Anisotropic case
238  multi1d<Real> cf=makeFermCoeffs(aniso);
239  for(int mu=0; mu < Nd; mu++) {
240  links_single[mu] *= cf[mu];
241  }
242  }
243 
244  // Now onto the inv param:
245  // Dslash type
246 
247  /****!!! FIXME: Before the final code remember to reset this to QUDA_CLOVER_WILSON_DSLASH */
248  QDPIO::cout << "Remember for production to reset quda_inv_param.dslash_typeto QUDA_CLOVER_WILSON_DSLASH" << std::endl;
249  quda_inv_param.dslash_type = QUDA_WILSON_DSLASH;
250  mg_inv_param.dslash_type = QUDA_WILSON_DSLASH;
251 
252  // Invert type:
253  switch( invParam.solverType ) {
254  case CG:
255  quda_inv_param.inv_type = QUDA_CG_INVERTER;
256  solver_string = "CG";
257  break;
258  case BICGSTAB:
259  quda_inv_param.inv_type = QUDA_BICGSTAB_INVERTER;
260  solver_string = "BICGSTAB";
261  break;
262  case GCR:
263  quda_inv_param.inv_type = QUDA_GCR_INVERTER;
264  solver_string = "GCR";
265  break;
266  default:
267  QDPIO::cerr << "Unknown Solver type" << std::endl;
268  QDP_abort(1);
269  break;
270  }
271  //Params added for now to get this to initialize.
272  mg_inv_param.inv_type = QUDA_GCR_INVERTER;
273  mg_inv_param.tol = 1e-10;
274  mg_inv_param.maxiter = 10000;
275  mg_inv_param.reliable_delta = 1e-10;
276  mg_inv_param.verbosity = QUDA_VERBOSE;
277  mg_inv_param.verbosity_precondition = QUDA_VERBOSE;
278 
279 
280 
281  Real diag_mass;
282  {
283  // auto is C++11 so I don't have to remember all the silly typenames
284  auto wlparams = invParam.WilsonParams;
285 
286  auto aniso = wlparams.anisoParam;
287 
288  Real ff = where(aniso.anisoP, aniso.nu / aniso.xi_0, Real(1));
289  diag_mass = 1 + (Nd-1)*ff + wlparams.Mass;
290  }
291 
292 
293  quda_inv_param.kappa = static_cast<double>(1)/(static_cast<double>(2)*toDouble(diag_mass));
294  /**** END FIXME XXX ***/
295 
296  quda_inv_param.tol = toDouble(invParam.RsdTarget);
297  quda_inv_param.maxiter = invParam.MaxIter;
298  quda_inv_param.reliable_delta = toDouble(invParam.Delta);
299  quda_inv_param.pipeline = invParam.Pipeline;
300 
301  // Solution type
302  //quda_inv_param.solution_type = QUDA_MATPC_SOLUTION;
303  //Taken from invert test.
304  quda_inv_param.solution_type = QUDA_MATPC_SOLUTION;
305  quda_inv_param.solve_type = QUDA_DIRECT_PC_SOLVE;
306 
307  QDPIO::cout << "Using Symmetric Linop: 1 - A^{-1}_oo D A^{-1}_ee D" << std::endl;
308  quda_inv_param.matpc_type = QUDA_MATPC_ODD_ODD;
309 
310  quda_inv_param.dagger = QUDA_DAG_NO;
311  quda_inv_param.mass_normalization = QUDA_KAPPA_NORMALIZATION;
312 
313  quda_inv_param.cpu_prec = cpu_prec;
314  quda_inv_param.cuda_prec = gpu_prec;
315  quda_inv_param.cuda_prec_sloppy = gpu_half_prec;
316  //Add some lines for mg_inv_param.
317  mg_inv_param.cpu_prec = cpu_prec;
318  mg_inv_param.cuda_prec = gpu_prec;
319  mg_inv_param.cuda_prec_sloppy = gpu_half_prec;
320 
321  //Clover stuff
322  // This doesn't hurt to leave in, to make sure these are in a well defined state
323  // But yucky!!!
324  mg_inv_param.clover_cpu_prec = cpu_prec;
325  mg_inv_param.clover_cuda_prec = gpu_prec;
326  mg_inv_param.clover_cuda_prec_sloppy = gpu_half_prec;
327  mg_inv_param.clover_cuda_prec_precondition = gpu_prec;
328  mg_inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
329  //
330  //Done...
331  quda_inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
332  quda_inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
333 
334 #ifndef BUILD_QUDA_DEVIFACE_SPINOR
335  quda_inv_param.dirac_order = QUDA_DIRAC_ORDER;
336 #else
337  QDPIO::cout << "MDAGM Using QDP-JIT spinor order" << std::endl;
338  quda_inv_param.dirac_order = QUDA_QDPJIT_DIRAC_ORDER;
339  quda_inv_param.input_location = QUDA_CUDA_FIELD_LOCATION;
340  quda_inv_param.output_location = QUDA_CUDA_FIELD_LOCATION;
341 #endif
342 
343  // Autotuning
344  if( invParam.tuneDslashP ) {
345  QDPIO::cout << "Enabling Dslash Autotuning" << std::endl;
346 
347  quda_inv_param.tune = QUDA_TUNE_YES;
348  mg_inv_param.tune = QUDA_TUNE_YES;
349  }
350  else {
351  QDPIO::cout << "Disabling Dslash Autotuning" << std::endl;
352 
353  quda_inv_param.tune = QUDA_TUNE_NO;
354  mg_inv_param.tune = QUDA_TUNE_NO;
355  }
356 
357 
358  // Setup padding
359  multi1d<int> face_size(4);
360  face_size[0] = latdims[1]*latdims[2]*latdims[3]/2;
361  face_size[1] = latdims[0]*latdims[2]*latdims[3]/2;
362  face_size[2] = latdims[0]*latdims[1]*latdims[3]/2;
363  face_size[3] = latdims[0]*latdims[1]*latdims[2]/2;
364 
365  int max_face = face_size[0];
366  for(int i=1; i <=3; i++) {
367  if ( face_size[i] > max_face ) {
368  max_face = face_size[i];
369  }
370  }
371 
372 
373  q_gauge_param.ga_pad = max_face;
374  // PADDING
375  quda_inv_param.sp_pad = 0;
376  quda_inv_param.cl_pad = 0;
377 
378  if( invParam.MULTIGRIDParamsP ) {
379  QDPIO::cout << "Setting MULTIGRID solver params" << std::endl;
380  // Dereference handle
381  const MULTIGRIDSolverParams& ip = *(invParam.MULTIGRIDParams);
382 
383  // Set preconditioner precision
384  switch( ip.prec ) {
385  case HALF:
386  mg_inv_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
387  quda_inv_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
388  mg_inv_param.clover_cuda_prec_precondition = QUDA_HALF_PRECISION;
389  quda_inv_param.clover_cuda_prec_precondition = QUDA_HALF_PRECISION;
390  q_gauge_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
391  break;
392 
393  case SINGLE:
394  mg_inv_param.cuda_prec_precondition = QUDA_SINGLE_PRECISION;
395  mg_inv_param.clover_cuda_prec_precondition = QUDA_SINGLE_PRECISION;
396  quda_inv_param.cuda_prec_precondition = QUDA_SINGLE_PRECISION;
397  quda_inv_param.clover_cuda_prec_precondition = QUDA_SINGLE_PRECISION;
398  q_gauge_param.cuda_prec_precondition = QUDA_SINGLE_PRECISION;
399  break;
400 
401  case DOUBLE:
402  mg_inv_param.cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
403  mg_inv_param.clover_cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
404  quda_inv_param.cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
405  quda_inv_param.clover_cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
406  q_gauge_param.cuda_prec_precondition = QUDA_DOUBLE_PRECISION;
407  break;
408  default:
409  mg_inv_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
410  mg_inv_param.clover_cuda_prec_precondition = QUDA_HALF_PRECISION;
411  quda_inv_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
412  quda_inv_param.clover_cuda_prec_precondition = QUDA_HALF_PRECISION;
413  q_gauge_param.cuda_prec_precondition = QUDA_HALF_PRECISION;
414  break;
415  }
416 
417  switch( ip.reconstruct ) {
418  case RECONS_NONE:
419  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_NO;
420  break;
421  case RECONS_8:
422  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_8;
423  break;
424  case RECONS_12:
425  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_12;
426  break;
427  default:
428  q_gauge_param.reconstruct_precondition = QUDA_RECONSTRUCT_12;
429  break;
430  };
431  }
432  // Set up the links
433  void* gauge[4];
434 
435 #ifndef BUILD_QUDA_DEVIFACE_GAUGE
436  for(int mu=0; mu < Nd; mu++) {
437  gauge[mu] = (void *)&(links_single[mu].elem(all.start()).elem().elem(0,0).real());
438  }
439 #else
440  GetMemoryPtrGauge(gauge,links_single);
441  //gauge[mu] = GetMemoryPtr( links_single[mu].getId() );
442  //QDPIO::cout << "MDAGM CUDA gauge[" << mu << "] in = " << gauge[mu] << "\n";
443 #endif
444 
445  loadGaugeQuda((void *)gauge, &q_gauge_param);
446 
447 
448  MULTIGRIDSolverParams ip = *(invParam.MULTIGRIDParams);
449  //
450  quda_inv_param.tol_precondition = toDouble(ip.tol[0]);
451  quda_inv_param.maxiter_precondition = ip.maxIterations[0];
452  quda_inv_param.gcrNkrylov = ip.outer_gcr_nkrylov;
453  mg_inv_param.gcrNkrylov = ip.precond_gcr_nkrylov;
454  //Replacing above with what's in the invert test.
455  switch( ip.schwarzType ) {
456  case ADDITIVE_SCHWARZ :
457  quda_inv_param.schwarz_type = QUDA_ADDITIVE_SCHWARZ;
458  break;
459  case MULTIPLICATIVE_SCHWARZ :
460  quda_inv_param.schwarz_type = QUDA_MULTIPLICATIVE_SCHWARZ;
461  break;
462  default:
463  quda_inv_param.schwarz_type = QUDA_ADDITIVE_SCHWARZ;
464  break;
465  }
466  quda_inv_param.precondition_cycle = 1;
467  //Invert test always sets this to 1.
468 
469  quda_inv_param.verbosity_precondition = QUDA_VERBOSE;
470 
471  //MG is the only option.
472  quda_inv_param.inv_type_precondition = QUDA_MG_INVERTER;
473  //New invert test changes here.
474 
475  mg_inv_param.sp_pad = 0;
476  mg_inv_param.cl_pad = 0;
477 
478  mg_inv_param.preserve_source = QUDA_PRESERVE_SOURCE_NO;
479  mg_inv_param.gamma_basis = QUDA_DEGRAND_ROSSI_GAMMA_BASIS;
480  mg_inv_param.dirac_order = QUDA_DIRAC_ORDER;
481 
482  mg_inv_param.input_location = QUDA_CPU_FIELD_LOCATION;
483  mg_inv_param.output_location = QUDA_CPU_FIELD_LOCATION;
484 
485  mg_inv_param.kappa = quda_inv_param.kappa;
486 
487  mg_inv_param.dagger = QUDA_DAG_NO;
488  mg_inv_param.mass_normalization = QUDA_KAPPA_NORMALIZATION;
489 
490  mg_inv_param.matpc_type = QUDA_MATPC_ODD_ODD;
491  mg_inv_param.solution_type = QUDA_MAT_SOLUTION;
492  mg_inv_param.solve_type = QUDA_DIRECT_SOLVE;
493 
494  mg_param.invert_param = &mg_inv_param;
495 
496  mg_inv_param.Ls = 1;
497  quda_inv_param.Ls = 1;
498 
499  // FIXME: Make this an XML Param
500  mg_param.run_verify = QUDA_BOOLEAN_NO;
501 
502  mg_param.n_level = ip.mg_levels;
503 
504  for (int i=0; i<mg_param.n_level-1; ++i) {
505  if( ip.setup_on_gpu[i] ) {
506  mg_param.setup_location[i] = QUDA_CUDA_FIELD_LOCATION;
507  }
508  else {
509  mg_param.setup_location[i] = QUDA_CPU_FIELD_LOCATION;
510  }
511 
512  }
513 
514  for (int i=0; i<mg_param.n_level; i++) {
515  for (int j=0; j< Nd; j++) {
516  if( i < mg_param.n_level-1 ) {
517  mg_param.geo_block_size[i][j] = ip.blocking[i][j];
518  }
519  else {
520  mg_param.geo_block_size[i][j] = 4;
521  }
522  }
523  mg_param.spin_block_size[i] = 1;
524  if( i < mg_param.n_level-1) {
525  mg_param.n_vec[i] = ip.nvec[i];
526  mg_param.nu_pre[i] = ip.nu_pre[i];
527  mg_param.nu_post[i] = ip.nu_post[i];
528  }
529  mg_param.smoother_tol[i] = toDouble(ip.smootherTol[i]);
530  mg_param.global_reduction[i] = QUDA_BOOLEAN_YES;
531 
532  switch( ip.smootherType[i] ) {
533  case MR:
534  mg_param.smoother[i] = QUDA_MR_INVERTER;
535  mg_param.omega[i] = 0.85;
536  break;
537  case CA_GCR:
538  mg_param.smoother[i] = QUDA_CA_GCR_INVERTER;
539  mg_param.omega[i] = 0.85;
540  break;
541  default:
542  QDPIO::cout << "Unknown or no smother type specified, no smoothing inverter will be used." << std::endl;
543  mg_param.smoother[i] = QUDA_INVALID_INVERTER;
544  break;
545  }
546  mg_param.location[i] = QUDA_CUDA_FIELD_LOCATION;
547  mg_param.smoother_solve_type[i] = QUDA_DIRECT_PC_SOLVE;
548  mg_param.coarse_grid_solution_type[i] = QUDA_MATPC_SOLUTION;
549  }
550 
551  // only coarsen the spin on the first restriction
552  mg_param.spin_block_size[0] = 2;
553 
554  // coarse grid solver is GCR
555  mg_param.smoother[ip.mg_levels-1] = QUDA_GCR_INVERTER;
556 
557  mg_param.compute_null_vector = ip.generate_nullspace ? QUDA_COMPUTE_NULL_VECTOR_YES
558  : QUDA_COMPUTE_NULL_VECTOR_NO;
559 
560  for(int l=0; l < ip.mg_levels; l++) {
561  mg_param.vec_infile[l][0] = '\0';
562  mg_param.vec_outfile[l][0] = '\0';
563  }
564 
565  QDPIO::cout<<"Basic MULTIGRID params copied."<<std::endl;
566  quda_inv_param.verbosity = QUDA_VERBOSE;
567 
568 
569 
570  // setup the multigrid solver
571  void *mg_preconditioner = newMultigridQuda(&mg_param);
572  QDPIO::cout<<"NewMultigridQuda state initialized."<<std::endl;
573  quda_inv_param.preconditioner = mg_preconditioner;
574  QDPIO::cout<<"MULTIGRID preconditioner set."<<std::endl;
575  //
576 
577 
578  }
579 
580 
581  //! Destructor is automatic
582  ~LinOpSysSolverQUDAMULTIGRIDWilson()
583  {
584  QDPIO::cout << "Destructing" << std::endl;
585  freeGaugeQuda();
586 
587  /* FIXME: I Don't use it. Should I still do this? */
588  // freeCloverQuda();
589  destroyMultigridQuda(quda_inv_param.preconditioner);
590  }
591 
592  //! Return the subset on which the operator acts
593  const Subset& subset() const {return A->subset();}
594 
595  //! Solver the linear system
596  /*!
597  * \param psi solution ( Modify )
598  * \param chi source ( Read )
599  * \return syssolver results
600  */
601  SystemSolverResults_t operator() (T& psi, const T& chi) const
602  {
603  SystemSolverResults_t res;
604 
605  START_CODE();
606  StopWatch swatch;
607  swatch.start();
608 
609  psi=zero;
610 
611  // T MdagChi;
612 
613  // This is a CGNE. So create new RHS
614  // (*A)(MdagChi, chi, MINUS);
615  // Handle< LinearOperator<T> > MM(new MdagMLinOp<T>(A));
616  if ( invParam.axialGaugeP ) {
617  T g_chi,g_psi;
618 
619  // Gauge Fix source and initial guess
620  QDPIO::cout << "Gauge Fixing source and initial guess" << std::endl;
621  g_chi[ rb[1] ] = GFixMat * chi;
622  g_psi[ rb[1] ] = GFixMat * psi;
623  QDPIO::cout << "Solving" << std::endl;
624  res = qudaInvert( g_chi,g_psi);
625  QDPIO::cout << "Untransforming solution." << std::endl;
626  psi[ rb[1]] = adj(GFixMat)*g_psi;
627 
628  }
629  else {
630  res = qudaInvert(chi,psi);
631  }
632 
633  swatch.stop();
634 
635 
636  {
637  T r;
638  r[A->subset()]=chi;
639  T tmp;
640  (*A)(tmp, psi, PLUS);
641  r[A->subset()] -= tmp;
642  res.resid = sqrt(norm2(r, A->subset()));
643  }
644 
645  Double rel_resid = res.resid/sqrt(norm2(chi,A->subset()));
646 
647  QDPIO::cout << "QUDA_MULTIGRID_"<< solver_string <<"_WILSON_SOLVER: " << res.n_count << " iterations. Rsd = " << res.resid << " Relative Rsd = " << rel_resid << std::endl;
648 
649  // Convergence Check/Blow Up
650  if ( ! invParam.SilentFailP ) {
651  if ( toBool( rel_resid > invParam.RsdToleranceFactor*invParam.RsdTarget) ) {
652  QDPIO::cerr << "ERROR: QUDA MULTIGRID Solver residuum is outside tolerance: QUDA resid="<< rel_resid << " Desired =" << invParam.RsdTarget << " Max Tolerated = " << invParam.RsdToleranceFactor*invParam.RsdTarget << std::endl;
653  QDP_abort(1);
654  }
655  }
656 
657  END_CODE();
658  return res;
659  }
660 
661 
662 private:
663  // Hide default constructor
664  LinOpSysSolverQUDAMULTIGRIDWilson() {}
665 
666 #if 1
667  Q links_orig;
668 #endif
669 
670  U GFixMat;
671  QudaPrecision_s cpu_prec;
672  QudaPrecision_s gpu_prec;
673  QudaPrecision_s gpu_half_prec;
674 
675  Handle< LinearOperator<T> > A;
676  const SysSolverQUDAMULTIGRIDWilsonParams invParam;
677  QudaGaugeParam q_gauge_param;
678  QudaInvertParam quda_inv_param;
679  QudaInvertParam mg_inv_param;
680  QudaMultigridParam mg_param;
681 
682  SystemSolverResults_t qudaInvert( const T& chi_s,
683  T& psi_s
684  )const ;
685 
686  std::string solver_string;
687 };
688 
689 
690 } // End namespace
691 
692 #endif // BUILD_QUDA
693 #endif
694 
aniso_io.h
Anisotropy parameters.
chroma_config.h
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.
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::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::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
l
int l
Definition: pade_trln_w.cc:111
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_multigrid_wilson_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