|
1 | 1 | #pragma once |
2 | 2 |
|
3 | | -#include <cuda_runtime.h> |
4 | 3 | #include "../utils/cuda_utils.cuh" |
| 4 | +#include <cuda_runtime.h> |
5 | 5 |
|
6 | 6 | /** |
7 | 7 | * Bank Conflict Free SGEMM Kernel |
8 | | - * |
| 8 | + * |
9 | 9 | * This implementation eliminates shared memory bank conflicts by adding padding |
10 | 10 | * to the shared memory arrays. |
11 | | - * |
| 11 | + * |
12 | 12 | * ============================================================================ |
13 | 13 | * Bank Conflict Explanation: |
14 | 14 | * ============================================================================ |
15 | | - * |
| 15 | + * |
16 | 16 | * Shared memory is divided into 32 banks (on modern GPUs). |
17 | 17 | * Each bank is 4 bytes wide (one float). |
18 | 18 | * Consecutive 4-byte words map to consecutive banks. |
19 | | - * |
| 19 | + * |
20 | 20 | * Bank assignment: bank_id = (address / 4) % 32 |
21 | | - * |
| 21 | + * |
22 | 22 | * For a 32x32 array stored row-major: |
23 | 23 | * - Element [i][j] is at address: (i * 32 + j) * 4 |
24 | 24 | * - Bank: (i * 32 + j) % 32 = j (since 32 % 32 = 0) |
25 | | - * |
| 25 | + * |
26 | 26 | * Problem: When threads in a warp access column j of the array, |
27 | 27 | * they all access the same bank j, causing a 32-way bank conflict! |
28 | | - * |
| 28 | + * |
29 | 29 | * Solution: Add 1 element of padding per row. |
30 | 30 | * - Element [i][j] is now at address: (i * 33 + j) * 4 |
31 | 31 | * - Bank: (i * 33 + j) % 32 = (i + j) % 32 |
32 | 32 | * - Now threads accessing column j get different banks! |
33 | | - * |
| 33 | + * |
34 | 34 | * ============================================================================ |
35 | | - * |
| 35 | + * |
36 | 36 | * C = A * B |
37 | 37 | * A: M x K (row-major) |
38 | 38 | * B: K x N (row-major) |
39 | 39 | * C: M x N (row-major) |
40 | 40 | */ |
41 | | -template<int TILE_SIZE> |
42 | | -__global__ void bank_conflict_free_sgemm_kernel( |
43 | | - const float* __restrict__ A, |
44 | | - const float* __restrict__ B, |
45 | | - float* __restrict__ C, |
46 | | - int M, int K, int N |
47 | | -) { |
48 | | - // Shared memory with padding to avoid bank conflicts |
49 | | - // Adding 1 to the second dimension shifts each row by 1 bank |
50 | | - // This ensures column accesses hit different banks |
51 | | - __shared__ float As[TILE_SIZE][TILE_SIZE + 1]; // +1 padding |
52 | | - __shared__ float Bs[TILE_SIZE][TILE_SIZE + 1]; // +1 padding |
53 | | - |
54 | | - int bx = blockIdx.x; |
55 | | - int by = blockIdx.y; |
56 | | - int tx = threadIdx.x; |
57 | | - int ty = threadIdx.y; |
58 | | - |
59 | | - int row = by * TILE_SIZE + ty; |
60 | | - int col = bx * TILE_SIZE + tx; |
61 | | - |
62 | | - float sum = 0.0f; |
63 | | - int numTiles = (K + TILE_SIZE - 1) / TILE_SIZE; |
64 | | - |
65 | | - for (int t = 0; t < numTiles; ++t) { |
66 | | - // Load tile of A into shared memory (coalesced access) |
67 | | - int aCol = t * TILE_SIZE + tx; |
68 | | - if (row < M && aCol < K) { |
69 | | - As[ty][tx] = A[row * K + aCol]; |
70 | | - } else { |
71 | | - As[ty][tx] = 0.0f; |
72 | | - } |
73 | | - |
74 | | - // Load tile of B into shared memory (coalesced access) |
75 | | - int bRow = t * TILE_SIZE + ty; |
76 | | - if (bRow < K && col < N) { |
77 | | - Bs[ty][tx] = B[bRow * N + col]; |
78 | | - } else { |
79 | | - Bs[ty][tx] = 0.0f; |
80 | | - } |
81 | | - |
82 | | - __syncthreads(); |
83 | | - |
84 | | - // Compute partial dot product |
85 | | - // Access pattern: As[ty][k] - row access (no conflict) |
86 | | - // Bs[k][tx] - column access (no conflict due to padding!) |
87 | | - #pragma unroll |
88 | | - for (int k = 0; k < TILE_SIZE; ++k) { |
89 | | - sum += As[ty][k] * Bs[k][tx]; |
90 | | - } |
91 | | - |
92 | | - __syncthreads(); |
| 41 | +template <int TILE_SIZE> |
| 42 | +__global__ void bank_conflict_free_sgemm_kernel(const float *__restrict__ A, |
| 43 | + const float *__restrict__ B, |
| 44 | + float *__restrict__ C, int M, |
| 45 | + int K, int N) { |
| 46 | + // Shared memory with padding to avoid bank conflicts |
| 47 | + // Adding 1 to the second dimension shifts each row by 1 bank |
| 48 | + // This ensures column accesses hit different banks |
| 49 | + __shared__ float As[TILE_SIZE][TILE_SIZE + 1]; // +1 padding |
| 50 | + __shared__ float Bs[TILE_SIZE][TILE_SIZE + 1]; // +1 padding |
| 51 | + |
| 52 | + int bx = blockIdx.x; |
| 53 | + int by = blockIdx.y; |
| 54 | + int tx = threadIdx.x; |
| 55 | + int ty = threadIdx.y; |
| 56 | + |
| 57 | + int row = by * TILE_SIZE + ty; |
| 58 | + int col = bx * TILE_SIZE + tx; |
| 59 | + |
| 60 | + float sum = 0.0f; |
| 61 | + int numTiles = (K + TILE_SIZE - 1) / TILE_SIZE; |
| 62 | + |
| 63 | + for (int t = 0; t < numTiles; ++t) { |
| 64 | + // Load tile of A into shared memory (coalesced access) |
| 65 | + int aCol = t * TILE_SIZE + tx; |
| 66 | + if (row < M && aCol < K) { |
| 67 | + As[ty][tx] = A[row * K + aCol]; |
| 68 | + } else { |
| 69 | + As[ty][tx] = 0.0f; |
93 | 70 | } |
94 | | - |
95 | | - if (row < M && col < N) { |
96 | | - C[row * N + col] = sum; |
| 71 | + |
| 72 | + // Load tile of B into shared memory (coalesced access) |
| 73 | + int bRow = t * TILE_SIZE + ty; |
| 74 | + if (bRow < K && col < N) { |
| 75 | + Bs[ty][tx] = B[bRow * N + col]; |
| 76 | + } else { |
| 77 | + Bs[ty][tx] = 0.0f; |
97 | 78 | } |
| 79 | + |
| 80 | + __syncthreads(); |
| 81 | + |
| 82 | +// Compute partial dot product |
| 83 | +// Access pattern: As[ty][k] - row access (no conflict) |
| 84 | +// Bs[k][tx] - column access (no conflict due to padding!) |
| 85 | +#pragma unroll |
| 86 | + for (int k = 0; k < TILE_SIZE; ++k) { |
| 87 | + sum += As[ty][k] * Bs[k][tx]; |
| 88 | + } |
| 89 | + |
| 90 | + __syncthreads(); |
| 91 | + } |
| 92 | + |
| 93 | + if (row < M && col < N) { |
| 94 | + C[row * N + col] = sum; |
| 95 | + } |
98 | 96 | } |
99 | 97 |
|
100 | 98 | /** |
101 | 99 | * Launch wrapper for bank conflict free SGEMM kernel |
102 | 100 | */ |
103 | | -template<int TILE_SIZE = 32> |
104 | | -void launch_bank_conflict_free_sgemm( |
105 | | - const float* A, |
106 | | - const float* B, |
107 | | - float* C, |
108 | | - int M, int K, int N, |
109 | | - cudaStream_t stream = 0 |
110 | | -) { |
111 | | - dim3 blockDim(TILE_SIZE, TILE_SIZE); |
112 | | - dim3 gridDim( |
113 | | - (N + TILE_SIZE - 1) / TILE_SIZE, |
114 | | - (M + TILE_SIZE - 1) / TILE_SIZE |
115 | | - ); |
116 | | - |
117 | | - bank_conflict_free_sgemm_kernel<TILE_SIZE><<<gridDim, blockDim, 0, stream>>>( |
118 | | - A, B, C, M, K, N |
119 | | - ); |
120 | | - |
121 | | - CUDA_CHECK(cudaGetLastError()); |
| 101 | +template <int TILE_SIZE = 32> |
| 102 | +void launch_bank_conflict_free_sgemm(const float *A, const float *B, float *C, |
| 103 | + int M, int K, int N, |
| 104 | + cudaStream_t stream = 0) { |
| 105 | + dim3 blockDim(TILE_SIZE, TILE_SIZE); |
| 106 | + dim3 gridDim((N + TILE_SIZE - 1) / TILE_SIZE, |
| 107 | + (M + TILE_SIZE - 1) / TILE_SIZE); |
| 108 | + |
| 109 | + bank_conflict_free_sgemm_kernel<TILE_SIZE> |
| 110 | + <<<gridDim, blockDim, 0, stream>>>(A, B, C, M, K, N); |
| 111 | + |
| 112 | + CUDA_CHECK(cudaGetLastError()); |
122 | 113 | } |
123 | 114 |
|
124 | 115 | /** |
125 | 116 | * Alternative: Transposed B storage to avoid bank conflicts |
126 | | - * |
| 117 | + * |
127 | 118 | * Instead of padding, we can store B transposed in shared memory. |
128 | 119 | * This changes the access pattern from column to row access. |
129 | 120 | */ |
130 | | -template<int TILE_SIZE> |
| 121 | +template <int TILE_SIZE> |
131 | 122 | __global__ void bank_conflict_free_transposed_sgemm_kernel( |
132 | | - const float* __restrict__ A, |
133 | | - const float* __restrict__ B, |
134 | | - float* __restrict__ C, |
135 | | - int M, int K, int N |
136 | | -) { |
137 | | - // No padding needed if we transpose B |
138 | | - __shared__ float As[TILE_SIZE][TILE_SIZE + 1]; |
139 | | - __shared__ float BsT[TILE_SIZE][TILE_SIZE + 1]; // B transposed |
140 | | - |
141 | | - int bx = blockIdx.x; |
142 | | - int by = blockIdx.y; |
143 | | - int tx = threadIdx.x; |
144 | | - int ty = threadIdx.y; |
145 | | - |
146 | | - int row = by * TILE_SIZE + ty; |
147 | | - int col = bx * TILE_SIZE + tx; |
148 | | - |
149 | | - float sum = 0.0f; |
150 | | - int numTiles = (K + TILE_SIZE - 1) / TILE_SIZE; |
151 | | - |
152 | | - for (int t = 0; t < numTiles; ++t) { |
153 | | - // Load A normally |
154 | | - int aCol = t * TILE_SIZE + tx; |
155 | | - if (row < M && aCol < K) { |
156 | | - As[ty][tx] = A[row * K + aCol]; |
157 | | - } else { |
158 | | - As[ty][tx] = 0.0f; |
159 | | - } |
160 | | - |
161 | | - // Load B transposed: BsT[tx][ty] instead of Bs[ty][tx] |
162 | | - int bRow = t * TILE_SIZE + ty; |
163 | | - if (bRow < K && col < N) { |
164 | | - BsT[tx][ty] = B[bRow * N + col]; // Note: indices swapped |
165 | | - } else { |
166 | | - BsT[tx][ty] = 0.0f; |
167 | | - } |
168 | | - |
169 | | - __syncthreads(); |
170 | | - |
171 | | - // Now both accesses are row-wise (no bank conflicts) |
172 | | - #pragma unroll |
173 | | - for (int k = 0; k < TILE_SIZE; ++k) { |
174 | | - sum += As[ty][k] * BsT[tx][k]; // Both row accesses |
175 | | - } |
176 | | - |
177 | | - __syncthreads(); |
| 123 | + const float *__restrict__ A, const float *__restrict__ B, |
| 124 | + float *__restrict__ C, int M, int K, int N) { |
| 125 | + // No padding needed if we transpose B |
| 126 | + __shared__ float As[TILE_SIZE][TILE_SIZE + 1]; |
| 127 | + __shared__ float BsT[TILE_SIZE][TILE_SIZE + 1]; // B transposed |
| 128 | + |
| 129 | + int bx = blockIdx.x; |
| 130 | + int by = blockIdx.y; |
| 131 | + int tx = threadIdx.x; |
| 132 | + int ty = threadIdx.y; |
| 133 | + |
| 134 | + int row = by * TILE_SIZE + ty; |
| 135 | + int col = bx * TILE_SIZE + tx; |
| 136 | + |
| 137 | + float sum = 0.0f; |
| 138 | + int numTiles = (K + TILE_SIZE - 1) / TILE_SIZE; |
| 139 | + |
| 140 | + for (int t = 0; t < numTiles; ++t) { |
| 141 | + // Load A normally |
| 142 | + int aCol = t * TILE_SIZE + tx; |
| 143 | + if (row < M && aCol < K) { |
| 144 | + As[ty][tx] = A[row * K + aCol]; |
| 145 | + } else { |
| 146 | + As[ty][tx] = 0.0f; |
178 | 147 | } |
179 | | - |
180 | | - if (row < M && col < N) { |
181 | | - C[row * N + col] = sum; |
| 148 | + |
| 149 | + // Load B transposed: BsT[tx][ty] instead of Bs[ty][tx] |
| 150 | + int bRow = t * TILE_SIZE + ty; |
| 151 | + if (bRow < K && col < N) { |
| 152 | + BsT[tx][ty] = B[bRow * N + col]; // Note: indices swapped |
| 153 | + } else { |
| 154 | + BsT[tx][ty] = 0.0f; |
182 | 155 | } |
| 156 | + |
| 157 | + __syncthreads(); |
| 158 | + |
| 159 | +// Now both accesses are row-wise (no bank conflicts) |
| 160 | +#pragma unroll |
| 161 | + for (int k = 0; k < TILE_SIZE; ++k) { |
| 162 | + sum += As[ty][k] * BsT[tx][k]; // Both row accesses |
| 163 | + } |
| 164 | + |
| 165 | + __syncthreads(); |
| 166 | + } |
| 167 | + |
| 168 | + if (row < M && col < N) { |
| 169 | + C[row * N + col] = sum; |
| 170 | + } |
183 | 171 | } |
184 | 172 |
|
185 | | -template<int TILE_SIZE = 32> |
186 | | -void launch_bank_conflict_free_transposed_sgemm( |
187 | | - const float* A, |
188 | | - const float* B, |
189 | | - float* C, |
190 | | - int M, int K, int N, |
191 | | - cudaStream_t stream = 0 |
192 | | -) { |
193 | | - dim3 blockDim(TILE_SIZE, TILE_SIZE); |
194 | | - dim3 gridDim( |
195 | | - (N + TILE_SIZE - 1) / TILE_SIZE, |
196 | | - (M + TILE_SIZE - 1) / TILE_SIZE |
197 | | - ); |
198 | | - |
199 | | - bank_conflict_free_transposed_sgemm_kernel<TILE_SIZE><<<gridDim, blockDim, 0, stream>>>( |
200 | | - A, B, C, M, K, N |
201 | | - ); |
202 | | - |
203 | | - CUDA_CHECK(cudaGetLastError()); |
| 173 | +template <int TILE_SIZE = 32> |
| 174 | +void launch_bank_conflict_free_transposed_sgemm(const float *A, const float *B, |
| 175 | + float *C, int M, int K, int N, |
| 176 | + cudaStream_t stream = 0) { |
| 177 | + dim3 blockDim(TILE_SIZE, TILE_SIZE); |
| 178 | + dim3 gridDim((N + TILE_SIZE - 1) / TILE_SIZE, |
| 179 | + (M + TILE_SIZE - 1) / TILE_SIZE); |
| 180 | + |
| 181 | + bank_conflict_free_transposed_sgemm_kernel<TILE_SIZE> |
| 182 | + <<<gridDim, blockDim, 0, stream>>>(A, B, C, M, K, N); |
| 183 | + |
| 184 | + CUDA_CHECK(cudaGetLastError()); |
204 | 185 | } |
0 commit comments