Hello Mandelbrot

This is a Hello World program for CUDA. It calculates a Mandelbrot fractal and saves it in a PPM image file, which can be opened in most image viewers.

Since the only reason to run software on the GPU is to speed it up, the example includes a timer that prints the execution time. No CUDA program should be without error checking, so that is included as well.

To run, save the code to hello.cu and compile it with nvcc -o hello hello.cu. The program will run for less than a second. When it exits, there will be a new file in the same directory, mandelbrot.ppm, which can be opened in an image viewer to see the fractal, shown here on the left.

#include <cuda.h>
#include <fstream>
#include <stdio.h>
#include <cassert>

const int bailout = 1000;
const int w = 1000;
const int h = 1000;

typedef unsigned int u32;
typedef unsigned char u8;

using namespace std;

#define check_cuda_call(ans) { _check_cuda_call((ans), __FILE__, __LINE__); }
inline void _check_cuda_call(cudaError_t code, const char* file, int line)
{
  if (code != cudaSuccess) {
    fprintf(stderr,"CUDA Error: %s %s %d\n", cudaGetErrorString(code), file, line);
    exit(code);
  }
}

// EventTimer by Mark Harris.
class EventTimer {
public:
  EventTimer() : mStarted(false), mStopped(false) {
    cudaEventCreate(&mStart);
    cudaEventCreate(&mStop);
  }
  ~EventTimer() {
    cudaEventDestroy(mStart);
    cudaEventDestroy(mStop);
  }
  void start(cudaStream_t s = 0) {
    cudaEventRecord(mStart, s);
    mStarted = true;
    mStopped = false;
  }
  void stop(cudaStream_t s = 0)  {
    assert(mStarted);
    cudaEventRecord(mStop, s);
    mStarted = false;
    mStopped = true;
  }
  float elapsed() {
    assert(mStopped);
    if (!mStopped) return 0;
    cudaEventSynchronize(mStop);
    float elapsed = 0;
    cudaEventElapsedTime(&elapsed, mStart, mStop);
    return elapsed;
  }

private:
  bool mStarted, mStopped;
  cudaEvent_t mStart, mStop;
};


__global__ void MandelbrotKernel(int* out, float cr1, float cr2, float ci1, float ci2)
{
  u32 x(blockIdx.x * blockDim.x + threadIdx.x);
  u32 y(blockIdx.y * blockDim.y + threadIdx.y);

  if (x >= w || y >= h) {
    return;
  }

  float cr = (x / (float)w) * (cr2 - cr1) + cr1;
  float ci = (y / (float)h) * (ci2 - ci1) + ci1;

  float zi = 0.0f, zr = 0.0f, zr2 = 0.0f, zi2 = 0.0f, zit;
  u32 iter = bailout;
  while(--iter && zr2 + zi2 < 4.0f) {
    zit = zr * zi;
    zi = zit + zit + ci;
    zr = (zr2 - zi2) + cr;
    zr2 = zr * zr;
    zi2 = zi * zi;
  }

  if (iter) {
    iter = bailout - iter;
  }

  out[x + y * w] = iter * 5.0f;
}


// Translate (center + zoom) to (upper left + lower right)
void GetTranslatedCoordinates(float* cr1, float* cr2, float* ci1, float* ci2, float center_r, float center_i, float zoom) {
  *cr1 = center_r - zoom;
  *cr2 = center_r + zoom;
  float aspect_ratio = (float)w / (float)h;
  *ci1 = center_i - (zoom / aspect_ratio);
  *ci2 = center_i + (zoom / aspect_ratio);
}


void writePpm(int* mandel_bailouts)
{
  ofstream ofs("mandelbrot.ppm", ios::binary);
  ofs << "P6" << "\n" << w << " " << h << " " << 255 << "\n";
  for (u32 y = 0; y < h; ++y) {
    for (u32 x = 0; x < w; ++x) {
      int v = mandel_bailouts[x + (y * w)];
      if (v > 255) {
        v = 255;
      }
      u8 vb = static_cast<u8>(v);
      ofs.write((char*)&vb, sizeof(u8));
      ofs.write((char*)&vb, sizeof(u8));
      ofs.write((char*)&vb, sizeof(u8));
    }
  }
}


int div_up(int a, int b) {
  return ((a % b) != 0) ? (a / b + 1) : (a / b);
}


int main(int argc, char *argv[])
{
  float cr1, cr2, ci1, ci2;
  float center_r = -0.5f, center_i = 0.0f;
  float zoom = 1.5f;
  GetTranslatedCoordinates(&cr1, &cr2, &ci1, &ci2, center_r, center_i, zoom);

  int* mandel_escape_times_d;
  check_cuda_call(cudaMalloc(&mandel_escape_times_d, w * h * sizeof(int)));

  dim3 threads_per_block(16, 16); // 16 * 16 = 256
  dim3 blocks(div_up(w, threads_per_block.x), div_up(h, threads_per_block.y));

  EventTimer timer;
  timer.start();
  MandelbrotKernel<<<blocks, threads_per_block>>>(mandel_escape_times_d, cr1, cr2, ci1, ci2);
  timer.stop();

  // Check that kernel ran successfully.
  check_cuda_call(cudaDeviceSynchronize());
  check_cuda_call(cudaPeekAtLastError());

  printf("Elapsed time: %fms\n", timer.elapsed());

  int* mandel_escape_times_h;
  check_cuda_call(cudaMallocHost(&mandel_escape_times_h, w * h * sizeof(int)));
  check_cuda_call(cudaMemcpy(mandel_escape_times_h, mandel_escape_times_d, w * h * sizeof(int), cudaMemcpyDeviceToHost));
  writePpm(mandel_escape_times_h);

  check_cuda_call(cudaFree(mandel_escape_times_d));
  check_cuda_call(cudaFreeHost(mandel_escape_times_h));

  return 0;
}