restructuring of gauge observables

Include/Core/global.h

@@ -98,6 +98,7 @@ GLB_VAR(int, *iup, = NULL);
 GLB_VAR(int, *idn, = NULL);
 GLB_VAR(char, *imask, = NULL);
 GLB_VAR(int, zerocoord[4], = { 0, 0, 0, 0 });
+GLB_VAR(int, *timeslices, = NULL);
 
 GLB_VAR(int, BLK_T, = 4);
 GLB_VAR(int, BLK_X, = 4);
@@ -154,6 +155,7 @@ GLB_VAR(int, *iup_gpu, = NULL);
 GLB_VAR(int, *idn_gpu, = NULL);
 GLB_VAR(char, *imask_gpu, = NULL);
 GLB_VAR(unsigned int, grid_size_max_gpu, = 65535);
+GLB_VAR(int, *timeslices_gpu, = NULL);
 #else
 #define STD_MEM_TYPE (CPU_MEM)
 #endif

Include/Update/gauge_observables.h

@@ -0,0 +1,58 @@
+/**
+ * @file avr_plaquette.h
+ * @brief Plaquette evaluation functions
+ */
+
+#ifndef AVR_PLAQUETTE_H
+#define AVR_PLAQUETTE_H
+
+#include "hr_complex.h"
+#include "Core/spinor_field.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// Local functions
+double plaq(suNg_field *gauge, int ix, int mu, int nu);
+void cplaq(hr_complex *ret, suNg_field *gauge, int ix, int mu, int nu);
+void clover_F(suNg_algebra_vector *F, suNg_field *V, int ix, int mu, int nu);
+double local_plaq(suNg_field *gauge, int ix);
+
+// Observables
+extern double (*avr_plaquette)();
+extern double (*avr_plaquette_suNg_field)(suNg_field *gauge);
+extern void (*avr_plaquette_time)(suNg_field *gauge, double *plaqt, double *plaqs);
+extern void (*full_plaquette)();
+extern void (*full_plaquette_suNg_field)(suNg_field *gauge);
+extern void (*local_plaquette)(suNg_field *gauge, scalar_field *s);
+extern double (*E)(suNg_field *V);
+extern void (*E_T)(double *E, suNg_field *V);
+extern double (*Esym)(suNg_field *V);
+extern void (*Esym_T)(double *E, suNg_field *V);
+extern double (*topo)(suNg_field *V);
+
+// Mostly necessary for testing
+double avr_plaquette_cpu();
+double avr_plaquette_suNg_field_cpu(suNg_field *u);
+void avr_plaquette_time_cpu(suNg_field *u, double *plaqt, double *plaqs);
+void full_plaquette_cpu();
+void full_plaquette_suNg_field_cpu(suNg_field *u);
+
+#ifdef WITH_GPU
+double avr_plaquette_gpu();
+double avr_plaquette_suNg_field_gpu(suNg_field *u);
+void avr_plaquette_time_gpu(suNg_field *u, double *plaqt, double *plaqs);
+void full_plaquette_gpu();
+void full_plaquette_suNg_field_gpu(suNg_field *u);
+#endif
+
+// Workspace functions
+void avr_ts_plaquette(void);
+void cplaq_wrk(hr_complex *ret, int ix, int mu, int nu);
+hr_complex avr_plaquette_wrk(void);
+
+#ifdef __cplusplus
+}
+#endif
+#endif

Include/update.h

@@ -6,7 +6,6 @@
 #ifndef UPDATE_H
 #define UPDATE_H
 
-#include "Update/avr_plaquette.h"
 #include "Update/cabmar.h"
 #include "Update/clover_tools.h"
 #include "Update/dirac.h"
@@ -23,5 +22,6 @@
 #include "Update/stout_smearing.h"
 #include "Update/update_algo.h"
 #include "Update/twisted_bc.h"
+#include "Update/gauge_observables.h"
 
 #endif

LibHR/Geometry/geometry_init_gpu.c

@@ -15,11 +15,13 @@ void init_neighbors_gpu() {
 
     CHECK_CUDA(cudaMalloc((void **)&iup_gpu, 4 * N * sizeof(int)));
     CHECK_CUDA(cudaMalloc((void **)&idn_gpu, 4 * N * sizeof(int)));
+    CHECK_CUDA(cudaMalloc((void **)&timeslices_gpu, N * sizeof(int)));
     CHECK_CUDA(cudaMalloc((void **)&imask_gpu, N * sizeof(char)));
     CHECK_CUDA(cudaMalloc((void **)&ipt_gpu,
                           (X + 2 * X_BORDER) * (Y + 2 * Y_BORDER) * (Z + 2 * Z_BORDER) * (T + 2 * T_BORDER) * sizeof(int)));
     CHECK_CUDA(cudaMemcpy(iup_gpu, iup, 4 * N * sizeof(int), cudaMemcpyHostToDevice));
     CHECK_CUDA(cudaMemcpy(idn_gpu, idn, 4 * N * sizeof(int), cudaMemcpyHostToDevice));
+    CHECK_CUDA(cudaMemcpy(timeslices_gpu, timeslices, N * sizeof(int), cudaMemcpyHostToDevice));
     CHECK_CUDA(cudaMemcpy(imask_gpu, imask, N * sizeof(*imask), cudaMemcpyHostToDevice));
     CHECK_CUDA(cudaMemcpy(ipt_gpu, ipt,
                           (X + 2 * X_BORDER) * (Y + 2 * Y_BORDER) * (Z + 2 * Z_BORDER) * (T + 2 * T_BORDER) * sizeof(int),

LibHR/Geometry/new_geom.c

@@ -410,6 +410,10 @@ static void index_alloc() {
     icoord = malloc((main_mem_volume + buf_mem_volume) * sizeof(*icoord));
     error((icoord == NULL), 1, __func__, "Cannot allocate memory for send/recv icoord");
     sb_icoord = icoord + main_mem_volume;
+
+    // memory for timeslice array
+    timeslices = malloc((main_mem_volume + buf_mem_volume) * sizeof(*icoord));
+    error((icoord == NULL), 1, __func__, "Cannot allocate memory for time slice array");
 }
 
 static void index_free() {
@@ -523,6 +527,7 @@ static void enumerate_lattice() {
     const int INVALID_POINT = -1;
     //NB: memset only works with 0 and -1 (assuming 2's complement numbers)
     memset(ipt, INVALID_POINT, T_EXT * X_EXT * Y_EXT * Z_EXT * sizeof(*ipt));
+    memset(timeslices, INVALID_POINT, T_EXT * X_EXT * Y_EXT * Z_EXT * sizeof(*timeslices));
 
     //enumerate points in each box
     //we assign an id to each point in a box and record it in ipt_ext and icoord
@@ -602,6 +607,18 @@ static void enumerate_lattice() {
         }
     }
 
+    for (int nt = 0; nt < T; nt++) {
+        for (int nx = 0; nx < X; nx++) {
+            for (int ny = 0; ny < Y; ny++) {
+                for (int nz = 0; nz < Z; nz++) {
+                    int ix = ipt(nt, nx, ny, nz);
+                    // Initialize the array that returns the local time slice for a given index
+                    timeslices[ix] = nt;
+                }
+            }
+        }
+    }
+
 #ifdef _OPENMP
     //Correct for NUMA affinity
     //we create a new copy of the arrays: iup, idn, imask and icoord

LibHR/Update/Dphi.c

@@ -608,7 +608,9 @@ void Dphi_cpu_(spinor_field *restrict out, spinor_field *restrict in) {
 #endif
 
         _SITE_FOR(out->type, ixp, ix) {
+#ifdef _OPENMP
             register int thread0 = hr_threadId();
+#endif
             int iy;
             suNf *up, *um;
             suNf_vector psi, chi, psi2, chi2;