source: CIVL/examples/translation/cuda/matMult1.cvl@ 764d472

1.23 2.0 main test-branch
Last change on this file since 764d472 was a86ac45, checked in by Andre Marianiello <andre.marianiello@…>, 12 years ago

Moved a few more files over from the civl-papers/examples repo

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@1187 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 9.0 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
30#define N 1024 /* size of square matrix */
31#define TILE_WIDTH 16
32
33
34/* MM kernel using global (not shared) memory. */
35void _kernel_myMM_global (dim3 gridDim, dim3 blockDim, const double * const A, const double * const B, double *C, int width) {
36
37 void _block(uint3 blockIdx) {
38
39 void _thread(uint3 threadIdx) {
40
41 /* Get row and column from block and thread IDs */
42 int row = (blockDim.y*blockIdx.y) + threadIdx.y;
43 int col = (blockDim.x*blockIdx.x) + threadIdx.x;
44
45 /* Initialize result of one element which one thread computes. */
46 double result=0.0;
47
48 /* Compute one element of the matrix product. */
49 for (int i = 0; i < width; ++i)
50 result += A[row*width + i] * B[i*width + col];
51
52 /* Store the result of one matrix element in matrix C. */
53 C[row * width + col] = result;
54 }
55 $gbarrier_destroy(_block_barrier);
56 }
57 _createProcs(gridDim, _block);
58}
59
60
61/* MM kernel using shared memory. */
62void _kernel_myMM_shared (dim3 gridDim, dim3 blockDim, const double * const A, const double * const B, double* C, int width) {
63
64 void _block(uint3 blockIdx) {
65 int _numThreads = blockDim.x * blockDim.y * blockDim.z;
66 $gbarrier _block_barrier = $gbarrier_create($here, _numThreads);
67
68 double A_shared[TILE_WIDTH][TILE_WIDTH];
69 double B_shared[TILE_WIDTH][TILE_WIDTH];
70
71 void _thread(uint3 threadIdx) {
72
73 int _tid = _index(blockDim, threadIdx);
74 $barrier _b = $barrier_create($here, _block_barrier, _tid);
75
76 int bx = blockIdx.x; int by = blockIdx.y;
77 int tx = threadIdx.x; int ty = threadIdx.y;
78
79 /* Identify the row and column of the C element to work on. */
80 int row = by * TILE_WIDTH + ty;
81 int col = bx * TILE_WIDTH + tx;
82
83 double result = 0.0;
84
85 /* Loop over the A and B tiles required to compute the C element. */
86 for (int phase = 0; phase < width/TILE_WIDTH; ++phase) {
87 /* Shared effort: loading of A and B tiles into shared memory. */
88
89 A_shared[ty][tx] = A[row*width + (phase*TILE_WIDTH + tx)];
90
91 B_shared[ty][tx] = B[col + (phase*TILE_WIDTH + ty)*width];
92
93 $barrier_call(_b);
94
95 for (int k = 0; k < TILE_WIDTH; ++k) {
96 result += A_shared[ty][k] * B_shared[k][tx];
97 }
98
99 $barrier_call(_b);
100
101 }
102 C[row*width+col] = result;
103 }
104 _createProcs(blockDim, _thread);
105 $gbarrier_destroy(_block_barrier);
106 }
107 _createProcs(gridDim, _block);
108}
109
110
111/************************************************************************/
112/************************************************************************/
113/************************************************************************/
114
115
116int main (int argc, char** argv) {
117 _host_scope = $here;
118
119 int _main ( void ) {
120
121 /* Set device based on input from command line */
122 if (argc > 1) {
123 if (cudaSetDevice(atoi(argv[1])) != cudaSuccess) {
124 int num_devices;
125 cudaGetDeviceCount(&num_devices);
126 fprintf(stderr, "Error initializing device %s,\
127 device value must be 0-%d\n", argv[1], (num_devices-1));
128 return 0;
129 }
130 } else {
131 fprintf(stderr, "Usage: %s gpu_device\n", argv[0]);
132 return 0;
133 }
134
135 /* Declare CPU arrays. */
136 double A[N*N],B[N*N],C[N*N]; /* linearized CPU double arrays */
137 int r,c;
138
139 /* Declare GPU arrays. */
140 double *G_A,*G_B,*G_C; /* linearized GPU double arrays */
141 size_t size_a,size_b,size_c; /* size of linearized array in bytes */
142 size_a = size_b = size_c = N*N;
143
144 /* Setup a clock. */
145 cudaEvent_t start, stop;
146 float CPU_elapsedtime, GPU_global_elapsedtime, GPU_shared_elapsedtime;
147 cudaEventCreate(&start);
148 cudaEventCreate(&stop);
149
150
151
152
153 /* 1) Initialize matrix operands as double-precision arrays on host (CPU). */
154 for (r=0;r<N;++r)
155 for (c=0;c<N;++c) {
156 A[r*N+c] = 1.0;
157 B[r*N+c] = 1.0;
158 }
159
160
161 /*-----------------------------------------------------------------------*/
162
163 /* MM on a CPU. */
164 cudaEventRecord(start,0);
165 for (int r = 0; r < N; ++r )
166 for (int c = 0; c < N; ++c )
167 for (int k = 0; k < N; ++k )
168 C[r*N+c] += A[r*N+c] * B[k*N+c];
169 cudaEventRecord(stop,0);
170 cudaEventSynchronize(stop);
171 cudaEventElapsedTime(&CPU_elapsedtime,start,stop);
172 printf(" speedup\n");
173 printf(" -------\n");
174 printf("Elapsed time in CPU: %7.1f milliseconds\n", CPU_elapsedtime);
175 /*-----------------------------------------------------------------------*/
176
177 /* MM on Global Memory of GPGPU. */
178 cudaEventRecord(start,0);
179
180 /* 2) Copy operands from CPU memory to GPGPU memory. */
181 cudaMalloc((void**)&G_A,size_a*sizeof(double)); /* alloc A in GPGPU */
182 cudaMalloc((void**)&G_B,size_b*sizeof(double)); /* alloc B in GPGPU */
183 cudaMalloc((void**)&G_C,size_c*sizeof(double)); /* alloc C in GPGPU */
184 cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
185 cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);
186
187 /* 3) Apply matrix operation to operands on GPGPU */
188 /* There is no partial final block in this example. */
189 dim3 block(TILE_WIDTH,TILE_WIDTH); /* using a 2D block: 16,16,1 */
190 dim3 grid(N/TILE_WIDTH,N/TILE_WIDTH); /* as many 16x16-thread blocks as needed: */
191 myMM_global<<< grid,block >>>(G_A,G_B,G_C,N); /* grid(16,16,1) */
192
193 /* 4) Copy result from GPGPU memory to CPU memory. */
194 cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);
195
196 /* Deallocate memory on GPGPU. */
197 cudaFree(G_A);
198 cudaFree(G_B);
199 cudaFree(G_C);
200
201 cudaEventRecord(stop,0);
202 cudaEventSynchronize(stop);
203 cudaEventElapsedTime(&GPU_global_elapsedtime,start,stop);
204 printf("Elapsed time in GPU (global memory): %7.1f milliseconds %5.1f\n",
205 GPU_global_elapsedtime,CPU_elapsedtime/GPU_global_elapsedtime);
206 /*
207 printf("\nGLOBAL MEMORY:\n");
208 for (r=0;r<10;++r)
209 for (c=0;c<10;++c) {
210 printf("%2d,%2d %g\n", r,c,C[r*N+c]);
211 }
212 */
213 /*-----------------------------------------------------------------------*/
214
215 /* MM on Shared Memory of GPGPU. */
216 cudaEventRecord(start,0);
217
218 /* 2) Copy operands from CPU memory to GPGPU memory. */
219 cudaMalloc((void**)&G_A,size_a*sizeof(double)); /* alloc A in GPGPU */
220 cudaMalloc((void**)&G_B,size_b*sizeof(double)); /* alloc B in GPGPU */
221 cudaMalloc((void**)&G_C,size_c*sizeof(double)); /* alloc C in GPGPU */
222 cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
223 cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);
224
225 /* 3) Apply matrix operation to operands on GPGPU */
226 /* There is not partial final block in this example. */
227 /* Use the same grid and block from the previous case. */
228 myMM_shared<<< grid,block >>>(G_A,G_B,G_C,N);
229
230 /* 4) Copy result from GPGPU memory to CPU memory. */
231 cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);
232
233 /* Deallocate memory on GPGPU. */
234 cudaFree(G_A);
235 cudaFree(G_B);
236 cudaFree(G_C);
237
238 cudaEventRecord(stop,0);
239 cudaEventSynchronize(stop);
240 cudaEventElapsedTime(&GPU_shared_elapsedtime,start,stop);
241 printf("Elapsed time in GPU (shared memory): %7.1f milliseconds %5.1f\n",
242 GPU_shared_elapsedtime,CPU_elapsedtime/GPU_shared_elapsedtime);
243 /*
244 printf("\nSHARED MEMORY:\n");
245 for (r=0;r<10;++r)
246 for (c=0;c<10;++c) {
247 printf("%2d,%2d %g\n", r,c,C[r*N+c]);
248 }
249 */
250 /*-----------------------------------------------------------------------*/
251
252 /* Deallocate the clock. */
253 cudaEventDestroy(start);
254 cudaEventDestroy(stop);
255
256 return 0;
257 }
258
259 _main();
260}
Note: See TracBrowser for help on using the repository browser.