summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/cluster.h8
-rw-r--r--lib/correlation.h36
-rw-r--r--lib/height.h62
-rw-r--r--lib/orthogonal.h3
-rw-r--r--lib/state.h2
-rw-r--r--lib/vector.h15
-rw-r--r--lib/wolff.h43
7 files changed, 123 insertions, 46 deletions
diff --git a/lib/cluster.h b/lib/cluster.h
index 5f770ad..8061d34 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -52,7 +52,7 @@ template <class T>
T subtract(T, T);
template <class R_t, class X_t>
-v_t 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> *state, v_t v0, R_t r, gsl_rng *rand) {
v_t nv = 0;
ll_t *stack = NULL; // create a new stack
@@ -104,8 +104,8 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
double dE = state->H(rs_old) - state->H(rs_new);
prob = 1.0 - exp(-dE / state->T);
- subtract (state->M, rs_old);
- add (state->M, rs_new);
+ subtract (&(state->M), rs_old);
+ add (&(state->M), rs_new);
state->E += dE;
free_spin (rs_old);
@@ -137,6 +137,6 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
free(marks);
- return nv;
+ state->last_cluster_size = nv;
}
diff --git a/lib/correlation.h b/lib/correlation.h
new file mode 100644
index 0000000..5722aab
--- /dev/null
+++ b/lib/correlation.h
@@ -0,0 +1,36 @@
+
+#pragma once
+
+#include "types.h"
+#include "state.h"
+
+#include <fftw3.h>
+
+template <class R_t, class X_t>
+double correlation_length(const state_t <R_t, X_t> *s) {
+ double *data = (double *)fftw_malloc(s->nv * sizeof(double));
+ int rank = s->D;
+ int *n = (int *)malloc(rank * sizeof(int));
+ fftw_r2r_kind *kind = (fftw_r2r_kind *)malloc(rank * sizeof(fftw_r2r_kind));
+ for (D_t i = 0; i < rank; i++) {
+ n[i] = s->L;
+ kind[i] = FFTW_R2HC;
+ }
+ fftw_plan plan = fftw_plan_r2r(rank, n, data, data, kind, 0);
+
+ for (v_t i = 0; i < s->nv; i++) {
+ data[i] = correlation_component(s->spins[i]);
+ }
+
+ fftw_execute(plan);
+
+ double length = pow(data[0], 2);
+
+ fftw_destroy_plan(plan);
+ fftw_free(data);
+ free(n);
+ free(kind);
+
+ return length;
+}
+
diff --git a/lib/height.h b/lib/height.h
new file mode 100644
index 0000000..16c4dc5
--- /dev/null
+++ b/lib/height.h
@@ -0,0 +1,62 @@
+
+#pragma once
+
+#include <cmath>
+
+#include "types.h"
+
+// object definition
+template <class T>
+struct height_t { T x; };
+
+// init, copy, add, subtract, scalar_multiple, free_spin, and
+// write_magnetization are necessary for the operation of wolff.h
+template <class T>
+void init(height_t *ptr) {
+ ptr->x = (T)0;
+}
+
+template <class T>
+height_t <T> copy (height_t h) {
+ return h;
+}
+
+template <class T>
+void add (height_t <T> *h1, height_t <T> h2) {
+ h1->x += h2.x;
+}
+
+template <class T>
+void subtract (height_t <T> *h1, height_T <T> h2) {
+ h1->x -= h2.x;
+}
+
+template <class T>
+height_t <T> scalar_multiple(v_t a, height_t <T> h) {
+ height_t <T> hm;
+ hm.x = a * h.x;
+
+ return hm;
+}
+
+template <class T>
+void free_spin (height_t <T> h) {
+}
+
+template <class T>
+void write_magnetization(height_t <T> M, FILE *outfile) {
+ fwrite(&(M.x), sizeof(T), 1, outfile);
+}
+
+template <class T>
+double correlation_component(height_t <T> h) {
+ return (double)h.x;
+}
+
+// below here are not necessary for operation
+
+template <class T>
+T dot(height_t <T> h1, height_t <T> h2) {
+ return (h1.x - h2.x) * (h1.x - h2.x);
+}
+
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index d8ad33d..340ee2c 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -184,9 +184,10 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state
}
double v2 = 0;
+ double factor = gsl_ran_ugaussian(r);
for (q_t i = 0; i < q; i++) {
- tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * gsl_ran_ugaussian(r);
+ tmp_v.x[i] = (tmp_v.x[i] - tmpM * s->M.x[i] / M2) + epsilon * factor * s->M.x[i] / sqrt(M2);
v2 += tmp_v.x[i] * tmp_v.x[i];
}
diff --git a/lib/state.h b/lib/state.h
index 8d4265b..1491938 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -19,6 +19,7 @@ class state_t {
R_t R;
double E;
X_t M; // the "sum" of the spins, like the total magnetization
+ v_t last_cluster_size;
std::function <double(X_t, X_t)> J;
std::function <double(X_t)> H;
@@ -36,6 +37,7 @@ class state_t {
init (&R);
E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
M = scalar_multiple (nv, spins[0]);
+ last_cluster_size = 0;
}
~state_t() {
diff --git a/lib/vector.h b/lib/vector.h
index 62ce59e..88934c1 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -30,16 +30,16 @@ vector_t <q, T> copy (vector_t <q, T> v) {
}
template <q_t q, class T>
-void add (vector_t <q, T> v1, vector_t <q, T> v2) {
+void add (vector_t <q, T> *v1, vector_t <q, T> v2) {
for (q_t i = 0; i < q; i++) {
- v1.x[i] += v2.x[i];
+ v1->x[i] += v2.x[i];
}
}
template <q_t q, class T>
-void subtract (vector_t <q, T> v1, vector_t <q, T> v2) {
+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];
+ v1->x[i] -= v2.x[i];
}
}
@@ -72,7 +72,7 @@ void free_spin (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(double), q, outfile);
+ fwrite(M.x, sizeof(T), q, outfile);
}
template <q_t q> // save some space and don't write whole doubles
@@ -83,3 +83,8 @@ void write_magnetization(vector_t <q, double> M, FILE *outfile) {
}
}
+template <q_t q, class T>
+double correlation_component(vector_t <q, T> v) {
+ return (double)v.x[0];
+}
+
diff --git a/lib/wolff.h b/lib/wolff.h
index 90470d5..2d7e8b2 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -1,7 +1,4 @@
-#include <time.h>
-#include <getopt.h>
-
#include "cluster.h"
#include "state.h"
@@ -13,7 +10,7 @@ double H_vector(vector_t <q, T> v1, T *H) {
}
template <class R_t, class X_t>
-void wolff(count_t N, D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, std::function <R_t(gsl_rng *, const state_t <R_t, X_t> *)> gen_R, unsigned long timestamp, bool silent) {
+void wolff(count_t N, D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, std::function <R_t(gsl_rng *, const state_t <R_t, X_t> *)> gen_R, unsigned int n_measurements, std::function <void(const state_t <R_t, X_t> *)> *measurements, bool silent) {
state_t <R_t, X_t> s(D, L, T, J, H);
@@ -21,52 +18,26 @@ void wolff(count_t N, D_t D, L_t L, double T, std::function <double(X_t, X_t)> J
gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
gsl_rng_set(r, rand_seed());
- char *filename_M = (char *)malloc(255 * sizeof(char));
- char *filename_E = (char *)malloc(255 * sizeof(char));
- char *filename_S = (char *)malloc(255 * sizeof(char));
-
- sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
- sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
- sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
-
- FILE *outfile_M = fopen(filename_M, "wb");
- FILE *outfile_E = fopen(filename_E, "wb");
- FILE *outfile_S = fopen(filename_S, "wb");
-
- free(filename_M);
- free(filename_E);
- free(filename_S);
-
- v_t cluster_size = 0;
-
if (!silent) printf("\n");
for (count_t steps = 0; steps < N; steps++) {
- if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, cluster_size);
+ if (!silent) printf("\033[F\033[JWOLFF: sweep %" 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);
-
- cluster_size = flip_cluster <R_t, X_t> (&s, v0, step, r);
-
+ flip_cluster <R_t, X_t> (&s, v0, step, r);
free_spin(step);
- {
- float smaller_E = (float)s.E;
- fwrite(&smaller_E, sizeof(float), 1, outfile_E);
+ for (unsigned int i = 0; i < n_measurements; i++) {
+ measurements[i](&s);
}
- write_magnetization(s.M, outfile_M);
- fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S);
}
if (!silent) {
printf("\033[F\033[J");
}
- printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, cluster_size);
+ printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, s.last_cluster_size);
- fclose(outfile_M);
- fclose(outfile_E);
- fclose(outfile_S);
+ fftw_cleanup();
gsl_rng_free(r);
}