CE_ArraySortImpl.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * NAME: CE_Array.h ( CE Library, C++)
 *
 * COMMENTS: UT_Array style class on GPU.
 */

#ifndef __CE_ArraySortImpl__
#define __CE_ArraySortImpl__

#include "CE_ArrayView.h"

#include <UT/UT_Swap.h>
#include <UT/UT_Tracing.h>
#include <UT/UT_WorkBuffer.h>

// Buffer for auxiliary allocation.
// Necessary because CE_ArrayView cannot depend on CE_BufferDevice.
template <typename T>
template <typename V>
class CE_ArrayView<T>::Buffer
{
public:
    exint size() { return mySize; }
    const cl::Buffer &buffer() { return myBuffer; }
    CE_ArrayView<V> view() { return CE_ArrayView<V>(myBuffer, mySize); }

    Buffer(exint size)
    {
        if (size > 0)
        {
            CE_Context *context = CE_Context::getContext();
            myBuffer = context->allocBuffer(size * sizeof(V));
            mySize = size;
        }
        else
        {
            mySize = 0;
        }
    }
    ~Buffer()
    {
        CE_Context *context = CE_Context::getContext();
        context->releaseBuffer(std::move(myBuffer));
    }

private:
    cl::Buffer  myBuffer;
    exint       mySize;
};

// Appends a type name we use as part of type defines in array.cl
template <typename T>
template <typename V>
void
CE_ArrayView<T>::appendElemType(UT_WorkBuffer &wb)
{
    if constexpr (std::is_same<V, uint8>::value)
        wb.append("UCHAR");
    if constexpr (std::is_same<V, int8>::value)
        wb.append("CHAR");
    if constexpr (std::is_same<V, uint16>::value)
        wb.append("USHORT");
    if constexpr (std::is_same<V, int16>::value)
        wb.append("SHORT");
    if constexpr (std::is_same<V, int32>::value)
        wb.append("INT");
    if constexpr (std::is_same<V, uint32>::value)
        wb.append("UINT");
    if constexpr (std::is_same<V, int64>::value)
        wb.append("EXINT");
    if constexpr (std::is_same<V, uint64>::value)
        wb.append("ULONG");
    if constexpr (std::is_same<V, fpreal16>::value)
        wb.append("HALF");
    if constexpr (std::is_same<V, fpreal32>::value)
        wb.append("FLOAT");
    if constexpr (std::is_same<V, fpreal64>::value)
        wb.append("DOUBLE");
}

template <typename T>
template <typename V>
void
CE_ArrayView<T>::convertFrom(const CE_ArrayView<V> &src,
                             int src_tuplesize, int dst_tuplesize,
                             exint src_offset, exint dst_offset,
                             exint nelements, T default_value)
{
    UT_ASSERT(src_tuplesize >= 1 && dst_tuplesize >= 1);
    if (!(src_tuplesize >= 1 && dst_tuplesize >= 1))
        return;

    exint src_nelem = src.size() / src_tuplesize;
    exint dst_nelem = size() / dst_tuplesize;
    if (nelements < 0)
        nelements = SYSmin(src_nelem - src_offset, dst_nelem - dst_offset);

    UT_ASSERT(src_offset >= 0 && dst_offset >= 0 &&
              (src_nelem - src_offset) >= nelements &&
              (dst_nelem - dst_offset) >= nelements);

    if (nelements > 0)
    {
        UT_WorkBuffer wb;
        wb.append(" -D CEARRAY_VALUE_");
        appendElemType<V>(wb);
        const char *opt = wb.buffer();
        cl::Kernel k = loadKernel("convertFrom", opt);

        CE_Context *context = CE_Context::getContext();
        cl::NDRange global_range, local_range;
        context->get1DRanges(k, nelements, global_range, local_range);
        cl::KernelFunctor convert_from = k.bind(context->getQueue(),
                                                global_range, local_range);

        convert_from(src.buffer(), buffer(),
                     src_tuplesize, dst_tuplesize,
                     src_offset, dst_offset,
                     nelements, scalarKernelArg(default_value));
    }
}

template <typename T>
template <typename V, typename I>
void
CE_ArrayView<T>::convertFromIndices(const CE_ArrayView<V> &src,
                                    const CE_ArrayView<I> &indices,
                                    int src_tuplesize, int dst_tuplesize,
                                    exint dst_offset,
                                    T default_value)
{
    SYS_STATIC_ASSERT(SYS_IsIntegral_v<I>);
    UT_ASSERT(src_tuplesize >= 1 && dst_tuplesize >= 1);
    if (!(src_tuplesize >= 1 && dst_tuplesize >= 1))
        return;
    
    SYS_MAYBE_UNUSED exint src_nelem = src.size() / src_tuplesize;
    SYS_MAYBE_UNUSED exint dst_nelem = size() / dst_tuplesize;
    UT_ASSERT((src_nelem - 1) <= std::numeric_limits<I>::max());
    UT_ASSERT(dst_offset >= 0 && (dst_nelem - dst_offset) >= indices.size());
    
    if (size() > 0)
    {
        UT_WorkBuffer wb;
        wb.append(" -D CEARRAY_VALUE_");
        appendElemType<V>(wb);
        wb.append(" -D CEARRAY_INDEX_");
        appendElemType<I>(wb);
        const char *opt = wb.buffer();
        cl::Kernel k = loadKernel("convertFromIndices", opt);

        CE_Context *context = CE_Context::getContext();
        cl::NDRange global_range, local_range;
        context->get1DRanges(k, indices.size(), global_range, local_range);
        cl::KernelFunctor convert_from = k.bind(context->getQueue(),
                                                global_range, local_range);

        convert_from(src.buffer(), indices.buffer(), buffer(),
                     src_tuplesize, dst_tuplesize, dst_offset,
                     indices.size(), scalarKernelArg(default_value));
    }
}

template <typename T>
template <typename I>
void
CE_ArrayView<T>::reorder(const CE_ArrayView<I> &order)
{
    // Index array must be an integer type with sufficient limits given size().
    SYS_STATIC_ASSERT(SYS_IsIntegral_v<I>);
    UT_ASSERT((size() - 1) <= std::numeric_limits<I>::max());
    UT_ASSERT(order.size() == size());
    if (!size()) return;

    UT_WorkBuffer wb;
    wb.append(" -D CEARRAY_INDEX_");
    appendElemType<I>(wb);
    const char *opt = wb.buffer();
    cl::Kernel k = loadKernel("buffer_reorder", opt);

    CE_Context *context = CE_Context::getContext();
    cl::NDRange global_range, local_range;
    context->get1DRanges(k, size(), global_range, local_range);
    cl::KernelFunctor op = k.bind(context->getQueue(),
                                  global_range, local_range);

    CE_BufferDevice<T> reordered(size());
    int elemsize = sizeof(T);
    op(buffer(), reordered.buffer(), size(), elemsize, order.buffer());

    // Copy reordered result to our buffer.
    context->copyBuffer(reordered.buffer(), buffer(), size() * sizeof(T));
}

template <typename T>
template <typename V>
void
CE_ArrayView<T>::sortInternal(
        CE_ArrayView<V> &vals,
        bool is_descending,
        int maxbits)
{
    utZoneScoped;
    exint nelem = size();
    utZoneValue(nelem);
    bool has_vals = !vals.isEmpty();
    CE_Context *context = CE_Context::getContext();

    // Number of compute units and local memory size.
    cl::Device dev = context->getDevice();
    size_t ncompute = dev.getInfo<CL_DEVICE_MAX_COMPUTE_UNITS>();
    size_t nlocal = dev.getInfo<CL_DEVICE_LOCAL_MEM_SIZE>();

    // Bits per radix.
    constexpr size_t radixbits = 8;
    // Radix value.
    constexpr size_t radix = 1 << radixbits;
    // How many total bits in the type.
    size_t totalbits = 8 * sizeof(T);
    // Limit the number of bits considered if provided.
    // We still have to do it in multiples of radixbits.
    if (maxbits > 0)
    {
        totalbits = SYSmin(SYSroundUpToMultipleOf(size_t(maxbits), radixbits),
                           totalbits);
    }
    size_t npasses =  totalbits / radixbits;
    // Max number of elements each workitem can process.
    constexpr size_t maxelemperworkitem = 256;
    // Workitems per group.
    size_t nworkitems = 16;
    // Usually 1K per workitem at radixbits = 8.
    size_t localsize = sizeof(uint32) * radix * nworkitems;
    // Give a little extra headroom in local memory.
    while (localsize > nlocal * 1.2)
    {
        nworkitems /= 2;
        localsize = sizeof(uint32) * radix * nworkitems;
    }
    // Elements processed per group.
    size_t maxelempergroup = maxelemperworkitem * nworkitems;
    // Just a heuristic, this is the maximum possible number of groups.
    // We seem to get decreasing performance above 128.
    size_t ngroups = SYSmin(size_t(128), ncompute * 4);
    // Clamp groups from 1 to the min number necessary to process all elements.
    size_t mingroups = SYSroundUpToMultipleOf(size_t(nelem), maxelempergroup) / maxelempergroup;
    // size_t SYSclamp is ambigious on mac so use int version.
    ngroups = SYSclamp((int)ngroups, 1, (int)mingroups);

    // TODO - write a single workgroup/single pass sort for when ngroups == 1.
    cl::NDRange globalRange {ngroups * nworkitems};
    cl::NDRange localRange {nworkitems};
    const char *opt = is_descending ? "-D SORT_DESCENDING" : nullptr;
    UT_WorkBuffer wb;
    if (has_vals)
    {        
        // Append a value flag and type.
        wb.append(opt);
        wb.append(" -D HAS_values -D CEARRAY_VALUE_");
        appendElemType<V>(wb);
        opt = wb.buffer();
    }
    cl::Kernel khist = loadKernel("radix_sort_histogram", opt);
    cl::KernelFunctor hist = khist.bind(context->getQueue(), globalRange, localRange);
    cl::Kernel kreorder = loadKernel("radix_sort_reorder", opt);
    cl::KernelFunctor reorder = kreorder.bind(context->getQueue(), globalRange, localRange);

    // Each group takes nworkitems Kb at radixbits = 8.
    Buffer<uint32> histogram(radix * nworkitems * ngroups);
    Buffer<uint32> histosums(histogram.size());
    auto localarg = cl::__local(localsize);

    CE_ArrayView<uint32> histogram_view = histogram.view();
    CE_ArrayView<uint32> histosums_view = histosums.view();

    // Set up secondary key and value (if needed) buffers.
    Buffer<T> dst_array(nelem);
    CE_ArrayView<T> dst_arrayview = dst_array.view();

    Buffer<V> dst_vals(has_vals ? nelem : 0);
    CE_ArrayView<V> dst_valsview = dst_vals.view();

    // Pointers to alternate src and dst arrays.
    CE_ArrayView<T> *src = this;
    CE_ArrayView<V> *v_src = &vals;
    CE_ArrayView<T> *dst = &dst_arrayview;
    CE_ArrayView<V> *v_dst = &dst_valsview;

    for(size_t pass=0; pass < npasses; pass++)
    {
        // Histogram.
        hist(src->buffer(), nelem, static_cast<int>(pass), histogram.buffer(), localarg);
        // Histogram prefix sum.
        histogram_view.prefixSum(histosums_view);
        // Sort, possibly with values.
        if (has_vals)
        {
            reorder(src->buffer(), nelem, static_cast<int>(pass), histosums.buffer(),
                    dst->buffer(), v_src->buffer(), v_dst->buffer(), localarg);
            UTswap(v_src, v_dst);
        }
        else
        {
            reorder(src->buffer(), nelem, static_cast<int>(pass),
                    histosums.buffer(), dst->buffer(), localarg);
        }
        // Swap buffers.
        UTswap(src, dst);
    }

    // Copy results if we had an odd number of passes
    if (npasses & 1)
    {
        context->copyBuffer(dst_array.buffer(), buffer(), nelem * sizeof(T));
        if (has_vals)
            context->copyBuffer(dst_vals.buffer(), vals.buffer(), nelem * sizeof(V));
    }
}

#endif