現代のコンパイラの凄さを体感してしまったというお話

この記事について

行列演算を実行するプログラムで通常のfor文を用いて計算する方法とSIMD命令を用いて計算する場合で実行時間を測定し、SIMD命令を用いたプログラムが通常のプログラムを圧倒するというストーリーを仕立てようとGitHub Copilot君の力も借りつつ検証してみました。しかし、結果は。。。現代のコンパイラの凄さが分かったという結論に至ったというお話です。

検証内容について

SIMD(Single Instruction Multiple Data)命令を上手く用いると行列演算プログラムの高速化が実現出来ると言われます。これはSIMD命令は1つの演算命令を複数のデータに適用できるため、同じ演算を複数の要素に適用する様な演算、つまり、行列の内積計算の様な演算に対して威力を発揮します。よって、NEON命令を用いて実装した行列演算プログラムとfor文を用いて実装した行列演算プロウグラムをそれぞれ作成して実行時間を比較してみることにしました。なお、今回の検証ではソースの作成はGitHub Copilotを用いて行いました。このプログラムを生成する過程で何度か生成結果を実際にビルドし、実行時間を測定して、NEONを用いたソースの性能 > for文を用いたソースの性能を達成するまで生成を繰り返しました。

検証環境

• Raspberry Pi 5
  • CPU : Cortex-A76
  • DRAM : 16GB
• Ubuntu 24.04.4 LTS

実際に作成した検証プログラム

GitHub Copilotと何度かのやり取りを繰り返して最終的に完成した検証プログラムは以下です。MacとRasPiでの実行を想定して切り分ける分岐をMakefileに入れていますが、Macでの検証は実行していません。さて、プログラムを生成する過程でまず4x4程度の行列で性能差を出そうとすると要素が整数固定値の行列を与える問題ではどう頑張っても通常のfor文を用いたプログラムが性能で勝ってしまいます。よって、疑似乱数で実数行列を作りその演算を行うプログラムになっています。また、ウォームアップを追加しています。この理由をCopilot君に質問すると。。。

ウォームアップを入れる理由は、最初の数回だけ遅くなりやすい「初期状態の揺れ」を本計測から外すためです。
主な要因は次です。
キャッシュ未温状態
最初は命令キャッシュ・データキャッシュ・TLBが冷えていて、メモリアクセスが遅くなりやすいです。
CPU周波数の立ち上がり（DVFS）
Linux環境だと、処理開始直後は低クロックで走り、負荷が続くとクロックが上がることがあります。最初の計測だけ不利になります。
分岐予測などの安定化
ループ実行が進むと予測器の状態が学習され、挙動が安定します。
実行環境の初回オーバーヘッド
初回のページ確保やライブラリ関連の遅延要素が混じることがあります。
このプログラムの目的は Neon と通常実装の純粋な計算性能比較なので、ウォームアップで環境を整えてから測ると、差がより再現性よく出ます。
実務では「ウォームアップ後の複数回計測で中央値を採用」が定番です。

との事です。またMakefileではfor文を利用したプログラムのビルド時にはコンパイラが自動的にベクトル演算命令を用いて最適化を実行しないように止めるオプションを追加しています。これをしないとコンパイラでビルドすると勝手に最適化が実行された実行ファイルが生成されます。

/* POSIX.1-2008 の API を有効化する
clock_gettime / CLOCK_MONOTONIC を使うために必要
標準ヘッダより前で定義すること */
#define _POSIX_C_SOURCE 200809L
#include <stdio.h>
#include <stdint.h>
#include <time.h>

#ifdef USE_NEON
#include <arm_neon.h>
#endif

/* 4x4 固定の行列積ベンチマーク */
#define SIZE 4

/* 入力行列プール数（2のべき乗にして idx 計算を軽くする） */
#define POOL_SIZE 256

/* CPUやキャッシュを温めるための試行回数 */
#define WARMUP_ITERS 20000UL

/* 1 trial あたりの計測回数 */
#define BENCH_ITERS 10000000UL

/* ノイズ低減のための試行回数（中央値を採用） */
#define TRIALS 7

/* 関数のインライン展開を抑制して、計測の偏りを減らす */
#if defined(__GNUC__) || defined(__clang__)
#define NOINLINE __attribute__((noinline))
#else
#define NOINLINE
#endif

/* 毎回同じ入力にならないよう、複数の入力行列を保持する */
static int32_t a_pool[POOL_SIZE][SIZE][SIZE];
static int32_t b_pool[POOL_SIZE][SIZE][SIZE];

/* timespec（sec+nsec）を単一のナノ秒値へ変換 */
static inline uint64_t timespec_to_ns(const struct timespec *ts)
{
    return (uint64_t)ts->tv_sec * 1000000000ULL + (uint64_t)ts->tv_nsec;
}

/* 値（ns）を見やすい単位へ自動変換する */
static const char *pick_time_unit(double *value_ns)
{
    if (*value_ns >= 1e9)
    {
        *value_ns /= 1e9;
        return "s";
    }
    if (*value_ns >= 1e6)
    {
        *value_ns /= 1e6;
        return "ms";
    }
    if (*value_ns >= 1e3)
    {
        *value_ns /= 1e3;
        return "us";
    }
    return "ns";
}

/* 結果確認用の行列表示 */
static void print_matrix(const char *name, const int32_t m[SIZE][SIZE])
{
    printf("%s:\n", name);
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            printf("%6d ", m[i][j]);
        }
        printf("\n");
    }
}

/* 再現性のある疑似乱数で入力プールを初期化 */
static void init_input_pool(void)
{
    uint32_t seed = 0x12345678u;

    for (int p = 0; p < POOL_SIZE; p++)
    {
        for (int i = 0; i < SIZE; i++)
        {
            for (int j = 0; j < SIZE; j++)
            {
                seed = seed * 1664525u + 1013904223u;
                a_pool[p][i][j] = (int32_t)((seed >> 27) - 16); /* -16..15 */
                seed = seed * 1664525u + 1013904223u;
                b_pool[p][i][j] = (int32_t)((seed >> 27) - 16); /* -16..15 */
            }
        }
    }
}

#ifndef USE_NEON
/* 非SIMD（スカラ）実装 */
static NOINLINE void matrix_multiply_standard(
    const int32_t a[SIZE][SIZE],
    const int32_t b[SIZE][SIZE],
    int32_t result[SIZE][SIZE])
{
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            int32_t sum = 0;
            for (int k = 0; k < SIZE; k++)
            {
                sum += a[i][k] * b[k][j];
            }
            result[i][j] = sum;
        }
    }
}
#endif

#ifdef USE_NEON
/* 手書きNEON実装（4列を同時処理） */
static NOINLINE void matrix_multiply_neon(
    const int32_t a[SIZE][SIZE],
    const int32_t b[SIZE][SIZE],
    int32_t result[SIZE][SIZE])
{
    /* Bの各行をベクトルとしてロード */
    int32x4_t b0 = vld1q_s32(b[0]);
    int32x4_t b1 = vld1q_s32(b[1]);
    int32x4_t b2 = vld1q_s32(b[2]);
    int32x4_t b3 = vld1q_s32(b[3]);

    for (int i = 0; i < SIZE; i++)
    {
        /* Aの1要素を4レーンへ複製して積和 */
        int32x4_t a0 = vdupq_n_s32(a[i][0]);
        int32x4_t a1 = vdupq_n_s32(a[i][1]);
        int32x4_t a2 = vdupq_n_s32(a[i][2]);
        int32x4_t a3 = vdupq_n_s32(a[i][3]);

        int32x4_t c = vmulq_s32(a0, b0);
        c = vmlaq_s32(c, a1, b1);
        c = vmlaq_s32(c, a2, b2);
        c = vmlaq_s32(c, a3, b3);

        vst1q_s32(result[i], c);
    }
}
#endif

/* 1 trial 分の計測を実行し、経過ナノ秒を返す */
static uint64_t run_once(volatile int32_t *checksum, int32_t last_result[SIZE][SIZE])
{
    struct timespec t0, t1;

    clock_gettime(CLOCK_MONOTONIC, &t0);

    for (uint64_t i = 0; i < BENCH_ITERS; i++)
    {
        /* 入力プールを巡回して、同一入力固定の最適化を避ける */
        uint64_t idx = i & (POOL_SIZE - 1);
#ifdef USE_NEON
        matrix_multiply_neon(a_pool[idx], b_pool[idx], last_result);
#else
        matrix_multiply_standard(a_pool[idx], b_pool[idx], last_result);
#endif
        /* 計算結果を消されないよう副作用を作る */
        *checksum += last_result[i & 3][(i >> 2) & 3];
    }

    clock_gettime(CLOCK_MONOTONIC, &t1);
    return timespec_to_ns(&t1) - timespec_to_ns(&t0);
}

/* 小さい配列なので挿入ソートで中央値を取得 */
static uint64_t median_ns(uint64_t a[TRIALS])
{
    for (int i = 1; i < TRIALS; i++)
    {
        uint64_t key = a[i];
        int j = i - 1;
        while (j >= 0 && a[j] > key)
        {
            a[j + 1] = a[j];
            j--;
        }
        a[j + 1] = key;
    }
    return a[TRIALS / 2];
}

int main(void)
{
    volatile int32_t checksum = 0;
    int32_t result[SIZE][SIZE] = {0};
    uint64_t trials[TRIALS];

    init_input_pool();

#ifdef USE_NEON
    printf("Mode: NEON\n");
#else
    printf("Mode: STANDARD\n");
#endif

    /* ウォームアップ（計測には含めない） */
    for (uint64_t i = 0; i < WARMUP_ITERS; i++)
    {
        uint64_t idx = i & (POOL_SIZE - 1);
#ifdef USE_NEON
        matrix_multiply_neon(a_pool[idx], b_pool[idx], result);
#else
        matrix_multiply_standard(a_pool[idx], b_pool[idx], result);
#endif
        checksum += result[i & 3][(i >> 2) & 3];
    }

    /* 複数 trial を実行して頑健な代表値を取る */
    for (int t = 0; t < TRIALS; t++)
    {
        trials[t] = run_once(&checksum, result);
    }

    uint64_t med_ns = median_ns(trials);
    double per_op_ns = (double)med_ns / (double)BENCH_ITERS;
    double ops_per_sec = 1e9 / per_op_ns;

    double total = (double)med_ns;
    const char *total_unit = pick_time_unit(&total);

    double per_op = per_op_ns;
    const char *per_op_unit = pick_time_unit(&per_op);

    printf("Median elapsed: %.3f %s (iters=%llu, trials=%d)\n",
           total, total_unit, (unsigned long long)BENCH_ITERS, TRIALS);
    printf("Time per multiply: %.3f %s/op\n", per_op, per_op_unit);
    printf("Throughput: %.2f Mops/s\n", ops_per_sec / 1e6);

    print_matrix("Result (last)", result);

    (void)checksum;
    return 0;
}

UNAME_S := $(shell uname -s)

SRC     := main.c
OUT_STD := matrix_standard.out
OUT_NEON:= matrix_neon.out

CFLAGS  := -O3 -Wall

ifeq ($(UNAME_S),Darwin)
CC           := clang
SCALAR_FLAGS := -fno-vectorize
NEON_FLAGS   := -DUSE_NEON
else
CC           := gcc
# 非SIMD比較用: 自動ベクトル化のみ抑制
SCALAR_FLAGS := -fno-tree-vectorize -fno-tree-slp-vectorize
# 手書きNEON用
NEON_FLAGS   := -march=armv8-a+simd -DUSE_NEON
endif

all: standard neon

standard: $(SRC)
	$(CC) $(CFLAGS) $(SCALAR_FLAGS) -o $(OUT_STD) $(SRC)

neon: $(SRC)
	$(CC) $(CFLAGS) $(NEON_FLAGS) -o $(OUT_NEON) $(SRC)

clean:
	rm -f *.out

実行時間を計測しました。確かにNEONを用いたプログラムの方が性能で勝っています。

$ ./matrix_neon.out                                                                    
Mode: NEON                                                                             
Median elapsed: 133.803 ms (iters=10000000, trials=7)                                   
Time per multiply: 13.380 ns/op                                                         
Throughput: 74.74 Mops/s                                                                
Result (last):                                                                          
    73    168   -240    176                                                             
  -243   -197     56    -16                                                             
   234    309    -48     86                                                             
   -96   -159    -36    -50 

$ ./matrix_standard.out
Mode: STANDARD                                                                      
Median elapsed: 278.375 ms (iters=10000000, trials=7)                               
Time per multiply: 27.838 ns/op                                                     
Throughput: 35.92 Mops/s                                                            
Result (last):                                                                      
    73    168   -240    176                                                         
  -243   -197     56    -16                                                         
   234    309    -48     86                                                         
   -96   -159    -36    -50

まとめ

最初は検算が出来る方が良いだろうと簡単な行列のプログラムを書いて検証出来ないかと検証してみたものの全く歯が立たず、コンパイラに色々とオプションを渡して、プログラムも複雑怪奇にしてやっと差が出せました。このソースを完成させるのに複数回試作を繰り返しました。これは問題設定に問題が有って、当初行列サイズが小さくなおかつ固定値の整数行列の内積だったので差が出にくい問題設定だったこと、そして、現代のコンパイラはこれまでの研究で実装者が何もしなくても様々な最適化を自動的に適用してくれるため、特別な命令を用いなくても(特殊な問題を除いて)高速に動いてしまうことが要因です。筆者はこれまであまり低レイヤーに関わってこなかったので、この辺りの知識が乏しいため、今回手元の環境で実験してみましたが、NEON命令を手書きで書く必要が有る場面は本当に限定的な特殊な場面で通常はコンパイラが頑張ってくれるということが分かった検証となりました。本当は失敗作のソースの生成過程も掲載したかったのですが、Copilot君とのやり取りを取り分けて引用するのが想像以上に煩雑だったため、今回は断念しました。但し、最初に作った初号機のソースだけ以下に掲載しておこうと思います。

#include <stdio.h>
#include <stdint.h>

#ifdef USE_NEON
#include <arm_neon.h>
#endif

#define SIZE 4

#ifndef USE_NEON
static void matrix_multiply_standard(
    const int32_t a[SIZE][SIZE],
    const int32_t b[SIZE][SIZE],
    int32_t result[SIZE][SIZE])
{
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            result[i][j] = 0;
            for (int k = 0; k < SIZE; k++)
            {
                result[i][j] += a[i][k] * b[k][j];
            }
        }
    }
}
#endif

#ifdef USE_NEON
static void matrix_multiply_neon(
    const int32_t a[SIZE][SIZE],
    const int32_t b[SIZE][SIZE],
    int32_t result[SIZE][SIZE])
{
    int32x4_t b0 = vld1q_s32(b[0]);
    int32x4_t b1 = vld1q_s32(b[1]);
    int32x4_t b2 = vld1q_s32(b[2]);
    int32x4_t b3 = vld1q_s32(b[3]);

    for (int i = 0; i < SIZE; i++)
    {
        int32x4_t a0 = vdupq_n_s32(a[i][0]);
        int32x4_t a1 = vdupq_n_s32(a[i][1]);
        int32x4_t a2 = vdupq_n_s32(a[i][2]);
        int32x4_t a3 = vdupq_n_s32(a[i][3]);

        int32x4_t c = vmulq_s32(a0, b0);
        c = vmlaq_s32(c, a1, b1);
        c = vmlaq_s32(c, a2, b2);
        c = vmlaq_s32(c, a3, b3);

        vst1q_s32(result[i], c);
    }
}
#endif

static void print_matrix(const char *name, const int32_t m[SIZE][SIZE])
{
    printf("%s:\n", name);
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            printf("%6d ", m[i][j]);
        }
        printf("\n");
    }
}

int main(void)
{
    int32_t a[SIZE][SIZE] = {
        {1, 2, 3, 4},
        {5, 6, 7, 8},
        {9, 10, 11, 12},
        {13, 14, 15, 16}};

    int32_t b[SIZE][SIZE] = {
        {17, 18, 19, 20},
        {21, 22, 23, 24},
        {25, 26, 27, 28},
        {29, 30, 31, 32}};

    int32_t result[SIZE][SIZE] = {0};

#ifdef USE_NEON
    matrix_multiply_neon(a, b, result);
    printf("Mode: NEON\n");
#else
    matrix_multiply_standard(a, b, result);
    printf("Mode: STANDARD\n");
#endif

    print_matrix("A", a);
    print_matrix("B", b);
    print_matrix("Result (A x B)", result);

    return 0;
}

UNAME_S := $(shell uname -s)

CFLAGS  := -O2 -Wall
SRC     := main.c

ifeq ($(UNAME_S),Darwin)
CC      := clang
NEON_ARCH :=
else
CC      := gcc
NEON_ARCH := -march=armv8-a+simd
endif

all: standard neon

standard: $(SRC)
	$(CC) $(CFLAGS) -o matrix_standard.out $(SRC)

neon: $(SRC)
	$(CC) $(CFLAGS) $(NEON_ARCH) -DUSE_NEON -o matrix_neon.out $(SRC)

clean:
	rm -f *.out