試行錯誤な日々: OpenCLで行列の掛け算をしてみた

背景

OpenCLとは、GPUを利用して並列処理プログラムを作れるライブラリです。
そのOpenCLを利用して、行列の掛け算をしてみました。

備忘録を兼ねて、自分が取り組んだ内容を共有します。

全体像

使ったもの
計算内容
プログラムの説明
行列を用意
OpenCLの変数を用意
デバイスIDを取得
並列処理（カーネル）を記述
並列処理（カーネル）を作成
プログラムの用意
メモリの確保
並列数の設定
実行
結果の取得
リソースの開放
比較のために、並列処理しないプログラムを作成
行列の大きを変えて実行速度を比較
まとめ

使ったもの

OpenCLのプログラムがコンパイルできる環境を使います。

下記の記事で作った環境を、この記事では利用しました。
UbuntuでIntelのCPUに内蔵されているGPUをOpenCLで使う方法

計算内容

幅aW、高さbWの行列Aと、幅bW。高さaWの行列Bの積を計算して、幅aW、高さbWの行列Rを取得します。

A x B = R

説明するプログラムは下記のリポジトリで管理しています。
asukiaaa/c_opencl_practice

下記のようなコマンドでダウンロードして実行できます。

git clone https://github.com/asukiaaa/c_opencl_practice.git
cd c_opencl_practice
gcc matrix_dot_matrix.c -lOpenCL
./a.out

行列を用意

行列Aの幅（wA）と高（hA）、行列Bの幅（wB）を指定して、行列A、Bと結果を入れるための行列Rをそれぞれ作成します。
wA、hA、wBを変えるたびにコンパイルし直すのは面倒なので、実行時の引数として指定できるようにしました。

計算に使うAとBは、下記のルールで値を設定しました。
行列A: A[x][y] = x*10 + y
行列B: xとyが1のところだけ2の単位行列
（行列Bを単位行列にしなかったのは、期待通りに行列計算できていることを確認するためです。）

int main(int argc, char *argv[]) {
  int wA, hA, wB;
  if (argc == 1) {
    wA = hA = wB = 10;
  } else if (argc == 2) {
    wA = hA = wB = atoi(argv[1]);
  } else if (argc == 4) {
    wA = atoi(argv[1]);
    hA = atoi(argv[2]);
    wB = atoi(argv[3]);
  } else {
    printf("No value or 1 or 3 values for wA, hA and wB are expected.\n");
    return EXIT_FAILURE;
  }
  int hB = wA, wR = wB, hR = hA;
  unsigned int matrixAMemSize = sizeof(float) * (unsigned int) (wA * hA);
  unsigned int matrixBMemSize = sizeof(float) * (unsigned int) (wB * hB);
  unsigned int matrixRMemSize = sizeof(float) * (unsigned int) (wR * hR);
  float* matrixA = (float*) malloc(matrixAMemSize);
  float* matrixB = (float*) malloc(matrixBMemSize);
  float* matrixR = (float*) malloc(matrixRMemSize);

  int i, j;
  for (i = 0; i < wA; i++) {
    for (j = 0; j < hA; j++) {
      matrixA[j*wA + i] = 10*i + j;
    }
  }
  for (i = 0; i < wB; i++) {
    for (j = 0; j < hB; j++) {
      int bValue = 0;
      if (i == j) {
        if (i == 1) {
          bValue = 2;
        } else {
          bValue = 1;
        }
      }
      matrixB[j*wB + i] = bValue;
    }
  }

  // ..
}

行列の内容を見たいので、出力関数も作りました。
行列が大きいと表示に時間がかかる上に人による読み取りが難しくなるので、リポジトリのプログラムはwAが10以下の時だけ行列の内容を出力するようにしました。

void printMatrix(float* matrix, int w, int h) {
  int i, j;
  printf("matrix\n");
  for (i = 0; i < h; i++) {
    for (j = 0; j < w; j++) {
      printf(" %6.2f", matrix[i*w + j]);
    }
    printf("\n");
  }
}

wA = 3、hA = 4、wB = 4のときの行列Aと行列Bは、このようになります。

行列A

printMatrix(matrixA, wA, hA);

matrix
   0.00  10.00  20.00
   1.00  11.00  21.00
   2.00  12.00  22.00
   3.00  13.00  23.00

行列B

printMatrix(matrixB, wB, hB);

matrix
   1.00   0.00   0.00   0.00
   0.00   2.00   0.00   0.00
   0.00   0.00   1.00   0.00

OpenCLの変数を用意

変数を準備しておきます。

  cl_device_id device_id = NULL;
  cl_context context = NULL;
  cl_command_queue command_queue = NULL;
  cl_mem matrixAMemObj = NULL, matrixBMemObj = NULL, matrixRMemObj = NULL;
  cl_program program = NULL;
  cl_kernel kernel = NULL;
  cl_platform_id platform_id = NULL;
  cl_uint ret_num_devices;
  cl_uint ret_num_platforms;
  cl_int ret;

ちなみに、retに関数の戻り値を保存する記述が多々あります。
retにはエラーを意味する数字が入ることがありまして、下記のヘッダファイルを見ると、その数字が何を意味するか分かります。

KhronosGroup/OpenCL-Headers/opencl22/CL/cl.h

デバイスIDを取得

OpenCLのデバイスIDを取得します。
ここで取得したdevice_idを、メモリの操作やプログラムを実行する関数の引数に利用します。

  /* Get Platform and Device Info */
  ret = clGetPlatformIDs(1, &platform_id, &ret_num_platforms);
  ret = clGetDeviceIDs(platform_id, CL_DEVICE_TYPE_DEFAULT, 1, &device_id, &ret_num_devices);
  if (ret != CL_SUCCESS) {
    printf("Error: Failed to create a device group!\n");
    return EXIT_FAILURE;
  }

並列処理（カーネル）を記述

カーネルの記述は、ファイルで管理する方法とcプログラムに文字列として記述する方法があります。
一長一短ですので、お好みの方法を選択してください。

ファイルで管理
利点: clファイルに対してエディタのシンタックスが効くなど、コードをきれいに管理できる
欠点: コンパイルしても実行ファイルと共にカーネルファイルも管理する必要がある

コードの中の文字列として管理
利点: コンパイル後は実行ファイルだけを管理すればよいという利点があります。
欠点: エディタのシンタックスが反映されないなど、カーネル記述の見た目が悪くなる

リポジトリのプログラムは、最初はファイルで管理していましたが、実行ファイルだけで動いてほしくなったので文字列での管理に変更しました。

ファイルで管理する方法

行列を処理するためのOpenCVのカーネル関数をファイルに定義します。

matrix_dot_matrix.cl

__kernel void matrix_dot_matrix(
  __global float *A,
  __global float *B,
  __global float *Result,
  int wA,
  int wB) {
  int tx = get_global_id(0);
  int ty = get_global_id(1);

  float value = 0;
  for (int k = 0; k < wA; ++k) {
    float elementA = A[ty * wA + k];
    float elementB = B[k * wB + tx];
    value += elementA * elementB;
  }

  Result[ty * wB + tx] = value;
}

上記のファイルを開き、文字列として保持します。

  FILE *fp;
  char fileName[] = "./matrix_dot_matrix.cl";
  char *source_str;
  size_t source_size;
  fp = fopen(fileName, "r");
  if (!fp) {
    fprintf(stderr, "Failed to load kernel.\n");
    exit(1);
  }
  source_str = (char*)malloc(MAX_SOURCE_SIZE);
  source_size = fread(source_str, 1, MAX_SOURCE_SIZE, fp);
  fclose(fp);

コード中の文字列として管理する方法

文字列として持ち、大きさを計ります。

  char *source_str =
    "__kernel void matrix_dot_matrix( \
       __global float *A, \
       __global float *B, \
       __global float *Result, \
       int wA, \
       int wB) { \
       int tx = get_global_id(0); \
       int ty = get_global_id(1); \
       float value = 0; \
       for (int k = 0; k < wA; ++k) { \
         float elementA = A[ty * wA + k]; \
         float elementB = B[k * wB + tx]; \
         value += elementA * elementB; \
       } \
       Result[ty * wB + tx] = value; \
     }";
  size_t source_size = strlen(source_str);

並列処理（カーネル）を作成

カーネルの文字列（source_str）と大きさ（source_size）、中間生成物を作りつつ、clCreateKernelでカーネルを作ります。
contextやcommand_queueなどの中間生成物が何に使えるのかはまだ把握していませんが、この記述で動作しました。
ここで作ったカーネルを呼び出すことで、並列処理を実行できます。

  /* Create OpenCL context */
  context = clCreateContext(NULL, 1, &device_id, NULL, NULL, &ret);

  /* Create Command Queue */
  command_queue = clCreateCommandQueueWithProperties(context, device_id, 0, &ret);

  /* Create Kernel Program from the source */
  program = clCreateProgramWithSource(context, 1, (const char **)&source_str, (const size_t *)&source_size, &ret);

  /* Build Kernel Program */
  ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);
  if (ret != CL_SUCCESS) {
    size_t len;
    char buffer[2048];
    printf("Error: Failed to build program executable!\n");
    clGetProgramBuildInfo(program, device_id, CL_PROGRAM_BUILD_LOG, sizeof(buffer), buffer, &len);
    printf("%s\n", buffer);
    exit(1);
  }

  /* Create OpenCL Kernel */
  kernel = clCreateKernel(program, "matrix_dot_matrix", &ret);
  if (!kernel || ret != CL_SUCCESS) {
    printf("Error: Failed to create compute kernel!\n");
    exit(1);
  }

メモリの確保

clCreateBufferでOpenCL用のバッファを確保します。
行列A、BからをREAD_ONLYなバッファを作りました。
行列Rは出力として使うので、WRITE_ONLYとして、バッファのザイズを渡しました。

  /* Create Memory Buffer */
  matrixAMemObj = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, matrixAMemSize, matrixA, &ret);
  matrixBMemObj = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, matrixBMemSize, matrixB, &ret);
  matrixRMemObj = clCreateBuffer(context, CL_MEM_WRITE_ONLY, matrixRMemSize, NULL, &ret);
  if (!matrixAMemObj || !matrixBMemObj || !matrixRMemObj) {
    printf("Error: Failed to allocate device memory!\n");
    exit(1);
  }

  /* Set OpenCL Kernel Parameters */
  ret = clSetKernelArg(kernel, 0, sizeof(cl_mem), (void *)&matrixAMemObj);
  ret |= clSetKernelArg(kernel, 1, sizeof(cl_mem), (void *)&matrixBMemObj);
  ret |= clSetKernelArg(kernel, 2, sizeof(cl_mem), (void *)&matrixRMemObj);
  ret |= clSetKernelArg(kernel, 3, sizeof(int), (void *)&wA);
  ret |= clSetKernelArg(kernel, 4, sizeof(int), (void *)&wB);
  if (ret != CL_SUCCESS) {
    printf("Error: Failed to set kernel arguments! %d\n", ret);
    exit(1);
  }

並列数の設定

個人的に最も時間を取られた設定です。
まずは全体像をご覧ください。

int getMaxCommonFactorOf2Pow(int target) {
  int commonFactor = 1;
  while (target % 2 == 0) {
    target /= 2;
    commonFactor *= 2;
  }
  return commonFactor;
}

  int workDim = 2;
  size_t globalWorkSize[workDim], localWorkSize[workDim];
  globalWorkSize[0] = wR;
  globalWorkSize[1] = hR;
  size_t localWR = wR, localHR = hR;
  while (maxLocalSizes[0] < localWR ||
         maxLocalSizes[1] < localHR ||
         localWR * localHR > maxWorkGroupSize) {
    if (maxLocalSizes[0] < localWR) {
      if (localWR % 2 == 0) {
        localWR /= 2;
      } else {
        localWR = getMaxCommonFactorOf2Pow(localWR);
      }
    } else if (maxLocalSizes[1] < localHR) {
      if (localWR % 2 == 0) {
        localHR /= 2;
      } else {
        localHR = getMaxCommonFactorOf2Pow(localHR);
      }
    } else if (localWR > localHR && localWR % 2 == 0) {
      localWR /= 2;
    } else if (localHR % 2 == 0) {
      localHR /= 2;
    } else if (localWR % 2 == 0) {
      localWR /= 2;
    } else if (localHR != 1) {
      localHR = getMaxCommonFactorOf2Pow(localHR);
    } else {
      localWR = getMaxCommonFactorOf2Pow(localWR);
    }
  }
  localWorkSize[0] = localWR;
  localWorkSize[1] = localHR;

workDim、globalWorkSize、localWorkSizeは、この後の「並列処理の実行」で使うclEnqueueNDRangeKernelという関数に渡すための値です。
globalWorkSizeは行いたい処理の総数、localWorkSizeは並列に行いたい処理の数をそれぞれ意味します。

globalWorkSizeとlocalWorkSizeには下記のような制約があるため、それを満たすために上記のような処理を行っています。

localWorkSizeの要素は、2の乗数（1,2,4,8..）をかけると、対応するglobalWorkSizeになること
例: localWorkSize[0] * pow(2,n) == globalWorkSize[0]
localWorkSizeの要素をかけ合わせた合計は、CL_DEVICE_MAX_WORK_GROUP_SIZEで取得する値以下になること
localWorkSizeはCL_DEVICE_MAX_WORK_ITEM_DIMENSIONSで取得した値以下になること

上記の制約があるため、計算対象となる行列によっては、並列計算できない場合があります。

例: CL_DEVICE_MAX_WORK_GROUP_SIZEが8192に対して8209x8209の行列を計算する場合。
8209は2で割ることができないため1x1が実行可能な要素の組み合わせになります。
この場合、処理が1つずつしか走らないので、GPUで実行する意味が無くなってしまいます。
（1x1の場合、CPUより処理に時間がかかる可能性が高いです。）

次元あたりの並列処理できる数が1にならないように、globalWorkSize（行列の幅や高さ）は2の倍数で設定するのが良さそうです。

並列処理を実行

これまでに用意した変数を利用して、下記の処理で並列処理を実行できます。

  ret = clEnqueueNDRangeKernel(command_queue, kernel, workDim, NULL, globalWorkSize, localWorkSize, 0, NULL, NULL);

結果の取得

計算した結果をmatrixRに配置します。

  ret = clEnqueueReadBuffer(command_queue, matrixRMemObj, CL_TRUE, 0, matrixRMemSize, matrixR, 0, NULL, NULL);

このような結果になりました。

printMatrix(matrixR, wR, hR);

matrix
   0.00  20.00  20.00   0.00
   1.00  22.00  21.00   0.00
   2.00  24.00  22.00   0.00
   3.00  26.00  23.00   0.00

期待通りに行列の掛け算を行えました。

リソースの開放

OpenCLのために確保した領域や、mallocした変数を開放します。

  ret = clFlush(command_queue);
  ret = clFinish(command_queue);
  ret = clReleaseKernel(kernel);
  ret = clReleaseProgram(program);
  ret = clReleaseMemObject(matrixAMemObj);
  ret = clReleaseMemObject(matrixBMemObj);
  ret = clReleaseMemObject(matrixRMemObj);
  ret = clReleaseCommandQueue(command_queue);
  ret = clReleaseContext(context);

  free(matrixA);
  free(matrixB);
  free(matrixR);
  free(maxLocalSizes);

比較のために、並列処理しないプログラムを作成

OpenCLで作ったプログラムがどのくらい早いのか比較するために、並列処理しないプログラムを作成しました。

matrix_dot_matrix_on_cpu.c

...
  for (i = 0; i < wA; i++) {
    for (j = 0; j < wB; j++) {
      float value = 0;
      for (k = 0; k < wA; k++) {
        float elementA = matrixA[j * wA + k];
        float elementB = matrixB[k * wB + i];
        value += elementA * elementB;
      }
      matrixR[j*wB + i] = value;
    }
  }
...

行列の大きさを変えて実行速度を比較

下記の環境で、実行速度を比較してみました。

分類	名前
モデル	Lenovo ThinkPad S1 Yoga
CPU	Intel Core i7-4600U CPU @ 2.10GHz × 4
GPU	Intel Haswell Mobile
OS	Ubuntu17.10

リポジトリにあるプログラムを利用する場合、このように実行できます。

プログラムをダウンロードして、その中に移動

git clone https://github.com/asukiaaa/c_opencl_practice.git
cd c_opencl_practice

2000x2000の行列の並列演算プログラムの実行

gcc matrix_dot_matrix.c -lOpenCL
./a.out 2000

2000x2000の行列の並列ではない演算プログラムの実行

gcc matrix_dot_matrix_on_cpu.c
./a.out 2000

結果はこのようになりました。

Program	2000x2000	4000x4000	10000x10000
並列処理しない版	62.44 秒	594.91 秒	4時間近くかかりそうなので試していません
並列処理版	18.72 秒	203.47 秒	4436.74 秒

並列処理すると3倍くらい早くなるようでした。
早くはなったのですが、10倍くらい早くなることを期待していたので、期待外れな結果になりました。

別の記事で実装したブロックマッチングはOpenCLを利用すると100倍位早くなったので、推測ですが、IntelのCPUに載っているGPUは掛け算や割り算の計算時間が長いのかもしれません。

まとめ

localWorkSizeを求めるのに時間を取られましたが、OpenCLを利用して行列の掛け算を並列処理できました。

何かの参考になれば嬉しいです。

参考

3.3 First OpenCL Program
How to compile OpenCL example in GCC?
Matrix multiplication in OpenCL

変更履歴

2018.02.18
コンパイルエラーを引き起こしていたconstを削除しました。
カーネルを文字列として管理する方法を追加しました。

2018.03.02
clCreateBufferで設定するメモリタイプを最適化しました。
（CL_MEM_READ_WRITEをCL_MEM_READ_ONLYかCL_MEM_WRITE_ONLYに変更しました。)

2018.03.07
ブロックマッチングでは100倍早くなったというコメントを追加しました。

試行錯誤な日々

ページ

2018年2月18日日曜日

OpenCLで行列の掛け算をしてみた

背景

全体像