summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/queue.c21
-rw-r--r--lib/queue.h18
-rw-r--r--lib/wolff.h44
-rw-r--r--lib/wolff_tools.c61
4 files changed, 63 insertions, 81 deletions
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 <inttypes.h>
+#include <stdbool.h>
+#include <stdlib.h>
+#include <string.h>
+
+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 <string.h>
-#include <math.h>
-#include <getopt.h>
+#pragma once
+
+#include <assert.h>
+#include <fftw3.h>
#include <float.h>
-#include <sys/types.h>
-#include <inttypes.h>
+#include <getopt.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
+#include <inttypes.h>
+#include <math.h>
#include <stdbool.h>
-#include <assert.h>
-#include <fftw3.h>
+#include <string.h>
+#include <sys/types.h>
#include <jst/graph.h>
#include <jst/rand.h>
+#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;
}
-