Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /trunk/test/MIP/mip_opencl.c
ViewVC logotype

Annotation of /trunk/test/MIP/mip_opencl.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 191 - (view) (download) (as text)

1 : lamonts 177 #include <OpenCL/OpenCl.h>
2 :     #include <assert.h>
3 :     #include <stdio.h>
4 :     #include <stdlib.h>
5 :     #include <sys/sysctl.h>
6 :     #include <sys/stat.h>
7 :    
8 :     #include <teem/nrrd.h>
9 :    
10 :     #define SIZE 200
11 :    
12 :     /*typedef float vec3[3];
13 :    
14 :     typedef struct {
15 :     int degree;
16 :     float coeff[];
17 :     } polynomial;
18 :    
19 :     typedef struct {
20 :     int support;
21 :     polynomial *segments[];
22 :     } kernel; */
23 :    
24 :    
25 :     int device_stats(cl_device_id device_id){
26 :    
27 :     int err,i;
28 :     size_t returned_size;
29 :    
30 :     // Report the device vendor and device name
31 :     //
32 :     cl_char vendor_name[1024] = {0};
33 :     cl_char device_name[1024] = {0};
34 :     cl_char device_profile[1024] = {0};
35 :     cl_char device_extensions[1024] = {0};
36 :     cl_device_local_mem_type local_mem_type;
37 :    
38 :     cl_ulong global_mem_size, global_mem_cache_size;
39 :     cl_ulong max_mem_alloc_size;
40 :    
41 :     cl_uint clock_frequency, vector_width, max_compute_units;
42 :    
43 :     size_t max_work_item_dims,max_work_group_size, max_work_item_sizes[3];
44 :    
45 :     cl_uint vector_types[] = {CL_DEVICE_PREFERRED_VECTOR_WIDTH_CHAR, CL_DEVICE_PREFERRED_VECTOR_WIDTH_SHORT, CL_DEVICE_PREFERRED_VECTOR_WIDTH_INT,CL_DEVICE_PREFERRED_VECTOR_WIDTH_LONG,CL_DEVICE_PREFERRED_VECTOR_WIDTH_FLOAT,CL_DEVICE_PREFERRED_VECTOR_WIDTH_DOUBLE};
46 :     char *vector_type_names[] = {"char","short","int","long","float","double"};
47 :    
48 :     err = clGetDeviceInfo(device_id, CL_DEVICE_VENDOR, sizeof(vendor_name), vendor_name, &returned_size);
49 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_NAME, sizeof(device_name), device_name, &returned_size);
50 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_PROFILE, sizeof(device_profile), device_profile, &returned_size);
51 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_EXTENSIONS, sizeof(device_extensions), device_extensions, &returned_size);
52 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_LOCAL_MEM_TYPE, sizeof(local_mem_type), &local_mem_type, &returned_size);
53 :    
54 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_GLOBAL_MEM_SIZE, sizeof(global_mem_size), &global_mem_size, &returned_size);
55 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_GLOBAL_MEM_CACHELINE_SIZE, sizeof(global_mem_cache_size), &global_mem_cache_size, &returned_size);
56 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_MAX_MEM_ALLOC_SIZE, sizeof(max_mem_alloc_size), &max_mem_alloc_size, &returned_size);
57 :    
58 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_MAX_CLOCK_FREQUENCY, sizeof(clock_frequency), &clock_frequency, &returned_size);
59 :    
60 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(max_work_group_size), &max_work_group_size, &returned_size);
61 :    
62 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS, sizeof(max_work_item_dims), &max_work_item_dims, &returned_size);
63 :    
64 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_MAX_WORK_ITEM_SIZES, sizeof(max_work_item_sizes), max_work_item_sizes, &returned_size);
65 :    
66 :     err|= clGetDeviceInfo(device_id, CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(max_compute_units), &max_compute_units, &returned_size);
67 :    
68 :     printf("Vendor: %s\n", vendor_name);
69 :     printf("Device Name: %s\n", device_name);
70 :     printf("Profile: %s\n", device_profile);
71 :     printf("Supported Extensions: %s\n\n", device_extensions);
72 :    
73 :     printf("Local Mem Type (Local=1, Global=2): %i\n",(int)local_mem_type);
74 :     printf("Global Mem Size (MB): %i\n",(int)global_mem_size/(1024*1024));
75 :     printf("Global Mem Cache Size (Bytes): %i\n",(int)global_mem_cache_size);
76 :     printf("Max Mem Alloc Size (MB): %ld\n",(long int)max_mem_alloc_size/(1024*1024));
77 :    
78 :     printf("Clock Frequency (MHz): %i\n\n",clock_frequency);
79 :    
80 :     for(i=0;i<6;i++){
81 :     err|= clGetDeviceInfo(device_id, vector_types[i], sizeof(clock_frequency), &vector_width, &returned_size);
82 :     printf("Vector type width for: %s = %i\n",vector_type_names[i],vector_width);
83 :     }
84 :    
85 :     printf("\nMax Work Group Size: %lu\n",max_work_group_size);
86 :     //printf("Max Work Item Dims: %lu\n",max_work_item_dims);
87 :     //for(size_t i=0;i<max_work_item_dims;i++)
88 :     // printf("Max Work Items in Dim %lu: %lu\n",(long unsigned)(i+1),(long unsigned)max_work_item_sizes[i]);
89 :    
90 :     printf("Max Compute Units: %i\n",max_compute_units);
91 :     printf("\n");
92 :    
93 :     return CL_SUCCESS;
94 :     }
95 :     //Loads the Kernel from a file
96 :     char * loadKernel (const char * filename)
97 :     {
98 :     struct stat statbuf;
99 :     FILE *fh;
100 :     char *source;
101 :    
102 :     fh = fopen(filename, "r");
103 :     if (fh == 0)
104 :     return 0;
105 :    
106 :     stat(filename, &statbuf);
107 :     source = (char *) malloc(statbuf.st_size + 1);
108 :     fread(source, statbuf.st_size, 1, fh);
109 :     source[statbuf.st_size] = '\0';
110 :    
111 :     return source;
112 :     }
113 :     float det3x3(float a, float b, float c, float d, float e, float f, float g, float h, float i)
114 :     {
115 :     return ( (a)*(e)*(i)
116 :     + (d)*(h)*(c)
117 :     + (g)*(b)*(f)
118 :     - (g)*(e)*(c)
119 :     - (d)*(b)*(i)
120 :     - (a)*(h)*(f));
121 :     }
122 :     float det4x4(cl_float16 m)
123 :     {
124 :     return (m[ 0] * det3x3(m[ 5], m[ 6], m[ 7],
125 :     m[ 9], m[10], m[11],
126 :     m[13], m[14], m[15])
127 :    
128 :     - m[ 1] * det3x3(m[ 4], m[ 6], m[ 7],
129 :     m[ 8], m[10], m[11],
130 :     m[12], m[14], m[15])
131 :     + m[ 2] * det3x3(m[ 4], m[ 5], m[ 7],
132 :     m[ 8], m[ 9], m[11],
133 :     m[12], m[13], m[15])
134 :    
135 :     - m[ 3] * det3x3(m[ 4], m[ 5], m[ 6],
136 :     m[ 8], m[ 9], m[10],
137 :     m[12], m[13], m[14]));
138 :    
139 :    
140 :    
141 :     }
142 :     void invMatrix(cl_float16 m, cl_float16 i)
143 :     {
144 :     float det = det4x4(m);
145 :    
146 :    
147 :     i[0] = det3x3(m[5],m[ 6],m[ 7],
148 :     m[ 9],m[10],m[11],
149 :     m[13],m[14],m[15])/det;
150 :    
151 :     i[ 1] = -det3x3(m[ 1],m[ 2],m[ 3],
152 :     m[ 9],m[10],m[11],
153 :     m[13],m[14],m[15])/det;
154 :    
155 :     i[ 2] = det3x3(m[ 1],m[ 2],m[ 3],
156 :     m[ 5],m[ 6],m[ 7],
157 :     m[13],m[14],m[15])/det;
158 :    
159 :     i[ 3] = -det3x3(m[ 1],m[ 2],m[ 3],
160 :     m[ 5],m[ 6],m[ 7],
161 :     m[ 9],m[10],m[11])/det;
162 :    
163 :     i[ 4] = -det3x3(m[ 4],m[ 6],m[ 7],
164 :     m[ 8],m[10],m[11],
165 :     m[12],m[14],m[15])/det;
166 :    
167 :     i[ 5] = det3x3(m[ 0],m[ 2],m[ 3],
168 :     m[ 8],m[10],m[11],
169 :     m[12],m[14],m[15])/det;
170 :    
171 :     i[ 6] = -det3x3(m[ 0],m[ 2],m[ 3],
172 :     m[ 4],m[ 6],m[ 7],
173 :     m[12],m[14],m[15])/det;
174 :    
175 :     i[ 7] = det3x3(m[ 0],m[ 2],m[ 3],
176 :     m[ 4],m[ 6],m[ 7],
177 :     m[ 8],m[10],m[11])/det;
178 :    
179 :     i[ 8] = det3x3(m[ 4],m[ 5],m[ 7],
180 :     m[ 8],m[ 9],m[11],
181 :     m[12],m[13],m[15])/det;
182 :    
183 :     i[ 9] = -det3x3(m[ 0],m[ 1],m[ 3],
184 :     m[ 8],m[ 9],m[11],
185 :     m[12],m[13],m[15])/det;
186 :    
187 :     i[10] = det3x3(m[ 0],m[ 1],m[ 3],
188 :     m[ 4],m[ 5],m[ 7],
189 :     m[12],m[13],m[15])/det;
190 :    
191 :     i[11] = -det3x3(m[ 0],m[ 1],m[ 3],
192 :     m[ 4],m[ 5],m[ 7],
193 :     m[ 8],m[ 9],m[11])/det;
194 :    
195 :     i[12] = -det3x3(m[ 4],m[ 5],m[ 6],
196 :     m[ 8],m[ 9],m[10],
197 :     m[12],m[13],m[14])/det;
198 :    
199 :     i[13] = det3x3(m[ 0],m[ 1],m[ 2],
200 :     m[ 8],m[ 9],m[10],
201 :     m[12],m[13],m[14])/det;
202 :    
203 :     i[14] = -det3x3(m[ 0],m[ 1],m[ 2],
204 :     m[ 4],m[ 5],m[ 6],
205 :     m[12],m[13],m[14])/det;
206 :    
207 :     i[15] = det3x3(m[ 0],m[ 1],m[ 2],
208 :     m[ 4],m[ 5],m[ 6],
209 :     m[ 8],m[ 9],m[10])/det;
210 :     }
211 :     void printMatrix(float * matrix, int rowSize)
212 :     {
213 :     int index = 0, end = 1, arraySize = rowSize * rowSize;
214 :    
215 :     for(index = 0; index < arraySize; index++)
216 :     {
217 :     if(end == rowSize)
218 :     {
219 :     printf(" %.2f\n",matrix[index]);
220 :     end = 1;
221 :     }
222 :     else
223 :     {
224 :     printf(" %.2f ",matrix[index]);
225 :     end++;
226 :     }
227 :     }
228 :     printf("\n");
229 :     }
230 :     void loadTransformMatrix(Nrrd * nin, cl_float16 transformMatrix)
231 :     {
232 :     int i,j, size = nin->spaceDim;
233 :     NrrdAxisInfo axisInfo;
234 :    
235 :     //Image axis Scaling and Rotation
236 :     for(i = 0; i < size; i++)
237 :     {
238 :     axisInfo = nin->axis[i];
239 :     for(j = 0; j < size; j++)
240 :     {
241 :     transformMatrix[ (size+ 1) * j + i] = axisInfo.spaceDirection[j];
242 :     }
243 :    
244 :     //Image Location
245 :     transformMatrix[ (i * (size + 1)) + size] = nin->spaceOrigin[i];
246 :    
247 :     //Bottom row of the Transform Matrix
248 :     transformMatrix[((size + 1) * (size)) + i ] = 0;
249 :     }
250 :     transformMatrix[((size + 1) * (size)) + size ] = 1;
251 :     }
252 :     Nrrd * loadNrrdFile(char * filename)
253 :     {
254 :     /* create a nrrd; at this point this is just an empty container */
255 :     Nrrd * nin;
256 :    
257 :     nin = nrrdNew();
258 :     char *err;
259 :    
260 :     /* read in the nrrd from file */
261 :     if (nrrdLoad(nin, filename, NULL)) {
262 :     err = biffGetDone(NRRD);
263 :     fprintf(stderr, "Mip: trouble reading \"%s\":\n%s", filename, err);
264 :     free(err);
265 :     return NULL;
266 :     }
267 :    
268 :     /* say something about the array
269 :     printf("Mip: \"%s\" is a %d-dimensional nrrd of type %d (%s)\n",
270 :     filename, nin->dim, nin->type,
271 :     airEnumStr(nrrdType, nin->type));
272 :     printf("Mip: the array contains %d elements, each %d bytes in size\n",
273 :     (int)nrrdElementNumber(nin), (int)nrrdElementSize(nin));*/
274 :    
275 :     return nin;
276 :    
277 :     }
278 :     int exe_MIP_Kernel(Nrrd * nin, float stepSize, cl_float4 eyeVec, cl_float4 origVec,
279 :     cl_float4 cVec, cl_float4 rVec, float * h1, float * h2, float * out)
280 :     {
281 :    
282 :     cl_program program;
283 :     cl_kernel kernel;
284 :    
285 :     cl_command_queue queue;
286 :     cl_context context;
287 :    
288 :     cl_device_id cpu = NULL, device = NULL;
289 :    
290 :     cl_int err = 0;
291 :    
292 :     cl_float16 transformMatrix;
293 :     cl_float16 inverseMatrix;
294 :    
295 :     int imageDataSize = (int)nrrdElementNumber(nin);
296 :    
297 : lamonts 191 float * data = (float *)nin->data;
298 :     printf("Data Image: %f\n", (float)data[4* nin->axis[1].size * nin->axis[2].size + 5 * nin->axis[2].size + 2]);
299 :    
300 : lamonts 177 cl_mem imageData_mem, out_mem, h1_mem, h2_mem;
301 :    
302 :    
303 :     /** Setup Device **/
304 :     err = clGetDeviceIDs(NULL,CL_DEVICE_TYPE_CPU,1,&cpu,NULL);
305 :     assert(err==CL_SUCCESS);
306 :    
307 :     err = clGetDeviceIDs(NULL,CL_DEVICE_TYPE_GPU,1,&device,NULL);
308 :     //if(err != CL_SUCCESS)
309 :     device = cpu;
310 :    
311 :     assert(device);
312 :    
313 :    
314 :     /** Retrieve Information about the device
315 :     cl_char vendor_name[1024] = {0};
316 :     cl_char device_name[1024] = {0};
317 :     err = clGetDeviceInfo(device, CL_DEVICE_VENDOR, sizeof(vendor_name), vendor_name, &returned_size);
318 :     err|= clGetDeviceInfo(device, CL_DEVICE_NAME, sizeof(device_name), device_name, &returned_size);
319 :     printf("Connecting to %s %s...\n", vendor_name, device_name);
320 :     device_stats(device); */
321 :    
322 :    
323 :     /* Setup Context and Command Queue */
324 :     context = clCreateContext(0,1,&device,NULL,NULL,&err);
325 :     assert(err == CL_SUCCESS);
326 :    
327 :     queue = clCreateCommandQueue(context,device,0,NULL);
328 :    
329 :     /** Load the Kernel and Program **/
330 :     const char * filename = "mip.cl";
331 :     char * kernel_source = loadKernel(filename);
332 :    
333 :     assert(kernel_source != 0);
334 :    
335 :     program = clCreateProgramWithSource(context,1,(const char **)&kernel_source,NULL,&err);
336 :    
337 :     assert(err == CL_SUCCESS);
338 :    
339 :     err = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
340 :    
341 :    
342 :     /** Retrieve information about the program build to check for any possible errors **/
343 :     char * build_log;
344 :     size_t log_size;
345 :    
346 :     // First call to know the proper size
347 :     clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
348 :     build_log = (char *) malloc(log_size+1);
349 :     // Second call to get the log
350 :     clGetProgramBuildInfo(program, device, CL_PROGRAM_BUILD_LOG, log_size, build_log, NULL);
351 :     build_log[log_size] = '\0';
352 :     printf("\nBuild Log:\n%s\n",build_log);
353 :     free(build_log);
354 :    
355 :     assert(err == CL_SUCCESS);
356 :    
357 :     kernel = clCreateKernel(program,"raycast",&err);
358 :    
359 :     assert(err == CL_SUCCESS);
360 :    
361 :     /** Memory Allocation for the Matrices **/
362 :    
363 :     h1_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(float) * 4,NULL,NULL);
364 :     err |= clEnqueueWriteBuffer(queue,h1_mem,CL_TRUE,0,sizeof(float) * 4,
365 :     (void *)h1 ,0,NULL,NULL);
366 :    
367 :     h2_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(float) * 4,NULL,NULL);
368 :     err |= clEnqueueWriteBuffer(queue,h2_mem,CL_TRUE,0,sizeof(float) * 4,
369 :     (void *)h2 ,0,NULL,NULL);
370 :    
371 :    
372 : lamonts 191 imageData_mem = clCreateBuffer(context,CL_MEM_READ_ONLY,sizeof(float) * imageDataSize,NULL,NULL);
373 :     err |= clEnqueueWriteBuffer(queue,imageData_mem,CL_TRUE,0,sizeof(float) * imageDataSize,
374 : lamonts 177 nin->data ,0,NULL,NULL);
375 :    
376 :     //Load the transformMatrix
377 :     loadTransformMatrix(nin,transformMatrix);
378 :     invMatrix(transformMatrix,inverseMatrix);
379 :    
380 :     err |= clEnqueueWriteBuffer(queue,imageData_mem,CL_TRUE,0,imageDataSize,
381 :     nin->data ,0,NULL,NULL);
382 :    
383 :     assert(err == CL_SUCCESS);
384 :    
385 :     out_mem = clCreateBuffer(context,CL_MEM_READ_WRITE, sizeof(float) * (SIZE *SIZE),NULL,NULL);
386 :    
387 :     clFinish(queue);
388 :    
389 :     size_t global_work_size[2], local_work_size[2];
390 :    
391 :     global_work_size[0] = 256;
392 :     global_work_size[1] = 256;
393 :    
394 :     local_work_size[0] = 1;
395 :     local_work_size[1] = 1;
396 :    
397 :     err =clSetKernelArg(kernel,0,sizeof(cl_mem), &imageData_mem);
398 :     err |=clSetKernelArg(kernel,1,sizeof(cl_mem), &h1_mem);
399 :     err |=clSetKernelArg(kernel,2,sizeof(cl_mem), &h2_mem);
400 :     err |=clSetKernelArg(kernel,3,sizeof(cl_mem), &out_mem);
401 :     err |=clSetKernelArg(kernel,4,sizeof(cl_float4), origVec);
402 :     err |=clSetKernelArg(kernel,5,sizeof(cl_float4), eyeVec);
403 :     err |=clSetKernelArg(kernel,6,sizeof(cl_float4), cVec);
404 :     err |=clSetKernelArg(kernel,7,sizeof(cl_float4), rVec);
405 :     err |=clSetKernelArg(kernel,8,sizeof(cl_float16), &inverseMatrix);
406 :     err |=clSetKernelArg(kernel,9,sizeof(float), &stepSize);
407 :     err |=clSetKernelArg(kernel,10,sizeof(int), &nin->axis[1].size);
408 :     err |=clSetKernelArg(kernel,11,sizeof(int), &nin->axis[2].size);
409 :    
410 :     assert(err == CL_SUCCESS);
411 :    
412 :     /** Retrieve the Recommend Work Group Size */
413 :     size_t thread_size;
414 :     clGetKernelWorkGroupInfo(kernel,device,CL_KERNEL_WORK_GROUP_SIZE,
415 :     sizeof(size_t),&thread_size,NULL);
416 :     printf("Recommended Size: %lu\n",thread_size);
417 :    
418 :    
419 :     err = clEnqueueNDRangeKernel(queue,kernel,2,NULL,global_work_size,
420 :     local_work_size,0,NULL,NULL);
421 :    
422 :     assert(err == CL_SUCCESS);
423 :    
424 :     clFinish(queue);
425 :    
426 :     err = clEnqueueReadBuffer(queue,out_mem,CL_TRUE,0, sizeof(float) * (SIZE *SIZE),out,0,NULL,NULL);
427 : lamonts 191
428 :     printMatrix(out,20);
429 : lamonts 177
430 :     clReleaseKernel(kernel);
431 :     clReleaseProgram(program);
432 :     clReleaseCommandQueue(queue);
433 :     clReleaseContext(context);
434 :    
435 :     clReleaseMemObject(imageData_mem);
436 :     clReleaseMemObject(h1_mem);
437 :     clReleaseMemObject(h2_mem);
438 :     clReleaseMemObject(out_mem);
439 :    
440 :     return CL_SUCCESS;
441 :     }
442 :     int main (int argc, char ** argv)
443 :     {
444 :     //Declaring and initializing input variables
445 :     Nrrd * nin;
446 :     char * dataFile = "txs.nrrd";
447 :     cl_float4 eyeVector = {25,15,10};
448 :     cl_float4 origVector = {8.83877,2.5911,7.65275};
449 :     cl_float4 cVector = {-0.0151831,0.0278357,0};
450 :     cl_float4 rVector = {0.0074887,0.00408474,-0.0305383};
451 :     float stepSize = 0.1;
452 :     float h1[] = {0.666667,0,-1,0.5};
453 :     float h2[] = {1.33333, -2, 1,-0.166667};
454 :     float * out;
455 :    
456 :     out = (float *) malloc(sizeof(float) * (SIZE * SIZE));
457 :    
458 :     nin = loadNrrdFile(dataFile);
459 :    
460 :     exe_MIP_Kernel(nin,stepSize,eyeVector,origVector,
461 :     cVector,rVector,h1,h2,out);
462 :    
463 :    
464 :     return 0;
465 :     }

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0