 1#include <stdio.h>
 2#include <stdint.h>
 3
 4#define MAX_ITER 100
 5#define DIM 12800
 6static uint32_t colors[MAX_ITER + 1];
 7static uint32_t data[DIM * DIM];
 8
 9__device__ uint32_t mandelbrot(double x, double y) {
10  double zr = 0, zi = 0, zrsqr = 0, zisqr = 0;
11  int i;
12
13  for (i = 0; i < MAX_ITER; i++) {
14    zi = zr * zi * 2 + y;
15    zr = zrsqr - zisqr + x;
16    zrsqr = zr * zr;
17    zisqr = zi * zi;
18    if (zrsqr + zisqr > 4.0) {
19      break; // SIMT threads diverges here!
20    }
21  }
22  
23  return i;
24}
25
26__global__ void mandelbrot_kernel(uint32_t *data, double xmin, double ymin, double step, uint32_t *colors) {
27  int pix_per_thread = DIM * DIM / (gridDim.x * blockDim.x);
28  int tId = blockDim.x * blockIdx.x + threadIdx.x;
29  int offset = pix_per_thread * tId;
30  for (int i = offset; i < offset + pix_per_thread; i++) {
31    int x = i % DIM;
32    int y = i / DIM;
33    double cr = xmin + x * step;
34    double ci = ymin + y * step;
35    data[y * DIM + x] = colors[mandelbrot(cr, ci)];
36  }
37  if (gridDim.x * blockDim.x * pix_per_thread < DIM * DIM
38      && tId < (DIM * DIM) - (blockDim.x * gridDim.x)) {
39    int i = blockDim.x * gridDim.x * pix_per_thread + tId;
40    int x = i % DIM;
41    int y = i / DIM;
42    double cr = xmin + x * step;
43    double ci = ymin + y * step;
44    data[y * DIM + x] = colors[mandelbrot(cr, ci)];
45  }
46}
47
48int main() {
49  float freq = 6.3 / MAX_ITER;
50  for (int i = 0; i < MAX_ITER; i++) {
51    char r = sin(freq * i + 3) * 127 + 128;
52    char g = sin(freq * i + 5) * 127 + 128;
53    char b = sin(freq * i + 1) * 127 + 128;
54    colors[i] = b + 256 * g + 256 * 256 * r;
55  }
56  colors[MAX_ITER] = 0;
57
58  uint32_t *dev_colors, *dev_data;
59  cudaMalloc((void**)&dev_colors, sizeof(colors));
60  cudaMalloc(&dev_data, sizeof(data));
61  cudaMemcpy(dev_colors, colors, sizeof(colors), cudaMemcpyHostToDevice);
62
63  double xcen = -0.5, ycen = 0, scale = 3;
64  mandelbrot_kernel<<<512, 512>>>(
65    dev_data,
66    xcen - (scale / 2),
67    ycen - (scale / 2),
68    scale / DIM,
69    dev_colors
70  );
71
72  cudaMemcpy(data, dev_data, sizeof(data), cudaMemcpyDeviceToHost);
73  cudaFree(dev_data);
74  cudaFree(dev_colors);
75
76  FILE *fp = fopen("mandelbrot.ppm", "w");
77  fprintf(fp, "P6\n%d %d 255\n", DIM, DIM);
78  for (int i = 0; i < DIM * DIM; i++) {
79    fputc((data[i] >> 16) & 0xff, fp);
80    fputc((data[i] >>  8) & 0xff, fp);
81    fputc((data[i] >>  0) & 0xff, fp);
82  }
83  
84  return 0;
85}
11. 真实世界的并发编程¶

Take-away Messages¶

课后习题/编程作业¶

1. 阅读材料¶
