efficient computation of the smallest fourier mode by doing a magnetization-style update anytime a bond with the external field changes

6 files changed, 78 insertions, 40 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 8061d34..2225e2f 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -94,12 +94,15 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
 
         if (is_ext) {
           X_t rs_old, rs_new;
+          v_t non_ghost;
           if (vn == state->nv) {
             rs_old = act_inverse (R_old, si_old);
             rs_new = act_inverse (R_old, si_new);
+            non_ghost = v;
           } else {
             rs_old = act_inverse (R_old, sj);
             rs_new = act_inverse (R_new, sj);
+            non_ghost = vn;
           }
           double dE = state->H(rs_old) - state->H(rs_new);
           prob = 1.0 - exp(-dE / state->T);
@@ -108,6 +111,15 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           add (&(state->M), rs_new);
           state->E += dE;
 
+          for (D_t i = 0; i < state->D; i++) {
+            double x = (double)((non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L) / (double)state->L; 
+            scalar_subtract(&(state->ReF[i]), cos(2 * M_PI * x), rs_old);
+            scalar_add(&(state->ReF[i]), cos(2 * M_PI * x), rs_new);
+
+            scalar_subtract(&(state->ImF[i]), sin(2 * M_PI * x), rs_old);
+            scalar_add(&(state->ImF[i]), sin(2 * M_PI * x), rs_new);
+          }
+
           free_spin (rs_old);
           free_spin (rs_new);
         } else {
diff --git a/lib/correlation.h b/lib/correlation.h
index b3b7d86..26c3a99 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -7,15 +7,13 @@
 #include <fftw3.h>
 
 template <class R_t, class X_t>
-double correlation_length(const state_t <R_t, X_t> *s, fftw_plan plan, double *in, double *out) {
-  for (v_t i = 0; i < s->nv; i++) {
-    in[i] = correlation_component(s->spins[i]);
-  }
-
-  fftw_execute(plan);
+double correlation_length(const state_t <R_t, X_t> *s) {
+  double total = 0;
 
-  double length = pow(out[0], 2);
+  for (D_t j = 0; j < s->D; j++) {
+    total += norm_squared(s->ReF[j]) + norm_squared(s->ImF[j]);
+  }
 
-  return length;
+  return total / s->D;
 }
 
diff --git a/lib/measure.h b/lib/measure.h
index 3474684..d20c081 100644
--- a/lib/measure.h
+++ b/lib/measure.h
@@ -36,9 +36,9 @@ std::function <void(const state_t <R_t, X_t> *)> measurement_magnetization_file(
 }
 
 template <class R_t, class X_t>
-std::function <void(const state_t <R_t, X_t> *)> measurement_fourier_file(FILE *file, fftw_plan plan, double *fftw_in, double *fftw_out) {
+std::function <void(const state_t <R_t, X_t> *)> measurement_fourier_file(FILE *file) {
   return [=](const state_t <R_t, X_t> *s) {
-    float smaller_X = (float)correlation_length(s, plan, fftw_in, fftw_out);
+    float smaller_X = (float)correlation_length(s);
     fwrite(&smaller_X, sizeof(float), 1, file);
   };
 }
diff --git a/lib/state.h b/lib/state.h
index 1491938..c70459f 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -20,6 +20,8 @@ class state_t {
     double E;
     X_t M; // the "sum" of the spins, like the total magnetization
     v_t last_cluster_size;
+    X_t *ReF;
+    X_t *ImF;
 
     std::function <double(X_t, X_t)> J;
     std::function <double(X_t)> H;
@@ -36,8 +38,14 @@ class state_t {
       }
       init (&R);
       E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
-      M = scalar_multiple (nv, spins[0]);
+      M = scalar_multiple(nv, spins[0]);
       last_cluster_size = 0;
+      ReF = (X_t *)malloc(D * sizeof(X_t));
+      ImF = (X_t *)malloc(D * sizeof(X_t));
+      for (D_t i = 0; i < D; i++) {
+        ReF[i] = scalar_multiple(0, spins[0]);
+        ImF[i] = scalar_multiple(0, spins[0]);
+      }
     }
 
     ~state_t() {
@@ -48,6 +56,12 @@ class state_t {
       free(spins);
       free_spin(R);
       free_spin(M);
+      for (D_t i = 0; i < D; i++) {
+        free_spin(ReF[i]);
+        free_spin(ImF[i]);
+      }
+      free(ReF);
+      free(ImF);
     }
 };
 
diff --git a/lib/vector.h b/lib/vector.h
index 2099d28..a1084e2 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -37,6 +37,20 @@ void add (vector_t <q, T> *v1, vector_t <q, T> v2) {
 }
 
 template <q_t q, class T>
+void scalar_add (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
+  for (q_t i = 0; i < q; i++) {
+    v1->x[i] += a * v2.x[i];
+  }
+}
+
+template <q_t q, class T>
+void scalar_subtract (vector_t <q, T> *v1, T a, vector_t <q, T> v2) {
+  for (q_t i = 0; i < q; i++) {
+    v1->x[i] -= a * v2.x[i];
+  }
+}
+
+template <q_t q, class T>
 void subtract (vector_t <q, T> *v1, vector_t <q, T> v2) {
   for (q_t i = 0; i < q; i++) {
     v1->x[i] -= v2.x[i];
@@ -44,6 +58,16 @@ void subtract (vector_t <q, T> *v1, vector_t <q, T> v2) {
 }
 
 template <q_t q, class T>
+T norm_squared (vector_t <q, T> v) {
+  T tmp = 0;
+  for (q_t i = 0; i < q; i++) {
+    tmp += v.x[i] * v.x[i];
+  }
+
+  return tmp;
+}
+
+template <q_t q, class T>
 vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
   vector_t <q, T> multiple;
   multiple.x = (T *)malloc(q * sizeof(T));
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 329f468..d718440 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -115,24 +115,33 @@ int main(int argc, char *argv[]) {
 
   FILE *outfile_info = fopen("wolff_metadata.txt", "a");
 
-  fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"%s\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, ON_strings[N_COMP], N_COMP, D, L, pow(L, D), D * pow(L, D), T);
-
-  for (q_t i = 0; i < N_COMP; i++) {
-    fprintf(outfile_info, "%.15f", H_vec[i]);
-    if (i < N_COMP - 1) {
-      fprintf(outfile_info, ", ");
+  fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"%s\", \"q\" -> %d, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"FIELD_TYPE\" -> ", timestamp, ON_strings[N_COMP], N_COMP, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T);
+  if (modulated_field) {
+    fprintf(outfile_info, "\"MODULATED\", \"ORDER\" -> %d, \"H\" -> %.15f, ", order, H_vec[0]);
+  } else {
+    fprintf(outfile_info, "\"VECTOR\", \"H\" -> {");
+    for (q_t i = 0; i < N_COMP; i++) {
+      fprintf(outfile_info, "%.15f", H_vec[i]);
+      if (i < N_COMP - 1) {
+        fprintf(outfile_info, ", ");
+      }
     }
+    fprintf(outfile_info, "}, ");
   }
 
-  fprintf(outfile_info, "}, \"GENERATOR\" -> \"%s\", \"EPS\" -> %g |>\n", pert_type, epsilon);
+  fprintf(outfile_info, "\"GENERATOR\" -> \"%s\"", pert_type);
+
+  if (use_pert) {
+    fprintf(outfile_info, ", \"EPS\" -> %g", epsilon);
+  }
+
+  fprintf(outfile_info, " |>\n");
 
   fclose(outfile_info);
 
   unsigned int n_measurements = 0;
   std::function <void(const On_t *)> *measurements = (std::function <void(const On_t *)> *)calloc(POSSIBLE_MEASUREMENTS, sizeof(std::function <void(const On_t *)>));
   FILE *outfile_M, *outfile_E, *outfile_S, *outfile_F;
-  double *fftw_in, *fftw_out;
-  fftw_plan plan;
 
   if (measurement_flags & measurement_energy) {
     char *filename_E = (char *)malloc(255 * sizeof(char));
@@ -166,22 +175,7 @@ int main(int argc, char *argv[]) {
     sprintf(filename_F, "wolff_%lu_F.dat", timestamp);
     outfile_F = fopen(filename_F, "wb");
     free(filename_F);
-
-    fftw_in = (double *)fftw_malloc(pow(L, D) * sizeof(double));
-    fftw_out = (double *)fftw_malloc(pow(L, D) * sizeof(double));
-    int rank = D;
-    int *n = (int *)malloc(rank * sizeof(int));
-    fftw_r2r_kind *kind = (fftw_r2r_kind *)malloc(rank * sizeof(fftw_r2r_kind));
-    for (D_t i = 0; i < rank; i++) {
-      n[i] = L;
-      kind[i] = FFTW_R2HC;
-    }
-    plan = fftw_plan_r2r(rank, n, fftw_in, fftw_out, kind, 0);
-
-    free(n);
-    free(kind);
-
-    measurements[n_measurements] = measurement_fourier_file<orthogonal_R_t, vector_R_t> (outfile_F, plan, fftw_in, fftw_out);
+    measurements[n_measurements] = measurement_fourier_file<orthogonal_R_t, vector_R_t> (outfile_F);
     n_measurements++;
   }
 
@@ -208,10 +202,6 @@ int main(int argc, char *argv[]) {
   }
   if (measurement_flags & measurement_fourierZero) {
     fclose(outfile_F);
-    fftw_destroy_plan(plan);
-    fftw_free(fftw_in);
-    fftw_free(fftw_out);
-    fftw_cleanup(); // fftw is only used if fourier modes are measured!
   }
 
   free(H_vec);