GPU并行计算版函数图像生成器（转）

http://www.cnblogs.com/Ninputer/archive/2011/08/18/2145045.html

前几天技术大牛Vczh同学开发了一个函数图像绘制程序，可以画出方程f(x,y)=0的图像。他的原理是用图像上每一点的坐标带入函数f得到针对x和y的两个方程，再用牛顿迭代法求解得到一组点集，然后画到图像上。用他的程序可以画出各种各样令人惊叹的方程图形。但是他的程序非常慢，因为对每一个点坐标都用牛顿迭代法求解是一项很费时的任务，即使采用了 Parallel.For，CPU算起来也很吃力。我研究了他的程序之后觉得可以用擅长并行计算的显卡来加速迭代法求解的过程。用OpenCL来完成这个任务再合适不过了。

整个过程还是相当顺利的，完全在Vczh原始程序的基础上改成。仅稍微改变了策略。步骤如下：

解析输入函数f之后，分别生成∂f/∂x和∂f/∂y两个偏导数，然后将这三个二元函数转化为合法的OpenCL表达式。
用OpenCL实现牛顿迭代法。
将图像上的每一点分派到一个OpenCL线程，然后由无数并行的OpenCL线程计算自己的点。

其中OpenCL代码如下：

fp_t func(fp_t x, fp_t y) 
{   
    return {动态生成}; 
}

fp_t df_dx(fp_t x, fp_t y)
{
    return {动态生成};
}

fp_t df_dy(fp_t x, fp_t y)
{
    return {动态生成};
}

fp_t solvex(fp_t start, const fp_t consty)
{
    for (int i = 0; i < MAX_ITER; ++i)
    {
        fp_t result = func(start, consty);
        
        if (result <= EPSILON && result >= -EPSILON)
        {
            return start;
        }

        fp_t d = df_dx(start, consty);
        if (d <= EPSILON && d >= -EPSILON)
        {
            return NAN;
        }
        else
        {
            start -= result / d;
        }
    }
    return NAN;
}

kernel void ComputeX(
    global write_only fp_t* points,
    int unit,
    int width,
    int cx,
    int cy,
    float origin_x,
    float origin_y)
{
    int gx = get_global_id(0);
    int gy = get_global_id(1);

    uint write_loc = gx + gy * width;
    
    fp_t py = origin_y + (fp_t)(gy + 1 - cy) / unit;
    fp_t px = origin_x + (fp_t)(cx - gx - 1) / unit;

    points[write_loc] = solvex(px, py);
}

这是求解f(x, a) = 0部分的代码，求解f(b, y) = 0的与之基本一样。其中fp_t是根据情况定义的typedef，可能是float或double。因为不是所有的OpenCL设备都支持双精度浮点，所以要写成通用类型的算法，用宏来控制。

在.NET中使用OpenCL，最容易的方法当然是使用Cloo库。Coo库完整地封装了OpenCL（1.1）的所有函数，并且是.NET非常容易使用的对象模型。我只使用了一次Cloo和OpenCL就再也不想忍受DirectX ComputeShader那麻烦的要死的类库。。

我的程序源代码已经完全上传到github。地址是：https://github.com/Ninputer/opencl-plot 点击Download即可打包下载所有代码。想直接运行的朋友可以点击此处下载二进制文件包。

要想运行这个程序，你必须安装了OpenCL的实现平台。目前Windows上的OpenCL主要有NVidia、AMD和Intel几家提供的实现平台。如果你有一块比较新的NVidia或AMD显卡，那么只需要装了最新的驱动包就自带了OpenCL。以下显卡支持双精度浮点：NVidia GeForce 200系列、400系列、500系列显卡；AMD Radeon HD 5800、5900、6900系列。其中Radeon 6900系列尚不支持官方双精度浮点数（cl_khr_fp64）扩展，所以本程序也支持cl_amd_fp64双精度浮点扩展，功能上完全一样。G80 和RV770等架构的显卡以及AMD中低端显卡只能支持单精度浮点数，绘制时的精度可能会略差。

没有支持OpenCL显卡的同学也可以用多核CPU来进行OpenCL计算，仍然要比原始的C#版本快一些。如果使用Intel Core i3、i5、i7系列CPU，可以使用Intel OpenCL SDK，下载地址：http://software.intel.com/en-us/articles/opencl-sdk/ 其他多核CPU都可以使用AMD APP SDK，下载地址：http://developer.amd.com/sdks/AMDAPPSDK/downloads/Pages/default.aspx

启动程序后允许先选择OpenCL计算平台和设备，如果安装了多个OpenCL平台可以任意选择。目前本程序暂时不支持多显卡并联技术（SLI、Crossfire）。NVidia CUDA平台界面示例：