matmul.c 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525
  1. /* StarPU --- Runtime system for heterogeneous multicore architectures.
  2. *
  3. * Copyright (C) 2010,2011, 2015, 2017 University of Bordeaux
  4. * Copyright (C) 2017 Inria
  5. *
  6. * StarPU is free software; you can redistribute it and/or modify
  7. * it under the terms of the GNU Lesser General Public License as published by
  8. * the Free Software Foundation; either version 2.1 of the License, or (at
  9. * your option) any later version.
  10. *
  11. * StarPU is distributed in the hope that it will be useful, but
  12. * WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  14. *
  15. * See the GNU Lesser General Public License in COPYING.LGPL for more details.
  16. */
  17. #ifndef STARPU_NON_BLOCKING_DRIVERS
  18. int main(int argc, const char** argv) {
  19. (void) argv;
  20. (void) argv;
  21. /* testcase does not seem to support blocking drivers */
  22. return 77;
  23. }
  24. #else
  25. #ifdef __APPLE_CC__
  26. #include <OpenCL/opencl.h>
  27. #else
  28. #include <CL/cl.h>
  29. #endif
  30. #include <stdio.h>
  31. #include <string.h>
  32. #include <stdlib.h>
  33. #include <stdint.h>
  34. #include <unistd.h>
  35. #include <assert.h>
  36. #include <math.h>
  37. #include <sys/time.h>
  38. #define error(...) do { fprintf(stderr, "Error: " __VA_ARGS__); exit(EXIT_FAILURE); } while(0)
  39. #define check(exp) do { err = exp; if(err != CL_SUCCESS) { fprintf(stderr, "OpenCL Error (%d): " #exp "\n", err); exit(EXIT_FAILURE); }} while(0)
  40. #define check2(exp) exp; if(err != CL_SUCCESS) { fprintf(stderr, "OpenCL Error (%d): " #exp "\n", err); exit(EXIT_FAILURE); }
  41. #define check3(exp, err) do { if(err != CL_SUCCESS) { fprintf(stderr, "OpenCL Error (%d): " #exp "\n", err); exit(EXIT_FAILURE); } } while(0)
  42. // Thread block size
  43. #define BLOCK_SIZE 16 // Kernel thread-block size
  44. #define WORK_SIZE 64 // Kernel global size in lines of A (or C)
  45. #define TYPE float
  46. // Basic Matrix dimensions
  47. #define WA (128L * BLOCK_SIZE) // Matrix A width
  48. #ifdef STARPU_QUICK_CHECK
  49. #define HA (128L * BLOCK_SIZE) // Matrix A height
  50. #else
  51. #define HA (512L * BLOCK_SIZE) // Matrix A height
  52. #endif
  53. #define WB (128L * BLOCK_SIZE) // Matrix B width
  54. #define HB WA // Matrix B height
  55. #define WC WB // Matrix C width
  56. #define HC HA // Matrix C height
  57. #define BLOCKS (HA / WORK_SIZE)
  58. ////////////////////////////////////////////////////////////////////////////////
  59. // declaration, forward
  60. void printDiff(TYPE*, TYPE*, int, int, int, TYPE);
  61. void computeReference(TYPE*, const TYPE*, const TYPE*, unsigned int, unsigned int, unsigned int);
  62. #define str(x) #x
  63. #define CODE "\
  64. #define TYPE float\n\
  65. __kernel void sgemmNN(int wa, int ha, int wb, __global TYPE* A, __global TYPE* B, __global TYPE* C) {\n\
  66. #define BS 16\n\
  67. #define BLOCK_SIZE 16\n\
  68. int bx = get_group_id(0);\n\
  69. int by = get_group_id(1);\n\
  70. \n\
  71. int tx = get_local_id(0);\n\
  72. int ty = get_local_id(1);\n\
  73. \n\
  74. int gx = get_global_id(0);\n\
  75. int gy = get_global_id(1);\n\
  76. __local float As[BS][BS+1];\
  77. __local float Bs[BS][BS+1];\
  78. \n\
  79. unsigned int block_w = min(wb - bx * BLOCK_SIZE, BLOCK_SIZE);\n\
  80. unsigned int block_h = min(ha - by * BLOCK_SIZE, BLOCK_SIZE);\n\
  81. \n\
  82. int valid = (gx < wb && gy < ha);\n\
  83. \n\
  84. TYPE Csub = (TYPE)0.0;\n\
  85. \n\
  86. int pos = 0;\n\
  87. while (pos < wa) {\n\
  88. unsigned int size = min(wa-pos, BLOCK_SIZE);\n\
  89. if (tx < size && gy < ha)\n\
  90. As[tx][ty] = A[pos + tx + wa * gy];\n\
  91. if (ty < size && gx < wb)\n\
  92. Bs[tx][ty] = B[gx + wb * (pos+ty)];\n\
  93. \n\
  94. barrier(CLK_LOCAL_MEM_FENCE);\n\
  95. \n\
  96. if (valid) {\n\
  97. for (int k = 0; k < size; ++k)\n\
  98. Csub += As[k][ty] * Bs[tx][k];\n\
  99. }\n\
  100. pos += size;\n\
  101. barrier(CLK_LOCAL_MEM_FENCE);\n\
  102. }\n\
  103. \n\
  104. if (valid)\n\
  105. C[wb * gy + gx] = Csub;\n\
  106. }"
  107. static char * code = CODE;
  108. int check = 0;
  109. static void __attribute__((unused)) parse_args(int argc, const char **argv)
  110. {
  111. int i;
  112. for (i = 1; i < argc; i++)
  113. {
  114. if (strcmp(argv[i], "-check") == 0)
  115. {
  116. check = 1;
  117. }
  118. if (strcmp(argv[i], "-h") == 0)
  119. {
  120. printf("usage : %s [-check]\n", argv[0]);
  121. }
  122. }
  123. }
  124. // Round Up Division function
  125. size_t roundUp(int group_size, int global_size) {
  126. int r = global_size % group_size;
  127. if(r == 0) {
  128. return global_size;
  129. } else {
  130. return global_size + group_size - r;
  131. }
  132. }
  133. void fillArray(TYPE* data, int size) {
  134. int i;
  135. const TYPE fScale = (TYPE)(1.0f / (float)RAND_MAX);
  136. for (i = 0; i < size; ++i) {
  137. data[i] = fScale * rand();
  138. }
  139. }
  140. void printArray(float* data, int size) {
  141. int i;
  142. for (i = 0; i < size; ++i) {
  143. printf("%d: %.3f\n", i, data[i]);
  144. }
  145. }
  146. /**
  147. * Compare two float arrays using L2-norm with an epsilon tolerance for equality
  148. * @return shrTRUE if \a reference and \a data are identical, otherwise shrFALSE
  149. * @param reference handle to the reference data / gold image
  150. * @param data handle to the computed data
  151. * @param len number of elements in reference and data
  152. * @param epsilon epsilon to use for the comparison
  153. */
  154. int shrCompareL2fe( const float* reference, const float* data, const unsigned int len, const float epsilon ) {
  155. assert(epsilon >= 0);
  156. float error = 0;
  157. float ref = 0;
  158. unsigned int i;
  159. for(i = 0; i < len; ++i) {
  160. float diff = reference[i] - data[i];
  161. error += diff * diff;
  162. ref += reference[i] * reference[i];
  163. }
  164. float normRef = sqrtf(ref);
  165. if (fabs(ref) < 1e-7) {
  166. #ifdef _DEBUG
  167. fprintf(stderr, "ERROR, reference l2-norm is 0\n");
  168. #endif
  169. return 0;
  170. }
  171. float normError = sqrtf(error);
  172. error = normError / normRef;
  173. int result = error < epsilon;
  174. #ifdef _DEBUG
  175. if( !result) {
  176. fprintf(stderr, "ERROR, l2-norm error %d is greater than epsilon %lf \n", error, epsilon);
  177. }
  178. #endif
  179. return result;
  180. }
  181. int main(int argc, const char** argv) {
  182. cl_uint platform_count;
  183. cl_platform_id platforms[5];
  184. cl_int err = CL_SUCCESS;
  185. unsigned int i, p;
  186. cl_device_type dev_type = CL_DEVICE_TYPE_ALL;
  187. void * ptrs[BLOCKS];
  188. cl_command_queue cqs[BLOCKS];
  189. cl_mem d_A[BLOCKS];
  190. cl_mem d_C[BLOCKS];
  191. cl_mem d_B[BLOCKS];
  192. cl_event GPUDone[BLOCKS];
  193. cl_event GPUExecution[BLOCKS];
  194. struct timeval start, end;
  195. int workOffset[BLOCKS];
  196. int workSize[BLOCKS];
  197. unsigned int sizePerGPU = HC / BLOCKS;
  198. unsigned int sizeMod = HC % BLOCKS;
  199. size_t A_size = WA * HA;
  200. size_t A_mem_size = sizeof(TYPE) * A_size;
  201. TYPE* A_data;
  202. size_t B_size = WB * HB;
  203. size_t B_mem_size = sizeof(TYPE) * B_size;
  204. TYPE* B_data;
  205. size_t C_size = WC * HC;
  206. size_t C_mem_size = sizeof(TYPE) * C_size;
  207. TYPE* C_data;
  208. parse_args(argc, argv);
  209. check(clGetPlatformIDs(5, platforms, &platform_count));
  210. if (platform_count == 0) {
  211. printf("No platform found\n");
  212. exit(77);
  213. }
  214. cl_uint device_count;
  215. cl_uint devs[platform_count];
  216. cl_device_id * devices[platform_count];
  217. cl_context ctx[platform_count];
  218. cl_command_queue * commandQueue[platform_count];
  219. device_count = 0;
  220. for (p=0; p<platform_count; p++) {
  221. cl_platform_id platform = platforms[p];
  222. err = clGetDeviceIDs(platform, dev_type, 0, NULL, &devs[p]);
  223. if (err == CL_DEVICE_NOT_FOUND) {
  224. devs[p] = 0;
  225. continue;
  226. }
  227. if (devs[p] == 0) {
  228. printf("No OpenCL device found\n");
  229. exit(77);
  230. }
  231. if (err != CL_SUCCESS) {
  232. fprintf(stderr, "OpenCL Error (%d) in clGetDeviceIDs()\n", err);
  233. exit(EXIT_FAILURE);
  234. }
  235. if (devs[p] == 0)
  236. continue;
  237. devices[p] = (cl_device_id*)malloc(sizeof(cl_device_id) * devs[p]);
  238. commandQueue[p] = (cl_command_queue*)malloc(sizeof(cl_command_queue) * devs[p]);
  239. check(clGetDeviceIDs(platform, dev_type, devs[p], devices[p], NULL));
  240. cl_context_properties properties[] = {CL_CONTEXT_PLATFORM, (cl_context_properties)platform, 0};
  241. check2(ctx[p] = clCreateContext(properties, devs[p], devices[p], NULL, NULL, &err));
  242. for(i = 0; i < devs[p]; ++i)
  243. {
  244. cl_device_id device = devices[p][i];
  245. char name[2048];
  246. name[0] = '\0';
  247. clGetDeviceInfo(device, CL_DEVICE_NAME, 2048, name, NULL);
  248. printf("Device %d: %s\n", i, name);
  249. commandQueue[p][i] = clCreateCommandQueue(ctx[p], device, CL_QUEUE_PROFILING_ENABLE | CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE, &err);
  250. if (err == CL_INVALID_VALUE) {
  251. fprintf(stderr, "Invalid property for clCreateCommandQueue\n");
  252. exit(77);
  253. }
  254. check3("clCreateCommandQueue", err);
  255. }
  256. device_count += devs[p];
  257. }
  258. if (device_count == 0)
  259. error("No device found\n");
  260. cl_kernel multiplicationKernel[platform_count];
  261. printf("\nUsing Matrix Sizes: A(%lu x %lu), B(%lu x %lu), C(%lu x %lu)\n",
  262. (unsigned long)WA, (unsigned long)HA, (unsigned long)WB, (unsigned long)HB, (unsigned long)WC, (unsigned long)HC);
  263. // allocate host memory for matrices A, B and C
  264. A_data = (TYPE*)malloc(A_mem_size);
  265. if (A_data == NULL) {
  266. perror("malloc");
  267. exit(-1);
  268. }
  269. B_data = (TYPE*)malloc(B_mem_size);
  270. if (B_data == NULL) {
  271. perror("malloc");
  272. exit(-1);
  273. }
  274. C_data = (TYPE*) malloc(C_mem_size);
  275. if (C_data == NULL) {
  276. perror("malloc");
  277. exit(-1);
  278. }
  279. cl_program program[platform_count];
  280. for (p=0; p<platform_count; p++) {
  281. if (devs[p] == 0)
  282. continue;
  283. check2(program[p] = clCreateProgramWithSource(ctx[p], 1, (const char **)&code, NULL, &err));
  284. check(clBuildProgram(program[p], 0, NULL, NULL, NULL, NULL));
  285. check2(multiplicationKernel[p] = clCreateKernel(program[p], "sgemmNN", &err));
  286. }
  287. printf("Initializing data...\n");
  288. srand(2008);
  289. fillArray(A_data, A_size);
  290. fillArray(B_data, B_size);
  291. memset(C_data, 0, C_size);
  292. printf("Computing...\n");
  293. workOffset[0] = 0;
  294. gettimeofday(&start, NULL);
  295. size_t localWorkSize[] = {BLOCK_SIZE, BLOCK_SIZE};
  296. int c = 0;
  297. for (p=0; p<platform_count;p++) {
  298. for (i=0; i<devs[p]; i++) {
  299. check2(d_B[c] = clCreateBuffer(ctx[p], CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, HB * WB * sizeof(TYPE), B_data, &err));
  300. c++;
  301. }
  302. }
  303. for(i=0; i < BLOCKS; ++i)
  304. {
  305. int d = i % device_count;
  306. cl_uint platform = 0;
  307. // determine device platform
  308. int dev = d;
  309. for (platform = 0; platform < platform_count; platform++) {
  310. if ((cl_int)(dev - devs[platform]) < 0)
  311. break;
  312. dev -= devs[platform];
  313. }
  314. assert(platform < platform_count);
  315. workSize[i] = (i < sizeMod) ? sizePerGPU+1 : sizePerGPU;
  316. check2(d_A[i] = clCreateBuffer(ctx[platform], CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, workSize[i] * WA * sizeof(TYPE), &A_data[workOffset[i] * WA], &err));
  317. check2(d_C[i] = clCreateBuffer(ctx[platform], CL_MEM_WRITE_ONLY | CL_MEM_USE_HOST_PTR, workSize[i] * WC * sizeof(TYPE), &C_data[workOffset[i] * WC], &err));
  318. check(clSetKernelArg(multiplicationKernel[platform], 0, sizeof(cl_int), &workSize[i]));
  319. check(clSetKernelArg(multiplicationKernel[platform], 1, sizeof(cl_int), &workSize[i]));
  320. check(clSetKernelArg(multiplicationKernel[platform], 2, sizeof(cl_int), &workSize[i]));
  321. check(clSetKernelArg(multiplicationKernel[platform], 3, sizeof(cl_mem), (void *) &d_A[i]));
  322. check(clSetKernelArg(multiplicationKernel[platform], 4, sizeof(cl_mem), (void *) &d_B[d]));
  323. check(clSetKernelArg(multiplicationKernel[platform], 5, sizeof(cl_mem), (void *) &d_C[i]));
  324. size_t globalWorkSize[] = {roundUp(BLOCK_SIZE,WC), roundUp(BLOCK_SIZE,workSize[i])};
  325. check(clEnqueueNDRangeKernel(commandQueue[platform][dev], multiplicationKernel[platform], 2, NULL, globalWorkSize, localWorkSize, 0, NULL, &GPUExecution[i]));
  326. // Non-blocking copy of result from device to host
  327. cqs[i] = commandQueue[platform][dev];
  328. check2(ptrs[i] = clEnqueueMapBuffer(cqs[i], d_C[i], CL_FALSE, CL_MAP_READ, 0, WC * sizeof(TYPE) * workSize[i], 1, &GPUExecution[i], &GPUDone[i], &err));
  329. if(i+1 < BLOCKS)
  330. workOffset[i + 1] = workOffset[i] + workSize[i];
  331. }
  332. // CPU sync with GPU
  333. for (p=0; p<platform_count;p++) {
  334. cl_uint dev;
  335. for (dev=0; dev<devs[p]; dev++) {
  336. clFinish(commandQueue[p][dev]);
  337. }
  338. }
  339. gettimeofday(&end, NULL);
  340. double timing = (double)((end.tv_sec - start.tv_sec)*1000000 + (end.tv_usec - start.tv_usec));
  341. double dSeconds = timing/1000/1000;
  342. double dNumOps = 2.0 * (double)WA * (double)HA * (double)WB;
  343. double gflops = 1.0e-9 * dNumOps/dSeconds;
  344. printf("Throughput = %.4f GFlops/s, Time = %.5f s, Size = %.0f, NumDevsUsed = %d, Blocks = %ld, Workgroup = %zu\n",
  345. gflops, dSeconds, dNumOps, device_count, BLOCKS, localWorkSize[0] * localWorkSize[1]);
  346. // compute reference solution
  347. if (check) {
  348. printf("Comparing results with CPU computation... ");
  349. TYPE* reference = (TYPE*)malloc(C_mem_size);
  350. computeReference(reference, A_data, B_data, HA, WA, WB);
  351. // check result
  352. int res = shrCompareL2fe(reference, C_data, C_size, 1.0e-6f);
  353. if (res == 0) {
  354. printf("\n\n");
  355. printDiff(reference, C_data, WC, HC, 100, 1.0e-5f);
  356. }
  357. else printf("PASSED\n\n");
  358. free(reference);
  359. }
  360. for(i = 0; i < BLOCKS; i++)
  361. {
  362. clEnqueueUnmapMemObject(cqs[i], d_C[i], ptrs[i], 0, NULL, NULL);
  363. }
  364. for(i = 0; i < BLOCKS; i++)
  365. {
  366. clFinish(cqs[i]);
  367. }
  368. for (i=0; i<device_count; i++) {
  369. clReleaseMemObject(d_B[i]);
  370. }
  371. for(i = 0; i < BLOCKS; i++)
  372. {
  373. clReleaseMemObject(d_A[i]);
  374. clReleaseMemObject(d_C[i]);
  375. clReleaseEvent(GPUExecution[i]);
  376. clReleaseEvent(GPUDone[i]);
  377. }
  378. for (p=0; p<platform_count;p++) {
  379. if (devs[p] == 0)
  380. continue;
  381. check(clReleaseKernel(multiplicationKernel[p]));
  382. check(clReleaseProgram(program[p]));
  383. check(clReleaseContext(ctx[p]));
  384. cl_uint k;
  385. for(k = 0; k < devs[p]; ++k)
  386. {
  387. check(clReleaseCommandQueue(commandQueue[p][k]));
  388. }
  389. }
  390. free(A_data);
  391. free(B_data);
  392. free(C_data);
  393. return 0;
  394. }
  395. void printDiff(TYPE *data1, TYPE *data2, int width, int height, int listLength, TYPE listTol) {
  396. printf("Listing first %d Differences > %.6f...\n", listLength, listTol);
  397. int i,j,k;
  398. int error_count=0;
  399. for (j = 0; j < height; j++) {
  400. if (error_count < listLength) {
  401. printf("\n Row %d:\n", j);
  402. }
  403. for (i = 0; i < width; i++) {
  404. k = j * width + i;
  405. float diff = fabs(data1[k] - data2[k]);
  406. if (diff > listTol) {
  407. if (error_count < listLength) {
  408. printf(" Loc(%d,%d)\tCPU=%.5f\tGPU=%.5f\tDiff=%.6f\n", i, j, data1[k], data2[k], diff);
  409. }
  410. error_count++;
  411. }
  412. }
  413. }
  414. printf(" \n Total Errors = %d\n\n", error_count);
  415. }
  416. /**
  417. * Compute reference data set
  418. * C = A * B
  419. * @param C reference data, computed but preallocated
  420. * @param A matrix A as provided to device
  421. * @param B matrix B as provided to device
  422. * @param hA height of matrix A
  423. * @param wB width of matrix B
  424. */
  425. void computeReference(TYPE* C, const TYPE* A, const TYPE* B, unsigned int hA, unsigned int wA, unsigned int wB) {
  426. unsigned int i,j,k;
  427. for (i = 0; i < hA; ++i)
  428. for (j = 0; j < wB; ++j) {
  429. double sum = 0;
  430. for (k = 0; k < wA; ++k) {
  431. double a = A[i * wA + k];
  432. double b = B[k * wB + j];
  433. sum += a * b;
  434. }
  435. C[i * wB + j] = (TYPE)sum;
  436. }
  437. }
  438. #endif /* STARPU_NON_BLOCKING_DRIVERS */