source: CIVL/examples/cuda/matMult1.cu@ 139c8d5

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

Got cuda-omp.cu example verified. Added matMult example from civl-papers directory. Updated CIVL2CudaTransformerTest to run matMult example. Minor changes to cuda.cvl

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