diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/cluster.h | 12 | ||||
| -rw-r--r-- | lib/correlation.h | 14 | ||||
| -rw-r--r-- | lib/measure.h | 4 | ||||
| -rw-r--r-- | lib/state.h | 16 | ||||
| -rw-r--r-- | lib/vector.h | 24 |
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)); |