1
# I found this example for PyCuda here:
2
# http://wiki.tiker.net/PyCuda/Examples/Mandelbrot
4
# I adapted it for PyOpenCL. Hopefully it is useful to someone.
5
# July 2010, HolgerRapp@gmx.net
7
# Original readme below these lines.
9
# Mandelbrot calculate using GPU, Serial numpy and faster numpy
10
# Use to show the speed difference between CPU and GPU calculations
11
# ian@ianozsvald.com March 2010
13
# Based on vegaseat's TKinter/numpy example code from 2006
14
# http://www.daniweb.com/code/snippet216851.html#
15
# with minor changes to move to numpy from the obsolete Numeric
21
import numpy.linalg as la
25
# You can choose a calculation routine below (calc_fractal), uncomment
26
# one of the three lines to test the three variations
27
# Speed notes are listed in the same place
29
# set width and height of window, more pixels take longer to calculate
33
def calc_fractal_opencl(q, maxiter):
34
ctx = cl.Context(cl.get_platforms()[0].get_devices())
35
queue = cl.CommandQueue(ctx)
37
output = np.empty(q.shape, dtype=np.uint64)# resize(np.array(0,), q.shape)
40
q_opencl = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
41
output_opencl = cl.Buffer(ctx, mf.WRITE_ONLY, output.nbytes)
43
prg = cl.Program(ctx, """
44
__kernel void mandelbrot(__global float2 *q,
45
__global long *output, long const maxiter)
47
int gid = get_global_id(0);
48
float nreal, real = 0;
50
for(int curiter = 0; curiter < maxiter; curiter++) {
51
nreal = real*real - imag*imag + q[gid][0];
52
imag = 2* real*imag + q[gid][1];
55
if (real*real + imag*imag > 4.) {
56
output[gid] = curiter;
63
prg.mandelbrot(queue, output.shape, None, q_opencl,
64
output_opencl, np.int32(maxiter))
66
cl.enqueue_read_buffer(queue, output_opencl, output).wait()
72
def calc_fractal_serial(q, maxiter):
73
# calculate z using numpy
74
# this routine unrolls calc_fractal_numpy as an intermediate
75
# step to the creation of calc_fractal_opencl
76
# it runs slower than calc_fractal_numpy
77
z = np.zeros(q.shape, np.complex64)
78
output = np.resize(np.array(0,), q.shape)
79
for i in range(len(q)):
80
for iter in range(maxiter):
81
z[i] = z[i]*z[i] + q[i]
88
def calc_fractal_numpy(q, maxiter):
89
# calculate z using numpy, this is the original
90
# routine from vegaseat's URL
91
output = np.resize(np.array(0,), q.shape)
92
z = np.zeros(q.shape, np.complex64)
94
for iter in range(maxiter):
96
done = np.greater(abs(z), 2.0)
97
q = np.where(done,0+0j, q)
98
z = np.where(done,0+0j, z)
99
output = np.where(done, iter, output)
102
# choose your calculation routine here by uncommenting one of the options
103
calc_fractal = calc_fractal_opencl
104
# calc_fractal = calc_fractal_serial
105
# calc_fractal = calc_fractal_numpy
107
if __name__ == '__main__':
113
class Mandelbrot(object):
117
self.root.title("Mandelbrot Set")
124
def draw(self, x1, x2, y1, y2, maxiter=300):
125
# draw the Mandelbrot set, from numpy example
126
xx = np.arange(x1, x2, (x2-x1)/w*2)
127
yy = np.arange(y2, y1, (y1-y2)/h*2) * 1j
128
q = np.ravel(xx+yy[:, np.newaxis]).astype(np.complex64)
130
start_main = time.time()
131
output = calc_fractal(q, maxiter)
132
end_main = time.time()
134
secs = end_main - start_main
135
print "Main took", secs
137
output = (output + (256*output) + (256**2)*output) * 8
138
# convert output to a string
139
self.mandel = output.tostring()
141
def create_image(self):
143
create the image from the draw() string
145
self.im = Image.new("RGB", (w/2, h/2))
146
# you can experiment with these x and y ranges
147
self.draw(-2.13, 0.77, -1.3, 1.3)
148
self.im.fromstring(self.mandel, "raw", "RGBX", 0, -1)
150
def create_label(self):
151
# put the image on a label widget
152
self.image = ImageTk.PhotoImage(self.im)
153
self.label = tk.Label(self.root, image=self.image)