Merge pull request #1366 from lattice/feature/quark-smearing.milc-int…

…erface MILC interface for two-link quark smearing
lattice · Aug 30, 2023 · d571552 · d571552
2 parents 1d6af2d + 073fd0f
commit d571552
Show file tree

Hide file tree

Showing 2 changed files with 87 additions and 0 deletions.
diff --git a/include/quda_milc_interface.h b/include/quda_milc_interface.h
@@ -118,6 +118,18 @@ extern "C" {
     int use_pinned_memory;   /** use page-locked memory in QUDA    */
   } QudaFatLinkArgs_t;
 
+  /**
+   * Parameters for two-link Gaussian quark smearing.
+   */
+  typedef struct {
+    int n_steps; /** Number of steps to apply **/
+    double width; /** The width of the Gaussian **/
+    int compute_2link; /** if nonzero then compute two-link, otherwise reuse gaugeSmeared **/
+    int delete_2link; /** if nonzero then delete two-link, otherwise keep two-link for future use **/
+    int t0; /** Set if the input spinor is on a time slice **/
+    int laplaceDim; /** Dimension of Laplacian **/
+  } QudaTwoLinkQuarkSmearArgs_t;
+
   /**
    * Optional: Set the MPI Comm Handle if it is not MPI_COMM_WORLD
    *
@@ -1061,6 +1073,17 @@ extern "C" {
     void* milc_sitelink
     );
 
+  /**
+   * @brief Perform two-link Gaussian smearing on a given spinor (for staggered fermions).
+   * @param[in] external_precision  Precision of host fields passed to QUDA (2 - double, 1 - single)
+   * @param[in] quda_precision  Precision for QUDA to use (2 - double, 1 - single)
+   * @param[in] h_gauge  Host gauge field
+   * @param[in,out] source  Spinor field to smear
+   * @param[in] qsmear_args  Struct setting some smearing metadata
+   */
+  void qudaTwoLinkGaussianSmear(int external_precision, int quda_precision, void * h_gauge, void * source,
+                                QudaTwoLinkQuarkSmearArgs_t qsmear_args);
+
   /* The below declarations are for removed functions from prior versions of QUDA. */
 
   /**

diff --git a/lib/milc_interface.cpp b/lib/milc_interface.cpp
@@ -3083,3 +3083,67 @@ void qudaGaugeFixingFFT( int precision,
   printfQuda("Time D2H: %lf\n", timeinfo[2]);
   printfQuda("Time all: %lf\n", timeinfo[0]+timeinfo[1]+timeinfo[2]);
 }
+
+void qudaTwoLinkGaussianSmear( int external_precision, int quda_precision, void * h_gauge, void * source, QudaTwoLinkQuarkSmearArgs_t qsmear_args )
+{
+  static const QudaVerbosity verbosity = getVerbosity();
+  qudamilc_called<true>(__func__, verbosity);
+
+  QudaPrecision cpu_prec = (external_precision == 2) ? QUDA_DOUBLE_PRECISION : QUDA_SINGLE_PRECISION;
+  QudaPrecision cuda_prec = (quda_precision == 2) ? QUDA_DOUBLE_PRECISION : QUDA_SINGLE_PRECISION;
+  QudaPrecision cuda_prec_sloppy = cuda_prec;
+
+  // inverter setup ---------------------
+  QudaInvertParam invertParam = newQudaInvertParam();
+
+  double mass = 0.0;
+  QudaParity parity = QUDA_EVEN_PARITY; // need to fix
+
+  // dummies
+  double target_residual = 1e-14;
+  double target_residual_hq = 1e-14;
+  int maxiter = 1;
+  double reliable_delta = 0.1;
+  QudaInverterType inverter = QUDA_CG_INVERTER;
+
+  setInvertParams( cpu_prec, cuda_prec, cuda_prec_sloppy, mass, target_residual, target_residual_hq, maxiter, reliable_delta, parity, verbosity, inverter, &invertParam );
+
+  invertParam.laplace3D = qsmear_args.laplaceDim;
+  //---------------------------- inverter setup
+
+  // gauge setup ---------------------------
+  QudaGaugeParam gaugeParam = newQudaGaugeParam();
+
+  int * dim = localDim;
+
+  // dummies
+  double tadpole = 0;
+  double naik_epsilon = 0;
+
+  setGaugeParams( gaugeParam, gaugeParam, nullptr, dim, cpu_prec, cuda_prec, cuda_prec_sloppy, tadpole, naik_epsilon );
+
+  gaugeParam.reconstruct = QUDA_RECONSTRUCT_NO; // need to fix
+  gaugeParam.staggered_phase_type = QUDA_STAGGERED_PHASE_NO;
+
+  gaugeParam.ga_pad = getLinkPadding( dim );
+  gaugeParam.mom_ga_pad = 0;
+  //--------------------------- gauge setup
+
+  // Load gauge field
+  if( qsmear_args.compute_2link == 0 ) gaugeParam.use_resident_gauge = 1;
+  loadGaugeQuda( const_cast<void *>(h_gauge), &gaugeParam );
+
+  // quark smearing parameters
+  QudaQuarkSmearParam qsmearParam;
+  qsmearParam.inv_param = &invertParam;
+  qsmearParam.n_steps = qsmear_args.n_steps;
+  qsmearParam.width = qsmear_args.width;
+  qsmearParam.compute_2link = qsmear_args.compute_2link;
+  qsmearParam.delete_2link = qsmear_args.delete_2link;
+  qsmearParam.t0 = qsmear_args.t0;
+
+  // run gaussian smearing
+  performTwoLinkGaussianSmearNStep( source, &qsmearParam );
+
+  qudamilc_called<false>(__func__, verbosity);  
+} //qudaTwoLinkGaussianSmear