|
| 1 | +/** |
| 2 | + * Copyright 1993-2012 NVIDIA Corporation. All rights reserved. |
| 3 | + * |
| 4 | + * Please refer to the NVIDIA end user license agreement (EULA) associated |
| 5 | + * with this source code for terms and conditions that govern your use of |
| 6 | + * this software. Any use, reproduction, disclosure, or distribution of |
| 7 | + * this software and related documentation outside the terms of the EULA |
| 8 | + * is strictly prohibited. |
| 9 | + * |
| 10 | + */ |
| 11 | + |
| 12 | +/** |
| 13 | + * Matrix multiplication: C = A * B. |
| 14 | + * Host code. |
| 15 | + * |
| 16 | + * This sample implements matrix multiplication as described in Chapter 3 |
| 17 | + * of the programming guide. |
| 18 | + * It has been written for clarity of exposition to illustrate various CUDA |
| 19 | + * programming principles, not with the goal of providing the most |
| 20 | + * performant generic kernel for matrix multiplication. |
| 21 | + * |
| 22 | + * See also: |
| 23 | + * V. Volkov and J. Demmel, "Benchmarking GPUs to tune dense linear algebra," |
| 24 | + * in Proc. 2008 ACM/IEEE Conf. on Superconducting (SC '08), |
| 25 | + * Piscataway, NJ: IEEE Press, 2008, pp. Art. 31:1-11. |
| 26 | + */ |
| 27 | + |
| 28 | +// System includes |
| 29 | +#include <stdio.h> |
| 30 | +#include <assert.h> |
| 31 | + |
| 32 | +// CUDA runtime |
| 33 | +#include <cuda_runtime.h> |
| 34 | + |
| 35 | +// Helper functions and utilities to work with CUDA |
| 36 | +//#include <helper_functions.h> |
| 37 | + |
| 38 | +/** |
| 39 | + * Matrix multiplication (CUDA Kernel) on the device: C = A * B |
| 40 | + * wA is A's width and wB is B's width |
| 41 | + */ |
| 42 | +template <int BLOCK_SIZE> __global__ void |
| 43 | +matrixMulCUDA(float *C, float *A, float *B, int wA, int wB) |
| 44 | +{ |
| 45 | + // Block index |
| 46 | + int bx = blockIdx.x; |
| 47 | + int by = blockIdx.y; |
| 48 | + |
| 49 | + // Thread index |
| 50 | + int tx = threadIdx.x; |
| 51 | + int ty = threadIdx.y; |
| 52 | + |
| 53 | + // Index of the first sub-matrix of A processed by the block |
| 54 | + int aBegin = wA * BLOCK_SIZE * by; |
| 55 | + |
| 56 | + // Index of the last sub-matrix of A processed by the block |
| 57 | + int aEnd = aBegin + wA - 1; |
| 58 | + |
| 59 | + // Step size used to iterate through the sub-matrices of A |
| 60 | + int aStep = BLOCK_SIZE; |
| 61 | + |
| 62 | + // Index of the first sub-matrix of B processed by the block |
| 63 | + int bBegin = BLOCK_SIZE * bx; |
| 64 | + |
| 65 | + // Step size used to iterate through the sub-matrices of B |
| 66 | + int bStep = BLOCK_SIZE * wB; |
| 67 | + |
| 68 | + // Csub is used to store the element of the block sub-matrix |
| 69 | + // that is computed by the thread |
| 70 | + float Csub = 0; |
| 71 | + |
| 72 | + // Loop over all the sub-matrices of A and B |
| 73 | + // required to compute the block sub-matrix |
| 74 | + for (int a = aBegin, b = bBegin; |
| 75 | + a <= aEnd; |
| 76 | + a += aStep, b += bStep) |
| 77 | + { |
| 78 | + |
| 79 | + // Declaration of the shared memory array As used to |
| 80 | + // store the sub-matrix of A |
| 81 | + __shared__ float As[BLOCK_SIZE][BLOCK_SIZE]; |
| 82 | + |
| 83 | + // Declaration of the shared memory array Bs used to |
| 84 | + // store the sub-matrix of B |
| 85 | + __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE]; |
| 86 | + |
| 87 | + // Load the matrices from device memory |
| 88 | + // to shared memory; each thread loads |
| 89 | + // one element of each matrix |
| 90 | + As[ty][tx] = A[a + wA * ty + tx]; |
| 91 | + Bs[ty][tx] = B[b + wB * ty + tx]; |
| 92 | + |
| 93 | + // Synchronize to make sure the matrices are loaded |
| 94 | + __syncthreads(); |
| 95 | + |
| 96 | + // Multiply the two matrices together; |
| 97 | + // each thread computes one element |
| 98 | + // of the block sub-matrix |
| 99 | +#pragma unroll |
| 100 | + |
| 101 | + for (int k = 0; k < BLOCK_SIZE; ++k) |
| 102 | + { |
| 103 | + Csub += As[ty][k] * Bs[k][tx]; |
| 104 | + } |
| 105 | + |
| 106 | + // Synchronize to make sure that the preceding |
| 107 | + // computation is done before loading two new |
| 108 | + // sub-matrices of A and B in the next iteration |
| 109 | + __syncthreads(); |
| 110 | + } |
| 111 | + |
| 112 | + // Write the block sub-matrix to device memory; |
| 113 | + // each thread writes one element |
| 114 | + int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx; |
| 115 | + C[c + wB * ty + tx] = Csub; |
| 116 | +} |
| 117 | + |
| 118 | +void constantInit(float *data, int size, float val) |
| 119 | +{ |
| 120 | + for (int i = 0; i < size; ++i) |
| 121 | + { |
| 122 | + data[i] = val; |
| 123 | + } |
| 124 | +} |
| 125 | + |
| 126 | +/** |
| 127 | + * Run a simple test of matrix multiplication using CUDA |
| 128 | + */ |
| 129 | +int matrixMultiply(dim3 &dimsA, dim3 &dimsB) |
| 130 | +{ |
| 131 | + int block_size = 32; |
| 132 | + // Allocate host memory for matrices A and B |
| 133 | + unsigned int size_A = dimsA.x * dimsA.y; |
| 134 | + unsigned int mem_size_A = sizeof(float) * size_A; |
| 135 | + float *h_A = (float *)malloc(mem_size_A); |
| 136 | + unsigned int size_B = dimsB.x * dimsB.y; |
| 137 | + unsigned int mem_size_B = sizeof(float) * size_B; |
| 138 | + float *h_B = (float *)malloc(mem_size_B); |
| 139 | + |
| 140 | + // Initialize host memory |
| 141 | + const float valB = 0.01f; |
| 142 | + constantInit(h_A, size_A, 1.0f); |
| 143 | + constantInit(h_B, size_B, valB); |
| 144 | + |
| 145 | + // Allocate device memory |
| 146 | + float *d_A, *d_B, *d_C; |
| 147 | + |
| 148 | + // Allocate host matrix C |
| 149 | + dim3 dimsC(dimsB.x, dimsA.y, 1); |
| 150 | + unsigned int mem_size_C = dimsC.x * dimsC.y * sizeof(float); |
| 151 | + float *h_C = (float *) malloc(mem_size_C); |
| 152 | + |
| 153 | + if (h_C == NULL) |
| 154 | + { |
| 155 | + fprintf(stderr, "Failed to allocate host matrix C!\n"); |
| 156 | + exit(EXIT_FAILURE); |
| 157 | + } |
| 158 | + |
| 159 | + cudaError_t error; |
| 160 | + |
| 161 | + error = cudaMalloc((void **) &d_A, mem_size_A); |
| 162 | + |
| 163 | + if (error != cudaSuccess) |
| 164 | + { |
| 165 | + printf("cudaMalloc d_A returned error code %d, line(%d)\n", error, __LINE__); |
| 166 | + exit(EXIT_FAILURE); |
| 167 | + } |
| 168 | + |
| 169 | + error = cudaMalloc((void **) &d_B, mem_size_B); |
| 170 | + |
| 171 | + if (error != cudaSuccess) |
| 172 | + { |
| 173 | + printf("cudaMalloc d_B returned error code %d, line(%d)\n", error, __LINE__); |
| 174 | + exit(EXIT_FAILURE); |
| 175 | + } |
| 176 | + |
| 177 | + error = cudaMalloc((void **) &d_C, mem_size_C); |
| 178 | + |
| 179 | + if (error != cudaSuccess) |
| 180 | + { |
| 181 | + printf("cudaMalloc d_C returned error code %d, line(%d)\n", error, __LINE__); |
| 182 | + exit(EXIT_FAILURE); |
| 183 | + } |
| 184 | + |
| 185 | + // copy host memory to device |
| 186 | + error = cudaMemcpy(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice); |
| 187 | + |
| 188 | + if (error != cudaSuccess) |
| 189 | + { |
| 190 | + printf("cudaMemcpy (d_A,h_A) returned error code %d, line(%d)\n", error, __LINE__); |
| 191 | + exit(EXIT_FAILURE); |
| 192 | + } |
| 193 | + |
| 194 | + error = cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice); |
| 195 | + |
| 196 | + if (error != cudaSuccess) |
| 197 | + { |
| 198 | + printf("cudaMemcpy (d_B,h_B) returned error code %d, line(%d)\n", error, __LINE__); |
| 199 | + exit(EXIT_FAILURE); |
| 200 | + } |
| 201 | + |
| 202 | + // Setup execution parameters |
| 203 | + dim3 threads(block_size, block_size); |
| 204 | + dim3 grid(dimsB.x / threads.x, dimsA.y / threads.y); |
| 205 | + |
| 206 | + // Create and start timer |
| 207 | + printf("\tRunning Matrix Multiplication CUDA Kernel...\n"); |
| 208 | + |
| 209 | + // Performs warmup operation using matrixMul CUDA kernel |
| 210 | + if (block_size == 16) |
| 211 | + { |
| 212 | + matrixMulCUDA<16><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x); |
| 213 | + } |
| 214 | + else |
| 215 | + { |
| 216 | + matrixMulCUDA<32><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x); |
| 217 | + } |
| 218 | + |
| 219 | + cudaDeviceSynchronize(); |
| 220 | + |
| 221 | + // Allocate CUDA events that we'll use for timing |
| 222 | + cudaEvent_t start; |
| 223 | + error = cudaEventCreate(&start); |
| 224 | + |
| 225 | + if (error != cudaSuccess) |
| 226 | + { |
| 227 | + fprintf(stderr, "Failed to create start event (error code %s)!\n", cudaGetErrorString(error)); |
| 228 | + exit(EXIT_FAILURE); |
| 229 | + } |
| 230 | + |
| 231 | + cudaEvent_t stop; |
| 232 | + error = cudaEventCreate(&stop); |
| 233 | + |
| 234 | + if (error != cudaSuccess) |
| 235 | + { |
| 236 | + fprintf(stderr, "Failed to create stop event (error code %s)!\n", cudaGetErrorString(error)); |
| 237 | + exit(EXIT_FAILURE); |
| 238 | + } |
| 239 | + |
| 240 | + // Record the start event |
| 241 | + error = cudaEventRecord(start, NULL); |
| 242 | + |
| 243 | + if (error != cudaSuccess) |
| 244 | + { |
| 245 | + fprintf(stderr, "Failed to record start event (error code %s)!\n", cudaGetErrorString(error)); |
| 246 | + exit(EXIT_FAILURE); |
| 247 | + } |
| 248 | + |
| 249 | + // Execute the kernel |
| 250 | + int nIter = 300; |
| 251 | + |
| 252 | + for (int j = 0; j < nIter; j++) |
| 253 | + { |
| 254 | + if (block_size == 16) |
| 255 | + { |
| 256 | + matrixMulCUDA<16><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x); |
| 257 | + } |
| 258 | + else |
| 259 | + { |
| 260 | + matrixMulCUDA<32><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x); |
| 261 | + } |
| 262 | + } |
| 263 | + |
| 264 | + // Record the stop event |
| 265 | + error = cudaEventRecord(stop, NULL); |
| 266 | + |
| 267 | + if (error != cudaSuccess) |
| 268 | + { |
| 269 | + fprintf(stderr, "Failed to record stop event (error code %s)!\n", cudaGetErrorString(error)); |
| 270 | + exit(EXIT_FAILURE); |
| 271 | + } |
| 272 | + |
| 273 | + // Wait for the stop event to complete |
| 274 | + error = cudaEventSynchronize(stop); |
| 275 | + |
| 276 | + if (error != cudaSuccess) |
| 277 | + { |
| 278 | + fprintf(stderr, "Failed to synchronize on the stop event (error code %s)!\n", cudaGetErrorString(error)); |
| 279 | + exit(EXIT_FAILURE); |
| 280 | + } |
| 281 | + |
| 282 | + float msecTotal = 0.0f; |
| 283 | + error = cudaEventElapsedTime(&msecTotal, start, stop); |
| 284 | + |
| 285 | + if (error != cudaSuccess) |
| 286 | + { |
| 287 | + fprintf(stderr, "Failed to get time elapsed between events (error code %s)!\n", cudaGetErrorString(error)); |
| 288 | + exit(EXIT_FAILURE); |
| 289 | + } |
| 290 | + |
| 291 | + // Compute and print the performance |
| 292 | + float msecPerMatrixMul = msecTotal / nIter; |
| 293 | + double flopsPerMatrixMul = 2.0 * (double)dimsA.x * (double)dimsA.y * (double)dimsB.x; |
| 294 | + double gigaFlops = (flopsPerMatrixMul * 1.0e-9f) / (msecPerMatrixMul / 1000.0f); |
| 295 | + printf( |
| 296 | + "\tPerformance= %.2f GFlop/s, Time= %.3f msec, Size= %.0f Ops\n", |
| 297 | + gigaFlops, |
| 298 | + msecPerMatrixMul, |
| 299 | + flopsPerMatrixMul); |
| 300 | + |
| 301 | + // Copy result from device to host |
| 302 | + error = cudaMemcpy(h_C, d_C, mem_size_C, cudaMemcpyDeviceToHost); |
| 303 | + |
| 304 | + if (error != cudaSuccess) |
| 305 | + { |
| 306 | + printf("cudaMemcpy (h_C,d_C) returned error code %d, line(%d)\n", error, __LINE__); |
| 307 | + exit(EXIT_FAILURE); |
| 308 | + } |
| 309 | + |
| 310 | + printf("\tChecking computed result for correctness: "); |
| 311 | + bool correct = true; |
| 312 | + |
| 313 | + for (int i = 0; i < (int)(dimsC.x * dimsC.y); i++) |
| 314 | + { |
| 315 | + if (fabs(h_C[i] - (dimsA.x * valB)) > 1e-5) |
| 316 | + { |
| 317 | + printf("Error! Matrix[%05d]=%.8f, ref=%.8f error term is > 1e-5\n", i, h_C[i], dimsA.x*valB); |
| 318 | + correct = false; |
| 319 | + } |
| 320 | + } |
| 321 | + |
| 322 | + printf("%s\n", correct ? "OK" : "FAIL"); |
| 323 | + |
| 324 | + // Clean up memory |
| 325 | + free(h_A); |
| 326 | + free(h_B); |
| 327 | + free(h_C); |
| 328 | + cudaFree(d_A); |
| 329 | + cudaFree(d_B); |
| 330 | + cudaFree(d_C); |
| 331 | + |
| 332 | + printf("\n\tNote: For peak performance, please refer to the matrixMulCUBLAS example.\n"); |
| 333 | + |
| 334 | + if (correct) |
| 335 | + { |
| 336 | + return EXIT_SUCCESS; |
| 337 | + } |
| 338 | + else |
| 339 | + { |
| 340 | + return EXIT_FAILURE; |
| 341 | + } |
| 342 | +} |
0 commit comments