source: CIVL/examples/cuda/matMult1.cu@ b9a6041

1.23 2.0 main test-branch
Last change on this file since b9a6041 was 47665c1, checked in by Andre Marianiello <andre.marianiello@…>, 11 years ago

Updated cuda.cvl and civl-cuda.cvl implementations to hide some struct definitions

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@2211 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 9.0 KB
Line 
1// TODO add assertions to check equivalence of implementation results
2#ifdef _CIVL
3#include <civlc.cvh>
4#endif
5/***********************************************************************
6* FILENAME: MM.cu
7* Matrix Multiplication
8* Matrix operands have row-major order.
9*
10* C = A * B
11* Multiplies two square matrices (NxN * NxN).
12* Matrix values have type double.
13*
14* A simple CUDA program has a basic workflow:
15* 1) Initialize matrix operands as double-precision arrays on host (CPU).
16* 2) Copy operands from host memory to GPU memory.
17* 3) Apply matrix operaton to operands on GPU
18* 4) Copy result from GPU memory to host memory.
19*
20*
21* CUDA C Programming Guide Version 4.2 (3.2.3, p.22):
22* http://developer.download.nvidia.com/compute/DevZone/docs/html/C/doc/CUDA_C_Programming_Guide.pdf
23*
24* MM with linearized matrix operands:
25* http://www.hpcwire.com/hpcwire/2008-10-08/compilers_and_more_programming_gpus_today.html
26*
27*************************************************************************/
28// online source: https://www.rcac.purdue.edu/userinfo/resources/carter/compile/MM.cu
29
30#include <stdio.h>
31#include <stdlib.h>
32#include "cuda.h"
33
34#ifdef _CIVL
35$input int N;
36$input int TILE_WIDTH;
37#else
38#define N 1024 /* size of square matrix */
39#define TILE_WIDTH 16
40#endif
41
42
43/* MM kernel using global (not shared) memory. */
44__global__
45void myMM_global (const double * const A, const double * const B, double *C, int width) {
46
47 /* Get row and column from block and thread IDs */
48 int row = (blockDim.y*blockIdx.y) + threadIdx.y;
49 int col = (blockDim.x*blockIdx.x) + threadIdx.x;
50
51 /* Initialize result of one element which one thread computes. */
52 double result=0.0;
53
54 /* Compute one element of the matrix product. */
55 for (int i = 0; i < width; ++i)
56 result += A[row*width + i] * B[i*width + col];
57
58 /* Store the result of one matrix element in matrix C. */
59 C[row * width + col] = result;
60}
61
62
63/* MM kernel using shared memory. */
64__global__
65void myMM_shared (const double * const A, const double * const B, double* C, int width) {
66 __shared__ double A_shared[TILE_WIDTH][TILE_WIDTH];
67 __shared__ double B_shared[TILE_WIDTH][TILE_WIDTH];
68
69 int bx = blockIdx.x; int by = blockIdx.y;
70 int tx = threadIdx.x; int ty = threadIdx.y;
71
72 /* Identify the row and column of the C element to work on. */
73 int row = by * TILE_WIDTH + ty;
74 int col = bx * TILE_WIDTH + tx;
75
76 double result = 0.0;
77
78 /* Loop over the A and B tiles required to compute the C element. */
79 for (int phase = 0; phase < width/TILE_WIDTH; ++phase) {
80 /* Shared effort: loading of A and B tiles into shared memory. */
81 A_shared[ty][tx] = A[row*width + (phase*TILE_WIDTH + tx)];
82 B_shared[ty][tx] = B[col + (phase*TILE_WIDTH + ty)*width];
83 __syncthreads();
84
85 for (int k = 0; k < TILE_WIDTH; ++k)
86 result += A_shared[ty][k] * B_shared[k][tx];
87 __syncthreads();
88
89 }
90 C[row*width+col] = result;
91}
92
93
94/************************************************************************/
95/************************************************************************/
96/************************************************************************/
97
98
99int main (int argc, char** argv) {
100
101#ifdef _CIVL
102 $assume(argc == 2);
103 $assume(atoi(argv[1]) == 0);
104#endif
105
106 /* Set device based on input from command line */
107 if (argc > 1) {
108 if (cudaSetDevice(atoi(argv[1])) != cudaSuccess) {
109 int num_devices;
110 cudaGetDeviceCount(&num_devices);
111 fprintf(stderr, "Error initializing device %s,\
112 device value must be 0-%d\n", argv[1], (num_devices-1));
113 return 0;
114 }
115 } else {
116 fprintf(stderr, "Usage: %s gpu_device\n", argv[0]);
117 return 0;
118 }
119
120 /* Declare CPU arrays. */
121#ifdef _CIVL
122 $input double A[N*N];
123 $input double B[N*N];
124#else
125 double A[N*N],B[N*N];
126#endif
127 double C[N*N]; /* linearized CPU double arrays */
128 int r,c;
129
130 /* Declare GPU arrays. */
131 double *G_A,*G_B,*G_C; /* linearized GPU double arrays */
132 size_t size_a,size_b,size_c; /* size of linearized array in bytes */
133 size_a = size_b = size_c = N*N;
134
135 /* Setup a clock. */
136 cudaEvent_t start, stop;
137 float CPU_elapsedtime, GPU_global_elapsedtime, GPU_shared_elapsedtime;
138 cudaEventCreate(&start);
139 cudaEventCreate(&stop);
140
141
142
143
144 /* 1) Initialize matrix operands as double-precision arrays on host (CPU). */
145#ifndef _CIVL
146 for (r=0;r<N;++r)
147 for (c=0;c<N;++c) {
148 A[r*N+c] = 1.0;
149 B[r*N+c] = 1.0;
150 }
151#endif
152
153
154/*-----------------------------------------------------------------------*/
155
156 /* MM on a CPU. */
157 cudaEventRecord(start,0);
158 for (int r = 0; r < N; ++r )
159 for (int c = 0; c < N; ++c )
160 for (int k = 0; k < N; ++k )
161 C[r*N+c] += A[r*N+c] * B[k*N+c];
162 cudaEventRecord(stop,0);
163 cudaEventSynchronize(stop);
164 cudaEventElapsedTime(&CPU_elapsedtime,start,stop);
165 printf(" speedup\n");
166 printf(" -------\n");
167 printf("Elapsed time in CPU: %7.1f milliseconds\n", CPU_elapsedtime);
168/*-----------------------------------------------------------------------*/
169
170 /* MM on Global Memory of GPGPU. */
171 cudaEventRecord(start,0);
172
173 /* 2) Copy operands from CPU memory to GPGPU memory. */
174 cudaMalloc((void**)&G_A,size_a*sizeof(double)); /* alloc A in GPGPU */
175 cudaMalloc((void**)&G_B,size_b*sizeof(double)); /* alloc B in GPGPU */
176 cudaMalloc((void**)&G_C,size_c*sizeof(double)); /* alloc C in GPGPU */
177 cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
178 cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);
179
180 /* 3) Apply matrix operation to operands on GPGPU */
181 /* There is no partial final block in this example. */
182 dim3 block = {TILE_WIDTH,TILE_WIDTH,1}; /* using a 2D block: 16,16,1 */
183 dim3 grid = {N/TILE_WIDTH,N/TILE_WIDTH,1}; /* as many 16x16-thread blocks as needed: */
184 myMM_global<<< grid,block >>>(G_A,G_B,G_C,N); /* grid(16,16,1) */
185
186 /* 4) Copy result from GPGPU memory to CPU memory. */
187 cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);
188
189 /* Deallocate memory on GPGPU. */
190 cudaFree(G_A);
191 cudaFree(G_B);
192 cudaFree(G_C);
193
194 cudaEventRecord(stop,0);
195 cudaEventSynchronize(stop);
196 cudaEventElapsedTime(&GPU_global_elapsedtime,start,stop);
197 printf("Elapsed time in GPU (global memory): %7.1f milliseconds %5.1f\n",
198 GPU_global_elapsedtime,CPU_elapsedtime/GPU_global_elapsedtime);
199//*
200 printf("\nGLOBAL MEMORY:\n");
201 for (r=0;r<N;++r)
202 for (c=0;c<N;++c) {
203 printf("%2d,%2d %g\n", r,c,C[r*N+c]);
204 }
205//*/
206/*-----------------------------------------------------------------------*/
207
208 /* MM on Shared Memory of GPGPU. */
209 cudaEventRecord(start,0);
210
211 /* 2) Copy operands from CPU memory to GPGPU memory. */
212 cudaMalloc((void**)&G_A,size_a*sizeof(double)); /* alloc A in GPGPU */
213 cudaMalloc((void**)&G_B,size_b*sizeof(double)); /* alloc B in GPGPU */
214 cudaMalloc((void**)&G_C,size_c*sizeof(double)); /* alloc C in GPGPU */
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 not partial final block in this example. */
220 /* Use the same grid and block from the previous case. */
221 myMM_shared<<< grid,block >>>(G_A,G_B,G_C,N);
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_shared_elapsedtime,start,stop);
234 printf("Elapsed time in GPU (shared memory): %7.1f milliseconds %5.1f\n",
235 GPU_shared_elapsedtime,CPU_elapsedtime/GPU_shared_elapsedtime);
236//*
237 printf("\nSHARED MEMORY:\n");
238 for (r=0;r<N;++r)
239 for (c=0;c<N;++c) {
240 printf("%2d,%2d %g\n", r,c,C[r*N+c]);
241 }
242//*/
243/*-----------------------------------------------------------------------*/
244
245 /* Deallocate the clock. */
246 cudaEventDestroy(start);
247 cudaEventDestroy(stop);
248
249 return 0;
250}
251
Note: See TracBrowser for help on using the repository browser.