1#include <gtest/gtest.h>
  2
  3#include <iostream>
  4#include <random>
  5
  6#include "../spqlios/reim4/reim4_arithmetic.h"
  7#include "test/testlib/reim4_elem.h"
  8
  9/// Actual tests
 10
 11typedef typeof(reim4_extract_1blk_from_reim_ref) reim4_extract_1blk_from_reim_f;
 12void test_reim4_extract_1blk_from_reim(reim4_extract_1blk_from_reim_f reim4_extract_1blk_from_reim) {
 13  static const uint64_t numtrials = 100;
 14  for (uint64_t m : {4, 8, 16, 1024, 4096, 32768}) {
 15    double* v = (double*)malloc(2 * m * sizeof(double));
 16    double* w = (double*)malloc(8 * sizeof(double));
 17    reim_view vv(m, v);
 18    for (uint64_t i = 0; i < numtrials; ++i) {
 19      reim4_elem el = gaussian_reim4();
 20      uint64_t blk = rand() % (m / 4);
 21      vv.set_blk(blk, el);
 22      reim4_extract_1blk_from_reim(m, blk, w, v);
 23      reim4_elem actual(w);
 24      ASSERT_EQ(el, actual);
 25    }
 26    free(v);
 27    free(w);
 28  }
 29}
 30
 31TEST(reim4_arithmetic, reim4_extract_1blk_from_reim_ref) {
 32  test_reim4_extract_1blk_from_reim(reim4_extract_1blk_from_reim_ref);
 33}
 34#ifdef __x86_64__
 35TEST(reim4_arithmetic, reim4_extract_1blk_from_reim_avx) {
 36  test_reim4_extract_1blk_from_reim(reim4_extract_1blk_from_reim_avx);
 37}
 38#endif
 39
 40typedef typeof(reim4_save_1blk_to_reim_ref) reim4_save_1blk_to_reim_f;
 41void test_reim4_save_1blk_to_reim(reim4_save_1blk_to_reim_f reim4_save_1blk_to_reim) {
 42  static const uint64_t numtrials = 100;
 43  for (uint64_t m : {4, 8, 16, 1024, 4096, 32768}) {
 44    double* v = (double*)malloc(2 * m * sizeof(double));
 45    double* w = (double*)malloc(8 * sizeof(double));
 46    reim_view vv(m, v);
 47    for (uint64_t i = 0; i < numtrials; ++i) {
 48      reim4_elem el = gaussian_reim4();
 49      el.save_as(w);
 50      uint64_t blk = rand() % (m / 4);
 51      reim4_save_1blk_to_reim_ref(m, blk, v, w);
 52      reim4_elem actual = vv.get_blk(blk);
 53      ASSERT_EQ(el, actual);
 54    }
 55    free(v);
 56    free(w);
 57  }
 58}
 59
 60TEST(reim4_arithmetic, reim4_save_1blk_to_reim_ref) { test_reim4_save_1blk_to_reim(reim4_save_1blk_to_reim_ref); }
 61#ifdef __x86_64__
 62TEST(reim4_arithmetic, reim4_save_1blk_to_reim_avx) { test_reim4_save_1blk_to_reim(reim4_save_1blk_to_reim_avx); }
 63#endif
 64
 65typedef typeof(reim4_extract_1blk_from_contiguous_reim_ref) reim4_extract_1blk_from_contiguous_reim_f;
 66void test_reim4_extract_1blk_from_contiguous_reim(
 67    reim4_extract_1blk_from_contiguous_reim_f reim4_extract_1blk_from_contiguous_reim) {
 68  static const uint64_t numtrials = 20;
 69  for (uint64_t m : {4, 8, 16, 1024, 4096, 32768}) {
 70    for (uint64_t nrows : {1, 2, 5, 128}) {
 71      double* v = (double*)malloc(2 * m * nrows * sizeof(double));
 72      double* w = (double*)malloc(8 * nrows * sizeof(double));
 73      reim_vector_view vv(m, nrows, v);
 74      reim4_array_view ww(nrows, w);
 75      for (uint64_t i = 0; i < numtrials; ++i) {
 76        uint64_t blk = rand() % (m / 4);
 77        for (uint64_t j = 0; j < nrows; ++j) {
 78          reim4_elem el = gaussian_reim4();
 79          vv.row(j).set_blk(blk, el);
 80        }
 81        reim4_extract_1blk_from_contiguous_reim_ref(m, nrows, blk, w, v);
 82        for (uint64_t j = 0; j < nrows; ++j) {
 83          reim4_elem el = vv.row(j).get_blk(blk);
 84          reim4_elem actual = ww.get(j);
 85          ASSERT_EQ(el, actual);
 86        }
 87      }
 88      free(v);
 89      free(w);
 90    }
 91  }
 92}
 93
 94TEST(reim4_arithmetic, reim4_extract_1blk_from_contiguous_reim_ref) {
 95  test_reim4_extract_1blk_from_contiguous_reim(reim4_extract_1blk_from_contiguous_reim_ref);
 96}
 97#ifdef __x86_64__
 98TEST(reim4_arithmetic, reim4_extract_1blk_from_contiguous_reim_avx) {
 99  test_reim4_extract_1blk_from_contiguous_reim(reim4_extract_1blk_from_contiguous_reim_avx);
100}
101#endif
102
103// test of basic arithmetic functions
104
105TEST(reim4_arithmetic, add) {
106  reim4_elem x = gaussian_reim4();
107  reim4_elem y = gaussian_reim4();
108  reim4_elem expect = x + y;
109  reim4_elem actual;
110  reim4_add(actual.value, x.value, y.value);
111  ASSERT_EQ(actual, expect);
112}
113
114TEST(reim4_arithmetic, mul) {
115  reim4_elem x = gaussian_reim4();
116  reim4_elem y = gaussian_reim4();
117  reim4_elem expect = x * y;
118  reim4_elem actual;
119  reim4_mul(actual.value, x.value, y.value);
120  ASSERT_EQ(actual, expect);
121}
122
123TEST(reim4_arithmetic, add_mul) {
124  reim4_elem x = gaussian_reim4();
125  reim4_elem y = gaussian_reim4();
126  reim4_elem z = gaussian_reim4();
127  reim4_elem expect = z;
128  reim4_elem actual = z;
129  expect += x * y;
130  reim4_add_mul(actual.value, x.value, y.value);
131  ASSERT_EQ(actual, expect) << infty_dist(expect, actual);
132}
133
134// test of dot products
135
136typedef typeof(reim4_vec_mat1col_product_ref) reim4_vec_mat1col_product_f;
137void test_reim4_vec_mat1col_product(reim4_vec_mat1col_product_f product) {
138  for (uint64_t ell : {1, 2, 5, 13, 69, 129}) {
139    std::vector<double> actual(8);
140    std::vector<double> a(ell * 8);
141    std::vector<double> b(ell * 8);
142    reim4_array_view va(ell, a.data());
143    reim4_array_view vb(ell, b.data());
144    reim4_array_view vactual(1, actual.data());
145    // initialize random values
146    for (uint64_t i = 0; i < ell; ++i) {
147      va.set(i, gaussian_reim4());
148      vb.set(i, gaussian_reim4());
149    }
150    // compute the mat1col product
151    reim4_elem expect;
152    for (uint64_t i = 0; i < ell; ++i) {
153      expect += va.get(i) * vb.get(i);
154    }
155    // compute the actual product
156    product(ell, actual.data(), a.data(), b.data());
157    // compare
158    ASSERT_LE(infty_dist(vactual.get(0), expect), 1e-10);
159  }
160}
161
162TEST(reim4_arithmetic, reim4_vec_mat1col_product_ref) { test_reim4_vec_mat1col_product(reim4_vec_mat1col_product_ref); }
163#ifdef __x86_64__
164TEST(reim4_arena, reim4_vec_mat1col_product_avx2) { test_reim4_vec_mat1col_product(reim4_vec_mat1col_product_avx2); }
165#endif
166
167typedef typeof(reim4_vec_mat2cols_product_ref) reim4_vec_mat2col_product_f;
168void test_reim4_vec_mat2cols_product(reim4_vec_mat2col_product_f product) {
169  for (uint64_t ell : {1, 2, 5, 13, 69, 129}) {
170    std::vector<double> actual(16);
171    std::vector<double> a(ell * 8);
172    std::vector<double> b(ell * 16);
173    reim4_array_view va(ell, a.data());
174    reim4_matrix_view vb(ell, 2, b.data());
175    reim4_array_view vactual(2, actual.data());
176    // initialize random values
177    for (uint64_t i = 0; i < ell; ++i) {
178      va.set(i, gaussian_reim4());
179      vb.set(i, 0, gaussian_reim4());
180      vb.set(i, 1, gaussian_reim4());
181    }
182    // compute the mat1col product
183    reim4_elem expect[2];
184    for (uint64_t i = 0; i < ell; ++i) {
185      expect[0] += va.get(i) * vb.get(i, 0);
186      expect[1] += va.get(i) * vb.get(i, 1);
187    }
188    // compute the actual product
189    product(ell, actual.data(), a.data(), b.data());
190    // compare
191    ASSERT_LE(infty_dist(vactual.get(0), expect[0]), 1e-10);
192    ASSERT_LE(infty_dist(vactual.get(1), expect[1]), 1e-10);
193  }
194}
195
196TEST(reim4_arithmetic, reim4_vec_mat2cols_product_ref) {
197  test_reim4_vec_mat2cols_product(reim4_vec_mat2cols_product_ref);
198}
199#ifdef __x86_64__
200TEST(reim4_arithmetic, reim4_vec_mat2cols_product_avx2) {
201  test_reim4_vec_mat2cols_product(reim4_vec_mat2cols_product_avx2);
202}
203#endif
204
205// for now, we do not need avx implementations,
206// so we will keep a single test function
207TEST(reim4_arithmetic, reim4_vec_convolution_ref) {
208  for (uint64_t sizea : {1, 2, 3, 5, 8}) {
209    for (uint64_t sizeb : {1, 3, 6, 9, 13}) {
210      std::vector<double> a(8 * sizea);
211      std::vector<double> b(8 * sizeb);
212      std::vector<double> expect(8 * (sizea + sizeb - 1));
213      std::vector<double> actual(8 * (sizea + sizeb - 1));
214      reim4_array_view va(sizea, a.data());
215      reim4_array_view vb(sizeb, b.data());
216      std::vector<reim4_elem> vexpect(sizea + sizeb + 3);
217      reim4_array_view vactual(sizea + sizeb - 1, actual.data());
218      for (uint64_t i = 0; i < sizea; ++i) {
219        va.set(i, gaussian_reim4());
220      }
221      for (uint64_t j = 0; j < sizeb; ++j) {
222        vb.set(j, gaussian_reim4());
223      }
224      // manual convolution
225      for (uint64_t i = 0; i < sizea; ++i) {
226        for (uint64_t j = 0; j < sizeb; ++j) {
227          vexpect[i + j] += va.get(i) * vb.get(j);
228        }
229      }
230      // partial convolution single coeff
231      for (uint64_t k = 0; k < sizea + sizeb + 3; ++k) {
232        double dest[8] = {0};
233        reim4_convolution_1coeff_ref(k, dest, a.data(), sizea, b.data(), sizeb);
234        ASSERT_LE(infty_dist(reim4_elem(dest), vexpect[k]), 1e-10);
235      }
236      // partial convolution dual coeff
237      for (uint64_t k = 0; k < sizea + sizeb + 2; ++k) {
238        double dest[16] = {0};
239        reim4_convolution_2coeff_ref(k, dest, a.data(), sizea, b.data(), sizeb);
240        ASSERT_LE(infty_dist(reim4_elem(dest), vexpect[k]), 1e-10);
241        ASSERT_LE(infty_dist(reim4_elem(dest + 8), vexpect[k + 1]), 1e-10);
242      }
243      // actual convolution
244      reim4_convolution_ref(actual.data(), sizea + sizeb - 1, 0, a.data(), sizea, b.data(), sizeb);
245      for (uint64_t k = 0; k < sizea + sizeb - 1; ++k) {
246        ASSERT_LE(infty_dist(vactual.get(k), vexpect[k]), 1e-10) << k;
247      }
248    }
249  }
250}
251
252EXPORT void reim4_convolution_ref(double* dest, uint64_t dest_size, uint64_t dest_offset, const double* a,
253                                  uint64_t sizea, const double* b, uint64_t sizeb);