【发布时间】：2017-04-08 19:43:17
【问题描述】：
我正在优化我在 pyOpenCL 中的 Mandelbrot 渲染器，并希望将迭代分成多个块，以便更好地利用我的 GPU。
最大迭代次数 = 1000 和 2 个“块”的示例：
1. 运行迭代 0-500 的 mandelbrot 逃逸算法。
2. 为迭代次数
第一个循环按预期工作，但之后的每个块都会导致错误的结果。我真的很想更具体一点，但我不知道真正的问题出在哪里（现在盯着代码超过 2 天）。
我怀疑从内核复制旧的 x,y（实部、虚部）部分时出了点问题，但我不知道如何调试它。
我在我的 GPU 和 CPU 上运行的结果相同，所以我猜它与 GPU 无关。
iterations=2000 和 10 个块的示例图像：
这几乎只是第一个块（加上一些“错误”像素）。
全部在一个块中完成（iterations=200 和 1 个块）：
iterations=2000 和 chunks = 1 的预期结果：
import pyopencl as cl
import numpy as np
from PIL import Image
from decimal import Decimal
def mandel(ctx, x, y, zoom, max_iter=1000, iter_steps=1, width=500, height=500, use_double=False):
mf = cl.mem_flags
cl_queue = cl.CommandQueue(ctx)
# build program
code = """
#if real_t == double
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#endif
kernel void mandel(
__global real_t *coords,
__global uint *output,
__global real_t *output_coord,
const uint max_iter,
const uint start_iter
){
uint id = get_global_id(0);
real_t2 my_coords = vload2(id, coords);
real_t x = my_coords.x;
real_t y = my_coords.y;
uint iter = 0;
for(iter=start_iter; iter<max_iter; ++iter){
if(x*x + y*y > 4.0f){
break;
}
real_t xtemp = x*x - y*y + my_coords.x;
y = 2*x*y + my_coords.y;
x = xtemp;
}
// copy the current x,y pair back
real_t2 val = (real_t2){x, y};
vstore2(val, id, output_coord);
output[id] = iter;
}
"""
_cltype, _nptype = ("double",np.float64) if use_double else ("float", np.float32)
prg = cl.Program(ctx, code).build("-cl-opt-disable -D real_t=%s -D real_t2=%s2" % (_cltype, _cltype))
# Calculate the "viewport".
x0 = x - ((Decimal(3) * zoom)/Decimal(2.))
y0 = y - ((Decimal(2) * zoom)/Decimal(2.))
x1 = x + ((Decimal(3) * zoom)/Decimal(2.))
y1 = y + ((Decimal(2) * zoom)/Decimal(2.))
# Create index map in x,y pairs
xx = np.arange(0, width, 1, dtype=np.uint32)
yy = np.arange(0, height, 1, dtype=np.uint32)
index_map = np.dstack(np.meshgrid(xx, yy))
# and local "coordinates" (real, imaginary parts)
coord_map = np.ndarray(index_map.shape, dtype=_nptype)
coord_map[:] = index_map
coord_map[:] *= (_nptype((x1-x0)/Decimal(width)), _nptype((y1-y0)/Decimal(height)))
coord_map[:] += (_nptype(x0), _nptype(y0))
coord_map = coord_map.flatten()
index_map = index_map.flatten().astype(dtype=np.uint32)
# Create input and output buffer
buffer_in_cl = cl.Buffer(ctx, mf.READ_ONLY, size=coord_map.nbytes)
buffer_out = np.zeros(width*height, dtype=np.uint32) # This will contain the iteration values of that run
buffer_out_cl = cl.Buffer(ctx, mf.WRITE_ONLY, size=buffer_out.nbytes)
buffer_out_coords = np.zeros(width*height*2, dtype=_nptype) # This the last x,y values
buffer_out_coords_cl = cl.Buffer(ctx, mf.WRITE_ONLY, size=buffer_out_coords.nbytes)
# 2D Buffer to collect the iterations needed per pixel
#iter_map = np.zeros(width*height, dtype=np.uint32).reshape((width, height)) #.reshape((height, width))
iter_map = np.zeros(width*height, dtype=np.uint32).reshape((height, width))
start_max_iter = 0
to_do = coord_map.size / 2
steps_size = int(max_iter / float(iter_steps))
while to_do > 0 and start_max_iter < max_iter:
end_max_iter = min(max_iter, start_max_iter + steps_size )
print "Iterations from iteration %i to %i for %i numbers" % (start_max_iter, end_max_iter, to_do)
# copy x/y pairs to device
cl.enqueue_copy(cl_queue, buffer_in_cl, coord_map[:to_do*2]).wait()
# and finally call the ocl function
prg.mandel(cl_queue, (to_do,), None,
buffer_in_cl,
buffer_out_cl,
buffer_out_coords_cl,
np.uint32(end_max_iter),
np.uint32(start_max_iter)
).wait()
# Copy the output back
cl.enqueue_copy(cl_queue, buffer_out_coords, buffer_out_coords_cl).wait()
cl.enqueue_copy(cl_queue, buffer_out, buffer_out_cl).wait()
# Get indices of "found" escapes
done = np.where(buffer_out[:to_do]<end_max_iter)[0]
# and write the iterations to the coresponding cell
index_reshaped = index_map[:to_do*2].reshape((to_do, 2))
tmp = index_reshaped[done]
iter_map[tmp[:,1], tmp[:,0]] = buffer_out[done]
#iter_map[tmp[:,0], tmp[:,1]] = buffer_out[done]
# Get the indices of non escapes
undone = np.where(buffer_out[:to_do]==end_max_iter)[0]
# and write them back to our "job" maps for the next loop
tmp = buffer_out_coords[:to_do*2].reshape((to_do, 2))
coord_map[:undone.size*2] = tmp[undone].flatten()
index_map[:undone.size*2] = index_reshaped[undone].flatten()
to_do = undone.size
start_max_iter = end_max_iter
print "%i done. %i unknown" % (done.size, undone.size)
# simple coloring by modulo 255 on the iter_map
return (iter_map % 255).astype(np.uint8).reshape((height, width))
if __name__ == '__main__':
ctx = cl.create_some_context(interactive=True)
img = mandel(ctx,
x=Decimal("-0.7546546453361122021732941811"),
y=Decimal("0.05020518634419688663435986387"),
zoom=Decimal("0.0002046859427855630601247281079"),
max_iter=2000,
iter_steps=1,
width=500,
height=400,
use_double=False
)
Image.fromarray(img).show()
编辑：Here 是实部/虚部永远不会离开 GPU 内存的另一个版本。
结果是一样的。
我完全没有想法。
【问题讨论】：
-
第三张图片，你说的对，我也觉得不对。它缺少您期望的细节（并且纵横比缩放不正确）。它是用 32 位浮点完成的吗（在 2000 次迭代时不够用）？在第一张图片下加上一些“错误”像素是否有线索？为什么会有错误的像素？
-
谢谢@WeatherVane。纵横比确实是错误的，这对这个测试应该没关系。我修复了缩放计算（仍然错误，但更好）。浮动图像：i.imgur.com/2rShDVl.png，双图像：i.imgur.com/SvwmhWg.png。相同（与另一个渲染器 guciek.github.io/… 类似（是的，我的渲染是垂直翻转的）。我会尽力为您提供一个更好的图像，说明计算出的错误像素。
-
但是恕我直言，真正的问题是我只是从内核中获取临时的 x,y（实数/虚数），如果该点没有逃脱（至少这是我尝试的）并再次提供它们去做）。所以我应该得到相同的图像，就像我在一次运行中计算出来的一样，还是我错了？
-
同意纵横比和图像翻转在这里不相关。我不使用 Python，但是
const uint start_iter是在哪里初始化的，或者设置为从之前的部分迭代中恢复？这似乎是一种非常奇怪的方式来进行并行处理的计算。在将每个点传递给下一个chunk之前，您是否检测到每个点是否已经逃逸？我原以为将一整行传递给每个线程会比恢复部分迭代更有效。当任何线程完成时，您可以将其放到下一行。 -
@WeatherVane 我在第 70 行和第 107 行设置/更新 start_iter。是的，我确实检测到了转义点，并且只继续使用非转义点。我的想法是尽量减少等待其他未逃脱的早期逃脱线程（据我所知，当使用工作组大小> 1时）。目标是减少计算高迭代点（200 万及以上）所需的时间。我在这里可能完全错了，但无法比较性能，因为我没有让它正常运行，）
标签： python numpy opencl mandelbrot pyopencl