summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/wolff.h16
-rw-r--r--lib/wolff_tools.c23
-rw-r--r--src/wolff.c127
3 files changed, 92 insertions, 74 deletions
diff --git a/lib/wolff.h b/lib/wolff.h
index 428ce85..82ccb9e 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -31,6 +31,17 @@ typedef struct {
int32_t dHb;
} cluster_t;
+typedef struct {
+ uint64_t n;
+ double x;
+ double dx;
+ double x2;
+ double m2;
+ double m4;
+ double c;
+ double dc;
+} meas_t;
+
int32_t sign(double x);
cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
@@ -40,3 +51,8 @@ graph_t *graph_add_ext(const graph_t *g);
uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
double *ps);
+
+void update_meas(meas_t *m, double x);
+
+double add_to_avg(double mx, double x, uint64_t n);
+
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 7886bab..59d04fc 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -152,3 +152,26 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
return n_flips;
}
+
+double add_to_avg(double mx, double x, uint64_t n) {
+ return mx * (n / (n + 1.)) + x * 1. / (n + 1.);
+}
+
+void update_meas(meas_t *m, double x) {
+ uint64_t n = m->n;
+
+ m->x = add_to_avg(m->x, x, n);
+ m->x2 = add_to_avg(m->x2, pow(x, 2), 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);
+
+ 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);
+ }
+
+ (m->n)++;
+}
diff --git a/src/wolff.c b/src/wolff.c
index 3b07a13..4c8d744 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -97,34 +97,23 @@ int main(int argc, char *argv[]) {
double diff = 1e31;
uint64_t n_runs = 0;
-
- double E1, E2, dE1, aM1, daM1, M1, M2, dM1, C, dC, aX, daX, X, dX, aMmu2, Mmu2, aMmu4, Mmu4, Emu2, Emu4;
double clust_per_sweep = 0;
- E1 = 0;
- E2 = 0;
- M1 = 0;
- aM1 = 0;
- M2 = 0;
- C = 0;
- dC = 0;
- X = 0;
- aX = 0;
- dX = 0;
- dE1 = 0;
- dM1 = 0;
- Mmu2 = 0;
- Mmu4 = 0;
- aMmu2 = 0;
- aMmu4 = 0;
- Emu2 = 0;
- Emu4 = 0;
+ meas_t *M, *aM, *eM, *mM, *E, *eE, *mE;
+
+ M = calloc(1, sizeof(meas_t));
+ aM = calloc(1, sizeof(meas_t));
+ eM = calloc(1, sizeof(meas_t));
+ mM = calloc(1, sizeof(meas_t));
+ E = calloc(1, sizeof(meas_t));
+ eE = calloc(1, sizeof(meas_t));
+ mE = calloc(1, sizeof(meas_t));
printf("\n");
while (diff > eps || diff == 0. || n_runs < min_runs) {
printf("\033[F\033[JWOLFF: sweep %" PRIu64
", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
- n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+ n_runs, fabs(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
uint32_t n_flips = 0;
uint32_t n_clust = 0;
@@ -134,69 +123,52 @@ int main(int argc, char *argv[]) {
n_clust++;
}
- E1 = E1 * (n_runs / (n_runs + 1.)) + s->H * 1. / (n_runs + 1.);
- M1 = M1 * (n_runs / (n_runs + 1.)) + s->M * 1. / (n_runs + 1.);
- aM1 = aM1 * (n_runs / (n_runs + 1.)) + abs(s->M) * 1. / (n_runs + 1.);
- E2 = E2 * (n_runs / (n_runs + 1.)) + pow(s->H, 2) * 1. / (n_runs + 1.);
- M2 = M2 * (n_runs / (n_runs + 1.)) + pow(s->M, 2) * 1. / (n_runs + 1.);
-
- Mmu2 = Mmu2 * (n_runs / (n_runs + 1.)) +
- pow(s->M - M1, 2) * 1. / (n_runs + 1.);
- Mmu4 = Mmu4 * (n_runs / (n_runs + 1.)) +
- pow(s->M - M1, 4) * 1. / (n_runs + 1.);
- aMmu2 = aMmu2 * (n_runs / (n_runs + 1.)) +
- pow(abs(s->M) - aM1, 2) * 1. / (n_runs + 1.);
- aMmu4 = aMmu4 * (n_runs / (n_runs + 1.)) +
- pow(abs(s->M) - aM1, 4) * 1. / (n_runs + 1.);
- Emu2 = Emu2 * (n_runs / (n_runs + 1.)) +
- pow(s->H - E1, 2) * 1. / (n_runs + 1.);
- Emu4 = Emu4 * (n_runs / (n_runs + 1.)) +
- pow(s->H - E1, 4) * 1. / (n_runs + 1.);
-
- if (n_runs > 1) {
- double Msigma2 = n_runs / (n_runs - 1) * (M2 - pow(M1, 2));
- double aMsigma2 = n_runs / (n_runs - 1) * (M2 - pow(aM1, 2));
- X = Msigma2 / T;
- aX = aMsigma2 / T;
- dX =
- sqrt((Mmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Mmu2, 2)) / n_runs) /
- T;
- daX =
- sqrt((aMmu4 - (n_runs - 3.) / (n_runs - 1.) * pow(aMmu2, 2)) / n_runs) /
- T;
-
-
- double Esigma2 = n_runs / (n_runs - 1) * (E2 - pow(E1, 2));
- C = Esigma2 / T;
- dC =
- sqrt((Emu4 - (n_runs - 3.) / (n_runs - 1.) * pow(Emu2, 2)) / n_runs) /
- T;
-
- dE1 = sqrt(Esigma2 / n_runs);
- dM1 = sqrt(Msigma2 / n_runs);
- daM1 = sqrt(aMsigma2 / n_runs);
-
- if (M_stop) {
- diff = fabs(dM1 / M1);
- } else {
- diff = fabs(dX / X);
- }
+ double ss = 1;
+ if (s->spins[s->g->nv - 1]) {
+ ss = -1;
+ }
+
+ double HH = 1;
+ if (H < 0) {
+ HH = -1;
+ }
+
+ update_meas(M, (double)(s->M));
+ update_meas(aM, HH * fabs((double)(s->M)));
+ update_meas(E, s->H);
+
+ if (HH * s->M * ss > 0) {
+ update_meas(eM, HH * fabs((double)(s->M)));
+ update_meas(eE, s->H);
+ } else {
+ update_meas(mM, - HH * fabs((double)(s->M)));
+ update_meas(mE, s->H);
+ }
+
+ if (M_stop) {
+ diff = fabs(M->dx / M->x);
+ } else {
+ diff = fabs(M->dc / M->c);
}
- clust_per_sweep = clust_per_sweep * (n_runs / (n_runs + 1.)) +
- (n_clust * 1. / N) * 1. / (n_runs + 1.);
+ clust_per_sweep = add_to_avg(clust_per_sweep, n_clust * 1. / N, n_runs);
n_runs++;
}
+
printf("\033[F\033[JWOLFF: sweep %" PRIu64
", dH/H = %.4f, dM/M = %.4f, dC/C = %.4f, dX/X = %.4f, cps: %.1f\n",
- n_runs, fabs(dE1 / E1), dM1 / M1, dC / C, dX / X, clust_per_sweep);
+ n_runs, fabs(E->dx / E->x), M->dx / M->x, E->dc / E->c, M->dc / M->c, clust_per_sweep);
FILE *outfile = fopen("out.dat", "a");
+
fprintf(outfile,
- "%u %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L,
- T, H, E1 / h->nv, dE1 / h->nv, M1 / h->nv, dM1 / h->nv, C / h->nv,
- dC / h->nv, X / h->nv, dX / h->nv, aM1 / h->nv, daM1 / h->nv, aX / h->nv, daX / h->nv);
+ "%u %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %" PRIu64 " %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", L,
+ T, H, n_runs,
+ E->x / h->nv, E->dx / h->nv, M->x / h->nv, M->dx / h->nv, E->c / h->nv, E->dc / h->nv, M->c / h->nv, M->dc / h->nv,
+ eE->n, eE->x / h->nv, eE->dx / h->nv, eM->x / h->nv, eM->dx / h->nv, eE->c / h->nv, eE->dc / h->nv, eM->c / h->nv, eM->dc / h->nv,
+ mE->n, mE->x / h->nv, mE->dx / h->nv, mM->x / h->nv, mM->dx / h->nv, mE->c / h->nv, mE->dc / h->nv, mM->c / h->nv, mM->dc / h->nv,
+ aM->x / h->nv, aM->dx / h->nv, aM->c / h->nv, aM->dc / h->nv);
fclose(outfile);
if (output_state) {
@@ -212,6 +184,13 @@ int main(int argc, char *argv[]) {
graph_free(s->g);
free(s->spins);
free(s);
+ free(M);
+ free(aM);
+ free(eM);
+ free(mM);
+ free(E);
+ free(eE);
+ free(mE);
free(bond_probs);
graph_free(h);