source: CIVL/examples/cuda/matMult1.cvl@ b2a2faa

1.23 2.0 main test-branch
Last change on this file since b2a2faa was 3ff27cf, checked in by Manchun Zheng <zmanchun@…>, 11 years ago

updated examples since $assert/$assume has been changed to functions; fixed the model builder for the new side-effect remover.

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