diff options
Diffstat (limited to 'lib/orthogonal.h')
-rw-r--r-- | lib/orthogonal.h | 27 |
1 files changed, 10 insertions, 17 deletions
diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 2a391ae..c0958f2 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -53,20 +53,19 @@ void free_spin (orthogonal_t <q, T> m) { template <q_t q, class T> vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) { vector_t <q, T> v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); if (m.is_reflection) { double prod = 0; for (q_t i = 0; i < q; i++) { - prod += v.x[i] * m.x[i]; + prod += v[i] * m.x[i]; } for (q_t i = 0; i < q; i++) { - v_rot.x[i] = v.x[i] - 2 * prod * m.x[i]; + v_rot[i] = v[i] - 2 * prod * m.x[i]; } } else { for (q_t i = 0; i < q; i++) { for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * i + j] * v.x[j]; + v_rot[i] += m.x[q * i + j] * v[j]; } } } @@ -112,11 +111,10 @@ vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) { return act(m, v); // reflections are their own inverse } else { vector_t <q, T> v_rot; - v_rot.x = (T *)calloc(q, sizeof(T)); for (q_t i = 0; i < q; i++) { for (q_t j = 0; j < q; j++) { - v_rot.x[i] += m.x[q * j + i] * v.x[j]; + v_rot[i] += m.x[q * j + i] * v[j]; } } @@ -176,33 +174,32 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q if (n > 1) { unsigned int rotation = gsl_rng_uniform_int(r, n); - v.x = (double *)malloc(q * sizeof(double)); double cosr = cos(2 * M_PI * rotation / (double)n / 2.0); double sinr = sin(2 * M_PI * rotation / (double)n / 2.0); - v.x[0] = v0.x[0] * cosr - v0.x[1] * sinr; - v.x[1] = v0.x[1] * cosr + v0.x[0] * sinr; + v[0] = v0[0] * cosr - v0[1] * sinr; + v[1] = v0[1] * cosr + v0[0] * sinr; for (q_t i = 2; i < q; i++) { - v.x[i] = v0.x[i]; + v[i] = v0[i]; } } else { - v.x = v0.x; + v = v0; } double m_dot_v = 0; for (q_t i = 0; i < q; i++) { m.x[i] = gsl_ran_ugaussian(r); - m_dot_v += m.x[i] * v.x[i]; + m_dot_v += m.x[i] * v[i]; } double v2 = 0; double factor = epsilon * gsl_ran_ugaussian(r); for (q_t i = 0; i < q; i++) { - m.x[i] = m.x[i] - m_dot_v * v.x[i] + factor * v.x[i]; + m.x[i] = m.x[i] - m_dot_v * v[i] + factor * v[i]; v2 += pow(m.x[i], 2); } @@ -212,10 +209,6 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q m.x[i] /= mag_v; } - if (n > 1) { - free(v.x); - } - return m; } |