異構計算（CPU + GPU）編程簡介
1. 概念
所謂異構計算，是指CPU+ GPU或者CPU+ 其它設備（如FPGA等）協同計算。一般我們的程序，是在CPU上計算。
但是，當大量的數據需要計算時，CPU顯得力不從心。那么，是否可以找尋其它的方法來解決計算速度呢？那就是異構計算。
例如可利用CPU（Central Processing Unit）、GPU（Graphic Processing Unit）、甚至APU(Accelerated Processing Units， CPU與GPU的融合)
等計算設備的計算能力從而來提高系統的速度。異構系統越來越普遍，對於支持這種環境的計算而言，也正受到越來越多的關注。
2. 異構計算的實現
目前異構計算使用最多的是利用GPU來加速。主流GPU都采用了統一架構單元，憑借強大的可編程流處理器陣容，
GPU在單精度浮點運算方面將CPU遠遠甩在身后。英特爾Core i7 965處理器，在默認情況下，
它的浮點計算能力只有NVIDIA GeForce GTX 280 的1/13，與AMD Radeon HD 4870相比差距就更大。
3.基於GPU編程
不同廠商通常僅僅提供對於自己設備編程的實現。對於異構系統一般很難用同種風格的編程語言來實現機構編程，
而且將不同的設備作為統一的計算單元來處理的難度也是非常大的。基於GPU編程的，目前主要兩大廠商提供：
一個是NVidia，其提供的GPU編程為CUDA，目前使用的CUDA SDK 4.2.
另一個是AMD，其提供的GPU編程為AMD APP （其前身是ATI Stream），目前最新版本 AMD APP 2.7。
這兩個東東是不兼容的，各自為政。作為軟件開發者而言，用CUDA開發的軟件只能在NVidia相應的顯卡上運行，
用AMD APP開發的軟件，只能在ATI相應的顯卡上運行。
4. OpenCL簡介
那么有沒有可能讓他們統一起來，簡化編程呢？有，那就是由蘋果公司發起並最后被業界認可的OpenCL，目前版本1.2。
開放式計算語言（Open Computing Language:OpenCL），旨在滿足這一重要需求。通過定義一套機制，
來實現硬件獨立的軟件開發環境。利用OpenCL可以充分利用設備的並行特性，支持不同級別的並行，
並且能有效映射到由CPU，GPU， FPGA（Field－Programmable Gate Array）和將來出現的設備所組成的同構或異構，
單設備或多設備的系統。OpenCL定義了運行時， 允許用來管理資源，將不同類型的硬件結合在同種執行環境中，
並且很有希望在不久的將來，以更加自然的方式支持動態地平衡計算、功耗和其他資源。
5. DirectCompute簡介
作為軟件行業的老大—微軟在這方面又做了什么呢？微軟也沒閑着，微軟推出DirectCompute，與OpenCL抗衡。
DirectCompute集成在DX中，目前版本DX11，其中就包括DirectCompute。由於微軟的地位，所以大多數廠商也都支持DirectCompute。
6. GPU計算模型
內核是執行模型的核心，能在設備上執行。當一個內核執行之前，需要指定一個 N-維的范圍（NDRange）。
一個NDRange是一個一維、二維或三維的索引空間。 還需要指定全局工作節點的數目，
工作組中節點的數目。如圖NDRange所示，全局工作節點的范圍為{12, 12}，工作組的節點范圍為{4, 4}，總共有9個工作組。
如果定義向量為1024維，特別地，我們可以定義全局工作節點為1024，工作組中節點為128，
則總共有8個組。定義工作組主要是為有些僅需在組內交換數據的程序提供方便。當然工作節點數目的多少要受到設備的限制。
如果一個設備有1024個處理節點，則1024維的向量，每個節點計算一次就能完成。而如果一個設備僅有128個處理節點，
那么每個節點需要計算8次。合理設置節點數目，工作組數目能提高程序的並行度。
7. 程序實例
不論是OpenCL還是DirectCompute，其編程風格都基本差不多，程序是分成兩部分的：
一部分是在設備上執行的（對於我們，是GPU），另一部分是在主機上運行的（對於我們，是CPU）。
在設備上執行的程序或許是你比較關注的。它是OpenCL和DirectCompute產生神奇力量的地方。
為了能在設備上執行代碼，OpenCL程序員需要寫一個特殊的函數（kernel函數）放在專用文件中（.cl），這個函數需要使用OpenCL語言編寫。
OpenCL語言采用了C語言的一部分加上一些約束、關鍵字和數據類型。在主機上運行的程序提供了API，所以可以管理你在設備上運行的程序。
主機程序可以用C或者C++編寫，它控制OpenCL的環境（上下文，指令隊列…）。
DirectCompute程序員需要寫Shader文件（.hlsl），在這個文件中寫函數。Shader文件的格式可以查MSDN。
在寫程序時，先初始化設備，然后編譯需要在GPU上運行的程序（運行在GPU上的程序是在應用程序運行時編譯的）。
然后映射需要在GPU上運行的函數名字，OpenCL 調用clEnqueueNDRangeKernel執行kernel函數，
DirectCompute調用ID3D11DeviceContext:: Dispatch執行Shader函數。函數是並發執行的。
運行在GPU上的函數一般都很簡單。以求和為例：
用CPU運算
- void vector_add_cpu (const float* fIn1,
- const float* fIn2,
- float* fOut,
- const int iNum)
- {
- for (int i = 0; i < iNum; i++)
- {
- fOut [i] = fIn1[i] + fIn2[i];
- }
- }
void vector_add_cpu (const float* fIn1,
const float* fIn2,
float* fOut,
const int iNum)
{
for (int i = 0; i < iNum; i++)
{
fOut [i] = fIn1[i] + fIn2[i];
}
}
以下是OPenCL的kernel函數
- //在GPU上，邏輯就會有一些不同。我們使每個線程計算一個元素的方法來代替cpu程序中的循環計算。每個線程的index與要計算的向量的index相同。
- __kernel void vector_add_gpu (__global const float* fIn1,
- __global const float* fIn2,
- __global float* fOut,
- const int iNum)
- {
- /* get_global_id(0) 返回正在執行的這個線程的ID。
- 許多線程會在同一時間開始執行同一個kernel，
- 每個線程都會收到一個不同的ID，所以必然會執行一個不同的計算。*/
- const int idx = get_global_id(0);
- /* 每個work-item都會檢查自己的id是否在向量數組的區間內。
- 如果在，work-item就會執行相應的計算。*/
- if (idx < iNum)
- {
- fOut [idx] = fIn1[idx] + fIn2[idx];
- }
- }
//在GPU上，邏輯就會有一些不同。我們使每個線程計算一個元素的方法來代替cpu程序中的循環計算。每個線程的index與要計算的向量的index相同。
__kernel void vector_add_gpu (__global const float* fIn1,
__global const float* fIn2,
__global float* fOut,
const int iNum)
{
/* get_global_id(0) 返回正在執行的這個線程的ID。
許多線程會在同一時間開始執行同一個kernel，
每個線程都會收到一個不同的ID，所以必然會執行一個不同的計算。*/
const int idx = get_global_id(0);
/* 每個work-item都會檢查自己的id是否在向量數組的區間內。
如果在，work-item就會執行相應的計算。*/
if (idx < iNum)
{
fOut [idx] = fIn1[idx] + fIn2[idx];
}
}
有一些需要注意的地方：
1. Kernel關鍵字定義了一個函數是kernel函數。Kernel函數必須返回void。
2. Global關鍵字位於參數前面。它定義了參數內存的存放位置。
另外，所有kernel都必須寫在“.cl”文件中，“.cl”文件必須只包含OpenCL代碼。
Shader函數
- #define NUM_THREAD 16
- StructuredBuffer<float> fInput1 : register( t0 );
- StructuredBuffer<float> fInput2 : register( t0 );
- StructuredBuffer<float> fOutput : register( u0 );
- [numthreads(NUM_THREAD, 1, 1)]
- void vector_add_gpu( uint3 Gid : SV_GroupID,
- uint3 DTid : SV_DispatchThreadID,
- uint3 GTid : SV_GroupThreadID,
- uint GI : SV_GroupIndex )
- {
- fOutput[DTid.x] = fInput1[DTid.x] + fInput2[DTid.x] ;
- }
#define NUM_THREAD 16
StructuredBuffer<float> fInput1 : register( t0 );
StructuredBuffer<float> fInput2 : register( t0 );
StructuredBuffer<float> fOutput : register( u0 );
[numthreads(NUM_THREAD, 1, 1)]
void vector_add_gpu( uint3 Gid : SV_GroupID,
uint3 DTid : SV_DispatchThreadID,
uint3 GTid : SV_GroupThreadID,
uint GI : SV_GroupIndex )
{
fOutput[DTid.x] = fInput1[DTid.x] + fInput2[DTid.x] ;
}
圖像旋轉是指把定義的圖像繞某一點以逆時針或順時針方向旋轉一定的角度，通常是指繞圖像的中心以逆時針方向旋轉。
假設圖像的左上角為(l, t), 右下角為(r, b)，則圖像上任意點(x, y) 繞其中心(xcenter, ycenter)逆時針旋轉θ角度后， 新的坐標位置(x',y')的計算公式為：
x′ = (x - xcenter) cosθ - (y － ycenter) sinθ + xcenter,
y′ = (x - xcenter) sinθ + (y － ycenter) cosθ + ycenter.
C代碼：
- void image_rotate(
- unsigned int* iInbuf,
- unsigned int* iOutbuf,
- int iWidth, int iHeight,
- float fSinTheta,
- float fCosTheta)
- {
- int i, j;
- int iXc = iWidth /2;
- int iYc = iHeight /2;
- for(i = 0; i < iHeight; i++)
- {
- for(j=0; j< iWidth; j++)
- {
- int iXpos = (j- iXc)*fCosTheta - (i - iYc) * fSinTheta + iXc;
- int iYpos = (j- iXc)*fSinTheta + (i - iYc) * fCosTheta + iYc;
- if(iXpos >=0&& iYpos >=0&& iXpos < iWidth && iYpos < iHeight)
- iOutbuf[iYpos * iWidth + iXpos] = iInbuf[i* iWidth +j];
- }
- }
- }
void image_rotate(
unsigned int* iInbuf,
unsigned int* iOutbuf,
int iWidth, int iHeight,
float fSinTheta,
float fCosTheta)
{
int i, j;
int iXc = iWidth /2;
int iYc = iHeight /2;
for(i = 0; i < iHeight; i++)
{
for(j=0; j< iWidth; j++)
{
int iXpos = (j- iXc)*fCosTheta - (i - iYc) * fSinTheta + iXc;
int iYpos = (j- iXc)*fSinTheta + (i - iYc) * fCosTheta + iYc;
if(iXpos >=0&& iYpos >=0&& iXpos < iWidth && iYpos < iHeight)
iOutbuf[iYpos * iWidth + iXpos] = iInbuf[i* iWidth +j];
}
}
}
CL代碼：
- __kernel void image_rotate(
- __global int * iInbuf,
- __global int * iOutbuf, //Data in global memory
- int iWidth , int iHeight, //Image Dimensions
- float fSinTheta, float fCosTheta ) //Rotation Parameters
- {
- const int ix = get_global_id(0);
- const int iy = get_global_id(1);
- int iXc = iWidth /2;
- int iYc = iHeight /2;
- int iXpos = ( ix- iXc)*fCosTheta - (iy- iYc)*fSinTheta+ iXc;
- int iYpos = (ix- iXc)*fSinTheta + ( iy- iYc)*fCosTheta+ iYc;
- if ((iXpos >=0) && (iXpos < iWidth ) && (iYpos >=0) && (iYpos < iHeight))
- iOutbuf[iYpos * iWidth + iXpos]= iInbuf[iy* iWidth +ix];
- }
__kernel void image_rotate(
__global int * iInbuf,
__global int * iOutbuf, //Data in global memory
int iWidth , int iHeight, //Image Dimensions
float fSinTheta, float fCosTheta ) //Rotation Parameters
{
const int ix = get_global_id(0);
const int iy = get_global_id(1);
int iXc = iWidth /2;
int iYc = iHeight /2;
int iXpos = ( ix- iXc)*fCosTheta - (iy- iYc)*fSinTheta+ iXc;
int iYpos = (ix- iXc)*fSinTheta + ( iy- iYc)*fCosTheta+ iYc;
if ((iXpos >=0) && (iXpos < iWidth ) && (iYpos >=0) && (iYpos < iHeight))
iOutbuf[iYpos * iWidth + iXpos]= iInbuf[iy* iWidth +ix];
}
不論是OpenCL還是 DirectCompute，其編程還是有些復雜，特別是對於設備的初始化，以及數據交換，非常麻煩。
對於初學者難度相當大。那么有沒有更簡單的編程方法呢？
8. C++AMP
還要提到微軟，因為我們基本上都使用微軟的東東。微軟也不錯，推出了C++AMP，這是個開放標准，嵌入到VS2012中（VS2012目前還是預覽版）
，在Win7和Win8系統中才能使用，其使用就簡單多了：
- void CppAmpMethod()
- {
- int aCPP[] = {1, 2, 3, 4, 5};
- int bCPP[] = {6, 7, 8, 9, 10};
- int sumCPP[size];
- // Create C++ AMP objects.
- array_view<const int, 1> a(size, aCPP);
- array_view<const int, 1> b(size, bCPP);
- array_view<int, 1> sum(size, sumCPP);
- sum.discard_data();
- parallel_for_each(
- // Define the compute domain, which is the set of threads that are created.
- sum.extent,
- // Define the code to run on each thread on the accelerator.
- [=](index<1> idx) restrict(amp)
- {
- sum[idx] = a[idx] + b[idx];
- }
- );
- }
void CppAmpMethod()
{
int aCPP[] = {1, 2, 3, 4, 5};
int bCPP[] = {6, 7, 8, 9, 10};
int sumCPP[size];
// Create C++ AMP objects.
array_view<const int, 1> a(size, aCPP);
array_view<const int, 1> b(size, bCPP);
array_view<int, 1> sum(size, sumCPP);
sum.discard_data();
parallel_for_each(
// Define the compute domain, which is the set of threads that are created.
sum.extent,
// Define the code to run on each thread on the accelerator.
[=](index<1> idx) restrict(amp)
{
sum[idx] = a[idx] + b[idx];
}
);
}
就這么簡單，只需要包含一個頭文件，使用一個命名空間，包含庫文件，一切就OK，在調用時，只是一個函數parallel_for_each。（有點類似OpenMP）。
9. 應用
MATLAB 2010b中Parallel Computing Toolbox與MATLAB Distributed Computing Server的最新版本可利用NVIDIA的CUDA並行計算架構
在NVIDIA計算能力1.3以上的GPU上處理數據，執行GPU加速的MATLAB運算，將用戶自己的CUDA Kernel函數集成到MATLAB應用程序當中。
另外，通過在台式機上使用Parallel Computing Toolbox以及在計算集群上使用MATLAB Distributed Computing Server來運行多個MATLAB worker程序，
從而可在多顆NVIDIA GPU上進行計算。AccelerEyes公司開發的Jacket插件也能夠使MATLAB利用GPU進行加速計算。Jacket不僅提供了GPU API（應用程序接口），而且還集成了GPU MEX功能。在一定程度說，Jacket是一個完全對用戶透明的系統，能夠自動的進行內存分配和自動優化。Jacket使用了一個叫“on-the- fly”的編譯系統，使MATLAB交互式格式的程序能夠在GPU上運行。目前，Jacket只是基於NVIDIA的CUDA技術，但能夠運行在各主流操作系統上。
從Photoshop CS4開始，Adobe將GPU通用計算技術引入到自家的產品中來。GPU可提供對圖像旋轉、縮放和放大平移這些常規瀏覽功能的加速，
還能夠實現2D/3D合成，高質量抗鋸齒，HDR高動態范圍貼圖，色彩轉換等。而在Photoshop CS5中，更多的算法和濾鏡也開始支持GPU加速。
另外，Adobe的其他產品如Adobe After Effects CS4、Adobe Premiere Pro CS4也開始使用GPU進行加速。這些軟件借助的也是NVIDIA的CUDA技術。
Windows 7 的核心組成部分包括了支持GPU通用計算的Directcompute API，為視頻處理、動態模擬等應用進行加速。
Windows 7借助Directcompute增加了對由GPU支持的高清播放的in-the-box支持，可以流暢觀看，同時CPU占用率很低。
Internet Explorer 9加入了對Directcompute技術的支持，可以調用GPU對網頁中的大計算量元素做加速計算；
Excel2010、Powerpoint2010也開始提供對Directcompute技術的支持。
比利時安特衛普大學，通用電氣醫療集團，西門子醫療，東芝中風研究中心和紐約州立大學水牛城分校的都針對GPU加速CT重建進行了各自的研究，
不僅如此，西門子醫療用GPU實現了加速MRI中的GRAPPA自動校准，完成MR重建，快速MRI網格化，
隨機擴散張量磁共振圖像（DT-MRI）連通繪圖等算法。其他的一些研究者則把醫學成像中非常重要的二維與三維圖像中器官分割（如Level Set算法），
不同來源圖像的配准，重建體積圖像的渲染等也移植到GPU上進行計算。
10 .不足
異構並行計算變得越來越普遍，然而對於現今存在的OpenCL和DirectCompute版本來說，的確還存在很多不足，例如編寫內核，
需要對問題的並行情況做較為深入的分析，對於內存的管理還是需要程序員來顯式地申明、顯式地在主存和設備的存儲器之間進行移動，
還不能完全交給系統自動完成。從這些方面，OpenCL和 DirectCompute的確還需加強，要使得人們能高效而又靈活地開發應用，還有很多工作要完成。