diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/cluster.h | 65 | ||||
-rw-r--r-- | lib/correlation.h | 8 | ||||
-rw-r--r-- | lib/measure.h | 10 | ||||
-rw-r--r-- | lib/wolff.h | 12 | ||||
-rw-r--r-- | lib/z2.h | 8 |
5 files changed, 53 insertions, 50 deletions
diff --git a/lib/cluster.h b/lib/cluster.h index f948586..f3271b0 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -12,45 +12,45 @@ #include "graph.h" template <class R_t, class X_t> -void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { +void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) { v_t nv = 0; std::stack<v_t> stack; stack.push(v0); - std::vector<bool> marks(state->g.nv, false); + std::vector<bool> marks(s.g.nv, false); while (!stack.empty()) { v_t v = stack.top(); stack.pop(); - if (!marks[v]) { + if (!marks[v]) { // don't reprocess anyone we've already visited! X_t si_old, si_new; R_t R_old, R_new; - R_old = state->R; + R_old = s.R; marks[v] = true; - bool v_is_ghost = (v == state->nv); + bool v_is_ghost = (v == s.nv); // ghost site has the last index if (v_is_ghost) { - R_new = r.act(R_old); + R_new = r.act(R_old); // if we are, then we're moving the transformation } else { - si_old = state->spins[v]; - si_new = r.act(si_old); + si_old = s.spins[v]; + si_new = r.act(si_old); // otherwise, we're moving the spin at our site } - for (const v_t &vn : state->g.v_adj[v]) { + for (const v_t &vn : s.g.v_adj[v]) { X_t sj; - bool vn_is_ghost = (vn == state->nv); + bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost if (!vn_is_ghost) { - sj = state->spins[vn]; + sj = s.spins[vn]; } double prob; - if (v_is_ghost || vn_is_ghost) { + if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond... X_t rs_old, rs_new; v_t non_ghost; if (vn_is_ghost) { @@ -63,51 +63,54 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { non_ghost = vn; } - double dE = state->H(rs_old) - state->H(rs_new); + double dE = s.H(rs_old) - s.H(rs_new); + #ifdef FINITE_STATES prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)]; #else - prob = 1.0 - exp(-dE / state->T); + prob = 1.0 - exp(-dE / s.T); #endif - state->M -= rs_old; - state->M += rs_new; + s.M -= rs_old; + s.M += rs_new; - state->E += dE; + s.E += dE; - for (D_t i = 0; i < state->D; i++) { - L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L; + for (D_t i = 0; i < s.D; i++) { + L_t x = (non_ghost / (v_t)pow(s.L, s.D - i - 1)) % s.L; - state->ReF[i] -= rs_old * state->precomputed_cos[x]; - state->ReF[i] += rs_new * state->precomputed_cos[x]; + s.ReF[i] -= rs_old * s.precomputed_cos[x]; + s.ReF[i] += rs_new * s.precomputed_cos[x]; - state->ImF[i] -= rs_old * state->precomputed_sin[x]; - state->ImF[i] += rs_new * state->precomputed_sin[x]; + s.ImF[i] -= rs_old * s.precomputed_sin[x]; + s.ImF[i] += rs_new * s.precomputed_sin[x]; } - } else { - double dE = state->J(si_old, sj) - state->J(si_new, sj); + } else { // otherwise, we're at a perfectly normal bond! + double dE = s.J(si_old, sj) - s.J(si_new, sj); + #ifdef FINITE_STATES prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)]; #else - prob = 1.0 - exp(-dE / state->T); + prob = 1.0 - exp(-dE / s.T); #endif - state->E += dE; + + s.E += dE; } - if (gsl_rng_uniform(rand) < prob) { // and with probability... + if (gsl_rng_uniform(rand) < prob) { stack.push(vn); // push the neighboring vertex to the stack } } if (v_is_ghost) { - state->R = R_new; + s.R = R_new; } else { - state->spins[v] = si_new; + s.spins[v] = si_new; nv++; } } } - state->last_cluster_size = nv; + s.last_cluster_size = nv; } diff --git a/lib/correlation.h b/lib/correlation.h index 49c6ff2..23b1bd9 100644 --- a/lib/correlation.h +++ b/lib/correlation.h @@ -7,13 +7,13 @@ #include <fftw3.h> template <class R_t, class X_t> -double correlation_length(const state_t <R_t, X_t> *s) { +double correlation_length(const state_t <R_t, X_t>& s) { double total = 0; - for (D_t j = 0; j < s->D; j++) { - total += norm_squared(s->ReF[j]) + norm_squared(s->ReF[j]); + for (D_t j = 0; j < s.D; j++) { + total += norm_squared(s.ReF[j]) + norm_squared(s.ReF[j]); } - return total / s->D; + return total / s.D; } diff --git a/lib/measure.h b/lib/measure.h index 8082dd2..2c5ffb7 100644 --- a/lib/measure.h +++ b/lib/measure.h @@ -29,17 +29,17 @@ FILE **measure_setup_files(unsigned char flags, unsigned long timestamp) { } template <class R_t, class X_t> -std::function <void(const state_t <R_t, X_t> *)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t> *)> other_f) { - return [=] (const state_t <R_t, X_t> *s) { +std::function <void(const state_t <R_t, X_t>&)> measure_function_write_files(unsigned char flags, FILE **files, std::function <void(const state_t <R_t, X_t>&)> other_f) { + return [=] (const state_t <R_t, X_t>& s) { if (flags & measurement_energy) { - float smaller_E = (float)s->E; + float smaller_E = (float)s.E; fwrite(&smaller_E, sizeof(float), 1, files[0]); } if (flags & measurement_clusterSize) { - fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, files[1]); + fwrite(&(s.last_cluster_size), sizeof(uint32_t), 1, files[1]); } if (flags & measurement_magnetization) { - write_magnetization(s->M, files[2]); + write_magnetization(s.M, files[2]); } if (flags & measurement_fourierZero) { float smaller_X = (float)correlation_length(s); diff --git a/lib/wolff.h b/lib/wolff.h index a4a663c..498f7f3 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -3,18 +3,18 @@ #include "state.h" template <class R_t, class X_t> -void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s0)> gen_R, std::function <void(const state_t <R_t, X_t> *)> measurements, gsl_rng *r, bool silent) { +void wolff(count_t N, state_t <R_t, X_t>& s, std::function <R_t(gsl_rng *, X_t s0)> gen_R, std::function <void(const state_t <R_t, X_t>&)> measurements, gsl_rng *r, bool silent) { #ifdef FINITE_STATES - initialize_probs(s->J, s->H, s->T); + initialize_probs(s.J, s.H, s.T); #endif if (!silent) printf("\n"); for (count_t steps = 0; steps < N; steps++) { - if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s->E, s->last_cluster_size); + if (!silent) printf("\033[F\033[JWOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, s.last_cluster_size); - v_t v0 = gsl_rng_uniform_int(r, s->nv); - R_t step = gen_R(r, s->spins[v0]); + v_t v0 = gsl_rng_uniform_int(r, s.nv); + R_t step = gen_R(r, s.spins[v0]); flip_cluster <R_t, X_t> (s, v0, step, r); measurements(s); @@ -23,7 +23,7 @@ void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s if (!silent) { printf("\033[F\033[J"); } - printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s->E, s->last_cluster_size); + printf("WOLFF: step %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size); } @@ -25,7 +25,7 @@ class z2_t { z2_t(bool x) : x(x) {} - ising_t act(const ising_t& s) { + ising_t act(const ising_t& s) const { if (x) { return ising_t(!s.x); } else { @@ -33,7 +33,7 @@ class z2_t { } } - z2_t act(const z2_t& r) { + z2_t act(const z2_t& r) const { if (x) { return z2_t(!r.x); } else { @@ -41,11 +41,11 @@ class z2_t { } } - ising_t act_inverse(const ising_t& s) { + ising_t act_inverse(const ising_t& s) const { return this->act(s); } - z2_t act_inverse(const z2_t& r) { + z2_t act_inverse(const z2_t& r) const { return this->act(r); } }; |