2
\$\begingroup\$

I was looking into SOA to better utilize SIMD instructions. I implemented a generic wrapper object that lets you chose the amount of individual vectors each structure consists of. I am using xsimd and boost libraries if you're interested in testing it out.

Next step would be to make initializing the buffer easier, like an emplace method that allows forwarding arguments to the nested vectors or smth., here your feedback would come into play :)

I implemented the stl-style transform algorithm on the buffer that lets you transform a buffer of size N into a buffer of size M. Here I took the pre-implemented transform algorithm from xsimd as reference. As an example, I calculate the cross-product and dot-product of some vectors. First I transform a buffer of size 3 into a buffer of size 3, second from 3 to 1.

Also added structured binding support by adding specializations for tuple_size and tuple_element, see here boost:: combine, range-based for and structured bindings

#include <chrono>
#include <iostream>
#include <vector>

#include <xsimd/xsimd.hpp>
#include <xsimd/stl/algorithms.hpp>

#include <boost/iterator/zip_iterator.hpp>

namespace voxel
{
    template <typename T, size_t N>
    struct Buffer final
    {
        std::array<std::vector<T, XSIMD_DEFAULT_ALLOCATOR(T)>, N> values;

        [[nodiscard]]
        auto begin() noexcept
        {
            return std::apply([](auto&&... xs) noexcept
            {
                return boost::make_zip_iterator(boost::make_tuple(std::begin(std::forward<decltype(xs)>(xs))...));
            }, values);
        }

        [[nodiscard]]
        auto cbegin() const noexcept
        {
            return std::apply([](auto&&... xs) noexcept
            {
                return boost::make_zip_iterator(boost::make_tuple(std::cbegin(std::forward<decltype(xs)>(xs))...));
            }, values);
        }

        [[nodiscard]]
        auto end() noexcept
        {
            return std::apply([](auto&&... xs) noexcept
            {
                return boost::make_zip_iterator(boost::make_tuple(std::end(std::forward<decltype(xs)>(xs))...));
            }, values);
        }

        [[nodiscard]]
        auto cend() const noexcept
        {
            return std::apply([](auto&&... xs) noexcept
            {
                return boost::make_zip_iterator(boost::make_tuple(std::cend(std::forward<decltype(xs)>(xs))...));
            }, values);
        }
    };
}

namespace std
{
    template <typename T, typename U>
    struct tuple_size<boost::tuples::cons<T, U>>
        : boost::tuples::length<boost::tuples::cons<T, U>>
    {
    };

    template <size_t I, typename T, typename U>
    struct tuple_element<I, boost::tuples::cons<T, U>>
        : boost::tuples::element<I, boost::tuples::cons<T, U>>
    {
    };
}

namespace voxel
{
    namespace detail
    {
        template <typename T, std::size_t I>
        struct RepeatHelper final
        {
            using Type = T;
        };
    }

    template <typename T, typename I>
    struct Repeat;

    template <typename T, size_t... Is>
    struct Repeat<T, std::index_sequence<Is...>> final
    {
        using Type = std::tuple<typename detail::RepeatHelper<T, Is>::Type...>;
    };

    namespace detail
    {
        template <typename Input, typename Output, typename Function, size_t... Is, size_t... Js>
        void transform(Input begin, Input end, Output out, Function&& function, std::index_sequence<Is...>,
                       std::index_sequence<Js...>)
        {
            using InputValueType = std::decay_t<decltype(begin->template get<0>())>;
            using OutputValueType = std::decay_t<decltype(out->template get<0>())>;

            static_assert(std::is_same_v<InputValueType, OutputValueType>,
                "Input and Output have to be of the same type!");

            using InputTraits = xsimd::simd_traits<InputValueType>;
            using OutputTraits = xsimd::simd_traits<OutputValueType>;

            const auto distance = static_cast<std::size_t>(std::distance(begin, end));
            const auto inputSize = InputTraits::size;
            const auto outputSize = OutputTraits::size;

            const auto alignInput =
                xsimd::get_alignment_offset(&(begin->template get<0>()), distance, inputSize);
            const auto alignOutput =
                xsimd::get_alignment_offset(&(out->template get<0>()), distance, outputSize);

            const auto alignEnd = alignInput + ((distance - alignInput) & ~(inputSize - 1));

            for (auto i = 0u; i < alignInput; ++i)
            {
                *(out + i) = function((begin + i)->template get<Is>()...);
            }

            typename Repeat<typename InputTraits::type, std::index_sequence<Is...>>::Type batches;
            if (alignInput == alignOutput)
            {
                for (auto i = alignInput; i < alignEnd; i += inputSize)
                {
                    (xsimd::load_aligned(&(begin + i)->template get<Is>(), std::get<Is>(batches)), ...);
                    const auto result = function(std::get<Is>(batches)...);
                    (xsimd::store_aligned(&(out + i)->template get<Js>(), result.template get<Js>()), ...);
                }
            }
            else
            {
                for (auto i = alignInput; i < alignEnd; i += inputSize)
                {
                    (xsimd::load_aligned(&(begin + i)->template get<Is>(), std::get<Is>(batches)), ...);
                    const auto result = function(std::get<Is>(batches)...);
                    (xsimd::store_unaligned(&(out + i)->template get<Js>(), result.template get<Js>()), ...);
                }
            }

            for (auto i = alignEnd; i < distance; ++i)
            {
                *(out + i) = function((begin + i)->template get<Is>()...);
            }
        }
    }

    template <typename Input, typename Output, typename Function>
    void transform(Input begin, Input end, Output out, Function&& function)
    {
        using InputSize = boost::tuples::length<std::decay_t<decltype(*begin)>>;
        using OutputSize = boost::tuples::length<std::decay_t<decltype(*out)>>;
        return detail::transform(begin, end, out, std::forward<Function>(function),
                                 std::make_index_sequence<InputSize::value>{},
                                 std::make_index_sequence<OutputSize::value>{});
    }
}

int main()
{
    voxel::Buffer<float, 3u> values;
    values.values[0].reserve(100U);
    values.values[1].reserve(100U);
    values.values[2].reserve(100U);

    for (auto y = 0; y < 10; ++y)
    {
        for (auto x = 0; x < 10; ++x)
        {
            values.values[0].emplace_back(x);
            values.values[1].emplace_back(y);
            values.values[2].emplace_back(x + y);
        }
    }

    std::cout << "cross products:\n";
    voxel::Buffer<float, 3u> crossProducts;
    crossProducts.values[0].resize(100U);
    crossProducts.values[1].resize(100U);
    crossProducts.values[2].resize(100U);

    const auto runs = 10U;
    const auto xx = 29.0F;
    const auto yy = 34.0F;
    const auto zz = 23.0F;

    auto accumulated = 0.0;

    for (auto i = 0u; i < runs; ++i)
    {
        const auto t0 = std::chrono::high_resolution_clock::now();

        voxel::transform(values.cbegin(), values.cend(), crossProducts.begin(),
                         [=](const auto& x, const auto& y, const auto& z) noexcept
                         {
                             return boost::make_tuple(yy * z - zz * y, zz * x - xx * z, xx * y - yy * x);
                         });

        const auto elapsed = std::chrono::duration_cast<std::chrono::duration<double>>(
            std::chrono::high_resolution_clock::now() - t0).count();
        accumulated += elapsed;
        std::cout << i << '.' << ' ' << elapsed << 's' << '\n';
    }
    std::cout << "average. " << accumulated / static_cast<double>(runs) << 's' << '\n' << '\n';

    for (const auto [x, y, z] : crossProducts)
    {
        std::cout << '(' << x << ',' << y << ',' << z << ')' << '\n';
    }

    std::cout << "\ndot products:\n";
    voxel::Buffer<float, 1u> dotProducts;
    dotProducts.values[0].resize(100U);

    accumulated = 0.0;

    for (auto i = 0u; i < runs; ++i)
    {
        const auto t0 = std::chrono::high_resolution_clock::now();

        voxel::transform(values.cbegin(), values.cend(), dotProducts.begin(),
                         [=](const auto& x, const auto& y, const auto& z) noexcept
                         {
                             return boost::make_tuple(xx * x + yy * y + zz * z);
                         });

        const auto elapsed = std::chrono::duration_cast<std::chrono::duration<double>>(
            std::chrono::high_resolution_clock::now() - t0).count();
        accumulated += elapsed;
        std::cout << i << '.' << ' ' << elapsed << 's' << '\n';
    }
    std::cout << "average. " << accumulated / static_cast<double>(runs) << 's' << '\n' << '\n';

    for (const auto [x] : dotProducts)
    {
        std::cout << '(' << x << ')' << '\n';
    }
}

EDIT: In the meantime I worked on implementing a join function, that lets you join 2 buffers into a new one, although it seems to work, I haven't fully tested it yet. As there is no answer at this time, I am going to add it to this post, instead of making a new one:

namespace detail
{
    template <typename Input1, typename Input2, typename Output, typename Function, size_t... Is, size_t... Js,
              size_t... Ks>
    void join(Input1 begin1, Input1 end1, Input2 begin2, Input2 end2, Output out, Function&& function,
              std::index_sequence<Is...>, std::index_sequence<Js...>, std::index_sequence<Ks...>)
    {
        using Input1ValueType = std::decay_t<decltype(begin1->template get<0>())>;
        using Input2ValueType = std::decay_t<decltype(begin2->template get<0>())>;
        using OutputValueType = std::decay_t<decltype(out->template get<0>())>;

        static_assert(std::conjunction_v<std::is_same<Input1ValueType, Input2ValueType>, std::is_same<
                                             Input1ValueType, OutputValueType>>,
            "Input1, Input2 and Output have to be the same type!");

        using Input1Traits = xsimd::simd_traits<Input1ValueType>;
        using Input2Traits = xsimd::simd_traits<Input2ValueType>;
        using OutputTraits = xsimd::simd_traits<OutputValueType>;

        const auto distance = static_cast<std::size_t>(std::min(std::distance(begin1, end1),
                                                                std::distance(begin2, end2)));
        const auto input1Size = Input1Traits::size;
        const auto input2Size = Input2Traits::size;
        const auto outputSize = OutputTraits::size;

        const auto alignBegin1 =
            xsimd::get_alignment_offset(&(begin1->template get<0>()), distance, input1Size);

        const auto alignBegin2 =
            xsimd::get_alignment_offset(&(begin2->template get<0>()), distance, input2Size);

        const auto alignOutput =
            xsimd::get_alignment_offset(&(out->template get<0>()), distance, outputSize);

        const auto alignEnd =
            alignBegin1 + ((distance - alignBegin1) & ~(input1Size - 1));

        for (auto i = 0U; i < alignBegin1; ++i)
        {
            *(out + i) = function((begin1 + i)->template get<Is>()..., (begin2 + i)->template get<Js>()...);
        }

        typename Repeat<typename Input1Traits::type, std::index_sequence<Is...>>::Type input1Batches;
        typename Repeat<typename Input2Traits::type, std::index_sequence<Js...>>::Type input2Batches;
        if (alignBegin1 == alignOutput && alignBegin1 == alignBegin2)
        {
            for (auto i = alignBegin1; i < alignEnd; i += input1Size)
            {
                (xsimd::load_aligned(&(begin1 + i)->template get<Is>(), std::get<Is>(input1Batches)), ...);
                (xsimd::load_aligned(&(begin2 + i)->template get<Js>(), std::get<Js>(input2Batches)), ...);
                const auto result = function(std::get<Is>(input1Batches)..., std::get<Js>(input2Batches)...);
                (xsimd::store_aligned(&(out + i)->template get<Ks>(), result.template get<Ks>()), ...);
            }
        }
        else if (alignBegin1 == alignOutput && alignBegin1 != alignBegin2)
        {
            for (auto i = alignBegin1; i < alignEnd; i += input1Size)
            {
                (xsimd::load_aligned(&(begin1 + i)->template get<Is>(), std::get<Is>(input1Batches)), ...);
                (xsimd::load_unaligned(&(begin2 + i)->template get<Js>(), std::get<Js>(input2Batches)), ...);
                const auto result = function(std::get<Is>(input1Batches)..., std::get<Js>(input2Batches)...);
                (xsimd::store_aligned(&(out + i)->template get<Ks>(), result.template get<Ks>()), ...);
            }
        }
        else if (alignBegin1 != alignOutput && alignBegin1 == alignBegin2)
        {
            for (auto i = alignBegin1; i < alignEnd; i += input1Size)
            {
                (xsimd::load_aligned(&(begin1 + i)->template get<Is>(), std::get<Is>(input1Batches)), ...);
                (xsimd::load_aligned(&(begin2 + i)->template get<Js>(), std::get<Js>(input2Batches)), ...);
                const auto result = function(std::get<Is>(input1Batches)..., std::get<Js>(input2Batches)...);
                (xsimd::store_unaligned(&(out + i)->template get<Ks>(), result.template get<Ks>()), ...);
            }
        }
        else
        {
            for (auto i = alignBegin1; i < alignEnd; i += input1Size)
            {
                (xsimd::load_aligned(&(begin1 + i)->template get<Is>(), std::get<Is>(input1Batches)), ...);
                (xsimd::load_unaligned(&(begin2 + i)->template get<Js>(), std::get<Js>(input2Batches)), ...);
                const auto result = function(std::get<Is>(input1Batches)..., std::get<Js>(input2Batches)...);
                (xsimd::store_unaligned(&(out + i)->template get<Ks>(), result.template get<Ks>()), ...);
            }
        }

        for (auto i = alignEnd; i < distance; ++i)
        {
            *(out + i) = function((begin1 + i)->template get<Is>()..., (begin2 + i)->template get<Js>()...);
        }
    }
}

template <typename Input1, typename Input2, typename Output, typename Function>
void join(Input1 begin1, Input1 end1, Input2 begin2, Input2 end2, Output out, Function&& function)
{
    using Input1Size = boost::tuples::length<std::decay_t<decltype(*begin1)>>;
    using Input2Size = boost::tuples::length<std::decay_t<decltype(*begin2)>>;
    using OutputSize = boost::tuples::length<std::decay_t<decltype(*out)>>;
    detail::join(begin1, end1, begin2, end2, out, std::forward<Function>(function),
                 std::make_index_sequence<Input1Size::value>{},
                 std::make_index_sequence<Input2Size::value>{},
                 std::make_index_sequence<OutputSize::value>{});
}
\$\endgroup\$

1 Answer 1

2
\$\begingroup\$

The code is clear and conforms to modern C++ programming guidelines.

Here are some small points to consider:

final

The point of final is to help the compiler devirtualize function calls. Since your classes (Buffer, RepeatHelper, etc.) has no virtual member functions, final does not help the compiler optimize. It is hence unreasonable to prevent inheriting from these classes (which is handy in many situations). Therefore, remove the final.

Template metaprogramming

index_sequence + get is handy, but std::apply is usually more readable. For example, compare

for (auto i = 0u; i < alignInput; ++i)
{
    *(out + i) = function((begin + i)->template get<Is>()...);
}

with

std::transform(begin, begin + i, out, [](auto&& element) {
    return std::apply(function, std::forward<decltype(element)>(element));
});

(assuming that the tuple from boost implements the standard tuple protocol).

Since you are already using boost, you can use Boost.Mp11 mp_repeat_c to produce a tuple of N copies of the same type:

boost::mp11::mp_repeat_c<std::tuple<T>, N> // produces std::tuple</* N copies of T */>

Also, this one from join:

std::conjunction_v<std::is_same<Input1ValueType, Input2ValueType>,
                   std::is_same<Input1ValueType, OutputValueType>>

can be replaced by:

boost::mp11::mp_same<Input1ValueType, Input2ValueType, OutputValueType>::value
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.