C++ AMP实战：绘制曼德勃罗特集图像

之前我写了一篇用GPU绘制曼德勃罗特（Mandelbrot）集图像的文章，里面使用的技术是与DirectX 11继承在一起的DirectCompute。DirectCompute执行在GPU上的kernel代码，必须用一种特殊的HLSL语言来编写。虽然这种语言有些类似于C，但一些特殊的细节使得没接触过DirectX的开发人员很不适应。相比于kernel代码，驱动HLSL所要进行的准备工作那简直麻烦得要命，所以我在那篇博客里索性略去了。如果要想要体会一下DirectCompute那麻烦又不直观的API，可以参考我翻译的另一篇博客。后来我转向了OpenCL，这是一种异构并行计算的行业标准。它的的kernel语言和API相比DirectCompute那可是简单了好几倍，又有.NET封装类库，无疑是个方便的选择。但是经过我的实践，至少在我的AMD 5850显卡上，OpenCL的性能没有DirectCompute那么好。能否鱼和熊掌兼得呢？现在这个选择出现了：C++ AMP

C++ AMP全名C++ Accelerated Massive Parallelism（加速大规模并行计算）。是微软提出的基于C++的异构化并行计算平台。它将随Visual Studio 11一起发布，目前为预览版本。所谓异构并行计算，主要的需求就来自于GPU通用计算的崛起。GPU非常适合大规模数据并行算法，即同一程序应多多组不同的数据进行并行运算。然而GPU的架构与主流CPU不同，而且常常更新换代，没法采用传统编程语言来编程。现有的GPU多数编程方案，如DirectCompute和OpenCL，都要使用不同的语言或编译器来编写运行于GPU上的kernel部分和运行在CPU上的host部分。C++ AMP统一了这两部分，可以用同一个编译器，同一种语法来编写kernel代码；无需任何编译器选项或设置。更甚至C++ AMP的API简单到了极致，比OpenCL的方便程度更上了一个层次。下面我们就用C++ AMP来实现经典的曼德勃罗特集图形绘制。

先来回顾一下曼德勃罗特集，这个是很多并行计算书都用来做例子的经典问题，因为它是一个典型的易并行算法。选任意复数c，让z从0开始采用以下公式

进行迭代。不同的c会产生很不一样的效果。有的c可以一直迭代下去（z的模都在半径为2的圆以内），而有的c会导致迭代若干次之后逃逸（其模会变得大于2，进而趋近于无穷）。不同的c迭代次数相差很大，如果把平面上的每个点坐标都当作c来进行迭代，然后把迭代次数绘制成不同的颜色，就会形成绚丽的图形。而不同c可以完全独立地进行独立，不会影响其他c值的计算结果。因此算法可以是完全并行的。

首先我们要开辟出一块空间来表示最后生成的图像。图像上的每一个像素将会映射到一个复数c上，这样a x b的图像，就表示a x b个不同的复数，分别进行上述迭代。C++ AMP提供了两种表示数据的方式：array模板类和array_view模板类。array模板类就像C++ 11的STL容器array一样，是一个固定尺寸、固定维度的数组。和STL array唯一的不同是C++ AMP的array直接分配在GPU的显存中。声明array时要传入两个模板参数——array元素的类型和维度。比如我们要构建一个unsigned int类型的二维数组，就可以这样写：

int a = 100, b = 100;
array<unsigned int, 2> buffer(a, b);

注意这个array模板是定义在Concurrency命名空间里的，需要提前using。

有时我们需要把CPU上的数据传送到GPU上，或者希望GPU计算的结果直接写入指定的CPU内存中，这时我们就可以用第二个模板类array_view。array_view本身不会分配任何内存，它只是一个CPU数据buffer的包装，其中数据buffer可以是一个指针、数组，也可以是诸如std::vector等各种常用容器。array_view仅仅当数据访问发生时才进行必要的数据拷贝。声明array_view和声明array一样，但必须传入它的基础容器：

int a = 100, b = 100;
std::vector<unsigned int> myarray(a * b);

array_view<unsigned int, 2> view(a, b, myarray);

接下来最重要的是让代码在GPU上执行，完成这个任务的入口方法是parallel_for_each方法。它非常像.NET 4里面的Parallel.ForEach方法，是一个通过lambda函数来分派任务的方法。我们先来完成一个简单的任务，现在有两个array_view，要生成一个array，其内容是两个array_view中相应元素的和。

std::vector<unsigned int> arr1(10);
std::vector<unsigned int> arr2(10);

array_view<unsigned int, 1> view1(10, arr1);
array_view<unsigned int, 1> view2(10, arr2);

array<unsigned int, 1> r(10);

parallel_for_each(r.grid, [view1, view2, &r](index<1> i) restrict(direct3d)
{
    r[i] = view1[i] + view2[i];
});

观察这里的parallel_for_each，有许多有趣的信息。首先它的第一个参数是一个grid类型的对象。grid是C++ AMP用来描述并行计算拓扑结构的对象。因为经过parallel_for_each方法分派的计算是在显卡计算的，我们需要事先指定有多少GPU线程来完成计算，以及这些线程的排列方式。这个任务当中，我们想要把两个array_view表示的数据想加，为了达到最大并行性，我们给数组的每一个元素都分配一个线程。代码里的r.grid就表示线程的排列方式由数组r自身的拓扑结构决定——在这个例子中就是一个1 x 10的一维结构。注意，GPU上通常可以开成百上千的线程，几乎没有任何代价，所以线程的拓扑结构应当尽量达到最大的并行度。

接下来那个方括号语法是C++的lambda表达式语法。C++ 11开始支持lambda表达式，和其他语言的lambda表达式一样，这是一种就地声明匿名函数的语法。lambda表达式在C++ 11中具有极其重要的地位。一开始的方括号是lambda表达式的捕获列表。与C#等语言的lambda表达式不同，C++得lambda表达式需要显式进行变量捕获，只有列在捕获列表中的变量才能在lambda函数体内使用。注意，捕获的时候直接写名字的变量（如上面例子中的view1和view2）是按值捕获的变量，在lambda函数体内修改它们的值，不会改变原始变量。而名字前加一个&符号的捕获变量是按引用捕获的变量，lambda函数对其修改会反应到原始变量中。在C++ AMP的parallel_for_each调用时有一个规定：所有变量都必须按值捕获，除了array对象，array对象必须按引用捕获。这就是上面例子为何要这样写的原因。

捕获列表之后是lambda函数的参数，定义方法和普通函数的参数方法是一样的。parallel_for_each要求传入的函数接受一个index类型的特殊参数。这个参数是实际运行时的线程编号。因为线程的拓扑结构已经由grid参数决定了，所以这个index运行时就会是grid表示范围中的一个值，维度和grid一样。比如上面的例子，index就会表示0-9号线程中的某一个，这些线程都是并行执行的。

参数表之后，出现了一个新的修饰符restrict(direct3d)。这是专为C++ AMP设计的新修饰符，用来表示函数的语法约束。众所周知GPU和CPU是有很多不同的，在GPU上执行的方法无法进行某些操作。比如GPU代码无法调用Windows API，也无法递归。restrict(direct3d)告诉编译器检查函数体中的代码，并且在编译时就能检测出所有不满足约束的语法。restrict修饰符还是函数签名的一部分，可以针对不同的restrict修饰符进行函数重载，这和const等C++原有的修饰符是一样的。除了restrict(direct3d)，还有一种restrict(cpu)修饰符。所有现有C++函数都默认带有restrict(cpu)修饰符。一个函数还可以同时约束为CPU和GPU代码，只要这么写即可：restrict(cpu, direct3d)。在restrict(direct3d)函数中调用其它的函数，会自动选择restrict(direct3d)的重载。稍后我们会看到，C++ AMP定义了一组restrict(direct3d)版本的数学函数，如log和sin等。在restrict(direct3d)代码中调用这些数学函数就会自动使用GPU版本，而普通的CPU代码则会去调用不带修饰的CPU版本。

最后是lambda函数的函数体，里面的代码非常直白，就是直接对当前线程index所对应的数组元素进行相加。当你运行上述代码，该代码就会在GPU上以并行计算的方式完成。没有任何额外的初始化代码，隐藏代码或者编译器选项设置——只要这么写就可以运行了。是不是很方便？

下面我们就来完成GPU上的曼德勃罗特集生成算法。先看头文件mandelbrot.h，里面include了amp.h——C++ AMP的主要头文件。还声明了一个fp_t类型，根据宏来控制它表示float或double。目前Windows 7上的C++ AMP还不支持双精度浮点运算。

#include "amp.h"

//#define FP64

#if defined FP64
#define _F(x) x
typedef double fp_t;
#else
#define _F(x) x##f
typedef float fp_t;
#endif

void generate_mandelbrot(
    Concurrency::array_view<unsigned int, 2> result,
    unsigned int max_iter,
    fp_t real_min,
    fp_t imag_min,
    fp_t real_max,
    fp_t imag_max );

接下来是实现的代码，还顺便生成了一个HSB（色调、饱和度、亮度）到RGB转换的GPU函数：

#include "stdafx.h"
#include "mandelbrot.h"

using namespace Concurrency;

unsigned int set_hsb (float hue, float saturate, float bright) restrict (direct3d);

void generate_mandelbrot(
    array_view<unsigned int, 2> result,
    unsigned int max_iter,
    fp_t real_min,
    fp_t imag_min,
    fp_t real_max,
    fp_t imag_max )
{
    int width = result.extent.get_x();
    int height = result.extent.get_y();

    fp_t scale_real = (real_max - real_min) / width;
    fp_t scale_imag = (imag_max - imag_min) / height;

    parallel_for_each(result.grid, [=](index<2> i) restrict(direct3d)
    {
        int gx = i.get_x();
        int gy = i.get_y();        

        fp_t cx = real_min + gx * scale_real;
        fp_t cy = imag_min + (height - gy) * scale_imag;

        fp_t zx = _F(0.0);
        fp_t zy = _F(0.0);

        fp_t temp;
        fp_t length_sqr;

        unsigned int count = 0;
        do
        {
            count++;
        
            temp = zx * zx - zy * zy + cx;
            zy = 2 * zx * zy + cy;
            zx = temp;

            length_sqr = zx * zx + zy * zy;
        }
        while((length_sqr < _F(4.0)) && (count < max_iter));
    
        //faster using multiplication than division
        float n = count * 0.0078125f; // n = count / 128.0f; 
        
        float h = 1.0f - 2.0f * fabs(0.5f - n + floor(n));

        //turn points at maximum iteration to black
        float bfactor = direct3d::clamp((float)(max_iter - count), 0.0f, 1.0f);

        result[i] = set_hsb(h, 0.7f, (1.0f - h * h * 0.83f) * bfactor);
    });
}

unsigned int set_hsb (float hue, float saturate, float bright) restrict (direct3d)
{   

    float red, green, blue;      
    float h = (hue * 256) / 60;  
    float p = bright * (1 - saturate);  
    float q = bright * (1 - saturate * (h - (int)h));  
    float t = bright * (1 - saturate * (1 - (h - (int)h)));  
    
    switch ((int)h) {  
    case 0:   
        red = bright,  green = t,  blue = p;  
        break;  
    case 1:  
        red = q,  green = bright,  blue = p;  
        break;  
    case 2:  
        red = p,  green = bright,  blue = t;  
        break;  
    case 3:  
        red = p,  green = q,  blue = bright;  
        break;  
    case 4:  
        red = t,  green = p,  blue = bright;  
        break;  
    case 5:  
    case 6:  
        red = bright,  green = p,  blue = q;  
        break;  
    }  

    unsigned int ired, igreen, iblue;
    ired = (unsigned int)(red * 255.0f);
    igreen = (unsigned int)(green * 255.0f);
    iblue = (unsigned int)(blue * 255.0f);
    
    return 0xff000000 | (ired << 16) | (igreen << 8) | iblue;  
 }

注意，在这段代码中我们使用的array_view是二维的，所以index也相应是一个二维的结构。我们可以从index的get_x()，get_y()等方法得到线程的编号。代码中还使用了若干小技巧，比如判断复数模的平方大于4而不是模大于2。因为开平方运算即使在GPU上也是个比较慢的操作，需要尽量避免。方法的最后部分使用一个特殊的算法将迭代数转化成了颜色。这个公式是我反复试验得到的，只是为了视觉上好看。大家自己实践可以根据自己喜好随便改。

C++ AMP的parallel_for_each调用是一种“准同步”的调用，当parallel_for_each完成时，GPU计算并不一定结束了，但代码会继续执行。只有当你访问计算的结果——比如包含结果的array_view时，它就会阻塞线程等待GPU计算的完成。我们可以调用array_view的synchornize方法来显式地等待。C++ AMP也提供真正异步接口，但用起来比较麻烦，在这个例子中就不采用了。下面的代码就是调用以及同步的代码。这里展示代码将结果保存在一个位图中：

int _tmain(int argc, _TCHAR* argv[])
{

    GdiplusStartupInput gdiplusStartupInput;
    ULONG_PTR           gdiplusToken;
    GdiplusStartup(&gdiplusToken, &gdiplusStartupInput, nullptr);


    const int edge = 1024;

    std::vector<unsigned int> a(edge * edge);
    array_view<unsigned int, 2> av(edge, edge, a);

    generate_mandelbrot(av, 512, -2, -2, 2, 2);

    av.synchronize();

    //construct a image    
    {
        Bitmap output(edge, edge, edge * 4, PixelFormat32bppARGB, reinterpret_cast<BYTE*>(a.data()));

        CLSID pngClsid;
        GetEncoderClsid(L"image/png", &pngClsid);

        output.Save(L"test.png", &pngClsid, nullptr);
    }
    
    GdiplusShutdown(gdiplusToken);
    return 0;
}

以上代码绘制了实部从-2到2，虚部从-2i到2i范围的曼德勃罗特集图形，绘制到一个边长1024的位图上，最大迭代次数512次。注意其中synchronize方法的调用。所得结果（缩小图）为：

曼德勃罗特集图形每一点放大都有无穷无尽美丽的花纹，为了能够任意探索图形的细节，我还编写了一个Windows界面的程序，使用Direct2D将所绘图像展示出来，而且加入了拖拽和鼠标滚轮的交互。实际运行效果非常酷，因为GPU并行计算的强大性能，所有交互都是平滑实时进行的。就好像是一张可以无限放大的图片一样。美中不足的是当前C++ AMP的Windows 7实现只支持单精度浮点数，所以放大到一定倍数就会模糊了。等到Windows 8的稳定版本出现，大家才能体验到双精度浮点的美妙结果。

我已经将所有源代码上传到了Github上，地址为：https://github.com/Ninputer/AMP-Demo

主要，要编译和运行此例子，需要以下软硬件环境：

1. Windows 7 （不支持XP和Vista）

2. 硬件支持DirectX 11所有特性的显卡。如Nvidia Geforce 400,500系列显卡，AMD Radeon HD5000,6000,7000系列显卡。无法运行于不支持DirectX 11的显卡上。

3. 需要Visual Studio 11 Developer Preview。下载地址

4. 需要Windows SDK。下载地址

我还翻译了一个Build大会关于C++ AMP的视频，是快速了解C++ AMP的最佳途径。观看地址：http://www.tudou.com/playlist/p/l13690258i105735621.html

欢迎关注我的微博http://weibo.com/ninputer，一起讨论C++ AMP和其他有趣的技术~