From 31f4244352b5e68eed770090419541d469f7f999 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 6 Jul 2018 14:51:29 -0400 Subject: split up some files --- lib/cluster.h | 177 +------------------------------------------------------ lib/orthogonal.c | 99 ------------------------------- lib/orthogonal.h | 116 +++++++++++++++++++++++++++++++++--- lib/vector.h | 72 ++++++++++++++++++++++ 4 files changed, 180 insertions(+), 284 deletions(-) delete mode 100644 lib/orthogonal.c create mode 100644 lib/vector.h diff --git a/lib/cluster.h b/lib/cluster.h index 29dd0cb..a20912e 100644 --- a/lib/cluster.h +++ b/lib/cluster.h @@ -20,6 +20,7 @@ #include "graph.h" #include "tree.h" #include "measurement.h" +#include "vector.h" #include "orthogonal.h" #include "dihedral.h" #include "dihinf.h" @@ -95,182 +96,6 @@ class state_t { } }; -template -struct vector_t { T *x; }; - -template -void init(vector_t *ptr) { - ptr->x = (T *)calloc(q, sizeof(T)); - - ptr->x[0] = (T)1; -} - -template -vector_t copy (vector_t v) { - vector_t v_copy; - - v_copy.x = (T *)calloc(q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - v_copy.x[i] = v.x[i]; - } - - return v_copy; -} - -template -void add (vector_t v1, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1.x[i] += v2.x[i]; - } -} - -template -void subtract (vector_t v1, vector_t v2) { - for (q_t i = 0; i < q; i++) { - v1.x[i] -= v2.x[i]; - } -} - -template -vector_t scalar_multiple(v_t a, vector_t v) { - vector_t multiple; - multiple.x = (T *)malloc(q * sizeof(T)); - for (q_t i = 0; i < q; i++) { - multiple.x[i] = a * v.x[i]; - } - - return multiple; -} - -template -T dot(vector_t v1, vector_t v2) { - T prod = 0; - - for (q_t i = 0; i < q; i++) { - prod += v1.x[i] * v2.x[i]; - } - - return prod; -} - -template -void free_spin (vector_t v) { - free(v.x); -} - -template -struct orthogonal_t { T *x; }; - -template -void init(orthogonal_t *ptr) { - ptr->x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - ptr->x[q * i + i] = (T)1; - } -} - -template -orthogonal_t copy (orthogonal_t m) { - orthogonal_t m_copy; - m_copy.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q * q; i++) { - m_copy.x[i] = m.x[i]; - } - - return m_copy; -} - -template -void free_spin (orthogonal_t m) { - free(m.x); -} - -template -vector_t act (orthogonal_t m, vector_t v) { - vector_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 * i + j] * v.x[j]; - } - } - - return v_rot; -} - - -template -orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { - orthogonal_t m2_rot; - m2_rot.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; - } - } - } - - return m2_rot; -} - -template -vector_t act_inverse (orthogonal_t m, vector_t v) { - vector_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]; - } - } - - return v_rot; -} - -template -orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { - orthogonal_t m2_rot; - m2_rot.x = (T *)calloc(q * q, sizeof(T)); - - for (q_t i = 0; i < q; i++) { - for (q_t j = 0; j < q; j++) { - for (q_t k = 0; k < q; k++) { - m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; - } - } - } - - return m2_rot; -} - -template -void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { - double *v = (double *)malloc(q * sizeof(double)); - double v2 = 0; - - for (q_t i = 0; i < q; i++) { - v[i] = gsl_ran_ugaussian(r); - v2 += v[i] * v[i]; - } - - ptr->x = (double *)calloc(q * q, sizeof(double)); - - for (q_t i = 0; i < q; i++) { - ptr->x[q * i + i] = 1.0; - for (q_t j = 0; j < q; j++) { - ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; - } - } - - free(v); -} - template v_t flip_cluster(state_t *state, v_t v0, R_t r, gsl_rng *rand) { v_t nv = 0; diff --git a/lib/orthogonal.c b/lib/orthogonal.c deleted file mode 100644 index 87569ae..0000000 --- a/lib/orthogonal.c +++ /dev/null @@ -1,99 +0,0 @@ - -#include "orthogonal.h" - -void vector_replace(q_t n, double *v1, const double *v2) { - // writes vector v2 of length n to memory located at v1 - for (q_t i = 0; i < n; i++) { - v1[i] = v2[i]; - } -} - -void vector_add(q_t n, double *v1, const double *v2) { - // adds vector v2 of length n to vector v1 - for (q_t i = 0; i < n; i++) { - v1[i] += v2[i]; - } -} - -void vector_subtract(q_t n, double *v1, const double *v2) { - // subtracts vector v2 of length n from vector v1 - for (q_t i = 0; i < n; i++) { - v1[i] -= v2[i]; - } -} - -double *vector_rotate(q_t n, const double *rot, const double *vec) { - // multiplies n by n rotation matrix rot to vector vec - double *rot_vec = (double *)malloc(n * sizeof(double)); - - double prod = 0.0; - for (q_t i = 0; i < n; i++) { - prod += rot[i] * vec[i]; - } - - for (q_t i = 0; i < n; i++) { - rot_vec[i] = vec[i] - 2 * prod * rot[i]; - } - - return rot_vec; -} - -double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) { - double *rot_vec = (double *)calloc(n, sizeof(double)); - - for (q_t i = 0; i < n; i++) { - for (q_t j = 0; j < n; j++) { - rot_vec[i] += rot[n * j + i] * vec[j]; - } - } - - return rot_vec; -} - -double vector_dot(q_t n, const double *v1, const double *v2) { - double dot = 0; - - for (q_t i = 0; i < n; i++) { - dot += v1[i] * v2[i]; - } - - return dot; -} - -double *orthogonal_rotate(q_t n, const double *r, const double *m) { - double *mul = (double *)calloc(n * n, sizeof(double)); - - for (q_t i = 0; i < n; i++) { - double akOki = 0; - - for (q_t k = 0; k < n; k++) { - akOki += r[k] * m[n * k + i]; - } - - for (q_t j = 0; j < n; j++) { - mul[n * j + i] = m[n * j + i] - 2 * akOki * r[j]; - } - } - - return mul; -} - -double *gen_rot(gsl_rng *r, q_t n) { - double *v = (double *)malloc(n * sizeof(double)); - - double v2 = 0; - - for (q_t i = 0; i < n; i++) { - v[i] = gsl_ran_ugaussian(r); - v2 += v[i] * v[i]; - } - - double magv = sqrt(v2); - - for (q_t i = 0; i < n; i++) { - v[i] /= magv; - } - - return v; -} - diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 60d5f49..0b2fdd5 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -1,24 +1,122 @@ +#pragma once + #include #include #include -#include +#include #include "types.h" -void vector_replace(q_t n, double *v1, const double *v2); +template +struct orthogonal_t { T *x; }; + +template +void init(orthogonal_t *ptr) { + ptr->x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = (T)1; + } +} + +template +orthogonal_t copy (orthogonal_t m) { + orthogonal_t m_copy; + m_copy.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q * q; i++) { + m_copy.x[i] = m.x[i]; + } + + return m_copy; +} + +template +void free_spin (orthogonal_t m) { + free(m.x); +} + +template +vector_t act (orthogonal_t m, vector_t v) { + vector_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 * i + j] * v.x[j]; + } + } + + return v_rot; +} + +template +orthogonal_t act (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k]; + } + } + } + + return m2_rot; +} + +template +vector_t act_inverse (orthogonal_t m, vector_t v) { + vector_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]; + } + } + + return v_rot; +} + +template +orthogonal_t act_inverse (orthogonal_t m1, orthogonal_t m2) { + orthogonal_t m2_rot; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); -void vector_add(q_t n, double *v1, const double *v2); + for (q_t i = 0; i < q; i++) { + for (q_t j = 0; j < q; j++) { + for (q_t k = 0; k < q; k++) { + m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k]; + } + } + } -void vector_subtract(q_t n, double *v1, const double *v2); + return m2_rot; +} -double *vector_rotate(q_t n, const double *rot, const double *vec); +template +void generate_rotation (gsl_rng *r, orthogonal_t *ptr) { + double *v = (double *)malloc(q * sizeof(double)); + double v2 = 0; -double *vector_rotate_inverse(q_t n, const double *rot, const double *vec); + for (q_t i = 0; i < q; i++) { + v[i] = gsl_ran_ugaussian(r); + v2 += v[i] * v[i]; + } -double vector_dot(q_t n, const double *v1, const double *v2); + ptr->x = (double *)calloc(q * q, sizeof(double)); + + for (q_t i = 0; i < q; i++) { + ptr->x[q * i + i] = 1.0; + for (q_t j = 0; j < q; j++) { + ptr->x[q * i + j] -= 2 * v[i] * v[j] / v2; + } + } -double *orthogonal_rotate(q_t n, const double *m1, const double *m2); + free(v); +} -double *gen_rot(gsl_rng *r, q_t n); diff --git a/lib/vector.h b/lib/vector.h new file mode 100644 index 0000000..c7f459c --- /dev/null +++ b/lib/vector.h @@ -0,0 +1,72 @@ + +#pragma once + +#include +#include + +#include "types.h" + +template +struct vector_t { T *x; }; + +template +void init(vector_t *ptr) { + ptr->x = (T *)calloc(q, sizeof(T)); + + ptr->x[0] = (T)1; +} + +template +vector_t copy (vector_t v) { + vector_t v_copy; + + v_copy.x = (T *)calloc(q, sizeof(T)); + + for (q_t i = 0; i < q; i++) { + v_copy.x[i] = v.x[i]; + } + + return v_copy; +} + +template +void add (vector_t v1, vector_t v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] += v2.x[i]; + } +} + +template +void subtract (vector_t v1, vector_t v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] -= v2.x[i]; + } +} + +template +vector_t scalar_multiple(v_t a, vector_t v) { + vector_t multiple; + multiple.x = (T *)malloc(q * sizeof(T)); + for (q_t i = 0; i < q; i++) { + multiple.x[i] = a * v.x[i]; + } + + return multiple; +} + +template +T dot(vector_t v1, vector_t v2) { + T prod = 0; + + for (q_t i = 0; i < q; i++) { + prod += v1.x[i] * v2.x[i]; + } + + return prod; +} + +template +void free_spin (vector_t v) { + free(v.x); +} + -- cgit v1.2.3-54-g00ecf