C++ AMP 实战：绘制曼德勃罗特集图像
之前我写了，里面使用的技术是与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) xtypedef double fp_t;#else
#define _F(x) x##ftypedef 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 Cunsigned 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) /
fp_t scale_imag = (imag_max - imag_min) /
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_
fp_t cy = imag_min + (height - gy) * scale_
fp_t zx = _F(0.0);
fp_t zy = _F(0.0);
fp_t length_
unsigned int count = 0;
temp = zx * zx - zy * zy +
zy = 2 * zx * zy +
length_sqr = zx * zx + 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,
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) {
red = bright,
green = t,
green = bright,
green = bright,
green = q,
green = p,
red = bright,
green = p,
unsigned int ired, igreen,
ired = (unsigned int)(red * 255.0f);
igreen = (unsigned int)(green * 255.0f);
iblue = (unsigned int)(blue * 255.0f);
return 0xff000000 | (ired && 16) | (igreen && 8) |
注意，在这段代码中我们使用的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 gdiplusStartupI
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 pngC
GetEncoderClsid(L&image/png&, &pngClsid);
output.Save(L&test.png&, &pngClsid, nullptr);
GdiplusShutdown(gdiplusToken);
以上代码绘制了实部从-2到2，虚部从-2i到2i范围的曼德勃罗特集图形，绘制到一个边长1024的位图上，最大迭代次数512次。注意其中synchronize方法的调用。所得结果（缩小图）为：
曼德勃罗特集图形每一点放大都有无穷无尽美丽的花纹，为了能够任意探索图形的细节，我还编写了一个Windows界面的程序，使用Direct2D 将所绘图像展示出来，而且加入了拖拽和鼠标滚轮的交互。实际运行效果非常酷，因为GPU并行计算的强大性能，所有交互都是平滑实时进行的。就好像是一张可以无限放大的图片一样。美中不足的是当前C++ AMP的Windows 7实现只支持单精度浮点数，所以放大到一定倍数就会模糊了。等到Windows 8的稳定版本出现，大家才能体验到双精度浮点的美妙结果。
我已经将所有源代码上传到了Github上，地址为：
主要，要编译和运行此例子，需要以下软硬件环境：
1. Windows 7 （不支持XP和Vista）
2. 硬件支持DirectX 11所有特性的显卡。如Nvidia Geforce 400,500系列显卡，AMD Radeon HD00系列显卡。无法运行于不支持DirectX 11的显卡上。
3. 需要Visual Studio 11 Developer Preview。
4. 需要Windows SDK。
我还翻译了一个Build大会关于C++ AMP的视频，是快速了解C++ AMP的最佳途径。观看地址：
欢迎关注我的微博，一起讨论C++ AMP和其他有趣的技术~
文章出处：
关注，看好C++ AMP。
正在下载VS11，准备试试……
请问，如果学习C++AMP的话有什么推荐的资料吗？
大哥，没有这个呀—— amp_math.h'&！
Error 1 error C1083: Cannot open include file: 'amp_math.h': No such file or directory d:\amp\ampdemo\mandelbrotviewer\mandelbrot.h 2 1 MandelbrotViewer百度拇指医生
&&&普通咨询
您的网络环境存在异常，
请输入验证码
验证码输入错误，请重新输入晏紫东 蒋梓乐 双人站 微博地址 /dongleclub/
粉丝：1413
视频地址复制
Flash地址复制
Html地址复制
离线看更方便
用或其他应用扫描二维码
花椒直播 日 初一 晏紫东花椒直播 预告
广播电视节目制作经营许可证：（沪）字第1248号
网络文化经营许可证：沪网文[6号
信息网络传播视听节目许可证：0910417
互联网ICP备案：沪ICP备号-3
沪ICP证：沪B2-
违法不良信息举报邮箱：
违法不良信息举报电话： 转 3系统提示信息
该信息不存在或已被删除！百度拇指医生
&&&普通咨询
您的网络环境存在异常，
请输入验证码
验证码输入错误，请重新输入