diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/wolff.h | 2 | ||||
| -rw-r--r-- | lib/wolff_tools.c | 50 |
2 files changed, 25 insertions, 27 deletions
diff --git a/lib/wolff.h b/lib/wolff.h index d3bc412..4270f33 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -36,7 +36,5 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, graph_t *graph_add_ext(const graph_t *g); -double *get_bond_probs(double T, double H, ising_state_t *s); - uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps); diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c index f45e1b5..74e5bd6 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -73,8 +73,10 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, x[v] = !x[v]; // flip the vertex for (uint32_t i = 0; i < nn; i++) { - uint32_t e, v1, v2, vn; + bool is_ext; + uint32_t e, v1, v2, vn, ind; int32_t *bond_counter; + double prob; e = g->ve[g->vei[v] + i]; // select the ith bond connected to site v1 = g->ev[2 * e]; @@ -82,14 +84,16 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself - bond_counter = - (v1 == g->nv - 1 || v2 == g->nv - 1) ? &(c->dHb) : &(c->dJb); + is_ext = (v1 == g->nv - 1 || v2 == g->nv - 1); + + bond_counter = is_ext ? &(c->dHb) : &(c->dJb); + prob = is_ext ? ps[1] : ps[0]; if (x[vn] == x0) { // if the neighboring site matches the flipping cluster... (*bond_counter)++; - if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]... + if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]... stack_push(&stack, vn); // push the neighboring vertex to the stack } } else { @@ -106,35 +110,23 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, return c; } -double *get_bond_probs(double T, double H, ising_state_t *s) { - double p = 1 - exp(-2 / T); - double q = 1 - exp(-2 * fabs(H) / T); - - double *ps = (double *)malloc(s->g->ne * sizeof(double)); - - for (uint32_t i = 0; i < s->g->ne; i++) { - uint32_t v1, v2; - v1 = s->g->ev[2 * i]; - v2 = s->g->ev[2 * i + 1]; - if (v1 == s->g->nv - 1 || v2 == s->g->nv - 1) { - ps[i] = q; - } else { - ps[i] = p; - } - } - - return ps; -} - uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps) { + bool no_r, no_ps; + no_r = false; + no_ps = false; + if (r == NULL) { + no_r = true; r = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(r, jst_rand_seed()); } if (ps == NULL) { - ps = get_bond_probs(T, H, s); + no_ps = true; + ps = (double *)malloc(2 * sizeof(double)); + ps[0] = 1 - exp(-2 / T); + ps[1] = 1 - exp(-2 * fabs(H) / T); } cluster_t *c = flip_cluster(s->g, ps, s->spins, r); @@ -146,5 +138,13 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, free(c); + if (no_ps) { + free(ps); + } + + if (no_r) { + gsl_rng_free(r); + } + return n_flips; } |