summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--CMakeLists.txt4
-rw-r--r--lib/cluster.h14
-rw-r--r--lib/ising.h80
-rw-r--r--lib/measure.h11
-rw-r--r--lib/measurement.c145
-rw-r--r--lib/measurement.h25
-rw-r--r--lib/state.h13
-rw-r--r--lib/vector.h98
-rw-r--r--lib/wolff.h21
-rw-r--r--lib/z2.h64
-rw-r--r--src/wolff_On.cpp39
-rw-r--r--src/wolff_ising.cpp160
12 files changed, 442 insertions, 232 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 98ef947..f9ee423 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,6 +14,7 @@ file(GLOB CSOURCES lib/*.c)
file(GLOB CPPSOURCES lib/*.cpp)
add_executable(wolff_finite src/wolff_finite.c ${CSOURCES})
+add_executable(wolff_ising src/wolff_ising.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(wolff_planar src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(wolff_heisenberg src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(analyze_correlations src/analyze_correlations.cpp ${CPPSOURCES} ${CSOURCES})
@@ -28,9 +29,10 @@ if (OPENMP_FOUND)
endif()
target_link_libraries(wolff_finite gsl cblas fftw3 m)
+target_link_libraries(wolff_ising gsl cblas fftw3 m)
target_link_libraries(wolff_heisenberg gsl cblas fftw3 m)
target_link_libraries(wolff_planar gsl cblas fftw3 m)
target_link_libraries(analyze_correlations gsl cblas fftw3 m)
-install(TARGETS wolff_finite wolff_heisenberg wolff_planar analyze_correlations DESTINATION bin)
+install(TARGETS wolff_finite wolff_ising wolff_heisenberg wolff_planar analyze_correlations DESTINATION bin)
diff --git a/lib/cluster.h b/lib/cluster.h
index 2225e2f..4b3c2f6 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -23,6 +23,8 @@
#include "measurement.h"
#include "vector.h"
#include "orthogonal.h"
+#include "ising.h"
+#include "z2.h"
#include "dihedral.h"
#include "dihinf.h"
#include "yule_walker.h"
@@ -107,17 +109,17 @@ void 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);
+ add(&(state->M), -1, rs_old);
+ add(&(state->M), 1, 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);
+ add(&(state->ReF[i]), -cos(2 * M_PI * x), rs_old);
+ 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);
+ add(&(state->ImF[i]), -sin(2 * M_PI * x), rs_old);
+ add(&(state->ImF[i]), sin(2 * M_PI * x), rs_new);
}
free_spin (rs_old);
diff --git a/lib/ising.h b/lib/ising.h
new file mode 100644
index 0000000..4e5164b
--- /dev/null
+++ b/lib/ising.h
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <cmath>
+#include <stdlib.h>
+
+#include "types.h"
+
+class ising_t {
+ public:
+ bool x;
+
+ typedef int M_t;
+ typedef double F_t;
+};
+
+void init(ising_t *p) {
+ p->x = false;
+}
+
+void free_spin(ising_t s) {
+ // do nothing!
+}
+
+void free_spin(int s) {
+ // do nothing
+}
+
+void free_spin(double s) {
+ // do nothing
+}
+
+ising_t copy(ising_t s) {
+ return s;
+}
+
+template <class T>
+void add(T *s1, T a, ising_t s2) {
+ if (s2.x) {
+ *s1 -= a;
+ } else {
+ *s1 += a;
+ }
+}
+
+int scalar_multiple(int factor, ising_t s) {
+ if (s.x) {
+ return -factor;
+ } else {
+ return factor;
+ }
+}
+
+
+double norm_squared(double s) {
+ return pow(s, 2);
+}
+
+template <class T>
+void write_magnetization(T M, FILE *outfile) {
+ fwrite(&M, sizeof(T), 1, outfile);
+}
+
+// below this line is unnecessary, but convenient
+
+double ising_dot(ising_t s1, ising_t s2) {
+ if (s1.x == s2.x) {
+ return 1.0;
+ } else {
+ return -1.0;
+ }
+}
+
+double scalar_field(ising_t s, double H) {
+ if (s.x) {
+ return -H;
+ } else {
+ return H;
+ }
+}
+
diff --git a/lib/measure.h b/lib/measure.h
index d20c081..52e43af 100644
--- a/lib/measure.h
+++ b/lib/measure.h
@@ -1,7 +1,9 @@
#pragma once
-#define POSSIBLE_MEASUREMENTS 4
+#include "measurement.h"
+
+#define POSSIBLE_MEASUREMENTS 5
const unsigned char measurement_energy = 1 << 0;
const unsigned char measurement_clusterSize = 1 << 1;
const unsigned char measurement_magnetization = 1 << 2;
@@ -43,6 +45,13 @@ std::function <void(const state_t <R_t, X_t> *)> measurement_fourier_file(FILE *
};
}
+template <class R_t, class X_t>
+std::function <void(const state_t <R_t, X_t> *)> measurement_average_cluster(meas_t *x) {
+ return [=](const state_t <R_t, X_t> *s) {
+ meas_update(x, s->last_cluster_size);
+ };
+}
+
#endif
diff --git a/lib/measurement.c b/lib/measurement.c
index b30cf6b..ad824f6 100644
--- a/lib/measurement.c
+++ b/lib/measurement.c
@@ -1,15 +1,6 @@
-#include "convex.h"
#include "measurement.h"
-meas_t *meas_initialize(count_t W) {
- meas_t *m = (meas_t *)calloc(1, sizeof(meas_t));
- m->W = W;
- m->xx = (double *)calloc(2 * W + 1, sizeof(double));
-
- return m;
-}
-
double add_to_avg(double mx, double x, count_t n) {
return mx * (n / (n + 1.0)) + x / (n + 1.0);
}
@@ -19,42 +10,24 @@ void meas_update(meas_t *m, double x) {
m->x = add_to_avg(m->x, x, n);
m->x2 = add_to_avg(m->x2, pow(x, 2), n);
- m->x4 = add_to_avg(m->x4, pow(x, 4), n);
m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n);
m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n);
- dll_t *tmp_window = m->x_window;
- dll_t *pos_save;
- count_t t = 0;
-
- while (tmp_window != NULL) {
- m->xx[t] = add_to_avg(m->xx[t], x * (tmp_window->x), m->n - t - 1);
- t++;
- if (t == 2 * (m->W)) {
- pos_save = tmp_window;
- }
- tmp_window = tmp_window->next;
+ /*
+ if (n > 1) {
+ double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2));
+ m->dx = sqrt(s2 / n);
+ m->c = s2;
+ m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n);
}
-
- if (t == 2 * (m->W) + 1) {
- if (2 * (m->W) + 1 == 1) {
- free(m->x_window);
- m->x_window = NULL;
- } else {
- free(pos_save->next);
- pos_save->next = NULL;
- }
- }
-
- stack_push_d(&(m->x_window), x);
+ */
(m->n)++;
}
double meas_dx(const meas_t *m) {
- return 2 * get_tau(m) * Cxx(m, 0) / m->n;
-// return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
+ return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
}
double meas_c(const meas_t *m) {
@@ -101,105 +74,3 @@ double rho(const autocorr_t *o, count_t i) {
return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
}
-double Cxx(const meas_t *m, count_t t) {
- return m->xx[t] - pow(m->x, 2);
-}
-
-double rho_m(const meas_t *m, count_t t) {
- return Cxx(m, t) / Cxx(m, 0);
-}
-
-double get_tau(const meas_t *m) {
- double *Gammas = (double *)malloc((m->W + 1) * sizeof(double));
-
- Gammas[0] = 1 + rho_m(m, 0);
- for (uint64_t i = 0; i < m->W; i++) {
- Gammas[1 + i] = rho_m(m, 2 * i + 1) + rho_m(m, 2 * i + 2);
- }
-
- uint64_t n;
- for (n = 0; n < m->W + 1; n++) {
- if (Gammas[n] <= 0) {
- break;
- }
- }
-
- double *conv_Gamma = get_convex_minorant(n, Gammas);
-
- double tau = - 0.5;
-
- for (uint64_t i = 0; i < n + 1; i++) {
- tau += conv_Gamma[i];
- }
-
- free(Gammas);
-
- return tau;
-}
-
-void print_meas(const meas_t *m, const char *sym, FILE *outfile) {
- fprintf(outfile, "%s-><|n->%" PRIcount ",x->%.15f,x^2->%.15f,x^4->%.15f,xx->{", sym, m->n, m->x, m->x2, m->x4);
- for (count_t i = 0; i < 2 * (m->W) + 1; i++) {
- fprintf(outfile, "%.15f", m->xx[i]);
- if (i < 2 * (m->W)) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "}|>");
-}
-
-void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile) {
- fprintf(outfile, "%s-><|n->{", sym);
- for (q_t i = 0; i < q; i++) {
- fprintf(outfile, "%" PRIcount, m[i]->n);
- if (i < q - 1) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "},x->{");
- for (q_t i = 0; i < q; i++) {
- fprintf(outfile, "%.15f", m[i]->x);
- if (i < q - 1) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "},x^2->{");
- for (q_t i = 0; i < q; i++) {
- fprintf(outfile, "%.15f", m[i]->x2);
- if (i < q - 1) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "},x^4->{");
- for (q_t i = 0; i < q; i++) {
- fprintf(outfile, "%.15f", m[i]->x4);
- if (i < q - 1) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "},xx->{");
- for (q_t i = 0; i < q; i++) {
- fprintf(outfile, "{");
- for (count_t j = 0; j < 2 * (m[i]->W) + 1; j++) {
- fprintf(outfile, "%.15f", m[i]->xx[j]);
- if (j < 2 * (m[i]->W)) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "}");
- if (i < q - 1) {
- fprintf(outfile, ",");
- }
- }
- fprintf(outfile, "}|>");
-}
-
-void free_meas(meas_t *m) {
- free(m->xx);
- while (m->x_window != NULL) {
- stack_pop_d(&(m->x_window));
- }
- free(m);
-}
-
-
diff --git a/lib/measurement.h b/lib/measurement.h
index d9bd52e..78fa51b 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -3,21 +3,20 @@
#include <math.h>
#include <stdlib.h>
-#include <stdio.h>
#include "types.h"
#include "stack.h"
+#ifdef __cplusplus
+extern "C" {
+#endif
+
typedef struct {
- count_t n;
+ uint64_t n;
double x;
double x2;
- double x4;
double m2;
double m4;
- count_t W;
- double *xx;
- dll_t *x_window;
} meas_t;
typedef struct {
@@ -41,14 +40,6 @@ void update_autocorr(autocorr_t *OO, double O);
double rho(const autocorr_t *o, uint64_t i);
-void print_meas(const meas_t *m, const char *sym, FILE *outfile);
-void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile);
-
-void free_meas(meas_t *m);
-
-meas_t *meas_initialize(count_t W);
-
-double get_tau(const meas_t *m);
-
-double Cxx(const meas_t *m, count_t t);
-
+#ifdef __cplusplus
+}
+#endif
diff --git a/lib/state.h b/lib/state.h
index c70459f..b40b85c 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -5,6 +5,7 @@
#include "types.h"
#include "graph.h"
+#include "ising.h"
template <class R_t, class X_t>
class state_t {
@@ -18,10 +19,10 @@ class state_t {
X_t *spins;
R_t R;
double E;
- X_t M; // the "sum" of the spins, like the total magnetization
+ typename X_t::M_t M; // the "sum" of the spins, like the total magnetization
v_t last_cluster_size;
- X_t *ReF;
- X_t *ImF;
+ typename X_t::F_t *ReF;
+ typename X_t::F_t *ImF;
std::function <double(X_t, X_t)> J;
std::function <double(X_t)> H;
@@ -38,10 +39,10 @@ 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((int)nv, spins[0]);
last_cluster_size = 0;
- ReF = (X_t *)malloc(D * sizeof(X_t));
- ImF = (X_t *)malloc(D * sizeof(X_t));
+ 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, spins[0]);
ImF[i] = scalar_multiple(0, spins[0]);
diff --git a/lib/vector.h b/lib/vector.h
index a1084e2..c478618 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -6,17 +6,49 @@
#include "types.h"
+/* The following is the minimum definition of a spin class.
+ *
+ * The class must contain an M_t and an F_t for holding the sum of an
+ * integer number of spins and a double-weighted number of spins,
+ * respectively.
+ *
+ * void init(X_t *p);
+ * void free_spin(X_t p);
+ * X_t copy(X_t x);
+ * void add(M_t *x1, int factor, X_t x2);
+ * void add(F_t *x1, double factor, X_t x2);
+ * M_t scalar_multiple(int factor, X_t x);
+ * double norm_squared(F_t x);
+ * void write_magnetization(M_t M, FILE *outfile);
+ *
+ */
+
template <q_t q, class T>
-struct vector_t { T *x; };
+class vector_t {
+ 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;
@@ -29,46 +61,15 @@ vector_t <q, T> copy (vector_t <q, T> v) {
return v_copy;
}
-template <q_t q, class T>
-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];
- }
-}
-
-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];
- }
-}
-
-template <q_t q, class T>
-T norm_squared (vector_t <q, T> v) {
- T tmp = 0;
+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++) {
- tmp += v.x[i] * v.x[i];
+ v1->x[i] += (U)(a * v2.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> 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++) {
@@ -79,19 +80,14 @@ vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
}
template <q_t q, class T>
-T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
- T prod = 0;
+double norm_squared (vector_t <q, T> v) {
+ double tmp = 0;
for (q_t i = 0; i < q; i++) {
- prod += v1.x[i] * v2.x[i];
+ tmp += pow(v.x[i], 2);
}
- return prod;
-}
-
-template <q_t q, class T>
-void free_spin (vector_t <q, T> v) {
- free(v.x);
+ return tmp;
}
template <q_t q, class T>
@@ -99,6 +95,8 @@ void write_magnetization(vector_t <q, T> M, FILE *outfile) {
fwrite(M.x, sizeof(T), q, outfile);
}
+// below functions and definitions are unnecessary for wolff.h but useful.
+
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++) {
@@ -108,8 +106,14 @@ 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];
+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];
+ }
+
+ return prod;
}
template <q_t q, class T>
diff --git a/lib/wolff.h b/lib/wolff.h
index eb357c7..21e88d9 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -3,33 +3,26 @@
#include "state.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 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);
-
- // initialize random number generator
- gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
- gsl_rng_set(r, rand_seed());
+void wolff(count_t N, state_t <R_t, X_t> *s, 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, gsl_rng *r, bool silent) {
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, 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);
- flip_cluster <R_t, X_t> (&s, v0, step, r);
+ v_t v0 = gsl_rng_uniform_int(r, s->nv);
+ R_t step = gen_R(r, s);
+ flip_cluster <R_t, X_t> (s, v0, step, r);
free_spin(step);
for (unsigned int i = 0; i < n_measurements; i++) {
- measurements[i](&s);
+ measurements[i](s);
}
}
if (!silent) {
printf("\033[F\033[J");
}
- printf("WOLFF: sweep %" 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);
- gsl_rng_free(r);
}
diff --git a/lib/z2.h b/lib/z2.h
new file mode 100644
index 0000000..599a6a5
--- /dev/null
+++ b/lib/z2.h
@@ -0,0 +1,64 @@
+
+#pragma once
+
+#include <gsl/gsl_rng.h>
+#include "types.h"
+#include "ising.h"
+#include "state.h"
+
+struct z2_t { bool x; };
+
+void init(z2_t *p) {
+ p->x = false;
+}
+
+void free_spin(z2_t p) {
+ // do nothing!
+}
+
+z2_t copy(z2_t x) {
+ return x;
+}
+
+ising_t act(z2_t r, ising_t s) {
+ ising_t rs;
+
+ if (r.x) {
+ rs.x = !s.x;
+ return rs;
+ } else {
+ rs.x = s.x;
+ return rs;
+ }
+}
+
+z2_t act(z2_t r1, z2_t r2) {
+ z2_t r3;
+
+ if (r1.x) {
+ r3.x = !r2.x;
+ return r3;
+ } else {
+ r3.x = r2.x;
+ return r3;
+ }
+}
+
+ising_t act_inverse(z2_t r, ising_t s) {
+ return act(r, s);
+}
+
+z2_t act_inverse(z2_t r1, z2_t r2) {
+ return act(r1, r2);
+}
+
+// these are all functions necessary for wolff.h
+
+z2_t generate_ising_rotation(gsl_rng *r, const state_t <z2_t, ising_t> *s) {
+ z2_t rot;
+ rot.x = true;
+ return rot;
+}
+
+
+
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index d718440..76831ba 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -40,6 +40,7 @@ int main(int argc, char *argv[]) {
bool silent = false;
bool use_pert = false;
+ bool N_is_sweeps = false;
bool modulated_field = false;
int order = 2;
@@ -51,7 +52,7 @@ int main(int argc, char *argv[]) {
unsigned char measurement_flags = measurement_energy | measurement_clusterSize;
- while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:")) != -1) {
+ while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:S")) != -1) {
switch (opt) {
case 'N': // number of steps
N = (count_t)atof(optarg);
@@ -82,11 +83,14 @@ int main(int argc, char *argv[]) {
modulated_field = true;
break;
case 'M':
- measurement_flags |= 1 << atoi(optarg);
+ measurement_flags ^= 1 << atoi(optarg);
break;
case 'o':
order = atoi(optarg);
break;
+ case 'S':
+ N_is_sweeps = true;
+ break;
default:
exit(EXIT_FAILURE);
}
@@ -179,6 +183,13 @@ int main(int argc, char *argv[]) {
n_measurements++;
}
+ meas_t *meas_sweeps;
+ if (N_is_sweeps) {
+ meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t));
+ measurements[n_measurements] = measurement_average_cluster<orthogonal_R_t, vector_R_t> (meas_sweeps);
+ n_measurements++;
+ }
+
std::function <double(vector_R_t)> H;
if (modulated_field) {
@@ -187,7 +198,24 @@ int main(int argc, char *argv[]) {
H = std::bind(H_vector <N_COMP, double>, std::placeholders::_1, H_vec);
}
- wolff <orthogonal_R_t, vector_R_t> (N, D, L, T, dot <N_COMP, double>, H, gen_R, n_measurements, measurements, silent);
+ // initialize random number generator
+ gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+ gsl_rng_set(r, rand_seed());
+
+ state_t <orthogonal_R_t, vector_R_t> s(D, L, T, dot <N_COMP, double>, H);
+
+ if (N_is_sweeps) {
+ count_t N_rounds = 0;
+ printf("\n");
+ while (N_rounds * N * meas_sweeps->x < N * s.nv) {
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+ wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, n_measurements, measurements, r, silent);
+ N_rounds++;
+ }
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+ } else {
+ wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, n_measurements, measurements, r, silent);
+ }
free(measurements);
@@ -204,7 +232,12 @@ int main(int argc, char *argv[]) {
fclose(outfile_F);
}
+ if (N_is_sweeps) {
+ free(meas_sweeps);
+ }
+
free(H_vec);
+ gsl_rng_free(r);
return 0;
}
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
new file mode 100644
index 0000000..4d4c791
--- /dev/null
+++ b/src/wolff_ising.cpp
@@ -0,0 +1,160 @@
+
+#include <getopt.h>
+
+#include <wolff.h>
+#include <correlation.h>
+#include <measure.h>
+
+int main(int argc, char *argv[]) {
+
+ count_t N = (count_t)1e7;
+
+ D_t D = 2;
+ L_t L = 128;
+ double T = 2.26918531421;
+ double H = 0.0;
+
+ bool silent = false;
+ bool N_is_sweeps = false;
+
+ int opt;
+
+ unsigned char measurement_flags = measurement_energy | measurement_clusterSize;
+
+ while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:sM:S")) != -1) {
+ switch (opt) {
+ case 'N': // number of steps
+ N = (count_t)atof(optarg);
+ break;
+ case 'D': // dimension
+ D = atoi(optarg);
+ break;
+ case 'L': // linear size
+ L = atoi(optarg);
+ break;
+ case 'T': // temperature
+ T = atof(optarg);
+ break;
+ case 'H': // external field. nth call couples to state n
+ H = atof(optarg);
+ break;
+ case 's': // don't print anything during simulation. speeds up slightly
+ silent = true;
+ break;
+ case 'M':
+ measurement_flags ^= 1 << atoi(optarg);
+ break;
+ case 'S':
+ N_is_sweeps = true;
+ break;
+ default:
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ unsigned long timestamp;
+
+ {
+ struct timespec spec;
+ clock_gettime(CLOCK_REALTIME, &spec);
+ timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec;
+ }
+
+ FILE *outfile_info = fopen("wolff_metadata.txt", "a");
+
+ fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"ISING\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> %.15f |>\n", timestamp, D, L, (v_t)pow(L, D), D * (v_t)pow(L, D), T, H);
+
+ fclose(outfile_info);
+
+ unsigned int n_measurements = 0;
+ std::function <void(const state_t <z2_t, ising_t> *)> *measurements = (std::function <void(const state_t <z2_t, ising_t> *)> *)calloc(POSSIBLE_MEASUREMENTS, sizeof(std::function <void(const state_t <z2_t, ising_t> *)>));
+ FILE *outfile_M, *outfile_E, *outfile_S, *outfile_F;
+
+ if (measurement_flags & measurement_energy) {
+ char *filename_E = (char *)malloc(255 * sizeof(char));
+ sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
+ outfile_E = fopen(filename_E, "wb");
+ free(filename_E);
+ measurements[n_measurements] = measurement_energy_file<z2_t, ising_t> (outfile_E);
+ n_measurements++;
+ }
+
+ if (measurement_flags & measurement_clusterSize) {
+ char *filename_S = (char *)malloc(255 * sizeof(char));
+ sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
+ outfile_S = fopen(filename_S, "wb");
+ free(filename_S);
+ measurements[n_measurements] = measurement_cluster_file<z2_t, ising_t> (outfile_S);
+ n_measurements++;
+ }
+
+ if (measurement_flags & measurement_magnetization) {
+ char *filename_M = (char *)malloc(255 * sizeof(char));
+ sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
+ outfile_M = fopen(filename_M, "wb");
+ free(filename_M);
+ measurements[n_measurements] = measurement_magnetization_file<z2_t, ising_t> (outfile_M);
+ n_measurements++;
+ }
+
+ if (measurement_flags & measurement_fourierZero) {
+ char *filename_F = (char *)malloc(255 * sizeof(char));
+ sprintf(filename_F, "wolff_%lu_F.dat", timestamp);
+ outfile_F = fopen(filename_F, "wb");
+ free(filename_F);
+ measurements[n_measurements] = measurement_fourier_file<z2_t, ising_t> (outfile_F);
+ n_measurements++;
+ }
+
+ meas_t *meas_sweeps;
+ if (N_is_sweeps) {
+ meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t));
+ measurements[n_measurements] = measurement_average_cluster<z2_t, ising_t> (meas_sweeps);
+ n_measurements++;
+ }
+
+ // initialize random number generator
+ gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+ gsl_rng_set(r, rand_seed());
+
+ state_t <z2_t, ising_t> s(D, L, T, ising_dot, std::bind(scalar_field, std::placeholders::_1, H));
+
+ std::function <z2_t(gsl_rng *, const state_t <z2_t, ising_t> *)> gen = generate_ising_rotation;
+
+ if (N_is_sweeps) {
+ count_t N_rounds = 0;
+ printf("\n");
+ while (N_rounds * N * meas_sweeps->x < N * s.nv) {
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+ wolff(N, &s, gen, n_measurements, measurements, r, silent);
+ N_rounds++;
+ }
+ printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+ } else {
+ wolff(N, &s, gen, n_measurements, measurements, r, silent);
+ }
+
+ free(measurements);
+
+ if (measurement_flags & measurement_energy) {
+ fclose(outfile_E);
+ }
+ if (measurement_flags & measurement_clusterSize) {
+ fclose(outfile_S);
+ }
+ if (measurement_flags & measurement_magnetization) {
+ fclose(outfile_M);
+ }
+ if (measurement_flags & measurement_fourierZero) {
+ fclose(outfile_F);
+ }
+
+ if (N_is_sweeps) {
+ free(meas_sweeps);
+ }
+
+ gsl_rng_free(r);
+
+ return 0;
+}
+