マンデルブロ集合をCUDAで描く（Mathematica） http://blog.unfindable.net/archives/3512

cuda-mandelbrot.nb

Needs["CUDALink`"]

CUDAQ[]

kernelCode = "
__device__ int check(Real_t x, Real_t y, int maxSteps) {
  Real_t re = 0, im = 0;
  int i = 0;
  while (i++ < maxSteps) {
    Real_t newRe = re * re - im * im + x;
    Real_t newIm = 2 * re * im + y;
    if (4 < newRe * newRe + newIm * newIm) return i;
    re = newRe;
    im = newIm;
  }
  return -1;
}

__global__ void mandelbrot(int* result, int maxSteps,
    int width, int height,
    Real_t xStart, Real_t xEnd, Real_t yStart, Real_t yEnd) {
  int xIndex = threadIdx.x + blockIdx.x * blockDim.x;
  int yIndex = threadIdx.y + blockIdx.y * blockDim.y;
  if (xIndex < width && yIndex < height) {
    Real_t x = xStart + (xEnd - xStart) * xIndex / width;
    Real_t y = yStart + (yEnd - yStart) * yIndex / height;
    result[xIndex + yIndex * width] = check(x, y, maxSteps);
  }
}";

kernel = CUDAFunctionLoad[kernelCode, "mandelbrot",
  {{"Integer32", _, "Output"}, "Integer32", "Integer32", "Integer32",
   _Real, _Real, _Real, _Real}, {16, 16}];

xMin = -2;
xMax = 1;
yMin = -3/2;
yMax = 3/2;
width = 300;
height = 300;

time = AbsoluteTiming[
    buffer = CUDAMemoryAllocate["Integer32", {height, width}]; 
    kernel[buffer, 100, width, height, N@xMin, N@xMax, N@yMin, N@yMax];
    result = CUDAMemoryGet[buffer];
    CUDAMemoryUnload[buffer];][[1]];
Print@time;
ArrayPlot[result]