投稿者: 管理者TeqStock

前回あまりにも Numba が遅いのでおかしいと思い、マルチコア、SIMD、CUDA を比較するためのテストプログラムを作ってみた。

結果は、以下の通り。シングル CPU の非SIMD演算が最も遅く、CUDA が最も速いという予想通りの結果になりました。

version: 201307
#20 processors, 20 threads
start parallel calculation.
<<single-nosimd>>
elapsed 1.115483 seconds

<<single-simd>>
elapsed 0.535289 seconds

<<multicore-nosimd>>
elapsed 0.189586 seconds

<<multicore-simd>>
elapsed 0.128129 seconds

<<cuda>>
elapsed 0.021330 seconds

1.234571 1.581539 1.109199 1.103452 0.831745 1.106268 0.878185 1.868425 1.353009 0.748571 
1.234571 1.581539 1.109199 1.103452 0.831745 1.106268 0.878185 1.868425 1.353009 0.748571 

ただし、cudaMallocManaged() による Unified Memory を使うと速度が極端に低下する現象が見られました。（20倍くらい遅くなる）

やはり PCIe バスの遅さが実行時間に影響していると考えてよいようです。Numba のように自動的にメモリ配置してしまうライブラリは却って速度低下の原因となることがわかりました。

というわけで、CUDA を使うときはなるべく Device Memory に閉じるように設計しましょう、という当たり前の結論となりました。

この下は今回使ったソースコードです。マルチコア化には OpenMP を、SIMDには AVX2 を使用しました。

[Nelem.h]

#define NELEM (400 * 1000 * 1000) # 400M elements

[sample_cuda.h]

#if !defined(SAMPLE_CUDA_H_)
#define SAMPLE_CUDA_H_

#if defined(__cplusplus)
extern "C"
#endif
void prepare(float *a, float *b, float *c);

#if defined(__cplusplus)
extern "C"
#endif
float *getResult(float *);

#if defined(__cplusplus)
extern "C"
#endif
void release();

#if defined(__cplusplus)
extern "C"
#endif
void wrapper_for_cuda(void);

#endif

[sample.c]

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <omp.h>
#include <time.h>
#include <x86intrin.h>
#include "Nelem.h"
#include "sample_cuda.h"

#define alignof(x) __alignof__(x)

float *a, *b, *c;

double elapsed(struct timespec start, struct timespec end) {
  time_t second = end.tv_sec - start.tv_sec;
  long nsecond = end.tv_nsec - start.tv_nsec;
  return (double)second + (double)nsecond / (1000.0 * 1000.0 * 1000.0);
}

typedef void (*CBFUNC)(void);

static inline void simd_add(int i) {
    __m256 ma = _mm256_load_ps(&a[i]);
    __m256 mb = _mm256_load_ps(&b[i]);
    __m256 mc = _mm256_add_ps(ma, mb);
    _mm256_store_ps(&c[i], mc);
}

static void single_nosimd_add(void) {
  int i;
  for (i = 0; i < NELEM; ++i) {
    c[i] = a[i] + b[i];
  }
}

static void single_simd_add(void) {
  int i;
  for (i = 0; i < NELEM; i += 8) {
    simd_add(i);
  }
}

static void multicore_nosimd_add(void) {
  int i;

  #pragma omp parallel for private(i)
  for (i = 0; i < NELEM; ++i) {
    c[i] = a[i] + b[i];
  }
}

static void multicore_simd_add(void) {
  int i;
  #pragma omp parallel for private(i)
  for (i = 0; i < NELEM; i += 8) {
    simd_add(i);
  }
}

char *cbtitle[] = { "single-nosimd", "single-simd",
                    "multicore-nosimd", "multicore-simd",
                    "cuda" };
CBFUNC cbfunc[] = {
  single_nosimd_add, single_simd_add,
  multicore_nosimd_add, multicore_simd_add,
  wrapper_for_cuda
};
#define NCBFUNC (sizeof(cbfunc)/sizeof(cbfunc[0]))

int main() {
  int i;
  float *f;
  struct timespec start, end;

#if defined(_OPENMP)
  {
    int np = omp_get_num_procs();
    int nth = omp_get_max_threads();
    printf("version: %d\n", _OPENMP);
    printf("#%d processors, %d threads\n", np, nth);
  }
#endif
 
  a = (float *)_mm_malloc(sizeof(float) * NELEM, 32);
  b = (float *)_mm_malloc(sizeof(float) * NELEM, 32);
  c = (float *)_mm_malloc(sizeof(float) * NELEM, 32);

  for (i = 0; i < NELEM; ++i) {
    a[i] = ((float)random() / (float)RAND_MAX);
    b[i] = ((float)random() / (float)RAND_MAX);
  }
  prepare(a, b, c);

  printf("start parallel calculation.\n");

  for (i = 0; i < NCBFUNC; ++i) {
    clock_gettime(CLOCK_REALTIME, &start);
    printf("<<%s>>\n", cbtitle[i]);
    (*cbfunc[i])();
    clock_gettime(CLOCK_REALTIME, &end);
    printf("elapsed %f seconds\n\n", elapsed(start, end));
  }
  

  for (i = 0; i < 10 && i < NELEM; ++i) {
    printf("%f ", c[i]);
  }
  fputs("\n", stdout);

  f = getResult(c);
  for (i = 0; i < 10 && i < NELEM; ++i) {
    printf("%f ", f[i]);
  }
  fputs("\n", stdout);
 
  _mm_free(a);
  _mm_free(b);
  _mm_free(c);
  release();
  return 0;
}

[sample_cuda.cu]

#include 
#include 
#include "Nelem.h"

#define NTHREAD (256)

float *pa, *pb, *pc;

__global__
void kernel(float *xc, float *a, float *b) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < NELEM) {
    xc[i] = a[i] + b[i];
  }
  __syncthreads();
}

__host__
void invoke_kernel(float *c, float *a, float *b) {
  int block = NELEM / NTHREAD;
  kernel<<<block, NTHREAD>>>(c, a, b);
}

extern "C" void prepare(float *a, float *b, float *c) {
  cudaMalloc((void**)&pa, sizeof(float) * NELEM);
  cudaMalloc((void**)&pb, sizeof(float) * NELEM);
  cudaMalloc((void**)&pc, sizeof(float) * NELEM);
  cudaMemcpy(pa, a, sizeof(float) * NELEM, cudaMemcpyHostToDevice);
  cudaMemcpy(pb, b, sizeof(float) * NELEM, cudaMemcpyHostToDevice);
}

extern "C" float *getResult(float *c) {
  cudaMemcpy(c, pc, sizeof(float) * NELEM, cudaMemcpyDeviceToHost);
  return c;
}

extern "C" void release(void) {
  cudaFree(pa);
  cudaFree(pb);
  cudaFree(pc);
}

extern "C" void wrapper_for_cuda(void) {
  invoke_kernel(pc, pa, pb);
  cudaDeviceSynchronize();
}

[Makefile]

SRC=sample.c
OBJ=$(SRC:.c=.o)
CUDASRC=sample_cuda.cu
CUDAOBJ=$(CUDASRC:.cu=.o)
TARGET_OPENMP=sample

all: $(TARGET_OPENMP)

$(TARGET_OPENMP): $(OBJ) $(CUDAOBJ)
        nvcc -g -O0 -o $(TARGET_OPENMP) --gpu-architecture=compute_61 --gpu-code=sm_61 --compiler-options="-fopenmp -DAVX2 -march=native" $(OBJ) $(CUDAOBJ)

$(CUDAOBJ): $(CUDASRC)
        nvcc --gpu-architecture=compute_61 --gpu-code=sm_61 -c $(CUDASRC)

.c.o:
        gcc -c -g -O0 -fopenmp -DAVX2 -march=native $<

clean:
        -rm $(TARGET_OPENMP) $(OBJ) $(CUDAOBJ)