From Florian, 3 Months ago, written in C++.
Embed
  1. #include "cuda_helper.cuh"
  2.  
  3. // C = A * B
  4.  
  5. #define C_ROWS 1000
  6. #define C_COLS 2000
  7. #define A_COLS 3000
  8.  
  9. #define A_ROWS C_ROWS
  10. #define B_ROWS A_COLS
  11. #define B_COLS C_COLS
  12.  
  13. #define BLOCK_DIM 32
  14.  
  15.  
  16. __global__
  17. void MatrixAddKernel(float *A, float *B, float *C) {
  18.         int row = blockIdx.x * blockDim.x + threadIdx.x;
  19.         int col = blockIdx.y * blockDim.y + threadIdx.y;
  20.         if (row >= C_ROWS || col >= C_COLS) {
  21.                 return;
  22.         }
  23.  
  24.         float sum = 0.0;
  25.         for (int k = 0; k < A_COLS; k++) {
  26.                 sum += A[row * A_COLS + k] * B[k * B_COLS + col];
  27.         }
  28.         C[row * C_COLS + col] = sum;
  29. }
  30.  
  31. void cudaMatrixMult(float *A, float *B, float *C, int repetitions, bool warmup) {
  32.         clock_t start = clock();
  33.  
  34.         size_t A_size = A_ROWS * A_COLS * sizeof(float);
  35.         size_t B_size = B_ROWS * B_COLS * sizeof(float);
  36.         size_t C_size = C_ROWS * C_COLS * sizeof(float);
  37.  
  38.         dim3 blockSize(BLOCK_DIM, BLOCK_DIM);
  39.         dim3 gridSize((C_ROWS + BLOCK_DIM - 1) / BLOCK_DIM,
  40.                                   (C_COLS + BLOCK_DIM - 1) / BLOCK_DIM);
  41.  
  42.         float *d_A, *d_B, *d_C;
  43.         cudaMalloc(&d_A, A_size);
  44.         cudaMalloc(&d_B, B_size);
  45.         cudaMalloc(&d_C, C_size);
  46.         cudaMemcpy(d_A, A, A_size, cudaMemcpyHostToDevice);
  47.         cudaMemcpy(d_B, B, B_size, cudaMemcpyHostToDevice);
  48.  
  49.         for (int i = 0; i < repetitions; i++)
  50.         {
  51.                 MatrixAddKernel<<<gridSize, blockSize>>>(d_A, d_B, d_C);
  52.         }
  53.  
  54.         cudaMemcpy(C, d_C, C_size, cudaMemcpyDeviceToHost);
  55.  
  56.         cudaFree(d_A);
  57.         cudaFree(d_B);
  58.         cudaFree(d_C);
  59.  
  60.         if (!warmup)
  61.         {
  62.                 float diff = float(clock() - start) / (CLOCKS_PER_SEC * repetitions);
  63.                 printf("CUDA: %.3lf seconds\n", diff);
  64.         }
  65. }
  66.  
  67. void fillRandomArray(float *A, int numElements) {
  68.         for (int i = 0; i < numElements; i++) {
  69.                 A[i] = rand() / (float)RAND_MAX;
  70.         }
  71. }
  72.  
  73. void verifyResults(float *A, float *B, float *C) {
  74.         for (int row = 0; row < C_ROWS; row++) {
  75.                 for (int col = 0; col < C_COLS; col++) {
  76.                         float sum = 0.0;
  77.                         for (int k = 0; k < A_COLS; k++) {
  78.                                 sum += A[row * A_COLS + k] * B[k * B_COLS + col];
  79.                         }
  80.                         if (fabs(C[row * C_COLS + col] - sum) > 1e-3) {
  81.                                 fprintf(stderr, "Result verification failed at element %d: %f vs. %f!\n", row, C[row * C_COLS + col], sum);
  82.                                 exit(EXIT_FAILURE);
  83.                         }
  84.                 }
  85.         }
  86. }
  87.  
  88. void sequentialMatrixMult(float *A, float *B, float *C) {
  89.         clock_t start = clock();
  90.  
  91.         for (int row = 0; row < C_ROWS; row++) {
  92.                 for (int col = 0; col < C_COLS; col++) {
  93.                         float sum = 0.0;
  94.                         for (int k = 0; k < A_COLS; k++) {
  95.                                 sum += A[row * A_COLS + k] * B[k * B_COLS + col];
  96.                         }
  97.                         C[row * C_COLS + col] = sum;
  98.                 }
  99.         }
  100.  
  101.         float diff = float(clock() - start) / CLOCKS_PER_SEC;
  102.         printf("Sequential: %.3lf seconds\n", diff);
  103. }
  104.  
  105. int main() {
  106.         int nofElemA = A_ROWS * A_COLS;
  107.         float *h_A = (float *)malloc(nofElemA * sizeof(float));
  108.         handleAllocationError(h_A);
  109.         fillRandomArray(h_A, nofElemA);
  110.  
  111.         int nofElemB = B_ROWS * B_COLS;
  112.         float *h_B = (float *)malloc(nofElemB * sizeof(float));
  113.         handleAllocationError(h_B);
  114.         fillRandomArray(h_B, nofElemB);
  115.        
  116.         int nofElemC = C_ROWS * C_COLS;
  117.         float *h_C = (float *)malloc(nofElemC * sizeof(float));
  118.         handleAllocationError(h_C);
  119.  
  120.         cudaMatrixMult(h_A, h_B, h_C, 2, true);
  121.         verifyResults(h_A, h_B, h_C);
  122.         cudaMatrixMult(h_A, h_B, h_C, 4, false);
  123.  
  124.         sequentialMatrixMult(h_A, h_B, h_C);
  125.  
  126.         free(h_A);
  127.         free(h_B);
  128.         free(h_C);
  129.  
  130.         return 0;
  131. }
  132.