summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-26 00:32:38 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-26 00:32:38 -0400
commitd3b3e39a525d0c3aa9663f824ba96e0c429a8313 (patch)
tree942a902af332a532ab76fccbdb51c83e98338d5e /lib
parent225905d1ad11a0f10d63ab0c6c71e6a5265e0894 (diff)
downloadc++-d3b3e39a525d0c3aa9663f824ba96e0c429a8313.tar.gz
c++-d3b3e39a525d0c3aa9663f824ba96e0c429a8313.tar.bz2
c++-d3b3e39a525d0c3aa9663f824ba96e0c429a8313.zip
partially class-ified, ising and On work but potts and height do not
Diffstat (limited to 'lib')
-rw-r--r--lib/cluster.h17
-rw-r--r--lib/correlation.h2
-rw-r--r--lib/ising.h65
-rw-r--r--lib/orthogonal.h27
-rw-r--r--lib/state.h25
-rw-r--r--lib/vector.h123
6 files changed, 110 insertions, 149 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);
-}
-
-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));
- for (q_t i = 0; i < q; i++) {
- v_copy.x[i] = v.x[i];
- }
-
- return v_copy;
-}
-
-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]);
- }
-}
-
-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 multiple;
-}
-
-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];
- }
+ vector_t() {
+ this->fill((T)0);
+ (*this)[1] = (T)1;
+ }
+
+ vector_t(const T *x) {
+ for (q_t i = 0; i < q; i++) {
+ (*this)[i] = 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;
+ }
+
+ 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;
+ }
+
+ 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];
+ }
+
+ return result;
+ }
+
+ 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];
+ }
+
+ 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));
}