文章連結： http://blog.csdn.net/u013108511/article/details/77096661
郵箱：[email protected]
  
  上一個部落格中我們瞭解到了在函式前面加__ global_ _修飾符來定義核函式並通過一個尖括號來呼叫就可以在GPU上執行這個函式的功能。方法簡單但是低效，因為經過優化的圖形處理器可以並行執行數百次的運算，然而在例項中只調用一個核函式且是序列執行。下面來了解運算的並行性。

##1、向量求和運算
  向量的求和運算，這裡對運算的模型進行簡單化，就理解成兩個陣列對應下標的資料元素相加，並將結果儲存到第三個陣列中。

###1.1、基於CPU的向量求和

  下面用傳統的標準的C方法實現兩個陣列相加的小程式.

#include <stdio.h>
#define N 10

void add(int *a, int *b, int *c) {
    int tid = 0;//第0個CPU，因此索引從0開始

    while(tid < N) {
        c[tid] = a[tid] + b[tid];
        tid ++;//單執行緒程式設計，由於只有一個CPU，因此每次遞增一
    }
}

int main( void ) {

    int a[N],b[N],c[N];

    //在CPU上位陣列a和b賦值
    for(int i = 0; i < N; i++) {
        a[i] = -i;
        b[i] = i * i;
    }

    add(a,b,c);

    //顯示結果
    for(int i = 0; i < N; i++) {
        printf("%d + %d = %d\n", a[i], b[i], c[i]);
    }

    return 0;
}
 

更簡單的相加函式可以通過for迴圈實現，同時如果想使用兩個CPU並行執行的話，策略是一個CPU執行下標為奇數的加法運算，一個CPU執行下標為偶數的加法運算。程式碼實現如下：

//第一個CPU核心
void add (int *a, int *b, int *c){
	int tid = 0;
	while(tid < N) {
		c[tid] = a[tid] + b[tid];
		tid += 2;
	}
}

//第二個CPU核心
void add (int *a, int *b, int *c){
	int tid = 1;
	while(tid < N) {
		c[tid] = a[tid] + b[tid];
		tid += 2;
	}
}
 

  然而即便如此，CPU按以上執行仍然需要增加更多的程式碼，如多執行緒的定義和控制等，但通常執行緒的實際執行情況並非如此。

##1.2、基於GPU的向量求和
  實現加法運算只需要將add（）編寫為裝置函式即可。

#include "../common/book.h"

#define N 65500

__global__ void add(int *a, int *b, int *c) {
    int tid = blockIdx.x;   
    if(tid < N)
        c[tid] = a[tid] + b[tid];
}

int main( void ) {
    int a[N], b[N], c[N];
    int *dev_a, *dev_b, *dev_c;
    
    //在GPU上分配記憶體
    HANDLE_ERROR(cudaMalloc ((void **)&dev_a, N*sizeof(int)));
    HANDLE_ERROR(cudaMalloc ((void **)&dev_b, N*sizeof(int)));
    HANDLE_ERROR(cudaMalloc ((void **)&dev_c, N*sizeof(int)));

    //CPU上位陣列a和b賦值
    for(int i = 0; i < N; i++) {
        a[i] = -i;
        b[i] =  i* 2;
    }

    //將a和b陣列複製到GPU中
    //注意：至於為什麼不採用直接在裝置上賦值以獲得更快的執行速度，這裡只是假設，後面會繼續討論
    HANDLE_ERROR(cudaMemcpy (dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
    HANDLE_ERROR(cudaMemcpy (dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));

    add<<<N, 1>>> (dev_a, dev_b, dev_c);

    //將陣列c從GPU拷貝到CPU
    HANDLE_ERROR(cudaMemcpy (c, dev_c, N *sizeof(int), cudaMemcpyDeviceToHost));

    //顯示結果
    for(int i = 0; i < N; i++) {
        printf("%d + %d = %d\n", a[i], b[i], c[i]);
    }

    //釋放在GPU上分配的記憶體
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return 0;
}

GPU程式設計有一些通用的模式，初學CUDA程式設計的時候熟悉模式可以快速上手程式設計，模式有一下幾點：
  1、GPU上的變數使用前都要為其分配記憶體，本例子中呼叫cudaMalloc（）函式為三個陣列分配記憶體，兩個陣列包含輸入值，一個數組包含結果；
  2、為了避免記憶體洩漏，使用結束後通過cudeFree（）函式釋放記憶體；
  3、通過cudaMemcpy()進行主機和記憶體之間資料的拷貝，通過指定引數cudaMemcpyHostToDevice將資料從主機拷貝到裝置中，通過制定引數cudaMemcpyDeviceToHost將資料從裝置拷貝到主機中。
  4、通過尖括號語法，在主機程式碼的main（）函式中執行add（）中的裝置程式碼。
  5、而編寫一個在裝置上執行的函式，只需要在函式的前面宣告一個修飾符__ global_ _
注意：至於為什麼不採用直接在裝置上賦值以獲得更快的執行速度，這裡只是假設，後面會繼續討論

上一篇部落格出現一下形式來啟動核函式：
  kernel<<<1,1>>>(param1，param2，…)
本例項中出現的核函式：
  add<<< N,1>>>（dev_a, dev_b, dev_c）

注意尖括號中的兩個引數，其中第一個引數表示裝置在執行核函式的時候使用的並行執行緒塊的數量。這個例項中指定這個引數為N。

###獲取程式碼當前執行的執行緒塊：變數blockIdx.x
  blockIdx是一個內建的變數，所以不需要定義，CUDA執行時已經預先定義，變數中包含的值就是當前執行裝置程式碼的執行緒塊的索引。使用blockIdx.x是因為CUDA支援二維的執行緒塊的陣列，使用二維的索引會帶來很大的便利。N個執行緒塊組成一個執行緒格（Grid），本例中一個一維的執行緒格包含N個執行緒塊，每個執行緒塊的blockIdx.x都不相同，第一個執行緒塊的blockIdx.x為0.最後一個執行緒塊的blockIdx.x為N-1。
  在啟動執行緒塊陣列時，陣列的每一維的最大數量都不能超過65535，這是一種硬體的限制，如果啟動的硬體的執行緒塊超過了這個限制，程式將會執行失敗。

#2、一個有趣的例子：繪製Julia集的曲線
  Julia集曲線是滿足某個複數計算函式的所有點構成的邊界。Julia集的基本演算法是，通過一個簡單的迭代等式對複平面中的點進行求值，如果早計算某一點時，迭代式的結果是發散的那麼這個點就不屬於Julia集合，相反如果迭代等式中計算的一系列的值位於某個邊界範圍以內，那麼這個點就是Julia集。
  具體的公式如下：
                   

##2.1、CPU的程式設計實現

#include "stdio.h"
#include "../common/cpu_bitmap.h"

//定義了畫素點長和寬的大小
#define DIM 1000

//複數的資料結構的定義以及對這個資料結構可以執行的操作。
struct cuComplex {
    float r;
    float i;
    cuComplex(float a, float b) {
        r = a;
        i = b;
    }
	//計算向量的長度的平方
    float magnitude2(void) {
        return r * r + i * i;
    }
	//乘法運算子過載
    cuComplex operator *(const cuComplex& a) {
        return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
    }
	//加法的運演算法的過載
    cuComplex operator +(const cuComplex& a) {
        return cuComplex(r+a.r, i+a.i);
    }
};

//判斷複數x+yi是否屬於Julia集，執行上面的操作
int julia(int x, int y) {
    const float scale = 1.5;

    float jx = scale * (float)(DIM/2 - x)/(DIM/2);

    float jy = scale * (float)(DIM/2 - y)/(DIM/2);

    cuComplex c(-0.8, 0.156);
    cuComplex a(jx, jy);

    int i = 0;
    for(i = 0; i < 200; i++) {
        a = a * a + c;
        if(a.magnitude2() > 1000)
            return 0;
    }
    return 1;
}

//核函式，對於設定的1000*1000畫素的區域，每個點是否是Julia集進行判斷
void kernel( unsigned char *ptr) {
    for(int y=0; y<DIM; y++) {
        for(int x=0; x<DIM; x++) {

            int offset = x + y * DIM;

            int juliaValue = julia(x, y);

            ptr[offset*4 + 0] = 255 *juliaValue;
            ptr[offset*4 + 1] = 0;
            ptr[offset*4 + 2] = 0;
            ptr[offset*4 + 3] = 255;
        }
    }
}

//主函式，由於是面向CPU的程式設計，程式全部在CPU上執行
int main( void ) {

    CPUBitmap bitmap( DIM, DIM );
    
    unsigned char *ptr = bitmap.get_ptr();

    kernel( ptr );
    
    bitmap.display_and_exit();

    return 0;
}

注意：由於一些庫和工具的引入，程式的編輯可能出現編譯錯誤的情況，CPU程式碼編譯指令如下：

g++ CPUJulia.c -lGL -lGLU -lglut -o CPUJulia

之所以用g++，應為長鬚中有一些c不支援的庫，g++編譯沒有問題。

##2.2、GPU的程式設計實現

#include "cuda.h"
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "../common/book.h"
#include "../common/cpu_bitmap.h"

#define DIM 1000

//複數的資料結構，__device__表示在裝置上執行的程式
//同時要注意的是:在裝置上執行的額程式在執行函式呼叫的時候也只能呼叫在裝置上執行的程式部分，而不能呼叫在主機上執行的函式部分
struct cuComplex {

    float r;
    float i;
    __device__ cuComplex(float a, float b) : r(a), i(b){}
    __device__ float magnitude2(void) {
        return r * r + i * i;
    }

    __device__ cuComplex operator *(const cuComplex& a) {
        return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
    }

    __device__ cuComplex operator +(const cuComplex& a) {
        return cuComplex(r+a.r, i+a.i);
    }
};

__device__ int julia(int x, int y) {
    const float scale = 1.5;

    float jx = scale * (float)(DIM/2 - x)/(DIM/2);

    float jy = scale * (float)(DIM/2 - y)/(DIM/2);

    cuComplex c(-0.8, 0.156);
    cuComplex a(jx, jy);

    int i = 0;
    for(i = 0; i < 200; i++) {
        a = a * a + c;
        if(a.magnitude2() > 1000)
            return 0;
    }
    return 1;
}

__global__ void kernel( unsigned char *ptr) {

    //將threadIdx/BlockIdx對映到畫素位置
    int x = blockIdx.x;
    int y = blockIdx.y;
    int offset = x + y * gridDim.x;
    
    //現在計算這個位置上的俄值
    int juliaValue = julia(x, y);
    ptr[offset*4 + 0] = 255 *juliaValue;
    ptr[offset*4 + 1] = 0;
    ptr[offset*4 + 2] = 0;
    ptr[offset*4 + 3] = 255;
}
struct DataBlock {  
    unsigned char   *dev_bitmap;  
}; 

int main( void ) {  
    DataBlock   data;  
    CPUBitmap bitmap( DIM, DIM, &data );  
    unsigned char    *dev_bitmap;  
  
    HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() ) );  
    data.dev_bitmap = dev_bitmap;  
  
    dim3    grid(DIM,DIM);  
    kernel<<<grid,1>>>( dev_bitmap );  
  
    HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,  
                              bitmap.image_size(),  
                              cudaMemcpyDeviceToHost ) );  
                                
    HANDLE_ERROR( cudaFree( dev_bitmap ) );  
                                
    bitmap.display_and_exit();  
}  

注意：由於一些庫和工具的引入，程式的編輯可能出現編譯錯誤的情況，GPU程式碼編譯指令如下：

nvcc GPUJulia.c -lGL -lGLU -lglut -o GPUJulia

同時還要注意，CPU高效能程式設計CUDA實戰中的程式碼編譯無法通過，要在cuComplex的建構函式前面加上_ _ device_ _

  對於上面的程式，前面對功能的理解完成之後，我們將學習重點放在GPU的並行程式設計的實現上面，上面的例子主要用到了GPU並行程式設計的一下知識點，逐漸的積累過程來加深對GPU並行程式設計的認識：
 1、整個程式執行的總的思路：main（）函式通過一個工具庫建立一個DIMDIM大小的點陣圖影象。
 2、對於10001000的點陣圖影象，每個畫素點代表一個對應的複數的值，每個值是否屬於Julia集使我們要在裝置上進行判斷的，也是我們程式的核心部分；
 3、程式中指定了多個執行緒塊來執行函式Kennel，由於每個畫素點的計算和其他的點的計算相互的獨立，因此每個計算的點都執行該函式的副本，這是對於DIM*DIM個點，需要二維的執行緒索引，這時就有：
 
                dim3 grid(DIM,DIM)
這是一個3維的陣列，但是當只用兩個值來初始化dim3的時候，第三個維度會自動的初始化為1。NVIDIA將來可能會支援三維的索引，目前只能這樣做來滿足API的要求。
 4、kernel函式被宣告為一個_ global 的函式，這樣這個函式可以再主機上呼叫並且在裝置上執行。
 5、其他的裝置上執行的函式與kernel不同的是：都**宣告為一個 device_ 函式，這表示程式碼在裝置上執行，另外與 _global _不同的是聲明瞭 device __ 的函式只能從其他的 _device _或者 _global __ 函式呼叫而無法通過主機函式的呼叫。**

##2.3、顯示結果

#3、總結
  通過本節，一個更復雜一點的例子，更深入的瞭解並行程式設計的實質。並行程式設計注意一下兩點：
（1）注意套用CUDA程式設計的模式，這樣可以將主要的注意力集中在並行執行緒的設計以及演算法的實現上來；
（2）注意程序的設定，加快程式理論執行的速率和效率。