diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/cluster.h | 17 | ||||
| -rw-r--r-- | lib/correlation.h | 2 | ||||
| -rw-r--r-- | lib/ising.h | 65 | ||||
| -rw-r--r-- | lib/orthogonal.h | 27 | ||||
| -rw-r--r-- | lib/state.h | 25 | ||||
| -rw-r--r-- | lib/vector.h | 109 |
6 files changed, 103 insertions, 142 deletions
diff --git a/lib/cluster.h b/lib/cluster.h index 8aae02a..427c3c8 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -68,22 +68,20 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { prob = 1.0 - exp(-dE / state->T); #endif - add(&(state->M), -1, rs_old); - add(&(state->M), 1, rs_new); + state->M -= rs_old; + state->M += rs_new; + state->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; - add(&(state->ReF[i]), -state->precomputed_cos[x], rs_old); - add(&(state->ReF[i]), state->precomputed_cos[x], rs_new); + state->ReF[i] -= rs_old * state->precomputed_cos[x]; + state->ReF[i] += rs_new * state->precomputed_cos[x]; - add(&(state->ImF[i]), -state->precomputed_sin[x], rs_old); - add(&(state->ImF[i]), state->precomputed_sin[x], rs_new); + state->ImF[i] -= rs_old * state->precomputed_sin[x]; + state->ImF[i] += rs_new * state->precomputed_sin[x]; } - - free_spin (rs_old); - free_spin (rs_new); } else { double dE = state->J(si_old, sj) - state->J(si_new, sj); #ifdef FINITE_STATES @@ -103,7 +101,6 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { free_spin(state->R); state->R = R_new; } else { - free_spin(state->spins[v]); state->spins[v] = si_new; } diff --git a/lib/correlation.h b/lib/correlation.h index 26c3a99..49c6ff2 100644 --- a/lib/correlation.h +++ b/lib/correlation.h @@ -11,7 +11,7 @@ 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->ImF[j]); + total += norm_squared(s->ReF[j]) + norm_squared(s->ReF[j]); } return total / s->D; diff --git a/lib/ising.h b/lib/ising.h index f08ae58..5932ecf 100644 --- a/lib/ising.h +++ b/lib/ising.h @@ -31,58 +31,49 @@ class ising_t { typedef int M_t; typedef double F_t; -}; -void init(ising_t *p) { - p->x = false; -} + ising_t() { + x = false; + } -void free_spin(ising_t s) { - // do nothing! -} + ising_t(bool x) : x(x) {} -void free_spin(int s) { - // do nothing -} + inline int operator*(v_t a) const { + if (x) { + return -(int)a; + } else { + return (int)a; + } + } -void free_spin(double s) { - // do nothing -} + inline double operator*(double a) const { + if (x) { + return -a; + } else { + return a; + } + } -ising_t copy(ising_t s) { - return s; -} +}; -void add(int *s1, int a, ising_t s2) { - if (s2.x) { - *s1 -= a; +inline int& operator+=(int& M, const ising_t &s) { + if (s.x) { + M--; } else { - *s1 += a; + M++; } -} -void add(double *s1, double a, ising_t s2) { - if (s2.x) { - *s1 -= a; - } else { - *s1 += a; - } + return M; } -int scalar_multiple(int factor, ising_t s) { +inline int& operator-=(int& M, const ising_t &s) { if (s.x) { - return -factor; + M++; } else { - return factor; + M--; } -} -double scalar_multiple(double factor, ising_t s) { - if (s.x) { - return -factor; - } else { - return factor; - } + return M; } double norm_squared(double s) { diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 2a391ae..c0958f2 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -53,20 +53,19 @@ void free_spin (orthogonal_t <q, T> m) { template <q_t q, class T> vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) { vector_t <q, T> v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); if (m.is_reflection) { double prod = 0; for (q_t i = 0; i < q; i++) { - prod += v.x[i] * m.x[i]; + prod += v[i] * m.x[i]; } for (q_t i = 0; i < q; i++) { - v_rot.x[i] = v.x[i] - 2 * prod * m.x[i]; + v_rot[i] = v[i] - 2 * prod * m.x[i]; } } else { for (q_t i = 0; i < q; i++) { for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * i + j] * v.x[j]; + v_rot[i] += m.x[q * i + j] * v[j]; } } } @@ -112,11 +111,10 @@ vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) { return act(m, v); // reflections are their own inverse } else { vector_t <q, T> v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); for (q_t i = 0; i < q; i++) { for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * j + i] * v.x[j]; + v_rot[i] += m.x[q * j + i] * v[j]; } } @@ -176,33 +174,32 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q if (n > 1) { unsigned int rotation = gsl_rng_uniform_int(r, n); - v.x = (double *)malloc(q * sizeof(double)); double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); - v.x[0] = v0.x[0] * cosr - v0.x[1] * sinr; - v.x[1] = v0.x[1] * cosr + v0.x[0] * sinr; + v[0] = v0[0] * cosr - v0[1] * sinr; + v[1] = v0[1] * cosr + v0[0] * sinr; for (q_t i = 2; i < q; i++) { - v.x[i] = v0.x[i]; + v[i] = v0[i]; } } else { - v.x = v0.x; + v = v0; } double m_dot_v = 0; for (q_t i = 0; i < q; i++) { m.x[i] = gsl_ran_ugaussian(r); - m_dot_v += m.x[i] * v.x[i]; + m_dot_v += m.x[i] * v[i]; } double v2 = 0; double factor = epsilon * gsl_ran_ugaussian(r); for (q_t i = 0; i < q; i++) { - m.x[i] = m.x[i] - m_dot_v * v.x[i] + factor * v.x[i]; + m.x[i] = m.x[i] - m_dot_v * v[i] + factor * v[i]; v2 += pow(m.x[i], 2); } @@ -212,10 +209,6 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q m.x[i] /= mag_v; } - if (n > 1) { - free(v.x); - } - return m; } diff --git a/lib/state.h b/lib/state.h index 3cef157..360e8f8 100644 --- a/lib/state.h +++ b/lib/state.h @@ -2,6 +2,7 @@ #pragma once #include <functional> +#include <vector> #include "types.h" #include "graph.h" @@ -15,7 +16,7 @@ class state_t { v_t ne; graph_t g; double T; - X_t *spins; + std::vector<X_t> spins; R_t R; double E; typename X_t::M_t M; // the "sum" of the spins, like the total magnetization @@ -29,23 +30,20 @@ class state_t { std::function <double(X_t, X_t)> J; std::function <double(X_t)> H; - state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), T(T), J(J), H(H), g(D, L) { + state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), g(D, L), T(T), J(J), H(H) { nv = g.nv; ne = g.ne; g.add_ext(); - spins = (X_t *)malloc(nv * sizeof(X_t)); - for (v_t i = 0; i < nv; i++) { - init (&(spins[i])); - } + spins.resize(nv); init (&R); E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]); - M = scalar_multiple((int)nv, spins[0]); + M = spins[0] * nv; last_cluster_size = 0; ReF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t)); ImF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t)); for (D_t i = 0; i < D; i++) { - ReF[i] = scalar_multiple(0.0, spins[0]); - ImF[i] = scalar_multiple(0.0, spins[0]); + ReF[i] = spins[0] * 0.0; + ImF[i] = spins[0] * 0.0; } precomputed_cos = (double *)malloc(L * sizeof(double)); precomputed_sin = (double *)malloc(L * sizeof(double)); @@ -56,16 +54,7 @@ class state_t { } ~state_t() { - for (v_t i = 0; i < nv; i++) { - free_spin(spins[i]); - } - 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); free(precomputed_sin); diff --git a/lib/vector.h b/lib/vector.h index 2fe6ab8..d5ced03 100644 --- a/lib/vector.h +++ b/lib/vector.h @@ -3,6 +3,7 @@ #include <stdlib.h> #include <cmath> +#include <array> #include "types.h" @@ -24,78 +25,67 @@ */ template <q_t q, class T> -class vector_t { +class vector_t : public std::array<T, q> { public: - T *x; // M_t needs to hold the sum of nv spins typedef vector_t <q, T> M_t; // F_t needs to hold the double-weighted sum of spins typedef vector_t <q, double> F_t; -}; - -template <q_t q, class T> -void init(vector_t <q, T> *ptr) { - ptr->x = (T *)calloc(q, sizeof(T)); - // initialize a unit vector - ptr->x[0] = (T)1; -} - -template <q_t q, class T> -void free_spin (vector_t <q, T> v) { - free(v.x); -} + vector_t() { + this->fill((T)0); + (*this)[1] = (T)1; + } -template <q_t q, class T> -vector_t <q, T> copy (vector_t <q, T> v) { - vector_t <q, T> v_copy; - - v_copy.x = (T *)calloc(q, sizeof(T)); + vector_t(const T *x) { + for (q_t i = 0; i < q; i++) { + (*this)[i] = x[i]; + } + } - for (q_t i = 0; i < q; i++) { - v_copy.x[i] = v.x[i]; - } + template <class U> + inline vector_t<q, T>& operator+=(const vector_t<q, U> &v) { + for (q_t i = 0; i < q; i++) { + (*this)[i] += (U)v[i]; + } + return *this; + } - return v_copy; -} + template <class U> + inline vector_t<q, T>& operator-=(const vector_t<q, U> &v) { + for (q_t i = 0; i < q; i++) { + (*this)[i] -= (U)v[i]; + } + return *this; + } -template <q_t q, class T, class U, class V> -void add(vector_t<q, U> *v1, V a, vector_t <q, T> v2) { - for (q_t i = 0; i < q; i++) { - v1->x[i] += (U)(a * v2.x[i]); - } -} + inline vector_t<q, T> operator*(v_t x) const { + vector_t<q, T> result; + for (q_t i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } -template <q_t q, class T> -vector_t <q, T> scalar_multiple(int a, vector_t <q, T> v) { - vector_t <q, T> multiple; - multiple.x = (T *)malloc(q * sizeof(T)); - for (q_t i = 0; i < q; i++) { - multiple.x[i] = a * v.x[i]; - } + return result; + } - return multiple; -} + inline vector_t<q, double> operator*(double x) const { + vector_t<q, double> result; + for (q_t i = 0; i < q; i++) { + result[i] = x * (*this)[i]; + } -template <q_t q, class T> -vector_t <q, T> scalar_multiple(double a, vector_t <q, T> v) { - vector_t <q, T> multiple; - multiple.x = (T *)malloc(q * sizeof(T)); - for (q_t i = 0; i < q; i++) { - multiple.x[i] = a * v.x[i]; - } + return result; + } +}; - return multiple; -} -template <q_t q, class T> -double norm_squared (vector_t <q, T> v) { +template<q_t q> +double norm_squared(vector_t<q, double> v) { double tmp = 0; - - for (q_t i = 0; i < q; i++) { - tmp += pow(v.x[i], 2); + for (double &x : v) { + tmp += pow(x, 2); } return tmp; @@ -103,7 +93,9 @@ double norm_squared (vector_t <q, T> v) { template <q_t q, class T> void write_magnetization(vector_t <q, T> M, FILE *outfile) { - fwrite(M.x, sizeof(T), q, outfile); + for (q_t i = 0; i < q; i++) { + fwrite(&(M[i]), sizeof(T), q, outfile); + } } // below functions and definitions are unnecessary for wolff.h but useful. @@ -111,7 +103,7 @@ void write_magnetization(vector_t <q, T> M, FILE *outfile) { template <q_t q> // save some space and don't write whole doubles void write_magnetization(vector_t <q, double> M, FILE *outfile) { for (q_t i = 0; i < q; i++) { - float M_tmp = (float)M.x[i]; + float M_tmp = (float)M[i]; fwrite(&M_tmp, sizeof(float), 1, outfile); } } @@ -121,7 +113,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) { T prod = 0; for (q_t i = 0; i < q; i++) { - prod += v1.x[i] * v2.x[i]; + prod += v1[i] * v2[i]; } return prod; @@ -129,8 +121,7 @@ T dot(vector_t <q, T> v1, vector_t <q, T> v2) { template <q_t q, class T> double H_vector(vector_t <q, T> v1, T *H) { - vector_t <q, T> H_vec; - H_vec.x = H; + vector_t <q, T> H_vec(H); return (double)(dot <q, T> (v1, H_vec)); } |