---
title: 第12回 配信講義　計算科学技術特論A （2023）
tags:  #計算科学 #a64fx #カーネルのチューニング  
author: [R-CCS　計算科学研究推進室](https://convert.docswell.com/user/2300203199)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/4JQ8MNP67P.jpg?width=480
description: 第12回7月6日 A64FXでの重力N体計算カーネルのチューニング 重力多体問題、いわゆるN体計算のカーネル部分のA64FXプロセッサ向けチューニングについて紹介する。SVE命令を活用すべく、ACLEの組込関数を利用した。IntelアーキテクチャのAVX-512命令を用いた同等のコードとの性能比較も行う。
published: July 05, 23
canonical: https://convert.docswell.com/s/2300203199/K24NJG-2023-07-05-132648
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/4JQ8MNP67P.jpg)

A64FXでの重⼒N体計算
カーネルのチューニング
計算科学技術特論A（2023）
第12回 7⽉6⽇
理化学研究所
運⽤技術部⾨
技師
計算科学研究センター（ R-CCS）
ソフトウェア開発技術ユニット
似⿃ 啓吾（にたどり けいご）


# Page. 2

![Page Image](https://bcdn.docswell.com/page/K742VGYMJ1.jpg)

概要
• A64FXプロセッサ向けの重⼒N体計算カーネルのチューニング
⼿法
• SVE (Scalable Vector Extensions)命令拡張を
ACLE (Arm C Language Extensions)のから⽤いた
• Intel AVX-512版も作成
コンパイラSIMD化機能やCPUアーキを⽐較
• 性能
• ⾃動SIMD化
• ACLE版
• +asm編集
: 16.5 cycle, 51.6%
: 13.5 cycle, 63.0%
: 13.0 cycle, 65.2%
(*) 理論限界として8.5 cycleを仮定


# Page. 3

![Page Image](https://bcdn.docswell.com/page/LJ1M3D6YEG.jpg)

公開コード、バージョン情報
• https://github.com/nitadori/nbody-a64fxにて公開
• 2021年から少しだけ整理
• 単なる作業場です
• 使って何か研究してみたい⽅はご⼀報ください
• 再現実験は断りなくどうぞ
• バージョン情報
• 富⼠通⾔語環境： tcsds-1.2.36, FCC 4.8.0 (trad mode)
• インテルコンパイラ（OneAPI）： icpc (ICC) 2021.2.0


# Page. 4

![Page Image](https://bcdn.docswell.com/page/GJW3PYW1J2.jpg)

構成
1.
2.
3.
4.
5.
イントロ
N体計算カーネルと逆数平⽅根
Intel AVX-512でのチューニング
A64FX Arm SVEでのチューニング
まとめ


# Page. 5

![Page Image](https://bcdn.docswell.com/page/4EZVKX5XJ3.jpg)

背景・ねらい
• 「富岳」のA64FXプロセッサ：
• 命令セットのがSPARC（＋富⼠通拡張）→ARM（SVE拡張）へ
• アーキテクチャレジスタ数が128→32に減少
• 命令レイテンシが9と⼤きい
• 演算器を埋めきるスケジューリングが難しい
• これまで：レジスタ数を活かしてコンパイル時の静的スケジューリン
グ（ソフトウェアパイプライニング）
• これから：CPUのOut-of-Order実⾏（動的スケジューリング）に頼ら
ざるを得ない
• 経験知：緩いSWPLを施したバイナリのほうが性能が出る
• N体計算カーネル
• 複雑すぎず簡単すぎない、メモリネックにならない
• マイクロアーキテクチャのスケジューリング評価にはうってつけ


# Page. 6

![Page Image](https://bcdn.docswell.com/page/Y765Q49PEV.jpg)

レイテンシとスループット
• 最近のCPUの中⾝は流れ作業（パイプライン）
• レイテンシ（Cycle per Instruction, CPI）
• 演算結果が利⽤可能になるまで待つ時間
• ⽐喩：⾞を発注してから納品されるまでの⽇数
• スループット（CPIもしくはIPC, Instruction per Cycle）
• その命令が演算器を専有する時間（CPI）
• ⽐喩：⼀⽇に何台の⾞を⽣産できるか（IPC）
• IPC（⼤きいほど良い）のほうが直感的なこともあるけどレイ
テンシと単位を揃えてCPI（⼩さい⽅が良い）で書くことが多
い
• 例：A64FXのFMA命令ならL-T=9-0.5
（9-1のユニットが2つあるとも）


# Page. 7

![Page Image](https://bcdn.docswell.com/page/G75YRQNP74.jpg)

例：Intel© Intinsic guideより


# Page. 8

![Page Image](https://bcdn.docswell.com/page/9J2GDP5ZJR.jpg)

「理論限界」の考え⽅
• A64FXには2本のSIMD演算パイプラインがある
• FLAとFLB、⼀部の命令はFLAのみが実⾏可能
• 演算をしないシャッフル系の命令やストア命令もFLAを消費すること
に注意
• ロード命令や整数命令はこれらと同時発⾏できる
• この2本で同時にFMA (fused multiply-add)命令が出たときにピーク
のFLOPSが出る（⾏列乗算など）
• 今回のN体カーネルは浮動⼩数点演算命令が17
• 内9がFMA
• 逆数平⽅根の実装によっては+1、後述の再計算では+3命令
• 理想的にレイテンシが隠されたときの値として、8.5 cycleを基準に効
率を求めた


# Page. 9

![Page Image](https://bcdn.docswell.com/page/DEYV25KNEM.jpg)

ループアンロールと
ソフトウェアパイプライニング
• ループの持つ並列度を利⽤して演算器を飽和させる
• レジスタ本数に対するプレッシャーは増加
for(int i=0; i&lt;n; i++){
A(i);
B(i);
C(i);
}
基本のループ（擬似コード）
for(int i=0; i&lt;n; i+=4){
A(i);
B(i);
C(i);
A(i+1); B(i+1); C(i+1);
A(i+2); B(i+2); C(i+2);
A(i+3); B(i+3); C(i+3);
}
ループアンローリング
for(int i=0; i&lt;n;
A(i); A(i+1);
B(i); B(i+1);
C(i); C(i+1);
}
i+=4){
A(i+2); A(i+3);
B(i+2); B(i+3);
C(i+2); C(i+3);
アンロール＋スケジューリング
（ストライピング）
A(0);
A(1);
for(int
A(i+2);
A(i+3);
A(i+4);
}
A(n-1);
B(0);
i=0; i&lt;n-3; i+=3){
B(i+1); C(i);
B(i+2); C(i+1);
B(i+3); C(i+2);
B(n-2); C(n-3);
B(n-1); C(n-2);
C(n-1);
ソフトウェアパイプライニング


# Page. 10

![Page Image](https://bcdn.docswell.com/page/VJNM8NRV78.jpg)

N体計算の基礎⽅程式
&lt;latexit sha1_base64=&quot;6nPV07uEIx++bQdcSfomXsZGFL0=&quot;&gt;AAADSHichVHLahRBFL3dMSbGaCa6UHBTOkQicYbqUVSEQNCFriQPJwlMzzTVbc2kJtWPVNeMxk7/gH6AC1cKLsTPcCNuxUU+QdwoEYLgwjvd7TMkuUVXnbp1z+lzuW4kRawp3TbMoSPDR0dGj40dHz9xcqI0eWo5DnvK43UvlKFadVnMpQh4XQst+WqkOPNdyVfc9duD95U+V7EIg/t6M+JNn3UC0RYe05hySk9tydt6msySit1WzEvsiCktmEx/I2K7PlGOSElW2yCVO8SOe76TdGettHWP5ETf6SIn3lA6sftM8SgWMgxaNTJDtnKFLqn80tpq1VLUU6KzppvFeckplWmVZkH2AqsAZShiPix9ABseQAge9MAHDgFoxBIYxLga0EbMYQPzl8ECChG+NyHBd4VIZLUcUhhDnR4yOFYwzK7j3sFbo8gGeB/oxxnbwz9K/BQyCUzRj/Q13aHv6Bv6if7YVyvJNAa+NvF0cy6PnIknZ5d2D2X5eGpY+8M60LPGzm9kXgV6j7LMoAsv5/cfP9tZurk4lVykL+ln9P+CbtO32EHQ/+a9WuCLzw/wE+D+EPUe/eVp/2oXq1IcrPX/GPeC5VrVula9snC1PHerGPEonIMLMI2zuw5zcBfmoY7qX40zBjHOm+/NL+au+T0vNY2Ccxr+iSHzJxNc0EU=&lt;/latexit&gt;
=
&lt;latexit sha1_base64=&quot;OA7lYFsX8Xqt787tKzsvFwk+wWE=&quot;&gt;AAADWHichVHLbtNAFL2Ogbbm0VA2SGxGRC2pgGicIkBIlSpYwAr1lbZS3Vj2dJJM6lfHk0CZ+gf4ARZdgcQC8Rls2CMW5Q8QyyKxYcGNY1QgapmRZ+7jnOMzun4SiFRRemiUzDNnz42NT1jnL1y8NFm+PLWWxj3JeIPFQSw3fC/lgYh4QwkV8I1Eci/0A77u7zwa9Nf7XKYijlbVXsK3Qq8diZZgnsKSWz5wWtJj2gk91ZEh2W7Ws+NEYUYcPyTSFWSePCZO2gtd3Z23s+ZTkjMtHbrdzNJOwFuqioBdqbTT9yRPUhHEUbNObpL9oUaX3P6ttp8rS9HuqNmmnssyqzqCmXXLFVqj+SKjgV0EFSjWYlz+BA5sQwwMehAChwgUxgF4kOLehBbGHHaxfgtsoJBgfws09iVGIsdyyMBCnR4yOCI8rO7g2cZss6hGmA/005zN8I8BfhKZBKbpZ/qOHtGP9D39Sn+eqKVzjYGvPbz9IZcn7uTLqys//ssK8VbQOWad6lnhy+/nXgV6T/LK4BVsyO+/eHW08mB5Ws/QN/Qb+n9ND+kHfEHU/87eLvHlg1P8RHg+Q73nf3g6Ge0jKsPB2v+OcTRYq9fsu7W5pTuVhYfFiMfhGlyHKs7uHizAE1iEBjDDMG4Y1LBLX0wwx8yJIbRkFJwr8Ncyp34B/2TQWg==&lt;/latexit&gt;
2
d
r8 = ⌧
2
dC
#
’
9=1
⇣p
m
m r8
&lt;9
Y 2 + |r 9
• εはsoftenin⻑といって⻑さ
の次元を持つパラメタ
• すべてのN 粒⼦への加速度を
計算するにはN 2回の相互作⽤
計算が必要
• 天体の計算だけでなく、分⼦
動⼒学法なども同じような式
のかたちをしている
r8 |2
&quot;
⌧
#
’
9=1
⌘ 3 (r 9
p
&lt;9
Y 2 + |r 9
r8 )
r8 |2
#!


# Page. 11

![Page Image](https://bcdn.docswell.com/page/YE9NKRMVE3.jpg)

include &lt;cmath&gt;
リファレンスコード
• 最内ループのメモリアクセ
スはbody[j]のみ
struct Body{
float x, y, z, m;
};
struct Acceleration{
float ax, ay, az;
};
void nbody(
const int n,
const float eps2,
const Body body[],
Acceleration acc[])
{
for(int i=0; i&lt;n; i++){
float xi=body[i].x, yi=body[i].y, zi=body[i].z;
float ax=0, ay=0, az=0;
• 演算をカウントしてみる
•
•
•
•
•
for(int j=0;
float dx
float dy
float dz
3-sub
3-mad
1-rsqrt
3-mul
3-mad
j&lt;n; j++){
= body[j].x - xi;
= body[j].y - yi;
= body[j].z - zi;
float r2 = eps2 + dx*dx;
r2 += dy*dy;
r2 += dz*dz;
float ri = 1.f / sqrtf(r2);
float mri = body[j].m * ri;
float ri2 = ri * ri;
float mri3 = mri * ri2;
• madは積和演算（1命令で2
演算できる）
• rsqrtは逆数平⽅根
ax += mri3 * dx;
ay += mri3 * dy;
az += mri3 * dz;
}
acc[i] = {ax, ay, az};
}
}
return;


# Page. 12

![Page Image](https://bcdn.docswell.com/page/GE8KRWL47D.jpg)

逆数平⽅根の重要性
• 1/ 𝑥の計算にはnaïveには除算と平⽅根が1回ずつ必要になる
• 除算も平⽅根も加減乗算に⽐べて桁で遅い
• 原理は⼈間の筆算や開平法と同じ
• 加減乗算と積和算は12命令しかないのでこのままでは確実にボトルネ
ック！
• 最近のプロセッサでのアプローチ
• 逆数1/𝑥と逆数平⽅根 1/ 𝑥の近似命令をサポート
（SIMDで1〜2サイクルスループット）
• 加減乗算の組み合わせで有効桁数を倍々に増やせる
• 初期値推定の有効ビット数はアーキテクチャごとにマチマチ
• SSE/AVX: 12-bit、AVX-512F: 14-bit、AVX-512ER (KNL): 28-bit、
Fujitsu (Sparc/Arm): 8-bit
• 収束補助の専⽤命令を備えるものも（3DNOW!、Arm SVE）


# Page. 13

![Page Image](https://bcdn.docswell.com/page/LELZKNLVJR.jpg)

収束公式の導出
H 0 :&#039; (1 + Y)G 1/=
() (⌘ := 1 GH 0 = ) = 1
&lt;latexit sha1_base64=&quot;F5DMtzbsYhv7gxiFzIVFMUqkK4g=&quot;&gt;AAAEC3ichVG9b9NAFH+u+Sjho2lZkFhORI0SlYRzqABFQqpgYUHqB2kr1U1ku5f4VPvs2peQYOUfQGJmYAKJATHChNgQEjtiKCssiLFILAw8OyZt+sVZvnv33vv93u/dM32Hh5LSbWVMPXHy1OnxM5mz585fmMhOTi2HXjuwWM3yHC9YNY2QOVywmuTSYat+wAzXdNiKuXk3jq90WBByTzyQPZ+tu0ZL8Ca3DImuxqRyP096DUqqeshdtkUK2ozeMQLmh9zxRJF061FJuyb6RNczOm82SZ7oDmvKArFJ9TbRSIl0SYQM/bogesBbtiwS9Jf2EcXRIcNICLMLWsku1qND6/TiMrsqhhVSzMAdd4DAuBUtztddQ9qBG0mjhw/QTxWPQP4xDaG67duGkJ6L3WgYH7xHCtXIDLGQ207OCrHrFbR0a8OT4V6qRjZHyzRZ5KChpUYO0jXvZT+DDhvggQVtcIGBAIm2AwaE+K1BE20GW+i/ChpQ8DG+DhHGA7R4ksugDxnkaSOCYYaB3k3cW3hbS70C7zF/mKAtrOjgHyCSwDT9Ql/RHfqJvqY/6J8juaKEI9bVw9McYJnfmHh8aen3f1EunhLsXdSxmiV2fivRylG7n3jiLqwBvvPo6c5SdXE6ytMX9Cfqf0636QfsQHR+WS8X2OKzY/QI3B8iX3ePpqOzTczq42C1/WM8aCxXytqN8vWF2dzcnXTE43AZrkABZ3cT5uAezEMNLOWt8lX5pnxXn6hv1Hfq+0HqmJJiLsLIUj/+BYj8Bq8=&lt;/latexit&gt;
(1 + Y) =
() 1 + Y = (1 ⌘) 1/=
⇣
⌘
1/=
() H := G
= (1 ⌘) 1/= H 0
⇣
⌘
H 1 := taylor (1 ⌘) 1/= H 0
&#039; 1 + 2 1 ⌘ + 2 2 ⌘2 + · · · H 0
⼀般にn乗根の逆数に対して公式が存在する
1. ⼩さな値hを求める
（正しく計算できていれば0）
2. 多項式𝑝 ℎ を計算
3. 近似値に𝑝 ℎ を掛ける
k次の多項式なら有効桁数は 𝑘 + 1 倍に


# Page. 14

![Page Image](https://bcdn.docswell.com/page/4JMV2XMNJW.jpg)

具体例
• 逆数
&lt;latexit sha1_base64=&quot;wiBiQREPOGPZUcJDJhsni7Jlx/c=&quot;&gt;AAADG3ichVHLahRREK20rzg+MupGcHNxSIhoxupEVAJC0I3LPJwkkE6a7s6d6Uv6le4744xDfsAfcCEuFFyInxEXIu4ki3yCupyAGxee7mnwERLr0nWrTtWpPpdyk0Blmnl/xDhx8tTpM6NnK+fOX7g4Vr10eTmL26knG14cxOmq62QyUJFsaKUDuZqk0gndQK64W4/y+kpHppmKoye6l8j10GlFqqk8RwOyqxs9m8XErJWpUG4L83ZXWFbFB/JAmGJKdEVeB9SzzQLM05tDMJBNPSl8pP7GdOFn4C1vM9aZsFLV8vUNu1rjOhcmDgdmGdSotPm4+oUs2qSYPGpTSJIi0ogDcijDWaMmYknbwG+RSUwJ6uvURz1FpIpeSTtUwZw2GBIdDtAt+BaytRKNkOfzs4Lt4Y8BvhRMQeO8x+94wB/5PX/ln0fO6hczcl093O6QKxN77PnVpR//ZYW4Nfm/Wcdq1nj5/UKrgvakQPJXeEN+59mLwdLs4nh/gt/wd+h/zfu8ixdEnQPv7YJcfHmMngj+KeZ1/9B0dLeLrh0s1vx3jYeD5em6ebc+s3CnNvewXPEoXaPrNInd3aM5ekzz1MD0D/SNBnRgvDJ2jU/G52GrMVJyrtBfZuz9Aj5kuf0=&lt;/latexit&gt;
参考： 逆数と平⽅根を求める⾼次収束アルゴリズム
http://www.finetune.co.jp/~lyuka/technote/fract/sqrt-jp.html
H 0 :&#039; 1/G
⌘ := 1
GH 0
H 1 := H 0 + H 0 ⌘ + ⌘2 + ⌘3 + · · ·
• 逆数平⽅根
&lt;latexit sha1_base64=&quot;zto1s2pYx8J+JNmHuC//zhtWe8Q=&quot;&gt;AAADQHichVHLahRBFL3d8ZGMjxnjJsFN4ZAQUceqGU1CQAi6cZmHkwTSsemu1EwX6ddU14yZNLNy5w+4cKXgQvwL3bgXF/kDgwuRCIK48HZPg4+QWEVXnXvuPbdPcd3Yl4mmdN8wR06dPnN2dKx07vyFi+XKpfG1JOoqLpo88iO14TqJ8GUomlpqX2zESjiB64t1d+d+ll/vCZXIKHyo+7HYCpx2KFuSOxopu/Kkb1MyvWAlMhAdwm5ZSUdpskssq+Qhf5cwchPDFMsGj+oZ3bdZnsiE1/PT8kVLzxCrpRzO6sRDOseNeeKhpojupGx2gEQjI/h2pBNiKdn29DW7UqU1mi9yFLACVKFYS1HlA1iwDRFw6EIAAkLQiH1wIMG9CS3EAjrI3wAGFGLMb0GKeYVI5rUCBlDCPl1UCKxwkN3Bs43RZsGGGGf9k1zN8Y8+fgqVBKboR/qaHtL39A09oD+P7ZXmPTJffbzdoVbEdvnpxOr3/6oCvDV4v1Unetb48vncq0Tvcc5kr+BDfW/v2eHqwspUOk1f0s/o/wXdp+/wBWHvG3+1LFaen+AnxPMx9tv9w9Px1S5WDXCw7N8xHgVr9RqbrTWWb1cX7xUjHoUrcBVmcHZzsAgPYAma2P3AKBsTxqT51vxkfjG/DktNo9Bchr+W+eMXjhfE6g==&lt;/latexit&gt;
p
H 0 :&#039; 1/ G
GH 0 2
✓
◆
1
3 2 5 3
H 1 := H 0 + H 0 ⌘ + ⌘ + ⌘ + · · ·
2
8
16
⌘ := 1


# Page. 15

![Page Image](https://bcdn.docswell.com/page/PJR2MNV4E9.jpg)

Break：
• 係数を求める程度ならブラ
ウザと⽇本語でOK
• 昔ならMaxima、最近なら
SymPyが便利かも
• Chat-GPTはWolframのプラ
グインがあればこのページ
と同じになる


# Page. 16

![Page Image](https://bcdn.docswell.com/page/PEXDVNZZ7X.jpg)

N体計算に有効な⽅法
• 逆数平⽅根を補正してから3乗するよりは、3乗してから補正す
るほうが少し有利
• 演算数が1つ減る（𝑟 !&quot; に再利⽤性がある）
• 直列な命令依存が減り命令レベル並列度が上がる（7段→5段）
• 補正多項式は 1 − ℎ !#/&quot; のテイラー展開
Ã
1
p
:&#039; 1/ A 2
Ã
2
:= Ã
&lt;latexit sha1_base64=&quot;FH2Elrj/31zpvm/zSUmyvKF3C4g=&quot;&gt;AAAESXichVHLbtNAFL2xC5TwaFo2ldgMjVohtQ3jtEAVCamCDcs+SFupbiPbmTSmfmU8CRTLPwAfwIIVSCwQn8EGiSVi0U9ArKBIIMSC60faxGnKWJ65c+ac43M9umeZvqD0MCfJI+fOXxi9mL90+crVscL4xIbvtrnBqoZruXxL13xmmQ6rClNYbMvjTLN1i23q+w+i880O477pOo/Egcd2bG3PMRumoQmEauO5FzNEFaZVZ4TvBvNKSCqqb9qsRZRbqt/iIuC75ZCoar6Ph1DlXkaIO5v5WTASNiOyQuYRLJ9C6/p7iWeDa4ZSJrMoO+ZGGAkWwqU5orbaWj3zlebQNM3YudtaljQ74JPKPHKsS1rtd0123Sx2j79de9xPSowSykJMOeH3EON/UCsUaYnGgwwWSloUIR0rbuEzqFAHFwxogw0MHBBYW6CBj882NLBm0EJ8DhSg4OH5DgR4zrEyYy6DEPLo00YFQ4aG6D7Oe7jbTlEH95G/H6sN/KKFL0clgWn6hb6jR/QjfU+/0r9DvYLYI8p1gKueaJlXG3s+uf7rvyobVwHNE9WZmQV2vhRnNTG7FyNRF0ai7zx7ebReWZsOZugb+g3zv6aH9AN24HR+Gm9X2dqrM/I4OD9Bv6c9mYazdWSFeLFK9hoHi41ySblTur26WFy+n17xKFyHKbiJd3cXluEhrEAVjNwPaVK6IU3Jn+Tv8m/5T0KVcqnmGvSNEfkf8Cka+g==&lt;/latexit&gt;
⌘ := 1
? :=
1
2
1
A 2 ⇥ Ã
+ ⌘ ⇥ 38 ,
A
1
:= Ã
1
A
2
:= A
1
&lt;A
3
⇥ Ã
+ Ã
⇥A
Ã
1
p
:&#039; 1/ A 2
1
Ã
2
:= Ã
2
⌘ := 1
&lt;latexit sha1_base64=&quot;ezYgX7H43PBpF7Q+6HNh4qu9tEo=&quot;&gt;AAAEtnichVG7b9NAGP/iGijh0bQsSCwnolZFbcPZLRBFqlTBwtgHaSvVTWQ7l9jUr9iXQLEysjDAyMAEEgPiz2BhRwz9ExBjkWBg4POjSeOk5Sz7vvvd93ucT/MsM+CUHuWECfHCxUuTl/NXrl67PlWYntkO3I6vs6ruWq6/q6kBs0yHVbnJLbbr+Uy1NYvtaAePov2dLvMD03We8EOP7dtqyzGbpq5yhOrTuT9zROGm1WDEr4VLUo9UlMC0WZtId5Wg7fPQr8k9oij5oT6EKqsZIq5sFgyDi0Rpd9QGsbMmqwk0QkicjKhBIksIymPa5IFu/enpjeVYecRtnEJiZGfYsfGoaCpgnNh6yembvqovy2QBaX2LCCOhdK9XPnEYBMu6LYzxT4W89JfHdorFmnye2GkTCmXiJfsSCvZDGf1FEsaoyYpvtgx+J53qhSIt0XiQ0UJKiyKkY90tfAMFGuCCDh2wgYEDHGsLVAjw2YMm1gzaiC+CBBQ83N+HEPd9rMy4l0EP8qjTQQbDDhXRA/y2cLWXog6uI/0gZuvoaOHrI5PALP1OP9Fj+pV+pj/o3zO1wlgjynWIs5ZwmVefenVz6/d/WTbOHIwB69zMHE9ejrOamN2LkegUesLvvnh7vFXZnA3n6Af6E/O/p0f0C57A6f7SP26wzXfn5HHw+wz1np/KdHa3hl09vFgpe42jxbZcku6XVjZWimsP0yuehFtwG+bx7h7AGjyGdaiCLtSFl8Jr4Y1YFmsiE1tJq5BLOTdgaIjeP7XcP1Q=&lt;/latexit&gt;
Ã
1
1
1
⌘ := Ã
⌘⇥?
,
&lt;A
1
1
⇥⌘
:= &lt; 9 ⇥ A
:= &lt;A 1 ⇥ A 2
7 stage, 9 ops
1
&lt;Ã
3
&lt;A
3
1
⇥ Ã
1
2
,
⌘ := &lt; 9 Ã
3
⇥ ⌘,
A 2 ⇥ Ã
:= &lt;Ã
&lt;A
3
3
,
&lt;Ã
&lt; 9 Ã
+ &lt;Ã
= &lt;Ã
3
1
3
:= &lt; ⇥ Ã
3
:= &lt;Ã
? :=
⌘⇥?
1 + 32 ⌘ +
3
2
1
⇥ Ã
+⌘⇥
15 2
8 ⌘
5 stage, 8 ops
1
15
8
2


# Page. 17

![Page Image](https://bcdn.docswell.com/page/3EKXQN8V7D.jpg)

外側ループと内側ループの並列度
• 外側ループ
• i-loop、i 並列
• ⼒を受ける粒⼦についての並
列度
• 可能ならこちらでのSIMD化
が望ましい
• j 粒⼦を1つ読んでくると
SIMD幅（16）の分再利⽤性
がある
void nbody(
const int n,
const float eps2,
const Body body[],
Acceleration acc[])
{
for(int i=0; i&lt;n; i++){
float xi=body[i].x, yi=body[i].y, zi=body[i].z;
float ax=0, ay=0, az=0;
for(int j=0;
float dx
float dy
float dz
j&lt;n; j++){
= body[j].x - xi;
= body[j].y - yi;
= body[j].z - zi;
float r2 = eps2 + dx*dx;
r2 += dy*dy;
r2 += dz*dz;
• 内側ループ
float ri = 1.f / sqrtf(r2);
• j-loop、j 並列
• ⼒を及ぼす粒⼦についての並
列度
• こちらでSIMD化すると、
ループの最後に縮約（総和）
演算が必要
• AoSからSoAへ変更しないと
性能を出しづらい
float mri = body[j].m * ri;
float ri2 = ri * ri;
float mri3 = mri * ri2;
ax += mri3 * dx;
ay += mri3 * dy;
az += mri3 * dz;
}
acc[i] = {ax, ay, az};
}
}
return;


# Page. 18

![Page Image](https://bcdn.docswell.com/page/L735PVMQE5.jpg)

コンパイラの対応状況
• Intel
• ディレクティブを付けることで外側（i ）ループのSIMD化にも対応
• 富⼠通
• 最内（j ）ループのみ⾃動的にSIMD化
• 最後の⽔平総和も出⼒してくれる
• 今回は双⽅のアーキテクチャに向けて外側ループの⼿動SIMD
化を⾏った


# Page. 19

![Page Image](https://bcdn.docswell.com/page/87DDL86W7G.jpg)

Intel編
• oneAPIとして今年から無償提供されるようになったコンパイ
ラを利⽤
• icpc -fast -qopt-zmm-usage=high -qopenmp nbody.cpp –S
icpx -fopenmp nbody.s
• デフォルトで256-bitのymmレジスタコードを出すので512-bitのzmm
レジスタを使うように
• OpenMPはHyper-Threadingのテスト⽤
• CascadeLake Xeonのクロックを2.0 GHzに固定して計測
• $ sudo cpupower frequency-set --min 2.0GHz --max 2.0GHz
• 外側もしくは内側もループにSIMD化ディレクティブを付けた
• #pragma omp simd
をループの前に付けるだけ


# Page. 20

![Page Image](https://bcdn.docswell.com/page/VJPX28YX78.jpg)

ソースコードの微変更
Before
After
#pragma omp simd
for(int i=0; i&lt;n; i++){
const float xi=body[i].x, yi=body[i].y,
zi=body[i].z;
float ax=0, ay=0, az=0;
for(int j=0; j&lt;n; j++){
float dx = body[j].x - xi;
float dy = body[j].y - yi;
float dz = body[j].z - zi;
...
ax += mri3 * dx;
ay += mri3 * dy;
az += mri3 * dz;
}
acc[i] = {ax, ay, az};
}
#pragma omp simd
for(int i=0; i&lt;n; i++){
const float xi=body[i].x, yi=body[i].y,
zi=body[i].z;
float ax=0, ay=0, az=0;
for(int j=0; j&lt;n; j++){
float dx = xi - body[j].x;
float dy = yi - body[j].y;
float dz = zi - body[j].z;
...
ax -= mri3 * dx;
ay -= mri3 * dy;
az -= mri3 * dz;
}
acc[i] = {ax, ay, az};
}
broadcastss (%rsi,%r14), %zmm1
vbroadcastss 4(%rsi,%r14), %zmm2
vbroadcastss 8(%rsi,%r14), %zmm3
vsubps
%zmm28, %zmm1, %zmm18
vsubps
%zmm25, %zmm2, %zmm17
vsubps
%zmm26, %zmm3, %zmm16
vsubps
vsubps
vsubps
(%rsi,%r14){1to16}, %zmm28, %zmm19
4(%rsi,%r14){1to16}, %zmm25, %zmm20
8(%rsi,%r14){1to16}, %zmm26, %zmm22
放送ロードがメモリオペランドになっている


# Page. 21

![Page Image](https://bcdn.docswell.com/page/2EVLDNG3JQ.jpg)

改善できそうな点
..B2.6:
# Preds ..B2.6 ..B2.5
# Execution count [1.25e+01]
vbroadcastss (%rsi,%r14), %zmm1
incq
%r12
vbroadcastss 4(%rsi,%r14), %zmm2
vbroadcastss 8(%rsi,%r14), %zmm3
vmovaps
%zmm20, %zmm0
vsubps
%zmm28, %zmm1, %zmm18
vsubps
%zmm25, %zmm2, %zmm17
vsubps
%zmm26, %zmm3, %zmm16
vfmadd231ps %zmm18, %zmm18, %zmm0
vfmadd231ps %zmm17, %zmm17, %zmm0
vfmadd231ps %zmm16, %zmm16, %zmm0
..___tag_value__Z11nbody_ipar0ifPK4BodyP12Acceleration.77:
call
*__svml_invsqrtf16_z0@GOTPCREL(%rip)
..___tag_value__Z11nbody_ipar0ifPK4BodyP12Acceleration.78:
vmulps
12(%rsi,%r14){1to16}, %zmm0, %zmm4
addq
$16, %r14
vmulps
%zmm0, %zmm0, %zmm5
vmulps
%zmm5, %zmm4, %zmm6
vfmadd231ps %zmm18, %zmm6, %zmm24
vfmadd231ps %zmm17, %zmm6, %zmm23
vfmadd231ps %zmm16, %zmm6, %zmm27
cmpq
%r15, %r12
jb
..B2.6
# Prob 82%
#24.15
#23.3
#25.15
#26.15
#28.25
#24.27
#25.27
#26.27
#28.25
#29.4
#30.4
#32.21
#34.28
#23.3
#35.21
#37.23
#39.4
#40.4
#41.4
#23.3
#23.3
• 逆数平⽅根が
SVML callにな
っている
• Short Vector
Math Library
• KNLのときは
vrsqrt28ps⼀発
だった


# Page. 22

![Page Image](https://bcdn.docswell.com/page/57GXYK8YEL.jpg)

Intrinsic版
{
const __m512 eps2 = _mm512_set1_ps(eps2_ss);
for(int i=0; i&lt;n; i+=16){
__m512 xi, yi, zi, mi;
transpose_4AoStoSoA(body+i, xi, yi, zi, mi);
__m512 ax, ay, az;
ax = ay = az = _mm512_set1_ps(0);
• ⼿動でSIMD化したもの
• rsqrtCubedは7命令、計
16命令
• GCCやClangもこのコー
ドを受け付ける（はず）
for(int j=0; j&lt;n; j++){
__m512 xj = _mm512_set1_ps(body[j].x);
__m512 yj = _mm512_set1_ps(body[j].y);
__m512 zj = _mm512_set1_ps(body[j].z);
__m512 mj = _mm512_set1_ps(body[j].m);
__m512 dx = _mm512_sub_ps(xi, xj);
__m512 dy = _mm512_sub_ps(yi, yj);
__m512 dz = _mm512_sub_ps(zi, zj);
• #include &lt;x86intrin.h&gt;
が必要
__m512 r2 = _mm512_fmadd_ps(dx, dx, eps2);
r2 = _mm512_fmadd_ps(dy, dy, r2);
r2 = _mm512_fmadd_ps(dz, dz, r2);
__m512 mri3 = rsqrtCubed(r2, mj);
ax = _mm512_fnmadd_ps(mri3, dx, ax);
ay = _mm512_fnmadd_ps(mri3, dy, ay);
az = _mm512_fnmadd_ps(mri3, dz, az);
}
transpose_3SoAtoAoS(ax, ay, az, acc+i);
}
}


# Page. 23

![Page Image](https://bcdn.docswell.com/page/4EQ8MN96JP.jpg)

性能値
• 内側ループSIMD化ではgathe命令が出ているが性能は良好
• 理論限界8.5 cycleに対して⼿動版は71%の効率
• Gflops値は相互作⽤数あたり「歴史的な係数38」を掛けたもの
、実際の演算数は24程度（ FMAを２演算それ以外を１演算）
cycle Gflops
Remark
i-par
16.9
71.9 外側ループSIMD
j-par
16.9
71.9 内側ループSIMD
intrin
11.9 102.6 ⼿動SIMD（外側）


# Page. 24

![Page Image](https://bcdn.docswell.com/page/KJ42VGXME1.jpg)

Hyper-Threading
• 1コアで2スレッド⾛らせることでスケジューリングを改善でき
る可能性がある
• 2スレッドがひとつのコアに当たるように以下の環境変数に
export KMP_AFFINITY=granularity=fine,compact
export OMP_NUM_THREADS=2
• 1割弱の改善はあったのだが、2スレッド実⾏環境でシングルス
レッド実⾏が遅くなってしまった
（裏で待機しているスレッドのせい？）
• 実アプリでは⾜を引っ張りそう
cycle(nt=1) cycle(nt=2)
ipar
16.9
21.4
jpar
16.9
21.4
intrin
11.9
15.0
intrin+omp
11.9
11.0


# Page. 25

![Page Image](https://bcdn.docswell.com/page/LE1M3DNY7G.jpg)

構造体アクセス命令とレジスタ内転置
• 外側ループであるが、
{x,y,z,m}[16] -&gt; {x[16],y[16],z[16],m[16]}
{ax[16],ay[16],az[16]} –&gt; {ax,ay,az}[16]
のような転置がある
• SVEにはld4、st3のようなstructure load/store命令がある
• AVX-512にそういうのはないのでレジスタ内で転置する
• 今回はボトルネックではないけどgather/scatteするのも悔しいので
• ⼀度作っておくと使い回しが効く
• 命令セットの格好の練習問題
• MDなど粒⼦をリストアクセスするケース、body[j]ではなく
body[idx[j]]が必要になるケースには有⽤


# Page. 26

![Page Image](https://bcdn.docswell.com/page/GEW3PY41E2.jpg)

Transpose_4AoStoSoA
• 2つのレジスタを連結さ
せて任意の要素を取り
出してくるという強⼒
な命令がある
c[i] = (a~b)[idx[i]]
• この命令でちょとした
テーブル参照もできる
r0
x0
y0
z0
m0
・・・
x3
y3
z3
m3
r1
x4
y4
z4
m4
・・・
x7
y7
z7
m7
r2
x8
x8
z8
m8
・・・
x11 y11 z11 m11
r3
x12 y12 z12 m12
・・・
x15 y11 z11 m15
xzlo x0
x4
z0
z4
・・・
x3
x7
z3
z7
unpcklo(r0,r1)
ymlo y0
y4
m0
m4
・・・
y3
y7
m3
m3
unpckhi(r0,r1)
xzhi x8
x12 z8
z12
・・・
x11 x11 z11 z11
unpcklo(r2,r3)
ymhi x8
x8
m8
m12
y11 y11 m11 m11
unpckhi(r2,r3)
xv
x0
x1
x2
x3
・・・
x12 x13 x14 x15
vpermt2ps(xzlo,xzhi)
yv
y0
y1
y2
y3
・・・
y12 y13 y14 y15
vpermt2ps(ymlo,ymhi)
zv
z0
z1
z2
z3
・・・
z12 z13 z14 z15
vpermt2ps(xzlo,xzhi)
mv
m0
m1
m2
m3
・・・
m12 m13 m14 m15
vpermt2ps(ymlo,ymhi)


# Page. 27

![Page Image](https://bcdn.docswell.com/page/47ZVKX2XE3.jpg)

Transpose_3SoAtoAoS
• 128-bit単位のシャッ
フル命令を使⽤
ax
x[0-3]
x[4-7]
x[8-11]
x[12-15]
ay
y[0-3]
x[4-7]
y[8-11]
y[12-15]
az
z[0-3]
z[4-7]
z[8-11]
z[12-15]
xylow
x[0-3]
x[4-7]
y[0-3]
x[4-7]
vshuff32x4(ax,ay)
xymid
x[4-7]
x[8-11]
x[4-7]
y[8-11]
vshuff32x4(ax,ay)
xyhig
x[8-11]
x[12-15]
y[8-11]
y[12-15]
vshuff32x4(ax,ay)
az
z0
z1
z2
a0
x0
y0
z0
a1
y5
z5
a2
z10
z3
・・・
・・・
・・・
・・・
z12 z13 z14 z15
x5
vpermt2ps(xlow,az)
x10 y10
vpermt2ps(xmid,az)
x15 y15 z15
vpermt2ps(xhig,az)


# Page. 28

![Page Image](https://bcdn.docswell.com/page/YJ65Q4VP7V.jpg)

Fujitsu ARM編
• とりあえずコンパイラの⾃動SIMD化から
• 性能が出るのはtradモード、clangモードはまだまだこれから
• 内側ループを⾃動SIMD化、gather命令が出てしまう
• Structure load命令は出してくれないようだ
（構造体をfloat[n][4]の配列に書き換えても出ない）
• SoA（x[n], y[n], z[n], m[n]）に並び替えた版も作ってみた
Gather
SoA
Cycle(noswp) Cycle(swp)
39.7
32.0
25.7
17.9
(2021年版）
• ソフトウェアパイプライニングが効いている
• コンパイルオプションは単に-Kfastもしくは-Kfast,noswp
• 18命令のループなのでコンパイラ任せで理論限界の半分


# Page. 29

![Page Image](https://bcdn.docswell.com/page/GJ5YRQDPJ4.jpg)

最新の値
• Gather版が少し速くなってた（コンパイラ改善？）
• ACLE版で効果のあった⼿法（unroll、recalc）もバックポート
してみた
• Gather
• SoA
• Unroll
• Recalc
: 28.8 cycles, 29.5%
: 17.9 cycles, 47.4%
: 17.1 cycles, 49.5%
: 16.5 cycles, 51.5%


# Page. 30

![Page Image](https://bcdn.docswell.com/page/9E2GDP8ZER.jpg)

ACLE版
void nbody_sve(
const int n,
const float eps2_ss,
const Body body[],
Acceleration acc[])
{
const svfloat32_t eps2 = svdup_f32(eps2_ss);
const svfloat32_t one
const svfloat32_t a
const svfloat32_t b
• Structure load/store命令を
外側ループに使⽤
• 外側ループを更にハンドア
ンロールしたバージョンも
作成した（ni32）
= svdup_f32(1.0);
= svdup_f32( 3./2.);
= svdup_f32(15./8.);
const svbool_t p0 = svptrue_b32();
for(int i=0; i&lt;n; i+=16){
svfloat32x4_t ibody = svld4_f32(p0, (const float *)(body+i));
svfloat32_t xi = svget4_f32(ibody, 0);
svfloat32_t yi = svget4_f32(ibody, 1);
svfloat32_t zi = svget4_f32(ibody, 2);
svfloat32_t ax, ay, az;
ax = ay = az = svdup_f32(0);
• ソフトウェアパイプライニン
グはそれでも効いた
for(int j=0; j&lt;n; j++){
svfloat32_t xj = svdup_f32(body[j].x);
svfloat32_t yj = svdup_f32(body[j].y);
svfloat32_t zj = svdup_f32(body[j].z);
svfloat32_t mj = svdup_f32(body[j].m);
cycle eff.
Gflops
ACLE 16.1 52.8
75.5
ni32
14.8 57.6
82.3
svfloat32_t dx = svsub_f32_x(p0, xj, xi);
svfloat32_t dy = svsub_f32_x(p0, yj, yi);
svfloat32_t dz = svsub_f32_x(p0, zj, zi);
svfloat32_t r2 = svmad_f32_x(p0, dx, dx, eps2);
r2 = svmad_f32_x(p0, dy, dy, r2);
r2 = svmad_f32_x(p0, dz, dz, r2);
理論限界を8.5 cycle/loop
相互作⽤あたり38 flopと仮定
svfloat32_t mri3 = rsqrtCubed(r2, mj, p0, one, a, b);
ax = svmla_f32_x(p0, ax, mri3, dx);
ay = svmla_f32_x(p0, ay, mri3, dy);
az = svmla_f32_x(p0, az, mri3, dz); // 17-ops
}
}
}
svfloat32x3_t acci = svcreate3_f32(ax, ay, az);
svst3_f32(p0, (float *)(acc+i), acci);


# Page. 31

![Page Image](https://bcdn.docswell.com/page/D7YV256N7M.jpg)

スケジューリングの
為の再計算
for(int j=0; j&lt;n; j++){
svfloat32_t xj = svdup_f32(body[j].x);
svfloat32_t yj = svdup_f32(body[j].y);
svfloat32_t zj = svdup_f32(body[j].z);
svfloat32_t mj = svdup_f32(body[j].m);
svfloat32_t dx_0 = svsub_f32_x(p0, xj, xi_0);
svfloat32_t dy_0 = svsub_f32_x(p0, yj, yi_0);
svfloat32_t dz_0 = svsub_f32_x(p0, zj, zi_0);
svfloat32_t dx_1 = svsub_f32_x(p0, xj, xi_1);
svfloat32_t dy_1 = svsub_f32_x(p0, yj, yi_1);
svfloat32_t dz_1 = svsub_f32_x(p0, zj, zi_1);
• 正直こんなことをやって速くな
るとは思わなかった
svfloat32_t r2_0 = svmad_f32_x(p0, dx_0, dx_0, eps2);
r2_0 = svmad_f32_x(p0, dy_0, dy_0, r2_0);
r2_0 = svmad_f32_x(p0, dz_0, dz_0, r2_0);
• ループ分割することになったらや
ろうと思っていたアイデア
svfloat32_t r2_1 = svmad_f32_x(p0, dx_1, dx_1, eps2);
r2_1 = svmad_f32_x(p0, dy_1, dy_1, r2_1);
r2_1 = svmad_f32_x(p0, dz_1, dz_1, r2_1);
• 座標の差分（dx, dy, dz）は⻑
寿命の変数
svfloat32_t mri3_0 = rsqrtCubed(r2_0, mj, p0, one, a, b);
svfloat32_t mri3_1 = rsqrtCubed(r2_1, mj, p0, one, a, b);
• 再計算すれば使い捨てにできる
• *bodyと*body2を別引数として渡
しておけばコンパイラに⾒抜かれ
ずに再計算させられる
//
xj
yj
zj
mj
// subtract again
dx_0 = svsub_f32_x(p0, xj, xi_0);
dy_0 = svsub_f32_x(p0, yj, yi_0);
dz_0 = svsub_f32_x(p0, zj, zi_0);
cycle eff. Gflops
ACLE
16.1 52.8
75.5
ni32
14.8 57.6
82.3
recalc 14.3 59.5
85.1
この段階で理論限界は10 cycleにな
っているが8.5 cycleで効率は計算
load again
= svdup_f32(body2[j].x);
= svdup_f32(body2[j].y);
= svdup_f32(body2[j].z);
= svdup_f32(body2[j].m);
dx_1 = svsub_f32_x(p0, xj, xi_1);
dy_1 = svsub_f32_x(p0, yj, yi_1);
dz_1 = svsub_f32_x(p0, zj, zi_1);
ax_0 = svmla_f32_x(p0, ax_0, mri3_0, dx_0);
ay_0 = svmla_f32_x(p0, ay_0, mri3_0, dy_0);
az_0 = svmla_f32_x(p0, az_0, mri3_0, dz_0);
}
ax_1 = svmla_f32_x(p0, ax_1, mri3_1, dx_1);
ay_1 = svmla_f32_x(p0, ay_1, mri3_1, dy_1);
az_1 = svmla_f32_x(p0, az_1, mri3_1, dz_1);


# Page. 32

![Page Image](https://bcdn.docswell.com/page/VENM8NZVJ8.jpg)

無駄命令の削減
• svdup_f32()に対して2種類の出⼒が混ざっ
ている
• 良い例：
• ld1rw
{z25.s}, p0/z, [x4]
• 放送ロード命令、演算パイプを消費しない
• 悪い例：
cycle eff. Gflops
ACLE
16.1 52.8
75.5
ni32
14.8 57.6
82.3
recalc
14.3 59.5
85.1
asmtune
13.6 62.6
89.5
• ldr
s30, [x4, 4]
dup
z0.s, z0.s[0]
• スカラロードとレジスタ間放送、演算パイプを消費する
• 仕⽅ないので⼿動で置換した
• ldrで取れた負のオフセットがld1rwでは取れない
のが⾯倒
• 18箇所修正、休憩と検証を⼊れて1箇所数分程度


# Page. 33

![Page Image](https://bcdn.docswell.com/page/Y79NKRZV73.jpg)

最新の値
• 逆数平⽅根の旧実装（コンパイラ相当）と⽐較してみた
• 2サイクル以上の効果あり、コンパイラからACLE
• Broadcastの無駄命令の効果も顕著
• j-並列に戻した版も作った（asm使わないものでは最速）
•
•
•
•
•
ni16 old rsqrt
ni16 new rsqrt
ni16 recalc
ni32 recalc
ni2 j-parallel
: 18.4 cycles, 46.2% ! new impl
: 16.1 cycles, 52.7%
: 16.4 cycles, 51.7%
: 14.3 cycles, 59.5%
: 13.5 cycles, 63.0% ! new impl
• ni16 recalcasm : 14.3 cycles, 59.5%
• ni32 recalcasm : 13.0 cycles, 65.2%


# Page. 34

![Page Image](https://bcdn.docswell.com/page/G78KRW14JD.jpg)

まとめ
• コンパイラ版
• とりあえず素直にSoAに
• ni2 nj16 recalc
: 17.9 cycles, 47.4%
: 16.5 cycles, 51.5%
• ACLE版
• ni16 new rsqrt
• ni2 nj16 recalc
: 16.1 cycles, 52.7%
: 13.5 cycles, 63.0%
• asm編集版
• ni32 recalc asmtune
: 13.0 cycles, 65.2%
• 逆数平⽅根の後置補正がよく効いている（2サイクル以上）
• Bcastの無駄命令もアンロール前で2サイクル、後で1サイクル程度
• アンロールや再計算の細かなテクニックで残り1、2サイクル程度


