第12回 配信講義 計算科学技術特論A （2023）
2.4K Views
July 05, 23
スライド概要
第12回7月6日 A64FXでの重力N体計算カーネルのチューニング
重力多体問題、いわゆるN体計算のカーネル部分のA64FXプロセッサ向けチューニングについて紹介する。SVE命令を活用すべく、ACLEの組込関数を利用した。IntelアーキテクチャのAVX-512命令を用いた同等のコードとの性能比較も行う。
R-CCS 計算科学研究推進室
関連スライド
各ページのテキスト
A64FXでの重⼒N体計算 カーネルのチューニング 計算科学技術特論A（2023） 第12回 7⽉6⽇ 理化学研究所 運⽤技術部⾨ 技師 計算科学研究センター（ R-CCS） ソフトウェア開発技術ユニット 似⿃ 啓吾（にたどり けいご）
概要 • 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を仮定
公開コード、バージョン情報 • https://github.com/nitadori/nbody-a64fxにて公開 • 2021年から少しだけ整理 • 単なる作業場です • 使って何か研究してみたい⽅はご⼀報ください • 再現実験は断りなくどうぞ • バージョン情報 • 富⼠通⾔語環境： tcsds-1.2.36, FCC 4.8.0 (trad mode) • インテルコンパイラ（OneAPI）： icpc (ICC) 2021.2.0
構成 1. 2. 3. 4. 5. イントロ N体計算カーネルと逆数平⽅根 Intel AVX-512でのチューニング A64FX Arm SVEでのチューニング まとめ
背景・ねらい • 「富岳」のA64FXプロセッサ： • 命令セットのがSPARC（＋富⼠通拡張）→ARM（SVE拡張）へ • アーキテクチャレジスタ数が128→32に減少 • 命令レイテンシが9と⼤きい • 演算器を埋めきるスケジューリングが難しい • これまで：レジスタ数を活かしてコンパイル時の静的スケジューリン グ（ソフトウェアパイプライニング） • これから：CPUのOut-of-Order実⾏（動的スケジューリング）に頼ら ざるを得ない • 経験知：緩いSWPLを施したバイナリのほうが性能が出る • N体計算カーネル • 複雑すぎず簡単すぎない、メモリネックにならない • マイクロアーキテクチャのスケジューリング評価にはうってつけ
レイテンシとスループット • 最近のCPUの中⾝は流れ作業（パイプライン） • レイテンシ（Cycle per Instruction, CPI） • 演算結果が利⽤可能になるまで待つ時間 • ⽐喩：⾞を発注してから納品されるまでの⽇数 • スループット（CPIもしくはIPC, Instruction per Cycle） • その命令が演算器を専有する時間（CPI） • ⽐喩：⼀⽇に何台の⾞を⽣産できるか（IPC） • IPC（⼤きいほど良い）のほうが直感的なこともあるけどレイ テンシと単位を揃えてCPI（⼩さい⽅が良い）で書くことが多 い • 例：A64FXのFMA命令ならL-T=9-0.5 （9-1のユニットが2つあるとも）
例：Intel© Intinsic guideより
「理論限界」の考え⽅ • A64FXには2本のSIMD演算パイプラインがある • FLAとFLB、⼀部の命令はFLAのみが実⾏可能 • 演算をしないシャッフル系の命令やストア命令もFLAを消費すること に注意 • ロード命令や整数命令はこれらと同時発⾏できる • この2本で同時にFMA (fused multiply-add)命令が出たときにピーク のFLOPSが出る（⾏列乗算など） • 今回のN体カーネルは浮動⼩数点演算命令が17 • 内9がFMA • 逆数平⽅根の実装によっては+1、後述の再計算では+3命令 • 理想的にレイテンシが隠されたときの値として、8.5 cycleを基準に効 率を求めた
ループアンロールと ソフトウェアパイプライニング • ループの持つ並列度を利⽤して演算器を飽和させる • レジスタ本数に対するプレッシャーは増加 for(int i=0; i<n; i++){ A(i); B(i); C(i); } 基本のループ（擬似コード） for(int i=0; i<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<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<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); ソフトウェアパイプライニング
N体計算の基礎⽅程式 <latexit sha1_base64="6nPV07uEIx++bQdcSfomXsZGFL0=">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=</latexit> = <latexit sha1_base64="OA7lYFsX8Xqt787tKzsvFwk+wWE=">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==</latexit> 2 d r8 = ⌧ 2 dC # ’ 9=1 ⇣p m m r8 <9 Y 2 + |r 9 • εはsoftenin⻑といって⻑さ の次元を持つパラメタ • すべてのN 粒⼦への加速度を 計算するにはN 2回の相互作⽤ 計算が必要 • 天体の計算だけでなく、分⼦ 動⼒学法なども同じような式 のかたちをしている r8 |2 " ⌧ # ’ 9=1 ⌘ 3 (r 9 p <9 Y 2 + |r 9 r8 ) r8 |2 #!
include <cmath>
リファレンスコード
• 最内ループのメモリアクセ
スは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<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<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;
逆数平⽅根の重要性 • 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）
収束公式の導出 H 0 :' (1 + Y)G 1/= () (⌘ := 1 GH 0 = ) = 1 <latexit sha1_base64="F5DMtzbsYhv7gxiFzIVFMUqkK4g=">AAAEC3ichVG9b9NAFH+u+Sjho2lZkFhORI0SlYRzqABFQqpgYUHqB2kr1U1ku5f4VPvs2peQYOUfQGJmYAKJATHChNgQEjtiKCssiLFILAw8OyZt+sVZvnv33vv93u/dM32Hh5LSbWVMPXHy1OnxM5mz585fmMhOTi2HXjuwWM3yHC9YNY2QOVywmuTSYat+wAzXdNiKuXk3jq90WBByTzyQPZ+tu0ZL8Ca3DImuxqRyP096DUqqeshdtkUK2ozeMQLmh9zxRJF061FJuyb6RNczOm82SZ7oDmvKArFJ9TbRSIl0SYQM/bogesBbtiwS9Jf2EcXRIcNICLMLWsku1qND6/TiMrsqhhVSzMAdd4DAuBUtztddQ9qBG0mjhw/QTxWPQP4xDaG67duGkJ6L3WgYH7xHCtXIDLGQ207OCrHrFbR0a8OT4V6qRjZHyzRZ5KChpUYO0jXvZT+DDhvggQVtcIGBAIm2AwaE+K1BE20GW+i/ChpQ8DG+DhHGA7R4ksugDxnkaSOCYYaB3k3cW3hbS70C7zF/mKAtrOjgHyCSwDT9Ql/RHfqJvqY/6J8juaKEI9bVw9McYJnfmHh8aen3f1EunhLsXdSxmiV2fivRylG7n3jiLqwBvvPo6c5SdXE6ytMX9Cfqf0636QfsQHR+WS8X2OKzY/QI3B8iX3ePpqOzTczq42C1/WM8aCxXytqN8vWF2dzcnXTE43AZrkABZ3cT5uAezEMNLOWt8lX5pnxXn6hv1Hfq+0HqmJJiLsLIUj/+BYj8Bq8=</latexit> (1 + Y) = () 1 + Y = (1 ⌘) 1/= ⇣ ⌘ 1/= () H := G = (1 ⌘) 1/= H 0 ⇣ ⌘ H 1 := taylor (1 ⌘) 1/= H 0 ' 1 + 2 1 ⌘ + 2 2 ⌘2 + · · · H 0 ⼀般にn乗根の逆数に対して公式が存在する 1. ⼩さな値hを求める （正しく計算できていれば0） 2. 多項式𝑝 ℎ を計算 3. 近似値に𝑝 ℎ を掛ける k次の多項式なら有効桁数は 𝑘 + 1 倍に
具体例
• 逆数
<latexit sha1_base64="wiBiQREPOGPZUcJDJhsni7Jlx/c=">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=</latexit>
参考： 逆数と平⽅根を求める⾼次収束アルゴリズム
http://www.finetune.co.jp/~lyuka/technote/fract/sqrt-jp.html
H 0 :' 1/G
⌘ := 1
GH 0
H 1 := H 0 + H 0 ⌘ + ⌘2 + ⌘3 + · · ·
• 逆数平⽅根
<latexit sha1_base64="zto1s2pYx8J+JNmHuC//zhtWe8Q=">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==</latexit>
p
H 0 :' 1/ G
GH 0 2
✓
◆
1
3 2 5 3
H 1 := H 0 + H 0 ⌘ + ⌘ + ⌘ + · · ·
2
8
16
⌘ := 1
Break： • 係数を求める程度ならブラ ウザと⽇本語でOK • 昔ならMaxima、最近なら SymPyが便利かも • Chat-GPTはWolframのプラ グインがあればこのページ と同じになる
N体計算に有効な⽅法
• 逆数平⽅根を補正してから3乗するよりは、3乗してから補正す
るほうが少し有利
• 演算数が1つ減る（𝑟 !" に再利⽤性がある）
• 直列な命令依存が減り命令レベル並列度が上がる（7段→5段）
• 補正多項式は 1 − ℎ !#/" のテイラー展開
Ã
1
p
:' 1/ A 2
Ã
2
:= Ã
<latexit sha1_base64="FH2Elrj/31zpvm/zSUmyvKF3C4g=">AAAESXichVHLbtNAFL2xC5TwaFo2ldgMjVohtQ3jtEAVCamCDcs+SFupbiPbmTSmfmU8CRTLPwAfwIIVSCwQn8EGiSVi0U9ArKBIIMSC60faxGnKWJ65c+ac43M9umeZvqD0MCfJI+fOXxi9mL90+crVscL4xIbvtrnBqoZruXxL13xmmQ6rClNYbMvjTLN1i23q+w+i880O477pOo/Egcd2bG3PMRumoQmEauO5FzNEFaZVZ4TvBvNKSCqqb9qsRZRbqt/iIuC75ZCoar6Ph1DlXkaIO5v5WTASNiOyQuYRLJ9C6/p7iWeDa4ZSJrMoO+ZGGAkWwqU5orbaWj3zlebQNM3YudtaljQ74JPKPHKsS1rtd0123Sx2j79de9xPSowSykJMOeH3EON/UCsUaYnGgwwWSloUIR0rbuEzqFAHFwxogw0MHBBYW6CBj882NLBm0EJ8DhSg4OH5DgR4zrEyYy6DEPLo00YFQ4aG6D7Oe7jbTlEH95G/H6sN/KKFL0clgWn6hb6jR/QjfU+/0r9DvYLYI8p1gKueaJlXG3s+uf7rvyobVwHNE9WZmQV2vhRnNTG7FyNRF0ai7zx7ebReWZsOZugb+g3zv6aH9AN24HR+Gm9X2dqrM/I4OD9Bv6c9mYazdWSFeLFK9hoHi41ySblTur26WFy+n17xKFyHKbiJd3cXluEhrEAVjNwPaVK6IU3Jn+Tv8m/5T0KVcqnmGvSNEfkf8Cka+g==</latexit>
⌘ := 1
? :=
1
2
1
A 2 ⇥ Ã
+ ⌘ ⇥ 38 ,
A
1
:= Ã
1
A
2
:= A
1
<A
3
⇥ Ã
+ Ã
⇥A
Ã
1
p
:' 1/ A 2
1
Ã
2
:= Ã
2
⌘ := 1
<latexit sha1_base64="ezYgX7H43PBpF7Q+6HNh4qu9tEo=">AAAEtnichVG7b9NAGP/iGijh0bQsSCwnolZFbcPZLRBFqlTBwtgHaSvVTWQ7l9jUr9iXQLEysjDAyMAEEgPiz2BhRwz9ExBjkWBg4POjSeOk5Sz7vvvd93ucT/MsM+CUHuWECfHCxUuTl/NXrl67PlWYntkO3I6vs6ruWq6/q6kBs0yHVbnJLbbr+Uy1NYvtaAePov2dLvMD03We8EOP7dtqyzGbpq5yhOrTuT9zROGm1WDEr4VLUo9UlMC0WZtId5Wg7fPQr8k9oij5oT6EKqsZIq5sFgyDi0Rpd9QGsbMmqwk0QkicjKhBIksIymPa5IFu/enpjeVYecRtnEJiZGfYsfGoaCpgnNh6yembvqovy2QBaX2LCCOhdK9XPnEYBMu6LYzxT4W89JfHdorFmnye2GkTCmXiJfsSCvZDGf1FEsaoyYpvtgx+J53qhSIt0XiQ0UJKiyKkY90tfAMFGuCCDh2wgYEDHGsLVAjw2YMm1gzaiC+CBBQ83N+HEPd9rMy4l0EP8qjTQQbDDhXRA/y2cLWXog6uI/0gZuvoaOHrI5PALP1OP9Fj+pV+pj/o3zO1wlgjynWIs5ZwmVefenVz6/d/WTbOHIwB69zMHE9ejrOamN2LkegUesLvvnh7vFXZnA3n6Af6E/O/p0f0C57A6f7SP26wzXfn5HHw+wz1np/KdHa3hl09vFgpe42jxbZcku6XVjZWimsP0yuehFtwG+bx7h7AGjyGdaiCLtSFl8Jr4Y1YFmsiE1tJq5BLOTdgaIjeP7XcP1Q=</latexit>
Ã
1
1
1
⌘ := Ã
⌘⇥?
,
<A
1
1
⇥⌘
:= < 9 ⇥ A
:= <A 1 ⇥ A 2
7 stage, 9 ops
1
<Ã
3
<A
3
1
⇥ Ã
1
2
,
⌘ := < 9 Ã
3
⇥ ⌘,
A 2 ⇥ Ã
:= <Ã
<A
3
3
,
<Ã
< 9 Ã
+ <Ã
= <Ã
3
1
3
:= < ⇥ Ã
3
:= <Ã
? :=
⌘⇥?
1 + 32 ⌘ +
3
2
1
⇥ Ã
+⌘⇥
15 2
8 ⌘
5 stage, 8 ops
1
15
8
2
外側ループと内側ループの並列度
• 外側ループ
• 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<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<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;
コンパイラの対応状況 • Intel • ディレクティブを付けることで外側（i ）ループのSIMD化にも対応 • 富⼠通 • 最内（j ）ループのみ⾃動的にSIMD化 • 最後の⽔平総和も出⼒してくれる • 今回は双⽅のアーキテクチャに向けて外側ループの⼿動SIMD 化を⾏った
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 をループの前に付けるだけ
ソースコードの微変更
Before
After
#pragma omp simd
for(int i=0; i<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<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<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<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
放送ロードがメモリオペランドになっている
改善できそうな点 ..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⼀発 だった
Intrinsic版
{
const __m512 eps2 = _mm512_set1_ps(eps2_ss);
for(int i=0; i<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<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 <x86intrin.h>
が必要
__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);
}
}
性能値 • 内側ループ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（外側）
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
構造体アクセス命令とレジスタ内転置 • 外側ループであるが、 {x,y,z,m}[16] -> {x[16],y[16],z[16],m[16]} {ax[16],ay[16],az[16]} –> {ax,ay,az}[16] のような転置がある • SVEにはld4、st3のようなstructure load/store命令がある • AVX-512にそういうのはないのでレジスタ内で転置する • 今回はボトルネックではないけどgather/scatteするのも悔しいので • ⼀度作っておくと使い回しが効く • 命令セットの格好の練習問題 • MDなど粒⼦をリストアクセスするケース、body[j]ではなく body[idx[j]]が必要になるケースには有⽤
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)
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)
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命令のループなのでコンパイラ任せで理論限界の半分
最新の値 • 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%
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<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<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);
スケジューリングの 為の再計算 for(int j=0; j<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);
無駄命令の削減 • 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箇所数分程度
最新の値 • 逆数平⽅根の旧実装（コンパイラ相当）と⽐較してみた • 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%
まとめ • コンパイラ版 • とりあえず素直に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サイクル程度