第1回 配信講義 計算科学技術特論A (2023)

37.8K Views

April 12, 23

スライド概要

第1回 4月13日 プログラム高速化の基礎
プログラム高速化の基礎知識、並列化プログラミング(MPI、OpenMP)の基礎知識、およびプログラム高速化の応用事例の座学を通して、計算科学で必要な高性能計算技術の基礎の習得を目指す。
https://www.r-ccs.riken.jp/outreach/schools/20230413-1/

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

関連スライド

各ページのテキスト
1.

内容に関する質問は [email protected] まで 第1回 プログラム高速化の基礎 名古屋大学情報基盤センター 1 2023年度 計算科学技術特論A 片桐孝洋

2.

本講義の位置づけ 2 2023年度 計算科学技術特論A

3.

講義日程と内容について  2023年度 計算科学技術特論A(木曜:13:00-14:30 )      3 第1回:プログラム高速化の基 礎、2023年4月13日  イントロダクション、ループアンローリング、キャッシュブロック化、 数値計算ライブラリの利用、その他 第2回:MPIの基礎、2023年4月20日  並列処理の基礎、MPIインターフェース、MPI通信の種類、その他 第3回:OpenMPの基礎、2023年4月27日  OpenMPの基礎、利用方法、その他 第4回:Hybrid並列化技法(MPIとOpenMPの応用)、2023年5月11日  背景、Hybrid並列化の適用事例、利用上の注意、その他 第5回:プログラム高速化実例と大規模学習への展開、2023年5月18日  プログラムの性能ボトルネック に関する考えかた(I/O、単体性能 (演算機ネック、メモリネック)、並列性能(バランス))、性能プロファイル、 機械学習におけるHPC、ほか 2023年度 計算科学技術特論A

4.

参考書  「計算科学のためのHPC技術1 」 下司雅章 (編集), 片桐孝洋 , 中田真秀, 渡辺宙志, 山 本有作, 吉井範行, Jaewoon Jung, 杉田 有治, 石村和 也, 大石進一, 関根晃太, 森倉悠介, 黒田久泰,著 出版社: 大阪大学出版会 (2017/4/3) ISBN-10: 4872595866, ISBN-13: 978-4872595864 発売日: 2017/4/3 【本書の特徴】  計算科学に必要なHPC技術について、基礎的な事 項を解説している  片桐担当(1章~5章)  プログラム高速化の基礎、MPIの基礎、OpenMP の基礎、Hybrid並列化技法(MPIとOpenMPの応 用)、プログラム高速化の応用      4 2023年度 計算科学技術特論A

5.

参考書(演習書)  「スパコンプログラミング入門 -並列処理とMPIの学習-」  片桐 孝洋 著、  東大出版会、ISBN978-4-13-062453-4、 発売日:2013年3月12日、判型:A5, 200頁  【本書の特徴】  C言語で解説  C言語、Fortran90言語のサンプルプログラムが付属  数値アルゴリズムは、図でわかりやすく説明  本講義の内容を全てカバー  内容は初級。初めて並列数値計算を学ぶ人向けの 入門書 5 2023年度 計算科学技術特論A

6.

参考書(演習書)  「並列プログラミング入門: サンプルプログラムで学ぶOpenMPとOpenACC」  片桐 孝洋 著  東大出版会、ISBN-10: 4130624563、 ISBN-13: 978-4130624565、発売日: 2015年5月25日  【本書の特徴】  C言語、Fortran90言語で解説  C言語、Fortran90言語の複数のサンプルプログラムが入手可能 (ダウンロード形式)  本講義の内容を全てカバー  Windows PC演習可能(Cygwin利用)。スパコンでも演習可能。  内容は初級。初めて並列プログラミングを学ぶ人向けの 入門書 6 2023年度 計算科学技術特論A

7.

参考書  「スパコンを知る: その基礎から最新の動向まで」  岩下武史、片桐孝洋、高橋大介 著  東大出版会、ISBN-10: 4130634550、 ISBN-13: 978-4130634557、 発売日:2015年2月20日、176頁  【本書の特徴】  スパコンの解説書です。以下を 分かりやすく解説します。  スパコンは何に使えるか  スパコンはどんな仕組みで、なぜ速く計算できるのか  最新技術、今後の課題と将来展望、など 7 2023年度 計算科学技術特論A

8.

参考書  「数値線形代数の数理とHPC (シリーズ応用数理 第 6巻)」  日本応用数理学会(監修)、櫻井鉄也、 松尾宇泰、片桐孝洋(著)  出版社: 共立出版 (2018/8/30)、 ISBN-10: 4320019555、発売日: 2018/8/30  【本書の特徴】  スパコンの解説書です。以下を 分かりやすく解説します。   8 前半:連立一次方程式の数値解法,行列の固有値問題および特異値問題 の数値解法,最小二乗問題の数値解法,行列関数の数値解法 後半:連立一次方程式の解法や,固有値および特異値の計算のスーパーコ ンピュータを利用する上で必要となる,データ分散・並列化・前処理・通信量 の削減の方法。HPCにおける計算手法や実装方法 2023年度 計算科学技術特論A

9.

参考書  「並列数値処理 - 高速化と性能向上のために -」  金田康正 東大教授 理博 編著、 片桐孝洋 東大特任准教授 博士(理学) 著、黒田久泰 愛媛大准教授 博士(理学) 著、山本有作 神戸大教授 博士(工学) 著、 五百木伸洋 ㈱日立製作所 著、  コロナ社、発行年月日:2010/04/30 , 判 型: A5, ページ数:272頁、 ISBN:978-4-339-02589-7, 定価:3,990円 (本体3,800円+税5%)  【本書の特徴】  Fortran言語で解説  数値アルゴリズムは、数式などで厳密に説明  本講義の内容に加えて、固有値問題の解法、疎行列反復解法、 FFT、ソート、など、主要な数値計算アルゴリズムをカバー  内容は中級~上級。専門として並列数値計算を学びたい 人向き 9 2023年度 計算科学技術特論A

10.

イントロダクション スパコンとは何か? 10 2023年度 計算科学技術特論A

11.

スーパーコンピュータとは  一般に「人工知能搭載のコンピュータ」とは言われていないが・・・  最近のスパコンはAIスパコン  ChatGPT :大規模言語モデル(LLM:Large Language Model)はGPU利用が必須  産総研「ABCI (AI Bridging Cloud Infrastructure)」  名古屋大学 スーパーコンピュータ「不老」TypeⅡサブシステム、など多数設置済み  明確な定義はない  現在の最高レベルの演算性能をもつ計算機のこと  経験的には、PCの1000倍以上 高速で、1000倍以上大容量メモリをもつ計算機  1000倍高速だと 世界が違う!  人の歩行速度:時速 約5km  1000倍だと 時速 5000km  ジェット旅客機:速くて 時速1000km  <ジェット旅客機の5倍の速度> と <歩く速さ>  最新鋭スパコンの能力はPCの10万倍以上高速 11 2023年度 計算科学技術特論A

12.

スーパーコンピュータとは  明確な定義はない  現在の最高レベルの演算性能をもつ計算機のこと  経験的には、PCの1000倍高速で、1000倍大容量なメモリをもつ計算機 輸出貿易管理令別表第一及び外国為替令別表の規定に 基づき貨物又は技術を定める省令(平成三年通商産業省令 第四十九号): 施行日: 令和四年十二月六日  (令和四年経済産業省令第七十八号による改正)  https://elaws.e-gov.go.jp/document?lawid=403M50000400049  第七条第三項ロ:デジタル電子計算機であって、 加重最高性能が七○実効テラ演算を超えるもの  現在、ほとんどすべてのスーパーコンピュータは並列計算機  名古屋大学情報基盤センターが所有するスーパーコンピュータ「不老」 も並列計算機 12 2023年度 計算科学技術特論A

13.

スーパーコンピュータで用いる単位 問)実効テラ演算とは・・・  答)TFLOPS(テラ・フロップス、 Tera Floating Point Operations Per Second)   1秒間に1回の演算能力(浮動小数点)が1FLOPS。       演算とは:足し算、引き算、かけ算、割り算、どれも 1回と計算する K(キロ)は1,000(千) M(メガ)は1,000,000(百万) G(ギガ)は1,000,000,000(十億) T(テラ)は1,000,000,000,000(一兆) 一秒間に一兆回の浮動小数点演算の能力がある こと。 13 2023年度 計算科学技術特論A

14.

スーパーコンピュータで用いる単位  PFLOPS(ぺタ・フロップス)    1秒間に0.1京(けい)回の浮動小数点演算能力があること。 「京」コンピュータ(2012年9月共用開始、11.2 PFLOPS、退役済) スーパーコンピュータ「富岳」(537 PFLOPS (倍精度、 ブーストモード)、2021年3月共用開始)  PCの演算能力は?  3.3GHz(1秒間に3.3G回のクロック周波数)として、 もし1クロックあたり1回の浮動小数点演算ができれば 3.3GFLOPS。  Intel Core i7 (Sandy Bridge)(2011年頃)では、 6コア、1クロックで8回の浮動小数計算ができるので、 3.3 GHz * 8回浮動小数点演算/Hz * 6コア = 158.4 GFLOPS  Cray-1は160MFLOPS。 1970年代のスパコンより、 PCの方が990倍以上高速! 14

15.

スーパーコンピュータ用語  理論性能(Theoretical Performance)    ハードウエア性能からはじき出した性能。 1クロックに実行できる浮動小数点回数から算出した FLOPS値を使うことが多い。 実効性能(Effective Performance)     15 何らかのベンチマークソフトウエアを実行して実行時間を計測。 そのベンチマークプログラムに使われている浮動小数点演算 を算出。 以上の値を基に算出したFLOPS値のこと。 連立一次方程式の求解ベンチマークであるLINPACKを 用いることが多い。 2023年度 計算科学技術特論A

16.

ムーアの法則  米Intel社の設立者ゴードン・ムーアが提唱した、半導体技術 の進歩に関する経験則。 「半導体チップの集積度は、およそ18ヵ月で2倍になる」  これから転じて、 「マイクロプロセッサの性能は、およそ18ヵ月で2倍になる」  上記によると、約5年で10倍となる。 16 2023年度 計算科学技術特論A

17.

主要なスーパーコンピュータ性能推移 Summit (DOE/SC/Oak Ridge) Frontier (理論性能) Sierra (DOE/NNSA/LLNL) (1 EXA FLOPS) (1.68 EXA FLOPS) Sunwei TaifuLight (NRCPC) Titan (DOE/SC/ORNL) Sequoia(DOE/NNSA/LLNL) スーパーコンピュータ K-Computer (RIKEN) 「富岳」 (552PFLOPS) Tianhe-1A(NUDT) Jaguar(ORNL) FX100(名大) TUBAME(東工大) T2K(東大) SR11000(東大) SX-8 FX1(JAXA) SR8000(東大) SX-7 E2S(地球Sim) SX-4 SX-6 地球シミュレータ SX-5 SX-4 VP-200 SX-3 SR2201(東大) SX-2 VP-2600 S-810 S-820 Cray-1 VP-200 ILLIAC-IV FACOM230 ENIAC 17 2023年度 計算科学技術特論A

18.

スーパーコンピュータのランキング  TOP500 Supercomputer Sites (http://www.top500.org/)    18 LINPACKの値から実効性能を算出した値の 500位までのランキング 米国オークリッジ国立研究所/テネシー大学 ノックスビル校の Jack Dongarra 教授が発案 毎年、6月、11月(米国の国際会議SC|xy) で発表 2023年度 計算科学技術特論A

19.

現在のランキング (2022年11月) https://www.top500.org/lists/top500/list/ 2022/11/ 19 2023年度 計算科学技術特論A

20.

1位 Frontier(米国)  1.68 EXA FLOPS (Theoretical) 2022.6- (ISC22)  米国エネルギー省(DOE)、科学局(SC)、オークリッジ国立研 究所  Power :21.10 MW Theoretical Peak:1.685 EXA FLOPS Linpack 1.110 EXA FLOPS (65% to theoretical peak ) 8,730,112 cores (約873万コア)⇒並列性 AMD Optimized 3rd Generation EPYC 64C (2GHz) 52.2 GFLOPS/Watt ⇒TOPレベルの電力性能      ※ Frontier TDS:62.6 GFLOPS/Watt 2023年度 計算科学技術特論A 20

21.

2位 スーパーコンピュータ「富岳」(日本)  537.21 PFLOPS (Theoretical) 2020.11- (SC20)        理化学研究所 計算科学研究センター(R-CCS) Power :29.89 MW Theoretical Peak:537.21 PFLOPS Linpack 442.01 PFLOPS ( 82% to theoretical peak ) 7,630,848 cores (約763万コア) ARM A64FX 48C (2.2GHz) 14.8 GFLOPS/Watt ※LINPACK稼働時の電力で算出ではない。 参考:15.4 GFLOPS/Watt (2020/11, Green500 List) Source: https://www.r-ccs.riken.jp/fugaku/system/ 2023年度 計算科学技術特論A 21

22.

3位 LUMI (フィンランド)  428.70PFLOPS (Theoretical) 2022.6- (ISC22)        EuroHPC/CSC Power :6.01 MW Theoretical Peak:428.70 PFLOPS Linpack 309.10 PFLOPS ( 72% to theoretical peak ) CPU: 2,220,288 cores (約222万コア) AMD Optimized 3rd Generation EPYC 64C (2GHz) 51.3 GFLOPS/Watt 2023年度 計算科学技術特論A 22

23.

4位 Leonardo (イタリア)  200.79PFLOPS (Theoretical) 2022.6- (ISC22)        EuroHPC/CINECA Power :5.61 MW Theoretical Peak:255.75 PFLOPS Linpack 174.7 PFLOPS ( 68% to theoretical peak ) CPU: 1,463,616 cores (約146万コア) CPU: Xeon Platinum 8358 32C (2.6GHz) + GPU: NVIDIA A100 SXM4 64 GB 31.1 GFLOPS/Watt 2023年度 計算科学技術特論A 23

24.

22位 AI Bridging Cloud Infrastructure (ABCI 2.0) (日本)  54.34 PFLOPS (Theoretical) 2018.6- (ISC18)          AI Bridging Cloud Infrastructure (ABCI) 2.0 産業技術総合研究所 Power :1.60MW Theoretical Peak:54.34 PFLOPS Linpack 22.21 PFLOPS (40% to theoretical peak) 504,000 cores (約54万コア) CPU: Xeon Gold 6148 20C (2.4GHz) + GPU: NVIDIA A100 (960基) + V100(4,352基) Interconnect: Infiniti Band HDRx4 13.8 GFLOPS/Watt Source: https://abci.ai/ja/ 2023年度 計算科学技術特論A 24

25.

数値計算+AIの融合スパコンの登場   CPUのみのサブマシンと、GPU主体のサブマシンを、 ファイルシステムで連結させて処理する「複合型スパコン」 数値シミュレーションの結果を用いて、AI処理をすることが可能    大規模言語モデル(LLM:Large Language Model)も 各サブシステム間で、シームレスなデータ移動 代表的なスパコン(国内)  名古屋大学 スーパーコンピュータ「不老」 (2020年7月1日~)  大阪大学 SQUID (2021年5月1日~)  東京大学 Wisteria / BDEC-01 (2021年5月14日~) 25 2023年度 計算科学技術特論A

26.

マルチコアとメニーコア  いわゆる、CPU (Central Processing Unit)  マルチコアCPU 例)マルチコアCPU (Intel Xeon Gold)  低電力化のため動作周波数を落として コア(CPU)をたくさん並べる  通常は8(PCレベル) ~ 32~80個 ☞スーパーコンピュータ「富岳」では48個/ノード  メニーコアCPU 例)メニーコアCPU (Intel Xeon Phi)  低電力化のため動作周波数をすごく落として コア(CPU)をもっとたくさん並べる  通常は60個以上、動作時には240並列以上 ☞生産停止 2023年現在、コア数的に、マルチコアと(旧)メニーコア の差がもうない 26 2023年度 計算科学技術特論A

27.

GPU (Graphics processing Unit) 例)NVIDIA A100、H100  ゲームとかで使われる グラフィックス用の演算加速器(GPU)を、 数値計算に使う  GPGPU (General Purpose GPU )  低電力化のため、すごく周波数が低い計算要素を、 すごく並べる  通常、1万~10万要素  単体では使えない  CPUと組み合わせて使う  そのため、演算加速器と呼ばれる  使うためには、専用言語が必要  NVIDIA CUDAなど 2023年度 計算科学技術特論A 27

28.

NVIDIA Hopper (H100)      (Source: https://www.nvidia.com/ja-jp/datacenter/h100/ https://www.nvidia.com/ja-jp/data-center/h100/ ) DP(FP64): 34 TFLOPS SP(FP32) : 67 TFLOPS Specialized for AI processing. NVIDIA Grace Hopper CPU+GPU Architecture. Half precision: 1.979 PFLOPS Transformer Engine(TE): FP16 + FP8 (newly added): for BERT and GPT-3  4x4 matrix-matrix-multiplications (for Tensor Core)      D = A B + C (Input: FP16, out: FP16 or FP32) Input (2x FP16), mult (Full precision), addition (FP32), Output (FP32) FP16 addition mode is supported. 14,592 cuda cores GPU memory: 80 GB (3.35 TB/s) 28 2023年度 計算科学技術特論A

29.

単体(CPU)最適化の方法 29 2023年度 計算科学技術特論A

30.

高速 最近の計算機のメモリ階層構造 O(1ナノ秒) O(10ナノ秒) レジスタ バイト Kバイト キャッシュ ~Mバイド O(100 ナノ秒) メインメモリ Mバイト ~Gバイド O(10 ミリ秒) ハードディスク Gバイト ~Tバイト 大容量 <メインメモリ>→<レジスタ>への転送コストは、 レジスタ上のデータ・アクセスコストの O(100)倍! 30 2023年度 計算科学技術特論A

31.

より直観的には… レジスタ キャッシュ メインメモリ 高性能(=速い)プログラミングをするには、 きわめて小容量のデータ範囲について 何度もアクセス(=局所アクセス)するように ループを書くしかない 31 2023年度 計算科学技術特論A

32.

Fujitsu FX1000のメモリ構成例 高速 レジスタ レベル1キャッシュ (64Kバイト/1コア) 11+ [TB/s] ●データ ●データ レベル2キャッシュ (8Mバイト/12コア 合計:32Mバイト) 3.6+ [TB/s] メインメモリ (32Gバイト/ノード) 32 ●データ 2023年度 計算科学技術特論A 大容量

33.

Fujitsu FX1000のメモリ構成例 高速 レジスタ レベル1キャッシュ (64Kバイト/1コア) 11+ [TB/s] ●データ ●データ レベル2キャッシュ (8Mバイト/12コア 合計:32Mバイト) 3.6+ [TB/s] データが L1キャッシュ上 にあれば、 速くアクセス可能 大容量 メインメモリ (32Gバイト/ノード) 33 2023年度 計算科学技術特論A

34.

Fujitsu FX1000のノードのメモリ構成例 ※12コア単位 × 4ソケット相当(CMG) コア0 コア1 コア2 コア3 コア0 コア1 コア2 コア3 コア0 コア1 コア2 コア0 コア1 コア2 コア3 コア3 コア コア コア コア コア コア コア コア コア コア コア コア 88コア 99コア 10 11 コア コア 10 11 8 9 10 11 8 9 10 L1 L1 L1 L111 … L1 L1 L1 L1 … L1 L1 L1 L1 L1 L1 L1 L1 … L1 L1 L1 L1 L1 L1 L1 L1L2… L1 L1 L1 L1 L1 L1 L1 L1 L2 L2 L2 メインメモリ ※階層メモリ構成となっている 34 2023年度 計算科学技術特論A

35.

Fujitsu FX1000全体メモリ構成 コ コ コ コ アコ アコ アコ アコ コ コ コ ア ア ア 0 1 2 3アコ 0アコ 1アコ 2アコ 3アコ L 0ア L 1ア L 2ア L 3ア L0 L1 L2 L3 コ コ コ コ アコ アコ アコ アコ コ ア 1 1アコ コ コ ア ア 8 9 1アコ 1アコ 8アコ 9アコ 0 1ア 1 1ア 0 1 … L 8ア L 9ア L 01 L 11 …1L 8 1L 9 1L 0 1L 1 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 コ コ コ コ アコ アコ アコ アコ コ コ コ ア ア ア 0 1 2 3アコ 0アコ 1アコ 2アコ 3アコ L 0ア L 1ア L 2ア L 3ア L0 L1 L2 L3 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 メインメモリ メインメモリ コ コ コ コ アコ アコ アコ アコ 0アコ 1アコ 2アコ 3アコ コ コ コ 0ア 1ア 2ア 3アコ 1ア L 2ア L 3ア L 0ア L L0 L1 L2 L3 コ コ コ コ アコ アコ アコ アコ 1アコ 1アコ 8アコ 9アコ コ 1ア 1アコ 8アコ 9アコ 0 1ア 1 1ア 0 1 ア ア … L 88 L 99 L 01 L 11 …1L 1L 1L 0 1L 1 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 コ コ コ コ アコ アコ アコ アコ コ ア 1 1アコ コ コ ア ア 8 9 1アコ 1アコ 8アコ 9アコ 0 1ア 1 1ア 0 1 … L 8ア L 9ア L 01 L 11 …1L 8 1L 9 1L 0 1L 1 コ コ コ コ アコ アコ アコ アコ 0アコ 1アコ 2アコ 3アコ コ コ コ 0ア 1ア 2ア 3アコ 1ア L 2ア L 3ア L 0ア L L0 L1 L2 L3 コ コ コ コ アコ アコ アコ アコ 1アコ 1アコ 8アコ 9アコ コ ア 1 1アコ 8アコ 9アコ 0 1ア 1 1ア 1 ア ア 0 8 9 … L 8 L 9 L 01 L 11 …1L 1L 1L 0 1L 1 メインメモリ コ コ コ コ アコ アコ アコ アコ 1アコ 1アコ 8アコ 9アコ コ ア 1 1アコ 8アコ 9アコ 0 1ア 1 1ア 1 ア ア 0 8 9 … L 8 L 9 L 01 L 11 …1L 1L 1L 0 1L 1 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 コ コ コ コ アコ アコ アコ アコ 0アコ 1アコ 2アコ 3アコ 0アコ 1アコ 2アコ 3アコ L 0ア L 1ア L 2ア L 3ア L0 L1 L2 L3 コ コ コ アコ アコ アコ 1アコ 8アコ 9アコ 1アコ 8アコ 9アコ 0 1ア ア ア 0 8 9 … L 8 L 9 L 01 …1L 1L 1L 0 コ アコ 1アコ 11アコ 11ア L 11 L 11 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 メインメモリ コ コ コ コ アコ アコ アコ アコ コ コ コ 0ア 1ア 2ア 3アコ 0アコ 1アコ 2アコ 3アコ L 0ア L 1ア L 2ア L 3ア L0 L1 L2 L3 コ コ コ コ アコ アコ アコ アコ コ 1ア 1アコ 8アコ 9アコ コ ア 1 1アコ 8アコ 9アコ 0 1ア 1 1ア 1 ア ア 0 8 9 … L 8 L 9 L 01 L 11 …1L 1L 1L 0 1L 1 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 メインメモリ メインメモリ 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 メインメモリ コ コ コ コ アコ アコ アコ アコ コ コ コ 0ア 1ア 2ア 3アコ コ コ コ ア ア ア 0 1 2 3アコ 1ア L 2ア L 3ア L 0ア L L0 L1 L2 L3 コ コ コ コ アコ アコ アコ アコ 0アコ 1アコ 2アコ 3アコ コ コ コ 0ア 1ア 2ア 3アコ 1ア L 2ア L 3ア L 0ア L L0 L1 L2 L3 コ コ コ アコ アコ アコ 1アコ 8アコ 9アコ 1アコ 8アコ 9アコ 0 1ア 9ア L0 1 … L 8ア L 0 …1L 8 1L 9 1L 0 コ アコ 1アコ 11アコ 11ア L 11 1L 1 1L 1L 1L 1L L L L L 1L 1L 1L 1L … …11L 11L 11L 11L 1 1 1 1 1 1 1 1L2 1 1 1 1 L2 L2 L2 メインメモリ TofuインターコネクトD メモリ階層が階層 35 (各ノードは周囲の隣接ノードへ同時に合計 40.8 GB/s×双方向 (1リンク当たり 6.8 GB/s × 双方向、6リンク同時通信可能)) 2023年度 計算科学技術特論A

36.

FX100計算ノードの構成 HMC 16GB TOFU2 Network 2ソケット相当、NUMA (Non Uniform Memory Access) Memory L2 (17コアで共有、12MB) L1 L1 L1 L1 Core #0 Core #1 Core #2 Core #3 : L1データ キャッシュ 64KB … L1 … Assist. Core … ICC L1 L1 L1 L1 Core #12 Core #13 Core #14 Core #15 Core #28 Core #29 Core #30 Core #31 L1 L1 L1 L1 ソケット0 (CMG(Core Memory Group)) Core #16 Core #17 Core #18 Core #19 L1 L1 L1 L1 ソケット1 (CMC) : L1データ キャッシュ … 64KB Assist. Core … … L1 L2 (17コアで共有、12MB) Memory HMC 16GB ノード内合計メモリ量:32GB 36 読込み:240GB/秒 書込み:240GB/秒=合計:480GB/秒

37.

FX100(名大)のCPU(SPARC64XIfx)の詳細情報 項目 値 アーキテクチャ名 HPC-ACE2 (SPARC-V9命令セット拡張仕様) 動作周波数 2.2 GHz L1キャッシュ 64 Kbytes (命令、データは分離) L2キャッシュ 24 Mbytes ソ フ ト ウ ェ ア 制 御 セクタキャッシュ キャッシュ 演算実行 2整数演算ユニット、8つの浮動小数点積和演算ユニット(FMA) SIMD命令実行 1命令で2つのFMAが動作 FMAは4つの浮動小数点演算(加算と乗算)を実行可能 レジスタ  浮動小数点レジスタ数:256本 その他 37 2023年度 計算科学技術特論A

38.

FX1000計算ノードの構成 4ソケット相当、NUMA (Non Uniform Memory Access) Tofu D Network Memory Memory HMC2 8GB L2 (12コアで共有、12MB) L2 (12コアで共有、12MB) L1 Core #0 L1 L1 L1 L1Core L1Core L1Core #1 #2 #3 Core Core Core #0 #1 #2 : L1データ L1 L1 キャッシュ … L1… Core : L1データ Assist. L1 64KBキャッシュ … Core Assist. … …#9 Core #3 64KB … Core ソケット0 (CMG(Core Memory Group)) ソケット2 (CMG) Assist. Core Core Core Core : L1データ Core Assist. #12 Core #13 Core #14 Core #15 Core … : L1データL1 Core L1 #12 L1 #13 L1 #14 L1 #15キャッシュ … L1 L1 ソケット1 (CMG) ソケット3 (CMG) HMC2 8GB L1 … キャッシュ L164KB … 64KB … L1 L1 38 L1 Core Core Core Core #20 Core #21 Core #22 Core #23 Core #20L1 #21 L1 #22 L1 #23 L1 … L1 ノード内合計メモリ量:32GB 読込み:512GB/秒 書込み:512GB/秒=合計:1024GB/秒 L1 L1Core L1Core L1Core L1 #9 #10 #11 Core Core Core Core #9 #9 #10 #11 L2 (12コアで共有、12MB) L2 (12コアで共有、12MB) Memory Memory ICC L1 L1 L1

39.

FX1000のCPU(A64FX)の詳細情報 項目 値 アーキテクチャ名 Armv8.2-A 動作周波数 2.2 GHz L1キャッシュ 64 Kbytes L2キャッシュ 32 Mbytes ソ フ ト ウ ェ ア 制 御 セクタキャッシュ キャッシュ  512 bit wide SIMD 演算実行 (Armv8-A Scalable Vector Extension, SVE) • 128 ~ 2048 bit 可変 • 倍精度、単精度、に加えて、半精度演算をサポート • 8, 16, 32, 64 bit 整数ベクトル演算対応 SIMD命令実行 512 wide SIMD レジスタ  32本 (Scalable Vector Register, 128~2048 bit)  16本 (Predicate Register,16~256 bit) その他  Gather/Scatter命令 (Source: FUJITSU Supercomputer PRIMEHPC FX1000 AI・エクサスケール時代を切り拓く HPC システム https://www.fujitsu.com/jp/products/computing/servers/supercomputer/downloads/ ) 39 2023年度 計算科学技術特論A

40.

FX100とFX1000のアーキテクチャ比較 FX100 FX1000(「不老」TypeI) 演算能力/ノード 倍精度: 1.011 TFLOPS 単精度: 2.022 TFLOPS 倍精度 3.3792 TFLOPS 単精度 6.7584 TFLOPS 半精度 13.5168 TFLOPS 演算コア数 32 48 アシスタントコア 2 2 SIMD幅 256 512 SIMD命令 整数演算、ストライド& 間接ロード/ストアを 強化 ・8/16/32ビット整数演算 ・Combined Gather (128バイト アラインブロック単位) L1Dキャッシュ/コア 64KB、4ウェイ 64KB、4ウェイ L2キャッシュ/ノード 24MB 32MB, 16ウェイ/CMG メモリバンド幅 480GB/秒 1024GB/秒 出典:https://www.ssken.gr.jp/MAINSITE/event/2015/20151028-sci/lecture-04/SSKEN_sci2015_miyoshi_presentation.pdf https://monoist.atmarkit.co.jp/mn/articles/1905/07/news013.html 40 2023年度 計算科学技術特論A

41.

演算パイプライン 演算の流れ作業 41 2023年度 計算科学技術特論A

42.

流れ作業   車を作る場合 1人の作業員1つの工程を担当(5名) フロント・バッ クガラスを つける 車体作成  外装 機能確認 上記工程が2ヶ月だとする(各工程は0.4ヶ月とする)    1台目 2台目 3台目 42 内装 2ヶ月後に1台できる 4ヶ月後に2台できる 2ヶ月/台 の効率 車体作成 フロント・ バックガ ラスをつ ける 内装 外装 • 各工程の作業員は、 0.4ヶ月働いて、 1.6ヶ月は休んでいる (=作業効率が低い) 機能確 認 車体作成 フロント・ バックガ ラスをつ ける 内装 外装 機能確 認 車体作成 フロント・ バックガ ラスをつ ける 内装 外装 2023年度 計算科学技術特論A 機能確 認 時間

43.

流れ作業   作業場所は、5ヶ所とれるとする 前の工程からくる車を待ち、担当工程が終わったら、 次の工程に速やかに送られるとする  ベルトコンベア 0.4ヶ月 車体作成 43 0.4ヶ月 フロント・バック ガラスをつける 0.4か月 内装 0.4か月 外装 2023年度 計算科学技術特論A 0.4か月 機能確認

44.

流れ作業  この方法では        1台目 2台目 3台目 4台目 5台目 44 2ヶ月後に、1台できる 2.4ヶ月後に、2台できる 2.8ヶ月後に、3台できる 3.2ヶ月後に、4台できる 3.4ヶ月後に、5台できる 3.8ヶ月後に、6台できる 0.63ヶ月/台 の効率 車体作成 フロント・ バックガ ラスをつ ける 車体作成 内装 外装 機能確 認 フロント・ バックガ ラスをつ ける 内装 外装 機能確 認 車体作成 フロント・ バックガ ラスをつ ける 内装 外装 車体作成 フロント・ バックガ ラスをつ ける 車体作成 機能確 認 内装 外装 機能確 認 フロント・ バックガ ラスをつ ける 内装 外装 機能確 認 •各作業員は、 十分に時間が立つと 0.4か月の単位時間あたり 休むことなく働いている (=作業効率が高い) •このような処理を、 <パイプライン処理> という 時間 2023年度 計算科学技術特論A

45.

計算機におけるパイプライン処理の形態 ハードウエア・パイプライニング 1. 計算機ハードウエアで行う 以下の形態が代表的   1. 2. ソフトウエア・パイプライニング 2. プログラムの書き方で行う 以下の形態が代表的   1. 2. 45 演算処理におけるパイプライン処理 メモリからのデータ(命令コード、データ)転送における パイプライン処理 コンパイラが行うパイプライン処理 (命令プリロード、データ・プリロード、データ・ポストストア) 人手によるコード改編によるパイプライン処理 (データ・プリロード、ループアンローリング) 2023年度 計算科学技術特論A

46.

演算器の場合  例:演算器の工程(注:実際の演算器の計算工程は異なる) データAを メモリから取る  データBを メモリから取る 演算 を行う 行列-ベクトル積の計算では for (j=0; j<n; j++) for (i=0; i<n; i++) { y[j] += A[j][i] * x[i] ; }  A[0][0]を メモリから取る 演算結果を 収納 演算器が稼働 する工程 パイプライン化しなければ以下のようになり無駄 x[0]をメモリから 取る A[0][0]* x[0] 結果 y[0]収納 A[0][1]を メモリから取る x[1]をメモリから 取る A[0][0]* x[1] 結果 y[0]収納 A[0][2]を メモリから取る 46 2023年度 計算科学技術特論A 時間 x[2]をメモリから 取る

47.

演算器の場合 これでは演算器は、4単位時間のうち、1単位時間しか 使われていないので無駄(=演算効率1/4=25%) 以下のようなパイプライン処理ができれば、 十分時間が経つと、毎単位時間で演算がなされる (=演算効率100%) 十分な時間とは、十分な   A[0][0]を メモリから取る x[0]をメモリから 取る A[0][1]を メモリから取る A[0][0]* x[0] 結果 y[0]収納 x[1]をメモリから 取る A[0][0]* x[1] 結果 y[0]収納 x[2]をメモリから 取る A[0][2]* x[2] A[0][2]を メモリから取る A[0][3]を メモリから取る x[3]をメモリから 取る A[0][4]を メモリから取る 結果 y[0]収納 ループ反復回数があること。 行列サイズNが大きいほど、 パイプラインが滞りなく流れ、 演算効率は良くなる。 →Nが小さいと演算効率 が悪い A[0][3]* x[3] 結果 y[0]収納 x[4]をメモリから 取る A[0][2]* x[4] … 47 2023年度 計算科学技術特論A 結果 y[0]収納 時間

48.

演算パイプラインのまとめ    演算器をフル稼働させるため(=高性能計算するため) に必要な概念 メインメモリからデータを取ってくる時間はとても大きい。 演算パイプラインをうまく組めば、メモリからデータを 取ってくる時間を<隠ぺい>できる (=毎単位時間、演算器が稼働した状態にできる) 実際は以下の要因があるので、そう簡単ではない 1. 2. 3. 4. 48 計算機アーキテクチャの構成による遅延(レジスタ数の制約、 メモリ→CPU・CPU→メモリへのデータ供給量制限、など)。 ループに必要な処理(ループ導入変数(i, j)の初期化と加算処理、 ループ終了判定処理) 配列データを参照するためのメモリアドレスの計算処理 コンパイラが正しくパイプライン化される命令を生成するか 2023年度 計算科学技術特論A

49.

実際のプロセッサの場合  実際のプロセッサでは 1. 2. 加減算 乗算 ごとに独立したパイプラインがある。  さらに、同時にパイプラインに流せる命令 (同時発行命令)が複数ある。  Intel Pentium4ではパイプライン段数が31段  演算器がフル稼働になるまでの時間が長い。  分岐命令、命令発行予測ミスなど、パイプラインを中断させる処理が多発 すると、演算効率がきわめて悪くなる。 近年の周波数の低い(低電力な)マルチコアCPU/メニーコアCPUでは、 パイプライン段数が少なくなりつつある(Xeon Phiは7段)  49 2023年度 計算科学技術特論A

50.

ループ内連続アクセス 50 2023年度 計算科学技術特論A

51.

単体最適化のポイント  配列のデータ格納方式を考慮して、連続アクセスすると速い (ループ内連続アクセス) NG  for (i=0; i<n; i++) { a[ i ][1] = b[ i ] * c[ i ]; } for (i=0; i<n; i++) { a[1][ i ] = b[ i ] * c[ i ]; OK } ループを細切れにし、データアクセス範囲をキャッシュ容量内 に収めると速い(ただしnが大きいとき)(キャッシュブロック化) NG 51 for (i=0; i<n; i++) { for (j=0; j<n; j++) { a[ i ][ j ] = b[ j ] * c[ j ]; }} for (jb=0; jb<n; jb+=m) for (i=0; i<n; i++) { for (j=jb; j<jb+m; j++) { OK a[ i ][ j ] = b[ j ] * c[ j ]; }}} 2023年度 計算科学技術特論A

52.

言語に依存した配列の格納方式の違い  Fortran言語の場合  C言語の場合 A(i, j) A[i][j] j j 1 2 3 4 1 5 9 13 5 6 7 8 2 6 10 14 9 10 11 12 3 7 11 15 13 14 15 16 4 8 12 16 格納方向 i i 52 格納方向 2023年度 計算科学技術特論A

53.

行列積コード例(C言語)  コード例 for (i=0; i<n; i++) for (j=0; j<n; j++) for (k=0; k<n; k++) C[i][j] += A[i][k] *B[k][j]; j k C i A i 53 j B k 2023年度 計算科学技術特論A

54.

行列の積 n  行列積 cij =  aik bkj (i, j = 1, 2, ..., n) k =1 の実装法は、次の二通りが知られている: ループ交換法 1.  ブロック化(タイリング)法 2.  54 連続アクセスの方向を変える目的で、行列-行列 積を実現する3重ループの順番を交換する キャッシュにあるデータを再利用する目的で、 あるまとまった行列の部分データを、何度も アクセスするように実装する 2023年度 計算科学技術特論A

55.

行列の積  ループ交換法  行列積のコードは、以下のような3重ループになる(C言語) for(i=0; i<n; i++) { for(j=0; j<n; j++) { for(k=0; k<n; k++) { c[ i ][ j ] = c[ i ][ j ] + a[ i ][ k ] * b[ k][ j ]; } } }  最内部の演算は、外側の3ループを交換しても、 計算結果が変わらない → 6通りの実現の方法がある 55 2023年度 計算科学技術特論A

56.

行列の積  ループ交換法  行列積のコードは、以下のような3重ループになる(Fortran言語) do i=1,n do j=1, n do k=1, n c( i , j ) = c( i, j) + a( i , k ) * b( k , j ) enddo enddo enddo  最内部の演算は、外側の3ループを交換しても、 計算結果が変わらない → 6通りの実現の方法がある 56 2023年度 計算科学技術特論A

57.

行列の積  行列データへのアクセスパターンから、 以下の3種類に分類できる 1. 内積形式 (inner-product form) 最内ループのアクセスパタンが <ベクトルの内積>と同等 2. 外積形式 (outer-product form) 最内ループのアクセスパタンが <ベクトルの外積>と同等 3. 中間積形式 (middle-product form) 内積と外積の中間 57 2023年度 計算科学技術特論A

58.
[beta]
行列の積


内積形式 (inner-product form)



ijk, jikループによる実現(C言語)
for (i=0; i<n; i++) {
for (j=0; j<n; j++) {
dc = 0.0;
for (k=0; k<n; k++) {
dc = dc + A[ i ][ k ] * B[ k ][ j ];
}
C[ i ][ j ]= dc;
}
}

※以降、最外のループからの変数の順番で実装法
を呼ぶ。たとえば上記のコードは<ijkループ>。
58

A

B
….

●行方向と列方向のアクセスあり
→行方向・列方向格納言語の
両方で性能低下要因
解決法:
A, Bどちらか一方を転置しておく
(ただし、データ構造の変更ができ
る場合)

2023年度 計算科学技術特論A

59.

行列の積  内積形式 (inner-product form)  ijk, jikループによる実現(Fortran言語)  do i=1, n do j=1, n dc = 0.0d0 do k=1, n dc = dc + A( i , k ) * B( k , j ) enddo C( i , j ) = dc enddo enddo ※以降、最外のループからの変数の順番で実装法 を呼ぶ。たとえば上記のコードは<ijkループ>。 59 A B …. ●行方向と列方向のアクセスあり →行方向・列方向格納言語の 両方で性能低下要因 解決法: A, Bどちらか一方を転置しておく (ただし、データ構造の変更ができ る場合) 2023年度 計算科学技術特論A

60.
[beta]
行列の積


外積形式 (outer-product form)



kij, kjiループによる実現(C言語)
for (i=0; i<n; i++) {
A
B
for (j=0; j<n; j++) {
C[ i ][ j ] = 0.0;
….
}
}
for (k=0; k<n; k++) {
for (j=0; j<n; j++) {
db = B[ k ][ j ];
●kjiループでは
for (i=0; i<n; i++) {
列方向アクセスがメイン
→列方向格納言語向き
C[ i ][ j ]= C[ i ][ j ]+ A[ i ][ k ]* db;
(Fortran言語)
}
}
}
60
2023年度 計算科学技術特論A

61.

行列の積  外積形式 (outer-product form) kij, kjiループによる実現(Fortran言語)  do i=1, n A B do j=1, n C( i , j ) = 0.0d0 …. enddo enddo do k=1, n do j=1, n db = B( k , j ) ●kjiループでは do i=1, n 列方向アクセスがメイン →列方向格納言語向き C( i , j ) = C( i , j )+ A( i , k ) * db (Fortran言語) enddo enddo 61 enddo 2023年度 計算科学技術特論A 

62.
[beta]
行列の積


中間積形式 (middle-product form)
ikj, jkiループによる実現(C言語)
for (j=0; j<n; j++) {
for (i=0; i<n; i++) {
C[ i ][ j ] = 0.0;
}
for (k=0; k<n; k++) {
db = B[ k ][ j ];
for (i=0; i<n; i++) {
C[ i ][ j ] = C[ i ][ j ] + A[ i ][ k ] * db;
}
}




}
62

A

B
.
.

●jkiループでは

全て列方向アクセス
→列方向格納言語に
最も向いている
(Fortran言語)

2023年度 計算科学技術特論A

63.

行列の積  中間積形式 (middle-product form)   ikj, jkiループによる実現(Fortran言語) do j=1, n do i=1, n C( i , j ) = 0.0d0 enddo do k=1, n db = B( k , j ) do i=1, n C( i , j ) = C( i , j ) + A( i , k ) * db enddo enddo enddo 63 A B . . ●jkiループでは 全て列方向アクセス →列方向格納言語に 最も向いている (Fortran言語) 2023年度 計算科学技術特論A

64.

ループアンローリング 64 2023年度 計算科学技術特論A

65.

ループアンローリング コンパイラが、 1. レジスタへのデータの割り当て; 2. パイプライニング; がよりできるようにするため、コードを書き 換えるチューニング技法  ループの刻み幅を、1ではなく、mにする  <m段アンローリング>とよぶ  65 2023年度 計算科学技術特論A

66.

ループアンローリングの例 (行列-行列積、C言語)  k-ループ2段展開 (nが2で割り切れる場合) for (i=0; i<n; i++) for (j=0; j<n; j++) for (k=0; k<n; k+=2) C[i][j] += A[i][k] *B[k][ j] + A[i][k+1]*B[k+1][ j];  k-ループのループ判定回数が1/2になる。 66 2023年度 計算科学技術特論A

67.

ループアンローリングの例 (行列-行列積、C言語)  j-ループ2段展開 (nが2で割り切れる場合) for (i=0; i<n; i++) for (j=0; j<n; j+=2) for (k=0; k<n; k++) { C[i][ j ] += A[i][k] *B[k][ j ]; C[i][ j+1] += A[i][k] *B[k][ j+1]; }  A[i][k]をレジスタに置き、高速にアクセスできるようになる。 一般に:演算式が増えることで、ビット幅が大きなSIMD化ができる ループ中の式が少ない場合、アンローリングして増やさないとSIMD化できない 67 2023年度 計算科学技術特論A

68.

ループアンローリングの例 (行列-行列積、C言語)  i-ループ2段展開 (nが2で割り切れる場合) for (i=0; i<n; i+=2) for (j=0; j<n; j++) for (k=0; k<n; k++) { C[i ][j] += A[i ][k] *B[k][j]; C[i+1][j] += A[i+1][k] *B[k][j]; }  B[i][j]をレジスタに置き、高速にアクセスできるようになる。 一般に:演算式が増えることで、ビット幅が大きなSIMD化ができる 68 2023年度 計算科学技術特論A

69.

ループアンローリングの例 (行列-行列積、C言語)  i-ループ、および j-ループ 2段展開 (nが2で割り切れる場合) for (i=0; i<n; i+=2) for (j=0; j<n; j+=2) for (k=0; k<n; k++) { C[i ][ j ] += A[i ][k] *B[k][ j ]; C[i ][ j+1] += A[i ][k] *B[k][ j+1]; C[i+1][ j ] += A[i+1][k] *B[k][ j ]; C[i+1][ j+1] += A[i+1][k] *B[k][ j +1]; }  69 A[i][j], A[i+1][k],B[k][j],B[k][j+1]をレジスタに置き、 高速にアクセスできるようになる。 2023年度 計算科学技術特論A

70.

ループアンローリングの例 (行列-行列積、C言語)  コンパイラにわからせるため、以下のように書く方がよい 場合がある  for (i=0; i<n; i+=2) for (j=0; j<n; j+=2) { dc00 = C[i ][ j ]; dc01 = C[i ][ j+1]; dc10 = C[i+1][ j ]; dc11 = C[i+1][ j+1] ; for (k=0; k<n; k++) { da0= A[i ][k] ; da1= A[i+1][k] ; db0= B[k][ j ]; db1= B[k][ j+1]; dc00 += da0 *db0; dc01 += da0 *db1; dc10 += da1 *db0; dc11 += da1 *db1; } C[i ][ j ] = dc00; C[i ][ j+1] = dc01; C[i+1][ j ] = dc10; C[i+1][ j+1] = dc11; } 70 2023年度 計算科学技術特論A

71.

ループアンローリングの例 (行列-行列積、Fortran言語)  k-ループ2段展開 (nが2で割り切れる場合) do i=1, n do j=1, n do k=1, n, 2 C(i, j) = C(i, j) +A(i, k) *B(k, j) + A(i, k+1)*B(k+1, j) enddo enddo enddo  k-ループのループ判定回数が1/2になる。 71 2023年度 計算科学技術特論A

72.

ループアンローリングの例 (行列-行列積、Fortran言語)  j-ループ2段展開 (nが2で割り切れる場合) do i=1, n do j=1, n, 2 do k=1, n C(i, j ) = C(i, j ) +A(i, k) * B(k, j ) C(i, j+1) = C(i, j+1) +A(i, k) * B(k, j+1) enddo enddo enddo  A(i, k)をレジスタに置き、高速にアクセスできるようになる。 72 2023年度 計算科学技術特論A

73.

ループアンローリングの例 (行列-行列積、Fortran言語)  i-ループ2段展開 (nが2で割り切れる場合) do i=1, n, 2 do j=1, n do k=1, n C(i , j) = C(i , j) +A(i , k) * B(k , j) C(i+1, j) = C(i+1, j) +A(i+1, k) * B(k , j) enddo enddo enddo  73 B(i, j)をレジスタに置き、高速にアクセスできるようになる。 2023年度 計算科学技術特論A

74.

ループアンローリングの例 (行列-行列積、Fortran言語)  i-ループ、および j-ループ 2段展開 (nが2で割り切れる場合) do i=1, n, 2 do j=1, n, 2 do k=1, n C(i , j ) = C(i , j ) +A(i , k) *B(k, j ) C(i , j+1) = C(i , j+1) +A(i , k) *B(k, j+1) C(i+1, j ) = C(i+1, j ) +A(i+1, k) *B(k, j ) C(i+1, j+1) =C(i+1, j+1) +A(i+1, k) *B(k, j +1) enddo; enddo; enddo;  A(i,j), A(i+1,k),B(k,j),B(k,j+1)をレジスタに置き、 高速にアクセスできるようになる。 74 2023年度 計算科学技術特論A

75.

ループアンローリング と レジスタスピル  レジスタにデータが載るようになり高速化される  ループアンローリング段数を増やす  レジスタ数の上限に達する  メモリにデータを書き戻す(レジスタスピル)⇒速度低下 一般に:ループ内の式が多いと、レジスタスピルが起きやすい ハードウェア上のレジスタ数が少ない場合も、レジスタスピルが起きやすい この場合、ループ分割 をして式を減らして レジスタスピルを防ぎ 高速化 75 2023年度 計算科学技術特論A

76.

ループアンローリングの例 (行列-行列積、Fortran言語)  コンパイラにわからせるため、以下のように書く方がよい 場合がある  do i=1, n, 2 do j=1, n, 2 dc00 = C(i ,j ); dc01 = C(i ,j+1) dc10 = C(i+1,j ); dc11 = C(i+1,j+1) do k=1, n da0= A(i ,k); da1= A(i+1, k) db0= B(k ,j ); db1= B(k, j+1) dc00 = dc00+da0 *db0; dc01 = dc01+da0 *db1; dc10 = dc10+da1 *db0; dc11 = dc11+da1 *db1; enddo C(i , j ) = dc00; C(i , j+1) = dc01 C(i+1, j ) = dc10; C(i+1, j+1) = dc11 enddo; enddo 76 2023年度 計算科学技術特論A

77.

キャッシュライン衝突 とびとびアクセスは弱い 77 2023年度 計算科学技術特論A

78.

不連続アクセスとは   C言語の場合 配列のデータ格納方式を考慮し 連続アクセスすると速い (ループ内連続アクセス) NG a[i][j] j for (i=0; i<n; i++) { a[ i ][1] = b[ i ] * c[ i ]; } 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 格納方向 i 間隔4での不連続アクセス 78 2023年度 計算科学技術特論A

79.

キャッシュメモリの構成 キャッシュライン CPU キャッシュメモリ (キャッシュ上のブロック) 演算器 演算 要求 10 演算 結果 0 6 2 14 レジスタ データ供給 データ供給 メインメモリ セット データ蓄積 (ブロック の並び) ブロック メインメモリ (記憶単位) 注)配列をアクセスすると、1要素分ではなく ブロック単位のデータ(例えば32バイト(倍精度4変数分) が同時にキャッシュに乗る(ブロックサイズと呼ぶ) 79 写像関数 データ蓄積 キャッシュメモリ ブロックと キャッシュラインの 対応 8 9 10 11 12 13 14 0 1 2 2023年度 計算科学技術特論A 3 4 6 7

80.

キャッシュとキャッシュライン  メインメモリ上とキャッシュ上のデータマッピング方式  読み出し: メインメモリ から キャッシュ へ    ダイレクト・マッピング方式: メモリバンクごとに直接的 セット・アソシアティブ方式: ハッシュ関数で写像(間接的) 書き込み: キャッシュ から メインメモリ へ   ストア・スルー方式: キャッシュ書き込み時に メインメモリと中身を一致させる ストア・イン方式: 対象となるキャッシュラインが 置き換え対象となったときに一致させる キャッシュメモリ キャッシュ ライン ライン0 ライン1 ライン2 ライン3 ライン4 ライン5 … 80 メインメモリ 写像関数 メモリブロック … 2023年度 計算科学技術特論A

81.

キャッシュライン衝突の例 直接メインメモリのアドレスをキャッシュに写像する、ダイレクト・マッピングを考える   物理結線は以下の通り マッピング間隔を、ここでは4とする   メインメモリ上のデータは、間隔4ごとに、同じキャッシュラインに乗る キャッシュラインは8バイト、メモリバンクも8バイトとする 配列aは 4×4の構成で、倍精度(8バイト)でメモリ確保されているとする   double a[4][4];  この前提で、格納方向と逆方向にアクセス(4とびのアクセス)する (=C言語の場合、i方向を連続アクセス) メインメモリ キャッシュ キャッシュメモリ 1 2 3 ライン ライン0 5 6 7 ライン1 9 10 11 ライン2 13 14 15 ライン3 物理結線 配列アクセス方向 81 2023年度 計算科学技術特論A … メモリ 連続方向 4 8 12 16

82.

キャッシュライン衝突の例  この前提の、<実際の配列構成>と<メモリブロック>の関係 実際は、以下のことがあるので、必ずしも、こうならないことに注意する  配列a[][]の物理メモリ上の配置はOSが動的に決定するので、ずれることがある  メモリブロックの容量は、8バイトより大きい  ダイレクト・マッピングではない  C言語の場合 配列a[i][j] メインメモリ上の ブロック構成 j 配列要素a[][] と メモリブロック構造と が完全一致 1 2 3 4 5 6 7 8 12 9 10 11 12 16 13 14 15 16 1 2 3 4 5 6 7 8 9 10 11 13 14 15 格納方向 i 82 2023年度 計算科学技術特論A …

83.

キャッシュライン衝突の例 1. a[0][0]があるブロック1がキャッシュライン0に乗る 2. すぐに、a[1][0]があるブロック5がアクセスされる 3. (物理結線先のキャッシュライン0に容量の空きがないので) キャッシュライン0のデータ(ブロック1の内容)を追い出さないといけない 4. ブロック5のデータがキャッシュライン0に乗る 5. すぐに、a[2][0]があるブロック9がアクセスされる 6. キャッシュライン0のデータ(ブロック5の内容)を追い出さないといけない …玉突きで、ライン1~3が空いていても、逐次的にキャッシュ上のデータが 追い出される キャッシュ ライン レジスタへ 83 キャッシュメモリ メモリ連続 メインメモリ 1 5 9 13 9 5 1 ライン0 ライン1 ライン2 ライン3 2 6 10 14 配列アクセス方向 2023年度 計算科学技術特論A … 3 7 11 15 4 8 12 16

84.

キャッシュライン衝突の例  1~6の状態が連続して発生する。 メモリ→キャッシュの回線が常に稼働    <回線お話し中>で、データが来るのが終わるまで、待たされる (回線レベルで並列にデータが持ってこれない) ストア・イン方式では、メモリにデータを書き戻すコストもかかる メモリからデータを逐次で読み出すのと同じ <キャッシュがない>のと同じ 演算器にデータが届かないので計算を中断。 演算器の利用効率が悪くなる 以上の現象を<キャッシュライン衝突>と呼ぶ 84 2023年度 計算科学技術特論A

85.

メモリ・インターリービング  物理的なメモリの格納方向に従いアクセスする時  データアクセス時、現在アクセス中のブロック上のデータは、 周辺ブロック上のデータも一括して(同時に)、別の キャッシュライン上に乗せるハードウェア機能がある キャッシュライン0のデータをアクセスしている最中に、 キャッシュライン1に近隣のブロック内データを(並列に) 持ってくることが可能 メモリの<インタリービング> 演算機から見たデータアクセス時間が短縮 演算器が待つ時間が減少(=演算効率が上がる) 物理的なデータ格納方向に連続アクセスするとよい 85 2023年度 計算科学技術特論A

86.

キャッシュライン衝突が起こる条件  メモリバンクのキャッシュラインへの割り付けは 2冪の間隔で行っていることが多い たとえば、32、64、128など 特定サイズの問題(たとえば1024次元)で、 性能が1/2~1/3、ときには1/10になる 場合、キャッシュライン衝突が生じている可能性あり   double a[1024][1024]; NG double precision a(1024, 1024) 実際は、OSやキャッシュ構成の影響で厳密な条件を見つけることは難しいが 2冪サイズでの配列確保は避けるべき 86 2023年度 計算科学技術特論A

87.

キャッシュライン衝突への対応  キャッシュライン衝突を防ぐ方法 パティング法: 配列に(2冪でない)余分な領域を確保 し確保配列の一部の領域を使う。 1.   2. 3. 87 余分な領域を確保して使う  例: double A[1024][1025]; で1024のサイズをアクセス コンパイラのオプションを使う データ圧縮法: 計算に必要なデータのみキャッシュ ライン衝突しないようにデータを確保し、かつ、必要な データをコピーする。 予測計算法: キャッシュライン衝突が起こる回数を 予測するルーチンを埋め込み、そのルーチンを配列 確保時に呼ぶ。 2023年度 計算科学技術特論A

88.

ブロック化 小さい範囲のデータ再利用 88 2023年度 計算科学技術特論A

89.

ブロック化によるアクセス局所化     89 キャッシュには大きさがあります。 この大きさを超えると、たとえ連続アクセスしても、 キャッシュからデータは追い出されます。 データが連続してキャッシュから追い出されると、 メモリから転送するのと同じとなり、高速な アクセス速度を誇るキャッシュの恩恵がなくなります。 そこで、高速化のためには、以下が必要です 1. キャッシュサイズ限界までデータを詰め込む 2. 詰め込んだキャッシュ上のデータを、何度も アクセスして再利用する 2023年度 計算科学技術特論A

90.

ブロック化によるキャッシュミスヒット 削減例 行列ー行列積  行列サイズ:8×8   double A[8][8]; キャッシュラインは4つ  1つのキャッシュラインに4つの行列要素が載る   キャッシュライン:4×8バイト(double)=32バイト 配列の連続アクセスは行方向(C言語)  キャッシュの追い出しアルゴリズム: Least Recently Used (LRU)  90 2023年度 計算科学技術特論A

91.

配列とキャッシュライン構成の関係  この前提の、<配列構成>と<キャッシュライン>の関係  ここでは、キャッシュライン衝突は考えません  C言語の場合 配列A[i][j]、B[i][j]、C[i][j] j 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 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64  1×4の配列要素が、 キャッシュラインに乗る  どのキャッシュラインに 乗るかは、<配列アクセス パターン> と <置き換え アルゴリズム>依存で決まる キャッシュラインの 構成 1 3 格納方向 i 91 2023年度 計算科学技術特論A 2 4

92.

行列-行列積の場合(ブロック化しない) LRU:直近で最もアクセス されていないラインの データを追い出す キャッシュミス① キャッシュミス② キャッシュミス③ キャッシュミス④ = C * キャッシュミス⑤ A B キャッシュライン ライン1 ライン2 ライン3 ライン4 92 ※キャッシュライン4つ、 置き換えアルゴリズム LRUの場合 2023年度 計算科学技術特論A

93.

行列-行列積の場合(ブロック化しない) キャッシュミス11 キャッシュミス⑥ = * キャッシュミス⑦ キャッシュミス⑧ キャッシュミス⑨ C キャッシュミス⑩ A B キャッシュライン ライン1 ライン2 ライン3 ライン4 93 ※キャッシュライン4つ、 置き換えアルゴリズム LRUの場合 2023年度 計算科学技術特論A

94.

行列-行列積の場合(ブロック化しない) ライン1 ライン2 ライン3 ライン4 キャッシュミス キャッシュミス キャッシュミス キャッシュミス キャッシュミス キャッシュミス = キャッシュミス キャッシュミス * キャッシュミス キャッシュミス キャッシュミス C A B キャッシュライン ※2要素計算するのに、 キャッシュミスヒット22回 94 ※キャッシュライン4つ、 置き換えアルゴリズム LRUの場合 2023年度 計算科学技術特論A

95.

行列-行列積の場合(ブロック化する:2要素) ライン1 ライン2 ライン3 ライン4 このブロック幅 単位で計算する キャッシュミス 12 キャッシュミス キャッシュミス 1 ① ② = C ① キャッシュミス キャッシュミス ② * A キャッシュミス B キャッシュライン ※キャッシュライン4つ、 置き換えアルゴリズム LRUの場合 95 2023年度 計算科学技術特論A

96.

行列-行列積の場合(ブロック化する:2要素) ライン1 ライン2 ライン3 ライン4 このブロック幅 単位で計算する 12 キャッシュミス キャッシュミス 1 ③ = ④ * キャッシュミス ③ キャッシュミス キャッシュミス ④ C キャッシュライン 96 キャッシュミス A ※2要素計算するのに、 キャッシュミスヒット10回 B ※キャッシュライン4つ、 置き換えアルゴリズム LRUの場合 2023年度 計算科学技術特論A

97.

行列積コード(C言語) :キャッシュブロック化なし  コード例 for (i=0; i<n; i++) for (j=0; j<n; j++) for (k=0; k<n; k++) C[i][j] += A[i][k] *B[k][j]; j k C i A i 97 j B k 2023年度 計算科学技術特論A

98.

行列-行列積のブロック化のコード (C言語)  nがブロック幅(ibl=16)で割り切れるとき、 以下のような6重ループのコードになる ibl = 16; for ( ib=0; ib<n; ib+=ibl ) { for ( jb=0; jb<n; jb+=ibl ) { for ( kb=0; kb<n; kb+=ibl ) { for ( i=ib; i<ib+ibl; i++ ) { for ( j=jb; j<jb+ibl; j++ ) { for ( k=kb; k<kb+ibl; k++ ) { C[i][j] += A[i][k] * B[k][j]; } } } } } } 98 2023年度 計算科学技術特論A

99.

行列-行列積のブロック化のコード (Fortran言語)  nがブロック幅(ibl=16)で割り切れるとき、 以下のような6重ループのコードになる ibl = 16 do ib=1, n, ibl do jb=1, n, ibl do kb=1, n, ibl do i=ib, ib+ibl-1 do j=jb, jb+ibl-1 do k=kb, kb+ibl-1 C(i, j) = C(i, j) + A(i, k) * B(k, j) enddo; enddo; enddo; enddo; enddo; enddo; 99 2023年度 計算科学技術特論A

100.

キャッシュブロック化時の データ・アクセスパターン ibl×iblの 小行列単位で 行列‐行列積 をする ibl ibl ibl ibl ibl ibl × = C 100 A 2023年度 計算科学技術特論A B

101.

キャッシュブロック化時の データ・アクセスパターン ibl×iblの 小行列単位で 行列‐行列積 をする ibl ibl ibl ibl ibl × = ibl C 101 A 2023年度 計算科学技術特論A B

102.

行列-行列積のブロック化のコードの アンローリング(C言語)   行列-行列積の6重ループのコードに加え、 さらに各6重ループにアンローリングを施すことができる。 i-ループ、およびj-ループ2段アンローリングは、以下のよ うなコードになる。(ブロック幅iblが2で割り切れる場合) ibl = 16; for (ib=0; ib<n; ib+=ibl) { for (jb=0; jb<n; jb+=ibl) { for (kb=0; kb<n; kb+=ibl) { for (i=ib; i<ib+ibl; i+=2) { for (j=jb; j<jb+ibl; j+=2) { for (k=kb; k<kb+ibl; k++) { C[i ][j ] += A[i ][k] * B[k][j ]; C[i+1][j ] += A[i+1][k] * B[k][j ]; C[i ][j+1] += A[i ][k] * B[k][j+1]; C[i+1][j+1] += A[i+1][k] * B[k][j+1]; } } } } } } 102 2023年度 計算科学技術特論A

103.

行列-行列積のブロック化のコードの アンローリング(Fortran言語)   行列-行列積の6重ループのコードに加え、 さらに各6重ループにアンローリングを施すことができる。 i-ループ、およびj-ループ2段アンローリングは、以下のよ うなコードになる。(ブロック幅iblが2で割り切れる場合) ibl = 16 do ib=1, n, ibl do jb=1, n, ibl do kb=1, n, ibl do i=ib, ib+ibl, 2 do j=jb, jb+ibl, 2 do k=kb, kb+ibl C(i , j ) = C(i , j ) + A(i , k) * B(k, j ) C(i+1, j ) = C(i+1, j ) + A(i+1, k) * B(k, j ) C(i , j+1) = C(i , j+1) + A(i , k) * B(k, j+1) C(i+1, j+1) = C(i+1, j+1) + A(i+1, k) * B(k, j+1) enddo; enddo; enddo; enddo; enddo; enddo; 103 2023年度 計算科学技術特論A

104.

その他の高速化技術 104 2023年度 計算科学技術特論A

105.

共通部分式の削除(1)  以下のプログラムは、冗長な部分がある。 d = a + b + c; f = d + a + b;  コンパイラがやる場合もあるが、以下のように書く方が 無難である。 temp = a + b; d = temp + c; f = d + temp; 105 2023年度 計算科学技術特論A

106.

共通部分式の削除(2)  配列のアクセスも、冗長な書き方をしないほうがよい。 for (i=0; i<n; i++) { xold[i] = x[i]; x[i] = x[i] + y[i]; }  以下のように書く。 for (i=0; i<n; i++) { dtemp = x[i]; xold[i] = dtemp; x[i] = dtemp + y[i]; } 106 2023年度 計算科学技術特論A

107.

コードの移動  割り算は演算時間がかかる。ループ中に書かない。 for (i=0; i<n; i++) { a[i] = a[i] / sqrt(dnorm); }  上記の例では、掛け算化して書く。 dtemp = 1.0 / sqrt(dnorm); for (i=0; i<n; i++) { a[i] = a[i] *dtemp; } 107 2023年度 計算科学技術特論A

108.

ループ中のIF文  なるべく、ループ中にIF文を書かない。 for (i=0; i<n; i++) { for (j=0; j<n; j++) { if ( i != j ) A[i][j] = B[i][j]; else A[i][j] = 1.0; } }  以下のように書く。 for (i=0; i<n; i++) { for (j=0; j<n; j++) { A[i][j] = B[i][j]; } } for (i=0; i<n; i++) A[i][i] = 1.0; 108 2023年度 計算科学技術特論A

109.

ソフトウェア・パイプライニングの強化  基のコード (2段のアンローリング) 定義-参照の距離が近い →ソフトウェア的には 何もできない  ソフトウェアパイプライニング を強化したコード (2段のアンローリング) 定義-参照の距離が遠い →ソフトウェアパイプライニング が適用できる機会が増加! 109 for (i=0; i<n; i+=2) { dtmpb0 = b[i]; dtmpc0 = c[i]; dtmpa0 = dtmpb0 + dtmpc0; a[i] = dtmpa0; dtmpb1 = b[i+1]; dtmpc1 = c[i+1]; dtmpa1 = dtmpb1 + dtmpc1; a[i+1] = dtmpa1; } for (i=0; i<n; i+=2) { dtmpb0 = b[i]; dtmpb1 = b[i+1]; dtmpc0 = c[i]; dtmpc1 = c[i+1]; dtmpa0 = dtmpb0 + dtmpc0; dtmpa1 = dtmpb1 + dtmpc1; a[i] = dtmpa0; a[i+1] = dtmpa1; } 2023年度 計算科学技術特論A

110.

数値計算ライブラリの利用 110 2023年度 計算科学技術特論A

111.

数値計算ライブラリ  密行列用ライブラリ      行列の要素に0がない(というデータ構造を扱う) 連立一次方程式の解法、固有値問題、FFT、その他 直接解法(反復解法もある) BLAS、LAPACK、ScaLAPACK、SuperLU、MUMPS、FFTW、など 疎行列用ライブラリ     行列の要素に0が多い 連立一次方程式の解法、固有値問題、その他 反復解法 PETSc、Xabclib、Lis、ARPACK、など 111 2023年度 計算科学技術特論A

112.

疎行列用ライブラリの特徴  疎行列を扱うアプリケーションはライブラリ化が難しい  疎行列データ形式の標準化が困難   カーネルの演算が微妙に違う、かつ、カーネルは広い範囲に分散   COO、CRS(CCS)、ELL、JDS、BCSR、・・・ 陽解法(差分法)を基にしたソフトウェア 数値ミドルウェアおよび領域特化型言語 (Domain Specific Language, DSL)    解くべき方程式や離散化方法に特化させることで、処理(対象となるプロ グラムの性質)を限定 以上の限定から、高度な最適化ができる言語(処理系)の作成(DSL)や、 ライブラリ化(数値ミドルウェア)ができる 数値ミドルウェアの例  112 ppOpen-HPC(東大)、PETSc(Argonne National Laboratory, USA.) 、 Trilinos (Sandia National Laboratory, USA)、など 2023年度 計算科学技術特論A

113.

BLAS  BLAS(Basic Linear Algebra Subprograms、 基本線形代数副プログラム集)    線形代数計算で用いられる、基本演算を標準化 (API化)したもの。 普通は、密行列用の線形代数計算用の基本演算 の副プログラムを指す。 疎行列の基本演算用の<スパースBLAS>というも のあるが、まだ定着していない。  113 スパースBLASはIntel MKL(Math Kernel Library)に入って いるが、広く使われているとは言えない。 2023年度 計算科学技術特論A

114.

BLAS  BLASでは、以下のように分類わけをして、 サブルーチンの命名規則を統一 1. 2. 3. 4.  演算対象のベクトルや行列の型(整数型、実数型、複素型) 行列形状(対称行列、三重対角行列) データ格納形式(帯行列を二次元に圧縮) 演算結果が何か(行列、ベクトル) 演算性能から、以下の3つに演算を分類  レベル1 BLAS: ベクトルとベクトルの演算  レベル2 BLAS: 行列とベクトルの演算  レベル3 BLAS: 行列と行列の演算 114 2023年度 計算科学技術特論A

115.

レベル1 BLAS  レベル1 BLAS     115 ベクトル内積、ベクトル定数倍の加算、など  例: y ← α x + y データの読み出し回数、演算回数がほほ同じ データの再利用(キャッシュに乗ったデータの再利用による データアクセス時間の短縮)がほとんどできない  実装による性能向上が、あまり期待できない  ほとんど、計算機ハードウエアの演算性能 レベル1BLASのみで演算を実装すると、演算が本来 持っているデータ再利用性がなくなる  例:行列-ベクトル積を、レベル1BLASで実装 2023年度 計算科学技術特論A

116.

レベル2 BLAS  レベル2 BLAS  行列-ベクトル積などの演算     例: y ← α A x + β y 前進/後退代入演算、T x = y (Tは三角行列)を xについて解く演算、を含む レベル1BLASのみの実装よる、データ再利用性の喪失 を回避する目的で提案 行列とベクトルデータに対して、データの再利用性あり   116 データアクセス時間を、実装法により短縮可能 (実装法により)性能向上がレベル1BLASに比べ しやすい(が十分でない) 2023年度 計算科学技術特論A

117.

レベル3 BLAS  レベル3 BLAS  行列-行列積などの演算 例: C ← α A B + β C 共有記憶型の並列ベクトル計算機では、レベル2 BLASでも 性能向上が達成できない。       117 並列化により1PE当たりのデータ量が減少する。 より大規模な演算をとり扱わないと、再利用の効果がない。 行列-行列積では、行列データ O(n2 ) に対して 演算は O (n 3 ) なので、データ再利用性が原理的に高い。 行列積は、アルゴリズムレベルでもブロック化できる。 さらにデータの局所性を高めることができる。 2023年度 計算科学技術特論A

118.

典型的なBLASの性能 性能 [FLOPS] 理論性能の限界 BLAS3 BLAS2 BLAS1 行列サイズ 118 2023年度 計算科学技術特論A

119.

BLAS利用例  倍精度演算BLAS3 C := alpha*op( A )*op( B ) + beta*C A: M*K; B:K*N; C:M*N; CALL DGEMM( ‘N’, ‘N’, n, n, n, ALPHA, A, N, B, N, BETA, C, N ) Aが転置しているか Bが転置しているか Mの大きさ Nの大きさ alpha の値 Aの アドレス Aの1次元目 の要素数 Bの アドレス Bの1次元目 の要素数 Kの大きさ 119 beta の値 2023年度 計算科学技術特論A Cの アドレス Cの1次元目 の要素数

120.

BLASの機能詳細 詳細はHP: http://www.netlib.org/blas/  命名規則: 関数名:XYYYY    X: データ型 S:単精度、D:倍精度、C:複素、Z:倍精度複素 YYYY: 計算の種類    120 レベル1: 例:AXPY:ベクトルをスカラー倍して加算 レベル2: 例:GEMV: 一般行列とベクトルの積 レベル3: 例:GEMM:一般行列どうしの積 2023年度 計算科学技術特論A

121.

GOTO BLASとは   後藤和茂 氏により開発された、ソースコードが 無償入手可能な、高性能BLASの実装(ライブラリ) 特徴   マルチコア対応がなされている 多くのコモディティハードウエア上の実装に特化    Intel Nehalem and Atom systems VIA Nanoprocessor AMD Shanghai and Istanbul 等  テキサス大学先進計算センター(TACC)で、 GOTO BLAS2として、ソースコードを配布している  121 https://www.tacc.utexas.edu/research-development/tacc-software/gotoblas2 2023年度 計算科学技術特論A

122.

LAPACK   密行列に対する、連立一次方程式の解法、 および固有値の解法の“標準”アルゴリズムルーチンを 無償で提供 その道の大学の専門家が集結    カリフォルニア大バークレー校: James Demmel教授 テネシー大ノックスビル校: Jack Dongarra教授 HP http://www.netlib.org/lapack/ 122 2023年度 計算科学技術特論A

123.

LAPACKの命名規則  命名規則: 関数名:XYYZZZ    123 X: データ型 S:単精度、D:倍精度、C:複素、Z:倍精度複素 YY: 行列の型 BD: 二重対角、DI:対角、GB:一般帯行列、GE:一般行列、 HE:複素エルミート、HP:複素エルミート圧縮形式、SY:対称 行列、…. ZZZ: 計算の種類 TRF: 行列の分解、TRS:行列の分解を使う、CON:条件数 の計算、RFS:計算解の誤差範囲を計算、TRI:三重対角行 列の分解、EQU:スケーリングの計算、… 2023年度 計算科学技術特論A

124.

インタフェース例:DGESV (1/3)  DGESV (N, NRHS, A, LDA, IPIVOT, B, LDB, INFO)      A X = B の解の行列Xを計算をする A * X = B、ここで A はN×N行列で、 X と B は N×NRHS行列 とする。 行交換の部分枢軸選択付きのLU分解 でA を A = P * L * U と分 解する。ここで、P は交換行列、L は下三角行列、Uは上三角行列 である。 分解されたA は、連立一次方程式A * X = Bを解くのに使われる。 引数  N (入力) - INTEGER  124 線形方程式の数。行列Aの次元数。 N >= 0。 2023年度 計算科学技術特論A

125.

インタフェース例:DGESV (2/3)  NRHS (入力) – INTEGER   A (入力/出力) – DOUBLE PRECISION, DIMENSION(:,:)   入力時は、N×Nの行列Aの係数を入れる。 出力時は、Aから分解された行列LとU = P*L*Uを圧縮して出力する。 Lの対角要素は1であるので、収納されていない。 LDA (入力) – INTEGER    右辺ベクトルの数。行列Bの次元数。 NRHS >= 0。 配列Aの最初の次元の大きさ。 LDA >= max(1,N)。 IPIVOT (出力) - DOUBLE PRECISION, DIMENSION(:)  125 交換行列Aを構成する枢軸のインデックス。 行列のi行がIPIVOT(i)行と交 換されている。 2023年度 計算科学技術特論A

126.

インタフェース例:DGESV (3/3)  B (入力/出力) – DOUBLE PRECISION, DIMENSION(:,:)    LDB (入力) -INTEGER   入力時は、右辺ベクトルの N×NRHS 行列Bを入れる。 出力時は、もし、INFO = 0 なら、N×NRHS行列である解行列Xが戻る。 配列Bの最初の次元の大きさ。 LDB >= max(1,N)。 INFO (出力) ーINTEGER    126 = 0: 正常終了 < 0: もし INFO = -i なら i-th 行の引数の値がおかしい。 > 0: もし INFO = i なら U(i,i) が厳密に0である。 分解は終わるが、 Uの分解は特異なため、解は計算されない。 2023年度 計算科学技術特論A

127.

ScaLAPACK    密行列に対する、連立一次方程式の解法、 および固有値の解法の“標準”アルゴリズムルーチンの 並列化版を無償で提供 ユーザインタフェースはLAPACKに<類似> ソフトウェアの<階層化>がされている     内部ルーチンはLAPACKを利用 並列インタフェースはBLACS データ分散方式に、2次元ブロック・サイクリック分散方式 を採用 (詳細は、「MPI」の講義で説明) HP: http://www.netlib.org/scalapack/ 127 2023年度 計算科学技術特論A

128.

ScaLAPACKのソフトウェア構成図 出典:http://www.netlib.org/scalapack/poster.html 分散メモリ用 アルゴリズムのライブラリ ScaLAPACK 分散メモリ用演算カーネル ライブラリ PBLAS キャッシュ 最適化アルゴリズム のライブラリ 大域アドレス 局所アドレス LAPACK 環境独立 環境依存 BLAS 演算カーネル ライブラリ 128 ScaLAPACK用 通信ライブラリ BLACS Message Passing Interface (MPI) 2023年度 計算科学技術特論A 汎用 通信ライブラリ

129.

BLACSとPBLAS  BLACS   ScaLAPACK中で使われる通信機能を関数化したもの。 通信ライブラリは、MPI、PVM、各社が提供する通信ライブラリを 想定し、ScaLAPACK内でコード修正せずに使うことを目的とする   現在、MPIがデファクトになったため、MPIで構築された BLACSのみ、現実的に利用されている。   いわゆる、通信ライブラリのラッパー的役割でScaLAPACK内で利用 なので、ScaLAPACKはMPIでコンパイルし、起動して利用する PBLAS   129 BLACSを用いてBLASと同等な機能を提供する関数群 並列版BLASといってよい。 2023年度 計算科学技術特論A

130.

ScaLAPACKの命名規則 原則: LAPACKの関数名の頭に“P”を付けたもの  そのほか、BLACS、PBLAS、データ分散を 制御するためのScaLAPACK用関数がある。  130 2023年度 計算科学技術特論A

131.

インタフェース例:PDGESV (1/4)  PDGESV ( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO )      sub(A) X = sub(B) の解の行列Xを計算をする ここで sub(A) はN×N行列を分散したA(IA:IA+N-1, JA:JA+N-1) の行列 X と B は N×NRHS行列を分散したB(IB:IB+N-1, JB:JB+NRHS-1) の行列 行交換の部分枢軸選択付きのLU分解 でsub(A) を sub(A) = P * L * U と分解する。ここで、P は交換行列、 L は下三角行列、Uは上三角行列である。 分解されたsub(A) は、連立一次方程式sub(A) * X = sub(B)を 解くのに使われる。 131 2023年度 計算科学技術特論A

132.

インタフェース例:PDGESV (2/4)  N (大域入力) – INTEGER   NRHS (大域入力) – INTEGER      右辺ベクトルの数。行列Bの次元数。 NRHS >= 0。 A (局所入力/出力) – DOUBLE PRECISION, DIMENSION(:,:)   線形方程式の数。行列Aの次元数。 N >= 0。 入力時は、N×Nの行列Aの局所化された係数を 配列A(LLD_A, LOCc( JA+N-1))を入れる。 出力時は、Aから分解された行列LとU = P*L*Uを圧縮して出力する。 Lの対角要素は1であるので、収納されていない。 IA(大域入力) -INTEGER :sub(A)の最初の行のインデックス JA(大域入力) -INTEGER :sub(A)の最初の列のインデックス DESCA (大域かつ局所入力) – INTEGER  132 分散された配列Aの記述子。 2023年度 計算科学技術特論A

133.

インタフェース例:PDGESV (3/4)  IPIVOT (局所出力) - DOUBLE PRECISION, DIMENSION(:)   B (局所入力/出力) – DOUBLE PRECISION, DIMENSION(:,:)    sub(B)の最初の行のインデックス JB(大域入力) -INTEGER   入力時は、右辺ベクトルの N×NRHSの行列Bの分散されたものを (LLD_B, LOCc(JB+NRHS-1)) に入れる。 出力時は、もし、INFO = 0 なら、N×NRHS行列である解行列Xが、 行列Bと同様の分散された状態で戻る。 IB(大域入力) -INTEGER   交換行列Aを構成する枢軸のインデックス。 行列のi行がIPIVOT(i)行と 交換されている。分散された配列( LOCr(M_A)+MB_A )として戻る。 sub(B)の最初の列のインデックス DESCB (大域かつ局所入力) – INTEGER  133 分散された配列Bの記述子。 2023年度 計算科学技術特論A

134.

インタフェース例:PDGESV (4/4)  INFO (大域出力) ーINTEGER    134 = 0: 正常終了 < 0:  もし i番目の要素が配列で、 そのj要素の値がおかしいなら、 INFO = -(i*100+j)となる。  もしi番目の要素がスカラーで、かつ、その値がおかしいなら、 INFO = -iとなる。 > 0: もし INFO = Kのとき U(IA+K-1, JA+K-1) が厳密に0である。 分解は完了するが、分解されたUは厳密に特異なので、 解は計算できない。 2023年度 計算科学技術特論A

135.

BLAS利用の注意  C言語からの利用   BLASライブラリは(たいてい)Fortranで書かれている 行列を1次元で確保する    Fortranに対して転置行列になるので、BLASの引数で転置を指定 引数は全てポインタで引き渡す 関数名の後に“_”をつける(BLASをコンパイルするコンパイラ依存) 例:dgemm_(...) 小さい行列は性能的に注意    キャッシュに載るようなサイズ(例えば100次元以下)の行列については、 BLASが高速であるとは限らない   全体の行列サイズは大きくても、利用スレッド数が多くなると、 スレッド当たりの行列サイズが小さくなるので注意!  135 BLASは、大規模行列で高性能になるように設計されている 例) N=8000でも272スレッド並列だと、スレッドあたり約480x480 まで小さくなる 2023年度 計算科学技術特論A

136.

その他のライブラリ(主に行列演算) 種類 問題 ライブラリ名 概要 密行列 BLAS MAGMA GPU、マルチコア、ヘテロジニ アス環境対応 疎行列 連立一次方程式 MUMPS 直接解法 SuperLU 直接解法 PETSc 反復解法、各種機能 Hypre 反復解法 連立一次方程式、 Lis 固有値ソルバ Xabclib 136 反復解法 (国産ライブラリ) 反復解法、自動チューニング (AT)機能 (国産ライブラリ) 2023年度 計算科学技術特論A

137.

その他のライブラリ(信号処理等) 種類 信号処理 問題 FFT ライブラリ名 FFTW 概要 FFTE 離散フーリエ変換 (国産ライブラリ) Spiral グラフ処理 137 グラフ分割 離散フーリエ変換、 AT機能 離散フーリエ変換、 AT機能 METIS、ParMETIS グラフ分割 SCOTCH、 グラフ分割 PT-SCOTCH 2023年度 計算科学技術特論A

138.

その他のライブラリ(フレームワーク) 種類 問題 ライブラリ名 概要 プログラミング 環境 マルチ フィジックス、 など Trilinos プログラミング フレームワークと 数値計算ライブラリ ステンシル 演算 Phisis ステンシル演算用 プログラミング フレームワーク (国産ライブラリ) FDM、FEM、DEM、 BEM、FVM ppOpen-HPC 5種の離散化手法に 基づくシミュレーション ソフトウェア、数値 ライブラリ、AT機能 (国産ライブラリ) 数値 ミドルウェア 138 2023年度 計算科学技術特論A

139.

レポート課題(その1)  問題レベルを以下に設定 問題のレベルに関する記述: •L00: きわめて簡単な問題。 •L10: ちょっと考えればわかる問題。 •L20: 標準的な問題。 •L30: 数時間程度必要とする問題。 •L40: 数週間程度必要とする問題。複雑な実装を必要とする。 •L50: 数か月程度必要とする問題。未解決問題を含む。 ※L40以上は、論文を出版するに値する問題。  教科書のサンプルプログラムは以下が利用可能     139 Sample-fx.tar Mat-Mat-noopt-fx.tar Mat-Vec-fx.tar Mat-Mat-fx.tar 2023年度 計算科学技術特論A

140.

レポート課題(その2) 1. 2. [L10] 利用できる計算機で、行列-行列積について、 メモリ連続アクセスとなる場合と、不連続となる場合 の性能を調査せよ。 [L15] 行列-行列積のアンローリングを、i, j, k ループについて施し、性能向上の度合いを調べよ。 どのアンローリング方式や段数が高速となるだろうか。 140 2023年度 計算科学技術特論A

141.

レポート課題(その3) 4. 5. 6. 7. [L15] 利用できる計算機で、ブロック化を行った 行列-行列積のコードに対し、アンローリングを 各ループについて施し性能を調査せよ。行列の大きさ (N)を変化させ、各Nに対して適切なアンローリング段 数を調査せよ。 [L5] 身近にある計算機の、キャッシュサイズと、 その構造を調べよ。 [L5] 身近にある計算機の、命令レベル並列性の 実装の仕組みを調べよ。 [L5] 本講義で取り扱っていないチューニング手法 を調べよ。 141 2023年度 計算科学技術特論A