1.   引言

隨著計算機圖形圖像技術的日趨進步，計算機圖形顯示處理器（Graphic Processor Unit）也在不斷的發展。目前先進的顯示處理器已經具備了以下的硬件條件^[1]：

1、    可進行自主編程（the Programmable Graphics Processor Unit）

2、    高并行性、多線程、多處理器（High Parallel、Multithreaded、Manycore processor）

3、    高存儲帶寬(High memory bandwidth)

硬件的發展同時也帶動了顯示處理器的運算速度的大幅提高。自2003年至2008年之間，以Nvidia公司的GForce系列顯卡為例，其每秒浮點數處理速度以及存儲帶寬相比CPU都有了相當大的差距。以下是Nvidia官方給出的硬件評測比較圖表：

圖1 GPU與CPU的浮點數運算速度比較                   圖2 GPU與CPU存儲帶寬比較

可以看出，目前顯示處理器的運算速度以及存儲容量均遠遠強于CPU。這是由于顯示處理器為了更好處理大型渲染計算，因此GPU的設計更傾向于數據計算處理(Data Processing)而不是像CPU那樣來進行數據的緩存(Data Cache)以及流控制(Flow Control)。^[1]因此在GPU的設計中，大部分晶體管被應用到了增加計算吞吐量的功能單元(ALU)。

 圖 3 CPU與GPU的構造區別

2.   通用GPU計算概念與原理

硬件上的區別使得GPU單元更適合大數據量的計算工作，但是從傳統定義上講，GPU的設計主要是用于來計算圖形數據的，因此如何將傳統計算映射到GPU上便成了關鍵的問題。

首先，計算機圖形計算的兩個主要屬性是數據并行性和獨立性，不僅是把相同或者相似的計算過程應用到很多的頂點和片斷組成的流中，而且在每個元素上的計算只有一少部分依賴于其它元素。這種計算模型被稱為離散計算模型。^[1]

其次，目前的GPU都具備了可進行自主編程的功能�？梢酝ㄟ^編寫頂點著色程序(Vertex Shading)和像素著色程序(Pixel Shading)來完成對計算邏輯的控制。^[1]

根據GPU的提供的計算資源，可以將其和CPU計算的資源進行類比，如下表所示：

表 1 GPU計算與CPU計算的映射關系表

通過使用GPU所提供的計算資源進行計算將能獲得使用CPU計算同樣的效果，多維向量以及多處理器則更能提高運算的并行程度與速度。利用GPU進行計算的步驟如下圖所示：

圖 4 通用GPU計算工作流程

利用GPU自身渲染流水線以及可自主變成特點資源進行數據計算的方式即為通用GPU計算的基本概念。^[1]

3.   通用GPU計算的具體實現方法

  雖然通用GPU計算概念不難被接受，但是由于GPU計算對于圖形學知識要求很高，因此目前并沒有得到廣泛的使用。本文將以具體的數學模型為實例，介紹通用GPU計算的實現方式。

3.1.      理想流體粒子計算數學模型

目前許多大型物理實驗都需要使用計算機來進行模擬計算，其中包括氣流模擬、流體模擬等。 這種模擬往往需要計算大量的粒子更新數據，由于粒子之間的計算是相對離散(Scatter Compute)的，因此可以使用通用GPU計算理念來處理這類的數學模型。從而解決傳統計算無法同時處理大密度離散數據處理的瓶頸。^[4]

 圖 5 流體模擬示意圖

粒子的位置(Position)和速度(Velocity)都是需要實時更新的，物體的速度需要由加速度(Acceleration)決定^[3][5],因此根據牛頓第二定律(Newtonian physics)可以使用以下公式：

l                   加速度計算公式：                                             (1)

a為加速度(Acceleration)

F為外力集合(Accumulated force)

M為單獨粒子質量(particle's mass)

l              速度計算公式                                       (2)

a為加速度(Acceleration)

V為更新后速度（Current Velocity）

 為先前速度（previous velocity）

l              粒子位置計算公式：                                     (3)

P為更新后位置（Current Position）

 為先前位置(previous Position) ^[5]

由于需要計算粒子所受力才可以確定粒子的更新，并且流體模擬中粒子都會受到全局外力的影響(Global Force)，外力往往會影響粒子按照共同的速度方向運動，并且對每個粒子的衰減影響均不相同，由于粒子可以近似作為球形體判斷，因此可以使用斯托克方程^[3][5] (Stoke's law of a drag force)來計算流體粘附力。

l              斯托克粘附力方程：                                  (4)

 為流體粘附力(Drag Force)

     為粘附系數(flow viscosity)

     為圓周率

R為粒子半徑(particle's Radius)

V為粒子速度(particle's Velocity)

根據上述公式，可以根據公式(4)計算粒子具體受力，利用(1)(2)來分別計算粒子速度，最后依據公式(3)進行粒子位移的更新。^[5]

3.2.      通用GPU計算模擬粒子的方法

通過CPU將所有粒子的生存周期和初始位置以及速度傳遞給GPU，并使用GPU編程計算粒子的每次更新，粒子每次更新的速度、位置、加速度值使用紋理進行存儲，并提供至GPU的片段處理器反復計算使用，并每次渲染顯示出粒子運動結果。^[2]

3.2.1.         將數據存儲到紋理(Texture)中

根據表格1所示，需要反復計算更新的數據都要存儲到紋理中。紋理具有四個存儲通道，分別為R通道（Red Channel）、G通道（Green Channel）、B通道（Blue Channel）以及Alpha通道（Alpha Channel），并且每一個通道都是一個浮點數(Float)。因此可以將計算需要使用的粒子位置信息，速度信息，加速度信息分別存入紋理的不同通道中，根據粒子數量K 可以創建三張   *  尺寸大小的紋理來存儲計算信息^[1]。如下圖所示：

圖 6 將計算數據存入紋理中

3.2.2.         GPU計算流程介紹

首先需要操作GPU讀入計算的粒子數據，之后通過頂點著色器(Vertex Shader)和像素著色器(Pixel Shader)程序來控制GPU的粒子數據計算更新，并將新數據存儲到紋理中，進行下一次GPU計算，同時顯示粒子更新效果。工作過程如下所示：

圖 7 工作流程圖示

3.3.      具體實現設計與開發

開發工具使用清單如表2所示

表 2 開發環境清單

由于需要進行GPU操作，因此需要直接進行顯卡計算指令的操作。本文采用了NVIDIA Cg語言進行編寫，最后通過編譯器轉換為顯卡指令集。^[2]

3.3.1.         主程序邏輯

    由于所有的粒子的數據計算均放到GPU中處理，因此CPU只需處理粒子的初始化，更新邏輯以及一些數據結構定義，如下圖：

圖 8 程序組成結構圖

表 3 程序結構說明

通過ParticleManager提供程序主入口完成3D繪制設備初始化并創建粒子系統對象，粒子系統對象在創建完成后會自我進行更新，這個過程會不斷通過OpenGl API調用GPU進行粒子計算。^[4]

int main(int argc, char **argv)

{

    glutInit(&argc, argv);

    glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);

    glutInitWindowSize(512, 512);

    (void) glutCreateWindow("gpu particle system");

    initGL();

    int w = 256, h = 256;

    psys = new ParticleSystem(w, h, 1, false);

    glutMainLoop();

    return 0;

}

3.3.2.         GPU計算

首先將初始化的粒子位置和速度錄入進用于存儲數據使用的紋理中。按照粒子的索引來向指定的貼圖通道中寫入粒子數據。^[4]

glBegin(GL_POINTS);

for (i = m_newParticles.begin(); i != m_newParticles.end(); i++)

{

        vec2f pos = IndexToPos(i->index);

        glTexCoord4fv(&i->pos[0]);

        glVertex2fv(&pos[0]);

    }

glEnd();

m_pos_buffer[m_previous]->Deactivate();

初始化結束后，通過編寫頂點著色器和像素著色器程序根據紋理中的位置和速度數據，在每次像素著色過程中重新計算粒子更新數據。以上操作使用Cg編寫完成。^[2]

float3 pos = texRECT(pos_tex, uv).xyz;

    float3 vel = texRECT(vel_tex, uv).xyz;

    float3 pos_next = pos + vel * timestep; 

    float3 force = gravity;

FloorCollide(pos_next,vel,0.0,0.5,forc)

vel += force * inv_mass * timestep;

 vel *= damping;

將更新的位置和速度信息重新存入紋理中，并循環至下次使用。同時調用GPU進行繪制。

// Vertex Shader

const float4 pointSize = 64.0;

    const float3 atten = float3(1.0, 0.0, 1.0);

    float4 pos_eye = mul(modelView, pos);

    float d = length(pos_eye.xyz);

psize = pointSize * sqrt(1.0 / (atten[0] + atten[2]*d*d));

// Pixel Shader

return tex2D(cloth_tex, texcoord0.xy) * color;

3.4.      測試結果與比較

最終的程序同時模擬了512 * 512個粒子的數據計算，具體測試環境如表4所示

表 4 測試環境清單

為了進行計算效率比較，分別采用了GPU計算和傳統CPU計算兩種方式來進行測試。結果如圖9所示。

圖 9 GPU計算與傳統計算速度比較

該圖表為計算一次粒子所需要時間顯示�？梢钥闯�，利用GPU進行離散類型數學計算速度要比傳統的CPU計算更高效，在多核并行計算的前提下則更能體現出GPU的計算優勢。利用GPU提高大規模計算速度也將會成為今后一段時間的主要手段。

3.5.      結束語

本論文主要闡述了GPU計算與傳統CPU計算的區別與比較，并以一個實例分析了GPU計算的可行性，此外還介紹了如何使用GPU計算的方法。通過模擬粒子流體的計算可以看出，通用GPU計算相對于傳統的CPU計算能夠更好的處理大規模離散性數據。

參考文獻(References)

[1] NVIDIA Corporation; GPU Gems 2: Programming Techniques for High-Performance Graphics and General-Purpose Computation

2005 http://developer.nvidia.com/object/gpu_gems_2_home.html

[2]Mark, William R.; Glanville, R. Steven; Akeley, Kurt; Kilgard, Mark J.; Cg: A System for Programming Graphics Hardware in a C-like Language. In SIGGRAPH Proceedings 2003

[3]Reeves, William T.; Particle Systems – Technique for Modeling a Class of Fuzzy

Objects. In SIGGRAPH Proceedings 1983

[4] NVIDIA Corporation; OpenGL Extension NV_pixel_data_range, 2002,

http://oss.sgi.com/projects/ogl-sample/registry/NV/pixel_data_range.txt

[5]Jakobsen, Thomas; Advanced Character Physics. In GDC Proceedings 2001

Lang2003: Lang, Hans W.; Odd-Even Merge Sort, 2003,

http://www.iti.fhflensburg.de/lang/algorithmen/sortieren/oemen.html

作者簡介：曹譞 24歲 2003年北京工業大學就讀于計算機學院并于2007年獲取本科學士，2007年就讀于北京工業大學計算機學院碩士。2007年代表中國參加羅馬尼亞舉辦Hard&Soft計算機比賽并獲取第一名，2007 Dream Build Play XNA 游戲開發比賽第三名成員。09834