Konstantin devel Wilson+Clover N_f=2+1 implementation

monomial/clover_trlog_monomial.c

@@ -54,7 +54,7 @@ void clover_trlog_heatbath(const int id, hamiltonian_field_t* const hf) {
   init_sw_fields();
   sw_term((const su3**)hf->gaugefield, mnl->kappa, mnl->c_sw);
   /*compute the contribution from the clover trlog term */
-  mnl->energy0 = -sw_trace(EO, mnl->mu);
+  mnl->energy0 = -sw_trace(EO, mnl->mu, mnl->single_flavor);
   if (g_proc_id == 0) {
     if (g_debug_level > 3) {
       printf("called clover_trlog_heatbath for id %d E = %e\n", id, mnl->energy0);
@@ -70,7 +70,7 @@ double clover_trlog_acc(const int id, hamiltonian_field_t* const hf) {
   mnl->energy1 = 0.;
   sw_term((const su3**)hf->gaugefield, mnl->kappa, mnl->c_sw);
   /*compute the contribution from the clover trlog term */
-  mnl->energy1 = -sw_trace(EO, mnl->mu);
+  mnl->energy1 = -sw_trace(EO, mnl->mu, mnl->single_flavor);
   if (g_proc_id == 0 && g_debug_level > 3) {
     if (g_debug_level > 3) {
       printf("called clover_trlog_acc for id %d dH = %1.10e\n", id, mnl->energy1 - mnl->energy0);

monomial/cloverdet_monomial.c

@@ -187,7 +187,7 @@ void cloverdet_derivative(const int id, hamiltonian_field_t *const hf) {
       // we again compute only the insertion matrices for S_det
       // the result is added to swp and swm
       // even sites only!
-      sw_deriv(EE, mnl->mu);
+      sw_deriv(EE, mnl->mu, 0);
     } else {
       /* \delta Q sandwitched by Y^\dagger and X */
       deriv_Sb_D_psi(mnl->w_fields[0], mnl->w_fields[1], hf, mnl->forcefactor);

monomial/monomial.c

@@ -96,6 +96,15 @@ int add_monomial(const int type) {
   monomial_list[no_monomials].forceprec = _default_g_eps_sq_force;
   monomial_list[no_monomials].maxiter = _default_max_solver_iterations;
   monomial_list[no_monomials].HB_maxiter = _default_max_solver_iterations;
+  if (monomial_list[no_monomials].type == RAT || monomial_list[no_monomials].type == RATCOR ||
+      monomial_list[no_monomials].type == CLOVERRAT ||
+      monomial_list[no_monomials].type == CLOVERRATCOR) {
+    monomial_list[no_monomials].single_flavor = 1;
+  } else {
+    monomial_list[no_monomials].single_flavor =
+        0;  // NOTE: the special case of CLOVERTRLOG requires knowledge of the "parent" monomial and
+            // is dealt with in init_monomials()
+  }
   if ((monomial_list[no_monomials].type == NDRAT) ||
       (monomial_list[no_monomials].type == NDRATCOR) ||
       (monomial_list[no_monomials].type == NDCLOVERRAT) ||
@@ -107,7 +116,7 @@ int add_monomial(const int type) {
   }
   monomial_list[no_monomials].solver_params.mcg_delta = _default_mixcg_innereps;
   monomial_list[no_monomials].solver_params.solution_type = TM_SOLUTION_M_MDAG;
-  // the defaut is 1 because the QPhiX interface is generalised in such a way
+  // the default is 1 because the QPhiX interface is generalised in such a way
   // that normal solves correspond to solves with one shift, this does not
   // affect the used parameters in any way!
   monomial_list[no_monomials].solver_params.no_shifts = 1;
@@ -578,6 +587,10 @@ int init_monomials(const int V, const int even_odd_flag) {
       monomial_list[no_monomials - 1].derivativefunction = NULL;
       monomial_list[no_monomials - 1].timescale = 0;
       monomial_list[no_monomials - 1].even_odd_flag = even_odd_flag;
+      if (monomial_list[clover_monomials[j]].type == CLOVERRAT)
+        monomial_list[no_monomials - 1].single_flavor =
+            1;  // only if the parent monomial is a CLOVERRAT (with "AddTrLog = yes") it is for a
+                // single flavor.
       if (g_proc_id == 0 && g_debug_level > 1) {
         printf("# Initialised monomial %s of type CLOVERTRLOG, no_monomials= %d\n",
                monomial_list[no_monomials - 1].name, no_monomials);

monomial/monomial.h

@@ -30,34 +30,36 @@
 
 #include <complex.h>
 
-#define DET 0
-#define DETRATIO 1
-#define GAUGE 2
-#define POLY 3
-#define NDPOLY 4
-#define SFGAUGE 5
-#define NDDETRATIO 6
-#define POLYDETRATIO 7
-#define CLOVERTRLOG 8
-#define CLOVERDET 9
-#define CLOVERDETRATIO 10
-#define NDCLOVER 11
-#define CLOVERNDTRLOG 12
-#define NDRAT 13
-#define NDCLOVERRAT 14
-#define NDRATCOR 15
-#define NDCLOVERRATCOR 16
-#define RAT 17
-#define RATCOR 18
-#define CLOVERRAT 19
-#define CLOVERRATCOR 20
-#define CLOVERDETRATIORW 21
-#define NDCLOVERDETRATIO 22
+typedef enum {
+  DET,
+  DETRATIO,
+  GAUGE,
+  POLY,
+  NDPOLY,
+  SFGAUGE,
+  NDDETRATIO,
+  POLYDETRATIO,
+  CLOVERTRLOG,
+  CLOVERDET,
+  CLOVERDETRATIO,
+  NDCLOVER,
+  CLOVERNDTRLOG,
+  NDRAT,
+  NDCLOVERRAT,
+  NDRATCOR,
+  NDCLOVERRATCOR,
+  RAT,
+  RATCOR,
+  CLOVERRAT,
+  CLOVERRATCOR,
+  CLOVERDETRATIORW,
+  NDCLOVERDETRATIO,
+} monomial_t;
 
 #define max_no_monomials 30
 
 typedef struct {
-  int type;
+  monomial_t type;
   int gtype;
   int initialised;
   int timescale;
@@ -121,6 +123,7 @@ typedef struct {
   double *MDPolyCoefs, *PtildeCoefs;
   /* rational approximation */
   rational_t rat;
+  int single_flavor;
   /* chronological solver fields */
   spinor **csg_field;
   spinor **csg_field2;

monomial/rat_monomial.c

@@ -143,7 +143,7 @@ void rat_derivative(const int id, hamiltonian_field_t* const hf) {
     }
   }
   if (mnl->type == CLOVERRAT && mnl->trlog) {
-    sw_deriv(EE, 0.);
+    sw_deriv(EE, 0., 1);
   }
   if (mnl->type == CLOVERRAT) {
     sw_all(hf, mnl->kappa, mnl->c_sw);
@@ -173,8 +173,14 @@ void rat_heatbath(const int id, hamiltonian_field_t* const hf) {
   }
   // we measure before the trajectory!
   if ((mnl->rec_ev != 0) && (hf->traj_counter % mnl->rec_ev == 0)) {
-    // if(mnl->type != CLOVERRAT) phmc_compute_ev(hf->traj_counter-1, id, &Qtm_pm_ndbipsi);
-    // else phmc_compute_ev(hf->traj_counter-1, id, &Qsw_pm_ndbipsi);
+    if (mnl->type != CLOVERRAT) {
+      phmc_compute_ev(hf->traj_counter - 1, id, &Qtm_pm_ndbipsi);
+    } else {
+      phmc_compute_ev(
+          hf->traj_counter - 1, id,
+          &Qsw_pm_ndbipsi);  // not that this one is technically not correct, but works anyways as
+                             // the eval computation is (only) implemented using QUDA.
+    }
   }
 
   // the Gaussian distributed random fields

operator/clover_deriv.c

@@ -60,7 +60,7 @@
 //
 // this function depends on mu
 
-void sw_deriv(const int ieo, const double mu) {
+void sw_deriv(const int ieo, const double mu, const int single_flavor) {
   tm_stopwatch_push(&g_timers, __func__, "");
 #ifdef TM_USE_OMP
 #pragma omp parallel
@@ -79,6 +79,9 @@ void sw_deriv(const int ieo, const double mu) {
       ioff = (VOLUME + RAND) / 2;
     }
     if (fabs(mu) > 0.) fac = 0.5;
+    if (single_flavor) fac = 0.5;
+    if (single_flavor && (fabs(mu) > 0.))
+      fatal_error("single_flavor flag and |mu|>0 are incompatible.", __func__);
 
 #ifndef TM_USE_OMP
     icy = 0;

operator/clover_det.c

@@ -100,13 +100,15 @@ void six_det(_Complex double* const rval, _Complex double a[6][6]) {
   *rval = det;
 }
 
-double sw_trace(const int ieo, const double mu) {
+double sw_trace(const int ieo, const double mu, const int single_flavor) {
   tm_stopwatch_push(&g_timers, __func__, "");
   double ALIGN res = 0.0;
 #ifdef TM_USE_MPI
   double ALIGN mres;
 #endif
 
+  const double fac = (single_flavor ? 0.5 : 1.);
+
 #ifdef TM_USE_OMP
 #pragma omp parallel
   {
@@ -171,10 +173,10 @@ double sw_trace(const int ieo, const double mu) {
 #ifdef TM_USE_MPI
   MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
   tm_stopwatch_pop(&g_timers, 0, 1, "");
-  return (mres);
+  return (fac * mres);
 #else
   tm_stopwatch_pop(&g_timers, 0, 1, "");
-  return (res);
+  return (fac * res);
 #endif
 }
 

operator/clover_leaf.h

@@ -29,13 +29,13 @@ extern su3 **swm, **swp;
 extern const double tiny_t;
 
 void sw_term(const su3** const gf, const double kappa, const double c_sw);
-double sw_trace(const int ieo, const double mu);
+double sw_trace(const int ieo, const double mu, const int single_flavor);
 double sw_trace_nd(const int ieo, const double mu, const double eps);
 void sw_invert(const int ieo, const double mu);
 void sw_invert_nd(const double mshift);
 void sw_invert_epsbar(const double epsbar);
 void sw_invert_mubar(const double mubar);
-void sw_deriv(const int ieo, const double mu);
+void sw_deriv(const int ieo, const double mu, const int single_flavor);
 void sw_deriv_nd(const int ieo);
 void sw_spinor_eo(const int ieo, const spinor* const kk, const spinor* const ll, const double fac);
 void sw_spinor(const spinor* const kk, const spinor* const ll, const double fac);

phmc.c

@@ -211,18 +211,16 @@ void phmc_compute_ev(const int trajectory_counter, const int id, matrix_mult_bi
   double atime, etime;
   _Complex double eval_min = 0.0;
   _Complex double eval_max = 0.0;
-  int max_iter_ev, no_eigenvalues;
+  int max_iter_ev = 1000;
+  int no_eigenvalues;
   char buf[100];
   char *phmcfilename = buf;
   FILE *countfile;
   monomial *mnl = &monomial_list[id];
-  ;
 
   sprintf(phmcfilename, "monomial-%.2d.data", id);
   atime = gettime();
 
-  max_iter_ev = 1000;
-
   if ((g_proc_id == 0) && (g_debug_level > 0)) {
     printf("# Computing eigenvalues for heavy doublet\n");
   }
@@ -235,7 +233,7 @@ void phmc_compute_ev(const int trajectory_counter, const int id, matrix_mult_bi
                  mnl->eig_n_kr, mnl->solver, g_relative_precision_flag,
                  1,  // we only support even-odd here
                  mnl->solver_params.refinement_precision, mnl->solver_params.sloppy_precision,
-                 mnl->solver_params.compression_type, 0);
+                 mnl->solver_params.compression_type, mnl->single_flavor);
     if (fabs(mnl->EVMax - 1) < 2 * DBL_EPSILON) {
       eval_min /= mnl->StildeMax;
     }
@@ -244,13 +242,24 @@ void phmc_compute_ev(const int trajectory_counter, const int id, matrix_mult_bi
       fprintf(stderr,
               "Error: Attempted to use QUDA eigensolver but this build was not configured for QUDA "
               "usage.\n");
+    }
+#ifdef TM_USE_MPI
+    MPI_Finalize();
+#endif
+    exit(-2);
+#endif
+  } else {
+    if (mnl->single_flavor) {
+      if (g_proc_id == 0) {
+        fprintf(stderr,
+                "Error: CPU version of eigenvalue computation not yet implemented for single quark "
+                "flavor.");
+      }
 #ifdef TM_USE_MPI
       MPI_Finalize();
 #endif
       exit(-2);
     }
-#endif
-  } else {
     eval_min = eigenvalues_bi(&no_eigenvalues, max_iter_ev, eigenvalue_precision, 0, Qsq);
   }
 
@@ -262,7 +271,7 @@ void phmc_compute_ev(const int trajectory_counter, const int id, matrix_mult_bi
                  mnl->eig_n_kr, mnl->solver, g_relative_precision_flag,
                  1,  // we only support even-odd here
                  mnl->solver_params.refinement_precision, mnl->solver_params.sloppy_precision,
-                 mnl->solver_params.compression_type, 0);
+                 mnl->solver_params.compression_type, mnl->single_flavor);
     if (fabs(mnl->EVMax - 1.) < 2 * DBL_EPSILON) {
       eval_max /= mnl->StildeMax;
     }
@@ -271,13 +280,24 @@ void phmc_compute_ev(const int trajectory_counter, const int id, matrix_mult_bi
       fprintf(stderr,
               "Error: Attempted to use QUDA eigensolver but this build was not configured for QUDA "
               "usage.\n");
+    }
+#ifdef TM_USE_MPI
+    MPI_Finalize();
+#endif
+    exit(-2);
+#endif
+  } else {
+    if (mnl->single_flavor) {
+      if (g_proc_id == 0) {
+        fprintf(stderr,
+                "Error: CPU version of eigenvalue computation not yet implemented for single quark "
+                "flavor.");
+      }
 #ifdef TM_USE_MPI
       MPI_Finalize();
 #endif
       exit(-2);
     }
-#endif
-  } else {
     eval_max = eigenvalues_bi(&no_eigenvalues, max_iter_ev, eigenvalue_precision, 1, Qsq);
   }
 

quda_interface.c

@@ -1423,10 +1423,13 @@ void _setOneFlavourSolverParam(const double kappa, const double c_sw, const doub
       inv_param.mu = 0.0;
       inv_param.tm_rho = -g_mu3 / 2. / kappa;
     } else {
+      inv_param.mu = 0.0;
+      inv_param.tm_rho = 0.0;
+      inv_param.clover_rho = 0.0;
       inv_param.twist_flavor = QUDA_TWIST_NO;
       inv_param.dslash_type = QUDA_CLOVER_WILSON_DSLASH;
     }
-    inv_param.matpc_type = QUDA_MATPC_EVEN_EVEN;
+    inv_param.matpc_type = QUDA_MATPC_ODD_ODD_ASYMMETRIC;
     inv_param.solution_type = QUDA_MAT_SOLUTION;
     inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
     inv_param.clover_coeff = c_sw * kappa;

reweighting_factor.c

@@ -65,14 +65,14 @@ void reweighting_factor(const int N, const int nstore) {
 
         sw_term((const su3**)hf.gaugefield, mnl->kappa, c_sw);
         if (mnl->type == CLOVERDETRATIORW) {
-          trlog[j] = -sw_trace(0, mnl->mu);
+          trlog[j] = -sw_trace(0, mnl->mu, 0);
         } else {
           trlog[j] = -sw_trace_nd(0, mnl->mubar, mnl->epsbar);
         }
 
         sw_term((const su3**)hf.gaugefield, mnl->kappa2, c_sw);
         if (mnl->type == CLOVERDETRATIORW) {
-          trlog[j] -= -sw_trace(0, mnl->mu2);
+          trlog[j] -= -sw_trace(0, mnl->mu2, 0);
         } else {
           trlog[j] -= -sw_trace_nd(0, mnl->mubar2, mnl->epsbar2);
         }