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_add_1blk_to_reim_ref) reim4_add_1blk_to_reim_f;
 66void test_reim4_add_1blk_to_reim(reim4_add_1blk_to_reim_f reim4_add_1blk_to_reim) {
 67  static const uint64_t numtrials = 100;
 68  for (uint64_t m : {4, 8, 16, 1024, 4096, 32768}) {
 69    double* v = (double*)malloc(2 * m * sizeof(double));
 70    double* w1 = (double*)malloc(8 * sizeof(double));
 71    double* w2 = (double*)malloc(8 * sizeof(double));
 72    // double* tmp = (double*)malloc(8 * sizeof(double));
 73    reim_view vv(m, v);
 74    for (uint64_t i = 0; i < numtrials; ++i) {
 75      reim4_elem el1 = gaussian_reim4();
 76      reim4_elem el2 = gaussian_reim4();
 77
 78      // v[blk] = w1
 79      uint64_t blk = rand() % (m / 4);
 80      el1.save_as(w1);
 81      reim4_save_1blk_to_reim_ref(m, blk, v, w1);
 82      el2.save_as(w2);
 83
 84      // v[blk] += w2
 85      reim4_add_1blk_to_reim_ref(m, blk, v, w2);
 86
 87      reim4_elem w1w2sum = vv.get_blk(blk);
 88      reim4_elem expected_sum = reim4_elem::zero();
 89      reim4_add(expected_sum.value, w1, w2);
 90      ASSERT_EQ(expected_sum, w1w2sum);
 91    }
 92    free(v);
 93    free(w1);
 94    free(w2);
 95  }
 96}
 97
 98TEST(reim4_arithmetic, reim4_add_1blk_to_reim_ref) { test_reim4_add_1blk_to_reim(reim4_add_1blk_to_reim_ref); }
 99#ifdef __x86_64__
100TEST(reim4_arithmetic, reim4_add_1blk_to_reim_avx) { test_reim4_add_1blk_to_reim(reim4_add_1blk_to_reim_ref); }
101#endif
102
103typedef typeof(reim4_extract_1blk_from_contiguous_reim_ref) reim4_extract_1blk_from_contiguous_reim_f;
104void test_reim4_extract_1blk_from_contiguous_reim(
105    reim4_extract_1blk_from_contiguous_reim_f reim4_extract_1blk_from_contiguous_reim) {
106  static const uint64_t numtrials = 20;
107  for (uint64_t m : {4, 8, 16, 1024, 4096, 32768}) {
108    for (uint64_t nrows : {1, 2, 5, 128}) {
109      double* v = (double*)malloc(2 * m * nrows * sizeof(double));
110      double* w = (double*)malloc(8 * nrows * sizeof(double));
111      reim_vector_view vv(m, nrows, v);
112      reim4_array_view ww(nrows, w);
113      for (uint64_t i = 0; i < numtrials; ++i) {
114        uint64_t blk = rand() % (m / 4);
115        for (uint64_t j = 0; j < nrows; ++j) {
116          reim4_elem el = gaussian_reim4();
117          vv.row(j).set_blk(blk, el);
118        }
119        reim4_extract_1blk_from_contiguous_reim_ref(m, nrows, blk, w, v);
120        for (uint64_t j = 0; j < nrows; ++j) {
121          reim4_elem el = vv.row(j).get_blk(blk);
122          reim4_elem actual = ww.get(j);
123          ASSERT_EQ(el, actual);
124        }
125      }
126      free(v);
127      free(w);
128    }
129  }
130}
131
132TEST(reim4_arithmetic, reim4_extract_1blk_from_contiguous_reim_ref) {
133  test_reim4_extract_1blk_from_contiguous_reim(reim4_extract_1blk_from_contiguous_reim_ref);
134}
135#ifdef __x86_64__
136TEST(reim4_arithmetic, reim4_extract_1blk_from_contiguous_reim_avx) {
137  test_reim4_extract_1blk_from_contiguous_reim(reim4_extract_1blk_from_contiguous_reim_avx);
138}
139#endif
140
141// test of basic arithmetic functions
142
143TEST(reim4_arithmetic, add) {
144  reim4_elem x = gaussian_reim4();
145  reim4_elem y = gaussian_reim4();
146  reim4_elem expect = x + y;
147  reim4_elem actual;
148  reim4_add(actual.value, x.value, y.value);
149  ASSERT_EQ(actual, expect);
150}
151
152TEST(reim4_arithmetic, mul) {
153  reim4_elem x = gaussian_reim4();
154  reim4_elem y = gaussian_reim4();
155  reim4_elem expect = x * y;
156  reim4_elem actual;
157  reim4_mul(actual.value, x.value, y.value);
158  ASSERT_EQ(actual, expect);
159}
160
161TEST(reim4_arithmetic, add_mul) {
162  reim4_elem x = gaussian_reim4();
163  reim4_elem y = gaussian_reim4();
164  reim4_elem z = gaussian_reim4();
165  reim4_elem expect = z;
166  reim4_elem actual = z;
167  expect += x * y;
168  reim4_add_mul(actual.value, x.value, y.value);
169  ASSERT_EQ(actual, expect) << infty_dist(expect, actual);
170}
171
172// test of dot products
173
174typedef typeof(reim4_vec_mat1col_product_ref) reim4_vec_mat1col_product_f;
175void test_reim4_vec_mat1col_product(reim4_vec_mat1col_product_f product) {
176  for (uint64_t ell : {1, 2, 5, 13, 69, 129}) {
177    std::vector<double> actual(8);
178    std::vector<double> a(ell * 8);
179    std::vector<double> b(ell * 8);
180    reim4_array_view va(ell, a.data());
181    reim4_array_view vb(ell, b.data());
182    reim4_array_view vactual(1, actual.data());
183    // initialize random values
184    for (uint64_t i = 0; i < ell; ++i) {
185      va.set(i, gaussian_reim4());
186      vb.set(i, gaussian_reim4());
187    }
188    // compute the mat1col product
189    reim4_elem expect;
190    for (uint64_t i = 0; i < ell; ++i) {
191      expect += va.get(i) * vb.get(i);
192    }
193    // compute the actual product
194    product(ell, actual.data(), a.data(), b.data());
195    // compare
196    ASSERT_LE(infty_dist(vactual.get(0), expect), 1e-10);
197  }
198}
199
200TEST(reim4_arithmetic, reim4_vec_mat1col_product_ref) { test_reim4_vec_mat1col_product(reim4_vec_mat1col_product_ref); }
201#ifdef __x86_64__
202TEST(reim4_arena, reim4_vec_mat1col_product_avx2) { test_reim4_vec_mat1col_product(reim4_vec_mat1col_product_avx2); }
203#endif
204
205typedef typeof(reim4_vec_mat2cols_product_ref) reim4_vec_mat2col_product_f;
206void test_reim4_vec_mat2cols_product(reim4_vec_mat2col_product_f product) {
207  for (uint64_t ell : {1, 2, 5, 13, 69, 129}) {
208    std::vector<double> actual(16);
209    std::vector<double> a(ell * 8);
210    std::vector<double> b(ell * 16);
211    reim4_array_view va(ell, a.data());
212    reim4_matrix_view vb(ell, 2, b.data());
213    reim4_array_view vactual(2, actual.data());
214    // initialize random values
215    for (uint64_t i = 0; i < ell; ++i) {
216      va.set(i, gaussian_reim4());
217      vb.set(i, 0, gaussian_reim4());
218      vb.set(i, 1, gaussian_reim4());
219    }
220    // compute the mat1col product
221    reim4_elem expect[2];
222    for (uint64_t i = 0; i < ell; ++i) {
223      expect[0] += va.get(i) * vb.get(i, 0);
224      expect[1] += va.get(i) * vb.get(i, 1);
225    }
226    // compute the actual product
227    product(ell, actual.data(), a.data(), b.data());
228    // compare
229    ASSERT_LE(infty_dist(vactual.get(0), expect[0]), 1e-10);
230    ASSERT_LE(infty_dist(vactual.get(1), expect[1]), 1e-10);
231  }
232}
233
234TEST(reim4_arithmetic, reim4_vec_mat2cols_product_ref) {
235  test_reim4_vec_mat2cols_product(reim4_vec_mat2cols_product_ref);
236}
237#ifdef __x86_64__
238TEST(reim4_arithmetic, reim4_vec_mat2cols_product_avx2) {
239  test_reim4_vec_mat2cols_product(reim4_vec_mat2cols_product_avx2);
240}
241#endif
242
243// for now, we do not need avx implementations,
244// so we will keep a single test function
245TEST(reim4_arithmetic, reim4_vec_convolution_ref) {
246  for (uint64_t sizea : {1, 2, 3, 5, 8}) {
247    for (uint64_t sizeb : {1, 3, 6, 9, 13}) {
248      std::vector<double> a(8 * sizea);
249      std::vector<double> b(8 * sizeb);
250      std::vector<double> expect(8 * (sizea + sizeb - 1));
251      std::vector<double> actual(8 * (sizea + sizeb - 1));
252      reim4_array_view va(sizea, a.data());
253      reim4_array_view vb(sizeb, b.data());
254      std::vector<reim4_elem> vexpect(sizea + sizeb + 3);
255      reim4_array_view vactual(sizea + sizeb - 1, actual.data());
256      for (uint64_t i = 0; i < sizea; ++i) {
257        va.set(i, gaussian_reim4());
258      }
259      for (uint64_t j = 0; j < sizeb; ++j) {
260        vb.set(j, gaussian_reim4());
261      }
262      // manual convolution
263      for (uint64_t i = 0; i < sizea; ++i) {
264        for (uint64_t j = 0; j < sizeb; ++j) {
265          vexpect[i + j] += va.get(i) * vb.get(j);
266        }
267      }
268      // partial convolution single coeff
269      for (uint64_t k = 0; k < sizea + sizeb + 3; ++k) {
270        double dest[8] = {0};
271        reim4_convolution_1coeff_ref(k, dest, a.data(), sizea, b.data(), sizeb);
272        ASSERT_LE(infty_dist(reim4_elem(dest), vexpect[k]), 1e-10);
273      }
274      // partial convolution dual coeff
275      for (uint64_t k = 0; k < sizea + sizeb + 2; ++k) {
276        double dest[16] = {0};
277        reim4_convolution_2coeff_ref(k, dest, a.data(), sizea, b.data(), sizeb);
278        ASSERT_LE(infty_dist(reim4_elem(dest), vexpect[k]), 1e-10);
279        ASSERT_LE(infty_dist(reim4_elem(dest + 8), vexpect[k + 1]), 1e-10);
280      }
281      // actual convolution
282      reim4_convolution_ref(actual.data(), sizea + sizeb - 1, 0, a.data(), sizea, b.data(), sizeb);
283      for (uint64_t k = 0; k < sizea + sizeb - 1; ++k) {
284        ASSERT_LE(infty_dist(vactual.get(k), vexpect[k]), 1e-10) << k;
285      }
286    }
287  }
288}
289
290EXPORT void reim4_convolution_ref(double* dest, uint64_t dest_size, uint64_t dest_offset, const double* a,
291                                  uint64_t sizea, const double* b, uint64_t sizeb);