ZEAM開発ログv0.1.3 AI/MLを爆速にしたい！ Flow のコードを OpenCL で書いてみる〜GPU編

（この記事は、「Elixir or Phoenix Advent Calendar 2017」の21日目です）

昨日は @twinbee さんの「Elixirから簡単にRustを呼び出せるRustler #1 準備編」でしたね。

「ZEAM開発ログ 目次」はこちら

おしらせ その1

「fukuoka.ex#11：DB／データサイエンスにコネクトするElixir」を6/22（金）19時に開催します！

私も「ZEAM（Zacky's Elixir Abstract Machine）開発ログ
　第3回〜AI/MLを爆速に！ElixirでGPUを駆動すると驚きのパフォーマンスに！」というタイトルでLTします。今，Qiitaで連載中のこの記事のまとめと今後の展望について熱く語ります！

おしらせ その2

第2回Web System Architecture研究会で論文発表しました！

タイトルは「Elixirの軽量コールバックスレッドの実装とPhoenixの同時セッション最大数・レイテンシ改善の構想」です。こちらもZEAMにつながるもう一つ別のメカニズム「軽量コールバックスレッド」とPhoenixへの応用について熱く語りました。

• 論文はこちら
• プレゼンテーションはこちら

さて本題〜はじめに

本連載の前回記事はこちら
|> ZEAM開発ログv0.1.0 Flow / GenStage による並列プログラミング入門
|> ZEAM開発ログv0.1.1 AI/MLを爆速にしたい！ Flow / GenStage でGPUを駆動できないの？
|> ZEAM開発ログv0.1.2 AI/MLを爆速にしたい！ Flow のコードを OpenCL で書いてみる〜CPU編

これらを一気にまとめると次のようになります。

(v0.1.0) Elixir Flow に沿ってプログラムを書くと，いい感じでマルチコア並列にしてくれます。stages を変化させたときに，マシンの論理コア数までは，コア数が増えるごとに実行速度が向上し，論理コア数以上にしたときには，実行速度が落ちるので，stagesを論理コア数と等しくするのが最もパフォーマンスが良くなります。これを裏付けるデータはこちら

stages(数)	benchmark(秒)	benchmark2(秒)	benchmark3(秒)	備考
1	52.795620	54.697525	44.455119
2	24.716176	25.297751	20.675610
4	15.016131	15.763084	13.610704	物理コア数
8	12.664873	13.366235	11.308742	HT込みコア数，最速
16	12.807277	13.611112	11.411827
32	12.841774	14.007026	11.714803
64	13.158978	14.013323	11.898896
128	13.217850	13.422258	11.914322
備考	mapの中で再帰ループ	mapを展開	mapの中に展開

(v0.1.1) Elixir Flow でコードはGPU向きのコードになっています。つまり，単純な構造で均質で大量にあるデータを，ほぼ同じような命令列で処理するコードになっています。

典型的な Flow のコードは下記の通りです。

list
|> Flow.from_enumerable
|> Flow.map(foo)
|> Flow.map(bar)
|> Flow.map(hoge)
|> Enum.to_list

list が単純な構造で均質で大量にあるデータであるリスト構造で，foo |> bar |> hogeという同じ命令列で処理します。

(v0.1.1) 並列化しないC言語とElixirで比較すると実行速度が12倍C言語の方が速いので，そもそも並列化以前にElixirの処理系には最適化の余地があります。

根拠となるデータは次の通りです。

Elixir(秒)	C言語(秒)
52.795620	4.232451

(v0.1.2) C言語ベースのOpenCL(CPU利用)をすることで，Elixirの同等プログラムと比べて7.64〜8.47倍の速度向上になりました。

根拠となるデータは次の通りです。

Elixir(秒)	Elixir(秒)	Elixir(秒)	C言語(秒)	OpenCL(秒)	OpenCL(秒)
1並列ループ	8並列ループ	8並列インライン展開	1並列ループ	8並列ループ	8並列インライン展開
52.795620	12.664873	11.308742	4.232451	1.496656	1.483530

(v0.1.2) OpenCLでCPUを駆動するよりも，じかに並列プログラミングした方がより高速になる可能性があります。根拠としては，OpenCLで駆動する時に必要となる配列の設定に多くの時間がかかっていたためです。

以上を踏まえると，ElixirをGPU対応すること，およびC言語レベルまで処理系を最適化することで，Elixirの性能が飛躍的に向上する期待があるというわけです。というわけで，今回こそは，いよいよGPUを駆動するプログラムを書いてみたいと思います！

OpenCLによるロジスティック写像のベンチマークプログラム(GPU対応版)

OpenCLでGPUを駆動するロジスティック写像のベンチマークプログラムを GitHub で公開しています。

https://github.com/zeam-vm/logistic_map_OpenCL

Mac の Xcode で動作します。ソースコードは次の通りです。

kernel.cl

#define INLINE

#ifdef INLINE
kernel void logisticsmap(
                         global int *mu,
                         global int *p,
                         global int *input,
                         global int *output)
{
    int x, m, pp;
    size_t i = get_global_id(0);
    x = input[i];
    m = mu[i];
    pp = p[i];
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    x = m * x * (x + 1) % pp;
    output[i] = x;
}
#else
kernel void logisticsmap(
                         global int *mu,
                         global int *p,
                         global int *input,
                         global int *output)
{
    size_t i = get_global_id(0);
    output[i] = input[i];
    for(int n = 0; n < 10; n++) {
        output[i] = mu[i] * output[i] * (output[i] + 1) % p[i];
    }
}
#endif

main.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

// This include pulls in everything you need to develop with OpenCL in OS X.
#include <OpenCL/opencl.h>

// Include the header file generated by Xcode.  This header file contains the
//  kernel block declaration.                                             // 1
#include "kernel.cl.h"

#define LOOP 10
#define P 6700417
#define MU 22

#define GPU

// Hard-coded number of values to test, for convenience.
#define NUM_VALUES 0x2000000


static int logisticsmap_calc(int x, int p, int mu) {
    return mu * x * (x + 1) % p;
}

static int logisticsmap_loopCalc(int num, int x, int p, int mu) {
    for(int i = 0; i < num; i++) {
        x = logisticsmap_calc(x, p, mu);
    }
    return x;
}

// A utility function that checks that our kernel execution performs the
// requested work over the entire range of data.
static int validate(cl_int *input, cl_int *output) {
    int i;
    for (i = 0; i < NUM_VALUES; i++) {
        int expected = logisticsmap_loopCalc(LOOP, input[i], P, MU);
        if ( output[i] != expected) {
            fprintf(stdout,
                    "Error: Element %d did not match expected output.\n", i);
            fprintf(stdout,
                    "       Saw %d, expected %d\n", output[i], expected);
            fflush(stdout);
            return 0;
        }
    }
    return 1;
}

int main (int argc, const char * argv[]) {
    int i;
    char name[128];

    struct timeval start_time;
    gettimeofday(&start_time, NULL);

    // First, try to obtain a dispatch queue that can send work to the
    // GPU in our system.                                             // 2
    dispatch_queue_t queue = NULL;

#ifdef GPU
    queue = gcl_create_dispatch_queue(CL_DEVICE_TYPE_GPU, NULL);
#endif

    // In the event that our system does NOT have an OpenCL-compatible GPU,
    // we can use the OpenCL CPU compute device instead.
    if (queue == NULL) {
        queue = gcl_create_dispatch_queue(CL_DEVICE_TYPE_CPU, NULL);
    }

    // This is not required, but let's print out the name of the device
    // we are using to do work.  We could use the same function,
    // clGetDeviceInfo, to obtain all manner of information about the device.
    cl_device_id gpu = gcl_get_device_id_with_dispatch_queue(queue);
    clGetDeviceInfo(gpu, CL_DEVICE_NAME, 128, name, NULL);
    fprintf(stdout, "Created a dispatch queue using the %s\n", name);

    struct timeval device_setting_time;
    gettimeofday(&device_setting_time, NULL);

    // Here we hardcode some test data.
    // Normally, when this application is running for real, data would come from
    // some REAL source, such as a camera, a sensor, or some compiled collection
    // of statistics—it just depends on the problem you want to solve.
    int *test_in = (int *)malloc(sizeof(cl_int) * NUM_VALUES);
    for (i = 0; i < NUM_VALUES; i++) {
        test_in[i] = (cl_int)i;
    }

    // Once the computation using CL is done, will have to read the results
    // back into our application's memory space.  Allocate some space for that.
    int *test_out = (int *)malloc(sizeof(cl_int) * NUM_VALUES);

    int *test_mu = (int *)malloc(sizeof(cl_int) * NUM_VALUES);
    for (i = 0; i < NUM_VALUES; i++) {
        test_mu[i] = (cl_int)MU;
    }

    int *test_p = (int *)malloc(sizeof(cl_int) * NUM_VALUES);
    for (i = 0; i < NUM_VALUES; i++) {
        test_p[i] = (cl_int)P;
    }

    // The test kernel takes two parameters: an input float array and an
    // output float array.  We can't send the application's buffers above, since
    // our CL device operates on its own memory space.  Therefore, we allocate
    // OpenCL memory for doing the work.  Notice that for the input array,
    // we specify CL_MEM_COPY_HOST_PTR and provide the fake input data we
    // created above.  This tells OpenCL to copy the data into its memory
    // space before it executes the kernel.                               // 3
    void* mem_mu  = gcl_malloc(sizeof(cl_int) * NUM_VALUES, test_mu,
                               CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR);
    void* mem_p  = gcl_malloc(sizeof(cl_int) * NUM_VALUES, test_p,
                               CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR);
    void* mem_in  = gcl_malloc(sizeof(cl_int) * NUM_VALUES, test_in,
                               CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR);

    // The output array is not initalized; we're going to fill it up when
    // we execute our kernel.                                             // 4
    void* mem_out =
    gcl_malloc(sizeof(cl_int) * NUM_VALUES, NULL, CL_MEM_WRITE_ONLY);

    struct timeval array_setting_time;
    gettimeofday(&array_setting_time, NULL);

    // Dispatch the kernel block using one of the dispatch_ commands and the
    // queue created earlier.                                            // 5

    dispatch_sync(queue, ^{
        // Although we could pass NULL as the workgroup size, which would tell
        // OpenCL to pick the one it thinks is best, we can also ask
        // OpenCL for the suggested size, and pass it ourselves.
        size_t wgs;
        gcl_get_kernel_block_workgroup_info(logisticsmap_kernel,
                                            CL_KERNEL_WORK_GROUP_SIZE,
                                            sizeof(wgs), &wgs, NULL);

        // The N-Dimensional Range over which we'd like to execute our
        // kernel.  In this case, we're operating on a 1D buffer, so
        // it makes sense that the range is 1D.
        cl_ndrange range = {                                              // 6
            1,                     // The number of dimensions to use.

            {0, 0, 0},             // The offset in each dimension.  To specify
            // that all the data is processed, this is 0
            // in the test case.                   // 7

            {NUM_VALUES, 0, 0},    // The global range—this is how many items
            // IN TOTAL in each dimension you want to
            // process.

            {wgs, 0, 0}            // The local size of each workgroup.  This
            // determines the number of work items per
            // workgroup.  It indirectly affects the
            // number of workgroups, since the global
            // size / local size yields the number of
            // workgroups.  In this test case, there are
            // NUM_VALUE / wgs workgroups.
        };
        // Calling the kernel is easy; simply call it like a function,
        // passing the ndrange as the first parameter, followed by the expected
        // kernel parameters.  Note that we case the 'void*' here to the
        // expected OpenCL types.  Remember, a 'float' in the
        // kernel, is a 'cl_float' from the application's perspective.   // 8

        logisticsmap_kernel(&range,(cl_int*)mem_mu, (cl_int*)mem_p, (cl_int*)mem_in, (cl_int*)mem_out);

        // Getting data out of the device's memory space is also easy;
        // use gcl_memcpy.  In this case, gcl_memcpy takes the output
        // computed by the kernel and copies it over to the
        // application's memory space.                                   // 9

        gcl_memcpy(test_out, mem_out, sizeof(cl_float) * NUM_VALUES);

    });

    struct timeval executing_time;
    gettimeofday(&executing_time, NULL);


    // Don't forget to free up the CL device's memory when you're done. // 10
    gcl_free(mem_in);
    gcl_free(mem_p);
    gcl_free(mem_mu);
    gcl_free(mem_out);

    // Finally, release your queue just as you would any GCD queue.    // 11
    dispatch_release(queue);

    struct timeval end_time;
    gettimeofday(&end_time, NULL);

    time_t devicediffsec = difftime(device_setting_time.tv_sec, start_time.tv_sec);
    suseconds_t devicediffsub = device_setting_time.tv_usec - start_time.tv_usec;
    double devicerealsec = devicediffsec + devicediffsub * 1e-6;
    printf("Device Setting: %f sec\n", devicerealsec);

    time_t arraydiffsec = difftime(array_setting_time.tv_sec, device_setting_time.tv_sec);
    suseconds_t arraydiffsub = array_setting_time.tv_usec - device_setting_time.tv_usec;
    double arrayrealsec = arraydiffsec + arraydiffsub * 1e-6;
    printf("Array Setting: %f sec\n", arrayrealsec);

    time_t executingdiffsec = difftime(executing_time.tv_sec, array_setting_time.tv_sec);
    suseconds_t executingdiffsub = executing_time.tv_usec - array_setting_time.tv_usec;
    double executingrealsec = executingdiffsec + executingdiffsub * 1e-6;
    printf("Executing: %f sec\n", executingrealsec);

    time_t terminationdiffsec = difftime(end_time.tv_sec, executing_time.tv_sec);
    suseconds_t terminationdiffsub = end_time.tv_usec - executing_time.tv_usec;
    double terminationrealsec = terminationdiffsec + terminationdiffsub * 1e-6;
    printf("Termination: %f sec\n", terminationrealsec);

    time_t totaldiffsec = difftime(end_time.tv_sec, start_time.tv_sec);
    suseconds_t totaldiffsub = end_time.tv_usec - start_time.tv_usec;
    double totalrealsec = totaldiffsec + totaldiffsub * 1e-6;
    printf("Total: %f sec\n", totalrealsec);


    // Check to see if the kernel did what it was supposed to:
    if ( validate(test_in, test_out)) {
        fprintf(stdout, "All values were OK.\n");
    }

    // And the same goes for system memory, as usual.
    free(test_in);
    free(test_out);
    free(test_p);
    free(test_mu);    
}

• kernel.cl で #define INLINE とするとインライン展開します。
• main.c で #define GPU とすると GPU を利用します。

実行結果

今回テストした環境は次の通りです。ちょっと古いGPUです。

Mac Pro (Mid 2010)
Processor 2.8GHz Quad-Core Intel Xeon
Memory 16GB
ATI Radeon HD 5770 1024MB

NVIDIA のグラフィックボードも試したのですが，OpenCLが対応していないようで，動作しませんでした。

実行結果はこんな感じです。

GPU (ループ)

Created a dispatch queue using the ATI Radeon HD 5770
Device Setting: 0.054858 sec
Array Setting: 0.534748 sec
Executing: 0.445068 sec
Termination: 0.037826 sec
Total: 1.072500 sec
All values were OK.


GPU (インライン展開)

Created a dispatch queue using the ATI Radeon HD 5770
Device Setting: 0.049323 sec
Array Setting: 0.542580 sec
Executing: 0.415137 sec
Termination: 0.040673 sec
Total: 1.047713 sec


だいたい1秒ジャストくらいでした！

配列の設定に0.5秒くらいの時間がかかっているのは相変わらずです。正味の実行時間は0.4秒ほどですから，すごいですね。

CPUと比べてGPUだとインライン展開したときに少しスピードアップしていますね。これはGPUは分岐予測などの機能が弱いからのように思います。

今までの実行時間をまとめると次の通りです。

Elixir(秒)	Elixir(秒)	Elixir(秒)	C言語(秒)	OpenCL(秒)	OpenCL(秒)	OpenCL(秒)	OpenCL(秒)
1並列ループ	8並列ループ	8並列インライン展開	1並列ループ	8並列ループ	8並列インライン展開	GPUループ	GPUインライン展開
52.795620	12.664873	11.308742	4.232451	1.496656	1.483530	1.072500	1.047713

• C言語の1並列ループ→OpenCLの8並列ループは2.82倍の速度向上でした
• C言語の1並列ループ→OpenCLのGPUループは3.95倍の速度向上でした。
• Elixirの1並列ループ→8並列ループは4.16倍の速度向上です。

GPUをもってしても並列処理による速度向上はElixirの方がありますね。OpenCLの場合は配列の設定に時間がかかっているので，その分のオーバーヘッドが大きいようです。

実際，正味の実行時間で計算してみると，次のようになります。CPU8並列でElixirと同等レベル，GPUだとさらに2倍以上の速度向上になります！

• C言語1並列ループ→OpenCL8並列ループで4.7倍の速度向上
• C言語1並列ループ→OpenCLGPUループで9.5倍の速度向上

マルチコアCPUでOpenCLを使わずに直にマルチコア並列プログラミングした場合の見積もり(正味の時間)と，OpenCLでGPUを用いた場合(配列設定時間を含む)を比較してみると，次のようになります。

• CPU(ループ) 0.896(秒) vs GPU(ループ)1.07(秒)
• CPU(インライン展開) 0.888(秒) vs GPU(インライン展開)1.05(秒)

このくらいの演算規模だとCPUでガリガリにチューニングしたほうが速そうです。CPUマルチコアかつSIMD命令とかAVX命令とかを使ったらGPU以上のパフォーマンスが出るかもしれませんね。もっとも，GPUに計算させている間にCPUで別の計算をするという使い方もできる可能性があるので，CPUに向いた計算とGPUに向いた計算でそれぞれ分業して負荷分散させるということをすると良さそうです。

• Elixirの8並列ループ→OpenCL CPUループ(合計)で8.47倍の速度向上です。
• Elixirの8並列ループ→OpenCL GPUループ(合計)で11.9倍の速度向上です。
• Elixirの8並列ループ→OpenCL CPUループ(正味)で14.2倍の速度向上です。
• Elixirの8並列インライン展開→OpenCL CPUインライン展開(合計)で7.64倍の速度向上です。
• Elixirの8並列インライン展開→OpenCL GPUインライン展開(合計)で10.8倍の速度向上です。
• Elixirの8並列インライン展開→OpenCL CPUインライン展開(正味)で12.7倍の速度向上です。

Elixirのコード生成をガリガリに最適化すると，このくらいのパフォーマンス向上が得られるという計算になります。正味と書いてある方が，CPUの場合でOpenCLを使わずに直にチューニングしたときにこのくらいの性能になるだろうという見込みです。

おわりに

1. OpenCLでGPUを利用した時にはC言語の1並列と比べて3.95倍の速度向上，Elixirの同等プログラムと比べて10.8〜11.9倍の速度向上になりました。
2. OpenCLを使わずにマルチコアかつSIMD命令やAVX命令を使った場合は，GPUの場合より高速になる可能性があります。見積もりではElixirの同等プログラムと比べて13倍前後の速度向上を期待できそうです。
3. GPUの場合は，データの転送に時間がかかっているので，データの転送量に比べて演算負荷が大きくなればなるほど，CPUよりGPUの方が有利になると思われます。
4. CPUとGPUで適性を見極めて適切に負荷分散をすること，さらにCPUとGPUを並列実行することで，さらなるパフォーマンスを引き出せる可能性があります。

ここまでで大体どのくらいの性能が引き出せるかが見積もれましたので，次のチャレンジとしては NIF(Native Implemented Functions)を使って，ElixirのプログラムからOpenCLのプログラムを呼び出すことを考えたいのですが，その前に小休止として，AI/MLのデファクトスタンダードであるPythonとElixir，OpenCLの比較をしてみましょう。

明日は @takasehideki さんの「ElixirでIoT#4：Nervesって何者？ラズパイでLチカできんの！？」です。こちらもお楽しみに！