1#include <memory.h>
  2
  3#include "zn_arithmetic_private.h"
  4
  5/** @brief size in bytes of a prepared matrix (for custom allocation) */
  6EXPORT uint64_t bytes_of_zn32_vmp_pmat(const MOD_Z* module,            // N
  7                                       uint64_t nrows, uint64_t ncols  // dimensions
  8) {
  9  return (nrows * ncols + 7) * sizeof(int32_t);
 10}
 11
 12/** @brief allocates a prepared matrix (release with delete_zn32_vmp_pmat) */
 13EXPORT ZN32_VMP_PMAT* new_zn32_vmp_pmat(const MOD_Z* module,  // N
 14                                        uint64_t nrows, uint64_t ncols) {
 15  return (ZN32_VMP_PMAT*)spqlios_alloc(bytes_of_zn32_vmp_pmat(module, nrows, ncols));
 16}
 17
 18/** @brief deletes a prepared matrix (release with free) */
 19EXPORT void delete_zn32_vmp_pmat(ZN32_VMP_PMAT* ptr) { spqlios_free(ptr); }
 20
 21/** @brief prepares a vmp matrix (contiguous row-major version) */
 22EXPORT void default_zn32_vmp_prepare_contiguous_ref(  //
 23    const MOD_Z* module,
 24    ZN32_VMP_PMAT* pmat,                                // output
 25    const int32_t* mat, uint64_t nrows, uint64_t ncols  // a
 26) {
 27  int32_t* const out = (int32_t*)pmat;
 28  const uint64_t nblk = ncols >> 5;
 29  const uint64_t ncols_rem = ncols & 31;
 30  const uint64_t final_elems = (8 - nrows * ncols) & 7;
 31  for (uint64_t blk = 0; blk < nblk; ++blk) {
 32    int32_t* outblk = out + blk * nrows * 32;
 33    const int32_t* srcblk = mat + blk * 32;
 34    for (uint64_t row = 0; row < nrows; ++row) {
 35      int32_t* dest = outblk + row * 32;
 36      const int32_t* src = srcblk + row * ncols;
 37      for (uint64_t i = 0; i < 32; ++i) {
 38        dest[i] = src[i];
 39      }
 40    }
 41  }
 42  // copy the last block if any
 43  if (ncols_rem) {
 44    int32_t* outblk = out + nblk * nrows * 32;
 45    const int32_t* srcblk = mat + nblk * 32;
 46    for (uint64_t row = 0; row < nrows; ++row) {
 47      int32_t* dest = outblk + row * ncols_rem;
 48      const int32_t* src = srcblk + row * ncols;
 49      for (uint64_t i = 0; i < ncols_rem; ++i) {
 50        dest[i] = src[i];
 51      }
 52    }
 53  }
 54  // zero-out the final elements that may be accessed
 55  if (final_elems) {
 56    int32_t* f = out + nrows * ncols;
 57    for (uint64_t i = 0; i < final_elems; ++i) {
 58      f[i] = 0;
 59    }
 60  }
 61}
 62
 63/** @brief prepares a vmp matrix (mat[row]+col*N points to the item) */
 64EXPORT void default_zn32_vmp_prepare_dblptr_ref(  //
 65    const MOD_Z* module,
 66    ZN32_VMP_PMAT* pmat,                                 // output
 67    const int32_t** mat, uint64_t nrows, uint64_t ncols  // a
 68) {
 69  for (uint64_t row_i = 0; row_i < nrows; ++row_i) {
 70    default_zn32_vmp_prepare_row_ref(module, pmat, mat[row_i], row_i, nrows, ncols);
 71  }
 72}
 73
 74/** @brief prepares the ith-row of a vmp matrix with nrows and ncols */
 75EXPORT void default_zn32_vmp_prepare_row_ref(  //
 76    const MOD_Z* module,
 77    ZN32_VMP_PMAT* pmat,                                                // output
 78    const int32_t* row, uint64_t row_i, uint64_t nrows, uint64_t ncols  // a
 79) {
 80  int32_t* const out = (int32_t*)pmat;
 81  const uint64_t nblk = ncols >> 5;
 82  const uint64_t ncols_rem = ncols & 31;
 83  const uint64_t final_elems = (row_i == nrows - 1) && (8 - nrows * ncols) & 7;
 84  for (uint64_t blk = 0; blk < nblk; ++blk) {
 85    int32_t* outblk = out + blk * nrows * 32;
 86    int32_t* dest = outblk + row_i * 32;
 87    const int32_t* src = row + blk * 32;
 88    for (uint64_t i = 0; i < 32; ++i) {
 89      dest[i] = src[i];
 90    }
 91  }
 92  // copy the last block if any
 93  if (ncols_rem) {
 94    int32_t* outblk = out + nblk * nrows * 32;
 95    int32_t* dest = outblk + row_i * ncols_rem;
 96    const int32_t* src = row + nblk * 32;
 97    for (uint64_t i = 0; i < ncols_rem; ++i) {
 98      dest[i] = src[i];
 99    }
100  }
101  // zero-out the final elements that may be accessed
102  if (final_elems) {
103    int32_t* f = out + nrows * ncols;
104    for (uint64_t i = 0; i < final_elems; ++i) {
105      f[i] = 0;
106    }
107  }
108}
109
110#if 0
111
112#define IMPL_zn32_vec_ixxx_matyyycols_ref(NCOLS) \
113  memset(res, 0, NCOLS * sizeof(int32_t));       \
114  for (uint64_t row = 0; row < nrows; ++row) {   \
115    int32_t ai = a[row];                         \
116    const int32_t* bb = b + row * b_sl;          \
117    for (uint64_t i = 0; i < NCOLS; ++i) {       \
118      res[i] += ai * bb[i];                      \
119    }                                            \
120  }
121
122#define IMPL_zn32_vec_ixxx_mat8cols_ref() IMPL_zn32_vec_ixxx_matyyycols_ref(8)
123#define IMPL_zn32_vec_ixxx_mat16cols_ref() IMPL_zn32_vec_ixxx_matyyycols_ref(16)
124#define IMPL_zn32_vec_ixxx_mat24cols_ref() IMPL_zn32_vec_ixxx_matyyycols_ref(24)
125#define IMPL_zn32_vec_ixxx_mat32cols_ref() IMPL_zn32_vec_ixxx_matyyycols_ref(32)
126
127void zn32_vec_i8_mat32cols_ref(uint64_t nrows, int32_t* res, const int8_t* a, const int32_t* b, uint64_t b_sl) {
128  IMPL_zn32_vec_ixxx_mat32cols_ref()
129}
130void zn32_vec_i16_mat32cols_ref(uint64_t nrows, int32_t* res, const int16_t* a, const int32_t* b, uint64_t b_sl) {
131  IMPL_zn32_vec_ixxx_mat32cols_ref()
132}
133
134void zn32_vec_i32_mat32cols_ref(uint64_t nrows, int32_t* res, const int32_t* a, const int32_t* b, uint64_t b_sl) {
135  IMPL_zn32_vec_ixxx_mat32cols_ref()
136}
137void zn32_vec_i32_mat24cols_ref(uint64_t nrows, int32_t* res, const int32_t* a, const int32_t* b, uint64_t b_sl) {
138  IMPL_zn32_vec_ixxx_mat24cols_ref()
139}
140void zn32_vec_i32_mat16cols_ref(uint64_t nrows, int32_t* res, const int32_t* a, const int32_t* b, uint64_t b_sl) {
141  IMPL_zn32_vec_ixxx_mat16cols_ref()
142}
143void zn32_vec_i32_mat8cols_ref(uint64_t nrows, int32_t* res, const int32_t* a, const int32_t* b, uint64_t b_sl) {
144  IMPL_zn32_vec_ixxx_mat8cols_ref()
145}
146typedef void (*zn32_vec_i32_mat8kcols_ref_f)(uint64_t nrows,                  //
147                                             int32_t* res,                    //
148                                             const int32_t* a,                //
149                                             const int32_t* b, uint64_t b_sl  //
150);
151zn32_vec_i32_mat8kcols_ref_f zn32_vec_i32_mat8kcols_ref[4] = {  //
152    zn32_vec_i32_mat8cols_ref, zn32_vec_i32_mat16cols_ref,      //
153    zn32_vec_i32_mat24cols_ref, zn32_vec_i32_mat32cols_ref};
154
155/** @brief applies a vmp product (int32_t* input) */
156EXPORT void default_zn32_vmp_apply_i32_ref(const MOD_Z* module,                //
157                                           int32_t* res, uint64_t res_size,    //
158                                           const int32_t* a, uint64_t a_size,  //
159                                           const ZN32_VMP_PMAT* pmat, uint64_t nrows, uint64_t ncols) {
160  const uint64_t rows = a_size < nrows ? a_size : nrows;
161  const uint64_t cols = res_size < ncols ? res_size : ncols;
162  const uint64_t ncolblk = cols >> 5;
163  const uint64_t ncolrem = cols & 31;
164  // copy the first full blocks
165  const uint32_t full_blk_size = nrows * 32;
166  const int32_t* mat = (int32_t*)pmat;
167  int32_t* rr = res;
168  for (uint64_t blk = 0;  //
169       blk < ncolblk;     //
170       ++blk, mat += full_blk_size, rr += 32) {
171    zn32_vec_i32_mat32cols_ref(rows, rr, a, mat, 32);
172  }
173  // last block
174  if (ncolrem) {
175    uint64_t orig_rem = ncols - (ncolblk << 5);
176    uint64_t b_sl = orig_rem >= 32 ? 32 : orig_rem;
177    int32_t tmp[32];
178    zn32_vec_i32_mat8kcols_ref[(ncolrem - 1) >> 3](rows, tmp, a, mat, b_sl);
179    memcpy(rr, tmp, ncolrem * sizeof(int32_t));
180  }
181  // trailing bytes
182  memset(res + cols, 0, (res_size - cols) * sizeof(int32_t));
183}
184
185#endif