1#include "reim4_elem.h"
  2
  3reim4_elem::reim4_elem(const double* re, const double* im) {
  4  for (uint64_t i = 0; i < 4; ++i) {
  5    value[i] = re[i];
  6    value[4 + i] = im[i];
  7  }
  8}
  9reim4_elem::reim4_elem(const double* layout) {
 10  for (uint64_t i = 0; i < 8; ++i) {
 11    value[i] = layout[i];
 12  }
 13}
 14reim4_elem::reim4_elem() {
 15  for (uint64_t i = 0; i < 8; ++i) {
 16    value[i] = 0.;
 17  }
 18}
 19void reim4_elem::save_re_im(double* re, double* im) const {
 20  for (uint64_t i = 0; i < 4; ++i) {
 21    re[i] = value[i];
 22    im[i] = value[4 + i];
 23  }
 24}
 25void reim4_elem::save_as(double* reim4) const {
 26  for (uint64_t i = 0; i < 8; ++i) {
 27    reim4[i] = value[i];
 28  }
 29}
 30reim4_elem reim4_elem::zero() { return reim4_elem(); }
 31
 32bool operator==(const reim4_elem& x, const reim4_elem& y) {
 33  for (uint64_t i = 0; i < 8; ++i) {
 34    if (x.value[i] != y.value[i]) return false;
 35  }
 36  return true;
 37}
 38
 39reim4_elem gaussian_reim4() {
 40  test_rng& gen = randgen();
 41  std::normal_distribution<double> dist(0, 1);
 42  reim4_elem res;
 43  for (uint64_t i = 0; i < 8; ++i) {
 44    res.value[i] = dist(gen);
 45  }
 46  return res;
 47}
 48
 49reim4_array_view::reim4_array_view(uint64_t size, double* data) : size(size), data(data) {}
 50reim4_elem reim4_array_view::get(uint64_t i) const {
 51  REQUIRE_DRAMATICALLY(i < size, "reim4 array overflow");
 52  return reim4_elem(data + 8 * i);
 53}
 54void reim4_array_view::set(uint64_t i, const reim4_elem& value) {
 55  REQUIRE_DRAMATICALLY(i < size, "reim4 array overflow");
 56  value.save_as(data + 8 * i);
 57}
 58
 59reim_view::reim_view(uint64_t m, double* data) : m(m), data(data) {}
 60reim4_elem reim_view::get_blk(uint64_t i) {
 61  REQUIRE_DRAMATICALLY(i < m / 4, "block overflow");
 62  return reim4_elem(data + 4 * i, data + m + 4 * i);
 63}
 64void reim_view::set_blk(uint64_t i, const reim4_elem& value) {
 65  REQUIRE_DRAMATICALLY(i < m / 4, "block overflow");
 66  value.save_re_im(data + 4 * i, data + m + 4 * i);
 67}
 68
 69reim_vector_view::reim_vector_view(uint64_t m, uint64_t nrows, double* data) : m(m), nrows(nrows), data(data) {}
 70reim_view reim_vector_view::row(uint64_t row) {
 71  REQUIRE_DRAMATICALLY(row < nrows, "row overflow");
 72  return reim_view(m, data + 2 * m * row);
 73}
 74
 75/** @brief addition */
 76reim4_elem operator+(const reim4_elem& x, const reim4_elem& y) {
 77  reim4_elem reps;
 78  for (uint64_t i = 0; i < 8; ++i) {
 79    reps.value[i] = x.value[i] + y.value[i];
 80  }
 81  return reps;
 82}
 83reim4_elem& operator+=(reim4_elem& x, const reim4_elem& y) {
 84  for (uint64_t i = 0; i < 8; ++i) {
 85    x.value[i] += y.value[i];
 86  }
 87  return x;
 88}
 89/** @brief subtraction */
 90reim4_elem operator-(const reim4_elem& x, const reim4_elem& y) {
 91  reim4_elem reps;
 92  for (uint64_t i = 0; i < 8; ++i) {
 93    reps.value[i] = x.value[i] + y.value[i];
 94  }
 95  return reps;
 96}
 97reim4_elem& operator-=(reim4_elem& x, const reim4_elem& y) {
 98  for (uint64_t i = 0; i < 8; ++i) {
 99    x.value[i] -= y.value[i];
100  }
101  return x;
102}
103/** @brief product */
104reim4_elem operator*(const reim4_elem& x, const reim4_elem& y) {
105  reim4_elem reps;
106  for (uint64_t i = 0; i < 4; ++i) {
107    double xre = x.value[i];
108    double yre = y.value[i];
109    double xim = x.value[i + 4];
110    double yim = y.value[i + 4];
111    reps.value[i] = xre * yre - xim * yim;
112    reps.value[i + 4] = xre * yim + xim * yre;
113  }
114  return reps;
115}
116/** @brief distance in infty norm */
117double infty_dist(const reim4_elem& x, const reim4_elem& y) {
118  double dist = 0;
119  for (uint64_t i = 0; i < 8; ++i) {
120    double d = fabs(x.value[i] - y.value[i]);
121    if (d > dist) dist = d;
122  }
123  return dist;
124}
125
126std::ostream& operator<<(std::ostream& out, const reim4_elem& x) {
127  out << "[\n";
128  for (uint64_t i = 0; i < 4; ++i) {
129    out << "  re=" << x.value[i] << ", im=" << x.value[i + 4] << "\n";
130  }
131  return out << "]";
132}
133
134reim4_matrix_view::reim4_matrix_view(uint64_t nrows, uint64_t ncols, double* data)
135    : nrows(nrows), ncols(ncols), data(data) {}
136reim4_elem reim4_matrix_view::get(uint64_t row, uint64_t col) const {
137  REQUIRE_DRAMATICALLY(row < nrows, "rows out of bounds" << row << " / " << nrows);
138  REQUIRE_DRAMATICALLY(col < ncols, "cols out of bounds" << col << " / " << ncols);
139  return reim4_elem(data + 8 * (row * ncols + col));
140}
141void reim4_matrix_view::set(uint64_t row, uint64_t col, const reim4_elem& value) {
142  REQUIRE_DRAMATICALLY(row < nrows, "rows out of bounds" << row << " / " << nrows);
143  REQUIRE_DRAMATICALLY(col < ncols, "cols out of bounds" << col << " / " << ncols);
144  value.save_as(data + 8 * (row * ncols + col));
145}