diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-07-10 12:37:02 -0400 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-07-10 12:37:02 -0400 |
commit | e53a4c09eb78e4c5a8365f1328a69ba7f9ff8992 (patch) | |
tree | 3c252af9ffafacab8392bf864270dcd034ed07ed /lib/orthogonal.h | |
parent | 609fb52b670d8ed74584a988b8c63da82d8d523b (diff) | |
parent | 1810103bc9ac4c9a8d432d113f5ca6eae6560fb4 (diff) | |
download | c++-e53a4c09eb78e4c5a8365f1328a69ba7f9ff8992.tar.gz c++-e53a4c09eb78e4c5a8365f1328a69ba7f9ff8992.tar.bz2 c++-e53a4c09eb78e4c5a8365f1328a69ba7f9ff8992.zip |
Merge branch 'master' of m5:/srv/git/wolff
Diffstat (limited to 'lib/orthogonal.h')
-rw-r--r-- | lib/orthogonal.h | 157 |
1 files changed, 148 insertions, 9 deletions
diff --git a/lib/orthogonal.h b/lib/orthogonal.h index 60d5f49..85f7a20 100644 --- a/lib/orthogonal.h +++ b/lib/orthogonal.h @@ -1,24 +1,163 @@ +#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 { bool is_reflection; T *x; }; + +template <q_t q, class T> +void init(orthogonal_t <q, T> *ptr) { + ptr->is_reflection = false; + 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.is_reflection = m.is_reflection; + + q_t size; + + if (m.is_reflection) { + size = q; + } else { + size = q * q; + } + + m_copy.x = (T *)calloc(size, sizeof(T)); + + for (q_t i = 0; i < size; 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)); + + if (m.is_reflection) { + double prod = 0; + for (q_t i = 0; i < q; i++) { + prod += v.x[i] * m.x[i]; + } + for (q_t i = 0; i < q; i++) { + v_rot.x[i] = v.x[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]; + } + } + } + + 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.is_reflection = false; + m2_rot.x = (T *)calloc(q * q, sizeof(T)); + + if (m1.is_reflection) { + for (q_t i = 0; i < q; i++) { + double akOki = 0; + + for (q_t k = 0; k < q; k++) { + akOki += m1.x[k] * m2.x[q * k + i]; + } + + for (q_t j = 0; j < q; j++) { + m2_rot.x[q * j + i] = m2.x[q * j + i] - 2 * akOki * m1.x[j]; + } + } + } else { + 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) { + if (m.is_reflection) { + 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]; + } + } + + 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) { + if (m1.is_reflection) { + return act(m1, m2); // reflections are their own inverse + } else { + 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) { + ptr->is_reflection = true; + ptr->x = (double *)calloc(q, sizeof(double)); -double *vector_rotate_inverse(q_t n, const double *rot, const double *vec); + double v2 = 0; -double vector_dot(q_t n, const double *v1, const double *v2); + for (q_t i = 0; i < q; i++) { + ptr->x[i] = gsl_ran_ugaussian(r); + v2 += ptr->x[i] * ptr->x[i]; + } -double *orthogonal_rotate(q_t n, const double *m1, const double *m2); + double mag_v = sqrt(v2); -double *gen_rot(gsl_rng *r, q_t n); + for (q_t i = 0; i < q; i++) { + ptr->x[i] /= mag_v; + } +} |