Trust简介

Thrust 是一个实现了众多基本并行算法的 C++ 模板库,类似于 C++ 的标准模板库(standard template library, STL)。该库自动包含在 CUDA 工具箱中。这是一个模板库,仅仅由一些头文件组成。在使用该库的某个功能时,包含需要的头文件即可。该库中的所有类型与函数都在命名空间thrust中定义,所以都以thrust::开头。用命名空间的目的是避免名称冲突。例如，Thrust中的thrust::sort和STL 中的 std::sort 就不会发生名称冲突。

数据结构

Thrust 中的数据结构主要是矢量容器(vector container)，类似于 STL中的std::vector。在 Thrust 中，有两种矢量:

(1)一种是存储于主机的矢量 thrust::host_vector<typename>。

(2)一种是存储于设备的矢量 thrust::device_vector<typename>。这里的 typename 可以是任何数据类型。例如，下面的语句定义了一个设备矢量x,元素类型为双精度浮点数(全部初始化为0)，长度为10:

thrust::device_vector<double>x(10，0);

要使用这两种矢量，需要分别包含如下头文件:

#incldue <thrust/host vector.h>

#incldue <thrust/device vector.h>

算法

Thrust 提供了5类常用算法，包括

(1)变换(transformation)。

(2)归约(reduction)。

(3)前缀和(prefxsum)。

(4)排序(sorting)与搜索(searching)。

(5)选择性复制、替换、移除、分区等重排(reordering)操作。

除了 thrust::copy，Thrust 算法的参数必须都来自于主机矢量或都来自于设备矢量。否则，编译器会报错。

---

实例分析

在了解 Thrust 库更多的细节之前，我们先分析Code1所示的程序，这个程序展示了Thrust库的一些显著特点。

Code1

#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>

#include <cuda_runtime.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <cstdlib>


int main()
{
    thrust::host_vector<int> h_vec(1 << 24);

    thrust::device_vector<int> d_vec = h_vec;

    thrust::generate(h_vec.begin(), h_vec.end(), rand);

    thrust::sort(d_vec.begin(), d_vec.end());

    thrust::copy(d_vec.begin(), d_vec.end(), h_vec.begin());

    return 0;
}

Code1分配了两个向量容器:host_vector与 device_vector。host_vector位于主机端，device_vector位于GPU设备端。Thrust 的向量容器与C++ STL中的向量容器类似，host_vector与 device_vector 是通用的容器(即可以存储任何数据类型)，可以动态调整大小。如Code1所示，容器可以自动分配和释放内存空间并且简化主机端和设备端之间的数据交换。

程序在向量容器上执行时，使用了generate、sort和copy算法。采用了STL中的迭代器进行遍历。在这个例子中，迭代器h_vec.beginO和h_vec.end()分别指向容器的第一个元素和最后一个元素的后一个位置（与STL一致左闭右开）。通过计算h_vec.end() – h_vec.beginO，我们可以得到容器的大小。

注意，在执行排序算法的时候，Thrust 会建议启动一个或多个CUDA kernel，但编程人员并不需要进行相关配置，因为Thrust的接口已经将这些细节抽象化了。对于性能敏感变量(比如 Thrust 库的网格和块大小)的选择，内存管理的细节，甚至排序算法的选择都留给具体实现的人自行决定。

迭代器和内存空间

虽然向量迭代器类似于数组的指针，但它们还包含了一些额外的信息。注意，我们不需要指定在 device_vector 元素上操作的sort算法，也不用暗示复制操作是从设备内存端到主机内存端。在Thrust库中，每个范围的内存空间可以通过迭代器参数自动推断，并调度合适的算法进行执行。

另外，关于内存空间，Thrust 的迭代器对大量信息进行隐式编码，这些信息可以用来指导进程调度。比如，Code1中sort的例子，它对基本的整型数据类型进行比较操作。在这个例子中，Thrust库中采用高度优化的基数排序(radix sort)算法，要比基于数据之间比较的排序算法(例如归并排序算法速度快很多。需要注意的是，这个调度过程并不会造成性能或存储开销:迭代器对元数据编码只存在于编译阶段并且它的调度策略已经确定。实际上，Thrust的静态调度策略可以利用迭代器类型的任何信息。

互操作性

Thrust库完全由CUDA C/C++实现，并且保持了与CUDA 生态系统其余部分的互操作性。互操作性是一个重要特性，因为没有一个单一的语言或库能够很好地解决所有问题。例如，尽管Thrust 算法在内部使用了像共享存储器的CUDA特性,但是并没有为用户提供机制通过 Thrust库直接使用共享存储器。因此，有时候应用程序需要直接访问CUDAC,实现一些特定的算法。Thrust和CUDA C之间的互操作性允许程序员只修改少量外围代码，就能用CUDA kerel函数替换Thrust kerel函数，反之亦然。

将Thrust转换成CUDA C很简单，类似于用标准C代码使用C++STL。外部库通过从向量中抽取“原始”指针，可以访问驻留在Thrust容器中的数据。Code2中的代码示例说明了使用原始指针转换，得到指向device_vector内容的整型指针。

Code2

//Thrust 与 CUDA C/C++的互操作

//Thrust dev To CUDA kernel
thrust::device_vector<int> d_vec(1 << 24);

thrust::device_vector<int> dev_Y;

reduction1<int> << <gridDim, threads, threads.x * sizeof(double) >> > (
            thrust::raw_pointer_cast(d_vec.data()),
            temp,
            thrust::raw_pointer_cast(dev_Y.data()));


//CUDA dev To Thrust dev

int* h_test = (int*)malloc((1 << 24) * sizeof(int));

int* d_test;

cudaMemcpy(d_test, h_test, (1 << 24) * sizeof(int),cudaMemcpyHostToDevice)

thrust::device_ptr<int> dev_ptr = thrust::device_pointer_cast(d_test);

thrust::sort(dev_ptr, dev_ptr + (1 << 24));

在Code2中,函数raw_pointer_cast()接受设备向量d_vec的元素0的地址（.data()与STL类似）作为参数,并且返回原始C指针raw_ptr。这个指针可用于调用CUDA C API函数(如cudaMemset()函数)，或者作为参数传递到CUDA C kerel函数中(reduction1函数)。

将 Thrust 算法应用到原始C指针也很简单。一旦原始指针经过 device_ptr 的包装，它便能作为普通的 Thrust迭代器。

Code2中，C指针raw_ptr 指向设备内存中由函数cudaMalloc()分配的一片内存。通过 device_pointer_cast()函数，它可以转换为指向设备向量的设备指针。转换后的指针提供了一些内存空间信息，以便Thrust库调用适当的算法实现，并且为从主机端访问设备存储器提供了方便的机制。在这个例子中，这些信息指明dev_ptr指向设备内存中的向量并且元素类型是整型。

Thrust的原生CUDA C的互操作性保证Thrust总是能作为CUDA C的很好补充,Thrust和CUDA C的结合使用通常比单独使用CUDA C或者Thrust效果好。事实上，即使能够完全使用 Thrust 函数编写完整的并行程序,但是在某些特定领域内直接使用CUDA C实现函数功能会取得更好的结果。原生CUDA C的抽象层次允许程序员能够细粒度地控制计算资源到特定问题的精确映射。在这个层次上编程给开发者提供了实现特定算法的灵活性。互操作性也有利于迭代开发策略:(1)使用Thrust库快速开发出并行应用的原型:(2)确定程序热点;(3)使用CUDA C实现特定算法并作必要优化。

Thrust性能分析

Code

耗时测试代码

#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>

#include <cuda_runtime.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <cstdlib>

#include "helper_cuda.h"
#include "error.cuh"

using namespace std;

const int FORTIME = 50;

template<typename T> __global__
void reduction1(T* X, uint32_t n, T* Y) {
    extern __shared__ uint8_t shared_mem[];
    T* partial_sum = reinterpret_cast<T*>(shared_mem);

    uint32_t tx = threadIdx.x;
    uint32_t i = blockIdx.x * blockDim.x + threadIdx.x;

    partial_sum[tx] = i < n ? X[i] : 0;
    __syncthreads();

    for (uint32_t stride = 1; stride < blockDim.x; stride <<= 1) {
        if (tx % (2 * stride) == 0)
            partial_sum[tx] += tx + stride < n ? partial_sum[tx + stride] : 0;
        __syncthreads();
    }

    if (tx == 0) Y[blockIdx.x] = partial_sum[0];
}


template<typename T>
void rand_array(T* array, size_t len) {

    for (int i = 0; i < len; ++i) {
        array[i] = ((T)rand()) / RAND_MAX;
    }
}

int main(int argc, char* argv[])
{
    thrust::host_vector<int> h_vec(1 << 24);
    cout <<"Test Mem :\t" << (1 << 24) * sizeof(int) / 1024 / 1024 << "MB" << endl;
    thrust::host_vector<int> h_vec1(5);
    thrust::generate(h_vec1.begin(), h_vec1.end(), rand);
    h_vec1[0] = 0;
    h_vec1[4] = 4;
    cout << "h_vec1[4] = \t" << h_vec1[4] << endl << "h_vec1.end() - 1 = \t" << *(h_vec1.end() - 1) << endl;

    thrust::generate(h_vec.begin(), h_vec.end(), rand);

    thrust::device_vector<int> d_vec(1 << 24);
      
    cudaEvent_t start, stop;
    float elapsed_time;
    checkCudaErrors(cudaEventCreate(&start));
    checkCudaErrors(cudaEventCreate(&stop));
    checkCudaErrors(cudaEventRecord(start));
    for (int i = 0; i < FORTIME; i++)
        d_vec = h_vec;

    checkCudaErrors(cudaEventRecord(stop));
    checkCudaErrors(cudaEventSynchronize(stop));
    checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
    std::cout << "thrust HostToDevice elapsed_time:" << elapsed_time / FORTIME << std::endl;

    thrust::sort(d_vec.begin(), d_vec.end());

    checkCudaErrors(cudaEventRecord(start));
    for (int i = 0; i < FORTIME; i++)
        thrust::copy(d_vec.begin(), d_vec.end(), h_vec.begin());
    checkCudaErrors(cudaEventRecord(stop));
    checkCudaErrors(cudaEventSynchronize(stop));
    checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
    std::cout << "thrust Copy DeviceToHost elapsed_time:" << elapsed_time / FORTIME << std::endl;

    checkCudaErrors(cudaEventRecord(start));
    for (int i = 0; i < FORTIME; i++)
        h_vec = d_vec;
    checkCudaErrors(cudaEventRecord(stop));
    checkCudaErrors(cudaEventSynchronize(stop));
    checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
    std::cout << "thrust DeviceToHost elapsed_time:" << elapsed_time / FORTIME << std::endl;

    //-------------------------------------------------------
    int* h_test = (int*)malloc((1 << 24) * sizeof(int));
    int* d_test;

    if (h_test == nullptr)
        return -1;

    rand_array(h_test, 1 << 24);


    checkCudaErrors(cudaMalloc((void**)&d_test, (1 << 24) * sizeof(int) ));

    checkCudaErrors(cudaEventRecord(start));
    for (int i = 0; i < FORTIME; i++)
        checkCudaErrors(cudaMemcpy(d_test, h_test, (1 << 24) * sizeof(int),cudaMemcpyHostToDevice));

    checkCudaErrors(cudaEventRecord(stop));
    checkCudaErrors(cudaEventSynchronize(stop));
    checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
    std::cout << "cudaMemcpy HostToDevice elapsed_time:" << elapsed_time / FORTIME << std::endl;

    checkCudaErrors(cudaEventRecord(start));
    for (int i = 0; i < FORTIME; i++)
        checkCudaErrors(cudaMemcpy(h_test, d_test, (1 << 24) * sizeof(int), cudaMemcpyDeviceToHost));

    checkCudaErrors(cudaEventRecord(stop));
    checkCudaErrors(cudaEventSynchronize(stop));
    checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
    std::cout << "cudaMemcpy DeviceToHost elapsed_time:" << elapsed_time / FORTIME << std::endl;


    //Thrust 与 CUDA C/C++的互操作



    thrust::device_ptr<int> dev_ptr = thrust::device_pointer_cast(d_test);
    thrust::sort(dev_ptr, dev_ptr + (1 << 24));
    
    thrust::device_vector<int> dev_Y;

    dim3 threads(1024);
    dim3 gridDim;
    uint32_t temp = 1 << 24; 
    int sumTime = 0;
    do {

        gridDim = dim3((temp + threads.x - 1) / threads.x);
        d_vec = dev_Y;
        dev_Y.resize(gridDim.x);

        checkCudaErrors(cudaEventRecord(start));

        reduction1<int> << <gridDim, threads, threads.x * sizeof(double) >> > (
            thrust::raw_pointer_cast(d_vec.data()),
            temp,
            thrust::raw_pointer_cast(dev_Y.data()));

        checkCudaErrors(cudaEventRecord(stop));
        checkCudaErrors(cudaEventSynchronize(stop));
        checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
        sumTime += elapsed_time;

        temp = gridDim.x;
    } while (temp > 1);




    free(h_test);
    cudaFree(d_test);

    return 0;
}

具体代码参考Code

可见Thrust的HostToDev、DevToHost和copy()耗时与CUDA C相似。

---

Reduction函数耗时分析：

Thrust虽然方便但是相对于固定优化的CUDA C耗时更长。其它Reduction函数请参考：【CUDA】 归约 Reduction

参考文献：

1、大规模并行处理器编程实战（第2版）

2、​​​CUDA C 编程：基础与实践