summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-30 06:32:51 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-30 06:32:51 -0400
commitfc242f52e835be85cc6030b6cae5619d18df7670 (patch)
treeeff9095238fbc9de89702b69863fca1f7bdcc29d
parent3223b527890e3090184384374f45a964cffa254a (diff)
downloadc++-fc242f52e835be85cc6030b6cae5619d18df7670.tar.gz
c++-fc242f52e835be85cc6030b6cae5619d18df7670.tar.bz2
c++-fc242f52e835be85cc6030b6cae5619d18df7670.zip
various changes
-rw-r--r--CMakeLists.txt5
-rw-r--r--lib/cluster.h68
-rw-r--r--lib/correlation.h4
-rw-r--r--lib/graph.cpp7
-rw-r--r--lib/graph.h1
-rw-r--r--lib/ising.h12
-rw-r--r--lib/state.h41
-rw-r--r--lib/vector.h6
-rw-r--r--src/wolff_On.cpp10
9 files changed, 105 insertions, 49 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fe564f7..568585f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -23,6 +23,7 @@ add_executable(wolff_7potts src/wolff_potts.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(wolff_3clock src/wolff_clock.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(wolff_5clock src/wolff_clock.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(wolff_planar src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES})
+add_executable(wolff_planar_2D src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(wolff_heisenberg src/wolff_On.cpp ${CPPSOURCES} ${CSOURCES})
add_executable(analyze_correlations src/analyze_correlations.cpp ${CPPSOURCES} ${CSOURCES})
@@ -32,6 +33,7 @@ SET_TARGET_PROPERTIES(wolff_7potts PROPERTIES COMPILE_FLAGS "-DPOTTSQ=7")
SET_TARGET_PROPERTIES(wolff_3clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=3")
SET_TARGET_PROPERTIES(wolff_5clock PROPERTIES COMPILE_FLAGS "-DPOTTSQ=5")
SET_TARGET_PROPERTIES(wolff_planar PROPERTIES COMPILE_FLAGS "-DN_COMP=2")
+SET_TARGET_PROPERTIES(wolff_planar_2D PROPERTIES COMPILE_FLAGS "-DN_COMP=2 -DDIMENSION=2")
SET_TARGET_PROPERTIES(wolff_heisenberg PROPERTIES COMPILE_FLAGS "-DN_COMP=3")
find_library(GSL REQUIRED NAMES gsl)
@@ -54,6 +56,7 @@ if (${GLUT} MATCHES "GLUT-NOTFOUND")
target_link_libraries(wolff_5clock cblas gsl m)
target_link_libraries(wolff_heisenberg cblas gsl m)
target_link_libraries(wolff_planar cblas gsl m)
+ target_link_libraries(wolff_planar_2D cblas gsl m)
else()
target_link_libraries(wolff_ising cblas gsl m glut GL GLU)
target_link_libraries(wolff_dgm cblas gsl m glut GL GLU)
@@ -65,6 +68,7 @@ else()
target_link_libraries(wolff_5clock cblas gsl m glut GL GLU)
target_link_libraries(wolff_heisenberg cblas gsl m glut GL GLU)
target_link_libraries(wolff_planar cblas gsl m glut GL GLU)
+ target_link_libraries(wolff_planar_2D cblas gsl m glut GL GLU)
target_compile_definitions(wolff_ising PUBLIC HAVE_GLUT)
target_compile_definitions(wolff_dgm PUBLIC HAVE_GLUT)
target_compile_definitions(wolff_cgm PUBLIC HAVE_GLUT)
@@ -74,6 +78,7 @@ else()
target_compile_definitions(wolff_3clock PUBLIC HAVE_GLUT)
target_compile_definitions(wolff_5clock PUBLIC HAVE_GLUT)
target_compile_definitions(wolff_planar PUBLIC HAVE_GLUT)
+ target_compile_definitions(wolff_planar_2D PUBLIC HAVE_GLUT)
target_compile_definitions(wolff_heisenberg PUBLIC HAVE_GLUT)
endif()
diff --git a/lib/cluster.h b/lib/cluster.h
index f3271b0..9dcf50d 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -25,78 +25,64 @@ void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) {
stack.pop();
if (!marks[v]) { // don't reprocess anyone we've already visited!
- X_t si_old, si_new;
- R_t R_old, R_new;
+ X_t si_new;
+ R_t R_new;
- R_old = s.R;
marks[v] = true;
bool v_is_ghost = (v == s.nv); // ghost site has the last index
if (v_is_ghost) {
- R_new = r.act(R_old); // if we are, then we're moving the transformation
+ R_new = r.act(s.R); // if we are, then we're moving the transformation
} else {
- si_old = s.spins[v];
- si_new = r.act(si_old); // otherwise, we're moving the spin at our site
+ si_new = r.act(s.spins[v]); // otherwise, we're moving the spin at our site
}
for (const v_t &vn : s.g.v_adj[v]) {
- X_t sj;
bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
- if (!vn_is_ghost) {
- sj = s.spins[vn];
- }
-
- double prob;
+ double dE, prob;
- if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond...
+ if (v_is_ghost || vn_is_ghost) { // this is a ghost-involved bond
X_t rs_old, rs_new;
v_t non_ghost;
+
if (vn_is_ghost) {
- rs_old = R_old.act_inverse(si_old);
- rs_new = R_old.act_inverse(si_new);
+ // if our neighbor is the ghost, the current site is a normal
+ // spin - rotate it back!
+ rs_old = s.R.act_inverse(s.spins[v]);
+ rs_new = s.R.act_inverse(si_new);
non_ghost = v;
} else {
- rs_old = R_old.act_inverse(sj);
- rs_new = R_new.act_inverse(sj);
+ /* if we're the ghost, we need to rotate our neighbor back in
+ both the old and new ways */
+ rs_old = s.R.act_inverse(s.spins[vn]);
+ rs_new = R_new.act_inverse(s.spins[vn]);
non_ghost = vn;
}
- double dE = s.H(rs_old) - s.H(rs_new);
+ dE = s.H(rs_old) - s.H(rs_new);
#ifdef FINITE_STATES
prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
-#else
- prob = 1.0 - exp(-dE / s.T);
#endif
- s.M -= rs_old;
- s.M += rs_new;
-
- s.E += dE;
-
- for (D_t i = 0; i < s.D; i++) {
- L_t x = (non_ghost / (v_t)pow(s.L, s.D - i - 1)) % s.L;
-
- s.ReF[i] -= rs_old * s.precomputed_cos[x];
- s.ReF[i] += rs_new * s.precomputed_cos[x];
-
- s.ImF[i] -= rs_old * s.precomputed_sin[x];
- s.ImF[i] += rs_new * s.precomputed_sin[x];
- }
- } else { // otherwise, we're at a perfectly normal bond!
- double dE = s.J(si_old, sj) - s.J(si_new, sj);
+ s.update_magnetization(rs_old, rs_new);
+ s.update_fourierZero(non_ghost, rs_old, rs_new);
+ } else { // this is a perfectly normal bond!
+ dE = s.J(s.spins[v], s.spins[vn]) - s.J(si_new, s.spins[vn]);
#ifdef FINITE_STATES
- prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
-#else
- prob = 1.0 - exp(-dE / s.T);
+ prob = J_probs[state_to_ind(s.spins[v])][state_to_ind(si_new)][state_to_ind(s.spins[vn])];
#endif
-
- s.E += dE;
}
+ s.update_energy(dE);
+
+#ifndef FINITE_STATES
+ prob = 1.0 - exp(-dE / s.T);
+#endif
+
if (gsl_rng_uniform(rand) < prob) {
stack.push(vn); // push the neighboring vertex to the stack
}
diff --git a/lib/correlation.h b/lib/correlation.h
index d08f713..29357a5 100644
--- a/lib/correlation.h
+++ b/lib/correlation.h
@@ -10,7 +10,11 @@ template <class R_t, class X_t>
double correlation_length(const state_t <R_t, X_t>& s) {
double total = 0;
+#ifdef DIMENSION
+ for (D_t j = 0; j < DIMENSION; j++) {
+#else
for (D_t j = 0; j < s.D; j++) {
+#endif
total += norm_squared(s.ReF[j]) + norm_squared(s.ImF[j]);
}
diff --git a/lib/graph.cpp b/lib/graph.cpp
index 8c97274..ca251f3 100644
--- a/lib/graph.cpp
+++ b/lib/graph.cpp
@@ -6,13 +6,17 @@ graph_t::graph_t(D_t D, L_t L) {
ne = D * nv;
v_adj.resize(nv);
+ coordinate.resize(nv);
for (std::vector<v_t> v_adj_i : v_adj) {
v_adj_i.reserve(2 * D);
}
for (v_t i = 0; i < nv; i++) {
+ coordinate[i].resize(D);
for (D_t j = 0; j < D; j++) {
+ coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
+
v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
v_adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
}
@@ -25,12 +29,15 @@ void graph_t::add_ext() {
}
v_adj.resize(nv + 1);
+ coordinate.resize(nv + 1);
v_adj[nv].reserve(nv);
for (v_t i = 0; i < nv; i++) {
v_adj[nv].push_back(i);
}
+ coordinate[nv].resize(coordinate[0].size());
+
ne += nv;
nv += 1;
}
diff --git a/lib/graph.h b/lib/graph.h
index a4732fb..de06924 100644
--- a/lib/graph.h
+++ b/lib/graph.h
@@ -13,6 +13,7 @@ class graph_t {
v_t ne;
v_t nv;
std::vector<std::vector<v_t>> v_adj;
+ std::vector<std::vector<double>> coordinate;
graph_t(D_t D, L_t L);
void add_ext();
diff --git a/lib/ising.h b/lib/ising.h
index 473a18a..aaa0c52 100644
--- a/lib/ising.h
+++ b/lib/ising.h
@@ -53,6 +53,18 @@ class ising_t {
return a;
}
}
+
+ inline int operator-(const ising_t &s) const {
+ if (x == s.x) {
+ return 0;
+ } else {
+ if (x) {
+ return -2;
+ } else {
+ return 2;
+ }
+ }
+ }
};
inline int& operator+=(int& M, const ising_t &s) {
diff --git a/lib/state.h b/lib/state.h
index 5abf65b..550d100 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -9,6 +9,10 @@
template <class R_t, class X_t>
class state_t {
+ private:
+ // updating fourier terms F requires many cos and sin calls, faster to do it beforehand.
+ std::vector<std::vector<double>> precomputed_cos;
+ std::vector<std::vector<double>> precomputed_sin;
public:
D_t D;
L_t L;
@@ -23,9 +27,6 @@ class state_t {
v_t last_cluster_size;
std::vector<typename X_t::F_t> ReF;
std::vector<typename X_t::F_t> ImF;
- // updating fourier terms F requires many cos and sin calls, faster to do it beforehand.
- std::vector<double> precomputed_cos;
- std::vector<double> precomputed_sin;
std::function <double(const X_t&, const X_t&)> J;
std::function <double(const X_t&)> H;
@@ -44,11 +45,35 @@ class state_t {
ReF[i] = spins[0] * 0.0;
ImF[i] = spins[0] * 0.0;
}
- precomputed_cos.resize(L);
- precomputed_sin.resize(L);
- for (L_t i = 0; i < L; i++) {
- precomputed_cos[i] = cos(2 * M_PI * (double)i / (double)L);
- precomputed_sin[i] = sin(2 * M_PI * (double)i / (double)L);
+ precomputed_cos.resize(nv);
+ precomputed_sin.resize(nv);
+ for (v_t i = 0; i < nv; i++) {
+ precomputed_cos[i].resize(D);
+ precomputed_sin[i].resize(D);
+ for (D_t j = 0; j < D; j++) {
+ precomputed_cos[i][j] = cos(2 * M_PI * g.coordinate[i][j] / (double)L);
+ precomputed_sin[i][j] = sin(2 * M_PI * g.coordinate[i][j] / (double)L);
+ }
+ }
+ }
+
+ void update_magnetization(const X_t& s_old, const X_t& s_new) {
+ M -= s_old;
+ M += s_new;
+ }
+
+ void update_energy(const double& dE) {
+ E += dE;
+ }
+
+ void update_fourierZero(v_t v, const X_t& s_old, const X_t& s_new) {
+#ifdef DIMENSION
+ for (D_t i = 0; i < DIMENSION; i++) {
+#else
+ for (D_t i = 0; i < D; i++) {
+#endif
+ ReF[i] += (s_new - s_old) * precomputed_cos[v][i];
+ ImF[i] += (s_new - s_old) * precomputed_sin[v][i];
}
}
};
diff --git a/lib/vector.h b/lib/vector.h
index 72633a8..7d0ee36 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -61,6 +61,12 @@ class vector_t : public std::array<T, q> {
return result;
}
+
+ inline vector_t<q, T> operator-(const vector_t<q, T>& v) const {
+ vector_t<q, T> diff = *this;
+ diff -= v;
+ return diff;
+ }
};
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 3fa5840..6cd8777 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -42,7 +42,11 @@ int main(int argc, char *argv[]) {
count_t N = (count_t)1e7;
+#ifdef DIMENSION
+ D_t D = DIMENSION;
+#else
D_t D = 2;
+#endif
L_t L = 128;
double T = 2.26918531421;
double *H_vec = (double *)calloc(MAX_Q, sizeof(double));
@@ -69,9 +73,15 @@ int main(int argc, char *argv[]) {
case 'N': // number of steps
N = (count_t)atof(optarg);
break;
+#ifdef DIMENSION
+ case 'D': // dimension
+ printf("Dimension was specified at compile time, you can't change it now!\n");
+ exit(EXIT_FAILURE);
+#else
case 'D': // dimension
D = atoi(optarg);
break;
+#endif
case 'L': // linear size
L = atoi(optarg);
break;