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