This time, as I promised, we will see an example problem which can be solved on GPU with 10x speedup. This is a calculation of widely known Mandelbrot set.
I've added Mandelbrot sample to FsGPU project so you can get it there.
Mandelbrot sample
There is a fascinating Expert F# book which I highly recommend to read, especially if you are new to functional programming world. It contains a lot of real world samples, and Mandelbrot is one of them. Here is book web site where you can get source code for all the samples in the book. And here is Mandelbrot source code that I selected for acceleration on GPU.
Originally, Mandelbrot sample uses straightforward F# function to calculate whether given point belongs to mandelbrot set M or not
let sqrMod (x:complex) = x.r * x.r + x.i * x.i
let rec mandel maxit (z:Complex) (c: Complex) count =
if (sqrMod(z) < 4.0) && (count < maxit) then
mandel maxit ((z * z) + c) c (count + 1)
else count
> mandel 50 Complex.Zero Complex.One 0;;
val it : int = 2
> mandel 50 Complex.Zero Complex.OneI 0;;
val it : int = 50
This function calculates the number of iterations before we move outside closed disk of radius 2 around the origin. If we can reach maxit iterations, then c belongs to set M. The less value mandel function returns the more quickly we move from the set M which is represented by different colors.The point is that this function is executed for each pixel on the screen. For 800x600 screen we have at least 800x600=480000 calculations and 800x600x50=24000000 when all screen pixels belong to M. If we want to implement real-time navigation and visualization we need to make this calculations as fast as possible.
Mandelbrot acceleration on GPU
If we look at mandel function above we notice that it can be executed for every pixel in parallel. This is ideal case for data parallel computations, so our implementation can scale up to thousands of independently executing parallel GPU threads.
That is how we need to change original sample:
- Introduce a new function CalcMandelbrot like this
let threadsX, threadsY = 32,16 (* only 512 threads per block allowed *)
(* GPU execution code for mandelbrot calculation.
NOTE: This method is not thread safe, so you need to call it from
appropriate thread. (see CudaExecutor class) *)let MandelbrotExecution(ctx, (width : int, height : int, crunch : int, xOff : float32, yOff: float32, scale : float32)) =
let exec = ctx.Execution
let num = width * height
let colors = Array.zeroCreate<int> num
exec.AddParameter("colors", colors, ParameterDirection.Out) |> ignore
exec.AddParameter("width", uint32(width)) |> ignore
exec.AddParameter("height", uint32(height)) |> ignore
exec.AddParameter("crunch", uint32(crunch)) |> ignore
exec.AddParameter("xOff", float32(xOff)) |> ignore
exec.AddParameter("yOff", float32(yOff)) |> ignore
exec.AddParameter("scale", float32(scale)) |> ignore
let blocksX, blocksY = iDivUp(width, threadsX), iDivUp(height, threadsY)
exec.Launch(blocksX, blocksY, threadsX, threadsY, 1) |> ignore exec.ReadData(colors, 0)
colors
(* Calculate mandelbrot set for specified parameters using GPU kernel execution.
NOTE: this is thread-safe invoker of corresponding GPU execution *)
let CalcMandelbrot =
let executor = new CudaExecutor<_, int[]>(0, "Mandelbrot_sm10.cubin", "Mandelbrot0_sm10", MandelbrotExecution)
executor.GetInvoker()
- Implement this function using GPU kernel which I borrowed from CUDA SDK
template<class T>
__device__ inline int CalcMandelbrot(const T xPos, const T yPos, const int crunch)
{
T y = yPos;
T x = xPos;
T yy = y * y;
T xx = x * x;
int i = crunch;
while (--i && (xx + yy < T(4.0))) {
y = x * y * T(2.0) + yPos;
x = xx - yy + xPos;
yy = y * y;
xx = x * x;
}
return i;
}
// The Mandelbrot CUDA GPU kernel function
extern "C" __global__
void Mandelbrot0_sm10(int *dst, const int imageW, const int imageH, const int crunch,
const float xOff, const float yOff, const float scale)
{
const int ix = blockDim.x * blockIdx.x + threadIdx.x;
const int iy = blockDim.y * blockIdx.y + threadIdx.y;
if ((ix < < imageW) && (iy < imageH)) // Calculate the location
const float xPos = xOff + (float)ix * scale;
const float yPos = yOff - (float)iy * scale;
// Calculate the Mandelbrot index for the current location
int m = CalcMandelbrot<float>(xPos, yPos, crunch);
m = m > 0 ? crunch - m : crunch;
// Output the pixel
int pixel = imageW * iy + ix;
dst[pixel] = m;
}
}
- Integrate this function into UI menu and use it to render Bitmap object on the screen
let calcColors (bw:int) (bh:int) maxit map =
let tl : Complex = map 0 0 // top left conner
let bl = map 0 (bh-1) // bottom left conner
let xOff, yOff = float32(tl.r), float32(tl.i)
let scale = (float32(tl.i) - float32(bl.i)) / float32(bh)
let num = bw * bh
CalcMandelbrot(bw, bh, maxit, xOff, yOff, scale)
let linearFillGPU (bw:int) (bh:int) maxit map =
let colors = calcColors bw bh maxit map
let rect = new Rectangle(0, 0, bw, bh);
lock bitmap (fun () ->
let mutable bmpData = null
try
bmpData <- bitmap.LockBits(rect, System.Drawing.Imaging.ImageLockMode.ReadWrite, Imaging.PixelFormat.Format24bppRgb);
let ptr = bmpData.Scan0
let bytes = bmpData.Stride * bh
let rgbValues = Array.create bytes 0uy
Marshal.Copy(ptr, rgbValues, 0, bytes)
for y = 0 to bh - 1 do
for x = 0 to bw - 1 do
let index = (y * bw + x)
let offset = index * 3
let c = colors.[index]
let color = pickColor maxit c
rgbValues.[offset] <- color.B
rgbValues.[offset + 1] <- color.G
rgbValues.[offset + 2] <- color.R
Marshal.Copy(rgbValues, 0, ptr, bytes)
finally
bitmap.UnlockBits(bmpData)
form.Invalidate()
)
You can compare Mandelbrot.fs file in my sample with original version from Expert F# to see all the changes.
Threading problem
One significant problem with CUDA interop should be mentioned here as it is closely relate to implementation of CalcMandelbrot sample above.
Here is how it is described in CUDA Programming Guide:
"... Several host threads can execute device code on the same device, but by design, a
host thread can execute device code on only one device. As a consequence, multiple
host threads are required to execute device code on multiple devices. Also, any
CUDA resources created through the runtime in one host thread cannot be used by
the runtime from another host thread. ..."
host thread can execute device code on only one device. As a consequence, multiple
host threads are required to execute device code on multiple devices. Also, any
CUDA resources created through the runtime in one host thread cannot be used by
the runtime from another host thread. ..."
Thus if we need to call GPU kernels from multiple threads, we cannot do it directly, insteed we need to wrap actual GPU execution code with some thread safe wrapper. When we call wrapper from any our working thread it will provide that GPU code will be called from appropriate thread.
That is exactly what CudaExecutor helper type allow us to do.
Summing up
When you start Mandelbrot sample it runs two tests to compare GPU and CPU performance.
Here is the output on my hardware
Stress test mandelbrot calculation on GPU
20 times: 0 min, 0 sec, 335 ms
Stress test mandelbrot calculation on CPU
20 times: 0 min, 16 sec, 646 ms
Feel the difference!
Posts
