argmin + xsimd

argminmax.cpp

    using vT = xsimd::simd_type<T>;
    using iT = xsimd::as_integer_t<T>;
    static const size_t vN = vT::size;
    const long n = elts.size();
    if(n >= std::numeric_limits<iT>::max()) {
      return _argminmax_seq<Op>(elts, minmax_elts);
    }

    auto viter = types::vectorizer_nobroadcast::vbegin(elts),
         vend = types::vectorizer_nobroadcast::vend(elts);

    const long bound = std::distance(viter, vend);
    long minmax_index = -1;
    if (bound > 0) {
      auto vacc = *viter;
      iT iota[vN] = {0};
      for (long i = 0; i < vN; ++i)
        iota[i] = i;
      auto curr = xsimd::load_unaligned(iota);
      xsimd::simd_type<iT> indices = curr;
      xsimd::simd_type<iT> step{vN};

      for (++viter; viter != vend; ++viter) {
        curr += step;
        auto c = *viter;
        vacc = typename Op::op{}(vacc, c);
        auto mask = c == vacc;
        indices = xsimd::select(xsimd::bool_cast(mask), curr, indices);
      }

      alignas(sizeof(vT)) T stored[vN];
      vacc.store_aligned(&stored[0]);
      alignas(sizeof(vT)) long indexed[vN];
      indices.store_aligned(&indexed[0]);

      for (size_t j = 0; j < vN; ++j) {
        if (Op::value(stored[j], minmax_elts)) {
          minmax_elts = stored[j];
          minmax_index = indexed[j];
        }
      }
    }