ESP32のC/C++の関数をアセンブリ言語化し手作業で最適化してみる
はじめに
前回の記事ではESP32を対象に、アセンブリ言語の読み方を説明しました。今回はもう少し具体的な例を挙げて説明をしていきます。
前提条件
-
XtensaコアのESP32シリーズ (ESP32/ESP32-S2/ESP32-S3) を使用していること。 -
ArduinoIDEまたはVSCode+PlatformIOでESP32用のプログラムを実行できること。 - 前回までの記事の内容までを把握していること。
※ ESP32-C3などはXtensaコアではないため、本記事の対象から外れます。
今回取り扱うサンプルコード
今回はC/C++で作成した関数をアセンブリ言語化し手作業で最適化する手法を説明します。
比較的手軽にアセンブリ言語化の恩恵が得やすい 『配列に対して何らかの演算を行い、結果を配列に格納する』 という処理例を用意します。
今回取り扱う例は以下のようにしました。
void test_func1(float* dst, const float* src, uint32_t len, float vol)
{
for (uint32_t i = 0; i < len; i++) {
dst[i] += src[i] * vol;
}
}
これは、float型の配列に入った波形データにボリューム係数を掛け、別の波形と加算して複数の波形を合成する場面…を想定しています。
引数は以下のようになります。
-
dst: 合成後の波形のfloat配列 -
src: 合成元の波形のfloat配列 -
len: 配列の要素数 -
vol:src波形に掛けるボリューム係数値
合成元の波形データとなる配列を4個用意し、それぞれに対してこの関数を呼び出すことで、dstに合成後の波形を作成します。
検証のためのコード全体像
最適化を始める前に、まずはC/C++で作成した関数の処理時間を確認しておきます。
ビルドするコードの全体像は以下の通りです。
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <freertos/FreeRTOS.h>
#include <freertos/task.h>
// アセンブリ言語で記述した関数を呼び出すための宣言
extern "C" {
void test_func1(float* dst, const float* src, uint32_t len, float level);
void test_func2(float* dst, const float* src, uint32_t len, float level);
};
// 波形データ用の配列要素数
static constexpr const uint32_t array_len = 1024;
static constexpr const uint32_t src_count = 4;
// 波形データ用の配列要素数 (SIMD命令を試すために16バイトアライメントしておく)
float src_data[src_count][array_len] __attribute__ ((aligned (16)));
float dst_data[array_len] __attribute__ ((aligned (16)));
// 最適化したい処理1
// weak指定しておき、同名の関数をアセンブリ言語で作成することで処理を上書き可能にしておく
__attribute((noinline, noclone, weak, optimize("-O3")))
void test_func1(float* dst, const float* src, uint32_t len, float vol) {
for (uint32_t i = 0; i < len; i++) {
dst[i] += src[i] * vol;
}
}
// 最適化したい処理2
// weak指定しておき、同名の関数をアセンブリ言語で作成することで処理を上書き可能にしておく
__attribute((noinline, noclone, weak, optimize("-O3")))
void test_func2(float* dst, const float* src, uint32_t len, float vol) {
for (uint32_t i = 0; i < len; i+=4) {
dst[i+0] += src[i+0] * vol;
dst[i+1] += src[i+1] * vol;
dst[i+2] += src[i+2] * vol;
dst[i+3] += src[i+3] * vol;
}
}
// 性能評価用の関数
uint32_t bench(void) {
// 合成後のデータをクリアしておく
memset(dst_data, 0, array_len * sizeof(float));
// 元データの用意。実用する際は波形データを用意するが、今回は乱数で埋める
for (int j = 0; j < src_count; ++j) {
for (int i = 0; i < array_len; ++i) {
src_data[j][i] = rand();
}
}
// 乱数を用いてボリューム係数を設定
float vol = rand() / 65536.0f + 1;
// 測定前にvTaskDelayを呼んでおくことで測定中に割り込みが入る可能性を低減
vTaskDelay(1);
uint32_t cycle = 0;
// 計測対象の処理を合計src_count回実行
for (int i = 0; i < src_count; ++i) {
uint32_t c0, c1;
// 処理前のCPUサイクル値を取得
__asm__ __volatile("rsr %0, ccount" : "=r"(c0) );
// ここで計測対象となる最適化したい処理を呼び出す
test_func1(dst_data, src_data[i], array_len, vol);
// test_func2(dst_data, src_data[i], array_len, vol);
// 処理後のCPUサイクル値を取得
__asm__ __volatile("rsr %0, ccount" : "=r"(c1) );
// 処理に要したCPUサイクル数を加算。20サイクル程度オーバーヘッドを差し引く
cycle += c1 - c0 - 20;
}
// 処理が正しく行われたか確認
for (int i = 0; i < array_len; ++i) {
float result = 0;
for (int j = 0; j < src_count; ++j) {
result += src_data[j][i] * vol;
}
// 計算結果が合わない要素があれば0を返す
// 最適化の状況によって浮動小数演算の誤差に差が生じるため、誤差を広めに許容する
if (fabsf(dst_data[i] - result) > (FLT_EPSILON * 2.0f) * fmaxf(1.f, fmaxf(fabsf(dst_data[i]), fabsf(result)))) { return 0; }
}
// 処理に要したCPUサイクルを返す
return cycle;
}
void loop(void) {
vTaskDelay(1000);
// 処理に要したサイクルカウントをシリアル出力
uint32_t result = bench();
uint32_t total_len = array_len * src_count;
uint32_t cycle = (result * 100 + (total_len>>1)) / total_len;
printf("total : %u / %u = %4.2f cycle\n", result, array_len * src_count, cycle / 100.0f);
}
void setup(void) {}
これをESP32で実行すると、処理に要したCPUサイクル数がシリアルモニタ表示されます。totalが処理に要したCPUサイクル数の合計で、これを配列の要素数 1024 × 4回 = 4096で割って、ループ一回あたりの所要CPUサイクル数を算出しています。
C/C++の最適化オプションの違いによる所要CPUサイクル数の確認
このソースコードのtest_func1の前行にある optimize("-O3") を書き換えて最適化レベルを変更し、それぞれの処理結果を求めた結果が以下のようになります。
| optimize | total | cycle |
|---|---|---|
| -Os | 61444 | 15 |
| -O1 | 49168 | 12 |
| -O2 | 36880 | 9 |
| -O3 | 36880 | 9 |
アセンブリ言語に変換する
test_func1を前回お伝えした Compiler Explorer を使って、test_func1をアセンブリ言語にしたソースコードを入手しましょう。optimize("-O3")での結果は以下の通りです。
test_func1(float*, float const*, unsigned long, float):
entry sp, 32 // 関数の先頭に必須のレジスタの退避処理
wfr f2, a5 // 浮動小数レジスタ f2 = vol の内容を書き込む
beqz.n a4, .L1 // ループ回数が0なら.L1(終了処理)へジャンプ
slli a8, a4, 2 // a8 = len a4 << 2;
addi a8, a8, -4 // a8 -= 4;
srli a8, a8, 2 // a8 >>= 2;
addi.n a8, a8, 1 // a8++;
loop a8, .L3_LEND // a8 の回数、.L3_LEND までの範囲をループ開始
lsi f1, a3, 0 // f1 = src[0];
lsi f0, a2, 0 // f0 = dst[0];
addi.n a3, a3, 4 // srcアドレスを4加算
madd.s f0, f1, f2 // f0 += f1 * f2
ssi f0, a2, 0 // dst[0] = f0;
addi.n a2, a2, 4 // dstアドレスを4加算
.L1: // このラベルの箇所に .L3_LEND も存在しているはず
retw.n // 関数終了
コメントを追加しましたので個別の説明は不要かと思いますが、前回までの記事に含まれていなかった命令が幾つかあるので少し解説します。
新しく登場した命令の説明
wfr 命令
wfr命令はアドレスレジスタから浮動小数レジスタへ値を移す命令です。
恐らく Write Float Register の略と思われますが、公式資料には Move AR to FR と記載されており、Wが何の略か不明でした。
この例では wfr f2, a5 となっています。第4引数 a5 の float vol の値を浮動小数演算用の f2 レジスタに移しています。
関数の引数はアドレスレジスタで受け渡しますが、float型も例外ではありません。しかし、アドレスレジスタは整数処理用のレジスタのため、浮動小数値として扱うことはできません。
浮動小数計算を扱う f0～f15 という専用のレジスタが用意されているので、そちらに移し替える必要があるのです。
beqz 命令
これは分岐命令 branch equal zero です。チェック対象のアドレスレジスタの値がゼロの場合は指定したラベルにジャンプします。
今回の例では beqz a4, .L1 となっています。a4 つまり第3引数 len の値がゼロなら .L1 の位置へジャンプしています。これにより、配列の長さがゼロなら何もせず終了します。
slli 命令 と srli 命令
左ビットシフト命令 Shift Left Logical Immediate と
右ビットシフト命令 Shift Right Logical Immediate です。
今回の例ではslli a8, a4, 2 は a8 = a4 << 2; 、srli a8, a8, 2 は a8 = a8 >> 2; の動作です。
loop 命令
これはループ命令です。Xtensaの特に便利な機能なので、これは後ほど詳しく解説します。ひとまず、指定されたレジスタ値の回数だけ、指定されたラベル位置までの処理を繰り返し実行するもの、とだけ説明しておきます。今回の例では loop a8, .L3_LEND ですから、a8の数だけ.L3_LENDまでの処理を繰り返し実行する動作になります。
lsi 命令
これはfloat値のロード命令 Load Single Immediate です。
今回の例ではlsi f1, a3, 0ですから、a3 つまり引数src[0] の位置からf1レジスタに読み出しています。
madd.s 命令
これはfloat値の掛け算と足し算を一度に行う命令 Multiply and Add Single です。
今回の例では madd.s f0, f1, f2ですから、f0 += f1 * f2 という処理を一度の命令で実行します。
ssi 命令
これはfloat値のストア命令 Store Single Immediate です。
今回の例では ssi f0, a2, 0ですから、a2 つまり引数 dst[0] の位置に f0 の値を書き込んでいます。
アセンブリ言語版の関数を使用する
次にこれを実機上で実行してみます。
前々回の記事では、C/C++のコード内にインライン・アセンブラで埋め込む記述方法を紹介しましたが、インライン・アセンブラでは長い処理を書くのに記述が煩雑なことや、デバッガで処理を追うことができないなどのデメリットがあります。そこで今回は独立したアセンブリ言語のソースファイルを用いる方法を紹介します。
以下のファイルをC/C++のソースと同じフォルダに配置してください。
ファイルの拡張子は大文字で .S です。
#if defined ( __XTENSA__ )
// void test_func1(float* dst, const float* src, uint32_t len, float vol);
// a2 = float* dst 出力先データ
// a3 = float* src 加算元データ
// a4 = uint32_t len ループ回数
// a5 = float vol 倍率
.global test_func1 // グローバル関数 test_func1 があることを示す
.section .text // .textセクションに配置することを示す
.align 4 // 4バイトアライメントを指定
test_func1: // 関数 test_func1 の先頭アドレスを示すラベル
entry sp, 32 // 関数の先頭に必須のレジスタの退避処理
beqz a4, FUNC1_LOOP_END // ループ回数が0ならFUNC1_LOOP_END(終了処理)へジャンプ
wfr f2, a5 // 浮動小数レジスタ f2 = vol の内容を書き込む
loop a4, FUNC1_LOOP_END // a4 の回数、FUNC1_LOOP_END までの範囲をループ開始
lsi f1, a3, 0 // f1 = src[0];
lsi f0, a2, 0 // f0 = dst[0];
addi a3, a3, 4 // srcアドレスを4加算
madd.s f0, f1, f2 // f0 += f1 * f2
ssi f0, a2, 0 // dst[0] = f0;
addi a2, a2, 4 // dstアドレスを4加算
FUNC1_LOOP_END:
retw
#endif
このファイルを先ほどのファイルと一緒にビルドすると、test_func1関数はC/C++版ではなくこのアセンブリ言語版のコードが実行されるようになります。
C/C++版のtest_func1関数にはweak属性を付与しているので、同名の関数が別に存在する場合はweakを指定していない方が使用されます。
アセンブリ言語版の動作状況を確認する
現時点でのアセンブリ言語版の実行結果はC/C++版の-O2,-O3と差がないので、本当にアセンブリ言語版が使用されているか不安になると思います。C/C++版の関数を削除したり、C/C++版の最適化オプションを-O1に変更するなどして実行結果を確認し、アセンブリ言語版が使用されていることを確かめてみてください。
ループ内で実行される命令数と実際に処理に要したサイクル数を比較する
さて先ほど確認したループ一回当たりのCPUサイクル数は-O3では9でした。
Xtensaの命令セットは基本的に1命令あたり1サイクルで実行されます。実際の命令数が9個あるか確認してみましょう。
// ～～ 省略 ～～
loop a4, FUNC1_LOOP_END // a4 の回数、FUNC1_LOOP_END までの範囲をループ開始
lsi f1, a3, 0 // f1 = src[0];
lsi f0, a2, 0 // f0 = dst[0];
addi a3, a3, 4 // srcアドレスを4加算
madd.s f0, f1, f2 // f0 += f1 * f2
ssi f0, a2, 0 // dst[0] = f0;
addi a2, a2, 4 // dstアドレスを4加算
FUNC1_LOOP_END:
// ～～ 省略 ～～
ループ処理に関わる部分は上記の範囲のみです。行数は9…ありませんね。命令数としては7個、ラベルの行を加えても8行です。そうするとループ処理のために追加のCPUサイクルが発生しているのでしょうか？
実際の命令数と所要サイクル数の関係について詳しく見て行きましょう。
loop命令 : Xtensaのゼロオーバーヘッドループ機能
ループに要する実行サイクル数を確認するために、ここでloop命令について詳しく説明します。
アセンブリ言語でループ処理を構築する場合、一般的には『ループカウンタの加減算処理』と『条件ジャンプ命令』という2個の処理を組み合わせます。レジスタを用いてループ回数をカウントし、ループの末尾にて条件分岐命令により回数に達していなければループ先頭へジャンプする、という流れです。
Xtensaでも同様の処理を行うことも出来ますが、単純なループを効率よく行える、loop命令が用意されています。loop命令では、最初にループの回数と終端アドレスを設定しておくと、以後は処理が終端に到達するたびに回数の判定とループ先頭へのジャンプが自動的に動作します。
今回の例では、FUNC1_LOOP_END:の前行 addi a2, a2, 4が実行された次のCPUサイクルでは、loop命令の次行にあるlsi f1, a3, 0が実行されます。
これはゼロオーバーヘッドループと呼ばれており、『ループカウンタのインクリメント処理』も『条件ジャンプ命令』も不要で、純粋にループ内にある命令のみを繰り返し実行できます。
さて、loop命令がゼロオーバーヘッドで実行される、ということは、今回のループ内の命令数は 6ということになります。すべての命令が1サイクルで処理できるのであれば、実行サイクル数も6であるべきですが、結果は9サイクルでした。どこかで余分に3サイクル使用していることになります。原因はどこにあるでしょうか。
floatの計算には時間が掛かる
結論を言うと、madd.s命令の実行結果が出るのに時間が掛かっています。より正確には、madd.s命令自体は1サイクルで次の処理に進めますが、演算結果の値はその時点ではまだ利用できないのです。
ソースコードを確認してみましょう。該当箇所は以下のようになっています。
// ～～ 省略 ～～
madd.s f0, f1, f2 // f0 += f1 * f2
ssi f0, a2, 0 // dst[0] = f0;
// ～～ 省略 ～～
madd.s命令でf0レジスタに演算結果が入りますが、次行ですぐにssi命令でf0の値を使おうとしています。このときssi命令はf0の値が利用可能になるまで待機させられてしまうのです。
本当でしょうか？ 試しに、何も処理を行わないnop命令を間に挟んで、実行結果がどう変化するか確認してみましょう。
// ～～ 省略 ～～
madd.s f0, f1, f2 // f0 += f1 * f2
nop
nop
nop
ssi f0, a2, 0 // dst[0] = f0;
// ～～ 省略 ～～
nop命令の個数が0～3個では実行結果が変わらず、4個目からは実行結果が伸びるはずです。
このことから、madd.s命令とssi命令の間には、3サイクルの待機時間が発生していることが確認できます。
手動で最適化を行う
ここまで前置きが長かったですが、ようやく今回の本題である最適化に着手します。
現状は6命令の実行に9サイクル掛かっています。madd.s命令の後に発生する3サイクルの待機を削減して6サイクルで動作するよう最適化することを目指してみましょう。実行結果が変わらない範囲で命令の実行順序などを調整し、madd.s命令とssi命令の間に他の命令を3個挟むことができれば最適化に成功したと言えます。
コードを眺めて影響の無さそうなところを探してみましょう。
最適化ポイント1:単純な実行順序の入れ替え
まず見つかるのは以下の箇所かと思います。
// ～～ 省略 ～～
loop a4, FUNC1_LOOP_END
lsi f1, a3, 0
lsi f0, a2, 0
addi a3, a3, 4 // ← これと…
madd.s f0, f1, f2 // ← これを入れ替えても良いのでは？
ssi f0, a2, 0
addi a2, a2, 4
FUNC1_LOOP_END:
// ～～ 省略 ～～
なるほどaddi a3, a3, 4の処理結果であるa3はループの最初でしか利用していませんから、後ろにずらしても影響は無さそうです。実際に入れ替えて実行結果を確かめてみましょう。
total : 32764 / 4096 = 8.00 cycle
結果は上記の通りです。まずは簡単に1サイクル削減できました。
あと2サイクル、いけるでしょうか。
最適化ポイント2:順序の入れ替えに伴う演算結果の相殺
他の命令は、単純に実行順序を入れ替えるだけでは結果が変わってしまいます。
では、以下のようにするのはどうでしょうか。
// ～～ 省略 ～～
beqz a4, FUNC1_LOOP_END
wfr f2, a5
// ←この辺でa2の値を-4しておき…
loop a4, FUNC1_LOOP_END
lsi f1, a3, 0
lsi f0, a2, 0 // ←ここのオフセットを+4に変更すれば…
madd.s f0, f1, f2
addi a3, a3, 4
ssi f0, a2, 0 // ←これと…
addi a2, a2, 4 // ←これを入れ替えられるのでは？
FUNC1_LOOP_END:
// ～～ 省略 ～～
ssi命令とその下のaddi命令を入れ替えると効果がありそうです。しかし単に入れ替えただけでは、ssi命令の保存先アドレスa2が+4された後になりますから、ループ開始前にaddi a2, a2, -4 を実行して影響を相殺しておきます。
そうすると lsi f0, a2, 0 命令の読出しアドレスが-4されてしまいますから、こちらも影響を相殺するべく、オフセット値を0から4に変更し lsi f0, a2, 4 とするのです。
変更後のソースは以下のようになります。
// ～～ 省略 ～～
beqz a4, FUNC1_LOOP_END
wfr f2, a5
addi a2, a2, -4 // ループに入る前に dstアドレスを4減算
loop a4, FUNC1_LOOP_END
lsi f1, a3, 0
lsi f0, a2, 4 // dstアドレスから値を読み出す際にオフセットを+4指定
madd.s f0, f1, f2
addi a3, a3, 4
addi a2, a2, 4 // ssi命令より前にdstアドレスを4加算
ssi f0, a2, 0
// ～～ 省略 ～～
FUNC1_LOOP_END:
これを実行すると、以下のようになりました。
total : 28676 / 4096 = 7.00 cycle
うまくいきました。
残りあと1サイクル、どうにかできるでしょうか。
最適化ポイント3:実行順序の変更とループ外での事前実行
ループの最初にある命令を後ろに移動するのはどうでしょうか。
// ～～ 省略 ～～
beqz a4, FUNC1_LOOP_END
wfr f2, a5
addi a2, a2, -4
// ←ここで事前にlsi f1,a3,0を実行しておけば…
loop a4, FUNC1_LOOP_END
lsi f1, a3, 0 // ←これを…
lsi f0, a2, 4
madd.s f0, f1, f2
addi a3, a3, 4
addi a2, a2, 4
// ←ここに移せるのでは？
ssi f0, a2, 0
FUNC1_LOOP_END:
// ～～ 省略 ～～
これは、次回のループで使うf1レジスタの値をssi命令の前の時点で用意しておくという考え方です。単に位置を変えただけでは初回のループで使うf1レジスタの値が不定になってしまうので、ループ開始前にもf1レジスタの値を用意するのです。
変更後のソースは以下のようになります。
test_func1:
// ～～ 省略 ～～
beqz a4, FUNC1_LOOP_END
wfr f2, a5
addi a2, a2, -4
lsi f1, a3, 0 // ループの前にも f1 = src[0]; を行う。
loop a4, FUNC1_LOOP_END
// ←ここにあったlsi命令をループの前とssiの前に配置
lsi f0, a2, 4
madd.s f0, f1, f2
addi a3, a3, 4
addi a2, a2, 4
lsi f1, a3, 0 // 次回ループで使用するための f1 = src[0]; を行う
ssi f0, a2, 0
FUNC1_LOOP_END:
// ～～ 省略 ～～
実行結果は以下のようになりました。
total : 24580 / 4096 = 6.00 cycle
うまく行きました。
6命令をきっちり6サイクルでループ処理できました。手作業での最適化は成功です！
もっと速くするには
ここまで見てきて気が付いた方もいると思いますが、ループ一回の処理の6命令の内訳を見ると、読書き演算で4命令、srcとdstの加算で2命令使っています。これを例えば一回のループで4回分の処理をし、ループの回数を1/4に減らせればどうでしょうか。
つまりループ内の内訳を以下のように…
- 読書き演算の 4命令×4回分 = 16命令
-
srcとdstの加算 = 2命令
合計で18命令を 1/4のループ回数で実行できるのではないでしょうか。
lsiとssiはオフセットを設定できますから、srcとdstの加算をせずに4個先までアクセスすることは問題ありません。浮動小数レジスタはf0～f15の合計16個ありますが現状では3個しか使用していませんから、レジスタの余裕もあります。
ループ一回あたりの処理要素数を増やす
ということで、変更後の関数が test_func2 になります。
__attribute((noinline, noclone, weak, optimize("-O3")))
void test_func2(float* dst, const float* src, uint32_t len, float vol) {
for (uint32_t i = 0; i < len; i+=4) {
dst[i+0] += src[i+0] * vol;
dst[i+1] += src[i+1] * vol;
dst[i+2] += src[i+2] * vol;
dst[i+3] += src[i+3] * vol;
}
}
ループ一回当たり配列4要素分を処理し、ループカウンタを+=4とすることでループ自体の回数が1/4に減っています。なおlenの値が4の倍数でないと端数の扱いに問題が生じますが、今回は呼び出し元が必ず4の倍数の長さの配列を用意するものとして、端数は考慮しないものとします。
計測処理の中を変更して、test_func1の代わりにtest_func2を呼び出すようにして実行結果を見てみましょう。
total : 28770 / 4096 = 7.00 cycle
C/C++版でも元の状態よりちょっと速くなりました。…ですが効果はいまひとつのようです。
ここで先ほどの最適化を思い出してください。処理速度を低下させていたのは、演算結果をすぐに配列に書き込もうとして待機が発生していたことでした。ということは、C/C++版でも以下のようにすれば余分な待機が発生しなくなるのではないでしょうか。
__attribute((noinline, noclone, weak, optimize("-O3")))
void test_func2(float* dst, const float* src, uint32_t len, float vol) {
for (uint32_t i = 0; i < len; i+=4) {
float result0 = dst[i+0] + src[i+0] * vol;
float result1 = dst[i+1] + src[i+1] * vol;
float result2 = dst[i+2] + src[i+2] * vol;
float result3 = dst[i+3] + src[i+3] * vol;
dst[i+0] = result0;
dst[i+1] = result1;
dst[i+2] = result2;
dst[i+3] = result3;
}
}
要点は、演算結果をすぐに使わず一旦変数にとっておき、配列4個分の演算を行った後でdst配列に書き込むことです。これにより、演算結果が利用可能になるまでの待ち時間の間に他の演算を行えるのではないでしょうか。
test_func2の中身を上記のものに変更して、実行結果を見てみましょう。
total : 18444 / 4096 = 4.50 cycle
予想通り、効果がありました。さきほど手動で最適化した結果よりも速くなっています。
test_func2も Compiler Explorer を使って、アセンブリ言語化されたソースコードを入手してみましょう。
興味がある方は、上記の変更の前後で出力結果がどう変わるかを試してみてください。
変更後のoptimize("-O3")での結果は以下の通りです。
test_func2(float*, float const*, unsigned long, float):
entry sp, 32
wfr f0, a5 // f0に引数volをセット
beqz a4, .L1
addi a8, a4, -1
slli a8, a8, 2
srli a8, a8, 4
addi a8, a8, 1
loop a8, .L3_LEND // ここからa8の回数、.L3_LENDまでをループ
lsi f8, a3, 0 // f8 = src[0];
lsi f7, a3, 4 // f7 = src[1];
lsi f6, a3, 8 // f6 = src[2];
lsi f5, a3, 12 // f5 = src[3];
lsi f4, a2, 0 // f4 = dst[0];
lsi f3, a2, 4 // f3 = dst[1];
lsi f2, a2, 8 // f2 = dst[2];
lsi f1, a2, 12 // f1 = dst[3];
madd.s f4, f8, f0 // f4 += src[0] * vol;
madd.s f3, f7, f0 // f3 += src[1] * vol;
madd.s f2, f6, f0 // f2 += src[2] * vol;
madd.s f1, f5, f0 // f1 += src[3] * vol;
ssi f4, a2, 0 // dst[0] = f4;
ssi f3, a2, 4 // dst[1] = f3;
ssi f2, a2, 8 // dst[2] = f2;
ssi f1, a2, 12 // dst[3] = f1;
addi a3, a3, 16 // srcのアドレスを16加算
addi a2, a2, 16 // dstのアドレスを16加算
.L1: // このラベルの箇所に .L3_LEND も存在しているはず
retw
madd.s命令とssi命令の間には3命令、他の処理が行われています。ループ内の命令数を数えてみると18命令あります。実行結果に表示されるサイクル数は4.5ですが、ループ回数が1/4になっているので実行結果を4倍し 4.5 × 4 = 18サイクルが実際のサイクル数となります。うまく最適化されていますね。余分な待機も発生していないので、これを手作業で最適化をする必要はなさそうです。
…しかし、手作業での最適化の出番がない、となると本記事としては目標がなくなって迷子になってしまいます。これをさらに速くするアプローチはあるでしょうか？
ESP32-S3のSIMD命令を使ってみる
というわけで、ESP32-S3にて使用可能になったSIMD命令を試してみることにします。
ここから先は普通のESP32では試すことができません。ESP32-S3にてお試しください。
SIMD命令は複数個のレジスタに同一の演算処理を同時に行うことができる命令セットです。今回の例にも適用できる『浮動小数レジスタ4個を対象にメモリの読書きをする命令』が用意されていますので、これを使用してみます。
なおC/C++のコードをESP32-S3用にコンパイルしても、SIMD命令を使ったコードは出てこないようです。今回はこちらで作成したコードを載せておきます。
#if __has_include (<sdkconfig.h>)
#include <sdkconfig.h>
#endif
#if CONFIG_IDF_TARGET_ESP32S3
// void test_func2(float* dst, const float* src, uint32_t len, float vol);
// a2 = float* dst 出力先データ
// a3 = float* src 加算元データ
// a4 = uint32_t len ループ回数
// a5 = float vol 倍率
.global test_func2
.section .text
.align 4
test_func2:
entry sp, 32
srli a4, a4, 2 // ループ回数 >>= 2 (1/4にする)
beqz a4, FUNC2_LOOP_END // ループ回数が0なら終了処理へジャンプ
wfr f0, a5 // 倍率を f0 レジスタにセット
loop a4, FUNC2_LOOP_END // ループ開始
ee.ldf.128.ip f5, f6, f7, f8, a3, 16 // 加算元から float 4個 読み、a3 アドレスを 16 加算
ee.ldf.128.ip f1, f2, f3, f4, a2, 0 // 出力先から float 4個 読み
madd.s f1, f5, f0 // f1 += f5 * f0
madd.s f2, f6, f0 // f2 += f6 * f0
madd.s f3, f7, f0 // f3 += f7 * f0
madd.s f4, f8, f0 // f4 += f8 * f0
ee.stf.128.ip f1, f2, f3, f4, a2, 16 // 出力先へ float 4個 書き、 a2 アドレスを 16 加算
FUNC2_LOOP_END:
retw
#endif
先ほど使ったアセンブリ言語ファイルesp32_asm_003_1.Sの末尾に上記のコードを追加してください。
C/C++版のtest_func2にはweak属性が指定されていますので、これを追加するだけでアセンブリ言語版のtest_func2が有効になります。
さて新たな命令が出てきたので少し説明します。
ee.ldf.128.ip 命令
float値をメモリから読み取る命令です。以下の動作を1サイクルで行います。
- 指定されたアドレスから16Byte読み込み、4個の浮動小数レジスタに格納する
- アドレス指定に使用したレジスタの値を加算する
今回の例の ee.ldf.128.ip f5, f6, f7, f8, a3, 16 であれば、f5 f6 f7 f8の4個にa3のアドレスから16Byte読み込み、さらにa3の値を +16 加算する、という動作になります。
ee.stf.128.ip 命令
float値をメモリに書き込む命令です。以下の動作を1サイクルで行います。
- 4個の浮動小数レジスタの内容を指定されたアドレスに16Byte書き込む
- アドレス指定に使用したレジスタの値を加算する
今回の例の ee.stf.128.ip f1, f2, f3, f4, a2, 16 であれば、f1 f2 f3 f4の4個の内容をa2のアドレスへ16Byte書き込み、さらにa2の値を +16 加算する、という動作になります。
浮動小数レジスタを4並列で同時に演算する命令は？
探してみましたが…見つけることができませんでした。おそらく用意されていないようです。
とはいえメモリの読書きだけでも4並列で実行できれば、実行命令数は大幅に削減できます。
SIMD命令版を実行してみる
実行結果は以下の通りです。
total : 10260 / 4096 = 2.50 cycle
かなり速くなりました。実際のループ一回あたりのサイクル数は、 2.5 × 4 = 10サイクルです。
しかし、ループ内の命令数を数えてみると、7命令しかありません。7命令に10サイクル掛かっているということは、つまり今回も madd.s 命令の直後に ee.stf.128.ip命令で演算結果の値を使おうとして3サイクル待機させられています。
先ほど手動で行った最適化手法が使えるでしょうか？ しかし、レジスタの読書きとアドレスの加算を同時に行う命令を使用しているため、順序を入れ替えられる命令がありません。1命令分は行けそうですが、3命令分となると難しいです。
1ループ中に処理する数をさらに増やす
ということで構造を変更します。1ループで処理する要素数を 8個に増やし、ループ回数を 1/8にします。この関数を呼び出す際には配列の要素数が8の倍数であることが条件となりますが、今回はそういう仕様ということにして進めます。
コードは以下のようになります。
test_func2:
entry sp, 32
srli a4, a4, 3 // ループ回数 >>= 3 (1/8にする)
beqz a4, FUNC2_LOOP_END // ループ回数が0なら終了処理へジャンプ
wfr f0, a5 // 倍率を f0 レジスタにセット
mov a15, a2 // 出力先を a15 レジスタにもセット
ee.ldf.128.ip f5, f6, f7, f8, a3, 16 // 加算元から float 4個 読み、a3 アドレスを 16 加算
ee.ldf.128.ip f1, f2, f3, f4, a2, 16 // 出力先から float 4個 読み、a2 アドレスを 16 加算
madd.s f1, f5, f0 // f1 += f5 * f0
madd.s f2, f6, f0 // f2 += f6 * f0
madd.s f3, f7, f0 // f3 += f7 * f0
madd.s f4, f8, f0 // f4 += f8 * f0
// ↑ここまでで、初回のストア命令に使用する演算結果をループ開始前に用意する
loop a4, FUNC2_LOOP_END // ループ開始
ee.ldf.128.ip f5, f6, f7, f8, a3, 16 // 加算元から float 4個 読み、a3 アドレスを 16 加算
ee.ldf.128.ip f9, f10,f11,f12,a2, 16 // 出力先から float 4個 読み、a2 アドレスを 16 加算
madd.s f9, f5, f0 // f9 += f5 * f0
madd.s f10,f6, f0 // f10+= f6 * f0
madd.s f11,f7, f0 // f11+= f7 * f0
madd.s f12,f8, f0 // f12+= f8 * f0
ee.stf.128.ip f1, f2, f3, f4, a15,16 // 出力先へ float 4個 書き、 a15 アドレスを 16 加算
// ↑このストア命令は、以前に求めておいた結果を使用している
ee.ldf.128.ip f5, f6, f7, f8, a3, 16 // 加算元から float 4個 読み、a3 アドレスを 16 加算
ee.ldf.128.ip f1, f2, f3, f4, a2, 16 // 出力先から float 4個 読み、a2 アドレスを 16 加算
madd.s f1, f5, f0 // f1 += f5 * f0
madd.s f2, f6, f0 // f2 += f6 * f0
madd.s f3, f7, f0 // f3 += f7 * f0
madd.s f4, f8, f0 // f4 += f8 * f0
ee.stf.128.ip f9, f10,f11,f12,a15,16 // 出力先へ float 4個 書き、 a15 アドレスを 16 加算
// ↑このストア命令は、以前に求めておいた結果を使用している
FUNC2_LOOP_END:
retw
行数がかなり増えましたが、要点としては以下の3点です。
- ループの中では、2回分の処理を行う
- 保存処理は1つ前のデータを扱う
-
loop命令の前に初回分のデータを用意する
こうすると、madd.s命令の結果が出るまでに他の処理を行えます。浮動小数レジスタを贅沢に使用していますが、デメリットは特にありませんので遠慮せず使って行きましょう。
実行結果は以下の通りです。
total : 7200 / 4096 = 1.76 cycle
結果の表示に端数が出て1.76になっていますが、これは1.75 cycleとして扱ってください。
ループ回数は 1/8 になっているので、 1.75 × 8 = 14 サイクルがループ一回あたりのサイクル数になります。ループ内の命令数も14ですから、余分な待機は発生しておらず、うまく最適化できています。
最適化前の test_func1 を -O3 でコンパイルした実行結果は 9 cycleでしたから、そこから比較すると処理時間を1/5にまで短縮できたことになります。
まとめ
以上、C/C++で作成した関数を、アセンブリ言語化して手作業で最適化する手法についての説明でした。
例として使用した test_func1とtest_func2の差や、test_func2の構造の工夫および__restrict__の指定については、C/C++の記述上だけでも処理時間を半分に短縮できるテクニックと言えますが、どうして実行時間に差が出るのかを正しく把握するためには、アセンブリ言語のソースコードを追って実際に何が起きているのか確認することが必要になってきます。今回取り上げた浮動小数計算の待ち時間などはC/C++の知識だけでは想像しづらいものですから、ぜひ皆さんもアセンブリ言語で色々試して体験してみてください。
次回は、引き続き最適化の具体的な手法を紹介しようと思います。
Discussion