CUDA lib GPUMLib 소개

요즘 CUDA에 관심이 있어서 CUDA를 설치도 해보고, 실행도 시켜 보았습니다.
현재 만들어진 좋은 library없을까 검색하던 중 CUDA 소스도 참고할만한 library가 있었습니다.

GPUMLib Download: http://sourceforge.net/projects/gpumlib
GPUMLib Homepage: http://gpumlib.sourceforge.net/

GPU machine learning 하는 open source입니다. C++ 과 CUDA로 되어있으며, 이미 BP,MBP 가 구현이 되어있어서 소스를 참고할 수 있어서 좋았습니다.

GPUMLib is an open source Graphic Processing Unit Machine Learning Library. This library aims to provide machine learning researchers and practitioners with a high performance library by taking advantage of the GPU enormous computational power. The library is developed in C++ and CUDA.

CUDA kernel 소스를 참고하자면 대표적으로 배열에 들어있는 합을 구할때 CUDA에서는 어떻게 하는게 효율적일까 고민을 좀 했습니다. BP 혹은 MLP 에서 이전 노드와 가중치를 곱한 모든 node의 합이 다음 노드값을 갖기 때문에 반듯이 필요한 부분입니다. Host 쪽에서 구현해도 되지만 node가 많아진다면 속도가 떨어지겠죠.
해당 소스에서는 다음과 같이 구현했습니다.

template <int blockSize> __device__ __forceinline__ void SumBeforeWarp(cudafloat * s) {
 if (blockSize >= 1024) {
  if (threadIdx.x < 512) s[threadIdx.x] += s[threadIdx.x + 512];
  __syncthreads();
 }

 if (blockSize >= 512) {
  if (threadIdx.x < 256) s[threadIdx.x] += s[threadIdx.x + 256];
  __syncthreads();
 }

 if (blockSize >= 256) {
  if (threadIdx.x < 128) s[threadIdx.x] += s[threadIdx.x + 128];
  __syncthreads();
 }

 if (blockSize >= 128) {
  if (threadIdx.x < 64) s[threadIdx.x] += s[threadIdx.x + 64];
  __syncthreads();
 }
}

template <int blockSize> __device__ __forceinline__ void SumWarp(volatile cudafloat * s) {
 if (blockSize >= 64) s[threadIdx.x] += s[threadIdx.x + 32];
 if (blockSize >= 32) s[threadIdx.x] += s[threadIdx.x + 16];
 if (blockSize >= 16) s[threadIdx.x] += s[threadIdx.x + 8];
 if (blockSize >= 8) s[threadIdx.x] += s[threadIdx.x + 4];
 if (blockSize >= 4) s[threadIdx.x] += s[threadIdx.x + 2];
 if (blockSize >= 2) s[threadIdx.x] += s[threadIdx.x + 1];
}

앞쪽 코드 SumBeforeWarp 함수에서는 절반씩 나누어서 합을 더하고 __syncthreads()를 호출하여 모든 thread들이 barrier에 멈추도록 합니다. 따라서 반반씩 하면 위와같은 코드가 나오는데 Index가 32 보다 작아지면 SumWarp라는 함수를 이용합니다.
아래는 호출하는 쪽 코드인데 threadIdx.x <32라는 조건이 있습니다. 32를 비교하는 이유는 NVIDIA GPU에서는 Warp단위로 스케줄링을 한다고 하는데요. Warp 단위가 32라고 합니다. 그래서 __syncthreads를 호출 안한다고 합니다. 문제는 Warp단위가 현재는 32이지만 나중에 바뀔 수도 있고요 __syncthreads 를 많이 호출하면 효율이 떨어진다고 합니다.

template <int blockSize> KERNEL FireLayerNeurons(cudafloat * inputs, cudafloat * weights, cudafloat * m, int mOffset, int totalNeuronsWithSelectiveActivation, cudafloat * outputs, int numInputs) {
 extern __shared__ cudafloat iw[];

 iw[threadIdx.x] = CUDA_VALUE(0.0);
 for(int i = threadIdx.x; i <= numInputs; i += blockDim.x) {
  cudafloat i_w = weights[NEURON * (numInputs + 1) + i];
  if (i > BIAS) i_w *= inputs[PATTERN * numInputs + (i - 1)];  
  iw[threadIdx.x] += i_w;
 }
 __syncthreads();

 SumBeforeWarp<blockSize>(iw);

 if (threadIdx.x < 32) {
  SumWarp<blockSize>(iw);
 
  if (threadIdx.x == 0) {
   cudafloat output = CUDA_SIGMOID(iw[0]);
   if (m != nullptr) output *= m[PATTERN * totalNeuronsWithSelectiveActivation + NEURON + mOffset];
   outputs[PATTERN * NUM_NEURONS + NEURON] = output;
  }
 }
}