CUDA Sample 01

https://ingun37.wordpress.com/2013/01/12/cuda-samplesv5-00_simplematrixmul-%EC%BD%94%EB%93%9C%EB%A6%AC%EB%B7%B0-2/
간단한 CUDA 샘플 소스코드

BY INGUN37 ON JANUARY 12, 2013 • ( LEAVE A COMMENT )

다음은 CUDA 툴킷이 제공하는 CUDA Samples에 포함된 matrixMul 프로젝트를 요약한것

이 코드의 목적은 Fermi이상의 gpu에서는 320*320크기의 행렬과 640*320행렬의 곱을, Fermi 아래에서는 160*160크기의 행렬과 320*160크기의 행렬의 곱을 300번씩 gpu를 이용해서 병렬로 처리하는 것이다. 즉 65536000번의 곱셈과 63488000의 덧셈을, Fermi 아래에서는 8192000번의 곱셈과 7680000번의 덧셈을 각각 300번씩 처리하는 셈이다.

1. block_size를 device가 Fermi 이상이면 32, 아래면 16으로 맞춘다.

1 2 3 4 5 6 7 8 9

int devID=0;cudaDeviceProp devicePropcudaGetDevice(&devID)cudaGetDeviceProperties(&deviceProp, devID);   if (deviceProp.computeMode == cudaComputeModeProhibited) {  fprintf(stderr, "Error: device is running in <Compute Mode Prohibited>, no threads can use ::cudaSetDevice().\n");  exit(EXIT_SUCCESS); }   int block_size = (deviceProp.major < 2) ? 16 : 32;

2. 연산할때 쓸 두 행렬과 결과값을 저장할 행렬의 크기만큼 host memory에 자리를잡고, 데이터값을 초기화한다음 device memory에복사한다.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

//연산할 행렬들을 모아놓을 테이블의 크기 를 정한다. dim3 dimsA(5*2*block_size, 5*2*block_size, 1); dim3 dimsB(5*4*block_size, 5*2*block_size, 1);   //dimsA와 dimsB를 host memory에 집어넣는다. unsigned int size_A = dimsA.x * dimsA.y; unsigned int mem_size_A = sizeof(float) * size_A; float *h_A = (float *)malloc(mem_size_A); unsigned int size_B = dimsB.x * dimsB.y; unsigned int mem_size_B = sizeof(float) * size_B; float *h_B = (float *)malloc(mem_size_B);   // 호스트 메모리 데이터값을 초기화한다. 아래 constantInit함수는 memcpy와 똑같은일을 한다. const float valB = 0.01f; constantInit(h_A, size_A, 1.0f); constantInit(h_B, size_B, valB);   // device memory float *d_A, *d_B, *d_C;   // 결과값을 받을 C행렬을 만든다. dim3 dimsC(dimsB.x, dimsA.y, 1); unsigned int mem_size_C = dimsC.x * dimsC.y * sizeof(float); float *h_C = (float *) malloc(mem_size_C);   //cudaMalloc과 cudaMemcpy는 device memory용 malloc, memcpy 와 같다. cudaMalloc((void **) &d_A, mem_size_A); cudaMalloc((void **) &d_B, mem_size_B); cudaMalloc((void **) &d_C, mem_size_C);   cudaMemcpy(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice); cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice);

3. 시험삼아 kernel을 한번 가동을 한다음, 퍼포먼스를 체크할 타이머를 작동하고, kernel을 300번 실행한다. 그리고 결과값이 올바른지 체크한 후 타이머 결과값과 함께 콘솔창에 출력한다.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92

//kernel함수를 호출할때 쓸 값들. 하나당 block_size^2개의 thread를 가진 (dimsB.x / threads.x, dimsA.y / threads.y)개의 block을 만든다. //(dimsB.x / threads.x, dimsA.y / threads.y) 이 값은 결과값을 저장할 C행렬의 크기와 같다. //즉 thread 하나당 C행렬의 element 하나를 담당하는셈 dim3 threads(block_size, block_size); dim3 grid(dimsB.x / threads.x, dimsA.y / threads.y);   // Performs warmup operation using matrixMul CUDA kernel //왜하는건진 잘 모르겠는데 아마 kernel자체가 워낙 무거워서 //사전체크없이 300번씩이나 돌리기에는 위험이 너무 커서 시험삼아 돌려보는것 같다. if (block_size == 16) {  matrixMulCUDA<16><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x); } else {  matrixMulCUDA<32><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x); }   cudaDeviceSynchronize();   // Allocate CUDA events that we'll use for timing cudaEvent_t start; cudaEventCreate(&start); cudaEventRecord(start, NULL);   int nIter = 300;   for (int j = 0; j < nIter; j++) {  if (block_size == 16)  {  matrixMulCUDA<16><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x);  }  else  {  matrixMulCUDA<32><<< grid, threads >>>(d_C, d_A, d_B, dimsA.x, dimsB.x);  } }   //Record the stop event cudaEventRecord(stop, NULL);   // Wait for the stop event to complete cudaEventSynchronize(stop);   float msecTotal = 0.0f; cudaEventElapsedTime(&msecTotal, start, stop);   // Compute and print the performance float msecPerMatrixMul = msecTotal / nIter; double flopsPerMatrixMul = 2.0 * (double)dimsA.x * (double)dimsA.y * (double)dimsB.x; double gigaFlops = (flopsPerMatrixMul * 1.0e-9f) / (msecPerMatrixMul / 1000.0f); printf(  "Performance= %.2f GFlop/s, Time= %.3f msec, Size= %.0f Ops, WorkgroupSize= %u threads/block\n",  gigaFlops,  msecPerMatrixMul,  flopsPerMatrixMul,  threads.x * threads.y);   // 결과값을 복사한다. cudaMemcpy(h_C, d_C, mem_size_C, cudaMemcpyDeviceToHost);   if (error != cudaSuccess) {  printf("cudaMemcpy (h_C,d_C) returned error code %d, line(%d)\n", error, __LINE__);  exit(EXIT_FAILURE); }   //결과값이 올바른지 체크하는부분. 위에서 A를 1, B를 0.01로 모든값을 초기화했으니 값이 하나라도 0.01*행렬길이가 아니면 틀린것 printf("Checking computed result for correctness: "); bool correct = true;   for (int i = 0; i < (int)(dimsC.x * dimsC.y); i++) {  if (fabs(h_C[i] - (dimsA.x * valB)) > 1e-5)  {  printf("Error! Matrix[%05d]=%.8f, ref=%.8f error term is > 1e-5\n", i, h_C[i], dimsA.x*valB);  correct = false;  } }   printf("%s\n", correct ? "OK" : "FAIL");   //memory 지워주는곳 free(h_A); free(h_B); free(h_C); cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);   cudaDeviceReset();

4. kernel 부분은 기본적인 행렬곱셈과 같은방법이나, thread들이 block별로 나눠져있는 관계로 threadIdx는 상대좌표로, blockIdx를 절대좌표로 이해한다. 현재 실행중인 kernel이 속해있는 blockIdx의 y좌표를 A행렬의 y번째 행의 첫번째 block을 시작점으로, x좌표를 B행렬의 x번째 열의 첫번째block을 시작점으로 잡은 후, kernel의 threadIdx값만큼 상대적으로 옮긴다음 루프를 돌면서 곱셈을한다. 더할때는 shared memory를 쓰게된다.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

//kernel이 속한 block의 index int bx = blockIdx.x; int by = blockIdx.y; //kernel이 속한 thread의 index int tx = threadIdx.x; int ty = threadIdx.y;   //주의해야할점 : 행렬은 2차원배열이지만 이 코드에서는 1차원으로 나열되있으므로 세로로 한줄 내려갈려면 행렬의 가로길이만큼 더해야한다. //A행렬은 행을 곱하므로 루프를 돌때마다 행의 시작점을 가르쳐야 하므로 행렬의 가로길이*block의 세로길이*kernel이 속한 block의 y int aBegin = wA * BLOCK_SIZE * by;//행렬가로줄전   //마지막 -1를 붙이지 않는다면 아래 for에서 a<=aEnd를 a<aEnd로 고쳐야 할듯(100%확실한건아님) //A행렬의 행의 끝점 int aEnd = aBegin + wA - 1; //A행렬은 가로방향이므로 block의 가로길이만큼 건너뛴다 int aStep = BLOCK_SIZE; //B행렬은 열을 곱하므로 루프를 돌때마다 열의 시작점을 가르쳐야 한다. 즉 block 의 가로길이 * kernel이 속한 block의 x int bBegin = BLOCK_SIZE * bx; //B행렬은 세로방향이므로 행렬의 가로길이만큼 건너뛴다. int bStep = BLOCK_SIZE * wB;   float Csub = 0;//결과값을 누적 //block단위로 도는 for문 //바로 곱하지 않고,2차원 배열의 shared memory로 옮긴다음 연산한다. for (int a = aBegin, b = bBegin;  a <= aEnd;  a += aStep, b += bStep) {  __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];   __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];   As[ty][tx] = A[a + wA * ty + tx];  Bs[ty][tx] = B[b + wB * ty + tx];   __syncthreads();   gma unroll   for (int k = 0; k < BLOCK_SIZE; ++k)  {  Csub += As[ty][k] * Bs[k][tx];  }  __syncthreads(); } int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx; C[c + wB * ty + tx] = Csub;