From 72301b3d5c3a91ff2e7fc6eedcad7bce8e647efa Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 19 Jul 2018 18:22:15 -0400
Subject: efficient computation of the smallest fourier mode by doing a
magnetization-style update anytime a bond with the external field changes
---
lib/cluster.h | 12 ++++++++++++
lib/correlation.h | 14 ++++++--------
lib/measure.h | 4 ++--
lib/state.h | 16 +++++++++++++++-
lib/vector.h | 24 ++++++++++++++++++++++++
5 files changed, 59 insertions(+), 11 deletions(-)
(limited to 'lib')
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
@@ -36,6 +36,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++) {
@@ -43,6 +57,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;
--
cgit v1.2.3-70-g09d2