【IT168专稿】Thrust是一个类似于STL的针对CUDA的C++模板库，能够使程序更简洁易读。Thrust提供与CUDA C完全兼容的接口，可以使我们高效地编写高性能并行程序。通过利用Thrust，程序员得以快速构建CUDA程序，并能够获得极高的稳定性和性能与精度，并行排序等例程的速度可提升5至100倍。
在之前文章中，我们给大家介绍了Thrust的快速入门的基础教程 (相关阅读：Thrust快速入门教程-基础)。今天给大家介绍的是Thrust的算法。
以下译文来自忆幽梦的博客，请参见：http://blog.cudachina.org/dreampursue/
Thrust提供了大量的常用并行算法。这些算法与STL的算法非常相似，于是我们使用了相同的名称(例如thrust::sort 与std::sort)。
所有的Thrust算法可以在主机端和设备端上使用。尤其是，当Thrust算法转入主机端迭代器时，将会调度主机端路径，同样，当使用设备端迭代器时将使用设备端实现。
thrust::copy是一个例外，他可以任意的拷贝主机端和设备端的数据。但是所有的迭代器参数必须符合Thrust算法的要求，要么都在主机端，要么都在设备端。当不能满足要求的时候，编译器会报错。
Transformations
Transformations算法作用是用来将目标容器赋上特定值(例如零)或者特定数列。之前的例子我们已经使用过thrust::fill，可以向所有元素赋特定值。此外transformations算法还包括thrust::sequence、thrust::replace、thrust::transform。完整的列表请参考文档。
下面的代码演示了几个transformation算法的用法。注意类似于C++中拥有的thrust::negate和thrust::modulus，Thrust在thrust/functional.h中也提供了，此外还有plus与multiplies等。
# include <thrust / transform .h>
# include <thrust / sequence .h>
# include <thrust / copy .h>
# include <thrust / fill .h>
# include <thrust / replace .h>
# include <thrust / functional .h>
# include <iostream >
int main ( void )
{
// allocate three device_vectors with 10 elements
thrust :: device_vector <int > X (10) ;
thrust :: device_vector <int > Y (10) ;
thrust :: device_vector <int > Z (10) ;
// initialize X to 0,1,2,3, ....
thrust :: sequence (X. begin () , X. end ());
// compute Y = -X
thrust :: transform (X. begin () , X.end () , Y. begin () , thrust :: negate <int >() );
// fill Z with twos
thrust :: fill (Z. begin () , Z. end () , 2);
// compute Y = X mod 2
thrust :: transform (X. begin () , X.end () , Z. begin () , Y. begin () , thrust :: modulus <int -
>() );
// replace all the ones in Y with tens
thrust :: replace (Y. begin () , Y. end () , 1, 10) ;
// print Y
thrust :: copy (Y. begin () , Y. end () , std :: ostream_iterator <int >( std :: cout , "\n"));
return 0;
}
thrust/fuctuional.h中的函数提供了大部分内置代数和比较运算，但是我们想提供更多出色的功能。比如，运算y < - a * x + y，x、y为向量，a为常数标量。这其实就是我们所熟知的由BLAS提供的SAXPY运算。
如果我们在thrust中实现SAXPY我们有几个选择。一个是，我们需要使用两个transformations(一个加和一个乘法)还有一个临时数则用于存储a乘后的值。另一更佳选择是使用一个单独的由用户自己定义函数的transformation，这才是我们真正先要的。我下面用源代码解释说明这两种方法。
{
const float a;
saxpy_functor ( float _a) : a(_a) {}
__host__ __device__
float operator ()( const float & x, const float & y) const {
return a * x + y;
}
};
void saxpy_fast ( float A, thrust :: device_vector <float >& X, thrust :: device_vector < -
float >& Y)
{
// Y <- A * X + Y
thrust :: transform (X. begin () , X.end () , Y. begin () , Y. begin () , saxpy_functor (A));
}
void saxpy_slow ( float A, thrust :: device_vector <float >& X, thrust :: device_vector < -
float >& Y)
{
thrust :: device_vector <float > temp (X. size ());
// temp <- A
thrust :: fill ( temp . begin () , temp . end () , A);
// temp <- A * X
thrust :: transform (X. begin () , X.end () , temp . begin () , temp . begin () , thrust :: -
multiplies <float >() );
// Y <- A * X + Y
thrust :: transform ( temp . begin () , temp . end () , Y. begin () , Y. begin () , thrust :: plus < -
float >() );
}
Saxpy_fast和saxpy_slow都是有效的SAXPY实现，尽管saxpy_fast会比saxpy_slow更快。忽略临时向量分配与代数运算的花费，其开销如下：
fast_saxpy：2N次读取和N次写入
slow_saxpy：4N次读取和3N写入
因为SAXPY受到内存约束(它的性能受限于内存的带宽，而不是浮点性能)更大量的读写操作使得saxpy_slow开销更加昂贵。而saxpy_fast执行速度与优化的BLAS实现中的SAXPY一样快。在类似SAXPY内存约束算法通常值得使用kernel融合(合并多个计算于单独的kernel)的方法以最小化内存的读写交换。
Thrust::transform只支持一个或者两个输入参数的transformations(例如f(x) -> y 和 f(x; y) -> z)。当transformation使用多于两个输入参数的时候需要使用其他方法了。例子arbitrary_transformation展示了使用thrust::zip_interator和thrust::for_each的解决方案。
Reductions
Reduction算法使用二元操作将输入序列规约为一个单值。例如，需要获得一数列的和，可以通过加运算规约此数组得到。相似的，数列的最大值，可以通过由两个输入值返回一个最大值的运算子规约得到。数列的求和的规约操作可以由thrust::reduce如下实现：
开始的两个参数定义了需要规约的数组，第三和第四个参数分别提供了初始值和相关的规约操作。实际上，通常使用的时候我们选择默认情况下没有初始值和不特别指出规约方法。所以下面三行代码是等同的：
int sum = thrust :: reduce (D. begin () , D. end () , ( int ) 0);
int sum = thrust :: reduce (D. begin () , D. end ())
虽然thrust::reduce能够有效的满足大部分的规约操作，但是，Thrust库依然提供了另外的一些函数以便使用(类似于STL)。例如，thrust::count能够返回给定序列的特定值的数量。
# include <thrust / device_vector .h>
...
// put three 1s in a device_vector
thrust :: device_vector <int > vec (5 ,0);
vec [1] = 1;
vec [3] = 1;
vec [4] = 1;
// count the 1s
int result = thrust :: count ( vec . begin () , vec .end () , 1);
// result is three
另一些规约操作，包括thrust::count_if、thrust::min_element、thrust::max_element、thrust::is_sorted、thrust::inner_product等，详细请参考documentation。
Transformations篇章中的SAXPY例子使用transformation内核展示了融合内核如何来减少内存交换。我们也可以使用thrust::transform_reduce实现融合内核来规约。下面的例子用来计算向量的模：
# include <thrust / functional .h>
# include <thrust / device_vector .h>
# include <thrust / host_vector .h>
# include <cmath >
// square <T> computes the square of a number f(x) -> x*x
template <typename T>
struct square
{
__host__ __device__
T operator ()( const T& x) const {
return x * x;
}
};
int main ( void )
{
// initialize host array
float x [4] = {1.0 , 2.0 , 3.0 , 4.0};
// transfer to device
thrust :: device_vector <float > d_x (x, x + 4);
// setup arguments
square <float > unary_op ;
thrust :: plus <float > binary_op ;
float init = 0;
// compute norm
float norm = std :: sqrt ( thrust :: transform_reduce ( d_x . begin () , d_x . end () , -
unary_op , init , binary_op ) );
std :: cout << norm << std :: endl ;
return 0;
}
这里有一个叫平方的一元操作，将输入序列的每个元素平方。然后，平方的和由标准的加法规约操作得到。类似较慢版本的SAXPY transformation，我们可以这样实现：首先将原数列同伙乘法转换成平方存储在一个临时数组，然后使用加法规约。但是显然这样做会带来不必要的浪费，速度降低。通过在规约内核中融合平方操作，我们就可以获得与自己编写内核一样的高性能。
Prefix-Sums
并行的前追求和，也叫scan操作，与压实流、基数排序等都是并行算法的重要模块。下面的源码将举例说明使用默认加法的inclusive scan：
int data [6] = {1, 0, 2, 2, 1, 3};
thrust :: inclusive_scan (data , data + 6, data ); // in - place scan
// data is now {1, 1, 3, 5, 6, 9}
Inclusive scan的每个输出元素为输入数列的相应部分和。例如，data[2] = data[0] + data[1] + data[2]。Exclusive scan类似，但是右移一个位置：
int data [6] = {1, 0, 2, 2, 1, 3};
thrust :: exclusive_scan (data , data + 6, data ); // in - place scan
// data is now {0, 1, 1, 3, 5, 6}
现在为data[2] = data[0] + data[1]。由例子可见，inclusive_sacn与exclusive_scan允许原址操作。Thrust也提供了函数transform_inclusive_scan与transform_exclusive_scan可以实现在scan操作前对输入数列进行一元操作。完整的scan变体说明请参见documentation。
以上译文来自忆幽梦的博客，请参见：http://blog.cudachina.org/dreampursue/
Reordering
Thrust通过下面的算法为partitioning和stream compaction提供支持：
1.copy_if：复制通过谓词测试的元素。
2. partition(分区)：根据谓词重新排列元素(真值在假值前)。
3. remove and remove_if：删除谓词测试失败的元素。
4. unique：在一系列参考文档中删除连续的副本，这些参考文档是关于函数重排和用法实例的完整列表。
Sorting
Thrust提供一系列函数根据既定的准则来分类数据和整理数据。thrust::sort和thrust::stable_sort函数是在STL(Standard Template Library 标准模板库)中的sort和stable_sort非常相似。
...
const int N = 6;
int A[N] = {1, 4, 2, 8, 5, 7};
thrust :: sort (A, A + N);
// A is now {1, 2, 4, 5, 7, 8}
另外，Thrust提供thrust::sort_by_key和thrust::stable_sort_by_key，这两个函数将 key-value 分为两类，分别存储在不同的位置。
...
const int N = 6;
int keys [N] = { 1, 4, 2, 8, 5, 7};
char values [N] = {'a', 'b', 'c', 'd', 'e', 'f'};
thrust :: sort_by_key (keys , keys + N, values );
// keys is now { 1, 2, 4, 5, 7, 8}
// values is now {'a ', 'c ', 'b ', 'e ', 'f ', 'd '}
跟STL类似，sorting函数同样接受用户自定义对比操作。
# include <thrust / functional .h>
...
const int N = 6;
int A[N] = {1, 4, 2, 8, 5, 7};
thrust :: stable_sort (A, A + N, thrust :: greater <int >());
// A is now {8, 7, 5, 4, 2, 1}
更多内容请点击：
CUDA专区：http://cuda.it168.com/
CUDA论坛：http://cudabbs.it168.com/
