summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/queue.c32
-rw-r--r--lib/wolff.h51
-rw-r--r--lib/wolff_tools.c179
3 files changed, 262 insertions, 0 deletions
diff --git a/lib/queue.c b/lib/queue.c
new file mode 100644
index 0000000..9a01741
--- /dev/null
+++ b/lib/queue.c
@@ -0,0 +1,32 @@
+
+#include "wolff.h"
+
+void stack_push(ll_t **q, uint32_t x) {
+ ll_t *nq = malloc(sizeof(ll_t));
+ nq->x = x;
+ nq->next = *q;
+
+ *q = nq;
+}
+
+uint32_t stack_pop(ll_t **q) {
+ ll_t *old_q = *q;
+
+ *q = old_q->next;
+ uint32_t x = old_q->x;
+
+ free(old_q);
+
+ return x;
+}
+
+bool stack_contains(const ll_t *q, uint32_t x) {
+ if (q == NULL) {
+ return false;
+ } else if (q->x == x) {
+ return true;
+ } else {
+ return stack_contains(q->next, x);
+ }
+}
+
diff --git a/lib/wolff.h b/lib/wolff.h
new file mode 100644
index 0000000..6beb03e
--- /dev/null
+++ b/lib/wolff.h
@@ -0,0 +1,51 @@
+
+#include <string.h>
+#include <math.h>
+#include <getopt.h>
+#include <float.h>
+#include <sys/types.h>
+#include <inttypes.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <stdbool.h>
+#include <assert.h>
+#include <omp.h>
+
+#include <jst/graph.h>
+#include <jst/rand.h>
+
+typedef struct {
+ graph_t *g;
+ bool *spins;
+ int32_t M;
+ double H;
+} ising_state_t;
+
+typedef struct ll_tag {
+ uint32_t x;
+ struct ll_tag *next;
+} ll_t;
+
+typedef struct {
+ int32_t nv;
+ double dH;
+} 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);
+
+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);
+
+int32_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
new file mode 100644
index 0000000..72cf9f7
--- /dev/null
+++ b/lib/wolff_tools.c
@@ -0,0 +1,179 @@
+
+#include "wolff.h"
+
+graph_t *graph_add_ext(const graph_t *g) {
+ graph_t *h = (graph_t *)calloc(1, sizeof(graph_t));
+
+ h->nv = g->nv + 1;
+ h->ne = g->ne + g->nv;
+
+ 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->vx = (double *)malloc(2 * h->nv * sizeof(double));
+ h->bq = (bool *)malloc(h->nv * sizeof(bool));
+
+ memcpy(h->ev, g->ev, 2 * g->ne * sizeof(uint32_t));
+ memcpy(h->vx, g->vx, 2 * g->nv * sizeof(double));
+ memcpy(h->bq, g->bq, g->nv * sizeof(bool));
+ h->vx[2 * g->nv] = -1;
+ h->vx[2 * g->nv + 1] = -0.5;
+ h->bq[g->nv] = false;
+
+ for (uint32_t i = 0; i < g->nv; i++) {
+ h->ev[2 * g->ne + 2 * i] = i;
+ h->ev[2 * g->ne + 2 * i + 1] = g->nv;
+ }
+
+ for (uint32_t i = 0; i < g->nv; i++) {
+ h->vei[i] = g->vei[i] + i;
+
+ for (uint32_t j = 0; j < g->vei[i + 1] - g->vei[i]; j++) {
+ h->ve[h->vei[i] + j] = g->ve[g->vei[i] + j];
+ }
+
+ h->ve[h->vei[i] + g->vei[i + 1] - g->vei[i]] = g->ne + i;
+ }
+
+ h->vei[g->nv] = g->vei[g->nv] + g->nv;
+ h->vei[g->nv + 1] = h->vei[g->nv] + g->nv;
+
+ for (uint32_t i = 0; i < g->nv; i++) {
+ h->ve[h->vei[g->nv] + i] = g->ne + i;
+ }
+
+ 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) {
+ uint32_t v0;
+ bool x0;
+ cluster_t *c;
+
+ v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex
+ x0 = x[v0]; // record its orientation
+
+ 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;
+ c->dH = 0;
+
+ while (stack != NULL) {
+ uint32_t v;
+ uint16_t nn;
+
+ v = stack_pop(&stack);
+ 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
+
+ for (uint16_t i = 0; i < nn; i++) {
+ uint32_t e, v1, v2, vn;
+
+ e = g->ve[g->vei[v] + i]; // select the ith bond connected to site
+ v1 = g->ev[2 * e];
+ v2 = g->ev[2 * e + 1];
+
+ vn = v == v1 ? v2 : v1; // distinguish neighboring site from site itself
+
+ if (x[vn] == x0) { // if the neighboring site matches the flipping cluster...
+ if (v1 == g->nv - 1 || v2 == g->nv - 1) {
+ c->dH += H;
+ } else {
+ c->dH += 1;
+ }
+
+ if (gsl_rng_uniform(r) < ps[e]) { // and with probability ps[e]...
+ stack_push(&stack, vn); // push the neighboring vertex to the stack
+ }
+ } else {
+ if (v1 == g->nv - 1 || v2 == g->nv - 1) {
+ c->dH -= H;
+ } else {
+ c->dH -= 1;
+ }
+ }
+ }
+
+ if (v != g->nv - 1) {
+ c->nv++;
+ }
+ }
+ }
+
+ if (x0) {
+ c->nv = -c->nv;
+ }
+
+ 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 * 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;
+}
+
+int32_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());
+ }
+
+ if (ps == NULL) {
+ ps = get_bond_probs(T, H, s);
+ }
+
+ cluster_t *c = flip_cluster(s->g, ps, H, s->spins, r);
+
+ s->M += 2 * c->nv;
+ s->H += c->dH;
+
+ int32_t n_flipped = c->nv;
+
+ free(c);
+
+ return n_flipped;
+}
+