Skip to content

Commit a3b2327

Browse files
harrismharrism
committed
Added CUDA C version of matrix transpose
1 parent 9c96d68 commit a3b2327

1 file changed

Lines changed: 300 additions & 0 deletions

File tree

‎series/cuda-cpp/transpose/transpose.cu‎

Lines changed: 300 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,300 @@
1+
// Copyright 2012 NVIDIA Corporation
2+
//
3+
// Licensed under the Apache License, Version 2.0 (the "License");
4+
// you may not use this file except in compliance with the License.
5+
// You may obtain a copy of the License at
6+
//
7+
// http://www.apache.org/licenses/LICENSE-2.0
8+
//
9+
// Unless required by applicable law or agreed to in writing, software
10+
// distributed under the License is distributed on an "AS IS" BASIS,
11+
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12+
// See the License for the specific language governing permissions and
13+
// limitations under the License.
14+
15+
#include <stdio.h>
16+
#include <assert.h>
17+
18+
// Convenience function for checking CUDA runtime API results
19+
// can be wrapped around any runtime API call. No-op in release builds.
20+
inline
21+
cudaError_t checkCuda(cudaError_t result)
22+
{
23+
#if defined(DEBUG) || defined(_DEBUG)
24+
if (result != cudaSuccess) {
25+
fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
26+
assert(result == cudaSuccess);
27+
}
28+
#endif
29+
return result;
30+
}
31+
32+
const int TILE_DIM = 32;
33+
const int BLOCK_ROWS = 8;
34+
const int NUM_REPS = 100;
35+
36+
// Check errors and print GB/s
37+
void postprocess(const float *ref, const float *res, int n, float ms)
38+
{
39+
bool passed = true;
40+
for (int i = 0; i < n; i++)
41+
if (res[i] != ref[i]) {
42+
printf("%d %f %f\n", i, res[i], ref[i]);
43+
printf("%25s\n", "*** FAILED ***");
44+
passed = false;
45+
break;
46+
}
47+
if (passed)
48+
printf("%20.2f\n", 2 * n * sizeof(float) * 1e-6 * NUM_REPS / ms );
49+
}
50+
51+
// simple copy kernel
52+
// Used as reference case representing best effective bandwidth.
53+
__global__ void copy(float *odata, const float *idata)
54+
{
55+
int x = blockIdx.x * TILE_DIM + threadIdx.x;
56+
int y = blockIdx.y * TILE_DIM + threadIdx.y;
57+
int width = gridDim.x * TILE_DIM;
58+
59+
for (int j = 0; j < TILE_DIM; j+= BLOCK_ROWS)
60+
odata[(y+j)*width + x] = idata[(y+j)*width + x];
61+
}
62+
63+
// copy kernel using shared memory
64+
// Also used as reference case, demonstrating effect of using shared memory.
65+
__global__ void copySharedMem(float *odata, const float *idata)
66+
{
67+
__shared__ float tile[TILE_DIM * TILE_DIM];
68+
69+
int x = blockIdx.x * TILE_DIM + threadIdx.x;
70+
int y = blockIdx.y * TILE_DIM + threadIdx.y;
71+
int width = gridDim.x * TILE_DIM;
72+
73+
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
74+
tile[(threadIdx.y+j)*TILE_DIM + threadIdx.x] = idata[(y+j)*width + x];
75+
76+
__syncthreads();
77+
78+
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
79+
odata[(y+j)*width + x] = tile[(threadIdx.y+j)*TILE_DIM + threadIdx.x];
80+
}
81+
82+
// naive transpose
83+
// Simplest transpose; doesn't use shared memory.
84+
// Global memory reads are coalesced but writes are not.
85+
__global__ void transposeNaive(float *odata, const float *idata)
86+
{
87+
int x = blockIdx.x * TILE_DIM + threadIdx.x;
88+
int y = blockIdx.y * TILE_DIM + threadIdx.y;
89+
int width = gridDim.x * TILE_DIM;
90+
91+
for (int j = 0; j < TILE_DIM; j+= BLOCK_ROWS)
92+
odata[x*width + (y+j)] = idata[(y+j)*width + x];
93+
}
94+
95+
// coalesced transpose
96+
// Uses shared memory to achieve coalesing in both reads and writes
97+
// Tile width == #banks causes shared memory bank conflicts.
98+
__global__ void transposeCoalesced(float *odata, const float *idata)
99+
{
100+
__shared__ float tile[TILE_DIM][TILE_DIM];
101+
102+
int x = blockIdx.x * TILE_DIM + threadIdx.x;
103+
int y = blockIdx.y * TILE_DIM + threadIdx.y;
104+
int width = gridDim.x * TILE_DIM;
105+
106+
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
107+
tile[threadIdx.y+j][threadIdx.x] = idata[(y+j)*width + x];
108+
109+
__syncthreads();
110+
111+
x = blockIdx.y * TILE_DIM + threadIdx.x; // transpose block offset
112+
y = blockIdx.x * TILE_DIM + threadIdx.y;
113+
114+
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
115+
odata[(y+j)*width + x] = tile[threadIdx.x][threadIdx.y + j];
116+
}
117+
118+
119+
// No bank-conflict transpose
120+
// Same as transposeCoalesced except the first tile dimension is padded
121+
// to avoid shared memory bank conflicts.
122+
__global__ void transposeNoBankConflicts(float *odata, const float *idata)
123+
{
124+
__shared__ float tile[TILE_DIM][TILE_DIM+1];
125+
126+
int x = blockIdx.x * TILE_DIM + threadIdx.x;
127+
int y = blockIdx.y * TILE_DIM + threadIdx.y;
128+
int width = gridDim.x * TILE_DIM;
129+
130+
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
131+
tile[threadIdx.y+j][threadIdx.x] = idata[(y+j)*width + x];
132+
133+
__syncthreads();
134+
135+
x = blockIdx.y * TILE_DIM + threadIdx.x; // transpose block offset
136+
y = blockIdx.x * TILE_DIM + threadIdx.y;
137+
138+
for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
139+
odata[(y+j)*width + x] = tile[threadIdx.x][threadIdx.y + j];
140+
}
141+
142+
int main(int argc, char **argv)
143+
{
144+
const int nx = 1024;
145+
const int ny = 1024;
146+
const int mem_size = nx*ny*sizeof(float);
147+
148+
dim3 dimGrid(nx/TILE_DIM, ny/TILE_DIM, 1);
149+
dim3 dimBlock(TILE_DIM, BLOCK_ROWS, 1);
150+
151+
int devId = 0;
152+
if (argc > 1) devId = atoi(argv[1]);
153+
154+
cudaDeviceProp prop;
155+
checkCuda( cudaGetDeviceProperties(&prop, devId));
156+
printf("\nDevice : %s\n", prop.name);
157+
printf("Matrix size: %d %d, Block size: %d %d, Tile size: %d %d\n",
158+
nx, ny, TILE_DIM, BLOCK_ROWS, TILE_DIM, TILE_DIM);
159+
printf("dimGrid: %d %d %d. dimBlock: %d %d %d\n",
160+
dimGrid.x, dimGrid.y, dimGrid.z, dimBlock.x, dimBlock.y, dimBlock.z);
161+
162+
checkCuda( cudaSetDevice(devId) );
163+
164+
float *h_idata = (float*)malloc(mem_size);
165+
float *h_cdata = (float*)malloc(mem_size);
166+
float *h_tdata = (float*)malloc(mem_size);
167+
float *gold = (float*)malloc(mem_size);
168+
169+
float *d_idata, *d_cdata, *d_tdata;
170+
checkCuda( cudaMalloc(&d_idata, mem_size) );
171+
checkCuda( cudaMalloc(&d_cdata, mem_size) );
172+
checkCuda( cudaMalloc(&d_tdata, mem_size) );
173+
174+
// check parameters and calculate execution configuration
175+
if (nx % TILE_DIM || ny % TILE_DIM) {
176+
printf("nx and ny must be a multiple of TILE_DIM\n");
177+
goto error_exit;
178+
}
179+
180+
if (TILE_DIM % BLOCK_ROWS) {
181+
printf("TILE_DIM must be a multiple of BLOCK_ROWS\n");
182+
goto error_exit;
183+
}
184+
185+
// host
186+
for (int j = 0; j < ny; j++)
187+
for (int i = 0; i < nx; i++)
188+
h_idata[j*nx + i] = j*nx + i;
189+
190+
// correct result for error checking
191+
for (int j = 0; j < ny; j++)
192+
for (int i = 0; i < nx; i++)
193+
gold[j*nx + i] = h_idata[i*nx + j];
194+
195+
// device
196+
checkCuda( cudaMemcpy(d_idata, h_idata, mem_size, cudaMemcpyHostToDevice) );
197+
198+
// events for timing
199+
cudaEvent_t startEvent, stopEvent;
200+
checkCuda( cudaEventCreate(&startEvent) );
201+
checkCuda( cudaEventCreate(&stopEvent) );
202+
float ms;
203+
204+
// ------------
205+
// time kernels
206+
// ------------
207+
printf("%25s%25s\n", "Routine", "Bandwidth (GB/s)");
208+
209+
// ----
210+
// copy
211+
// ----
212+
printf("%25s", "copy");
213+
checkCuda( cudaMemset(d_cdata, 0, mem_size) );
214+
// warm up
215+
copy<<<dimGrid, dimBlock>>>(d_cdata, d_idata);
216+
checkCuda( cudaEventRecord(startEvent, 0) );
217+
for (int i = 0; i < NUM_REPS; i++)
218+
copy<<<dimGrid, dimBlock>>>(d_cdata, d_idata);
219+
checkCuda( cudaEventRecord(stopEvent, 0) );
220+
checkCuda( cudaEventSynchronize(stopEvent) );
221+
checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
222+
checkCuda( cudaMemcpy(h_cdata, d_cdata, mem_size, cudaMemcpyDeviceToHost) );
223+
postprocess(h_idata, h_cdata, nx*ny, ms);
224+
225+
// -------------
226+
// copySharedMem
227+
// -------------
228+
printf("%25s", "shared memory copy");
229+
checkCuda( cudaMemset(d_cdata, 0, mem_size) );
230+
// warm up
231+
copySharedMem<<<dimGrid, dimBlock>>>(d_cdata, d_idata);
232+
checkCuda( cudaEventRecord(startEvent, 0) );
233+
for (int i = 0; i < NUM_REPS; i++)
234+
copySharedMem<<<dimGrid, dimBlock>>>(d_cdata, d_idata);
235+
checkCuda( cudaEventRecord(stopEvent, 0) );
236+
checkCuda( cudaEventSynchronize(stopEvent) );
237+
checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
238+
checkCuda( cudaMemcpy(h_cdata, d_cdata, mem_size, cudaMemcpyDeviceToHost) );
239+
postprocess(h_idata, h_cdata, nx * ny, ms);
240+
241+
// --------------
242+
// transposeNaive
243+
// --------------
244+
printf("%25s", "naive transpose");
245+
checkCuda( cudaMemset(d_tdata, 0, mem_size) );
246+
// warmup
247+
transposeNaive<<<dimGrid, dimBlock>>>(d_tdata, d_idata);
248+
checkCuda( cudaEventRecord(startEvent, 0) );
249+
for (int i = 0; i < NUM_REPS; i++)
250+
transposeNaive<<<dimGrid, dimBlock>>>(d_tdata, d_idata);
251+
checkCuda( cudaEventRecord(stopEvent, 0) );
252+
checkCuda( cudaEventSynchronize(stopEvent) );
253+
checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
254+
checkCuda( cudaMemcpy(h_tdata, d_tdata, mem_size, cudaMemcpyDeviceToHost) );
255+
postprocess(gold, h_tdata, nx * ny, ms);
256+
257+
// ------------------
258+
// transposeCoalesced
259+
// ------------------
260+
printf("%25s", "coalesced transpose");
261+
checkCuda( cudaMemset(d_tdata, 0, mem_size) );
262+
// warmup
263+
transposeCoalesced<<<dimGrid, dimBlock>>>(d_tdata, d_idata);
264+
checkCuda( cudaEventRecord(startEvent, 0) );
265+
for (int i = 0; i < NUM_REPS; i++)
266+
transposeCoalesced<<<dimGrid, dimBlock>>>(d_tdata, d_idata);
267+
checkCuda( cudaEventRecord(stopEvent, 0) );
268+
checkCuda( cudaEventSynchronize(stopEvent) );
269+
checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
270+
checkCuda( cudaMemcpy(h_tdata, d_tdata, mem_size, cudaMemcpyDeviceToHost) );
271+
postprocess(gold, h_tdata, nx * ny, ms);
272+
273+
// ------------------------
274+
// transposeNoBankConflicts
275+
// ------------------------
276+
printf("%25s", "conflict-free transpose");
277+
checkCuda( cudaMemset(d_tdata, 0, mem_size) );
278+
// warmup
279+
transposeNoBankConflicts<<<dimGrid, dimBlock>>>(d_tdata, d_idata);
280+
checkCuda( cudaEventRecord(startEvent, 0) );
281+
for (int i = 0; i < NUM_REPS; i++)
282+
transposeNoBankConflicts<<<dimGrid, dimBlock>>>(d_tdata, d_idata);
283+
checkCuda( cudaEventRecord(stopEvent, 0) );
284+
checkCuda( cudaEventSynchronize(stopEvent) );
285+
checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
286+
checkCuda( cudaMemcpy(h_tdata, d_tdata, mem_size, cudaMemcpyDeviceToHost) );
287+
postprocess(gold, h_tdata, nx * ny, ms);
288+
289+
error_exit:
290+
// cleanup
291+
checkCuda( cudaEventDestroy(startEvent) );
292+
checkCuda( cudaEventDestroy(stopEvent) );
293+
checkCuda( cudaFree(d_tdata) );
294+
checkCuda( cudaFree(d_cdata) );
295+
checkCuda( cudaFree(d_idata) );
296+
free(h_idata);
297+
free(h_tdata);
298+
free(h_cdata);
299+
free(gold);
300+
}

0 commit comments

Comments
 (0)