source: CIVL/examples/cuda/matMult1.cvl@ 5c27aa5

1.23 2.0 main test-branch
Last change on this file since 5c27aa5 was e3151da, checked in by Ziqing Luo <ziqing@…>, 11 years ago

re-organized example directory

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@1763 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 11.3 KB
Line 
1/***********************************************************************
2* FILENAME: MM.cu
3* Matrix Multiplication
4* Matrix operands have row-major order.
5*
6* C = A * B
7* Multiplies two square matrices (NxN * NxN).
8* Matrix values have type double.
9*
10* A simple CUDA program has a basic workflow:
11* 1) Initialize matrix operands as double-precision arrays on host (CPU).
12* 2) Copy operands from host memory to GPU memory.
13* 3) Apply matrix operaton to operands on GPU
14* 4) Copy result from GPU memory to host memory.
15*
16*
17* CUDA C Programming Guide Version 4.2 (3.2.3, p.22):
18* http://developer.download.nvidia.com/compute/DevZone/docs/html/C/doc/CUDA_C_Programming_Guide.pdf
19*
20* MM with linearized matrix operands:
21* http://www.hpcwire.com/hpcwire/2008-10-08/compilers_and_more_programming_gpus_today.html
22*
23*************************************************************************/
24// online source: https://www.rcac.purdue.edu/userinfo/resources/carter/compile/MM.cu
25
26#include "civlc.h"
27#include "civlc-cuda.cvl"
28#include <stdio.h>
29#include <stdlib.h>
30
31$input int N; /* size of square matrix */
32$input int TILE_WIDTH;
33
34
35/* MM kernel using global (not shared) memory. */
36void _kernel_myMM_global (dim3 gridDim, dim3 blockDim, cudaStream_t s, const double * const A, const double * const B, double *C, int width) {
37
38 void _kernel (_kernelInstance *this, cudaEvent_t e) {
39
40 _waitInQueue(this, e);
41
42 void _block(uint3 blockIdx) {
43 int _numThreads = blockDim.x * blockDim.y * blockDim.z;
44 $gbarrier _block_barrier = $gbarrier_create($here, _numThreads);
45
46 void _thread(uint3 threadIdx) {
47
48 int _tid = _index(blockDim, threadIdx);
49 $barrier _b = $barrier_create($here, _block_barrier, _tid);
50
51 /* Get row and column from block and thread IDs */
52 int row = (blockDim.y*blockIdx.y) + threadIdx.y;
53 int col = (blockDim.x*blockIdx.x) + threadIdx.x;
54
55 /* Initialize result of one element which one thread computes. */
56 double result=0.0;
57
58 /* Compute one element of the matrix product. */
59 for (int i = 0; i < width; ++i)
60 result += A[row*width + i] * B[i*width + col];
61
62 /* Store the result of one matrix element in matrix C. */
63 C[row * width + col] = result;
64 $barrier_destroy(_b);
65 }
66 _runProcs(blockDim, _thread);
67 $gbarrier_destroy(_block_barrier);
68 }
69 _runProcs(gridDim, _block);
70 _kernelFinish(this);
71 }
72 _enqueueKernel(s, _kernel);
73}
74
75
76/* MM kernel using shared memory. */
77void _kernel_myMM_shared (dim3 gridDim, dim3 blockDim, cudaStream_t s, const double * const A, const double * const B, double* C, int width) {
78
79 void _kernel (_kernelInstance *this, cudaEvent_t e) {
80
81 _waitInQueue(this, e);
82
83 void _block(uint3 blockIdx) {
84 int _numThreads = blockDim.x * blockDim.y * blockDim.z;
85 $gbarrier _block_barrier = $gbarrier_create($here, _numThreads);
86
87 double A_shared[TILE_WIDTH][TILE_WIDTH];
88 double B_shared[TILE_WIDTH][TILE_WIDTH];
89
90 void _thread(uint3 threadIdx) {
91
92 int _tid = _index(blockDim, threadIdx);
93 $barrier _b = $barrier_create($here, _block_barrier, _tid);
94
95 int bx = blockIdx.x; int by = blockIdx.y;
96 int tx = threadIdx.x; int ty = threadIdx.y;
97
98 /* Identify the row and column of the C element to work on. */
99 int row = by * TILE_WIDTH + ty;
100 int col = bx * TILE_WIDTH + tx;
101
102 double result = 0.0;
103
104 /* Loop over the A and B tiles required to compute the C element. */
105 for (int phase = 0; phase < width/TILE_WIDTH; ++phase) {
106 /* Shared effort: loading of A and B tiles into shared memory. */
107
108 A_shared[ty][tx] = A[row*width + (phase*TILE_WIDTH + tx)];
109
110 B_shared[ty][tx] = B[col + (phase*TILE_WIDTH + ty)*width];
111
112 $barrier_call(_b);
113
114 for (int k = 0; k < TILE_WIDTH; ++k) {
115 result += A_shared[ty][k] * B_shared[k][tx];
116 }
117
118 $barrier_call(_b);
119
120 }
121 C[row*width+col] = result;
122 $barrier_destroy(_b);
123 }
124 _runProcs(blockDim, _thread);
125 $gbarrier_destroy(_block_barrier);
126 }
127 _runProcs(gridDim, _block);
128 _kernelFinish(this);
129 }
130 _enqueueKernel(s, _kernel);
131}
132
133
134/************************************************************************/
135/************************************************************************/
136/************************************************************************/
137
138
139int main (int argc, char** argv) {
140
141 int _main ( void ) {
142
143 /* Set device based on input from command line */
144 /*
145 if (argc > 1) {
146 if (cudaSetDevice(atoi(argv[1])) != cudaSuccess) {
147 int num_devices;
148 cudaGetDeviceCount(&num_devices);
149 fprintf(stderr, "Error initializing device %s,\
150 device value must be 0-%d\n", argv[1], (num_devices-1));
151 return 0;
152 }
153 } else {
154 fprintf(stderr, "Usage: %s gpu_device\n", argv[0]);
155 return 0;
156 }
157 */
158
159 /* Declare CPU arrays. */
160 double A[N*N],B[N*N],C[N*N]; /* linearized CPU double arrays */
161 double cpuResult[N*N], gpuGlobalResult[N*N], gpuSharedResult[N*N];
162 int r,c;
163
164 /* Declare GPU arrays. */
165 double *G_A,*G_B,*G_C; /* linearized GPU double arrays */
166 size_t size_a,size_b,size_c; /* size of linearized array in bytes */
167 size_a = size_b = size_c = N*N;
168
169 /* Setup a clock. */
170 cudaEvent_t start, stop;
171 float CPU_elapsedtime, GPU_global_elapsedtime, GPU_shared_elapsedtime;
172 cudaEventCreate(&start);
173 cudaEventCreate(&stop);
174
175
176
177
178 /* 1) Initialize matrix operands as double-precision arrays on host (CPU). */
179 for (r=0;r<N;++r)
180 for (c=0;c<N;++c) {
181 A[r*N+c] = 1.0;
182 B[r*N+c] = 1.0;
183 C[r*N+c] = 0.0;
184 }
185
186
187 /*-----------------------------------------------------------------------*/
188
189 /* MM on a CPU. */
190 cudaEventRecord(start,0);
191 for (int r = 0; r < N; ++r )
192 for (int c = 0; c < N; ++c )
193 for (int k = 0; k < N; ++k )
194 C[r*N+c] += A[r*N+c] * B[k*N+c];
195 cudaEventRecord(stop,0);
196 cudaEventSynchronize(stop);
197 cudaEventElapsedTime(&CPU_elapsedtime,start,stop);
198 memcpy(cpuResult, C, size_c*sizeof(double));
199 printf(" speedup\n");
200 printf(" -------\n");
201 printf("Elapsed time in CPU: %7.1f milliseconds\n", CPU_elapsedtime);
202 /*-----------------------------------------------------------------------*/
203
204 /* MM on Global Memory of GPGPU. */
205 cudaEventRecord(start,0);
206
207 /* 2) Copy operands from CPU memory to GPGPU memory. */
208 //cudaMalloc((void**)&G_A,size_a*sizeof(double)); /* alloc A in GPGPU */
209 G_A = (double*)$malloc($root, size_a*sizeof(double));
210 //cudaMalloc((void**)&G_B,size_b*sizeof(double)); /* alloc B in GPGPU */
211 G_B = (double*)$malloc($root, size_b*sizeof(double));
212 //cudaMalloc((void**)&G_C,size_c*sizeof(double)); /* alloc C in GPGPU */
213 G_C = (double*)$malloc($root, size_c*sizeof(double));
214 cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
215 cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);
216
217 /* 3) Apply matrix operation to operands on GPGPU */
218 /* There is no partial final block in this example. */
219 dim3 block = { TILE_WIDTH, TILE_WIDTH, 1 }; /* using a 2D block: 16,16,1 */
220 dim3 grid = { N/TILE_WIDTH, N/TILE_WIDTH, 1 }; /* as many 16x16-thread blocks as needed: */
221 _kernel_myMM_global(grid, block, 0, G_A, G_B, G_C, N); /* grid(16,16,1) */
222
223 /* 4) Copy result from GPGPU memory to CPU memory. */
224 cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);
225
226 /* Deallocate memory on GPGPU. */
227 cudaFree(G_A);
228 cudaFree(G_B);
229 cudaFree(G_C);
230
231 cudaEventRecord(stop,0);
232 cudaEventSynchronize(stop);
233 cudaEventElapsedTime(&GPU_global_elapsedtime,start,stop);
234 memcpy(gpuGlobalResult, C, size_c*sizeof(double));
235 printf("Elapsed time in GPU (global memory): %7.1f milliseconds %5.1f\n",
236 GPU_global_elapsedtime,CPU_elapsedtime/GPU_global_elapsedtime);
237 /*
238 printf("\nGLOBAL MEMORY:\n");
239 for (r=0;r<10;++r)
240 for (c=0;c<10;++c) {
241 printf("%2d,%2d %g\n", r,c,C[r*N+c]);
242 }
243 */
244 /*-----------------------------------------------------------------------*/
245
246 /* MM on Shared Memory of GPGPU. */
247 cudaEventRecord(start,0);
248
249 /* 2) Copy operands from CPU memory to GPGPU memory. */
250 G_A = (double*)$malloc($root, size_a*sizeof(double));
251 G_B = (double*)$malloc($root, size_b*sizeof(double));
252 G_C = (double*)$malloc($root, size_c*sizeof(double));
253 cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
254 cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);
255
256 /* 3) Apply matrix operation to operands on GPGPU */
257 /* There is not partial final block in this example. */
258 /* Use the same grid and block from the previous case. */
259 printf("a\n");
260 _kernel_myMM_shared(grid, block, 0, G_A, G_B, G_C, N);
261 printf("b\n");
262
263 /* 4) Copy result from GPGPU memory to CPU memory. */
264 cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);
265
266 /* Deallocate memory on GPGPU. */
267 cudaFree(G_A);
268 cudaFree(G_B);
269 cudaFree(G_C);
270
271 cudaEventRecord(stop,0);
272 cudaEventSynchronize(stop);
273 cudaEventElapsedTime(&GPU_shared_elapsedtime,start,stop);
274 memcpy(gpuSharedResult, C, size_c*sizeof(double));
275 printf("Elapsed time in GPU (shared memory): %7.1f milliseconds %5.1f\n",
276 GPU_shared_elapsedtime,CPU_elapsedtime/GPU_shared_elapsedtime);
277 /*
278 printf("\nSHARED MEMORY:\n");
279 for (r=0;r<10;++r)
280 for (c=0;c<10;++c) {
281 printf("%2d,%2d %g\n", r,c,C[r*N+c]);
282 }
283 */
284 /*-----------------------------------------------------------------------*/
285
286 /* Deallocate the clock. */
287 cudaEventDestroy(start);
288 cudaEventDestroy(stop);
289 //$assert($forall {int i | i >= 0 && i < N*N} cpuResult[i] == gpuGlobalResult[i] && cpuResult[i] == gpuSharedResult[i]);
290 for (int i = 0; i < N*N; i ++) {
291 printf("i = %d\n", i);
292 printf("cpu = %f, gpuGlobal = %f, gpuShared = %f\n", cpuResult[i], gpuGlobalResult[i], gpuSharedResult[i]);
293 $assert(cpuResult[i] == gpuGlobalResult[i]);
294 $assert(cpuResult[i] == gpuSharedResult[i]);
295 }
296
297 return 0;
298 }
299
300 printf("1\n");
301 _cudaInit();
302 printf("2\n");
303 _main();
304 printf("3\n");
305 _cudaFinalize();
306 printf("4\n");
307 return 0;
308}
Note: See TracBrowser for help on using the repository browser.