Elixirで並列レイトレーサーを作った話 その2

こちら、
- Elixirで並列レイトレーサーを作った話 その1

の続きです。今回も長文で分けるか迷いましたが、OpenMP/CLは参考程度なので、特にOpenCLは手順とコードを書いてはいますが、太字の手順だけ流し読みして、どれくらい手間がかかるかだけ確認していただければと思います。

【2020/03/20：注釈】 Pythonよりだいぶ速かったのですが、素の算術演算でなく、numpyでの行列演算になると、Pythonはめっちゃ速くなります。numpyすごい(実際動いてるのCだけど)。

本記事の概要

今回はElixirの並列化とC++のCPU、GPUそれぞれの並列化について簡単に説明した上でベンチマークの結果、という流れになります。ベンチマークに使うレイトレーシングのシーンはこちらです。

これをベンチマークに使ったのは、ある程度時間のかかる処理でないと差がわかりにくいからです。単純な処理だと、特に比較対象のC++でのCPUの並列化とGPUの並列化の差がわからなくなってしまいます。

Elixirでの並列化

自力でspawnするという手もありますが、Flowを使うのが簡単です。
FlowはGenStageをラップした並列化のパッケージで、Enum.mapとかEnum.reduceを置き換えるだけで自動的に並列化を行ってくれます。

Flowによる並列化は
@zacky1972 さんの

• ZEAM開発ログv0.1.1 AI/MLを爆速にしたい！ Flow / GenStage でGPUを駆動できないの？
• ZEAM開発ログv0.1.2 AI/MLを爆速にしたい！ Flow のコードを OpenCL で書いてみる〜CPU編

辺りがわかりやすいと思います。というか、基本的にここでやっていることも大差ありません。
レイトレーサーでCPUバウンドな状況を作ってベンチマーク取ることで、CGでもElixir結構いけるんじゃない？という確認がしたかったこというというのがあります。自分の会社でもElixirと言えばフロントエンドなイメージがあるようでWebの技術カテゴリに入れられてましたが、いやいや、並列化できるならCPUバウンドなプログラムも十分早くなるでしょ、というところを調べたかったというのも自分の中でありました。また、ここから機械学習につなげて、GPUと絡む流れは楽しいなぁ、と。

Flowの導入

リンク貼るだけではちょっと雑なので、簡単に説明します。

まず、Elixir側ですが、依存性にFlowを追加します。

mix.exs

  # Run "mix help deps" to learn about dependencies.
  defp deps do
    [
      {:erlport, "~> 0.10.0" }, 
      {:flow, "~> 0.14.3"},
    ]
  end

いつものようにmix deps.getで依存性を解決します。

EnumをFlowに置き換える

EnumをFlowで置き換えるには、まず、Flowモジュールにリストを変換する必要があります。
これにはリストに対してFlow.from_enumerableを呼べばOKです。

Rangeはそのままではfrom_enumerable呼べないって記事を見た気がしますが、新しめのバージョンだといけますね。

注意点として、Flow.mapは順序の保証がありません。

遅延評価だと思っていたのですが、バックスラッシュで改行した時はちゃんと並ぶようなので、もしかすると裏で動いてる？

重要なパラメータ

上のリンクした @zacky1972 さんの記事でもありますが、max_demandとstagesが重要なパラメータになります。stagesは3個上の画像では自動で4になっており、このMacBook Proの論理コア数は4なので、自動で最適化されてそうです（違ったらごめんなさい)。この辺あまり言及されてない気がします。
max_demandはリスト分割の最大数です。ある程度の塊として並列化したプロセスに渡すので、データが少ない時はこれを設定しないと分割されない可能性があります。前項最後の例だと2個ずつ分割したために順番がバラけているのがわかります。していない場合、1プロセスで処理させて順番そのままだと予想できます。

前回のレイトレーサーを並列処理にしてみます。

raytracer_ex.ex

    scene.ny-1..0
    |> Enum.to_list
    |> Flow.from_enumerable(max_demand: 4)
    |> Flow.map(fn y->
      row = 
      0..scene.nx-1
        |> Enum.map(fn x ->
          color =
          1..scene.ns
          |> Enum.reduce(%Vec3{}, fn (_x, c) ->
            u = (x + :rand.uniform()) / scene.nx
            v = (y + :rand.uniform()) / scene.ny
            ray = Camera.get_ray(scene.camera, u, v)
            Vec3.add(c, _color(ray, scene.objects, 0))
          end)
          Vec3.to_list(Vec3.div(color, scene.ns)) |> Enum.map(&(:math.sqrt(&1) * 255))
        end)
      IO.puts("#{scene.ny - y} / #{scene.ny}")
      {y, row}
      end)
    |> Enum.sort(&(elem(&1, 0) > elem(&2, 0)))
    |> Enum.map(&elem(&1, 1))
    |> show

ここでは行単位で並列化しています。並列化の粒度は実験してみればわかると思いますが、細か過ぎるとデータ分割やプロセス生成のオーバーヘッドで並列化の恩恵が相殺されてしまいます。ご多分に漏れず、今回も列単位で分割した場合は速度はほとんど高速化されませんでした。

Elixirによる並列化はこれだけです。

少し説明すると、{y, row}部分で行数と行の列リストのタプルを作って、Enum.sortでタプルの1個目の行数でソートして、Enum.mapで2個めの要素、列リストを取り出し絵t、行x列（xRGB)の多次元配列にしています。
結果が単純なソートじゃない場合にEnum.sortを使う際はEnum.with_indexを使ったりする方法もありますが、行数(ピクセルのy要素)でループをしてる場合、それをそのまま使えます。

これはElixirによるCPU並列化です。
次に比較対象のC++でのCPU/GPU並列化を簡単に説明します。

C/C++でCPU並列化

CやC++でCPUの並列化をするにはOpenMP(Multi Processing)を使うのが楽です。

OpneMPの利用

インストールは環境に合わせてして下さい。
あとは

src/raytracer.cpp

#include <omp.h>

でincludeして、

src/raytracer.cpp

  #pragma omp parallel for schedule(static, 1)
  for (int y = ny-1; y >= 0; --y) {

for文の前にディレクティブ#pragma omp parallel forを書いてあげればOKです。
for文以外の並列化の手法もありますが、今回はfor文の並列化のみです。scheduleに関してはなくても平気です。1つ目の引数がstaticかdynamic、2つ目でデータの割当数(チャンクサイズ）を指定します。staticは分割時にどのスレッドがどのデータを処理するか決まり、dynamicの場合、終わったスレッドに次のチャンクを割り当てます。割当のオーバーヘッドがあるので、データ毎の処理時間がほぼ同じ場合、staticの方がパフォーマンスが出る傾向にあります。scheduleはオプションなので書かなくても大丈夫です。

また、並列数を指定する場合、

src/raytracer.cpp

#ifdef _OPENMP   
  int num_processors = omp_get_num_procs( );
  omp_set_num_threads(num_processors);
#endif

のようにomp_set_num_threads()で指定します。
OpenMPが使える場合、_OPENMPが宣言されるので上記のようにプリプロセッサでインストールされてる時のみ並列化、みたいなこともできます。上記ではomp_get_num_procs()で論理スレッド数を取得して設定しています。

C/C++でGPUによる並列化

GPUでの並列化の一般的な目的での利用(GPGPU)ではOpenCLがあります。

OpenCLの導入

インストールは環境ごとに調べて実行して下さい。
パスも環境で違ったりするので、ご注意下さい。

src/raytracer.cpp

#include <OpenCL/cl.h>

簡単に説明しますが、これくらい手間がかかる、という意味で書くので、流し読みして下さい。省略していますが、本来、各処理の後にstatusをチェックして処理が成功しているか確認すべきです。
1. platformの取得

src/raytracer.cpp

  cl_int status;

  cl_platform_id platform;
  status = clGetPlatformIDs(1, &platform, NULL); 

2. deviceの取得

src/raytracer.cpp

  cl_device_id device;
  status = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);

3. CPU側のバッファの準備

src/raytracer.cpp

  float *h_results = new float[kColorSize];

4. contextの作成

src/raytracer.cpp

  cl_context context = clCreateContext(NULL, 1, &device, NULL, NULL, &status);

5. command queueの生成

src/raytracer.cpp

  cl_command_queue cmd_queue = clCreateCommandQueue(context, device, 0, &status);

6. GPU側バッファの準備

src/raytracer.cpp

  cl_mem d_results = clCreateBuffer(context, CL_MEM_WRITE_ONLY, data_size, NULL, 

7. CPU側のデータをGPU側のバッファに書き込む(引数渡す場合)

src/raytracer.cpp

  status = clEnqueueWriteBuffer(cmd_queue, d_spheres, CL_FALSE, 0, sphere_size, h_spheres, 0, NULL, NULL);

8. kernel(GPUで実行する命令)の文書読み込み

src/raytracer.cpp

  cl_file.seekg (0, cl_file.end);
  int length = cl_file.tellg();
  cl_file.seekg (0, cl_file.beg);

  char *cl_program_text = new char[length + 1];
  cl_file.read(cl_program_text, length);
  cl_program_text[length] = '\0';

  cl_file.close();

  char *strings[1];
  strings[0] = cl_program_text;

9. programの作成及びビルド

src/raytracer.cpp

  cl_program program = clCreateProgramWithSource(context, 1, (const char**)strings, NULL, &status);
  delete[] cl_program_text;

  // compile and link the kernel code
  char *options = {""};
  status = clBuildProgram(program, 1, &device, options, NULL, NULL);

10. Kernelの生成

src/raytracer.cpp

  // create kernel object
  cl_kernel kernel = clCreateKernel(program, args[1], &status);

11. 引数の設定

src/raytracer.cpp

  status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_results);

12. 処理の実行

src/raytracer.cpp

  status = clEnqueueNDRangeKernel(cmd_queue, kernel, 3, NULL, global_worker_size, NULL, 0, NULL, NULL);

13. 結果の読み込み

src/raytracer.cpp

  status = clEnqueueReadBuffer(cmd_queue, d_results, CL_TRUE, 0, data_size, h_results, 0, NULL, NULL);

14. 後処理

src/raytracer.cpp

  clReleaseKernel(kernel);
  clReleaseProgram(program);
  clReleaseCommandQueue(cmd_queue);
  clReleaseMemObject(d_results);
  clReleaseMemObject(d_spheres);

  delete[] h_spheres;

みたいな流れです。カーネル（GPUでの処理）の文書はこんな感じです。

cl/raytracer.cl

kernel
void Draw(global float *d_results, constant Sphere *s, int n_sphere, constant int *types,
                constant float4 *materials, constant float3 *uvs)
{
  const int x = get_global_id(0);
  const int y = get_global_id(1);
  const int sn = get_global_id(2);

  const int idx = ((y * get_global_size(0) + x) * get_global_size(2) + sn) * 3;

  float3 lookfrom = (float3)(11.0, 2.0, 3.0);
  float3 lookat = (float3)(0.0, 0.0, 0.0);
  float dist_to_focus = fast_length(lookfrom - lookat);
  float aperture = 0.1;

  const Camera camera = SetupCamera(lookfrom, lookat, (float3)(0.0, 1.0, 0.0), 20.0,
                                    float(get_global_size(0))/ float(get_global_size(1)), aperture, dist_to_focus);

  const float u = float(x + uvs[sn].x) / get_global_size(0);
  const float v = float(y + uvs[sn].y) / get_global_size(1);

  Ray ray = GetRay(camera, u, v, uvs);

  float3 color = Color(ray, s, n_sphere, types, materials, uvs, 0);

  d_results[idx] = color.x;
  d_results[idx + 1] = color.y;
  d_results[idx + 2] = color.z;
}

関数切って意外といろんなことができますが、ランダム性がない(同時に実行されるため？seedも効かず、drand48()を使っても1より大きい数が返ってきたり、ちょっと普通のC/C++とOpenCL C/C++は違うみたいです。今回はCPU側からランダムの配列を渡して利用しました。

性能比較結果

これにpythonでのプログラムを含めての性能比較結果は

Elixir	Elixir	Python	C++	C++	C++
非並列	4並列(Flow)	非並列	非並列	CPU並列	GPU並列(OpenMP)
1,949,963	1,087,444	9,976,854	43,587	25,479	767
1.79	1.0	9.17	0.035	0.02	0.0007

マシン性能
MacBook Pro2017
OSX10.13.6 High Sierra
2.3GHz Intel Core i5(2core + Hyper Threading=4)
16GB LPDDR3

使ったコード（github)は下記になります。READMEとかテストちゃんと書いてないので、時間ある時に書きますが…C++はautotoolsを使ってます。

• Elixir
• Python
• C++(OpenMP/OpenCL含む)

まとめ

• Elixirの論理コア数での並列はPythonの9倍強速かった
• C++はElixirの約30倍速かった
• C++のCPU並列化はElixir並列の約50倍速かった
• C++のGPU並列化はめちゃくちゃ速かった

Elixirを加速させるための今後の展望

性能比較の通り、ElixirはPythonのようなハイレベル言語にはかなりのアドバンテージで速いと言えそうですが、C/C++などのローレベル言語にはまだまだ太刀打ちできません。ただ、ElixirからGPUでの処理が可能になればかなりのスピードアップが期待できると言えると思います。
参加させていただいているfukuoka.exではHastega (ヘイスガ)という名前のElixirでのGPGPUでの処理高速化に向けた取り組みを行っています(下記をご参照下さい)。
WebGL / WebGPU + Hastega / Elixir / Phoenix で分散／エッジ・コンピューティング
上記でわかる通り、OpenCLでのGPGPUの利用は結構なコード量であり、OpenCL C/C++でのカーネルの記述は処理上の制限があったりして、結構使い勝手が悪いのですが、Hastegaでは簡単な記述でのGPGPUが期待でき、簡潔で高速なプログラムを実現できるようになる、はずです。

私の本職は自然言語処理のエンジニアでCGやDeep Learningによる画像認識は趣味なのですが、これらのCPUバウンドな分野にもHastegaによるGPUでの加速が伴えばElixirは十分勝負できると思いますし、そこに貢献していきたいと思っています。今年はその辺りに注力していければ。

長文にお付き合いいただきありがとうございました。