使用 CUDA 和 OpenCL 进行排列
4.57/5 (6投票s)
在 GPU 上查找字典序排列
引言
CUDA 和 OpenCL 使任何应用程序都能访问图形处理单元进行非图形计算，并能利用 GPU 进行大规模通用计算。在 CUDA 和 OpenCL 出现之前，GPGPU 是通过图形 API 着色器（HLSL 和 GLSL）完成的，这些着色器在处理通用计算任务时显得笨拙而不自然，因为程序员必须“弯曲”着色器来完成通用编程。CUDA 可用于 Nvidia 图形卡和 x86 CPU。OpenCL 可用于 Nvidia 和 AMD 图形卡、x86 CPU 和 ARM CPU。而 DirectCompute 仅在 MS Vista 和 Windows 7 平台以及 DirectX 10/11 级图形卡上可用。要阅读本文，读者应相当熟悉 CUDA 和 OpenCL 的工作原理、其内存模型以及内核（kernel）、主机（host）和设备（device）等技术术语。
算法
本文使用的算法与我之前的 《C++ 中的排列，第 2 部分》 文章中的算法基本相同。该算法利用 阶乘 分解来查找字典序排列，并且是从零开始计数的。正如读者从我上一篇文章的评分中可以看到，我对算法的解释不太擅长，读者可能希望在其他地方阅读解释。为了解释这个算法，我们将使用两组元素，第一组称为未选定集（unpicked set），最初存储所有元素；第二组称为已选定集（picked set），最初为空，但最终将包含最终的排列。元素逐个地从未选定集移动到已选定集。让我们举一个例子。假设我们有一个未选定集 {0,1,2}。通过计算 3 的阶乘，我们知道 3 个元素总共有 6 种不同的排列。注意：阶乘运算通常用 ! 符号表示，例如 2 的阶乘简单地写成 2!。1 的阶乘定义为 1。
0,1,2 0,2,1 1,0,2 1,2,0 2,0,1 2,1,0
从（零基）字典序排列列表中，我们知道第一个排列（第 0 个）是 0,1,2，最后一个排列（第 5 个）是 2,1,0。让我们用阶乘分解来找出第 4 个排列（我们从上面的列表中得知，它是 2,0,1）。
2! * 0 = 0 2! * 1 = 2 2! * 2 = 4
该算法要求我们找到大于或等于我们的第 n 个索引（这里是 4）的最高乘积，所以我们知道第一个元素是 2，因为 4 等于第 3 个乘积，而未选定集 {0,1,2} 中的第 3 个元素是 2。让我们继续找出第二个元素。在此之前，我们必须从我们的第 n 个索引中减去找到的乘积，即 4 - 4 = 0。{2,?,?} 是我们当前的已选定集，未选定集中剩下 {0,1}。
1! * 0 = 0 1! * 1 = 1
同样，我们找到大于或等于我们减去的第 n 个索引（这里是 0）的最高数字，所以我们知道第二个元素是 0，因为 0 等于第 1 个乘积，而未选定集 {0,1} 中的第 1 个元素是 0。{2,0,?} 是我们当前部分找到的排列，未选定集中剩下 {1}。我们无需进一步计算，因为未选定集中只剩下 1。所以第 4 个排列是 {2,0,1}。让我们再试一个例子！这次，我们将找出 3 个元素 {0,1,2} 的第 3 个排列。
2! * 0 = 0 2! * 1 = 2 2! * 2 = 4
正如读者现在应该已经知道的，我们需要找到大于或等于我们的第 n 个索引（这里是 3）的最高乘积。我们知道第一个元素是 1，因为 3 小于第 3 个乘积但大于第 2 个乘积，而未选定集 {0,1,2} 中的第 2 个元素是 1。让我们继续找出已选定集的第二个元素。首先，我们必须从我们的第 n 个索引中减去找到的乘积，即 3 - 2 = 1。注意：{1,?,?} 是我们当前已选定集，未选定集中剩下 {0,2}。
1! * 0 = 0 1! * 1 = 1
同样，我们找到大于或等于我们减去的第 n 个索引（这里是 1）的最高数字，所以我们知道第二个元素是 2，因为 2 是未选定集 {0,2} 中的第 2 个元素。{1,2,?} 是我们当前部分找到的排列，未选定集中剩下 {0}。因为未选定集中只剩下 0，所以第 3 个排列是 {1,2,0}。
CUDA 源代码
我将展示 CUDA 内核的源代码。内核函数名为 Permute。arrDest 参数是一个未初始化的数组，将在 CUDA 操作完成后填充结果。offset 是要开始计算的第 n 个排列，Max 是要计算的最大排列。作为说明 offset 和 Max 参数用法的简单示例，假设您有 1000 个排列要查找，但您的 GPU 只有足够的内存一次计算 500 个排列，因此您必须将 CUDA 操作分成两部分。第一个操作计算第 0 个到第 499 个排列（offset=0，Max=500），第二个操作计算第 500 个到第 999 个排列（offset=500，Max=500）。算法中使用的第 n 个排列是从零开始计数的。常量数组 arr[i][j] 存储了阶乘(i) * j 的预计算值，其中 i 和 j 都是从零开始计数的，并且假定阶乘(0) 为零而不是通常的 1，以使代码正常工作。CUDA 中的常量内存访问比全局内存快，但它们通常限于 64KB 并且不能写入。当我们从 n 个元素的集合中查找第 k 个元素时，我们使用 (n-k) 的阶乘。arr 的每个元素都是一个 64 位整数，最多可以存储到阶乘 20。next_permutation 是从 Visual C++ 10 STL 复制和修改而来的。NEXT_PERM_LOOP 指定了您希望使用 next_permutation 的次数。
#define NUMBER_OF_ELEMENTS 5
#define BLOCK_DIM 1024
#define OFFSET 0
// When MAX_PERM = 0, means find all permutations
#define MAX_PERM 0
#define NEXT_PERM_LOOP 1
__constant__ long long arr[20][20] = { /*Not shown here to save space*/ };
// function to swap character
// a - the character to swap with b
// b - the character to swap with a
__device__ void swap(
char* a,
char* b)
{
char tmp = *a;
*a = *b;
*b = tmp;
}
// function to reverse the array (sub array in array)
// first - 1st character in the array (sub-array in array)
// last - 1 character past the last character
__device__ void reverse(
char* first,
char* last)
{
for (; first != last && first != --last; ++first)
swap(first, last);
}
// function to find the next permutation (sub array in array)
// first - 1st character in the array (sub-array in array)
// last - 1 character past the last character
__device__ void next_permutation(
char* first,
char* last)
{
char* next = last;
--next;
if(first == last || first == next)
return;
while(true)
{
char* next1 = next;
--next;
if(*next < *next1)
{
char* mid = last;
--mid;
for(; !(*next < *mid); --mid)
;
swap(next, mid);
reverse(next1, last);
return;
}
if(next == first)
{
reverse(first, last);
return;
}
}
}
__global__ void PermuteHybrid(char* arrDest, long long* offset, long long* Max)
{
long long index = threadIdx.x + blockIdx.x * blockDim.x;
if(index >= (*Max/(NEXT_PERM_LOOP+1)))
return;
index *= NEXT_PERM_LOOP+1;
long long tmpindex = index;
index += *offset;
char arrSrc[NUMBER_OF_ELEMENTS];
char arrTaken[NUMBER_OF_ELEMENTS];
for(char i=0; i<NUMBER_OF_ELEMENTS; ++i)
{
arrSrc[i] = i;
arrTaken[i] = 0;
}
char size = NUMBER_OF_ELEMENTS;
for(char i=NUMBER_OF_ELEMENTS-1; i>=0; --i)
{
for(char j=i; j>=0; --j)
{
if(index >= arr[i][j])
{
char foundcnt = 0;
index = index - arr[i][j];
for(char k=0;k<NUMBER_OF_ELEMENTS; ++k)
{
if(arrTaken[k]==0) // not taken
{
if(foundcnt==j)
{
arrTaken[k] = 1; // set to taken
arrDest[ (tmpindex*NUMBER_OF_ELEMENTS) + (NUMBER_OF_ELEMENTS-size) ] = arrSrc[k];
break;
}
foundcnt++;
}
}
break;
}
}
--size;
}
long long idx = tmpindex*NUMBER_OF_ELEMENTS;
for(char a=1; a<NEXT_PERM_LOOP+1; ++a)
{
long long idx2 = a*NUMBER_OF_ELEMENTS;
for(char i=0; i<NUMBER_OF_ELEMENTS; ++i)
{
arrDest[ idx + idx2 + i ] =
arrDest[ idx + ((a-1)*NUMBER_OF_ELEMENTS) + i ];
}
next_permutation(arrDest + idx + idx2,
arrDest+idx + idx2 + NUMBER_OF_ELEMENTS);
}
}
要控制操作参数，我们需要相应地设置宏。例如，下面的宏从第 0 个到第 29 个排列查找 5 个元素的排列。如果内核无法运行，读者可能需要尝试修改 BLOCK_DIM 的值，使其为 256 的倍数。
#define NUMBER_OF_ELEMENTS 5
#define BLOCK_DIM 1024
#define OFFSET 0
// When MAX_PERM = 0, means find all permutations
#define MAX_PERM 30
#define NEXT_PERM_LOOP 1
我包含了一个 check 函数来将生成的排列与 next_permutation 进行比较，还有一个 display 函数来显示结果。对于大量排列，读者可能希望注释掉 display 函数。
void check(char* arrDest, long long Max)
{
printf("\nChecking...\n");
char check[NUMBER_OF_ELEMENTS];
for(int i=0; i<NUMBER_OF_ELEMENTS ;++i)
{
check[i] = i;
}
if(OFFSET!=0)
{
for(int i=0; i<OFFSET; ++i)
{
std::next_permutation(check, check+NUMBER_OF_ELEMENTS);
}
}
for(int i=0; i<Max ;++i)
{
for(int j=0;j<NUMBER_OF_ELEMENTS;++j)
{
if(arrDest[i*NUMBER_OF_ELEMENTS+j] != check[j])
{
fprintf(stderr, "Diff check failed at %d!", i);
return;
}
}
std::next_permutation(check, check+NUMBER_OF_ELEMENTS);
}
}
void display(char* arrDest, long long Max)
{
for(int i=0; i<Max ;++i)
{
for(int j=0;j<NUMBER_OF_ELEMENTS;++j)
printf("%d ", arrDest[i*NUMBER_OF_ELEMENTS+j]);
printf("\n");
}
}
OpenCL 源代码
下面展示了 OpenCL 内核函数。OpenCL 代码与 CUDA 代码类似，除了主机代码和内核代码写在两个单独的源文件中。这导致需要将 NUMBER_OF_ELEMENTS 和 NEXT_PERM_LOOP 宏设置两次（一次在主机文件中，一次在内核源文件中）。arr 是一个一维数组，与 CUDA 中的二维数组不同，因为 OpenCL 中不可能定义二维数组。
#pragma OPENCL EXTENSION cl_khr_byte_addressable_store: enable
#define NUMBER_OF_ELEMENTS 5
#define NEXT_PERM_LOOP 1
__constant long arr[400] = { /*Not shown here to save space*/ };
// function to swap character
// a - the character to swap with b
// b - the character to swap with a
void swap(
__global char* a,
__global char* b)
{
char tmp = *a;
*a = *b;
*b = tmp;
}
// function to reverse the array (sub array in array)
// first - 1st character in the array (sub-array in array)
// last - 1 character past the last character
void reverse(
__global char* first,
__global char* last)
{
for (; first != last && first != --last; ++first)
swap(first, last);
}
// function to find the next permutation (sub array in array)
// first - 1st character in the array (sub-array in array)
// last - 1 character past the last character
void next_permutation(
__global char* first,
__global char* last)
{
__global char* next = last;
--next;
if(first == last || first == next)
return;
while(true)
{
__global char* next1 = next;
--next;
if(*next < *next1)
{
__global char* mid = last;
--mid;
for(; !(*next < *mid); --mid)
;
swap(next, mid);
reverse(next1, last);
return;
}
if(next == first)
{
reverse(first, last);
return;
}
}
}
__kernel void PermuteHybrid(__global char* arrDest, __global long* offset, __global long* Max)
{
long index = get_global_id(0);
if(index>=(*Max/(NEXT_PERM_LOOP+1)))
return;
index *= NEXT_PERM_LOOP+1;
long tmpindex = index;
index += *offset;
char arrSrc[NUMBER_OF_ELEMENTS];
char arrTaken[NUMBER_OF_ELEMENTS];
for(int i=0; i<NUMBER_OF_ELEMENTS; ++i)
{
arrSrc[i] = i;
arrTaken[i] = 0;
}
int size = NUMBER_OF_ELEMENTS;
for(char i=NUMBER_OF_ELEMENTS-1; i>=0; --i)
{
for(char j=i; j>=0; --j)
{
if(index >= arr[i*20+j])
{
char foundcnt = 0;
index = index - arr[i*20+j];
for(char k=0;k<NUMBER_OF_ELEMENTS; ++k)
{
if(arrTaken[k]==0) // not taken
{
if(foundcnt==j)
{
arrTaken[k] = 1; // set to taken
arrDest[ (tmpindex*NUMBER_OF_ELEMENTS) + (NUMBER_OF_ELEMENTS-size) ] = arrSrc[k];
break;
}
foundcnt++;
}
}
break;
}
}
--size;
}
long idx = tmpindex*NUMBER_OF_ELEMENTS;
for(char a=1; a<NEXT_PERM_LOOP+1; ++a)
{
long idx2 = a*NUMBER_OF_ELEMENTS;
for(char i=0; i<NUMBER_OF_ELEMENTS; ++i)
{
arrDest[ idx + idx2 + i ] =
arrDest[ idx + ((a-1)*NUMBER_OF_ELEMENTS) + i ];
}
next_permutation(arrDest + idx + idx2,
arrDest+idx + idx2 + NUMBER_OF_ELEMENTS);
}
}
OpenCL 版本也有类似的 check 和 display 函数。应根据需要注释掉它们的调用。
void check(char* arrDest, long long Max)
{
std::cout << std::endl;
std::cout << "Checking..." << std::endl;
char check[NUMBER_OF_ELEMENTS];
for(int i=0; i<NUMBER_OF_ELEMENTS ;++i)
{
check[i] = i;
}
if(OFFSET!=0)
{
for(int i=0; i<OFFSET; ++i)
{
std::next_permutation(check, check+NUMBER_OF_ELEMENTS);
}
}
for(int i=0; i<Max ;++i)
{
for(int j=0;j<NUMBER_OF_ELEMENTS;++j)
{
if(arrDest[i*NUMBER_OF_ELEMENTS+j] != check[j])
{
printf("Diff check failed at %d!", i);
return;
}
}
std::next_permutation(check, check+NUMBER_OF_ELEMENTS);
}
}
void display(char* arrDest, long long Max)
{
for(int i=0; i<Max ;++i)
{
for(int j=0;j<NUMBER_OF_ELEMENTS;++j)
std::cout << (int)(arrDest[i*NUMBER_OF_ELEMENTS+j]);
std::cout << std::endl;
}
}
基准测试
我对 CUDA 和 CPU 应用程序 进行了基准测试。基准测试中使用的 CPU 和 GPU 分别是 Intel i7 870（8 核，2.93Ghz）和 Nvidia Geforce 460。CPU 应用程序充分利用了 8 个核心来查找排列。CPU 应用程序使用阶乘分解将第 n 个排列分配给不同的 CPU 核心，在每个工作线程中，使用 STL next_permutation 来查找从第 n 个排列到每个连续排列。计算 11 个元素排列的结果如下表所示。找到的 11 个元素的总排列数为 39,916,800。存储结果所需的数组大小为 39,916,800 x 11 = 439,084,800。这是我的 1GB 内存 GPU 能够存储的最大排列数。
CPU : 550ms Version 1 with pure factorial decomposition (Average timing) CUDA : 550ms OpenCL : 581ms Version 2 with 1 next_permutation per factorial decomposition (Average timing) CUDA : 317ms OpenCL : 373ms Version 3 with 9 next_permutation per factorial decomposition (Average timing) CUDA : 681ms OpenCL : 456ms
8 核 CPU 和 GPU 的基准测试结果大致相同。这是否意味着 GPU 排列毫无用处？GPU 排列绝非毫无用处，因为它允许将排列和自定义内核完全在 GPU 上完成，而不是采用 CPU/GPU 混合方法。阶乘分解的复杂度为 O(n^2)，而 STL next_permutation 的复杂度为 O(n)，但我们不能在 GPU 上纯粹使用 next_permutation，因为 next_permutation 依赖于当前排列来生成下一个。然而，阶乘分解可以与 GPU 上的 next_permutation 结合使用以提高性能（类似于我为 CPU 基准测试所做的）。next_permutation 的使用次数越多，对性能的边际效益似乎就越小，因为必须有大量的 workitem 才能使 GPU 工作单元保持忙碌；300ms 是 2、3 或 4 次 next_permutation 的最低时间，之后性能会下降。版本 1（纯阶乘分解）很容易修改，以便从自定义内核调用。想象一下您使用的是版本 3，它可以在一次内核调用中为您提供不同数量的 next_permutation 的灵活性，如果您的自定义内核只需要 1 个排列输入，那么程序员将如何处理 next_permutation 生成的其余排列？
大规模排列基准测试
我决定修改代码来基准测试查找 12 个元素的排列。总排列数为 479,001,600。存储所有排列所需的总字节数为 479,001,600 * 12 = 5,748,019,200（约 6GB），这远远超出了我的 1GB 显卡能够存储的范围。我将计算分成 12 个批次，以适应我的显卡内存。您选择拆分批次的数字必须是总排列数的因子，这意味着总排列数可以被该数字整除而没有余数。第一批使用阶乘分解计算，其余批次使用纯 next_permutation。结果如下。
Timing: 649ms Checking... next_permutation Timing: 18ms Checking... next_permutation Timing: 29ms Checking... next_permutation Timing: 18ms Checking... next_permutation Timing: 29ms Checking... next_permutation Timing: 18ms Checking... next_permutation Timing: 32ms Checking... next_permutation Timing: 18ms Checking... next_permutation Timing: 28ms Checking... next_permutation Timing: 18ms Checking... next_permutation Timing: 29ms Checking... next_permutation Timing: 19ms Checking... Total timing: 905ms Executed program successfully.
第一次使用阶乘分解计时为 649 毫秒。总计时为 905 毫秒。CPU 版本需要 6.5 秒。GPU 性能比 CPU 版本提高了 7 倍！读者可以在 PermCudaBenchmark 文件夹下找到此基准测试源代码。读者可能希望注释掉耗时很长的检查函数，该函数需要检查所有生成的排列。
限制
本文提出的算法最多只能查找 20 个元素的排列。这是因为 64 位整数最多只能存储到阶乘(20)。注意：拥有 1GB 显卡可以存储 11 个元素的排列（约 0.5GB）。对于 12 个或更多元素的排列，程序员必须将操作分解成批次（最好是第一批使用阶乘分解，其余批次使用 next_permutation）。如 wiki 文档 所述，oclHashcat-plus 的排列限制为 11。对于不熟悉的人来说，oclHashcat-plus 是在 GPU 上运行的高级密码恢复软件。此外，该算法不处理重复元素。谨防：OpenCL 代码仅在 Nvidia 显卡上进行了测试，AMD 用户必须启用检查函数（默认启用）以确保排列的正确性。如中所述，OpenCL 代码未在任何 x86 OpenCL 实现上进行测试。
关注点
我对排列的热爱是如何产生的？故事始于 1998 年，当时当地广播电台 FM 93.3 举办了一个有奖游戏。广播电台会给出 6 位数字，听众必须将这 6 位数字排列成 2 个三位数，这两个数字加起来必须等于 933。现金奖非常诱人，而且每次听众未能给出正确答案时都会累积。我采用了蛮力方法来解决这个问题。花了几天时间才在 Visual Basic 6 中实现排列。有趣的是，我从未成功拨通广播电台的热线电话给出我的答案。我无助地听着，看着 1000 多美元的现金奖从我身边溜走。
结论
在本文中，我们研究了如何在 GPU 上实现字典序排列的生成，使用了 Nvidia 的 CUDA。基准测试似乎表明 CPU 和 GPU 具有相似的性能。我们通过阶乘分解和 next_permutation 的混合方法将性能提高了 40%。GPU 的优势在于可以将自定义操作主要在 GPU 上运行，而不是采用 CPU/GPU 混合方法（如果您预先在 CPU 上计算了排列）。要构建和运行 CUDA 源代码，读者必须下载最新的 CUDA SDK。源代码托管在 Codeplex。非常欢迎对如何进一步优化 CUDA 和 OpenCL 代码发表评论。
参考文献
- 《并行处理器编程：实践方法》（Programming Massively Parallel Processors: A Hands-on Approach），作者：David B. Kirk, Wen-mei W. Hwu
- 《CUDA 入门》（CUDA by Example），作者：Jason Sanders 和 Edward Kandrot
- 《OpenCL 编程指南》（OpenCL Programming Guide），作者：Aaftab Munshi, Benedict Gaster, Timothy G. Mattson, James Fung, Dan Ginsburg
- 《CUDA 应用程序设计与开发》（CUDA Application Design and Development），作者：Rob Farber
历史
- 2012-05-15：添加了 限制 部分。
- 2012-05-13：包含大规模排列基准测试，性能比 CPU 提高了 7 倍。
- 2012-05-09：使用 STL
next_permutation进行优化。 - 2012-05-09：添加了 OpenCL 部分和源代码。
- 2012-05-08：初次发布。