diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/cluster.h | 177 | ||||
| -rw-r--r-- | lib/orthogonal.c | 99 | ||||
| -rw-r--r-- | lib/orthogonal.h | 116 | ||||
| -rw-r--r-- | lib/vector.h | 72 |
4 files changed, 180 insertions, 284 deletions
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 <q_t q, class T> -struct vector_t { T *x; }; - -template <q_t q, class T> -void init(vector_t <q, T> *ptr) { - ptr->x = (T *)calloc(q, sizeof(T)); - - ptr->x[0] = (T)1; -} - -template <q_t q, class T> -vector_t <q, T> copy (vector_t <q, T> v) { - vector_t <q, 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 <q_t q, class T> -void add (vector_t <q, T> v1, vector_t <q, T> v2) { - for (q_t i = 0; i < q; i++) { - v1.x[i] += v2.x[i]; - } -} - -template <q_t q, class T> -void subtract (vector_t <q, T> v1, vector_t <q, T> v2) { - for (q_t i = 0; i < q; i++) { - v1.x[i] -= v2.x[i]; - } -} - -template <q_t q, class T> -vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) { - vector_t <q, 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 <q_t q, class T> -T dot(vector_t <q, T> v1, vector_t <q, T> v2) { - T prod = 0; - - for (q_t i = 0; i < q; i++) { - prod += v1.x[i] * v2.x[i]; - } - - return prod; -} - -template <q_t q, class T> -void free_spin (vector_t <q, T> v) { - free(v.x); -} - -template <q_t q, class T> -struct orthogonal_t { T *x; }; - -template <q_t q, class T> -void init(orthogonal_t <q, 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 <q_t q, class T> -orthogonal_t <q, T> copy (orthogonal_t <q, T> m) { - orthogonal_t <q, 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 <q_t q, class T> -void free_spin (orthogonal_t <q, T> m) { - free(m.x); -} - -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)); - - 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 <q_t q, class T> -orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) { - orthogonal_t <q, 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 <q_t q, class T> -vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) { - 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]; - } - } - - return v_rot; -} - -template <q_t q, class T> -orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) { - orthogonal_t <q, 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 <q_t q> -void generate_rotation (gsl_rng *r, orthogonal_t <q, double> *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 <class R_t, class X_t> v_t flip_cluster(state_t <R_t, X_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 <stdlib.h> #include <gsl/gsl_rng.h> #include <gsl/gsl_randist.h> -#include <math.h> +#include <cmath> #include "types.h" -void vector_replace(q_t n, double *v1, const double *v2); +template <q_t q, class T> +struct orthogonal_t { T *x; }; + +template <q_t q, class T> +void init(orthogonal_t <q, 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 <q_t q, class T> +orthogonal_t <q, T> copy (orthogonal_t <q, T> m) { + orthogonal_t <q, 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 <q_t q, class T> +void free_spin (orthogonal_t <q, T> m) { + free(m.x); +} + +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)); + + 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 <q_t q, class T> +orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) { + orthogonal_t <q, 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 <q_t q, class T> +vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) { + 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]; + } + } + + return v_rot; +} + +template <q_t q, class T> +orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) { + orthogonal_t <q, 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 <q_t q> +void generate_rotation (gsl_rng *r, orthogonal_t <q, double> *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 <stdlib.h> +#include <cmath> + +#include "types.h" + +template <q_t q, class T> +struct vector_t { T *x; }; + +template <q_t q, class T> +void init(vector_t <q, T> *ptr) { + ptr->x = (T *)calloc(q, sizeof(T)); + + ptr->x[0] = (T)1; +} + +template <q_t q, class T> +vector_t <q, T> copy (vector_t <q, T> v) { + vector_t <q, 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 <q_t q, class T> +void add (vector_t <q, T> v1, vector_t <q, T> v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] += v2.x[i]; + } +} + +template <q_t q, class T> +void subtract (vector_t <q, T> v1, vector_t <q, T> v2) { + for (q_t i = 0; i < q; i++) { + v1.x[i] -= v2.x[i]; + } +} + +template <q_t q, class T> +vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) { + vector_t <q, 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 <q_t q, class T> +T dot(vector_t <q, T> v1, vector_t <q, T> v2) { + T prod = 0; + + for (q_t i = 0; i < q; i++) { + prod += v1.x[i] * v2.x[i]; + } + + return prod; +} + +template <q_t q, class T> +void free_spin (vector_t <q, T> v) { + free(v.x); +} + |