Merge pull request #577 from etmc/compare_derivative

move compare_derivative into its own compilation unit and make the strictness variable, specify a 'name' argument

Makefile.in

@@ -58,7 +58,7 @@ MODULES = read_input gamma measure_gauge_action start \
 	little_D block operator \
 	spinor_fft X_psi P_M_eta \
 	jacobi fatal_error invert_clover_eo gettime \
-	tm_debug_printf \
+	tm_debug_printf compare_derivative \
         @SPI_FILES@ @QUDA_INTERFACE@ @DDalphaAMG_INTERFACE@
 
 CXXMODULES = @QPHIX_INTERFACE@

compare_derivative.c

@@ -0,0 +1,67 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2024 Bartosz Kostrzewa
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * tmLQCD is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ ***********************************************************************/
+
+#ifdef HAVE_CONFIG_H
+# include<tmlqcd_config.h>
+#endif
+#ifdef TM_USE_OMP 
+# include <omp.h>
+#endif
+#include <stdio.h>
+#include "global.h"
+#include "monomial/monomial.h"
+
+/* this function compares two derivatives calculated by an external library and tmLQCD */
+void compare_derivative(monomial *mnl, su3adj **ext_lib, su3adj **native, 
+    const double threshold, const char * name){
+  int n_diff = 0;
+
+  for(int ix = 0; ix < VOLUME; ix++){
+    for(int mu=0; mu<4; mu++){
+      double *ext=&(ext_lib[ix][mu].d1);
+      double *nat=&(native[ix][mu].d1);
+      for(int j=0; j<8; ++j){
+        double diff=ext[j]-nat[j];
+        if (sqrt(diff*diff) > threshold || isnan( ext[j] ) || isinf(ext[j]) ){
+            n_diff++;
+            printf("derivative at (t,x,y,z,mu,j) %d,%d,%d,%d,%d,%d,"
+                   " ext: %-14e, native: %-14e ratio: %-14g diff %-14g  on proc_id %d\n", 
+                  g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu, j,
+                  ext[j], nat[j], ext[j]/nat[j], ext[j]-nat[j], g_proc_id);
+        }
+      }
+    }
+  }
+  if(n_diff > 0){
+    printf("%s: the deviation between tmLQCD and the external library "
+           "exceeds the threshold %.1e in %d case(s) for parameters: c0=%e c1=%e g_beta=%e on proc_id: %d\n",
+           name,
+           threshold,
+           n_diff,
+           mnl->c0,
+           mnl->c1,
+           mnl->beta,
+           g_proc_id);
+
+    if(g_strict_residual_check) fatal_error("Difference between external library and tmLQCD-native function!", 
+                                            name);
+  }
+}
+

compare_derivative.h

@@ -0,0 +1,29 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2024 Bartosz Kostrzewa
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * tmLQCD is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ ***********************************************************************/
+
+#ifndef COMPARE_DERIVATIVE_H
+#define COMPARE_DERIVATIVE_H
+
+#include "monomial/monomial.h"
+#include "su3adj.h"
+
+void compare_derivative(monomial *mnl, su3adj **ext_lib, su3adj **native, const double threshold, const char * name);
+
+#endif

monomial/cloverdet_monomial.c

@@ -51,7 +51,7 @@
 #include "operator/clovertm_operators_32.h"
 #include "cloverdet_monomial.h"
 #include "xchange/xchange_deri.h"
-#include "monomial/gauge_monomial.h"
+#include "compare_derivative.h"
 #ifdef TM_USE_QUDA
 #  include "quda_interface.h"
 #endif
@@ -143,7 +143,7 @@ void cloverdet_derivative(const int id, hamiltonian_field_t * const hf) {
       #ifdef TM_USE_MPI
         xchange_deri(hf->derivative);// this function use ddummy inside
       #endif
-      compare_derivative(mnl, given, hf->derivative);
+      compare_derivative(mnl, given, hf->derivative, 1e-9, "cloverdet_derivative");
       mnl->external_library = QUDA_LIB;
       hf->derivative = given;
     }

monomial/cloverdetratio_monomial.c

@@ -48,6 +48,7 @@
 #include "cloverdetratio_monomial.h"
 #include "xchange/xchange_deri.h"
 #include "monomial/gauge_monomial.h"
+#include "compare_derivative.h"
 #ifdef TM_USE_QUDA
 #  include "quda_interface.h"
 #endif
@@ -253,7 +254,7 @@ void cloverdetratio_derivative(const int no, hamiltonian_field_t * const hf) {
       #ifdef TM_USE_MPI
         xchange_deri(hf->derivative);// this function use ddummy inside
       #endif
-      compare_derivative(mnl, given, hf->derivative);
+      compare_derivative(mnl, given, hf->derivative, 1e-9, "cloverdetratio_derivative");
       mnl->external_library = QUDA_LIB;
       hf->derivative = given;
     }

monomial/gauge_monomial.c

@@ -44,49 +44,12 @@
 #include "monomial/monomial.h"
 #include "hamiltonian_field.h"
 #include "gauge_monomial.h"
+#include "compare_derivative.h"
 #include "fatal_error.h"
 #ifdef TM_USE_QUDA
 #include "quda_interface.h"
 #endif
 
-/* this function compares the gauge derivative calculated by an external library and tmLQCD */
-void compare_derivative(monomial *mnl, su3adj **ext_lib, su3adj **native ){
-  const double threshold = 1e-7;
-  int n_diff = 0;
-
-  for(int ix = 0; ix < VOLUME; ix++){
-    for(int mu=0; mu<4; mu++){
-      double *ext=&(ext_lib[ix][mu].d1);
-      double *nat=&(native[ix][mu].d1);
-      for(int j=0; j<8; ++j){
-        double diff;
-        if (fabs(nat[j])>1e-7) diff=(ext[j]-nat[j])/nat[j];
-        else diff=ext[j]-nat[j];
-        if (sqrt(diff*diff) > threshold || isnan( ext[j] ) || isinf(ext[j]) ){
-            n_diff++;
-            printf("derivative at (t,x,y,z,mu,j) %d,%d,%d,%d,%d,%d, ext: %-14e, native: %-14e ratio: %-14g diff %-14g  on proc_id %d\n", 
-                  g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu, j,
-                  ext[j], nat[j], ext[j]/nat[j], ext[j]-nat[j], g_proc_id);
-        }
-      }
-    }
-  }
-  if(n_diff > 0){
-    printf("gauge_derivative: the relative deviation between tmLQCD and the external library "
-           "exceeds the threshold %.1e in %d case(s) for parameters: c0=%e c1=%e g_beta=%e on proc_id: %d\n",
-           threshold,
-           n_diff,
-           mnl->c0,
-           mnl->c1,
-           mnl->beta,
-           g_proc_id);
-
-    if(g_strict_residual_check) fatal_error("Difference between external library and tmLQCD-native function!", 
-                                            "gauge_derivative");
-  }
-}
-
-
 /* this function calculates the derivative of the momenta: equation 13 of Gottlieb */
 void gauge_derivative(const int id, hamiltonian_field_t * const hf) {
   monomial * mnl = &monomial_list[id];
@@ -107,7 +70,7 @@ void gauge_derivative(const int id, hamiltonian_field_t * const hf) {
      hf->derivative=ddummy;
      mnl->external_library=NO_EXT_LIB;
      gauge_derivative(id, hf);
-     compare_derivative(mnl,given, ddummy);
+     compare_derivative(mnl, given, ddummy, 1e-9, "gauge_derivative");
      mnl->external_library=QUDA_LIB;
      hf->derivative=given;
     }

monomial/gauge_monomial.h

@@ -26,6 +26,5 @@ void gauge_derivative(const int id, hamiltonian_field_t * const hf);
 void gauge_EMderivative(const int id, hamiltonian_field_t * const hf);
 void gauge_heatbath(const int id, hamiltonian_field_t * const hf);
 double gauge_acc(const int id, hamiltonian_field_t * const hf);
-void compare_derivative(monomial *mnl, su3adj **ext_lib, su3adj **native );
 
 #endif