gistfile1.cpp

#include <stdio.h>
#include <pgm.h>

#define MAXITER 255

#define CUDA_SAFE_CALL(E) do {                                          \
    cudaError_t e = (E);                                                \
    if (e != cudaSuccess) {                                             \
      printf("line %d: CUDA error: %s\n", __LINE__, cudaGetErrorString(e)); \
      exit(-2);                                                         \
    }                                                                   \
  } while (false)

static __global__ void
kernel(unsigned char *d_data, unsigned int height, unsigned int width)
{
  unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
  unsigned int j = blockIdx.y * blockDim.y + threadIdx.y;
  if (i >= width || j >= height) {
    return;
  }
  float2 c = make_float2(4.0f * i / width - 2.0f, 2.0f - 4.0f * j / height);
  float2 z = c;
  unsigned int iter;
  for (iter = 0; iter < MAXITER; ++iter) {
    if (z.x * z.x + z.y * z.y > 4.0f) {
      break;
    }
    z = make_float2(z.x * z.x - z.y * z.y + c.x, 2.0f * z.x * z.y + c.y);
  }
  d_data[j * width + i] = iter;
}

static void
writepgm(unsigned char *data, size_t height, size_t width)
{
  gray **array = pgm_allocarray(width, height);
  for (size_t y = 0; y < height; ++y) {
    for (size_t x = 0; x < width; ++x) {
      array[y][x] = data[y * width + x];
    }
  }
  FILE *file = fopen("mandelbrot.pgm", "w");
  pgm_writepgm(file, array, width, height, MAXITER, 0);
  fclose(file);
  pgm_freearray(array, height);
}

int
main(int argc, char **argv)
{
  pgm_init(&argc, argv);
  if (argc < 3) {
    printf("usage: %s height width\n", argv[0]);
    return -1;
  }
  size_t height = atoi(argv[1]);
  size_t width = atoi(argv[2]);
  printf("Height: %zu, Width: %zu\n", height, width);
  cudaEvent_t start, end;
  CUDA_SAFE_CALL(cudaEventCreate(&start));
  CUDA_SAFE_CALL(cudaEventCreate(&end));
  unsigned char *d_data;
  CUDA_SAFE_CALL(cudaMalloc((void **)&d_data,
                            height * width * sizeof(unsigned char)));
  CUDA_SAFE_CALL(cudaEventRecord(start));
  kernel<<<dim3((height + 15) / 16, (width + 15) / 16), dim3(16, 16)>>>
    (d_data, height, width);
  CUDA_SAFE_CALL(cudaEventRecord(end));
  CUDA_SAFE_CALL(cudaEventSynchronize(end));
  float ms;
  CUDA_SAFE_CALL(cudaEventElapsedTime(&ms, start, end));
  printf("%f milliseconds\n%f kilopixels/second\n", ms, height * width / ms);
  unsigned char *h_data = new unsigned char[height * width];
  CUDA_SAFE_CALL(cudaMemcpy(h_data, d_data,
                            height * width * sizeof(unsigned char),
                            cudaMemcpyDeviceToHost));
  CUDA_SAFE_CALL(cudaFree(d_data));
  writepgm(h_data, height, width);
  delete[] h_data;
  return 0;
}