merge with quda_work

default_input_values.h

@@ -147,6 +147,10 @@
 #define _default_phmc_pure_phmc 0
 #define _default_stilde_max 3.
 #define _default_stilde_min 0.01
+#define _default_eig_polydeg 128
+#define _default_eig_amin 0.001
+#define _default_eig_amax 4
+#define _default_eig_n_kr 96
 #define _default_degree_of_p 48
 #define _default_propagator_splitted 1
 #define _default_source_splitted 1
@@ -199,12 +203,14 @@
 
 #define _default_external_inverter 0
 
+#define _default_external_eigsolver 0
+
 #define _default_external_library 0
 
 #define _default_subprocess_flag 0
 #define _default_lowmem_flag 0
 
-#define _default_g_barrier_monomials_convergence 0
+#define _default_g_barrier_monomials_convergence 1
 
 /* default input values for QUDA interface */
 /* These follow the recommendations of https://github.com/lattice/quda/wiki/Multigrid-Solver */

doc/bibliography.bib

@@ -7904,3 +7904,17 @@ @inbook{Joo2013
 	year = "2013"
 }
 
+@article{Kostrzewa:2022hsv,
+    author = "Kostrzewa, Bartosz and Bacchio, Simone and Finkenrath, Jacob and Garofalo, Marco and Pittler, Ferenc and Romiti, Simone and Urbach, Carsten",
+    collaboration = "ETM",
+    title = "{Twisted mass ensemble generation on GPU machines}",
+    eprint = "2212.06635",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-lat",
+    doi = "10.22323/1.430.0340",
+    journal = "PoS",
+    volume = "LATTICE2022",
+    pages = "340",
+    year = "2023"
+}
+

include/tmLQCD.h

@@ -83,6 +83,12 @@ int tmLQCD_read_gauge(const int nconfig);
 int tmLQCD_invert(double *const propagator, double *const source, const int op_id,
                   const int write_prop);
 
+// invert with source and propagator provided in TXYZ spin colour complex lexicographic order
+// propagator has kappa normalisation
+// the two propagators and sources correspond to two flavours
+int tmLQCD_invert_doublet(double* const propagator0, double* const propagator1, double* const source0, 
+                  double* const source1, const int op_id, const int write_prop);
+
 // invert on odd part of lattice with prepared source
 int tmLQCD_invert_eo(double *const Odd_out, double *const Odd_in, const int op_id);
 

misc_types.h

@@ -89,6 +89,12 @@ typedef enum ExternalInverter_s {
   QPHIX_INVERTER
 } ExternalInverter;
 
+/* enumeration type for the external eigensolver */
+typedef enum ExternalEigSolver_s {
+  NO_EXT_EIGSOLVER = 0,
+  QUDA_EIGSOLVER
+} ExternalEigSolver;
+
 /* enumeration type for the external inverter */
 typedef enum ExternalLibrary_s {
   NO_EXT_LIB = 0,

monomial/monomial.c

@@ -116,6 +116,7 @@ int add_monomial(const int type) {
   monomial_list[no_monomials].solver_params.external_inverter = _default_external_inverter;
   monomial_list[no_monomials].solver_params.sloppy_precision = _default_operator_sloppy_precision_flag;
   monomial_list[no_monomials].external_library = _default_external_library;
+  monomial_list[no_monomials].external_eigsolver = _default_external_eigsolver;
   monomial_list[no_monomials].solver_params.refinement_precision = _default_operator_sloppy_precision_flag;
   monomial_list[no_monomials].HB_solver_params = monomial_list[no_monomials].solver_params;
   monomial_list[no_monomials].even_odd_flag = _default_even_odd_flag;
@@ -143,6 +144,11 @@ int add_monomial(const int type) {
   monomial_list[no_monomials].PrecisionHfinal = _default_g_acc_Hfin;
   monomial_list[no_monomials].PrecisionPtilde = _default_g_acc_Ptilde;
 
+  monomial_list[no_monomials].eig_polydeg = _default_eig_polydeg;
+  monomial_list[no_monomials].eig_amin = _default_eig_amin;
+  monomial_list[no_monomials].eig_amax = _default_eig_amax;
+  monomial_list[no_monomials].eig_n_kr = _default_eig_n_kr;
+
   monomial_list[no_monomials].rat.order = 12;
   monomial_list[no_monomials].rat.range[0] = _default_stilde_min;
   monomial_list[no_monomials].rat.range[1] = _default_stilde_max;

monomial/monomial.h

@@ -114,7 +114,10 @@ typedef struct {
   double PrecisionHfinal;
   double StildeMin, StildeMax;
   double EVMin, EVMax, EVMaxInv;
+  int eig_polydeg, eig_n_kr;
+  double eig_amin, eig_amax;
   ExternalLibrary external_library;
+  ExternalEigSolver external_eigsolver;
   double * MDPolyCoefs, * PtildeCoefs;
   /* rational approximation */
   rational_t rat;

operator.c

@@ -410,6 +410,12 @@ void op_invert(const int op_id, const int index_start, const int write_prop) {
         break;
     }
   } else if(optr->type == DBTMWILSON || optr->type == DBCLOVER) {
+    // there is no default for this set anywhere other than prepare_source.c
+    // such that calling this branch from an external program via tmLQCD_invert_doublet
+    // would result in undefined behaviour
+    // we thus set a default here unless it's explicitly set to 2
+    if( SourceInfo.no_flavours != 2 ) SourceInfo.no_flavours = 1;
+
     if(optr->type == DBCLOVER) {
       if (g_cart_id == 0 && g_debug_level > 1) {
         printf("#\n# csw = %e, computing clover leafs\n", g_c_sw);

phmc.c

@@ -25,6 +25,7 @@
 #include <math.h>
 #include <string.h>
 #include <time.h>
+#include <float.h>
 
 #include "global.h"
 
@@ -40,6 +41,10 @@
 #include "solver/matrix_mult_typedef_bi.h"
 #include "gettime.h"
 
+#ifdef TM_USE_QUDA
+#  include "quda_interface.h"
+#endif
+
 //                                          --> in  monomial
 double phmc_Cpol;                        // --> MDPolyLocNormConst
 double phmc_cheb_evmin, phmc_cheb_evmax; // --> EVMin, EVMax
@@ -206,7 +211,9 @@ void init_phmc() {
 void phmc_compute_ev(const int trajectory_counter,
 		     const int id,
 		     matrix_mult_bi Qsq) {
-  double atime, etime, temp=0., temp2=0.;
+  double atime, etime;
+  _Complex double eval_min = 0.0;
+  _Complex double eval_max = 0.0;
   int max_iter_ev, no_eigenvalues;
   char buf[100];
   char * phmcfilename = buf;
@@ -223,28 +230,75 @@ void phmc_compute_ev(const int trajectory_counter,
   }
 
   no_eigenvalues = 1;
-
-  temp = eigenvalues_bi(&no_eigenvalues, max_iter_ev, eigenvalue_precision, 0, Qsq);
+  if(mnl->external_eigsolver == QUDA_EIGSOLVER) {
+  #ifdef TM_USE_QUDA
+    eigsolveQuda(&eval_min, no_eigenvalues, eigenvalue_precision, 1, 0, max_iter_ev, 0,
+                 mnl->accprec, mnl->maxiter, mnl->eig_polydeg, mnl->eig_amin, 
+                 mnl->eig_amax, 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);
+    if( fabs(mnl->EVMax - 1) < 2*DBL_EPSILON ) {
+      eval_min /= mnl->StildeMax;
+    }
+  #else
+    if(g_proc_id == 0) {
+      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 {
+    eval_min = eigenvalues_bi(&no_eigenvalues, max_iter_ev, eigenvalue_precision, 0, Qsq);
+  }
+
 
   no_eigenvalues = 1;
-  temp2 = eigenvalues_bi(&no_eigenvalues, max_iter_ev, eigenvalue_precision, 1, Qsq);
+  if(mnl->external_eigsolver == QUDA_EIGSOLVER) {
+  #ifdef TM_USE_QUDA
+    eigsolveQuda(&eval_max, no_eigenvalues, eigenvalue_precision, 1, 0, max_iter_ev, 1,
+                 mnl->accprec, mnl->maxiter, mnl->eig_polydeg, mnl->eig_amin, 
+                 mnl->eig_amax, 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);
+    if( fabs(mnl->EVMax - 1.) < 2*DBL_EPSILON ) {
+      eval_max /= mnl->StildeMax;
+    }
+  #else
+    if(g_proc_id == 0) {
+      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 {
+    eval_max = eigenvalues_bi(&no_eigenvalues, max_iter_ev, eigenvalue_precision, 1, Qsq);
+  }
+
 
   if((g_proc_id == 0) && (g_debug_level > 1)) {
     printf("# %s: lowest eigenvalue end of trajectory %d = %e\n", 
-	   mnl->name, trajectory_counter, temp);
+	   mnl->name, trajectory_counter, creal(eval_min));
     printf("# %s: maximal eigenvalue end of trajectory %d = %e\n", 
-	   mnl->name, trajectory_counter, temp2);
+	   mnl->name, trajectory_counter, creal(eval_max));
   }
   if(g_proc_id == 0) {
-    if(temp2 > mnl->EVMax) {
-      fprintf(stderr, "\nWarning: largest eigenvalue for monomial %s larger than upper bound!\n\n", mnl->name);
+    if(creal(eval_max) > mnl->EVMax) {
+      fprintf(stderr, "\nWarning: largest eigenvalue for monomial %s: %.6f is larger than upper bound: %.6f\n\n", mnl->name, creal(eval_max), mnl->EVMax);
     }
-    if(temp < mnl->EVMin) {
-      fprintf(stderr, "\nWarning: smallest eigenvalue for monomial %s smaller than lower bound!\n\n", mnl->name);
+    if(creal(eval_min) < mnl->EVMin) {
+      fprintf(stderr, "\nWarning: smallest eigenvalue for monomial %s: %.6f is smaller than lower bound: %.6f\n\n", mnl->name, creal(eval_min), mnl->EVMin);
     }
     countfile = fopen(phmcfilename, "a");
     fprintf(countfile, "%.8d %1.5e %1.5e %1.5e %1.5e\n", 
-	    trajectory_counter, temp, temp2, mnl->EVMin, mnl->EVMax);
+	    trajectory_counter, creal(eval_min), creal(eval_max), mnl->EVMin, mnl->EVMax);
     fclose(countfile);
   }
   etime = gettime();

quda_dummy_types.h

@@ -39,6 +39,11 @@ typedef enum QudaInverterType_s {
  QUDA_CG_INVERTER = CG,
  QUDA_MR_INVERTER = MR,
  QUDA_GCR_INVERTER = GCR,
+ QUDA_CGNE_INVERTER = CGNE,
+ QUDA_CGNR_INVERTER = CGNR,
+ QUDA_CA_CG_INVERTER = CA_CG,
+ QUDA_CA_CGNE_INVERTER = CA_CGNE, 
+ QUDA_CA_CGNR_INVERTER = CA_CGNR,
  QUDA_CA_GCR_INVERTER = CA_GCR
 } QudaInverterType;