CUDA implementation of a totalistic cellular automaton. http://imgur.com/Yp6YwPU 50000 steps

makefile

all:
	nvcc -o celular totalistic-celular-automaton.cu -lX11

clean:
	rm ./celular

totalistic-celular-automaton.cu

#include <X11/Xlib.h> // must precede most other headers!
#include <stdlib.h>
#include <limits.h>
#include <time.h>
#include <sys/time.h>
#include <iostream>
#include <cstring>
#include <stdio.h>

void update_screen();

#define THREADS 20
#define HEIGHT 50
#define WIDTH 100
#define ENV_ITEMS (WIDTH * HEIGHT * THREADS * THREADS)
#define ENV_SIZE (ENV_ITEMS * sizeof(double))

#define cudaCheckErrors(msg) \
  do { \
    cudaError_t __err = cudaGetLastError(); \
    if (__err != cudaSuccess) { \
      fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
        msg, cudaGetErrorString(__err), \
        __FILE__, __LINE__); \
      fprintf(stderr, "*** FAILED - ABORTING\n"); \
      exit(1); \
    } \
  } while (0)

// window
Display *dsp;
Window win;
Atom wmDelete;
XEvent evt;
GC gc;

// env

__global__ void gpuNext(double* env, int row) {
  // int x = threadIdx.x;
  // int y = threadIdx.y + row;
  // env[y * WIDTH + x] = 1;

  int height = gridDim.y * blockDim.y;
  int width = gridDim.x * blockDim.x;

  int x = blockDim.x * blockIdx.x + threadIdx.x;
  int y = blockIdx.y * blockDim.y + threadIdx.y;

  if (y < row || y > row) return;

  int idx = gridDim.x * blockDim.x * (blockIdx.y * blockDim.y + threadIdx.y) +
    blockDim.x * blockIdx.x + threadIdx.x;


  int wrapNorth = ((height + y - 1) % height) * width;
  int wrapSouth = ((height + y + 1) % height) * width;

  int wrapEast = (width + x + 1) % width;
  int wrapWest = (width + x - 1) % width;

  double eastMiddle = env[y * width + wrapEast];
  double westMiddle = env[y * width + wrapWest];

  double eastNorth = env[wrapNorth + wrapEast];
  double westNorth = env[wrapNorth + wrapWest];

  double eastSouth = env[wrapSouth + wrapEast];
  double westSouth = env[wrapSouth + wrapWest];

  double middleNorth = env[wrapNorth + x];
  double middleSouth = env[wrapSouth + x];

  __syncthreads();

  double average = (westNorth + middleNorth + eastNorth + 0.4 * 3) / 3;
  env[idx] = average - (long int) average;

  // if (westNorth == 6 && middleNorth == 6 && eastNorth == 6)
  //   env[idx] = 3;
  //
  // if (westNorth == 3 && middleNorth == 3 && eastNorth == 3)
  //   env[idx] = 3;
  //
  // if (westNorth == 2 && middleNorth == 2 && eastNorth == 2)
  //   env[idx] = 2;

  // ###
  // if (westNorth && middleNorth && eastNorth)
  //   env[idx] = 0;
  //
  // // ##X
  // if (westNorth && middleNorth && !eastNorth)
  //   env[idx] = 1;
  //
  // // #X#
  // if (westNorth && !middleNorth && eastNorth)
  //   env[idx] = 0;
  //
  // // #XX
  // if (westNorth && !middleNorth && !eastNorth)
  //   env[idx] = 1;
  //
  // // X##
  // if (!westNorth && middleNorth && eastNorth)
  //   env[idx] = 1;
  //
  // // X##
  // if (!westNorth && middleNorth && !eastNorth)
  //   env[idx] = 0;
  //
  // // XX#
  // if (!westNorth && !middleNorth && eastNorth)
  //   env[idx] = 1;
  //
  // // XXX
  // if (!westNorth && !middleNorth && !eastNorth)
  //   env[idx] = 0;
}

unsigned int createRGB(int r, int g, int b) {
  return ((r & 0xff) << 16) + ((g & 0xff) << 8) + (b & 0xff);
}

void drawEnv(double* env) {
  XClearWindow(dsp,win);

  for(int i = 0; i < ENV_ITEMS; i++) {
    XSetForeground(dsp, gc, 0x000000);
    if (!env[i]) continue;

    // int weight = (int) (255 + 255 * env[i]) % 255;
    // unsigned int color = createRGB(weight, weight, weight);
    unsigned int color = INT_MAX * (1 - env[i]);
    XSetForeground(dsp, gc, color);

    XDrawPoint(dsp, win, gc, i % (WIDTH * THREADS), (int) i / (WIDTH * THREADS));
  }
}

__global__ void clearEnv(double* env) {
  env[gridDim.x * blockDim.x * (blockIdx.y * blockDim.y + threadIdx.y) +
    blockDim.x * blockIdx.x + threadIdx.x] = 0;

  env[(int) ((blockDim.x * gridDim.x) / 2)] = 1;
}

void initWindow() {
  dsp = XOpenDisplay(NULL);

  if (!dsp) {
    return;
  }

  int screen = DefaultScreen(dsp);
  unsigned int white = WhitePixel(dsp, screen);
  unsigned int black = BlackPixel(dsp, screen);

  win = XCreateSimpleWindow(dsp,
    DefaultRootWindow(dsp),
    0, 0, // origin
    THREADS * WIDTH, THREADS * HEIGHT, // size
    0, black, // border width/clr
    white); // backgrd clr

  wmDelete = XInternAtom(dsp, "WM_DELETE_WINDOW", True);
  XSetWMProtocols(dsp, win, &wmDelete, 1);

  gc = XCreateGC(dsp, win,
    0, // mask of values
    NULL); // array of values

  XSetForeground(dsp, gc, black);

  long eventMask = StructureNotifyMask;
  eventMask |= ExposureMask | ButtonPressMask | ButtonReleaseMask | KeyPressMask | KeyReleaseMask;
  XSelectInput(dsp, win, eventMask);

  KeyCode keyQ = XKeysymToKeycode(dsp, XStringToKeysym("Q"));

  XMapWindow(dsp, win);

  // wait until window appears
  // while (evt.type != MapNotify) {
  //   XNextEvent(dsp, &evt);
  // }
}

void sleepcp(int milliseconds) {
  clock_t time_end;
  time_end = clock() + milliseconds * CLOCKS_PER_SEC / 1000;
  while (clock() < time_end) {}
}

int main(){
  initWindow();

  dim3 blockDim(THREADS, THREADS, 1);
  dim3 gridDim(WIDTH, HEIGHT, 1);

  double* env = (double*) malloc(ENV_SIZE);
  double* dEnv;

  cudaMalloc((void**) &dEnv, ENV_SIZE);
  clearEnv<<<gridDim, blockDim>>>(dEnv);

  srand(time(0));

  bool draw = true;
  int i = 0;
  while (draw) {
    // XNextEvent(dsp, &evt);
    // XDrawRectangle(dsp, win, WhitePixel(dsp, gc), 1, 1, 497, 497);
    // XDrawRectangle(dsp, win, WhitePixel(dsp, gc), 50, 50, 398, 398);

    for (int j = 0; j < 200; j++) {
      gpuNext<<<gridDim, blockDim>>>(dEnv, (1 + i * 200 + j) % (HEIGHT * THREADS));
    }
    // gpuNext<<<gridDim, blockDim>>>(dEnv, (1 + i) % (HEIGHT * THREADS));

    cudaMemcpy(env, dEnv, ENV_SIZE, cudaMemcpyDeviceToHost);
    // std::cout << env;
    cudaCheckErrors("cudamalloc fail");

    drawEnv(env);
    // cudaCheckErrors("cudamalloc fail");

    i++;
    // sleepcp(1000);
  }

  // for (int j = 0; j < 50000; j++) {
  //   gpuNext<<<gridDim, blockDim>>>(dEnv, (1 + i * 50000 + j) % (HEIGHT * THREADS));
  // }
  // // gpuNext<<<gridDim, blockDim>>>(dEnv, (1 + i) % (HEIGHT * THREADS));
  //
  // cudaMemcpy(env, dEnv, ENV_SIZE, cudaMemcpyDeviceToHost);
  // // std::cout << env;
  // cudaCheckErrors("cudamalloc fail");
  //
  // drawEnv(env);
  //
  // while (draw) {}

  // XDestroyWindow(dsp, win);
  // XCloseDisplay(dsp);

  return 0;
}

// void update_screen()
// {
//   XClearWindow(dsp, win);
//
//   long i;
//   for (i = 0; i < XRES; i++) {
//     // XDrawPoint(dsp, win, gc, rand() % XRES, rand() % YRES);
//     pts[i].x = rand() % XRES;
//     pts[i].y = rand() % YRES;
//   }
//
//
//
//   return;
// }