diff options
| -rw-r--r-- | lib/cluster.h | 6 | ||||
| -rw-r--r-- | lib/orthogonal.h | 22 |
2 files changed, 20 insertions, 8 deletions
diff --git a/lib/cluster.h b/lib/cluster.h index ccd8725..5f770ad 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -58,7 +58,7 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { ll_t *stack = NULL; // create a new stack stack_push(&stack, v0); // push the initial vertex to the stack - bool *marks = (bool *)calloc(state->nv, sizeof(bool)); + bool *marks = (bool *)calloc(state->nv + 1, sizeof(bool)); while (stack != NULL) { v_t v = stack_pop(&stack); @@ -67,10 +67,10 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { X_t si_old, si_new; R_t R_old, R_new; + R_old = state->R; marks[v] = true; if (v == state->nv) { - R_old = state->R; R_new = act (r, R_old); } else { si_old = state->spins[v]; @@ -94,7 +94,7 @@ v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) { if (is_ext) { X_t rs_old, rs_new; - if (vn == state->g->nv - 1) { + if (vn == state->nv) { rs_old = act_inverse (R_old, si_old); rs_new = act_inverse (R_old, si_new); } else { diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 6166223..d8ad33d 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -167,9 +167,9 @@ orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, const state_t <o template <q_t q> orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state_t <orthogonal_t <q, double>, vector_t <q, double>> *s, double epsilon) { - orthogonal_t <q, double> ptr; + orthogonal_t <q, double> m; vector_t <q, double> tmp_v; - ptr.is_reflection = true; + m.is_reflection = true; tmp_v.x = (double *)malloc(q * sizeof(double)); @@ -196,10 +196,22 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const state tmp_v.x[i] /= mag_v; } - vector_t <q, double> v = act(s->R, tmp_v); + vector_t <q, double> v = act_inverse(s->R, tmp_v); free(tmp_v.x); - ptr.x = v.x; + m.x = v.x; - return ptr; + v2 = 0; + + for (q_t i = 0; i < q; i++) { + v2 += m.x[i] * m.x[i]; + } + + mag_v = sqrt(v2); + + for (q_t i = 0; i < q; i++) { + m.x[i] /= mag_v; + } + + return m; } |