使用CUDA計算遞歸數組

Question

我有一個遞歸數組定義。 讓它成爲

A(x, y + 1) = f(A(x - 1, y), A(x, y), A(x + 1, y)) 

第一層被初始化

A(x, 0) = g(x) 

我想通過層中使用CUDA來計算這樣的陣列層。問題是做這種事情的首選方式是什麼。單核應該計算一個數組A(tid, y)爲y在[1, height)同步每一步？還是應該計算單點，但多次調用？或者，最好將問題分解成更大的獨立部分？例如。這個數組可以被rhombs分解，使得如果前一層rhombs完成，每個整個菱形可以被獨立地計算（在菱形內部沒有同步）。

如果圖層是2D而不是1D，情況會不同嗎？

我打算計算寬度〜10000（可能更少就足夠了）和高度44100每秒一個數組。如果有問題，問題實際上是3D（200x50x44100）。爲了簡單起見，我只是製作了2D。

Answer 1

直接的方法可能是剛開始您介紹了這裏什麼：

「如果單核計算[1 y的數組A（TID，Y），高度）上的每一步同步？ 「

這應該很容易實現。

x的「寬度」爲10,000，以保持GPU與該多個線程合理忙碌。

對於複雜的f()函數，能夠每秒執行44100次迭代（平均迭代時間約22 us）可能會很有挑戰性。不過，對於一個相當簡單的f()函數，似乎可能基於下面的快速測試。我們受益於這樣的事實：通過像這樣迭代地啓動內核，隱藏了內核啓動開銷的大部分內容。

這裏是寫在推力展示概念證明樣本代碼：

$ cat t708.cu 
#include <thrust/device_vector.h> 
#include <thrust/host_vector.h> 
#include <thrust/for_each.h> 
#include <thrust/iterator/zip_iterator.h> 
#include <thrust/copy.h> 
#include <stdlib.h> 
#include <iostream> 

#define DSIZE 10000 
#define YSIZE 2 
#define NUM_ITER 44100 
#define AVG_SIZE 3 
#define DISP_WIDTH 5 

struct f 
{ 
    template <typename T> 
    __host__ __device__ 
    void operator()(T t) { 

    thrust::get<AVG_SIZE>(t) = thrust::get<0>(t); 
    thrust::get<AVG_SIZE>(t) += thrust::get<1>(t); 
    thrust::get<AVG_SIZE>(t) += thrust::get<2>(t); 

    thrust::get<AVG_SIZE>(t) /= AVG_SIZE;} 
}; 

int main(){ 

    thrust::host_vector<float> h_A(DSIZE); 
    for (int i =0; i < DSIZE; i++) h_A[i] = rand()/(float)RAND_MAX; // A(x, 0) = g(x) 
    thrust::device_vector<float> d_A[YSIZE]; 
    d_A[0].resize(h_A.size()); 
    d_A[1].resize(h_A.size()); 
    thrust::copy(h_A.begin(), h_A.end(), d_A[0].begin()); 
    thrust::copy(h_A.begin(), h_A.end(), d_A[1].begin()); 
    std::cout << "input left end: " << std::endl; 
    thrust::copy(d_A[0].begin(), d_A[0].begin()+DISP_WIDTH, std::ostream_iterator<float>(std::cout, ",")); 
    std::cout << std::endl << "input right end: " << std::endl; 
    thrust::copy(d_A[0].end() - DISP_WIDTH, d_A[0].end(), std::ostream_iterator<float>(std::cout, ",")); 
    std::cout << std::endl; 

    cudaEvent_t start, stop; 
    cudaEventCreate(&start); cudaEventCreate(&stop); 
    int cur = 0; 
    int nxt = 1; 
    cudaEventRecord(start, 0); 
    for (int i = 0; i < NUM_ITER; i++){ 
    thrust::for_each(thrust::make_zip_iterator(thrust::make_tuple(d_A[cur].begin(), d_A[cur].begin()+1, d_A[cur].begin()+2, d_A[nxt].begin()+1)), thrust::make_zip_iterator(thrust::make_tuple(d_A[cur].end()-2, d_A[cur].end()-1, d_A[cur].end(), d_A[nxt].end()-1)), f()); 
    cur = (cur==0) ? 1:0; // modify for a full storage in y 
    nxt = (nxt==0) ? 1:0;} 
    cudaDeviceSynchronize(); 
    cudaEventRecord(stop, 0); 
    cudaEventSynchronize(stop); 
    float et; 
    cudaEventElapsedTime(&et, start, stop); 
    std::cout << "elapsed time: " << et << "ms" << std::endl << "output left end: " << std::endl; 
    thrust::copy(d_A[cur].begin(), d_A[cur].begin()+DISP_WIDTH, std::ostream_iterator<float>(std::cout, ",")); 
    std::cout << std::endl << "output right end: " << std::endl; 
    thrust::copy(d_A[cur].end() - DISP_WIDTH, d_A[cur].end(), std::ostream_iterator<float>(std::cout, ",")); 
    std::cout << std::endl; 

    return 0; 
} 

$ nvcc -O3 -o t708 t708.cu 
$ ./t708 
input left end: 
0.840188,0.394383,0.783099,0.79844,0.911647, 
input right end: 
0.865333,0.828169,0.311025,0.373209,0.888766, 
elapsed time: 368.337ms 
output left end: 
0.840188,0.838681,0.837174,0.835667,0.83416, 
output right end: 
0.881355,0.883207,0.88506,0.886913,0.888766, 
$ 

注：

1. 整體執行時間爲〜370ms的44100次反覆提示的周圍的平均循環時間8US。
2. 這是在Fedora20系統上運行的，具有CUDA 7和Quadro5000 GPU。
3. 我包括輸入和輸出數據集來驗證「平均」 f()函數的左側和右側的顯示（函子，在這種情況下）。這實際上是一種「放鬆」形式，所以我們可以預期數據集的左端收斂到左端值（不變），數據集的右端收斂到右端值，之間的近似直線。
4. 我只保留Y（n）和Y（N-1）的數據集，用於計算和定時目的。如果你需要的所有Y（0..44100）它的代碼的一個相當簡單的修改。
5. 如果你不熟悉推力，quick start guide可能會讓你加速，你會發現覆蓋了一個有點相似的算法this question，它應該讓你知道這個推力方法可以轉換成一個相當於cuda內核方法。
6. 如果不是很明顯，那麼在調用thrust::for_each時隱含的cuda內核調用將完成「每步同步」，該調用是設備範圍的同步。