疎行列-ベクトル積性能評価: CUDA、MIC、そしてCPU

出典: トータル・ディスクロージャ・サイト(事実をありのままに)

いくつかの高速化ソリューションの案件で、疎行列-ベクトル積(SpMV)を高速化したい、という話が上がってきている。CPUで実行するにも、MKLをはじめとする数値計算ライブラリの高速化がよく効いていてまずまず速く、CUDA付属のCUSPARSEライブラリを使えば、GPGPUを使って更なる性能向上を享受できる。

ところが、Xeon Phi上でMIC向けMKLを使ってSpMVを計算させようとすると、どうにも期待を裏切るように思われる。MKL_MIC_ENABLE環境変数を設定しても、MKLは実行をオフロードせずにCPU上で計算するようであり、またNative実行でMKLを走らせても、期待したような性能を発揮しない。今後のチューニング作業の出発点として、また読者の皆様の知恵を拝借したく、Wikiの記事としてメモを残すこととする。

ベンチマーク用ソースコード

#include <cstdio>
#include <cstdlib>

#include <algorithm>
#include <numeric>
#include <vector>

#ifdef CUDA
#include <cuda_runtime.h>
#include <cusparse_v2.h>
#else
#include <mkl.h>
#endif
#include <sys/time.h>

const size_t n = 50000;
const size_t m = 50000;

static size_t
zrandom()
{
  return (size_t) random() * ((size_t) RAND_MAX + (size_t) 1) + (size_t) random();
}

static double
drandom()
{
  return (double) zrandom() / (double) (((size_t) RAND_MAX + (size_t) 1) * ((size_t) RAND_MAX + (size_t) 1));
}

int
main()
{
  double *b(new double[n]), *x(new double[m]);
  std::vector<bool> *iarray;
  srandom(0);
  for (size_t i = 0; i < n; ++i) {
    b[i] = drandom();
  }
#ifdef CUDA
  cusparseHandle_t handle;
  cusparseCreate(&handle);
  double *d_b, *d_x;
  cudaMalloc((void **)&d_b, n * sizeof(double));
  cudaMalloc((void **)&d_x, m * sizeof(double));
#endif
  for (size_t percent = 1; percent <= 10; ++percent) {
    iarray = new std::vector<bool>(m * n, false);
    size_t percentn = m * n * percent / 100;
    std::fill(iarray->begin(), iarray->begin() + percentn, true);
    for (size_t i = 0; i < percentn; ++i) {
      size_t j = zrandom() % (m * n);
      bool tmp = (*iarray)[i];
      (*iarray)[i] = (*iarray)[j];
      (*iarray)[j] = tmp;
    }
    int *ia = new int[m + 1];
    int *ja = new int[percentn];
    {
      size_t percenti = 0;
      for (size_t i = 0; i < m; ++i) {
	ia[i] = percenti;
	for (size_t j = 0; j < n; ++j) {
	  if ((*iarray)[i * n + j]) {
	    ja[percenti++] = j;
	  }
	}
      }
      ia[m] = percenti;
    }
    delete iarray;
    double *a = new double[percentn];
    for (size_t i = 0; i < percentn; ++i) {
      a[i] = drandom();
    }
    timeval start, end;
#ifdef CUDA
    cusparseMatDescr_t descrA;
    cusparseCreateMatDescr(&descrA);
    double *d_a;
    cudaMalloc((void **)&d_a, percentn * sizeof(double));
    cudaMemcpy(d_a, a, percentn * sizeof(double),
	       cudaMemcpyHostToDevice);
    int *d_ia, *d_ja;
    cudaMalloc((void **)&d_ia, (m + 1) * sizeof(int));
    cudaMemcpy(d_ia, ia, (m + 1) * sizeof(int),
	       cudaMemcpyHostToDevice);
    cudaMalloc((void **)&d_ja, percentn * sizeof(int));
    cudaMemcpy(d_ja, ja, percentn * sizeof(int),
	       cudaMemcpyHostToDevice);
#else
#pragma offload_transfer target(mic:0) in(a[0:percentn], ia[0:m + 1], ja[0:percentn] : alloc_if(1) free_if(0))
#endif
    gettimeofday(&start, 0);
#ifdef CUDA
    cudaMemcpy(d_b, b, n * sizeof(double),
	       cudaMemcpyHostToDevice);
    const double alpha = 1.0;
    const double beta = 0.0;
    cusparseDcsrmv(handle,
		   CUSPARSE_OPERATION_NON_TRANSPOSE,
		   m, n, percentn,
		   &alpha, descrA,
		   d_a, d_ia, d_ja,
		   d_b, &beta, d_x);
    cudaMemcpy(x, d_x, m * sizeof(double),
	       cudaMemcpyDeviceToHost);
#else
#pragma offload target(mic:0) in(b[0:n]) out(x[0:m]) nocopy(a[0:percentn], ia[0:m + 1], ja[0:percentn])
    {
      char transa = 'N';
      int m_ = m;
      int k_ = n;
      double alpha = 1.0;
      char matdescra[] = "G**C**";
      double beta = 0.0;
      mkl_dcsrmv(&transa, &m_, &k_, &alpha, matdescra, a, ja, ia, ia + 1, b, &beta, x);
    }
#endif
    gettimeofday(&end, 0);
    long long us;
    us = (end.tv_sec - start.tv_sec) * 1000000 + end.tv_usec - start.tv_usec;
    printf("%2zu %% (%10zu elements):\n"
           "%lld usec\n"
           "%e MFLOPS\n"
           "CHECKSUM: %e\n",
           percent, percentn, us, 2.0 * percentn / us,
           std::accumulate(x, x + m, 0.0));
#ifdef CUDA
    cudaFree(d_ja);
    cudaFree(d_ia);
    cudaFree(d_a);
#else
#pragma offload_transfer target(mic:0) nocopy(a[0:percentn], ia[0:m + 1], ja[0:percentn] : alloc_if(0) free_if(1))
#endif
    delete[] a;
    delete[] ja;
    delete[] ia;
  }
  return 0;
}

50000 × 50000行列のうち、1/2/3/4/5/6/7/8/9/10 %の要素を埋めた疎行列に対して、ベクトル積計算を行い、ベクトルの転送との合計の処理時間(行列の転送は含まない)を測定するプログラムである。mkl_dcsrmvはMKLのSpMVルーチン、cusparseDcsrmvはCUSPARSEのSpMVルーチンである。どちらもCompressed Sparse Row (CSR)形式を使うことができ、相互の置き換えに便利なようにAPIが設計されている。

ベンチマーク環境

ベンチマーク環境は以下のとおりである。

CPU: Xeon E5-2670 v2 (2.5 GHz, 10 Cores) × 2
GPGPU: Tesla K20
MIC: Xeon Phi 7120P
OS: CentOS 6.4
C++: Intel C++ Compiler XE 13.1.3.192.20130607
MKL: Math Kernel Library 11.0
CUDA: CUDA 5.5
NVIDIAドライバ: 319.37

ベンチマーク

要素数 CPU (MFLOPS) GPGPU (MFLOPS) MIC (MFLOPS)
25000000 (1 %) 404.8255 12863.39 129.6317
50000000 (2 %) 5176.787 15573.90 1486.370
75000000 (3 %) 5134.173 16933.85 1589.842
100000000 (4 %) 5075.627 17528.48 1689.146
125000000 (5 %) 4996.502 17857.14 1691.555
150000000 (6 %) 4965.819 18260.39 1746.938
175000000 (7 %) 4957.648 18600.20 1768.096
200000000 (8 %) 4976.300 18920.58 1758.698
225000000 (9 %) 5035.867 19210.25 1804.345
250000000 (10 %) 5066.318 19487.86 1796.132

MKLの結果で最初の1 %目の性能が低いのは、MKLの初期化オーバーヘッドであろうか。ともかく傾向としては、Xeonがおよそ5 GFLOPS、Teslaが10 GFLOPS台後半なのに対して、Xeon Phiが1 GFLOPS後半しか性能を出せていないというものである。


この記事へのコメントをお寄せください

  • サイトへの書き込みに差し支えございましたら トータルディスクロージャーサイトサポート係へメールをお送りください
  • トータル・ディスクロージャ・サイトに投稿された文章と画像は、すべてその著作権がHPCシステムズ株式会社に帰属し、HPCシステムズ株式会社が著作権を所有することに同意してください。
  • あなたの文章が他人によって自由に編集、配布されることを望まない場合は、投稿を控えてください。
  • コメントを書き込む場合は名前にひらがなを織り交ぜてください。
  • あなたの投稿する文章と画像はあなた自身によって書かれたものであるか、パブリック・ドメインかそれに類する自由なリソースからの複製であることを約束してください。あなたが著作権を保持していない作品を許諾なしに投稿してはいけません!

<comments hideform="false" />


Comments

ノート:疎行列-ベクトル積性能評価: CUDA、MIC、そしてCPU

個人用ツール