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

 

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

  1. 解析输入函数f之后,分别生成∂f/∂x和∂f/∂y两个偏导数,然后将这三个二元函数转化为合法的OpenCL表达式。
  2. 用OpenCL实现牛顿迭代法。
  3. 将图像上的每一点分派到一个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