之前我写了一篇用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开始采用以下公式

 

clip_image002

 

进行迭代。不同的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方法的调用。所得结果(缩小图)为:

test

 

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

image

 

我已经将所有源代码上传到了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和其他有趣的技术~

 posted on 2012-01-03 15:02  装配脑袋  阅读(11048)  评论(11编辑  收藏  举报