diff options
author | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-07-06 14:51:29 -0400 |
---|---|---|
committer | Jaron Kent-Dobias <jaron@kent-dobias.com> | 2018-07-06 14:51:29 -0400 |
commit | 31f4244352b5e68eed770090419541d469f7f999 (patch) | |
tree | c8fbe8259618833cab736593e263f0984bdb01ab /lib/orthogonal.h | |
parent | 2d8fcebf2f56efd1c3913ba49eaff6520ffdb33d (diff) | |
download | c++-31f4244352b5e68eed770090419541d469f7f999.tar.gz c++-31f4244352b5e68eed770090419541d469f7f999.tar.bz2 c++-31f4244352b5e68eed770090419541d469f7f999.zip |
split up some files
Diffstat (limited to 'lib/orthogonal.h')
-rw-r--r-- | lib/orthogonal.h | 116 |
1 files changed, 107 insertions, 9 deletions
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); |