MandelBrot

The Mandelbrot set is the set of complex numbers 'c' for which the function f(z) = z² + c doesn't diverge when iterated from z = 0.

This example illustrates a slightly more complex code, and presents usage of CUDA 2D grid in C# code.

We compute a rendered image on the square [-2.0,2.0] [-2.0, 2.0]. Since the program computes a static image (and is not interactive), parameters are stored in constant values:

  const int maxiter = 256; // iteration upper bound 
  const int N = 1024; // pixel count in rendered image
  const float fromX = -2.0f; // x lower bound
  const float fromY = -2.0f; // y lower bound
  const float size = 4.0f;  
  const float h = size / (float)N; // grid element size

The IterCount function is used to compute the number of iterations before sequence gets out the control disk. If we reach maxiter, it's assumed the point is part of the set.

  [Kernel]
  public static int IterCount(float cx, float cy)
  {
      int result = 0;
      float x = 0.0f;
      float y = 0.0f;
      float xx = 0.0f, yy = 0.0f;
      while (xx + yy <= 4.0f && result < maxiter)
      {
          xx = x * x;
          yy = y * y;
          float xtmp = xx - yy + cx;
          y = 2.0f * x * y + cy;
          x = xtmp;
          result++;
      }

      return result;
  }

The function is plain scalar C# code, decorated with [Kernel] to ensure it will generate a device function.

In order to get best performances, we need to perform computations on a 2D grid (to save registers and distribute work evenly among multiprocessors):

    [EntryPoint("run")]
    public static void Run(int[] light, int lineFrom, int lineTo)
    {
        for (int line = lineFrom + threadIdx.y + blockDim.y * blockIdx.y; line < lineTo; line += gridDim.y * blockDim.y)
        {
            for (int j = threadIdx.x + blockIdx.x * blockDim.x; j < N; j += blockDim.x * gridDim.x)
            {
                float x = fromX + line * h;
                float y = fromY + j * h;
                light[line * N + j] = IterCount(x, y);
            }
        }
    }

As usual, we wrap the parent type and invoke this entry point:

    HybRunner runner = HybRunner.Cuda("Man
    dynamic wrapper = runner.Wrap(new Program());
    wrapper.Run(light, 0, N);

The C# code (not accelerated), can still be called in a Parallel.For:

Parallel.For(0, N, (line) =>
{
    Run(light, line, line + 1);
});

Performance results of various versions can be found on our website.

Hello world

Home

Printf