マルチコア vs SIMD vs CUDA

前回あまりにも 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)

 

Anaconda の Numba を CUDA 9.0 で使う

Volta 対応の CUDA 9.0 が出ましたが、Anaconda の cudatoolkit は 8.0 までしか対応していないため、そのままでは Python  の Numba パッケージが使えません。

というわけで、悪戦苦闘した結果を書いておきます。

環境変数の設定

Numba を利用するためには以下の環境変数を設定する必要があります。

  • NUMBAPRO_NVVM
  • NUMBAPRO_LIBDEVICE

CUDA 9.0 をインストールすると /usr/local/cuda/nvvm/lib64/libnvvm.so.* というライブラリがが作成されますので、このパスを NUMBAPRO_NVVM に設定します。(ファイル名まで含む。)

また、/usr/local/cuda/nvvm/libdevice の下に libdevice.10.bc というファイルがありますので、シンボリックリンクで
libdevice.compute_XX.10.bc というリンクを貼っておきます。(XX は使用している GPU の Compute Capability を指定します。例えば GTX 1080 なら compute_61 となります。)NUMBAPRO_LIBDEVICE には /usr/local/cuda/nvvm/libdevice を指定します。

compute_XX の追加

Anaconda 5.0.0 の Numba には Compute Capability が 5.0 までしか登録されていないため、GPU によっては実行エラーが出てしまいます。ちょっと汚いやり方ですが、Numba パッケージに新しい Compute Capability を追加する方法で回避しました。

GTX 1080 用なら、具体的には ~/anaconda3/lib/python3.6/site-packages/numba/cuda/cudadrv/nvvm.py に以下の2点を追加します。

変数やクラス 修正内容
SUPPORTED_CC リストの最後に (6, 1) を追加
LibDevice クラスの _known_arch 変数 “compute_61” をリストの最後に追加

nVidia のデモプログラムへの修正

NVidia のサイトで説明しているサンプルは numbapro を使ったものなので、numba 用に修正が必要となります。

import numpy as np
from numpy.random import rand
from numba import vectorize

@vectorize(['float32(float32, float32)'], target='cuda')
def VectorAdd(a, b):
    return a + b

def main():
    N = 100000000 # 100M
    print("wait for initialization...")
    A = rand(N).astype(np.float32)
    B = rand(N).astype(np.float32)
    C = np.zeros(N, dtype=np.float32)
    print("OK. Let's start!")

    C = VectorAdd(A, B)
    print("C[:5] = " + str(C[:5]))

if __name__ == "__main__":
    main()

ちなみに、Anaconda3 は Intel の MKL を抱えているっぽいので、この程度の演算だと GPU の方が CPU より遅いという結果となりました。PCIe 上のメモリ転送がボトルネックを引き起こしているためと思われます。

デバイスメモリに結果をキャッシュするような例なら速くなるはず。