mandelbrot.c 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512
  1. /* StarPU --- Runtime system for heterogeneous multicore architectures.
  2. *
  3. * Copyright (C) 2010-2021 Université de Bordeaux, CNRS (LaBRI UMR 5800), Inria
  4. *
  5. * StarPU is free software; you can redistribute it and/or modify
  6. * it under the terms of the GNU Lesser General Public License as published by
  7. * the Free Software Foundation; either version 2.1 of the License, or (at
  8. * your option) any later version.
  9. *
  10. * StarPU is distributed in the hope that it will be useful, but
  11. * WITHOUT ANY WARRANTY; without even the implied warranty of
  12. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  13. *
  14. * See the GNU Lesser General Public License in COPYING.LGPL for more details.
  15. */
  16. #include <stdio.h>
  17. #include <stdlib.h>
  18. #include <string.h>
  19. #include <unistd.h>
  20. /* Uncomment this to activate X11 display */
  21. //#define USE_X11
  22. #define SHORT_LOG 1
  23. #define ROUND_ROBIN
  24. #ifdef USE_X11
  25. #include <X11/Xlib.h>
  26. #include <X11/Xutil.h>
  27. int use_x11 = 1;
  28. #else
  29. int use_x11 = 0;
  30. #endif
  31. int demo = 0;
  32. int frames = -1;
  33. #include <pthread.h>
  34. #include <assert.h>
  35. #include <sys/time.h>
  36. #ifdef __APPLE_CC__
  37. #include <OpenCL/opencl.h>
  38. #else
  39. #include <CL/cl.h>
  40. #endif
  41. #define error(...) do { fprintf(stderr, "Error: " __VA_ARGS__); exit(EXIT_FAILURE); } while(0)
  42. #define check(err, str) do { if(err != CL_SUCCESS) { fprintf(stderr, "OpenCL Error (%d): %s\n",err, str); exit(EXIT_FAILURE); }} while(0)
  43. #ifdef UNUSED
  44. #elif defined(__GNUC__)
  45. # define UNUSED(x) UNUSED_ ## x __attribute__((unused))
  46. #else
  47. # define UNUSED(x) x
  48. #endif
  49. const char * kernel_src = "\
  50. #pragma OPENCL EXTENSION cl_khr_fp64 : enable\n\
  51. #define TYPE double \n\
  52. #define MIN(a,b) (((a)<(b))? (a) : (b))\n\
  53. __kernel void mandelbrot_kernel(__global uint * a,\n\
  54. TYPE leftX, TYPE topY,\n\
  55. TYPE stepX, TYPE stepY,\n\
  56. uint maxIt, uint iby, uint block_size)\n\
  57. {\n\
  58. TYPE xc = leftX + get_global_id(0) * stepX;\n\
  59. TYPE yc = iby*block_size*stepY + topY + get_global_id(1) * stepY;\n\
  60. int it;\n\
  61. TYPE x,y;\n\
  62. x = y = (TYPE)0.0;\n\
  63. for (it=0;it<maxIt;it++)\n\
  64. {\n\
  65. TYPE x2 = x*x;\n\
  66. TYPE y2 = y*y;\n\
  67. if (x2+y2 > (TYPE)4) break; \n\
  68. TYPE twoxy = (TYPE)2*x*y;\n\
  69. x = x2 - y2 + xc;\n\
  70. y = twoxy + yc;\n\
  71. }\n\
  72. uint v = MIN((1024*((float)(it)/(2000))), 256);\n\
  73. a[get_global_id(0) + get_global_id(1)*get_global_size(0)] = (v<<16|(255-v)<<8); \n\
  74. }";
  75. static cl_uint nblocks = 8;
  76. static cl_uint height = 768;
  77. static cl_uint width = 1024;
  78. static cl_uint maxIt = 20000;
  79. static cl_uint group_size = 64;
  80. static double leftX = -0.745;
  81. static double rightX = -0.74375;
  82. static double topY = .15;
  83. static double bottomY = .14875;
  84. #ifdef USE_X11
  85. /* X11 data */
  86. static Display *dpy;
  87. static Window win;
  88. static XImage *bitmap;
  89. static GC gc;
  90. static KeySym Left=-1, Right, Down, Up, Alt ;
  91. static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
  92. static void exit_x11(void)
  93. {
  94. XDestroyImage(bitmap);
  95. XDestroyWindow(dpy, win);
  96. XCloseDisplay(dpy);
  97. }
  98. static void init_x11(int width, int height, cl_uint *buffer)
  99. {
  100. /* Attempt to open the display */
  101. dpy = XOpenDisplay(NULL);
  102. /* Failure */
  103. if (!dpy)
  104. exit(0);
  105. unsigned long white = WhitePixel(dpy,DefaultScreen(dpy));
  106. unsigned long black = BlackPixel(dpy,DefaultScreen(dpy));
  107. win = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy), 0, 0,
  108. width, height, 0, black, white);
  109. /* We want to be notified when the window appears */
  110. XSelectInput(dpy, win, StructureNotifyMask);
  111. /* Make it appear */
  112. XMapWindow(dpy, win);
  113. XTextProperty tp;
  114. char name[128] = "Mandelbrot";
  115. char *n = name;
  116. Status st = XStringListToTextProperty(&n, 1, &tp);
  117. if (st)
  118. XSetWMName(dpy, win, &tp);
  119. /* Wait for the MapNotify event */
  120. XFlush(dpy);
  121. int depth = DefaultDepth(dpy, DefaultScreen(dpy));
  122. Visual *visual = DefaultVisual(dpy, DefaultScreen(dpy));
  123. /* Make bitmap */
  124. bitmap = XCreateImage(dpy, visual, depth,
  125. ZPixmap, 0, (char *)buffer,
  126. width, height, 32, 0);
  127. /* Init GC */
  128. gc = XCreateGC(dpy, win, 0, NULL);
  129. XSetForeground(dpy, gc, black);
  130. XSelectInput(dpy, win, ExposureMask | KeyPressMask | StructureNotifyMask);
  131. Atom wmDeleteMessage;
  132. wmDeleteMessage = XInternAtom(dpy, "WM_DELETE_WINDOW", False);
  133. XSetWMProtocols(dpy, win, &wmDeleteMessage, 1);
  134. Left = XStringToKeysym ("Left");
  135. Right = XStringToKeysym ("Right");
  136. Up = XStringToKeysym ("Up");
  137. Down = XStringToKeysym ("Down");
  138. Alt = XStringToKeysym ("Alt");
  139. }
  140. static int handle_events(void)
  141. {
  142. XEvent event;
  143. XNextEvent(dpy, &event);
  144. KeySym key;
  145. char text[255];
  146. double coef = 0.05;
  147. if (event.type == KeyPress)
  148. {
  149. XLookupString(&event.xkey,text,255,&key,0);
  150. if (key == Left)
  151. {
  152. double widthX = rightX - leftX;
  153. leftX -= coef*widthX;
  154. rightX -= coef*widthX;
  155. }
  156. else if (key == Right)
  157. {
  158. double widthX = rightX - leftX;
  159. leftX += coef*widthX;
  160. rightX += coef*widthX;
  161. }
  162. else if (key == Down)
  163. {
  164. double heightY = topY - bottomY;
  165. topY += coef*heightY;
  166. bottomY += coef*heightY;
  167. }
  168. else if (key == Up)
  169. {
  170. double heightY = topY - bottomY;
  171. topY -= coef*heightY;
  172. bottomY -= coef*heightY;
  173. }
  174. else {
  175. double widthX = rightX - leftX;
  176. double heightY = topY - bottomY;
  177. if (text[0] == '-')
  178. {
  179. /* Zoom out */
  180. leftX -= (coef/2)*widthX;
  181. rightX += (coef/2)*widthX;
  182. topY += (coef/2)*heightY;
  183. bottomY -= (coef/2)*heightY;
  184. }
  185. else if (text[0] == '+')
  186. {
  187. /* Zoom in */
  188. leftX += (coef/2)*widthX;
  189. rightX -= (coef/2)*widthX;
  190. topY -= (coef/2)*heightY;
  191. bottomY += (coef/2)*heightY;
  192. }
  193. }
  194. if (text[0]=='q') {
  195. return -1;
  196. }
  197. }
  198. if (event.type==ButtonPress) {
  199. /* tell where the mouse Button was Pressed */
  200. printf("You pressed a button at (%i,%i)\n",
  201. event.xbutton.x,event.xbutton.y);
  202. }
  203. return 0;
  204. }
  205. #endif //USE_X11
  206. static void parse_args(int argc, char **argv)
  207. {
  208. int i;
  209. for (i = 1; i < argc; i++) {
  210. if (strcmp(argv[i], "-h") == 0) {
  211. fprintf(stderr, "Usage: %s [-h] [ -width 1024] [-height 768] [-nblocks 16] [-group_size 64] [-no-x11] [-demo] [-frames N] [-pos leftx:rightx:bottomy:topy]\n", argv[0]);
  212. exit(-1);
  213. }
  214. if (strcmp(argv[i], "-width") == 0) {
  215. char *argptr;
  216. width = strtol(argv[++i], &argptr, 10);
  217. }
  218. if (strcmp(argv[i], "-frames") == 0) {
  219. char *argptr;
  220. frames = strtol(argv[++i], &argptr, 10);
  221. }
  222. if (strcmp(argv[i], "-height") == 0) {
  223. char *argptr;
  224. height = strtol(argv[++i], &argptr, 10);
  225. }
  226. if (strcmp(argv[i], "-group_size") == 0) {
  227. char *argptr;
  228. group_size = strtol(argv[++i], &argptr, 10);
  229. }
  230. if (strcmp(argv[i], "-nblocks") == 0) {
  231. char *argptr;
  232. nblocks = strtol(argv[++i], &argptr, 10);
  233. }
  234. if (strcmp(argv[i], "-pos") == 0) {
  235. int ret = sscanf(argv[++i], "%lf:%lf:%lf:%lf", &leftX, &rightX, &bottomY, &topY);
  236. assert(ret == 4);
  237. }
  238. if (strcmp(argv[i], "-demo") == 0) {
  239. demo = 1;
  240. leftX = -50.22749575062760;
  241. rightX = 48.73874621262927;
  242. topY = -49.35016705749115;
  243. bottomY = 49.64891691946615;
  244. }
  245. if (strcmp(argv[i], "-no-x11") == 0) {
  246. #ifdef USE_X11
  247. use_x11 = 0;
  248. #endif
  249. }
  250. }
  251. }
  252. int main(int argc, char **argv) {
  253. #define MAX_DEVICES 20
  254. cl_platform_id platforms[15];
  255. cl_uint num_platforms;
  256. cl_device_id devices[15];
  257. cl_uint num_devices;
  258. cl_context context;
  259. cl_program program;
  260. cl_kernel kernel;
  261. cl_command_queue cq[MAX_DEVICES];
  262. cl_int err;
  263. cl_uint i;
  264. parse_args(argc, argv);
  265. cl_uint block_size = height/nblocks;
  266. assert((height % nblocks) == 0);
  267. assert((width % group_size) == 0);
  268. clGetPlatformIDs(0, NULL, &num_platforms);
  269. if (num_platforms == 0) {
  270. printf("No OpenCL platform found\n");
  271. exit(0);
  272. }
  273. err = clGetPlatformIDs(sizeof(platforms)/sizeof(cl_platform_id), platforms, NULL);
  274. check(err, "clGetPlatformIDs");
  275. unsigned int platform_idx;
  276. for (platform_idx=0; platform_idx<num_platforms; platform_idx++) {
  277. err = clGetDeviceIDs(platforms[platform_idx], CL_DEVICE_TYPE_GPU, sizeof(devices)/sizeof(cl_device_id), devices, &num_devices);
  278. check(err, "clGetDeviceIDs");
  279. if (num_devices != 0)
  280. break;
  281. }
  282. if (num_devices == 0)
  283. error("No OpenCL device found\n");
  284. cl_context_properties properties[] = {CL_CONTEXT_PLATFORM, (cl_context_properties)platforms[platform_idx], 0};
  285. context = clCreateContext(properties, num_devices, devices, NULL, NULL, &err);
  286. check(err, "clCreateContext");
  287. program = clCreateProgramWithSource(context, 1, &kernel_src, NULL, &err);
  288. check(err, "clCreateProgram");
  289. err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
  290. check(err, "clBuildProgram");
  291. kernel = clCreateKernel(program, "mandelbrot_kernel", &err);
  292. check(err, "clCreateKernel");
  293. for (i=0; i<num_devices; i++)
  294. cq[i] = clCreateCommandQueue(context, devices[i], CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE, &err);
  295. check(err, "clCreateCommandQueue");
  296. cl_uint *buffer;
  297. buffer = malloc(height*width*sizeof(cl_uint));
  298. #ifdef USE_X11
  299. if (use_x11)
  300. init_x11(width, height, buffer);
  301. #endif // USE_X11
  302. cl_mem block_handles[nblocks];
  303. cl_uint iby;
  304. for (iby = 0; iby < nblocks; iby++) {
  305. cl_uint *data = &buffer[iby*block_size*width];
  306. block_handles[iby] = clCreateBuffer(context, CL_MEM_WRITE_ONLY | CL_MEM_USE_HOST_PTR, block_size*width*sizeof(cl_uint), data, &err);
  307. }
  308. int stop = 0;
  309. int frame = 0;
  310. while (!stop) {
  311. struct timeval start, end;
  312. gettimeofday(&start, NULL);
  313. if (frames != -1) {
  314. frame++;
  315. stop = (frame == frames);
  316. }
  317. double stepX = (rightX - leftX)/width;
  318. double stepY = (topY - bottomY)/height;
  319. cl_event ker_events[nblocks];
  320. void * ptrs[nblocks];
  321. for (iby = 0; iby < nblocks; iby++) {
  322. err = clSetKernelArg(kernel, 0, sizeof(cl_mem), &block_handles[iby]);
  323. check(err, "clSetKernelArg out");
  324. err = clSetKernelArg(kernel, 1, sizeof(cl_double), &leftX);
  325. check(err, "clSetKernelArg leftX");
  326. err = clSetKernelArg(kernel, 2, sizeof(cl_double), &topY);
  327. check(err, "clSetKernelArg topY");
  328. err = clSetKernelArg(kernel, 3, sizeof(cl_double), &stepX);
  329. check(err, "clSetKernelArg leftX");
  330. err = clSetKernelArg(kernel, 4, sizeof(cl_double), &stepY);
  331. check(err, "clSetKernelArg topY");
  332. err = clSetKernelArg(kernel, 5, sizeof(cl_uint), &maxIt);
  333. check(err, "clSetKernelArg maxIt");
  334. err = clSetKernelArg(kernel, 6, sizeof(cl_uint), &iby);
  335. check(err, "clSetKernelArg iby");
  336. err = clSetKernelArg(kernel, 7, sizeof(cl_uint), &block_size);
  337. check(err, "clSetKernelArg block_size");
  338. size_t local[3] = {group_size, 1, 1};
  339. size_t global[3] = {width, block_size, 1};
  340. #ifdef ROUND_ROBIN
  341. int dev = iby % num_devices;
  342. #else
  343. int dev = 0;
  344. #endif
  345. err = clEnqueueNDRangeKernel(cq[dev], kernel, 3, NULL, global, local, 0, NULL, &ker_events[iby]);
  346. check(err, "clEnqueueNDRangeKernel");
  347. }
  348. for (iby = 0; iby < nblocks; iby++) {
  349. #ifdef ROUND_ROBIN
  350. int dev = iby % num_devices;
  351. #else
  352. int dev = 0;
  353. #endif
  354. ptrs[iby] = clEnqueueMapBuffer(cq[dev], block_handles[iby], CL_FALSE,CL_MAP_READ, 0, block_size*width*sizeof(cl_uint), 1, &ker_events[iby], NULL, NULL);
  355. }
  356. #ifdef ROUND_ROBIN
  357. for (i = 0; i < num_devices; i++)
  358. clFinish(cq[i]);
  359. #else
  360. clFinish(cq[0]);
  361. #endif
  362. gettimeofday(&end, NULL);
  363. double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
  364. #ifdef SHORT_LOG
  365. fprintf(stderr, "%f\n", timing/1000.0);
  366. #else
  367. fprintf(stderr, "Time to generate frame : %f ms\n", timing/1000.0);
  368. fprintf(stderr, "%14.14f:%14.14f:%14.14f:%14.14f\n", leftX, rightX, bottomY, topY);
  369. #endif
  370. #ifdef USE_X11
  371. if (use_x11) {
  372. for (iby = 0; iby < nblocks; iby++) {
  373. pthread_mutex_lock(&mutex);
  374. XPutImage(dpy, win, gc, bitmap,
  375. 0, iby*block_size,
  376. 0, iby*block_size,
  377. width, block_size);
  378. pthread_mutex_unlock(&mutex);
  379. }
  380. }
  381. #endif
  382. for (iby = 0; iby < nblocks; iby++) {
  383. #ifdef ROUND_ROBIN
  384. int dev = iby % num_devices;
  385. #else
  386. int dev = 0;
  387. #endif
  388. clEnqueueUnmapMemObject(cq[dev], block_handles[iby], ptrs[iby], 0, NULL, NULL);
  389. clReleaseEvent(ker_events[iby]);
  390. }
  391. if (demo) {
  392. /* Zoom in */
  393. double zoom_factor = 0.05;
  394. double widthX = rightX - leftX;
  395. double heightY = topY - bottomY;
  396. leftX += (zoom_factor/2)*widthX;
  397. rightX -= (zoom_factor/2)*widthX;
  398. topY -= (zoom_factor/2)*heightY;
  399. bottomY += (zoom_factor/2)*heightY;
  400. }
  401. else {
  402. #ifdef USE_X11
  403. if (use_x11) {
  404. handle_events();
  405. }
  406. #else
  407. stop = 1;
  408. #endif
  409. }
  410. }
  411. #ifdef USE_X11
  412. if (use_x11)
  413. exit_x11();
  414. #endif
  415. for (iby = 0; iby < nblocks; iby++) {
  416. clReleaseMemObject(block_handles[iby]);
  417. }
  418. for (i=0; i<num_devices; i++)
  419. clReleaseCommandQueue(cq[i]);
  420. clReleaseKernel(kernel);
  421. clReleaseProgram(program);
  422. clReleaseContext(context);
  423. return 0;
  424. }