<commit_before>#include <cstring>

#if defined _OPENMP
#include <omp.h>
#else
#define omp_get_thread_num()  0
#endif

#ifdef _MSC_VER
#pragma warning(disable:4127)
#endif
#define __DYN_GRID_SIZE
#include "common.h"
#include "metrix.h"
#include "OskarBinReader.h"
#include "aligned_malloc.h"

#define as256p(p) (reinterpret_cast<__m256d*>(p))
#define as256pc(p) (reinterpret_cast<const __m256d*>(p))

inline
void addGrids(
    Double4c dst[]
  , const Double4c srcs[]
  , int nthreads
  , int grid_pitch
  , int grid_size
  )
{
  int siz = grid_size*grid_pitch;
#pragma omp parallel for
  for (unsigned int i = 0; i < siz*sizeof(Double4c)/(256/8); i++) {
    __m256d sum = as256pc(srcs)[i];
    // __m256d sum = _mm256_loadu_pd(reinterpret_cast<const double*>(as256pc(srcs)+i));

    for (int g = 1; g < nthreads; g ++)
      sum = _mm256_add_pd(sum, as256pc(srcs + g * siz)[i]);

    as256p(dst)[i] = sum;
  }
}

// We could simply use pointer-to-function template
// but most C++ compilers seem to produce worse code
// in such a case. Thus we wrap it in a class.
template <
    int over
  , bool do_mirror
  , typename Inp
  > struct cvt {};

template <
    int over
  , bool do_mirror
  > struct cvt<over, do_mirror, Pregridded> {
  static void pre(double, double, Pregridded inp, Pregridded & outpr, int) {outpr = inp;}
};

template <
    int over
  , bool do_mirror
  > struct cvt<over, do_mirror, Double3> {
  static void pre(double scale, double wstep, Double3 inp, Pregridded & outpr, int grid_size) {
    pregridPoint<over, do_mirror>(scale, wstep, inp, outpr, grid_size);
  }
};

template <
    int over
  , bool is_half_gcf
  , bool use_permutations

  , typename Inp
  >
// grid must be initialized to 0s.
void gridKernel_scatter(
    double scale
  , double wstep
  , int baselines
  , const BlWMap permutations[/* baselines */]
  , Double4c grids[]
    // We have a [w_planes][over][over]-shaped array of pointers to
    // variable-sized gcf layers, but we precompute (in pregrid)
    // exact index into this array, thus we use plain pointer here
  , const complexd * gcf[]
  , const Inp _uvw[]
  , const Double4c _vis[]
  , int ts_ch
  , int grid_pitch
  , int grid_size
  ) {
  __ACC(Inp, uvw, ts_ch);
  int siz = grid_size*grid_pitch;
#pragma omp parallel
  {
    Double4c * _grid = grids + omp_get_thread_num() * siz;
    memset(_grid, 0, sizeof(Double4c) * siz);
    __ACC(Double4c, grid, grid_pitch);

#pragma omp for schedule(dynamic)
    for(int bl0 = 0; bl0 < baselines; bl0++) {
      int bl;
      if (use_permutations) bl = permutations[bl0].bl;
      else bl = bl0;
      int max_supp_here = get_supp(permutations[bl].wp);

      // Rotation
      // VLA requires "--std=gnu..." extension
      Double4c vis[ts_ch];
      int off = bl*ts_ch;
      for(int n=0; n<ts_ch; n++){
        #define __ROT_N_COPY(pol) vis[n].pol = rotw(_vis[off + n].pol, _uvw[off + n].w);
        __ROT_N_COPY(XX)
        __ROT_N_COPY(XY)
        __ROT_N_COPY(YX)
        __ROT_N_COPY(YY)
      }
      
      for (int su = 0; su < max_supp_here; su++) { // Moved from 2-levels below according to Romein
        for (int i = 0; i < ts_ch; i++) {
            Pregridded p;
            cvt<over, is_half_gcf, Inp>::pre(scale, wstep, uvw[bl][i], p, grid_size);

#ifdef __AVX__
          // We port Romein CPU code to doubles here (for MS2 we didn't)
          // vis0 covers XX and XY, vis1 -- YX and YY
          __m256d vis0, vis1;
          vis0 = _mm256_load_pd((const double *) &vis[i].XX);
          vis1 = _mm256_load_pd((const double *) &vis[i].YX);
          // vis0 = _mm256_loadu_pd((const double *) &vis[bl][i].XX);
          // vis1 = _mm256_loadu_pd((const double *) &vis[bl][i].YX);
#endif
          for (int sv = 0; sv < max_supp_here; sv++) {
            // Don't forget our u v are already translated by -max_supp_here/2
            int gsu, gsv;
            gsu = p.u + su;
            gsv = p.v + sv;

            complexd supportPixel;
            #define __layeroff su * max_supp_here + sv
            if (is_half_gcf) {
              int index;
              index = p.gcf_layer_index;
              // Negative index indicates that original w was mirrored
              // and we shall negate the index to obtain correct
              // offset *and* conjugate the result.
              if (index < 0) {
                supportPixel = conj(gcf[-index][__layeroff]);
              } else {
                supportPixel = gcf[index][__layeroff];
              }
            } else {
                supportPixel = gcf[p.gcf_layer_index][__layeroff];
            }

#ifdef __AVX__
            __m256d weight_r, weight_i;
            weight_r = _mm256_set1_pd(supportPixel.real());
            weight_i = _mm256_set1_pd(supportPixel.imag());

            /* _mm256_permute_pd(x, 5) swaps adjacent doubles pairs of x:
               d0, d1, d2, d3 goes to d1, d0, d3, d2
             */
            #define __DO(p,f)                                     \
            __m256d * gridptr##p = (__m256d *) &grid[gsu][gsv].f; \
            __m256d tr##p = _mm256_mul_pd(weight_r, vis##p);      \
            __m256d ti##p = _mm256_mul_pd(weight_i, vis##p);      \
            __m256d tp##p = _mm256_permute_pd(ti##p, 5);          \
            __m256d t##p = _mm256_addsub_pd(tr##p, tp##p);        \
            gridptr##p[p] = _mm256_add_pd(gridptr##p[p], t##p)
            //_mm256_storeu_pd(reinterpret_cast<double*>(gridptr##p + p), _mm256_add_pd(gridptr##p[p], t##p));

            __DO(0, XX);
            __DO(1, YX);
#else
            #define __GRID_POL(pol) grid[gsu][gsv].pol += vis[bl][i].pol * supportPixel
            __GRID_POL(XX);
            __GRID_POL(XY);
            __GRID_POL(YX);
            __GRID_POL(YY);
#endif
          }
        }
      }
    }
  }
}

<commit_msg>Fix non-avx CPU gridder variant.<commit_after>#include <cstring>

#if defined _OPENMP
#include <omp.h>
#else
#define omp_get_thread_num()  0
#endif

#ifdef _MSC_VER
#pragma warning(disable:4127)
#endif
#define __DYN_GRID_SIZE
#include "common.h"
#include "metrix.h"
#include "OskarBinReader.h"
#include "aligned_malloc.h"

#define as256p(p) (reinterpret_cast<__m256d*>(p))
#define as256pc(p) (reinterpret_cast<const __m256d*>(p))

inline
void addGrids(
    Double4c dst[]
  , const Double4c srcs[]
  , int nthreads
  , int grid_pitch
  , int grid_size
  )
{
  int siz = grid_size*grid_pitch;
#pragma omp parallel for
  for (unsigned int i = 0; i < siz*sizeof(Double4c)/(256/8); i++) {
    __m256d sum = as256pc(srcs)[i];
    // __m256d sum = _mm256_loadu_pd(reinterpret_cast<const double*>(as256pc(srcs)+i));

    for (int g = 1; g < nthreads; g ++)
      sum = _mm256_add_pd(sum, as256pc(srcs + g * siz)[i]);

    as256p(dst)[i] = sum;
  }
}

// We could simply use pointer-to-function template
// but most C++ compilers seem to produce worse code
// in such a case. Thus we wrap it in a class.
template <
    int over
  , bool do_mirror
  , typename Inp
  > struct cvt {};

template <
    int over
  , bool do_mirror
  > struct cvt<over, do_mirror, Pregridded> {
  static void pre(double, double, Pregridded inp, Pregridded & outpr, int) {outpr = inp;}
};

template <
    int over
  , bool do_mirror
  > struct cvt<over, do_mirror, Double3> {
  static void pre(double scale, double wstep, Double3 inp, Pregridded & outpr, int grid_size) {
    pregridPoint<over, do_mirror>(scale, wstep, inp, outpr, grid_size);
  }
};

template <
    int over
  , bool is_half_gcf
  , bool use_permutations

  , typename Inp
  >
// grid must be initialized to 0s.
void gridKernel_scatter(
    double scale
  , double wstep
  , int baselines
  , const BlWMap permutations[/* baselines */]
  , Double4c grids[]
    // We have a [w_planes][over][over]-shaped array of pointers to
    // variable-sized gcf layers, but we precompute (in pregrid)
    // exact index into this array, thus we use plain pointer here
  , const complexd * gcf[]
  , const Inp _uvw[]
  , const Double4c _vis[]
  , int ts_ch
  , int grid_pitch
  , int grid_size
  ) {
  __ACC(Inp, uvw, ts_ch);
  int siz = grid_size*grid_pitch;
#pragma omp parallel
  {
    Double4c * _grid = grids + omp_get_thread_num() * siz;
    memset(_grid, 0, sizeof(Double4c) * siz);
    __ACC(Double4c, grid, grid_pitch);

#pragma omp for schedule(dynamic)
    for(int bl0 = 0; bl0 < baselines; bl0++) {
      int bl;
      if (use_permutations) bl = permutations[bl0].bl;
      else bl = bl0;
      int max_supp_here = get_supp(permutations[bl].wp);

      // Rotation
      // VLA requires "--std=gnu..." extension
      Double4c vis[ts_ch];
      int off = bl*ts_ch;
      for(int n=0; n<ts_ch; n++){
        #define __ROT_N_COPY(pol) vis[n].pol = rotw(_vis[off + n].pol, _uvw[off + n].w);
        __ROT_N_COPY(XX)
        __ROT_N_COPY(XY)
        __ROT_N_COPY(YX)
        __ROT_N_COPY(YY)
      }
      
      for (int su = 0; su < max_supp_here; su++) { // Moved from 2-levels below according to Romein
        for (int i = 0; i < ts_ch; i++) {
            Pregridded p;
            cvt<over, is_half_gcf, Inp>::pre(scale, wstep, uvw[bl][i], p, grid_size);

#ifdef __AVX__
          // We port Romein CPU code to doubles here (for MS2 we didn't)
          // vis0 covers XX and XY, vis1 -- YX and YY
          __m256d vis0, vis1;
          vis0 = _mm256_load_pd((const double *) &vis[i].XX);
          vis1 = _mm256_load_pd((const double *) &vis[i].YX);
          // vis0 = _mm256_loadu_pd((const double *) &vis[bl][i].XX);
          // vis1 = _mm256_loadu_pd((const double *) &vis[bl][i].YX);
#endif
          for (int sv = 0; sv < max_supp_here; sv++) {
            // Don't forget our u v are already translated by -max_supp_here/2
            int gsu, gsv;
            gsu = p.u + su;
            gsv = p.v + sv;

            complexd supportPixel;
            #define __layeroff su * max_supp_here + sv
            if (is_half_gcf) {
              int index;
              index = p.gcf_layer_index;
              // Negative index indicates that original w was mirrored
              // and we shall negate the index to obtain correct
              // offset *and* conjugate the result.
              if (index < 0) {
                supportPixel = conj(gcf[-index][__layeroff]);
              } else {
                supportPixel = gcf[index][__layeroff];
              }
            } else {
                supportPixel = gcf[p.gcf_layer_index][__layeroff];
            }

#ifdef __AVX__
            __m256d weight_r, weight_i;
            weight_r = _mm256_set1_pd(supportPixel.real());
            weight_i = _mm256_set1_pd(supportPixel.imag());

            /* _mm256_permute_pd(x, 5) swaps adjacent doubles pairs of x:
               d0, d1, d2, d3 goes to d1, d0, d3, d2
             */
            #define __DO(p,f)                                     \
            __m256d * gridptr##p = (__m256d *) &grid[gsu][gsv].f; \
            __m256d tr##p = _mm256_mul_pd(weight_r, vis##p);      \
            __m256d ti##p = _mm256_mul_pd(weight_i, vis##p);      \
            __m256d tp##p = _mm256_permute_pd(ti##p, 5);          \
            __m256d t##p = _mm256_addsub_pd(tr##p, tp##p);        \
            gridptr##p[p] = _mm256_add_pd(gridptr##p[p], t##p)
            //_mm256_storeu_pd(reinterpret_cast<double*>(gridptr##p + p), _mm256_add_pd(gridptr##p[p], t##p));

            __DO(0, XX);
            __DO(1, YX);
#else
            #define __GRID_POL(pol) grid[gsu][gsv].pol += vis[i].pol * supportPixel
            __GRID_POL(XX);
            __GRID_POL(XY);
            __GRID_POL(YX);
            __GRID_POL(YY);
#endif
          }
        }
      }
    }
  }
}


{'fraction_non_alphanumeric': 0.12533310410040402, 'fraction_numerical': 0.027851800911200895, 'mean_word_length': 2.2488131806757887, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5136777', 'n_tokens_mistral': 4487, 'n_tokens_neox': 4231, 'n_words': 1046}