From 9db6ee734df8477a2529f56e4a6f4b1784bf941b Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 22 Jun 2017 16:05:05 -0400 Subject: many changes, simplification of some functions, removal of unneeded ones --- lib/queue.c | 21 +++++++++---------- lib/queue.h | 18 ++++++++++++++++ lib/wolff.h | 44 ++++++++++++++++----------------------- lib/wolff_tools.c | 61 ++++++++++++++++--------------------------------------- 4 files changed, 63 insertions(+), 81 deletions(-) create mode 100644 lib/queue.h (limited to 'lib') diff --git a/lib/queue.c b/lib/queue.c index 9a01741..aebbb2a 100644 --- a/lib/queue.c +++ b/lib/queue.c @@ -1,23 +1,23 @@ -#include "wolff.h" +#include "queue.h" void stack_push(ll_t **q, uint32_t x) { - ll_t *nq = malloc(sizeof(ll_t)); - nq->x = x; - nq->next = *q; + ll_t *nq = malloc(sizeof(ll_t)); + nq->x = x; + nq->next = *q; - *q = nq; + *q = nq; } uint32_t stack_pop(ll_t **q) { - ll_t *old_q = *q; + ll_t *old_q = *q; - *q = old_q->next; - uint32_t x = old_q->x; + *q = old_q->next; + uint32_t x = old_q->x; - free(old_q); + free(old_q); - return x; + return x; } bool stack_contains(const ll_t *q, uint32_t x) { @@ -29,4 +29,3 @@ bool stack_contains(const ll_t *q, uint32_t x) { return stack_contains(q->next, x); } } - diff --git a/lib/queue.h b/lib/queue.h new file mode 100644 index 0000000..c7acec1 --- /dev/null +++ b/lib/queue.h @@ -0,0 +1,18 @@ + +#pragma once + +#include +#include +#include +#include + +typedef struct ll_tag { + uint32_t x; + struct ll_tag *next; +} ll_t; + +void stack_push(ll_t **q, uint32_t x); + +uint32_t stack_pop(ll_t **q); + +bool stack_contains(const ll_t *q, uint32_t x); diff --git a/lib/wolff.h b/lib/wolff.h index cec9ee3..d3bc412 100644 --- a/lib/wolff.h +++ b/lib/wolff.h @@ -1,19 +1,23 @@ -#include -#include -#include +#pragma once + +#include +#include #include -#include -#include +#include #include #include +#include +#include #include -#include -#include +#include +#include #include #include +#include "queue.h" + typedef struct { graph_t *g; bool *spins; @@ -21,32 +25,18 @@ typedef struct { double H; } ising_state_t; -typedef struct ll_tag { - uint32_t x; - struct ll_tag *next; -} ll_t; - typedef struct { uint32_t nv; - double dH; - int32_t dM; + int32_t dJb; + int32_t dHb; } cluster_t; -double get_hamiltonian(graph_t *g, double *coupling, bool *x); - -void stack_push(ll_t **q, uint32_t x); - -uint32_t stack_pop(ll_t **q); - -bool stack_contains(const ll_t *q, uint32_t x); - -cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r); +cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, + gsl_rng *r); graph_t *graph_add_ext(const graph_t *g); -double hh(double th); - 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); - +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 016e19b..f45e1b5 100644 --- a/lib/wolff_tools.c +++ b/lib/wolff_tools.c @@ -1,4 +1,5 @@ +#include "queue.h" #include "wolff.h" graph_t *graph_add_ext(const graph_t *g) { @@ -9,7 +10,7 @@ graph_t *graph_add_ext(const graph_t *g) { h->ev = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t)); h->vei = (uint32_t *)malloc((h->nv + 1) * sizeof(uint32_t)); - h->ve = (uint32_t *) malloc(2 * h->ne * sizeof(uint32_t)); + h->ve = (uint32_t *)malloc(2 * h->ne * sizeof(uint32_t)); h->vx = (double *)malloc(2 * h->nv * sizeof(double)); h->bq = (bool *)malloc(h->nv * sizeof(bool)); @@ -45,42 +46,21 @@ graph_t *graph_add_ext(const graph_t *g) { return h; } -double get_hamiltonian(graph_t *g, double *coupling, bool *x) { - double hamiltonian = 0; - - for (uint32_t i = 0; i < g->ne; i++) { - uint32_t v1, v2; - - v1 = g->ev[2 * i]; - v2 = g->ev[2 * i + 1]; - - if (x[v1] == x[v2]) { - hamiltonian -= coupling[i]; - } else { - hamiltonian += coupling[i]; - } - } - - return hamiltonian; -} - -cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, gsl_rng *r) { +cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, + gsl_rng *r) { uint32_t v0; int32_t n_h_bonds, n_bonds; bool x0; cluster_t *c; - + v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex - x0 = x[v0]; // record its orientation + x0 = x[v0]; // record its orientation - ll_t *stack = NULL; // create a new stack + ll_t *stack = NULL; // create a new stack stack_push(&stack, v0); // push the initial vertex to the stack // initiate the data structure for returning flip information c = (cluster_t *)calloc(1, sizeof(cluster_t)); - c->nv = 0; - - n_h_bonds = 0; n_bonds = 0; while (stack != NULL) { uint32_t v; @@ -90,7 +70,7 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g nn = g->vei[v + 1] - g->vei[v]; if (x[v] == x0) { // if the vertex hasn't already been flipped - x[v] = !x[v]; // flip the vertex + x[v] = !x[v]; // flip the vertex for (uint32_t i = 0; i < nn; i++) { uint32_t e, v1, v2, vn; @@ -102,9 +82,11 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself - bond_counter = (v1 == g->nv - 1 || v2 == g->nv - 1) ? &n_h_bonds : &n_bonds; + bond_counter = + (v1 == g->nv - 1 || v2 == g->nv - 1) ? &(c->dHb) : &(c->dJb); - if (x[vn] == x0) { // if the neighboring site matches the flipping cluster... + 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]... @@ -115,22 +97,15 @@ cluster_t *flip_cluster(const graph_t *g, const double *ps, double H, bool *x, g } } - if (v != g->nv - 1) { + if (v != g->nv - 1) { // count the number of non-external sites that flip c->nv++; } } } - c->dH = n_bonds + H * n_h_bonds; - c->dM = n_h_bonds; - return c; } -double hh(double th) { - return (th - pow(th, 3) / 1.16951) * (1 - 0.222389 * pow(th, 2) - 0.043547 * pow(th, 4) - 0.014809 * pow(th, 6) - 0.007168 * pow(th, 8)); -} - 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); @@ -151,7 +126,8 @@ double *get_bond_probs(double T, double H, ising_state_t *s) { return ps; } -uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps) { +uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, + double *ps) { if (r == NULL) { r = gsl_rng_alloc(gsl_rng_mt19937); gsl_rng_set(r, jst_rand_seed()); @@ -161,10 +137,10 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps ps = get_bond_probs(T, H, s); } - cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r); + cluster_t *c = flip_cluster(s->g, ps, s->spins, r); - s->M += -2 * c->dM; - s->H += 2 * c->dH; + s->M += -2 * c->dHb; + s->H += 2 * (c->dJb + H * c->dHb); uint32_t n_flips = c->nv; @@ -172,4 +148,3 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r, double *ps return n_flips; } - -- cgit v1.2.3-54-g00ecf