用GPU通用并行计算绘制曼德勃罗特集图形 下篇

上篇中我们用DirectX Compute Shader在显卡上编写了一个并行算法来计算好看的曼德勃罗特集迭代数图形。那么使用显卡进行通用计算到底有多少优势呢？我们本次就来比较一番。首先我们在CPU上也实现该算法。为了方便起见我们设计了一个类：

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中的数据，现在放在了类的成员处。注意我们这个类可以计算自定义复平面区间的曼德勃罗特集。rmin和rmax表示复平面的实数轴范围而imin、imax则表示复平面的虚数轴范围。这些参数的意义和上次HLSL中用的参数是一样的，如果想自己实现该程序可以参考一下。

下面则是类的实现。我们用几乎和HLSL一模一样的做法来实现。其中某些HLSL的内置方法采用类似的C++实现代替：

#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);
        }
}

最后我们增加了一个驱动运算的程序：Calculate()成员函数。它的实现中采用了OpenMP指令（#pragma omp）。OpenMP是C++的一个扩展，用于实现统一地址空间的并行算法。这里采用的是一种静态任务分配的做法，将for转化为多个线程并行执行。这个方法对曼德勃罗特集来说有一个弊端，一会我们再详细讨论。

影响这个程序性能的主要有三点：1、输出像素的尺寸（参数中的width和height控制）；2、最大迭代次数（常数MAX_ITER控制）；3、所选的复平面区域（参数中的rmin、rmax、imin、imax控制）。因为复平面中各个点的迭代次数都不一样，所以无法确定算法的复杂度。按不超过最大迭代数来记，是一个系数很大的O(N)算法。我们这次测试固定所选的复平面区间为实数轴[-1.101，-1.099]以及虚数轴[2.229i，2.231i]的范围。它的图形是上篇中最后演示迭代次数那一组图。该区间有相当大的运算量。然后我们分别固定最大迭代数和输出像素数，并变动另外一个参数进行多次测量，比较CPU和GPU进行运算的性能。

这次测试采用的CPU是Intel Core i7 920，具有四个核心，默认主频2.66GHz，搭配6GB DDR3-1333内存。它还具有超线程技术，可以同时运行8个线程。GPU是AMD Ati HD5850显卡，默认频率725MHz（本次超频775MHz）搭配1GB 1250MHz GDDR5显存。该显卡新片具有18组SIMD处理器，共有1440个Stream Core运算单元。

首先我们固定最大迭代次数为512次，然后依次输出像素512x512、1024x1024、2048x2048、4096x4096、8192x8192、16384x16384像素的图片。下面是结果（时间单位为毫秒）：

输出像素	CPU成绩	GPU成绩	速度比（G:C）
512 x 512	213	23	9.26
1024 x 1024	635	83	7.65
2048 x 2048	2403	312	7.70
4096 x 4096	9279	1227	7.56
8192 x 8192	37287	4894	7.61
16384 x 16384	152015	35793	4.24

前五项数据的图表：

我们可以看到，GPU具有非常巨大的性能优势。即使是8个线程同时运转的酷睿i7处理器也完全不能同GPU相比。还可以观察到一些事实：在像素增大的时候GPU和CPU以类似的速度变慢。GPU:CPU的速度比在很大范围内都保持在7.6倍左右。而当像素数增大到16384x16384之后GPU的性能突然下降了，速度比降低到了4.24。这是因为在如此像素数下，显卡的显存已经不够用了，CPU分配了内存给显卡做虚拟显存使用。这个过程消耗了不少时间，导致显卡的性能大大受损。在处理较大数据的时候显存容量显得非常重要。因此nVidia和AMD都推出了通用计算专用版显卡（Tesla和FireStream），这些专用计算卡的一大特点就是具有更大的显存。

接下来我们将输出像素数固定在4096x4096上，然后测试不同的最大迭代数。从512一直增大到16384。下面是测试结果（时间单位为毫秒）：

最大迭代数	CPU成绩	GPU成绩	速度比（G:C）
512	9363	1247	7.51
1024	18042	1513	11.92
2048	35132	2058	17.07
4096	68398	2897	23.61
8192	135074	4347	31.07
16384	266547	7082	37.63

 

前五项数据的图表：

这次我们看到GPU出奇的优势。在迭代数增加到16384时，GPU竟然比CPU快出37倍之多！为什么会这样呢？除了GPU本身并行计算确实强劲之外，我们的CPU算法也有一个问题。在曼德勃罗特集的运算当中，不是每个点都能达到最大迭代数。相当多的点在不到最大迭代数之前就已经计算完了。如果我们将复平面的点均匀分给多个线程的话，那就会有一些线程先计算完成，有一些线程后计算完成的问题。如果我们观察运算过程中的CPU占用率就会发现，差不多一半的时间里CPU都不能够100%地充分利用：

当迭代数增大的时候，这种现象就更加明显了。而GPU中我们的任务非常细（每个线程仅负责一个坐标点的计算），DirectX实现了动态线程分配，使得GPU中一旦有空闲的运算单元就可以分配新的线程进行运算。因此我们CPU程序吃了很大的亏。不过，若我们重新编写CPU算法，也采用动态分配的话，就可以显著提高CPU利用率，缩小这一差距。这个任务就留给读者来完成了，试试看你能否写出让CPU尽可能保持100%的曼德勃罗特集算法。