目次: 自宅サーバー
Dockerを使っていると要らなくなったイメージ、コンテナ、ビルドキャッシュがたまってきて、/varディレクトリ以下が肥大化していることがあります。いつも忘れてしまうお掃除用のコマンドをメモしておきます。
### 終了しているcontainerの削除 $ docker container prune ### 確認 $ docker container ls -a ### タグもなく使われていないimageの削除 $ docker image prune ### 確認 $ docker image ls -a ### build cacheの削除 $ docker builder prune ### 確認 $ docker builder ls
Dockerがディスク容量をどの程度使用しているのかについてはsystem dfが便利(docker system dfのマニュアル)です。
$ docker system df TYPE TOTAL ACTIVE SIZE RECLAIMABLE Images 1 1 14.26GB 0B (0%) Containers 2 1 0B 0B Local Volumes 0 0 0B 0B Build Cache 17 0 0B 0B
昔はコマンドの名前に一貫性がなかった記憶がありますが、今はxxx lsとすれば大抵の場合は一覧が出るため統一感があります。私のようなライトユーザーがやりたいと思う程度の機能は、大抵既に存在しており良くできていてありがたいです。
目次: ベンチマーク
お手軽最適化のメモ、昨日の続きです。行列の掛け算を題材にします。前回は行列の掛け算と素朴な実装のコードを紹介しました。今回はお手軽最適化を紹介します。
スカラー処理だと遅いけれど、お手軽に最適化(数倍程度)がしたいときの参考になれば幸いです。
素朴版のコードですとi, j, kの順でループになっていて、kを最内ループにしていました。ループ内の計算は、
c[i * nn + j] += a[i * kk + k] * b[k * nn + j]
でした。このときメモリアクセスのパターンは、
です。GCCの自動ベクトル化ですと、このアクセスパターンをうまく最適化できないようです。対象が最内ループのみなのかもしれません。ループを入れ替えjを最内にして、BとCを行方向に読むようにします。するとメモリアクセスのパターンは、
です。ループを入れ替えるとループ内でCの0初期化ができないので、ループの外に追い出して最終的に下記のようなコードになります。
void sgemm_inner(const float *a, const float *b, float *c, int mm, int nn, int kk)
{
for (int i = 0; i < mm; i++) {
for (int j = 0; j < nn; j++) {
c[i * nn + j] = 0.0f;
}
}
for (int i = 0; i < mm; i++) {
for (int k = 0; k < kk; k++) {
for (int j = 0; j < nn; j++) { //★★jが最内ループ★★
c[i * nn + j] += a[i * kk + k] * b[k * nn + j];
}
}
}
}
$ gcc -Wall -g -O2 -fno-tree-vectorize -static -march=znver3 sgemm.c $ ./a.out matrix size: M:1519, N:1517, K:1523 time: 1.528314 (参考: 素朴版の実行時間) time: 2.277758 (参考: OpenBLASシングルスレッドの実行時間) $ export OPENBLAS_NUM_THREADS=1 ----- use CBLAS verify: 0.052149
ループ入れ替えで倍くらい速くなっていますが、この最適化の本領はコンパイラの自動ベクトル化です。GCCならば -ftree-vectorizeオプションを指定すると、行列Bと行列CへのアクセスにSIMD命令を使うようになります。
Ryzen 7 5700Xの場合はAVX2命令を使えます。他のCPUをお使いの場合は -marchを適宜変更してください。
$ gcc -Wall -g -O2 -ftree-vectorize -static -march=znver3 sgemm.c $ ./a.out matrix size: M:1519, N:1517, K:1523 time: 0.181133 (参考: OpenBLASシングルスレッドの実行時間) $ export OPENBLAS_NUM_THREADS=1 ----- use CBLAS verify: 0.052149
素朴版とOpenBLASでは1/43もの差がありましたが、ループ入れ替えと自動ベクトル化によってOpenBLASの1/3.5程度まで近づきました。GEMMが計算偏重の処理で最適化の効果が出やすい、という点を考慮する必要はあるものの僅かな書き換えで得られる効果にしては割と良いのではないでしょうか。
ソースコードを書き換える元気があれば、別の最適化方法もあります。GEMM特有の話に近付いてしまい、汎用的な話から遠ざかりますが、最適化ポイントの例という意味では参考になるはず……です。たぶん。
今回はSIMD命令でjの方向に一気に読むまでは同じですが、jの方向に進めるのではなく、kの方向に進めて、計算結果をCに足していく戦略です。具体的に言えばAi,kとBk,j〜Bk,j+7の8要素を一気に掛け算してCi,j〜Ci,j+7へ一気に足します。なぜ8要素かというとAVX/AVX2のレジスタ長(256bit)を使うと、32bit長のfloatを一度に8要素処理できるためです。
SIMD命令にはAi,kをSIMDレジスタの全要素に配る命令(set1_psから生成されるbroadcast命令)や、掛け算と足し算を一度に行うfmadd命令など、この計算順に最適な命令が揃っています。
AVX/AVX2のIntrinsicsの詳細についてはIntelのサイトなどを見ていただくとして(Intel Intrinsics Guide)、コードは下記のようになります。
#include <immintrin.h>
void sgemm_avx(const float *a, const float *b, float *c, int mm, int nn, int kk)
{
for (int j = 0; j < nn;) {
if (nn - j >= 8) {
for (int i = 0; i < mm; i++) {
__m256 vc = _mm256_set1_ps(0.0f);
for (int k = 0; k < kk; k++) {
__m256 va = _mm256_set1_ps(a[i * kk + k]);
__m256 vb = _mm256_loadu_ps(&b[k * nn + j]);
vc = _mm256_fmadd_ps(va, vb, vc);
}
_mm256_storeu_ps(&c[i * nn + j], vc);
}
j += 8;
} else {
for (int i = 0; i < mm; i++) {
c[i * nn + j] = 0.0f;
for (int k = 0; k < kk; k++) {
c[i * nn + j] += a[i * kk + k] * b[k * nn + j];
}
}
j++;
}
}
}
$ gcc -Wall -g -O2 -static -march=znver3 sgemm.c $ ./a.out matrix size: M:1519, N:1517, K:1523 time: 0.384861 (参考: 素朴版の実行時間) time: 2.277758 (参考: ループ入れ替え版+自動ベクトル化の実行時間) time: 0.181133 (参考: OpenBLASシングルスレッドの実行時間) $ export OPENBLAS_NUM_THREADS=1 ----- use CBLAS verify: 0.052149
素朴版と比べると6倍速いですが、ループ入れ替え+自動ベクトル化には負けています。
先程のコードはSIMDレジスタを3個しか同時に使っていませんでした。AVX/AVX2のYMMレジスタは16個もあるのに3個しか使わないのはもったいですから、iのループを8要素ずつアンローリングしてSIMDレジスタを同時にたくさん使いましょう。レジスタをうまく使いまわせば12要素のアンローリング(Bの保持に1個、Cの保持に12個、Aの保持に1個、計14個)まではできそうです。たぶん。
#include <immintrin.h>
void sgemm_avx_unroll8(const float *a, const float *b, float *c, int mm, int nn, int kk)
{
for (int j = 0; j < nn;) {
if (nn - j >= 8) {
int i = 0;
for (; i < (mm & ~7); i += 8) {
__m256 vc0 = _mm256_set1_ps(0.0f);
__m256 vc1 = _mm256_set1_ps(0.0f);
__m256 vc2 = _mm256_set1_ps(0.0f);
__m256 vc3 = _mm256_set1_ps(0.0f);
__m256 vc4 = _mm256_set1_ps(0.0f);
__m256 vc5 = _mm256_set1_ps(0.0f);
__m256 vc6 = _mm256_set1_ps(0.0f);
__m256 vc7 = _mm256_set1_ps(0.0f);
for (int k = 0; k < kk; k++) {
__m256 vb = _mm256_loadu_ps(&b[k * nn + j]);
__m256 va0 = _mm256_set1_ps(a[(i+0) * kk + k]);
__m256 va1 = _mm256_set1_ps(a[(i+1) * kk + k]);
__m256 va2 = _mm256_set1_ps(a[(i+2) * kk + k]);
__m256 va3 = _mm256_set1_ps(a[(i+3) * kk + k]);
__m256 va4 = _mm256_set1_ps(a[(i+4) * kk + k]);
__m256 va5 = _mm256_set1_ps(a[(i+5) * kk + k]);
__m256 va6 = _mm256_set1_ps(a[(i+6) * kk + k]);
__m256 va7 = _mm256_set1_ps(a[(i+7) * kk + k]);
vc0 = _mm256_fmadd_ps(va0, vb, vc0);
vc1 = _mm256_fmadd_ps(va1, vb, vc1);
vc2 = _mm256_fmadd_ps(va2, vb, vc2);
vc3 = _mm256_fmadd_ps(va3, vb, vc3);
vc4 = _mm256_fmadd_ps(va4, vb, vc4);
vc5 = _mm256_fmadd_ps(va5, vb, vc5);
vc6 = _mm256_fmadd_ps(va6, vb, vc6);
vc7 = _mm256_fmadd_ps(va7, vb, vc7);
}
_mm256_storeu_ps(&c[(i+0) * nn + j], vc0);
_mm256_storeu_ps(&c[(i+1) * nn + j], vc1);
_mm256_storeu_ps(&c[(i+2) * nn + j], vc2);
_mm256_storeu_ps(&c[(i+3) * nn + j], vc3);
_mm256_storeu_ps(&c[(i+4) * nn + j], vc4);
_mm256_storeu_ps(&c[(i+5) * nn + j], vc5);
_mm256_storeu_ps(&c[(i+6) * nn + j], vc6);
_mm256_storeu_ps(&c[(i+7) * nn + j], vc7);
}
for (; i < mm; i++) {
__m256 vc = _mm256_set1_ps(0.0f);
for (int k = 0; k < kk; k++) {
__m256 vb = _mm256_loadu_ps(&b[k * nn + j]);
__m256 va = _mm256_broadcast_ss(&a[(i+0) * kk + k]);
vc = _mm256_fmadd_ps(va, vb, vc);
}
_mm256_storeu_ps(&c[i * nn + j], vc);
}
j += 8;
} else {
for (int i = 0; i < mm; i++) {
c[i * nn + j] = 0.0f;
for (int k = 0; k < kk; k++) {
c[i * nn + j] += a[i * kk + k] * b[k * nn + j];
}
}
j++;
}
}
}
$ ./a.out matrix size: M:1519, N:1517, K:1523 time: 0.108420 (参考: ループ入れ替え版+自動ベクトル化の実行時間) time: 0.181133 (参考: OpenBLASシングルスレッドの実行時間) $ export OPENBLAS_NUM_THREADS=1 ----- use CBLAS verify: 0.052149
もはや最適化前のコードの原型がありませんが、ループ入れ替え版+自動ベクトル化の1.7倍くらいの速度になりました。OpenBLASの1/2程度まで迫っています。この最適化手法が汎用的か?と聞かれると何とも言えないですが、SIMDレジスタを同時にたくさん使う、最内ループ以外もアンローリング(最内ループはコンパイラがやってくれる)辺りは割と汎用的なアイデアです。
あと前回言った通り、GEMMは最適化の題材として取り上げただけなので、実際にGEMMを計算する場合はこのコードや自分で書いたコードを使うのではなく、信頼と実績のOpenBLASを使ってくださいませ。
目次: ベンチマーク
お手軽最適化のメモです。行列の掛け算を題材にします。
最初にお断りしておくとGEMMのような汎用処理の場合は、自分で最適化せずにOpenBLASを使ってください(素人が最適化しても勝てません)。しかしOpenBLASのような限界まで最適化されたコードは誰でも簡単には書ける、とは言えません。
スカラー処理だと遅いけれど、お手軽に最適化(数倍程度)がしたいときの参考になれば幸いです。
GEMMはGeneral Matrix Multiplyの略で、高校数学辺りでやった(はず)の行列の掛け算のことです。floatの場合はSGEMMと呼ばれ、doubleの場合はDGEMMと呼ばれます。最適化の題材はどちらでも良いんですけど、今回はSGEMMを使います。
忘れている方のために2行3列の行列Aと3行2列Bの掛け算A x B = Cだとこんな感じです。
Aの列数とBの行数は一致していなければなりません。行数と列数の関係を表すと、行列A(M行K列)x行列B(K行N列)= 行列C(M行N列)となります。
Cの1要素を計算するには、Aの1行とBの1列が必要です。式、および、視覚的に示すと下記のようになります。
説明はこれくらいにしてコードを見ましょう。
SGEMMを素直にコードにするとこんな感じです。行方向にデータを格納(Row-major orderといいます)しているので、N列の行列Cのi行j列(以降Ci,jと書く)にアクセスする際はc[i * N + j] とします。
void sgemm_naive(const float *a, const float *b, float *c, int mm, int nn, int kk)
{
for (int i = 0; i < mm; i++) {
for (int j = 0; j < nn; j++) {
c[i * nn + j] = 0.0f;
for (int k = 0; k < kk; k++) {
c[i * nn + j] += a[i * kk + k] * b[k * nn + j];
}
}
}
}
行列のサイズを適当に設定(M = 1519, N = 1517, K = 1523)して、実行時間を測ります。CPUはRyzen 7 5700Xです。
$ gcc -Wall -g -O2 -static sgemm.c $ ./a.out matrix size: M:1519, N:1517, K:1523 time: 2.277758
実行時間は実行するたびに変わりますが、大体2.27秒くらいでしょうか。OpenBLASのシングルスレッド(環境変数OPENBLAS_NUM_THREADS=1にするとシングルスレッド動作になります)で計算した時間を見ると、
c_ex = malloc(m * n * sizeof(float) * 2);
// C = alpha AB + beta C
float alpha = 1.0f, beta = 0.0f;
int lda = k, ldb = n, ldc = n;
printf("----- use CBLAS\n");
gettimeofday(&st, NULL);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
m, n, k, alpha, a, lda, b, ldb, beta, c_ex, ldc);
gettimeofday(&ed, NULL);
timersub(&ed, &st, &ela);
printf("verify: %d.%06d\n", (int)ela.tv_sec, (int)ela.tv_usec);
$ export OPENBLAS_NUM_THREADS=1 $ gcc -Wall -g -O2 -static -L path/to/openblas/ sgemm.c -lopenblas $ ./a.out matrix size: M:1519, N:1517, K:1523 ----- use CBLAS verify: 0.052149
わずか0.05秒、実に43倍という驚異のスピードです。すごいですね……。
行列の掛け算の説明でほぼ終わってしまいました。お手軽最適化はまた次に。
東京の道はようわからん通称(明治通り、みたいなやつ)が付いています。東京育ちの人にはなじみ深い名前だと思うんですが、都外出身者からすると、東京の道の通称と国道何号線、都道何号線の区分が全く合っていないので、知らない場所の道の通称を言われると結構困ります。
なぜかというと国道何号線、都道何号線という表記は地図で省かれることはない一方で、東京の道の通称は高速道路や地下鉄と重なったときに省かれてしまうことがあるためです。地図を見ても「〇〇通り?何それ?どこ??」と困惑します。
東京の道の通称は東京都が決めています(東京都通称道路名〜道路のわかりやすく親しみやすい名称〜 - 東京都建設局)。道案内の青看板に載るような名前と思ってもらえばわかりやすいでしょう。
東京都がまとめてくれているのはありがたいんですが「東京都が公式に定めた通称」というのは不思議な響きで、それは公称ではなかろうか?通称とは一体……??
東京の道の通称には「年号」+「通り」という名前がいくつかあります。このうちなぜか大正通りだけは存在しません。不思議ですね。
調べてみると割と有名な話らしいです(「明治通り」「昭和通り」はあるのに、なぜか「大正通り」はない東京のちょっとした謎 - アーバン ライフ メトロ)。詳しくは記事を読んでいただくとして、簡単に言えば、
ということみたいです。
目次: RISC-V
以前(2022年12月22日の日記参照)の日記でCoreMarkのスコアを測って表にしました。実はCoreMarkはOfastのみでは最速にはならず、コンパイルオプションをガチガチにチューンすると結構差が出ます。実際にNSITEXE NS31Aの測定結果でお見せしたいと思います。
どうしてNS31Aかというと、非常にシンプルなCPU&自分の会社で作っているので素性が明確であるためです。複雑なCPU、中身の分からないCPUになればなるほど、総当たりでオプションの組み合わせを試す不毛な作業になりがちです。今回はそういう組み合わせ問題を解きたいわけじゃないんで、簡単な奴で行きます。
まずはベースとなるOfastの結果です。実はO3でも結果は同じです。
2K performance run parameters for coremark. CoreMark Size : 666 Total ticks : 18912 Total time (secs): 18.912000 Iterations/Sec : 58.164129 Iterations : 1100 Compiler version : GCC12.2.0 Compiler flags : -Ofast -gdwarf-4 -march=rv32im -mabi=ilp32 -mcmodel=medany Memory location : Please put data memory location here (e.g. code in flash, data on heap etc) seedcrc : 0xe9f5 [0]crclist : 0xe714 [0]crcmatrix : 0x1fd7 [0]crcstate : 0x8e3a [0]crcfinal : 0x33ff Correct operation validated. See README.md for run and reporting rules. CoreMark 1.0 : 58.164129 / GCC12.2.0 -Ofast -gdwarf-4 -march=rv32im -mabi=ilp32 -mcmodel=medany / Heap
動作周波数は25MHzですので、58.164129 / 25 = 2.326 CM/MHz です。
最適化の基本となる、ループアンローリング、インライン化(-funroll-all-loops, -finline-functions)を足します。
キャッシュラインが32バイトなので、関数の先頭を32バイト境界に配置します(-falign-functions=32)。関数先頭で命令キャッシュミスヒットが発生したときに、同じキャッシュラインに後続の命令(1ラインに32 / 4 = 8命令)が載ります。後続の命令フェッチがキャッシュヒットすれば、最初のミスヒットを挽回できるだろうという目的です。
ジャンプやループの際に実行しない命令が中途半端にキャッシュに取り込まれないよう(= 利用効率の向上)、ジャンプやループの位置は8バイト境界に配置します(-falign-jumps=8 -falign-loops=8)。これも32バイト境界にすべきかと思いましたが、コード領域が散逸しすぎるためか逆に遅いです。
基本的に関数はインライン化した方がcall, retを省略、レジスタ共用など全体的に最適化できて速いです。しかしNS31Aは命令キャッシュが小さめ(FPGA向けコンフィグでは16KB)なので、無差別に関数をインライン化すると命令キャッシュがあふれてキャッシュミスヒットが発生してしまい、逆に遅くなります。
従ってあまりにも大きな関数はインライン化しないように設定します(-finline-limit=300)。デフォルト値600の1/2にしています(※)。
2K performance run parameters for coremark. CoreMark Size : 666 Total ticks : 15819 Total time (secs): 15.819000 Iterations/Sec : 69.536633 Iterations : 1100 Compiler version : GCC12.2.0 Compiler flags : -Ofast -gdwarf-4 -march=rv32im -mabi=ilp32 -mcmodel=medany -funroll-all-loops -finline-functions -finline-limit=300 -falign-functions=32 -falign-jumps=8 -falign-loops=8 Memory location : Please put data memory location here (e.g. code in flash, data on heap etc) seedcrc : 0xe9f5 [0]crclist : 0xe714 [0]crcmatrix : 0x1fd7 [0]crcstate : 0x8e3a [0]crcfinal : 0x33ff Correct operation validated. See README.md for run and reporting rules. CoreMark 1.0 : 69.536633 / GCC12.2.0 -Ofast -gdwarf-4 -march=rv32im -mabi=ilp32 -mcmodel=medany -funroll-all-loops -finline-functions -finline-limit=300 -falign-functions=32 -falign-jumps=8 -falign-loops=8 / Heap
動作周波数は25MHzですので、69.536633 / 25 = 2.781 CM/MHzです。ハードウェアは何も変えていませんが、性能1.2倍です。コンパイルオプションの威力恐るべし。
(※)この数値はGCC内部で使う仮想命令のライン数らしく、300が本当に適切か示すのは不可能です。マニュアルを見ると1/2や1/4に調整することが多いようなので、それに倣っています(参考: GCCのマニュアル)。
Twitterで見かけて面白かった問題、解法をメモしておきます。
100階建てのビルから卵を落とします。卵はある階よりも低ければ割れませんが、ある階よりも高いと割れます。卵を2つ持っているとき、卵が何階で割れるか調べる方法を提示してください。そのとき卵を落とす最大の回数をできる限り小さくしてください。
という問題です。
1階から順に落とせば確実にわかりますが、最大の投擲回数が100回(卵が割れる階数が100Fのとき)になり、最適な方法ではありません。
卵が2個あることを利用し、1個目は複数階を一気に飛ばして(例10F, 20F, 30F, ...)投擲します。例えば10F飛ばしで試すとすると10F, 20F, 30Fで割れず40Fで割れたら、2個目は「最後に割れなかった階の1つ上の階」すなわち31Fから投擲します。このように試すともっと早く卵が割れる階数がわかります。
卵の投擲回数の最大値を計算すると、例えば10F飛ばしであれば1個目が最大10回(10F, 20F, 30F, ..., 90F, 100F)、2個目が最大で9回(x1F〜x9F)なので、19回(卵が割れる階数が99Fのとき)です。100回と比べるとだいぶ少なくなりましたが、まだ無駄が残っています。
固定階数を飛ばす方法だと1個目と2個目の最大の投擲回数の和を見たとき、上の階になるほど悪化します。例えば、
1個目の投擲方法を変えて最初は10F飛ばし、次は9F飛ばし、その次は8F飛ばし、のようにすると投擲回数の最大値が悪化しないことに気づくと思います。
こんな表で示すとわかりやすいでしょうか。
1個目を投擲する階数 | 次に何階飛ばすか | 1個目の投擲回数 | 2個目の投擲開始階 | 2個目の最大の投擲回数 | 最大の投擲回数 |
---|---|---|---|---|---|
100 | 12 | 98 | 2 | 14 | |
97 | 3 | 11 | 94 | 3 | 14 |
93 | 4 | 10 | 89 | 4 | 14 |
88 | 5 | 9 | 83 | 5 | 14 |
82 | 6 | 8 | 76 | 6 | 14 |
75 | 7 | 7 | 68 | 7 | 14 |
67 | 8 | 6 | 59 | 8 | 14 |
58 | 9 | 5 | 49 | 9 | 14 |
48 | 10 | 4 | 38 | 10 | 14 |
37 | 11 | 3 | 26 | 11 | 14 |
25 | 12 | 2 | 13 | 12 | 14 |
12 | 13 | 1 | 1 | 11 | 12 |
1個目は12F, 25F, 37F, 48F, 58F, 67F, 75F, 82F, 88F, 93F, 97F, 100Fの順に投擲します。最大で12回です。割れたら最後に成功した階の1つ上の階から投擲すると、最大14回で卵が割れる階数がわかります。
最大の回数になるケースは12通りですが、一例として卵が割れる階数が99Fの場合を示します。12, 25, ..., 100(1個目が割れる), 98, 99(2個目が割れる)の14回になります。
他には1個目を投げる階数を下記のようにしても良いです。
いずれも投擲回数は最大14回となります。
< | 2023 | > | ||||
<< | < | 03 | > | >> | ||
日 | 月 | 火 | 水 | 木 | 金 | 土 |
- | - | - | 1 | 2 | 3 | 4 |
5 | 6 | 7 | 8 | 9 | 10 | 11 |
12 | 13 | 14 | 15 | 16 | 17 | 18 |
19 | 20 | 21 | 22 | 23 | 24 | 25 |
26 | 27 | 28 | 29 | 30 | 31 | - |
合計:
本日: