我正在解决一个问题,它涉及到快速排序10个数字(int32)。我的应用程序需要以尽可能快的速度对10个数字进行百万次排序。我正在对一个包含数十亿个元素的数据集进行采样,每次我都需要从中挑选10个数字(简化)并对它们进行排序(并从排序的10个元素列表中得出结论)。

目前我正在使用插入排序,但我想我可以为我的10个数字的特定问题实现一个非常快速的自定义排序算法,这将击败插入排序。

我该如何处理这个问题呢?


当前回答

这个问题并没有说这是某种基于web的应用程序。有一件事引起了我的注意:

我正在对一个包含数十亿个元素的数据集进行采样,每次我都需要从中挑选10个数字(简化)并对它们进行排序(并从排序的10个元素列表中得出结论)。

As a software and hardware engineer this absolutely screams FPGA to me. I don't know what kind of conclusions you need to draw from the sorted set of numbers or where the data comes from, but I know it would be almost trivial to process somewhere between one hundred million and a billion of these "sort-and-analyze" operations per second. I've done FPGA-assisted DNA sequencing work in the past. It is nearly impossible to beat the massive processing power of FPGAs when the problem is well suited for that type of a solution.

在某种程度上,唯一的限制因素是将数据铲入FPGA的速度有多快,以及取出数据的速度有多快。

As a point of reference, I designed a high performance real-time image processor that received 32 bit RGB image data at a rate of about 300 million pixels per second. The data streamed through FIR filters, matrix multipliers, lookup tables, spatial edge detection blocks and a number of other operations before coming out the other end. All of this on a relatively small Xilinx Virtex2 FPGA with internal clocking spanning from about 33 MHz to, if I remember correctly, 400 MHz. Oh, yes, it also had a DDR2 controller implementation and ran two banks of DDR2 memory.

当工作在数百MHz时,FPGA可以在每次时钟转换中输出10个32位数字。当数据填满处理管道时,操作开始时会有短暂的延迟。在此之后,您应该能够在每个时钟获得一个结果。如果可以通过复制排序和分析管道使处理并行化,则会更多。原则上,解决方案几乎是微不足道的。

关键在于:如果应用程序不是pc绑定的,并且数据流和处理与FPGA解决方案“兼容”(无论是独立的还是作为机器中的协处理器卡),那么无论使用哪种算法,用任何语言编写的软件都无法击败可达到的性能水平。

我只是快速搜索了一下,找到了一篇可能对你有用的论文。看起来要追溯到2012年。在今天(甚至在过去),您可以在性能上做得更好。下面就是:

fpga上的排序网络

其他回答

您可以完全展开插入排序。

为了简化这一点,可以在没有函数开销的情况下使用递归模板。因为int已经是模板,所以它也可以是模板形参。这也使得编码10以外的数组大小变得很简单。

注意,要排序int x[10],调用insert_sort<int, 9>::sort(x);因为类使用最后一项的索引。这可以被包装,但是需要阅读更多的代码。

template <class T, int NUM>
class insert_sort;

template <class T>
class insert_sort<T,0>
// Stop template recursion
// Sorting one item is a no operation 
{
public:
    static void place(T *x) {}
    static void sort(T * x) {}
};

template <class T, int NUM>
class insert_sort
// Use template recursion to do insertion sort.
// NUM is the index of the last item, e.g. for x[10] call <9>
{
public:
    static void place(T *x)
    {
        T t1=x[NUM-1];
        T t2=x[NUM];
        if (t1 > t2)
        {
            x[NUM-1]=t2;
            x[NUM]=t1;
            insert_sort<T,NUM-1>::place(x);
        }
    }
    static void sort(T * x)
    {
        insert_sort<T,NUM-1>::sort(x); // Sort everything before
        place(x);                    // Put this item in
    }
};

在我的测试中,这比排序网络的例子要快。

既然可以移动,为什么要交换?一条x86高速缓存线有足够的额外内存供您执行归并排序。

我可能会分别插入排序索引0-1、2-4、5-6、7-9,或者更好的是在插入时保持这些小组的排序,这样每次插入最多需要一到两次移位。

然后合并5、6和7-9 -> 10-14,合并0-1和2-4 -> 5-9,最后合并5-9和10-14 -> 0-9

使用排序网络,以4为一组进行比较,这样就可以在SIMD寄存器中进行比较。一对打包的最小/最大指令实现了打包比较器函数。对不起,我现在没有时间去找我记得看到过的关于这个的页面,但希望在SIMD或SSE排序网络上搜索会找到一些东西。

x86 SSE确实为4个32位整型的向量提供了打包的32位整型的min和max指令。AVX2 (Haswell及后续版本)具有相同的功能,但用于256b的8个整型向量。还有高效的洗牌指令。

如果你有很多独立的小排序,用向量并行地做4到8个排序是可能的。特别是,如果你随机选择元素(所以要排序的数据在内存中不会连续),你可以避免打乱,只需要按照你需要的顺序进行比较。10个寄存器保存来自4个(AVX2: 8) 10个整数列表的所有数据,仍然留下6个reg作为临时空间。

如果还需要对相关数据进行排序,则向量排序网络的效率较低。在这种情况下,最有效的方法似乎是使用wrapped -compare来获得更改元素的掩码,并使用该掩码来混合相关数据的(引用)向量。

以下是运行在10个CUDA线程上的CUDA内核(秩排序算法),在42毫秒内对1000个数组排序1000次,每次排序42纳秒或每次排序~70个周期(1.7 GHz):

inline
__device__ int findOrder(const int mask, const int data0)
{
    const int order1 = data0>__shfl_sync(mask,data0,0);
    const int order2 = data0>__shfl_sync(mask,data0,1);
    const int order3 = data0>__shfl_sync(mask,data0,2);
    const int order4 = data0>__shfl_sync(mask,data0,3);
    const int order5 = data0>__shfl_sync(mask,data0,4);
    const int order6 = data0>__shfl_sync(mask,data0,5);
    const int order7 = data0>__shfl_sync(mask,data0,6);
    const int order8 = data0>__shfl_sync(mask,data0,7);
    const int order9 = data0>__shfl_sync(mask,data0,8);
    const int order10 = data0>__shfl_sync(mask,data0,9);
    return order1 + order2 + order3 + order4 + order5 + order6 + order7 + order8 + order9 + order10;
}

// launch this with 10 CUDA threads in 1 block in 1 grid
// sorts 10 elements using only SIMT registers
__global__ void rankSort(int * __restrict__ buffer)
{    
    const int id  = threadIdx.x;

    // enable first 10 lanes of warp for shuffling 
    const int mask = __activemask();

    __shared__ int data[10000];

    for(int i=0;i<1000;i++)
    {
        data[id+i*10] = buffer[id+i*10];
    }
    __syncwarp();
    // benchmark 1000 times
    for(int k=0;k<1000;k++)
    {

        // sort 1000 arrays
        for(int j=0;j<1000;j+=5)
        {
            // sort 5 arrays at once to hide latency
            const int data1 = data[id+j*10];
            const int data2 = data[id+(j+1)*10];
            const int data3 = data[id+(j+2)*10];
            const int data4 = data[id+(j+3)*10];
            const int data5 = data[id+(j+4)*10];

            const int order1 = findOrder(mask,data1);        
            const int order2 = findOrder(mask,data2);
            const int order3 = findOrder(mask,data3);
            const int order4 = findOrder(mask,data4);
            const int order5 = findOrder(mask,data5);

            data[order1+j*10]=data1;         
            data[order2+(j+1)*10]=data2;           
            data[order3+(j+2)*10]=data3;
            data[order4+(j+3)*10]=data4;
            data[order5+(j+4)*10]=data5;
        }

    }
    __syncwarp();
    for(int i=0;i<1000;i++)
    {
        buffer[id+i*10] = data[id+i*10];
    }
}   

由于所有10个线程都在同一个SIMT单元上处理,它类似于运行在矢量寄存器上的AVX512优化版本,但除了更多的寄存器空间以隐藏更多的延迟之外。此外,SIMT单元是32宽的,因此它可以运行线性时间复杂度直到32个元素。

该算法假设元素是唯一的。对于重复的数据,需要一个额外的缩减步骤来将重复的顺序值解压缩为10个元素。首先,它计算每个元素的秩,然后直接将它们复制到数组中作为索引的秩。顺序值需要蛮力(O(N x N))次比较,为了减少延迟,使用warp-shuffles从(向量寄存器的)不同的warp-lanes收集数据。

出于类似于我在这里描述的原因,以下排序函数sort6_iterator()和sort10_iterator_local()应该能很好地执行,其中排序网络是从这里取的:

template<class IterType> 
inline void sort10_iterator(IterType it) 
{
#define SORT2(x,y) {if(data##x>data##y)std::swap(data##x,data##y);}
#define DD1(a)   auto data##a=*(data+a);
#define DD2(a,b) auto data##a=*(data+a), data##b=*(data+b);
#define CB1(a)   *(data+a)=data##a;
#define CB2(a,b) *(data+a)=data##a;*(data+b)=data##b;
  DD2(1,4) SORT2(1,4) DD2(7,8) SORT2(7,8) DD2(2,3) SORT2(2,3) DD2(5,6) SORT2(5,6) DD2(0,9) SORT2(0,9) 
  SORT2(2,5) SORT2(0,7) SORT2(8,9) SORT2(3,6) 
  SORT2(4,9) SORT2(0,1) 
  SORT2(0,2) CB1(0) SORT2(6,9) CB1(9) SORT2(3,5) SORT2(4,7) SORT2(1,8) 
  SORT2(3,4) SORT2(5,8) SORT2(6,7) SORT2(1,2) 
  SORT2(7,8) CB1(8) SORT2(1,3) CB1(1) SORT2(2,5) SORT2(4,6) 
  SORT2(2,3) CB1(2) SORT2(6,7) CB1(7) SORT2(4,5) 
  SORT2(3,4) CB2(3,4) SORT2(5,6) CB2(5,6) 
#undef CB1
#undef CB2
#undef DD1
#undef DD2
#undef SORT2
}

为了调用这个函数,我给它传递了一个std::vector迭代器。