コグノスケ

RISC-V 64bit Linux とシステムコールと 32bit の引数 その 1 - libc の実装

目次: RISC-V - まとめリンク

タイトルに反して実は RISC-V はあまり関係ないですが……。GNU libc と musl libc の実装を見ていたら、64bit 環境で reboot システムコールの引数 magic に、

• 0x00000000_fee1deadを渡す派（musl libc）
• 0xffffffff_fee1deadを渡す派（GNU libc）

がいることに気づきました。

musl libc と GNU libc の実装

C ライブラリの reboot の関数宣言は int reboot(int cmd) となっています。musl libc の場合 0xfee1dead をキャストせず long 型の引数に渡します。符号拡張は行われず 0x00000000_fee1dead がシステムコールの引数に渡されます。

// musl/src/linux/reboot.c

int reboot(int type)
{
	return syscall(SYS_reboot, 0xfee1dead, 672274793, type);
}


// musl/arch/riscv64/syscall_arch.h

static inline long __syscall3(long n, long a, long b, long c)
{
	register long a7 __asm__("a7") = n;
	register long a0 __asm__("a0") = a;
	register long a1 __asm__("a1") = b;
	register long a2 __asm__("a2") = c;
	__asm_syscall("r"(a7), "0"(a0), "r"(a1), "r"(a2))
}

一方の GNU libc の場合は、0xfee1dead を int にキャストしてからシステムコールに渡します。この値は負の数なので符号拡張が行われて 0xffffffff_fee1dead がシステムコールの引数に渡されます。

// glibc/sysdeps/unix/sysv/linux/reboot.c

/* Call kernel with additional two arguments the syscall requires.  */
int
reboot (int howto)
{
  return INLINE_SYSCALL (reboot, 3, (int) 0xfee1dead, 672274793, howto);
}


// glibc/sysdeps/unix/sysv/linux/riscv/sysdep.h

# define internal_syscall3(number, arg0, arg1, arg2)      		\
({ 									\
	long int _sys_result;						\
	long int _arg0 = (long int) (arg0);				\
	long int _arg1 = (long int) (arg1);				\
	long int _arg2 = (long int) (arg2);				\
									\
	{								\
	register long int __a7 asm ("a7") = number;			\
	register long int __a0 asm ("a0") = _arg0;			\
	register long int __a1 asm ("a1") = _arg1;			\
	register long int __a2 asm ("a2") = _arg2;			\
	__asm__ volatile ( 						\
	"scall\n\t" 							\
	: "+r" (__a0)							\
	: "r" (__a7), "r" (__a1), "r" (__a2)				\
	: __SYSCALL_CLOBBERS); 						\
	_sys_result = __a0;						\
	}								\
	_sys_result;							\
})

最初に気づいたときはどちらかがバグっている……？と勘違いしましたが、当たり前ですがどちらも正常に動きます。すなわち Linux 側で何か対応しているはずです。続きはまた今度。

RISC-V 64bit Linux とシステムコールと 32bit の引数 その 2 - Linux の動作

目次: RISC-V - まとめリンク

昨日の続きです。64bit Linux 環境の reboot() は magic の符号拡張をしてもしなくても正常に動くという話をしました。

この話は reboot には限らないので、もっと一般化すると「64bit 環境においてシステムコールの引数がレジスタ長（= 64bit）以下だった場合（例えば int が典型例）、どう扱うか？」となります。

動作を確認するにはデバッガで追えば早いでしょう。RISC-V 64bit 向けに、

• musl libc を使うツールチェーン（後述）
• GNU libc を使うツールチェーン（後述）
• linux
• busybox
• QEMU

を用意して起動します。符号拡張を行わない musl libc 版のツールチェーンで busybox をビルドして QEMU 起動、GDB でアタッチします。

$ qemu-system-riscv64 \
  -machine virt \
  -net none \
  -nographic \
  -chardev stdio,id=con,mux=on \
  -serial chardev:con -mon chardev=con,mode=readline \
  -kernel linux/arch/riscv/boot/Image \
  -initrd busybox/initramfs.cpio \
  -s


$ riscv64-unknown-linux-gnu-gdb linux/vmlinux

(gdb) target remote :1234

...

Linux のシステムコールの実装は下記のようになります。SYSCALL_DEFINE() とはなんぞや？というところが気になると思いますが、今はひとまずシステムコール名の頭に __do_sys_ を付けた関数が定義されると理解しておけば大丈夫です。

//linux/kernel/reboot.c

SYSCALL_DEFINE4(reboot, int, magic1, int, magic2, unsigned int, cmd,
		void __user *, arg)
{

...

ブレークポイントを __do_sys_reboot に仕掛けて QEMU 側のコンソールで busybox poweroff コマンドを実行します。

(gdb) b __do_sys_reboot
Breakpoint 1 at 0xffffffff8002d860: file kernel/reboot.c, line 702.

(gdb) c
Continuing.

Breakpoint 1, __do_sys_reboot (magic1=-18751827, magic2=672274793,
    cmd=1126301404, arg=0x4321fedc) at kernel/reboot.c:702
702     {

(gdb) p/x $a0
$1 = 0xfffffffffee1dead

システムコールの引数は負の数（つまり 0xffffffff_fee1dead）となっています。この関数に辿り着く前に符号拡張されるようです。

(gdb) bt
#0  __do_sys_reboot (magic1=-18751827, magic2=672274793, cmd=1126301404,
    arg=0x4321fedc) at kernel/reboot.c:702
#1  0xffffffff8002dab2 in __se_sys_reboot (magic1=<optimized out>,
    magic2=<optimized out>, cmd=<optimized out>, arg=<optimized out>)
    at kernel/reboot.c:700
#2  0xffffffff80002fe2 in handle_exception () at arch/riscv/kernel/entry.S:231
Backtrace stopped: frame did not save the PC

バックトレースを見ると __se_sys_reboot() という関数も呼ばれています。この関数を調べるためにブレークを掛けて再度実行します。

(gdb) b __se_sys_reboot
Breakpoint 1 at 0xffffffff8002daa0: file kernel/reboot.c, line 700.

(gdb) c
Continuing.

Breakpoint 1, __se_sys_reboot (magic1=4276215469, magic2=672274793,
    cmd=1126301404, arg=1126301404) at kernel/reboot.c:700
700     SYSCALL_DEFINE4(reboot, int, magic1, int, magic2, unsigned int, cmd,

(gdb) p/x $a0
$1 = 0xfee1dead

システムコール例外ハンドラに渡ってきた引数 magic はレジスタ a0 に入っていて、値は 0xfee1dead です。__do_sys_reboot() や __se_sys_reboot() 関数を定義する SYSCALL_DEFINE() はマクロのため、デバッガから見てもソースコードが表示されません。

不思議な SYSCALL_DEFINE() マクロの仕組みは後日調べるとして、とりあえず今は逆アセンブルします。

(gdb) disas
Dump of assembler code for function __se_sys_reboot:
=> 0xffffffff8002daa0 <+0>:     addi    sp,sp,-16
   0xffffffff8002daa2 <+2>:     sd      s0,0(sp)
   0xffffffff8002daa4 <+4>:     sd      ra,8(sp)
   0xffffffff8002daa6 <+6>:     addi    s0,sp,16
   0xffffffff8002daa8 <+8>:     sext.w  a2,a2    ★
   0xffffffff8002daaa <+10>:    sext.w  a1,a1    ★
   0xffffffff8002daac <+12>:    sext.w  a0,a0    ★
   0xffffffff8002daae <+14>:    jal     ra,0xffffffff8002d860 <__do_sys_reboot>
   0xffffffff8002dab2 <+18>:    ld      ra,8(sp)
   0xffffffff8002dab4 <+20>:    ld      s0,0(sp)
   0xffffffff8002dab6 <+22>:    addi    sp,sp,16
   0xffffffff8002dab8 <+24>:    ret

RISC-V のアセンブラを見たことない方のために補足すると、sext.w という命令は RISC-V の符号拡張命令です。また引数は a0, a1, a2, a3, a4, a5 レジスタを経由して渡されます。すなわち __se_sys_reboot() 関数は a0, a1, a2 レジスタ（引数 magic1, magic2, cmd に相当）を符号拡張して __do_sys_reboot() 関数に渡しています。

ここで当初の疑問が解消しますね。疑問の内容は下記の通りで、

• magic1 に 0xffffffff_fee1deadを渡す派（GNU libc）→ 符号拡張するが値は変わらず
• magic1 に 0x00000000_fee1deadを渡す派（musl libc）→ 符号拡張して 0xffffffff_fee1dead にする

答えは途中で符号拡張しているのでどちらでも問題なかった、というわけです。なるほどね。

補足: ツールチェーンの準備

ツールチェーンのビルドは crosstool-NG などでもできますし、面倒な場合は GitHub に Ubuntu 20.04 用の *.deb ファイルを置いたのでお使いください。

GNU libc の場合は、

• riscv64-linux-glibc-binutils.deb
• riscv64-linux-glibc-headers.deb
• riscv64-linux-glibc.deb
• riscv64-linux-glibc-gcc.deb

また musl libc の場合は、

• riscv64-linux-musl-binutils.deb
• riscv64-linux-musl-headers.deb
• riscv64-linux-musl.deb
• riscv64-linux-musl-gcc.deb

をお使いください。

RISC-V 64bit Linux とシステムコールと 32bit の引数 その 3 - SYSCALL_DEFINE マクロ

目次: RISC-V - まとめリンク

昨日の続きです。「64bit 環境においてシステムコールの引数がレジスタ長（= 64bit）以下だった場合（例えば int が典型例）、どう扱うか？」コードの解説編です。

どうして符号拡張されるのか？逆アセンブルで説明したつもりになってはいけない、逃げてはダメだ！という熱心なあなたのために、マクロを展開しながら説明します。今回の解説は SYSCALL_DEFINE4() たった 1行です。

// linux/kernel/reboot.c

SYSCALL_DEFINE4(reboot, int, magic1, int, magic2, unsigned int, cmd,
		void __user *, arg)
{

...


// linux/include/linux/syscalls.h

#define SYSCALL_DEFINE4(name, ...) SYSCALL_DEFINEx(4, _##name, __VA_ARGS__)

#define SYSCALL_DEFINEx(x, sname, ...)				\
	SYSCALL_METADATA(sname, x, __VA_ARGS__)			\
	__SYSCALL_DEFINEx(x, sname, __VA_ARGS__)

下記のように展開されます。

SYSCALL_METADATA(_reboot, 4, int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)
__SYSCALL_DEFINEx(4, _reboot, int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)

SYSCALL_METADATA() の方は関数の実体と無関係なので今回は無視して、__SYSCALL_DEFINEx() を見ます。

// linux/include/linux/syscalls.h

/*
 * The asmlinkage stub is aliased to a function named __se_sys_*() which
 * sign-extends 32-bit ints to longs whenever needed. The actual work is
 * done within __do_sys_*().
 */
#ifndef __SYSCALL_DEFINEx
#define __SYSCALL_DEFINEx(x, name, ...)					\
	__diag_push();							\
	__diag_ignore(GCC, 8, "-Wattribute-alias",			\
		      "Type aliasing is used to sanitize syscall arguments");\
	asmlinkage long sys##name(__MAP(x,__SC_DECL,__VA_ARGS__))	\
		__attribute__((alias(__stringify(__se_sys##name))));	\
	ALLOW_ERROR_INJECTION(sys##name, ERRNO);			\
	static inline long __do_sys##name(__MAP(x,__SC_DECL,__VA_ARGS__));\
	asmlinkage long __se_sys##name(__MAP(x,__SC_LONG,__VA_ARGS__));	\
	asmlinkage long __se_sys##name(__MAP(x,__SC_LONG,__VA_ARGS__))	\
	{								\
		long ret = __do_sys##name(__MAP(x,__SC_CAST,__VA_ARGS__));\
		__MAP(x,__SC_TEST,__VA_ARGS__);				\
		__PROTECT(x, ret,__MAP(x,__SC_ARGS,__VA_ARGS__));	\
		return ret;						\
	}								\
	__diag_pop();							\
	static inline long __do_sys##name(__MAP(x,__SC_DECL,__VA_ARGS__))
#endif /* __SYSCALL_DEFINEx */

長いマクロですね。reboot システムコールの宣言部分 __SYSCALL_DEFINE4(reboot, ...) は下記のように展開されます。

SYSCALL_DEFINE4(reboot, int, magic1, int, magic2, unsigned int, cmd,
		void __user *, arg)

↓

__diag_push();

__diag_ignore(GCC, 8, "-Wattribute-alias", "Type aliasing is used to sanitize syscall arguments");

asmlinkage long sys_reboot(__MAP(4,__SC_DECL,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)) __attribute__((alias(__stringify(__se_sys_reboot))));

ALLOW_ERROR_INJECTION(sys_reboot, ERRNO);

static inline long __do_sys_reboot(__MAP(4,__SC_DECL,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));

asmlinkage long __se_sys_reboot(__MAP(4,__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));

asmlinkage long __se_sys_reboot(__MAP(4,__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg))
{
	long ret = __do_sys_reboot(__MAP(4,__SC_CAST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));
	__MAP(4,__SC_TEST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg);
	__PROTECT(4, ret,__MAP(4,__SC_ARGS,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));
	return ret;
}

__diag_pop();

static inline long __do_sys_reboot(__MAP(4,__SC_DECL,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));

色々興味深い部分（__diag_push, __diag_ignore, などなど）はありますが、関数と関係ないので無視して、符号拡張を行う __se_sys_reboot() の実装部分を見ます。

asmlinkage long __se_sys_reboot(__MAP(4,__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg))
{
	long ret = __do_sys_reboot(__MAP(4,__SC_CAST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));
	__MAP(4,__SC_TEST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg);
	__PROTECT(4, ret,__MAP(4,__SC_ARGS,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));
	return ret;
}

この関数本体の部分がどう展開されるかを見ます。続きはまた今度。

RISC-V 64bit Linux とシステムコールと 32bit の引数 その 4 - __MAP(), __SC_X() マクロの実装

目次: RISC-V - まとめリンク

昨日の続きです。「64bit 環境においてシステムコールの引数がレジスタ長（= 64bit）以下だった場合（例えば int が典型例）、どう扱うか？」コードの解説編です。

どうして符号拡張されるのか？逆アセンブルで説明したつもりになってはいけない、逃げてはダメだ！という熱心なあなたのために、マクロを展開しながら説明します。今回の解説は __MAP(), __SC_LONG(), __SC_CAST() マクロの意味や実装です。

任意の引数に変換を掛ける __MAP() マクロ

鍵となるのは __MAP(4, ...) の展開後に出てくる __MAP4(m, t, a, ...) というマクロです。

// linux/include/linux/syscalls.h

/*
 * __MAP - apply a macro to syscall arguments
 * __MAP(n, m, t1, a1, t2, a2, ..., tn, an) will expand to
 *    m(t1, a1), m(t2, a2), ..., m(tn, an)
 * The first argument must be equal to the amount of type/name
 * pairs given.  Note that this list of pairs (i.e. the arguments
 * of __MAP starting at the third one) is in the same format as
 * for SYSCALL_DEFINE<n>/COMPAT_SYSCALL_DEFINE<n>
 */
#define __MAP0(m,...)
#define __MAP1(m,t,a,...) m(t,a)
#define __MAP2(m,t,a,...) m(t,a), __MAP1(m,__VA_ARGS__)
#define __MAP3(m,t,a,...) m(t,a), __MAP2(m,__VA_ARGS__)
#define __MAP4(m,t,a,...) m(t,a), __MAP3(m,__VA_ARGS__)
#define __MAP5(m,t,a,...) m(t,a), __MAP4(m,__VA_ARGS__)
#define __MAP6(m,t,a,...) m(t,a), __MAP5(m,__VA_ARGS__)
#define __MAP(n,...) __MAP##n(__VA_ARGS__)

マクロに渡された型 t と引数 a リストのペアを 1つずつ m(t, a) のように m で変換を掛けるマクロですが、イメージしづらいと思うので展開結果を載せます。

__MAP(4,__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)

↓

__MAP4(__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)

↓

__SC_LONG(int, magic1),
__SC_LONG(int, magic2),
__SC_LONG(unsigned int, cmd),
__SC_LONG(void __user *, arg)

実装はぱっと見では良くわかりませんが、展開結果はわかりやすいですね。

キャストの要、__SC_LONG(), __SC_CAST() マクロ

先頭の展開結果 __SC_LONG(int, magic1) に着目します。

// linux/include/linux/syscalls.h

#define __SC_LONG(t, a) __typeof(__builtin_choose_expr(__TYPE_IS_LL(t), 0LL, 0L)) a

#define __TYPE_AS(t, v)	__same_type((__force t)0, v)
#define __TYPE_IS_LL(t) (__TYPE_AS(t, 0LL) || __TYPE_AS(t, 0ULL))


// linux/include/linux/compiler_types.h

#define __same_type(a, b) __builtin_types_compatible_p(typeof(a), typeof(b))

# define __force

ここまでの実装を踏まえて展開すると下記のようになります。

__SC_LONG(int, magic1),

↓ __TYPE_IS_LL() を展開

__typeof(
  __builtin_choose_expr(
    (__same_type((int)0, 0LL) ||
     __same_type((int)0, 0ULL)), 0LL, 0L))
      magic1,

↓ __same_type を展開

__typeof(
  __builtin_choose_expr(
    (__builtin_types_compatible_p(typeof((int)0), typeof(0LL)) ||
     __builtin_types_compatible_p(typeof((int)0), typeof(0ULL))), 0LL, 0L))
       magic1,

最内側の __same_type() は引数に渡した値の型を比較するマクロです。マクロで使用する __builtin_types_compatible_p() は 1つ目と 2つ目の型が同じかどうか？を int で返す組み込み関数、typeof() は引数の型を返すキーワードです（関数ではないので値は返しません）。

• その他の組み込み関数のマニュアル（__builtin_types_compatible_p() も載っている）
• typeof のマニュアル

従って __same_type((int)0, 0LL) は int と 0LL の型、int 型と long long int 型が等しいか比較しています。当然、等しくありません。ちなみに int は __SC_LONG(int, magic1) の最初の引数 int から来たことを思い出していただけると嬉しいです。例えば __SC_LONG(long long, magic1) であれば比較結果は等しいと判断されるでしょう。

最内側の条件式 __same_type((int)0, 0LL) || __same_type((int)0, 0ULL)) の部分は型が long long もしくは unsigned long long かどうか？を調べています。

判断結果を受け取る __builtin_choose_expr(const_exp, exp1, exp2) は、const_exp が真なら、exp1 を返し、偽なら exp2 を返す組み込み関数です。もし const_exp がコンパイル時に判定できなければ、コンパイルエラーです。

以上をまとめると __SC_LONG(AAA, BBB) は、

• AAA が long long もしくは unsigned long long なら long long BBB
• AAA がそれ以外の型ならば long BBB

という働きをします。また __SC_CAST() はその名の通りキャストを生成します。

// linux/include/linux/syscalls.h

#define __SC_CAST(t, a)	(__force t) a

です。長いですね。続きはまた今度。

RISC-V 64bit Linux とシステムコールと 32bit の引数 その 5 - __MAP(), __SC_X() マクロの展開

目次: RISC-V - まとめリンク

昨日の続きです。「64bit 環境においてシステムコールの引数がレジスタ長（= 64bit）以下だった場合（例えば int が典型例）、どう扱うか？」コードの解説編です。

どうして符号拡張されるのか？逆アセンブルで説明したつもりになってはいけない、逃げてはダメだ！という熱心なあなたのために、マクロを展開しながら説明します。今回の解説は __MAP(), __SC_LONG(), __SC_CAST() マクロの展開結果です。

__MAP(), __SC_LONG(), __SC_CAST() マクロの合わせ技

今まで説明してきた __MAP(), __SC_LONG(), __SC_CAST() マクロを総動員して全部展開すると、

__MAP(4,__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)

↓

__SC_LONG(int, magic1),
__SC_LONG(int, magic2),
__SC_LONG(unsigned int, cmd),
__SC_LONG(void __user *, arg)

↓

long magic1,
long magic2,
long cmd,
long arg


__MAP(4,__SC_CAST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg)

↓

__SC_CAST(int, magic1),
__SC_CAST(int, magic2),
__SC_CAST(unsigned int, cmd),
__SC_CAST(void __user *, arg)

↓

(int)magic1,
(int)magic2,
(unsigned int)cmd,
(void __user *)arg

になります。他のマクロ __SC_TEST() と __SC_ARGS() と __PROTECT() は本筋ではないので詳細は省きますが、下記のような実装です。

// linux/include/linux/syscalls.h

#define __SC_TEST(t, a) (void)BUILD_BUG_ON_ZERO(!__TYPE_IS_LL(t) && sizeof(t) > sizeof(long))

#define __SC_ARGS(t, a)	a

#define __PROTECT(...) asmlinkage_protect(__VA_ARGS__)


// linux/include/linux/build_bug.h

/*
 * Force a compilation error if condition is true, but also produce a
 * result (of value 0 and type int), so the expression can be used
 * e.g. in a structure initializer (or where-ever else comma expressions
 * aren't permitted).
 */
#define BUILD_BUG_ON_ZERO(e) ((int)(sizeof(struct { int:(-!!(e)); })))


// linux/include/linux/linkage.h

#ifndef __ASSEMBLY__
#ifndef asmlinkage_protect
# define asmlinkage_protect(n, ret, args...)	do { } while (0)
#endif
#endif

展開するとこうなります。

__SC_TEST(AAA, BBB)

// AAA が long long でも unsigned long long でもない型で、
// long 型より大きい型（構造体とかが該当する）の場合に、
// コンパイルエラーにする。


__SC_ARG(AAA, BBB)

↓

BBB


__PROTECT(n, ret, args...)

↓

// m68k アーキテクチャでは引数をスタックに置かないようにする役目のコードになるらしい。
// 他のアーキテクチャでは下記のようになって何も意味はない。
do { } while (0)

これで不明な要素は全てなくなったはずです。続きはまた今度。

RISC-V 64bit Linux とシステムコールと 32bit の引数 その 6 - SYSCALL_DEFINE() を全部展開

目次: RISC-V - まとめリンク

昨日の続きです。「64bit 環境においてシステムコールの引数がレジスタ長（= 64bit）以下だった場合（例えば int が典型例）、どう扱うか？」コードの解説編です。

どうして符号拡張されるのか？逆アセンブルで説明したつもりになってはいけない、逃げてはダメだ！という熱心なあなたのために、マクロを展開しながら説明します。今回の解説は __MAP(), __SC_LONG(), __SC_CAST() マクロの意味や実装です。

SYSCALL_DEFINE() の展開

今までの SYSCALL_DEFINE(), __MAP(), __SC_LONG(), __SC_CAST() マクロを全部合体させて展開して、組み込み関数などを外すとこうなるはずです。

asmlinkage long __se_sys_reboot(__MAP(4,__SC_LONG,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg))
{
	long ret = __do_sys_reboot(__MAP(4,__SC_CAST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));
	__MAP(4,__SC_TEST,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg);
	__PROTECT(4, ret,__MAP(4,__SC_ARGS,int, magic1, int, magic2, unsigned int, cmd, void __user *, arg));
	return ret;
}

↓こうなる

asmlinkage long __se_sys_reboot(long magic1, long magic2, long cmd, long arg)
{
	long ret = __do_sys_reboot((int)magic1, (int)magic2, (unsigned int)cmd, (void __user *)arg));
	//__SC_TEST() による型のチェック、チェックに引っかかったらコンパイルエラー
	//__PROTECT(magic1, magic2, cmd, arg) に展開されるが RISC-V 64bit では何も意味はない
	return ret;
}

つまり __se_sys_reboot() ではシステムコールからの引数は long long, unsigned long long を除いて常に long 型の引数として扱い（__SC_LONG マクロの働き）、システムコールの本体 __do_sys_reboot() にはシステムコールが定義する本来の型にキャストして渡します（__SC_CAST マクロの働き）。

例えば musl libc の reboot() の magic1 にあたるレジスタ a0 には 0x00000000_fee1dead が渡されます。long の値として解釈すると正の数 4276215469 ですが、__do_sys_reboot() の呼び出し時に int にキャストして渡すため 0xfee1dead = -18751827 であると解釈されます（たぶん、大抵は）。

キャストされた負の int 型の値を __do_sys_reboot() は long 型の引数で受け取ります。このとき long 型で -18751827 と解釈される値（= 0xfee1dead を 64bit に符号拡張した 0xffffffff_fee1dead）をレジスタ a0 に格納して、関数呼び出しを行います。

(gdb) b __do_sys_reboot
Breakpoint 1 at 0xffffffff8002d860: file kernel/reboot.c, line 702.

(gdb) c
Continuing.

Breakpoint 1, __do_sys_reboot (magic1=-18751827, magic2=672274793,
    cmd=1126301404, arg=0x4321fedc) at kernel/reboot.c:702
702     {

(gdb) p/x $a0
$1 = 0xfffffffffee1dead

なので __do_sys_reboot() でブレークしてレジスタ a0 の内容を表示すると、0xffffffff_fee1dead になっていたわけです。

いやー長かった。マクロの魔術ですね。

C 言語で long を int にキャストしたらどうなるの？

目次: C 言語と libc - まとめリンク

最近 RISC-V 64bit 向けの Linux システムコールの実装を調べていました（2023年 2月 1日の日記など）。そのなかで long 型を int にキャストし、負の値として解釈している部分がありました。

C 言語に詳しい方は「あれ？」と思ったかもしれません。一見して問題なさそうな操作ですが、C 言語の仕様によれば結果は実装依存（implementation-defined）です。アーキテクチャやコンパイラ次第で結果が変わるということです。

条件は long より int の方が大きい（sizeof(long) > sizeof(int) になる）ときです。現代では 64bit 環境が該当するでしょう。32bit 環境は大抵 long と int が同じ値域なので問題は発生しません。

概要をお伝えしたところで、言語仕様を見てみましょう。

6.3.1.3 Signed and unsigned integers

When a value with integer type is converted to another integer type other than
_Bool, if the value can be represented by the new type, it is unchanged.

Otherwise, if the new type is unsigned, the value is converted by repeatedly
adding or subtracting one more than the maximum value that can be represented
in the new type until the value is in the range of the new type.

Otherwise, the new type is signed and the value cannot be represented in it;
either the result is implementation-defined or an implementation-defined signal
is raised.


（ざっくり和訳）

整数型の値が _Bool 以外の整数型に変換される場合、もし値が新しい型で表現できるなら
値は変化しません。

新しい型が符号なしで表現できないならば、新しい型の範囲に入るまで、新しい型が表現できる
最大値より 1 多い値を繰り返し加算または減算して値を変換します。

（訳注: パディング等も認められているのでややこしい言い方になっているが、16bit 整数で
最大値 = 0xffff だとすると、最大値 + 1 = 0x1_0000 となる。16bit 以上の上位ビットの
0 クリア処理に相当する）

新しい型が符号付きで値を表現できないならば、結果は実装定義になるか、もしくは実装定義の
シグナルが生成されます。

Linux は対象とするアーキテクチャとコンパイラが限定（GCC か Clang）されており、実装依存の動きを把握した上で使っていると思われます。承知の上というやつですね。

どの実装で使われるかわからないコードを書くなら、long から int へのキャストは常に結果が同じとは限らないので、キャストの結果で何が起きるか理解して使う必要があります。

気にしすぎても仕方がない

建前は上記の通りですが「全ての C コンパイラで動作するコード」を書く機会って、ほぼないと思います。基本的には細かい仕様を気にしすぎて何も進まなくなるより、時代で受け入れられる常識的な制約を付けて進めて、後で困ったら直せば良いんじゃない？と思います。

一方で 2038年問題のように、プログラムが書かれた時代は常識的な制約だったしみんな使っていたけど、何十年も経ってから大問題になるケースもあります。「常識的な制約 = ずっと無罪放免」かどうかは誰にもわかりません。

今の 64bit アーキテクチャに由来する時間や容量の制約は十分余裕があると思われますが、100年後の世界なんてわかりません。「誰が 64bit 決め打ちの実装にした！？100年前の奴らだと……？？」って泣きながら直している人がいても不思議じゃないですよね。

卵が割れる階数は？問題

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個目の最大の投擲回数の和を見たとき、上の階になるほど悪化します。例えば、

• 卵が割れる階数が 39F: 10F 〜 30F までで 3回、31F 〜 39F で最大 9回 = 12回
• 卵が割れる階数が 99F: 10F 〜 90F までで 9回、91F 〜 99F で最大 9回 = 18回

1個目の投擲方法を変えて最初は 10F 飛ばし、次は 9F 飛ばし、その次は 8F 飛ばし、のようにすると投擲回数の最大値が悪化しないことに気づくと思います。

こんな表で示すとわかりやすいでしょうか。

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個目を投げる階数を下記のようにしても良いです。

• 9F, 22F, 34F, 45F, 55F, 64F, 72F, 79F, 85F, 90F, 94F, 97F, 99F, 100F
• 10F, 23F, 35F, 48F, 56F, 65F, 73F, 80F, 86F, 91F, 95F, 98F, 100F

いずれも投擲回数は最大 14回となります。

CoreMark のコンパイルオプションをチューンする

目次: 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 のマニュアル）。

東京の道の名前

東京の道はようわからん通称（明治通り、みたいなやつ）が付いています。東京育ちの人にはなじみ深い名前だと思うんですが、都外出身者からすると、東京の道の通称と国道何号線、都道何号線の区分が全く合っていないので、知らない場所の道の通称を言われると結構困ります。

なぜかというと国道何号線、都道何号線という表記は地図で省かれることはない一方で、東京の道の通称は高速道路や地下鉄と重なったときに省かれてしまうことがあるためです。地図を見ても「〇〇通り？何それ？どこ？？」と困惑します。

通称とは一体？

東京の道の通称は東京都が決めています（東京都通称道路名〜道路のわかりやすく親しみやすい名称〜 - 東京都建設局）。道案内の青看板に載るような名前と思ってもらえばわかりやすいでしょう。

東京都がまとめてくれているのはありがたいんですが「東京都が公式に定めた通称」というのは不思議な響きで、それは公称ではなかろうか？通称とは一体……？？

なぜか存在しない大正通り

東京の道の通称には「年号」＋「通り」という名前がいくつかあります。このうちなぜか大正通りだけは存在しません。不思議ですね。

• 江戸通り（No.26: 国道6号など、千代田区大手町二丁目〜台東区花川戸二丁目）
• 明治通り（No.3: 都道416, 306号など、港区南麻布二丁目〜江東区夢の島）
• 昭和通り（No.24: 国道4号など、港区新橋一丁目〜台東区根岸五丁目）

調べてみると割と有名な話らしいです（「明治通り」「昭和通り」はあるのに、なぜか「大正通り」はない東京のちょっとした謎 - アーバン ライフ メトロ）。詳しくは記事を読んでいただくとして、簡単に言えば、

• 東京日日新聞の公募で大正通りと呼ぼうとした（今の靖国通り）が定着せず
• 東京都の公募で靖国通りが選ばれた、東京都建設局のまとめる一覧に大正通りは入ってない
• 東京に「大正通り」はある（武蔵野市吉祥寺）
• 「平成通り」もある（中央区兜町〜築地二丁目）
• 「令和通り」はまだない

ということみたいです。

SIMD を使ったお手軽最適化 - その 1

お手軽最適化のメモです。行列の掛け算を題材にします。

最初にお断りしておくと GEMM のような汎用処理の場合は、自分で最適化せずに OpenBLAS を使ってください（素人が最適化しても勝てません）。しかし OpenBLAS のような限界まで最適化されたコードは誰でも簡単には書ける、とは言えません。

スカラー処理だと遅いけれど、お手軽に最適化（数倍程度）がしたいときの参考になれば幸いです。

GEMM ってなんですか？

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列が必要です。式、および、視覚的に示すと下記のようになります。

行列の掛け算（式）

行列の掛け算（A の行と B の列の関係）

説明はこれくらいにしてコードを見ましょう。

基本コース - 素朴に演算

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倍という驚異のスピードです。すごいですね……。

行列の掛け算の説明でほぼ終わってしまいました。お手軽最適化はまた次に。

SIMD を使ったお手軽最適化 - その 2

お手軽最適化のメモ、昨日の続きです。行列の掛け算を題材にします。前回は行列の掛け算と素朴な実装のコードを紹介しました。今回はお手軽最適化を紹介します。

スカラー処理だと遅いけれど、お手軽に最適化（数倍程度）がしたいときの参考になれば幸いです。

お手軽コース - 自動ベクトル化

素朴版のコードですと i, j, k の順でループになっていて、k を最内ループにしていました。ループ内の計算は、

c[i * nn + j] += a[i * kk + k] * b[k * nn + j]

でした。このときメモリアクセスのパターンは、

• A: 行方向に連続（メモリアクセスパターンは連続 = 自動ベクトル化できる）
• B: 列方向に連続（メモリアクセスパターンは飛び飛びアクセス = 自動ベクトル化できない）
• C: 同一要素を何度もアクセス

素朴版の計算順（i, j, k ループ）

です。GCC の自動ベクトル化ですと、このアクセスパターンをうまく最適化できないようです。対象が最内ループのみなのかもしれません。ループを入れ替え j を最内にして、B と C を行方向に読むようにします。するとメモリアクセスのパターンは、

• A: 同一要素を何度もアクセス
• B: 行方向に連続（メモリアクセスパターンは連続 = 自動ベクトル化できる）
• C: 行方向に連続（メモリアクセスパターンは連続 = 自動ベクトル化できる）

ループ入れ替え版の計算順（i, k, j ループ）

です。ループを入れ替えるとループ内で 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 が計算偏重の処理で最適化の効果が出やすい、という点を考慮する必要はあるものの僅かな書き換えで得られる効果にしては割と良いのではないでしょうか。

もう少し頑張るコース - intrinsics

ソースコードを書き換える元気があれば、別の最適化方法もあります。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要素処理できるためです。

Intrinsics 版の計算順

SIMD 命令には Ai,k を SIMD レジスタの全要素に配る命令（set1_ps から生成される broadcast 命令）や、掛け算と足し算を一度に行う fmadd 命令など、この計算順に最適な命令が揃っています。

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 を使ってくださいませ。