用GPU通用并行計(jì)算繪制曼德勃羅特集圖形 下篇

上篇中我們用DirectX Compute Shader在顯卡上編寫了一個(gè)并行算法來計(jì)算好看的曼德勃羅特集迭代數(shù)圖形。那么使用顯卡進(jìn)行通用計(jì)算到底有多少優(yōu)勢(shì)呢？我們本次就來比較一番。首先我們?cè)贑PU上也實(shí)現(xiàn)該算法。為了方便起見我們?cè)O(shè)計(jì)了一個(gè)類：

class CPUCalc
{
private:
    int m_stride;
    int m_width;
    int m_height;
    float m_realMin;
    float m_imagMin;
    float m_scaleReal;
    float m_scaleImag;
    unsigned char* m_pData;

    void CalculatePoint(unsigned int x, unsigned int y);

public:
    CPUCalc(int stride, int width, int height, float rmin, float rmax, float imin, float imax, unsigned char* pData) 
        : m_stride(stride), m_width(width), m_height(height), m_realMin(rmin), m_imagMin(imin), m_scaleReal(0), m_scaleImag(0), m_pData(pData) 
    {
        m_scaleReal = (rmax - rmin) / width;
        m_scaleImag = (imax - imin) / height;
    }

    void Calculate();

};

在HLSL代碼中放在Constant Buffer中的數(shù)據(jù)，現(xiàn)在放在了類的成員處。注意我們這個(gè)類可以計(jì)算自定義復(fù)平面區(qū)間的曼德勃羅特集。rmin和rmax表示復(fù)平面的實(shí)數(shù)軸范圍而imin、imax則表示復(fù)平面的虛數(shù)軸范圍。這些參數(shù)的意義和上次HLSL中用的參數(shù)是一樣的，如果想自己實(shí)現(xiàn)該程序可以參考一下。

下面則是類的實(shí)現(xiàn)。我們用幾乎和HLSL一模一樣的做法來實(shí)現(xiàn)。其中某些HLSL的內(nèi)置方法采用類似的C++實(shí)現(xiàn)代替：

#include <algorithm>
#include <math.h>
#include "CPUCalc.h"

using std::max;

typedef unsigned int uint;
const uint MAX_ITER = 4096;

struct float2
{
    float x;
    float y;
};

inline float smoothstep(const float minv, const float maxv, const float v)
{
    if (v < minv) 
        return 0.0f;
    else if (v > maxv)
        return 1.0f;
    else
        return (v - minv) / (maxv - minv);
}

inline uint ComposeColor(uint index)
{    
    if (index == MAX_ITER) return 0xff000000;

    uint red, green, blue;
    
    float phase = index * 3.0f / MAX_ITER;
    red = (uint)(max(0.0f, phase - 2.0f) * 255.0f);
    green = (uint)(smoothstep(0.0f, 1.0f, phase - 1.3f) * 255.0f);
    blue = (uint)(max(0.0f, 1.0f - abs(phase - 1.0f)) * 255.0f);
    
    return 0xff000000 | (red << 16) | (green << 8) | blue;
}

void CPUCalc::CalculatePoint(uint x, uint y)
{
    float2 c;
    c.x = m_realMin + (x * m_scaleReal);
    c.y = m_imagMin + ((m_width - y) * m_scaleImag);
    
    float2 z;
    z.x = 0.0f;
    z.y = 0.0f;
    
    float temp, lengthSqr;
    uint count = 0;
    
    do
    {
        temp = z.x * z.x - z.y * z.y + c.x;
        z.y = 2 * z.x * z.y + c.y;
        z.x = temp;
        
        lengthSqr = z.x * z.x + z.y * z.y;
        count++;
    }
    while ((lengthSqr < 4.0f) && (count < MAX_ITER));    
    
    //write to result
    uint currentIndex = x * 4 + y * m_stride;
    uint& pPoint = *reinterpret_cast<uint*>(m_pData + currentIndex);
    
    pPoint = ComposeColor(static_cast<uint>(log((float)count) / log((float)MAX_ITER) * MAX_ITER));
}

void CPUCalc::Calculate()
{
    #pragma omp parallel for
    for (int y = 0; y < m_height; y++)
        for (int x = 0; x < m_width; x++)
        {
            CalculatePoint(x, y);
        }
}

最后我們?cè)黾恿艘粋€(gè)驅(qū)動(dòng)運(yùn)算的程序：Calculate()成員函數(shù)。它的實(shí)現(xiàn)中采用了OpenMP指令（#pragma omp）。OpenMP是C++的一個(gè)擴(kuò)展，用于實(shí)現(xiàn)統(tǒng)一地址空間的并行算法。這里采用的是一種靜態(tài)任務(wù)分配的做法，將for轉(zhuǎn)化為多個(gè)線程并行執(zhí)行。這個(gè)方法對(duì)曼德勃羅特集來說有一個(gè)弊端，一會(huì)我們?cè)僭敿?xì)討論。

影響這個(gè)程序性能的主要有三點(diǎn)：1、輸出像素的尺寸（參數(shù)中的width和height控制）；2、最大迭代次數(shù)（常數(shù)MAX_ITER控制）；3、所選的復(fù)平面區(qū)域（參數(shù)中的rmin、rmax、imin、imax控制）。因?yàn)閺?fù)平面中各個(gè)點(diǎn)的迭代次數(shù)都不一樣，所以無法確定算法的復(fù)雜度。按不超過最大迭代數(shù)來記，是一個(gè)系數(shù)很大的O(N)算法。我們這次測(cè)試固定所選的復(fù)平面區(qū)間為實(shí)數(shù)軸[-1.101，-1.099]以及虛數(shù)軸[2.229i，2.231i]的范圍。它的圖形是上篇中最后演示迭代次數(shù)那一組圖。該區(qū)間有相當(dāng)大的運(yùn)算量。然后我們分別固定最大迭代數(shù)和輸出像素?cái)?shù)，并變動(dòng)另外一個(gè)參數(shù)進(jìn)行多次測(cè)量，比較CPU和GPU進(jìn)行運(yùn)算的性能。

這次測(cè)試采用的CPU是Intel Core i7 920，具有四個(gè)核心，默認(rèn)主頻2.66GHz，搭配6GB DDR3-1333內(nèi)存。它還具有超線程技術(shù)，可以同時(shí)運(yùn)行8個(gè)線程。GPU是AMD Ati HD5850顯卡，默認(rèn)頻率725MHz（本次超頻775MHz）搭配1GB 1250MHz GDDR5顯存。該顯卡新片具有18組SIMD處理器，共有1440個(gè)Stream Core運(yùn)算單元。

首先我們固定最大迭代次數(shù)為512次，然后依次輸出像素512x512、1024x1024、2048x2048、4096x4096、8192x8192、16384x16384像素的圖片。下面是結(jié)果（時(shí)間單位為毫秒）：

前五項(xiàng)數(shù)據(jù)的圖表：

我們可以看到，GPU具有非常巨大的性能優(yōu)勢(shì)。即使是8個(gè)線程同時(shí)運(yùn)轉(zhuǎn)的酷睿i7處理器也完全不能同GPU相比。還可以觀察到一些事實(shí)：在像素增大的時(shí)候GPU和CPU以類似的速度變慢。GPU:CPU的速度比在很大范圍內(nèi)都保持在7.6倍左右。而當(dāng)像素?cái)?shù)增大到16384x16384之后GPU的性能突然下降了，速度比降低到了4.24。這是因?yàn)樵谌绱讼袼財(cái)?shù)下，顯卡的顯存已經(jīng)不夠用了，CPU分配了內(nèi)存給顯卡做虛擬顯存使用。這個(gè)過程消耗了不少時(shí)間，導(dǎo)致顯卡的性能大大受損。在處理較大數(shù)據(jù)的時(shí)候顯存容量顯得非常重要。因此nVidia和AMD都推出了通用計(jì)算專用版顯卡（Tesla和FireStream），這些專用計(jì)算卡的一大特點(diǎn)就是具有更大的顯存。

接下來我們將輸出像素?cái)?shù)固定在4096x4096上，然后測(cè)試不同的最大迭代數(shù)。從512一直增大到16384。下面是測(cè)試結(jié)果（時(shí)間單位為毫秒）：

前五項(xiàng)數(shù)據(jù)的圖表：

這次我們看到GPU出奇的優(yōu)勢(shì)。在迭代數(shù)增加到16384時(shí)，GPU竟然比CPU快出37倍之多！為什么會(huì)這樣呢？除了GPU本身并行計(jì)算確實(shí)強(qiáng)勁之外，我們的CPU算法也有一個(gè)問題。在曼德勃羅特集的運(yùn)算當(dāng)中，不是每個(gè)點(diǎn)都能達(dá)到最大迭代數(shù)。相當(dāng)多的點(diǎn)在不到最大迭代數(shù)之前就已經(jīng)計(jì)算完了。如果我們將復(fù)平面的點(diǎn)均勻分給多個(gè)線程的話，那就會(huì)有一些線程先計(jì)算完成，有一些線程后計(jì)算完成的問題。如果我們觀察運(yùn)算過程中的CPU占用率就會(huì)發(fā)現(xiàn)，差不多一半的時(shí)間里CPU都不能夠100%地充分利用：

當(dāng)?shù)鷶?shù)增大的時(shí)候，這種現(xiàn)象就更加明顯了。而GPU中我們的任務(wù)非常細(xì)（每個(gè)線程僅負(fù)責(zé)一個(gè)坐標(biāo)點(diǎn)的計(jì)算），DirectX實(shí)現(xiàn)了動(dòng)態(tài)線程分配，使得GPU中一旦有空閑的運(yùn)算單元就可以分配新的線程進(jìn)行運(yùn)算。因此我們CPU程序吃了很大的虧。不過，若我們重新編寫CPU算法，也采用動(dòng)態(tài)分配的話，就可以顯著提高CPU利用率，縮小這一差距。這個(gè)任務(wù)就留給讀者來完成了，試試看你能否寫出讓CPU盡可能保持100%的曼德勃羅特集算法。

原文：http://www.cnblogs.com/Ninputer/archive/2009/11/25/1610079.html