【CUDA】 Trust基本特性介紹及性能分析

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 devint* 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;
}