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

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

昨日は @twinbee さんの「Elixir並列処理「Flow」の2段ステージ構造を理解する」でしたね。

「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の同時セッション最大数・レイテンシ改善の構想」です。Qiitaで連載中のGPUの話とはまた別の，ZEAMにつながる新しいメカニズムの構想について熱く語ります！

さて，本題〜はじめに

前回のZEAM開発ログv0.1.1 AI/MLを爆速にしたい！ Flow / GenStage でGPUを駆動できないの？では，次のようにまとめました。

1. **GPUはSIMD(シムディー)**というモデルで動作します。これに対し，**CPUはMIMD(ミムディー)**というモデルで動作します。
2. GPUはSIMDなので，単純な構造で均質で大量にあるデータを，同じような命令列で処理する場合に効果を発揮します。
3. ElixirのFlowを手がかりに最適化すると，GPUに向いたプログラムに変換することができます。
4. C言語はElixirよりざっと12倍は速い。Elixirにはまだまだ高速化の余地があります。

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

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

今回は，いよいよGPUを駆動するプログラムを書いてみたいと思いますが，その前段階として，OpenCLを使ってCPUのポテンシャルを使い切ってみたいと思います。

OpenCLとは？

OpenCLは，アップルによって提案され，現在は標準化団体クロノス・グループの OpenCL Working Group によって策定されている標準規格で，マルチコアCPUやGPU，DSP，FPGAなど，異なる計算資源を活用するコンピュータシステム(ヘテロジーニアス環境とか，ヘテロジーニアス・コンピューティングとか呼びます)でのクロスプラットフォームな並列プログラミングフレームワークです。

• ヘテロジーニアス環境ということは，ほぼ同一のプログラムで，マルチコアCPUもGPUもDSPもFPGAも駆動できます。
• クロスプラットフォームということは，ほぼ同一のプログラムで，MacでもLinuxでもWindowsでも動作します。

Mac の OpenCL の対応状況はこんな感じです。

OpenGL および OpenCL グラフィックスを扱う Mac コンピュータ

Xcode で OpenCL のプログラミングを簡単に始められます。アップル公式のプログラミングガイドはこちらです。今回は，このプログラミングガイドを参考に，ほぼ丸写ししています。

OpenCL Programming Guide for Mac

OpenCLによるロジスティック写像のベンチマークプログラム

では，OpenCLでロジスティック写像のベンチマークプログラムを書いてみましょう。(今回はあえてまだGitHubレポジトリをご紹介しません！)

kernel.cl (ループ)

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];
    }
}

kernel.cl (インライン展開)

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;
}

main.c

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

#include <OpenCL/opencl.h>
#include "kernel.cl.h"

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

#define GPU

#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;
}

// 計算結果が合っているかの確認
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);
    
    // CPUデバイスを開く
    dispatch_queue_t queue;
    queue = gcl_create_dispatch_queue(CL_DEVICE_TYPE_CPU, NULL);
    cl_device_id cpu = gcl_get_device_id_with_dispatch_queue(queue);
    clGetDeviceInfo(cpu, 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);

    // 値の入力 in 用のホストのメモリを確保して初期化する
    int *test_in = (int *)malloc(sizeof(cl_int) * NUM_VALUES);
    for (i = 0; i < NUM_VALUES; i++) {
        test_in[i] = (cl_int)i;
    }
    
    // 結果の出力 out 用のホストのメモリを確保して初期化する
    int *test_out = (int *)malloc(sizeof(cl_int) * NUM_VALUES);

    // 入力 MU 用のホストのメモリを確保して初期化する    
    int *test_mu = (int *)malloc(sizeof(cl_int) * NUM_VALUES);
    for (i = 0; i < NUM_VALUES; i++) {
        test_mu[i] = (cl_int)MU;
    }

    // 入力 P 用のホストのメモリを確保して初期化する    
    int *test_p = (int *)malloc(sizeof(cl_int) * NUM_VALUES);
    for (i = 0; i < NUM_VALUES; i++) {
        test_p[i] = (cl_int)P;
    }

    // 入力について，デバイスのメモリを確保してホストからコピーする
    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);

    // 出力について，デバイス側のメモリを確保する
    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);
    
    // kernel.cl で書かれたプログラムを呼び出す
    dispatch_sync(queue, ^{
        size_t wgs;
        gcl_get_kernel_block_workgroup_info(logisticsmap_kernel,
                                            CL_KERNEL_WORK_GROUP_SIZE,
                                            sizeof(wgs), &wgs, NULL);
        
        // kernelで実行するN次元の範囲
        // この場合は，1次元
        cl_ndrange range = {
            1,                     // 次元の数
            
            {0, 0, 0},             // それぞれの次元のオフセット値。この場合はすべてのデータを処理したいので0にセットする
            
            {NUM_VALUES, 0, 0},    // 各次元にいくつの値が含まれるか
            
            {wgs, 0, 0}            // ワークグループごとのローカルサイズ
                                   // これによりワークグループごとのワーク項目の数が決まる
                                   // これは間接的にワークグループの数に影響を及ぼす
        };
        // カーネルの呼び出しは単純である。単に関数のように呼び出す。
        // 第1引数として ndrange，第2引数以下でそれぞれの引数を渡す。
        // 整数型の場合は cl_int を使う        
        logisticsmap_kernel(&range,(cl_int*)mem_mu, (cl_int*)mem_p, (cl_int*)mem_in, (cl_int*)mem_out);
        
        // 出力された結果を，デバイスからホストにコピーする
        
        gcl_memcpy(test_out, mem_out, sizeof(cl_float) * NUM_VALUES);
        
    });

    // 時間計測(実行時間)
    struct timeval executing_time;
    gettimeofday(&executing_time, NULL);

    
    // OpenCL のデバイスのメモリの解放を忘れずに
    gcl_free(mem_in);
    gcl_free(mem_p);
    gcl_free(mem_mu);
    gcl_free(mem_out);

    // キューを解放する
    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);

    
    // 値の整合性チェック
    if ( validate(test_in, test_out)) {
        fprintf(stdout, "All values were OK.\n");
    }
    
    // システムメモリの解放
    free(test_in);
    free(test_out);
    free(test_p);
    free(test_mu);
}

実行結果

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

ループ

CPU (non-INLINE)
Created a dispatch queue using the Intel(R) Xeon(R) CPU           W3530  @ 2.80GHz
Device Setting: 0.047694 sec
Array Setting: 0.522358 sec
Executing: 0.896386 sec
Termination: 0.030218 sec
Total: 1.496656 sec
All values were OK.

インライン展開

CPU (INLINE)
Created a dispatch queue using the Intel(R) Xeon(R) CPU           W3530  @ 2.80GHz
Device Setting: 0.048654 sec
Array Setting: 0.518187 sec
Executing: 0.887887 sec
Termination: 0.028802 sec
Total: 1.483530 sec
All values were OK.

配列の設定に意外と時間がかかっていますね。

またインライン展開してもそれほど速度差はないみたいです。これはCPUは分岐予測などの高度な機能が備わっているからかもしれません。分岐予測などが整っていないGPUの場合にどうなるかが興味が湧くところですね！

今までの実行時間をまとめるとこんな感じです。

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

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

並列処理による速度向上はElixirの方がありますね。OpenCLの場合は配列の設定に時間がかかっているので，その分のオーバーヘッドだと考えられます。実際，正味の実行時間で計算してみると，C言語1並列ループ→OpenCL8並列ループで4.7倍の速度向上なので，Elixirと同等レベルなので，当たっていそうです。おそらくマルチコアCPUをC言語化するときはOpenCLを使わずに直にマルチコア並列プログラミングした方が速いと思われます。

• Elixirの8並列ループ→OpenCL8並列ループで8.47倍の速度向上です。
• Elixirの8並列インライン展開→OpenCL8並列インライン展開で7.64倍の速度向上です。

Elixirのコード生成をガリガリに最適化すると，このくらいのパフォーマンス向上が得られるという計算になります。

おわりに

1. OpenGLはマルチコアCPUやGPUをはじめ，いろいろな環境で並列プログラミングをすることができます。
2. OpenCLでは配列の設定に意外と時間がかかるので，Elixirで並列数を増やしたときの速度向上より，C言語の1並列からOpenCLのCPU8並列の速度向上の方が小さかったです。おそらくマルチコアCPUをC言語化するときはOpenCLを使わずに直にマルチコア並列プログラミングした方が速いと思われます。
3. OpenCLでCPUを利用した時にはElixirの同等プログラムと比べて7.64〜8.47倍の速度向上になりました。直にマルチコア並列プログラミングしたときの性能が楽しみです！

今回でGPUを駆動したデータを紹介するつもりでしたが，思ったよりコード量があったので，マルチコアCPUだけにとどまりました。次回「ZEAM開発ログv0.1.3 AI/MLを爆速にしたい！ Flow のコードを OpenCL で書いてみる〜GPU編」
です！ お楽しみに！

明日は， @takasehideki さんの「ElixirでIoT#3：IoTボードで動いた！Phoenixが立った！性能評価と考察」です。こちらもお楽しみに！