user input for poly acc settings

etmc · Jun 17, 2023 · 3d67264 · 3d67264
1 parent 4e6033f
commit 3d67264
Show file tree

Hide file tree

Showing 7 changed files with 45 additions and 9 deletions.
diff --git a/default_input_values.h b/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

diff --git a/monomial/monomial.c b/monomial/monomial.c
@@ -144,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;

diff --git a/monomial/monomial.h b/monomial/monomial.h
@@ -114,6 +114,8 @@ 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;

diff --git a/phmc.c b/phmc.c
@@ -231,7 +231,8 @@ void phmc_compute_ev(const int trajectory_counter,
   if(mnl->external_eigsolver == QUDA_EIGSOLVER) {
   #ifdef TM_USE_QUDA
     temp = eigsolveQuda(no_eigenvalues, eigenvalue_precision, 1, 0, max_iter_ev, 0,
-                        mnl->accprec, mnl->maxiter, mnl->solver, g_relative_precision_flag,
+                        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,
@@ -257,7 +258,8 @@ void phmc_compute_ev(const int trajectory_counter,
   if(mnl->external_eigsolver == QUDA_EIGSOLVER) {
   #ifdef TM_USE_QUDA
     temp2 = eigsolveQuda(no_eigenvalues, eigenvalue_precision, 1, 0, max_iter_ev, 1,
-                         mnl->accprec, mnl->maxiter, mnl->solver, g_relative_precision_flag,
+                         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,

diff --git a/quda_interface.c b/quda_interface.c
@@ -2598,7 +2598,8 @@ Interface function for Eigensolver on Quda
 
 
 double eigsolveQuda(int n, double tol, int blksize, int blkwise, int max_iterations, int maxmin,
-                    const double precision, const int max_iter, const int solver_flag, const int rel_prec,
+                    const double precision, const int max_iter, const int polydeg, const double amin, 
+                    const double amax, const int n_kr, const int solver_flag, const int rel_prec,
                     const int even_odd_flag, const SloppyPrecision refinement_precision,
                     SloppyPrecision sloppy_precision, CompressionType compression) {
 
@@ -2646,6 +2647,7 @@ double eigsolveQuda(int n, double tol, int blksize, int blkwise, int max_iterati
   QudaInvertParam eig_invert_param = newQudaInvertParam();
   eig_invert_param = inv_param;
   eig_param.invert_param = &eig_invert_param;
+  eig_param.invert_param->verbosity = QUDA_VERBOSE;
   /* AS The following two are set to cuda_prec, otherwise 
    * it gives an error. Such high precision might not be
    * necessary. But have not found a way to consistently set
@@ -2689,10 +2691,10 @@ double eigsolveQuda(int n, double tol, int blksize, int blkwise, int max_iterati
   // at a relatively coarse lattice spacing for the eigenvalues
   // of the twisted-clover ND operator with values of musigma / mudelta
   // reproducing physical sea strange and charm quark masses
-  eig_param.use_poly_acc = maxmin == 1 ? QUDA_BOOLEAN_FALSE : QUDA_BOOLEAN_TRUE;
-  eig_param.poly_deg = 128;
-  eig_param.a_min = 1e-3;
-  eig_param.a_max = 4;
+  eig_param.use_poly_acc = (maxmin == 1) && (polydeg != 0) ? QUDA_BOOLEAN_FALSE : QUDA_BOOLEAN_TRUE;
+  eig_param.poly_deg = polydeg;
+  eig_param.a_min = amin;
+  eig_param.a_max = amax;
 
   /* Daggers the operator. Not necessary for 
    * most cases. */
@@ -2730,7 +2732,7 @@ double eigsolveQuda(int n, double tol, int blksize, int blkwise, int max_iterati
    * From my understanding, QUDA automatically scales
    * this search space, however more testing on this 
    * might be necessary */
-  eig_param.n_kr = 96;
+  eig_param.n_kr = n_kr;
 
   eig_param.max_restarts = max_iterations;
 

diff --git a/quda_interface.h b/quda_interface.h
@@ -176,7 +176,8 @@ void compute_WFlow_quda(const double eps ,const double tmax, const int traj, FIL
 
 
 double eigsolveQuda(int n, double tol, int blksize, int blkwise, int max_iterations, int maxmin,
-                    const double precision, const int max_iter, const int solver_flag, const int rel_prec,
+                    const double precision, const int max_iter, const int polydeg, const double amin, 
+                    const double amax, const int n_kr, const int solver_flag, const int rel_prec,
                     const int even_odd_flag, const SloppyPrecision refinement_precision,
                     SloppyPrecision sloppy_precision, CompressionType compression);
 

diff --git a/read_input.l b/read_input.l
@@ -2607,6 +2607,26 @@ static inline double fltlist_next_token(int * const list_end){
     if(myverbose) printf("  Do not use external eigensolver line %d monomial %d\n", line_of_file, current_monomial);
     mnl->external_eigsolver = NO_EXT_EIGSOLVER;
   }
+  {SPC}*EigAmin{EQL}{FLT} {
+    sscanf(yytext, " %[a-zA-Z] = %lf", name, &c);
+    mnl->eig_amin = c;
+    if(myverbose!=0) printf("  eig_amin set to %e line %d monomial %d\n", c, line_of_file, current_monomial);
+  }
+  {SPC}*EigAmax{EQL}{FLT} {
+    sscanf(yytext, " %[a-zA-Z] = %lf", name, &c);
+    mnl->eig_amax = c;
+    if(myverbose!=0) printf("  eig_amax set to %e line %d monomial %d\n", c, line_of_file, current_monomial);
+  }
+  {SPC}*EigPolyDeg{EQL}{DIGIT}+ {
+    sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+    mnl->eig_polydeg = a;
+    if(myverbose!=0) printf("  eig_polydeg set to %d line %d monomial %d\n", a, line_of_file, current_monomial);
+  }
+  {SPC}*EigNkr{EQL}{DIGIT}+ {
+    sscanf(yytext, " %[a-zA-Z] = %d", name, &a);
+    mnl->eig_n_kr = a;
+    if(myverbose!=0) printf("  eig_n_kr set to %d line %d monomial %d\n", a, line_of_file, current_monomial);
+  }
 }
 <NDPOLYMONOMIAL,CLPOLYMONOMIAL>{
   {SPC}*MaxPtildeDegree{EQL}{DIGIT}+ {