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

5 files changed, 59 insertions, 11 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));