内容に関する質問はkatagiri@cc.nagoya-u.ac.jp
まで

第１回
プログラム高速化の基礎
名古屋大学情報基盤センター 片桐孝洋
12017年度

計算科学技術特論A

本講義の位置づけ2

2017年度

計算科学技術特論A

講義日程と内容について
 ２０１７年度

計算科学技術特論Ａ（１学期：木曜3限 ）第１回：プログラム高速化の基礎、2017年4月13日
 イントロダクション、ループアンローリング、キャッシュブロック化、
数値計算ライブラリの利用、その他
 第２回：MPIの基礎、2017年4月20日
 並列処理の基礎、MPIインターフェース、MPI通信の種類、その他
 第３回：OpenMPの基礎、2017年4月27日
 OpenMPの基礎、利用方法、その他
 第４回：Hybrid並列化技法(MPIとOpenMPの応用)、2017年5月11日
 背景、Hybrid並列化の適用事例、利用上の注意、その他
 第５回：プログラム高速化の応用、2017年5月18日
 プログラムの性能ボトルネック に関する考えかた（I/O、単体性能(演算機ネック、メモリネック)、並列性能(バランス))、性能プロファイル、
その他


3

2017年度

計算科学技術特論A

参考書
 「計算科学のためのHPC技術1 」

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

2017年度

計算科学技術特論A

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

2017年度

計算科学技術特論A

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

6

2017年度

計算科学技術特論A

参考書
 「スパコンを知る:

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

2017年度

計算科学技術特論A

参考書


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

2017年度

計算科学技術特論A

イントロダクション
スパコンとは何か？

9

2017年度

計算科学技術特論A

スーパーコンピュータとは 「人工知能搭載のコンピュータではない」が・・・・AIスパコンが導入(2017-)されつつあり、近年のスパコンはAI搭載となりつつある

産総研「ABCI (AI Bridging Cloud Infrastructure)」

理研「ディープラーニング解析システム」

 明確な定義はない
 現在の最高レベルの演算性能をもつ計算機のこと
 経験的には、PCの1000倍以上 高速で、1000倍以上 大容量なメモリをもつ計算機

 1000倍高速だと 世界が違う！
 人の歩行速度：時速 約5km
 １０００倍だと 時速 5000km
 ジェット旅客機：速くて 時速1000km
 <ジェット旅客機の5倍の速度> と <歩く速さ> 最新鋭スパコンの能力はPCの10万倍以上高速
10

2017年度

計算科学技術特論A

スーパーコンピュータとは
 外為法安全保障貿易管理の外国為替及び外国貿易法の法令（平成２６年８月１４日公布、９月１５日施行）の規制対象
デジタル電子計算機
 第７条第三項ハ：

デジタル電子計算機であって、加重最高性能が八・〇実効テラ演算を超えるもの
 現在、ほとんどすべてのスーパーコンピュータは並列計算機
 名古屋大学情報基盤センターが所有する

FUJITSU PRIMEHPC FX100 も 並列計算機11

2017年度

計算科学技術特論A

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

計算科学技術特論A

スーパーコンピュータで用いる単位 PFLOPS（ぺタ・フロップス）
 １秒間に0.１京（けい）回の浮動小数点演算能力がある。
 「京コンピュータ」（2012年9月共用開始、11.2PFLOPS）
 PCの演算能力は？
 3.3GHｚ（１秒間に3.3G回のクロック周波数）として、もし１クロックあたり１回の浮動小数点演算ができれば3.3GFLOPS。
 Intel Core i7 (Sandy Bridge)では、6コア、１クロックで8回の浮動小数計算ができるので、3.3 GHz * 8回浮動小数点演算/Hz * 6コア = 158.4 GFLOPS
 Cray-１は160MFLOPS。 １９７０年代のスパコンより、PCの方が990倍以上高速！
13スーパーコンピュータ用語
 理論性能（Theoretical

Performance）

ハードウエア性能からはじき出した性能。
 １クロックに実行できる浮動小数点回数から算出した
FLOPS値を使うことが多い。
 実効性能（Effective Performance）

何らかのベンチマークソフトウエアを実行して実行時間を計測。
 そのベンチマークプログラムに使われている浮動小数点演算を算出。
 以上の値を基に算出したFLOPS値のこと。
 連立一次方程式の求解ベンチマークであるLINPACKを用いることが多い。


142017年度

計算科学技術特論A

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

15

2017年度

計算科学技術特論A

スーパーコンピュータ性能推移（主に日本製、理論性能）
Tianhe-2 (NUDT)
Titan (DOE/SC/ORNL)
Sequoia(DOE/NNSA/LLNL)
K-Computer (RIKEN)
Tianhe-1A(NUDT)FX100（名大）Jaguar(ORNL)
TUBAME（東工大）
SR11000(東大)T2K（東大）SX-8
FX1（JAXA）
SR8000(東大)SX-７ E2S（地球Sim）
SX-4
SX-６ 地球シミュレータ
SX-5
SX-4
SX-3
SR2201(東大)SX-2
VP-2600
S-810
S-820
VP-200(100PFLOPS)

Cray-1ILLIAC-IV

FACOM23016

2017年度

計算科学技術特論A

スーパーコンピュータのランキング
 TOP５００

Supercomputer Sites（http://www.top500.org/）
 LINPACKの値から実効性能を算出した値の５００位までのランキング
 米国オークリッジ国立研究所／テネシー大学
ノックスビル校の Jack Dongarra 教授が発案
 毎年、６月、１１月（米国の国際会議SC｜ｘｙ）に発表

17

2017年度

計算科学技術特論A

現在のランキング（2016年11月）

https://www.top500.org/list/2016/11/

18

2017年度

計算科学技術特論A

世界一のスパコン（中国） 中国の125PFLOPSスパコン2016.6- (ISC16)
 Sunway TaihuLight 神威太湖之光

National Supercomputing Center in Wuxi無錫（むしゃく）
国立スパコンセンター
 電力：15.3MW
 全体性能：125PFLOPS
 LINPACK: 93 PFLOPS（74%対ピーク）
 10,649,600コア（約1千万コア）
 Sunway SW26010 260C 1.45GHz


19

2017年度

計算科学技術特論A

日本最高速スパコン Oakforest-PACS(JCAHPC、東京大学情報基盤センター＆筑波大学計算科学研究センター)(Source: 全体構成 http://www.cc.u-tokyo.ac.jp/image/Oakforest-PACS.jpg)
内容仕様

Total Theoretical Peak

25.004 PFLOPS

Total Number of Nodes

8208

Total Memory Amount

897 TByte

Total Storage Amount

26 PB

Total Storage Bandwidth

500 GB/sec

Total Fast Cache Amount

940 TB

Total Fast Cache Bandwidth

1,560 GB/sec

8,208ノード (558,144 コア)
 ノード構成

CPUIntel® Xeon Phi ™ 7250 (Knights Landing, KNL)

#Processors

1 (68 Physical Cores )Frequency

1.4 GHz

Theoretical Peak for each Node3.0464 TFLOPS (double) / 6.0928 TFLOPS (single)

Memory96 GB(DDR4)＋ 16 GB(MCDRAM)

Inter20Connect

Intel ® Omni-Path Network (100 Gbps)

Intel Xeon Phi
KNLのアーキテクチャ (Knights Landing)メニーコア 各CPUあたり68物理コア
 HTで最大272スレッド

速い
遅い⇒高スレッド並列性の考慮

Non Uniform
Memory Access(NUMA)⇒NUMA affinityの考慮

512ビッドのSIMD演算⇒ループ中の演算量の考慮（アンローリング等）
21 Source: Knights Landing (KNL): 2nd Generation Intel® Xeon Phi™ Processor, Avinash

Sodani KNL Chief Architect Senior Principal Engineer, Intel Corp.

I/O技術の進展 2016年11月稼働の東大・筑波大(JCAHPC)の25PFLOPSスパコン
Oakforest-PACS
 ファイルシステム（26PB）の1PB分にDDN社のIME® (Infinite Memory Engine) (Source:
http://www.cc.u-tokyo.ac.jp/image/Oakforestを採用
PACS.jpg)
 I/OベンチマークIORにて
FPP(File Per Process : 並列プロセスがそれぞれ独立したファイルI/Oを行う)
SSF(Single Shared File : 全ての並列プロセスが単一共有ファイルにI/Oを行う)の異なるI/Oアクセスパターンにおいていずれも1TB/秒を達成(Source:
http://prtimes.jp/main/html/rd/p/000000098.000005769.html)
22

2017年度

計算科学技術特論A

AIスパコンの登場
 産総研「ABCI (AI Bridging Cloud Infrastructure)」 Source: http://www.hpcwire.jp/archives/11617
 2018年初頭までに130～200 PFLOPS (半精度、AI FLOPS) 電力：3MW未満
 環境温水冷却、大型リチウムイオン電池と高効率電源
 東京大学の柏キャンパスに設置予定
 理研「ディープラーニング解析システム」
 Source: http://pr.fujitsu.com/jp/news/2017/03/6.html
 2017年4月に稼働予定、富士通（株）が落札
 4PFLOPS

計算サーバ



Source:http://pr.fujitsu.com/jp/
news/2017/03/6.html

NVIDIA® Tesla® P100アクセラレーターを8基搭載のNVIDIA社「DGX-1」 24台「FujitsuPRIMERGY RX2530 M2」 32台

ストレージシステム



PCサーバ「FUJITSU Server PRIMERGY RX2540 M2」 6台
ストレージシステム「FUJITSU Storage ETERNUS（エターナス） DX200 S3」 8台「FUJITSU Storage ETERNUS DX100 S3」 1台
FUJITSU Software FEFS

23

2017年度

計算科学技術特論A

単体（CPU）最適化の方法

24

2017年度

計算科学技術特論A

高速

最近の計算機のメモリ階層構造
Ｏ（1ナノ秒）Ｏ（1０ナノ秒）

レジスタ

バイト

Ｋバイト
キャッシュ ～Ｍバイド

Ｏ（1００
ナノ秒）

メインメモリ

Ｍバイト
～Ｇバイド

Ｏ（1０
ミリ秒）

ハードディスク

Ｇバイト
～Ｔバイト
大容量

<メインメモリ>→<レジスタ>への転送コストは、レジスタ上のデータ・アクセスコストの Ｏ（100）倍！
25

2017年度

計算科学技術特論A

より直観的には…
レジスタ
キャッシュ

メインメモリ高性能（＝速い）プログラミングをするには、
きわめて小容量のデータ範囲について何度もアクセス（＝局所アクセス）するようにループを書くしかない26

2017年度

計算科学技術特論A

Fujitsu FX10のメモリ構成例
高速

レジスタ

●データ

レベル１キャッシュ
（32Ｋバイト/１コア）●データ

レベル２キャッシュ
（12Ｍバイト/16コア）

メインメモリ（３２Ｇバイト／ノード）
27

●データ

2017年度

大容量
計算科学技術特論A

Fujitsu FX10のメモリ構成例
レジスタ

高速

レベル１キャッシュ
（32Ｋバイト/１コア）

●データ

データが
Ｌ１キャッシュ上にあれば、速くアクセス可能●データ

レベル２キャッシュ
（12Ｍバイト/16コア）

大容量

メインメモリ（３２Ｇバイト／ノード）28

2017年度

計算科学技術特論A

Fujitsu FX10のノードのメモリ構成例
コア
１２

コア
１３

コア
１４

コア
１５

Ｌ１ … Ｌ１

Ｌ１

Ｌ１

Ｌ１

コア０ コア１ コア２ コア３

Ｌ１

Ｌ１

Ｌ１

Ｌ２
メインメモリ※階層メモリ構成となっている

29

2017年度

計算科学技術特論A

Fujitsu FX10全体メモリ構成
コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ２

Ｌ２

メインメモリ

メインメモリ…

Ｌ２

Ｌ２

メインメモリ

メインメモリ

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

コ
ア
０

コ
ア
１

コ
ア
２

コ
ア
３

コ
ア
１
２

コ
ア
１
３

コ
ア
１
４

コ
ア
１
５

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

… Ｌ
１

Ｌ
１

Ｌ
１

Ｌ
１

Ｌ２

Ｌ２

メインメモリ

メインメモリ…

Ｌ２

Ｌ２

メインメモリ

メインメモリＴＯＦＵネットワーク（５Ｇバイト／秒
×双方向）

メモリ階層が階層30

2017年度

計算科学技術特論A

FX10計算ノードの構成

TOFU
Network

１ソケットのみ

各ＣＰＵの内部構成

20GB/秒

Core
#0

Core
#1

Core
#2

Core
#3

…

Core
#12

Core
#13

Core
#14

Core
#15

L1

L1

L1

L1

: L1データキャッシュ32KB

L1

L1

L1

L1

ICC

L2 (16コアで共有、12MB)

85GB/秒=(8Byte×1333MHz
×8 channel)

MemoryMemoryMemoryMemory4GB ×2枚

4GB ×2枚

4GB ×2枚

4GB ×2枚

ノード内合計メモリ量：8GB×4＝32GB
31

2017年度

計算科学技術特論A

DDR3 DIMM

Fujitsu FX10のCPU(SPARC64IXfx)の詳細情報
項目値
アーキテクチャ名
HPC-ACE (SPARC-V9命令セット拡張仕様)
動作周波数
1.848GHz
L1キャッシュ
32 Kbytes (命令、データは分離)
L2キャッシュ
12 Mbytes
ソ フ ト ウ ェ ア 制 御 セクタキャッシュ
キャッシュ
演算実行
２整数演算ユニット、４つの浮動小数点積和演算ユニット（FMA）
SIMD命令実行
1命令で2つのFMAが動作
FMAは2つの浮動小数点演算（加算と乗算）を実行可能
レジスタ 浮動小数点レジスタ数：２５６本その他

 三角関数sin, cosの専用命令
 条件付き実行命令
 除算、平方根近似命令

32

2017年度

計算科学技術特論A

FX100計算ノードの構成
HMC
16GB

2ソケット、NUMA(Non Uniform Memory Access)

TOFU2Network

MemoryL2 (17コアで共有、12MB)
L1

L1

L1

L1

Core
#0

Core
#1

Core
#2

Core
#3

: L1データ
キャッシュ
64KB …

ICC

L1

L1

L1

L1

L1

Assist.
Core

Core
#12

Core
#13

Core
#14

Core
#15

Assist.
Core

… …

Core
#28

Core
#29

Core
#30

Core
#31

L1

L1

L1

L1

L1

… …

ソケット0 (CMG(Core Memory Group))
Core
#16

Core
#17

Core
#18

Core
#19

L1

L1

L1

L1

ソケット１ (CMC): L1データ
キャッシュ
…
64KB

L2 (17コアで共有、12MB)

MemoryHMC
16GB

ノード内合計メモリ量：32GB
33

読込み：240GB/秒
書込み：240GB/秒=合計：480GB/秒

FX10とFX100のアーキテクチャ比較
FX10

FX100

演算能力／ノード

倍精度／単精度：
236 GFLOPS

倍精度：1.011 TFLOPS
単精度：2.022 TFLOPS

演算コア数

16

32

アシスタントコア

なし

2

SIMD幅

128 ビット

256 ビット

SIMD命令浮動小数点演算、連続 右に加え、整数演算、ストライ
ロード／ストア
ド＆間接ロード／ストア

L1Dキャッシュ／コア

32KB、2ウェイ

64KB、4ウェイ

L2キャッシュ／ノード

12MB

24MB

メモリバンド幅

85GB/秒

480GB/秒 (HMC)

出典：https://www.ssken.gr.jp/MAINSITE/event/2015/20151028-sci/
lecture-04/SSKEN_sci2015_miyoshi_presentation.pdf
34

2017年度

計算科学技術特論A

FX100(名大)のCPU(SPARC64XIfx)の詳細情報
項目値
アーキテクチャ名
HPC-ACE2 (SPARC-V9命令セット拡張仕様)
動作周波数
2.2 GHz
L1キャッシュ64 Kbytes (命令、データは分離)
L2キャッシュ
24 Mbytes
ソ フ ト ウ ェ ア 制 御 セクタキャッシュ
キャッシュ
演算実行
２整数演算ユニット、８つの浮動小数点積和演算ユニット（FMA）
SIMD命令実行
1命令で２つのFMAが動作
FMAは４つの浮動小数点演算（加算と乗算）を実行可能
レジスタ 浮動小数点レジスタ数：２５６本その他

35

2017年度

計算科学技術特論A

演算パイプライン
演算の流れ作業36

2017年度

計算科学技術特論A

流れ作業
 車を作る場合
 １人の作業員１つの工程を担当（５名）

フロント・バッ
クガラスを
つける

車体作成

内装

外装

機能確認

 上記工程が２ヶ月だとする（各工程は0.4ヶ月とする）

２ヶ月後に１台できる
 ４ヶ月後に２台できる
 ２ヶ月／台 の効率


１台目
２台目
３台目
37

車体作成

フロント・
バックガ
ラスをつ
ける

内装

外装

• 各工程の作業員は、

０．４ヶ月働いて、
１．６ヶ月は休んでいる（＝作業効率が低い）

機能確
認
車体作成

フロント・
バックガ
ラスをつ
ける

内装

外装

機能確
認
車体作成

2017年度

フロント・
バックガ
ラスをつ
ける

内装

外装

計算科学技術特論A

機能確
認

時間

流れ作業
 作業場所は、５ヶ所とれるとする
 前の工程からくる車を待ち、担当工程が終わったら、次の工程に速やかに送られるとする

ベルトコンベア
０．４ヶ月
車体作成

38

０．４ヶ月
フロント・バック
ガラスをつける

０．４か月
内装

2017年度

０．４か月
外装

計算科学技術特論A

０．４か月
機能確認

流れ作業
 この方法では２ヶ月後に、１台できる
 ２．４ヶ月後に、２台できる
 ２．８ヶ月後に、３台できる
 ３．２ヶ月後に、４台できる
 ３．４ヶ月後に、５台できる
 ３．８ヶ月後に、６台できる
 ０．６３ヶ月／台 の効率


１台目
２台目
３台目
４台目
５台目
39

車体作成

フロント・
バックガ
ラスをつ
ける
車体作成

内装

外装

機能確
認

フロント・
バックガ
ラスをつ
ける

内装

外装

機能確
認

フロント・
バックガ
ラスをつ
ける

内装

外装

車体作成

車体作成

フロント・
バックガ
ラスをつ
ける
車体作成

機能確
認

内装

外装

機能確
認

フロント・
バックガ
ラスをつ
ける

内装

外装

機能確
認

2017年度

•各作業員は、十分に時間が立つと０．４か月の単位時間あたり休むことなく働いている（＝作業効率が高い）•このような処理を、＜パイプライン処理＞という時間
計算科学技術特論A

計算機におけるパイプライン処理の形態
ハードウエア・パイプライニング

1.

計算機ハードウエアで行う
以下の形態が代表的



1.
2.

ソフトウエア・パイプライニング

2.

プログラムの書き方で行う
以下の形態が代表的



1.
2.

40

演算処理におけるパイプライン処理
メモリからのデータ（命令コード、データ）転送におけるパイプライン処理

コンパイラが行うパイプライン処理（命令プリロード、データ・プリロード、データ・ポストストア）
人手によるコード改編によるパイプライン処理（データ・プリロード、ループアンローリング）2017年度

計算科学技術特論A

演算器の場合
 例：演算器の工程（注：実際の演算器の計算工程は異なる）
データＡをメモリから取る

データＢをメモリから取る

演算
を行う

演算結果を
収納

 行列-ベクトル積の計算では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]をメモリから取る41

2017年度

計算科学技術特論A

時間x[2]をメモリから取る

演算器の場合
 これでは演算器は、４単位時間のうち、１単位時間しか使われていないので無駄（＝演算効率１／４＝２５％）
 以下のようなパイプライン処理ができれば、十分時間が経つと、毎単位時間で演算がなされる（＝演算効率１００％）
十分な時間とは、十分な
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]を

x[2]をメモリから取る

A[0][2]*x[2]メモリから取るメモリから取るA[0][3]をメモリから取る

結果
y[0]収納ループ反復回数があること。
行列サイズNが大きいほど、
パイプラインが滞りなく流れ、演算効率は良くなる。
→Nが小さいと演算効率が悪い

x[3]をメモリから取る

A[0][3]*x[3]結果y[0]収納A[0][4]を

x[4]をメモリから取る

A[0][2]*x[4]メモリから取る…
42

2017年度

計算科学技術特論A

結果
y[0]収納時間

演算パイプラインのまとめ 演算器をフル稼働させるため（＝高性能計算するため）に必要な概念 メインメモリからデータを取ってくる時間はとても大きい。
演算パイプラインをうまく組めば、メモリからデータを取ってくる時間を＜隠ぺい＞できる（＝毎単位時間、演算器が稼働した状態にできる）
 実際は以下の要因があるので、そう簡単ではない
1.
2.
3.
4.
43

計算機アーキテクチャの構成による遅延（レジスタ数の制約、
メモリ→CPU・CPU→メモリへのデータ供給量制限、など）。
ループに必要な処理（ループ導入変数（i, j）の初期化と加算処理、
ループ終了判定処理）
配列データを参照するためのメモリアドレスの計算処理
コンパイラが正しくパイプライン化される命令を生成するか
2017年度

計算科学技術特論A

実際のプロセッサの場合
 実際のプロセッサでは
1.
2.

加減算
乗算

ごとに独立したパイプラインがある。
 さらに、同時にパイプラインに流せる命令（同時発行命令）が複数ある。
 Intel Pentium4ではパイプライン段数が３１段


演算器がフル稼働になるまでの時間が長い。

分岐命令、命令発行予測ミスなど、パイプラインを中断させる処理が多発すると、演算効率がきわめて悪くなる。
 近年の周波数の低い（低電力な）マルチコアCPU／メニーコアCPUでは、
パイプライン段数が少なくなりつつある（Xeon Phiは7段）


44

2017年度

計算科学技術特論A

FX10のハードウエア情報
 １クロックあたり、8回の演算ができる

浮動小数点積和演算ユニット（ＦＭＡ）あたり、乗算および加算が２つ（４つの浮動小数点演算）
 １クロックで、２つのＦＭＡが動作
 ４浮動小数点演算×2FMA＝８浮動小数点演算／クロック
 １コア当たり1.848ＧＨｚのクロックなので、
 理論最大演算は、

1.848 GHz* 8回 = 14.784 GFLOPS / コア １ノード16コアでは、
14.784 * 16コア = 236.5 GFLOPS / ノード
 レジスタ数（浮動小数点演算用）
 256個 / コア45

2017年度

計算科学技術特論A

ループ内連続アクセス46

2017年度

計算科学技術特論A

単体最適化のポイント
 配列のデータ格納方式を考慮して、連続アクセスすると速い（ループ内連続アクセス）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
47

for (i=0; i<n; i++) {for (j=0; j<n; j++) {a[ i ][ j ] = b[ j ] * c[ j ];}}2017年度for (jb=0; jb<n; jb+=m)for (i=0; i<n; i++) {for (j=jb; j<jb+m; j++) {OKa[ i ][ j ] = b[ j ] * c[ j ];}}}計算科学技術特論A

言語に依存した配列の格納方式の違い
 Ｆｏｒｔｒａｎ言語の場合

 Ｃ言語の場合

Ａ（i, j）

Ａ［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

48

2017年度格納方向

計算科学技術特論A

行列積コード例（Ｃ言語） コード例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

B
k

i
49

j

2017年度

計算科学技術特論A

行列の積
n

 行列積

cij   aik bkj (i, j  1, 2, ..., n)k 1の実装法は、次の二通りが知られている：

ループ交換法1.


ブロック化（タイリング）法

2.


50

連続アクセスの方向を変える目的で、行列-行列積を実現する３重ループの順番を交換する
キャッシュにあるデータを再利用する目的で、あるまとまった行列の部分データを、何度もアクセスするように実装する
2017年度

計算科学技術特論A

行列の積
 ループ交換法
 行列積のコードは、以下のような３重ループになる（Ｃ言語）

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 ];}}}
 最内部の演算は、外側の３ループを交換しても、計算結果が変わらない→ ６通りの実現の方法がある
51

2017年度

計算科学技術特論A

行列の積
 ループ交換法
 行列積のコードは、以下のような３重ループになる（Fortran言語）

do i=１，ndo j=１, n
do k=１, nc( i , j ) = c( i, j) + a( i , k ) * b( k , j )
enddo
enddo
enddo
 最内部の演算は、外側の３ループを交換しても、計算結果が変わらない→ ６通りの実現の方法がある
52

2017年度

計算科学技術特論A

行列の積
 行列データへのアクセスパターンから、以下の３種類に分類できる
1. 内積形式 (inner-product form)
最内ループのアクセスパタンが＜ベクトルの内積＞と同等
2. 外積形式 (outer-product form)
最内ループのアクセスパタンが＜ベクトルの外積＞と同等
3. 中間積形式 (middle-product form)
内積と外積の中間
53

2017年度

計算科学技術特論A

行列の積
 内積形式 (inner-product form）

ijk, jikループによる実現（Ｃ言語） 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;
}
}

A



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

2017年度B….

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

計算科学技術特論A

行列の積
 内積形式 (inner-product form）


ijk, jikループによる実現（Fortran言語）



do i=１, n
do j=１, n
dc = 0.0d0
do k=１, n
dc = dc + A( i , k ) * B( k , j )enddo
C( i , j ) = dc
enddo
enddo

2017年度B….

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

55

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

計算科学技術特論A

行列の積
 外積形式 (outer-product form）

kij, kjiループによる実現（Ｃ言語） for (i=0; i<n; i++) {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 ];for (i=0; i<n; i++) {C[ i ][ j ]= C[ i ][ j ]+ A[ i ][ k ]* db;
}
}
}


A

B
….

●kjiループでは
列方向アクセスがメイン→列方向格納言語向き（Ｆｏｒｔｒａｎ言語）

行列の積
 外積形式 (outer-product form）

kij, kjiループによる実現（Fortran言語）
 do i=１, n
do j=１, n
C( i , j ) = 0.0d0
enddo
enddo
do k=１, n
do j=１, n
db = B( k , j )do i=１, n
C( i , j ) = C( i , j )+ A( i , k ) * db
enddo
enddo
enddo


A

B
….

●kjiループでは
列方向アクセスがメイン→列方向格納言語向き（Ｆｏｒｔｒａｎ言語）

行列の積
 中間積形式 (middle-product form）

ikj, jkiループによる実現（Ｃ言語） 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;
}
}
}


58

2017年度AB.
.

●jkiループでは

全て列方向アクセス→列方向格納言語に最も向いている（Ｆｏｒｔｒａｎ言語）

計算科学技術特論A

行列の積
 中間積形式 (middle-product form）

ikj, jkiループによる実現（Fortran言語）
 do j=１, n
do i=１, n
C( i , j ) = 0.0d0
enddo
do k=１, n
db = B( k , j )do i=１, n
C( i , j ) = C( i , j ) + A( i , k ) * db
enddo
enddo
enddo


59

2017年度AB.
.

●jkiループでは

全て列方向アクセス→列方向格納言語に最も向いている（Ｆｏｒｔｒａｎ言語）

計算科学技術特論A

ループアンローリング60

2017年度

計算科学技術特論A

ループアンローリング コンパイラが、レジスタへのデータの割り当て；
2. パイプライニング；がよりできるようにするため、コードを書き換えるチューニング技法
 ループの刻み幅を、１ではなく、ｍにする
 ＜ｍ段アンローリング＞とよぶ
1.

61

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Ｃ言語）
 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+１]*B[k+１][ j]; k-ループのループ判定回数が１/2になる。

62

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Ｃ言語）
 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+１] += A[i][k] *B[k][ j+１];}
 A[i][k]をレジスタに置き、高速にアクセスできるようになる。

一般に：演算式が増えることで、ビット幅が大きなSIMD化ができる
63

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Ｃ言語）
 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+１][j] += A[i+１][k] *B[k][j];} B[i][j]をレジスタに置き、高速にアクセスできるようになる。

一般に：演算式が増えることで、ビット幅が大きなSIMD化ができる
64

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Ｃ言語）
 i-ループ、および j-ループ 2段展開

(nが２で割り切れる場合)

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+１] += A[i ][k] *B[k][ j+１];
C[i+１][ j ] += A[i+１][k] *B[k][ j ];C[i+１][ j+１] += A[i+１][k] *B[k][ j +１];}
 A[i][j], A[i+１][k],B[k][j],B[k][j+１]をレジスタに置き、高速にアクセスできるようになる。
65

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Ｃ言語）
 コンパイラにわからせるため、以下のように書く方がよい

場合がある for (i=0; i<n; i+=2)
for (j=0; j<n; j+=2) {dc00 = C[i ][ j ]; dc01 = C[i ][ j+１];
dc10 = C[i+１][ j ]; dc11 = C[i+１][ j+１] ;for (k=0; k<n; k++) {da0= A[i ][k] ; da1= A[i+１][k] ;db0= B[k][ j ]; db1= B[k][ j+１];
dc00 += da0 *db0; dc01 += da0 *db1;
dc10 += da1 *db0; dc11 += da1 *db1;}
C[i ][ j ] = dc00; C[i ][ j+１] = dc01;
C[i+１][ j ] = dc10; C[i+１][ j+１] = dc11;}66

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Fortran言語）
 k-ループ2段展開(nが2で割り切れる場合)

do i=１, n
do j=１, n
do k=１, n, 2C(i, j) = C(i, j) +A(i, k) *B(k, j) + A(i, k+１)*B(k+１, j)
enddo
enddo
enddo
 k-ループのループ判定回数が１/2になる。
67

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Fortran言語）
 j-ループ2段展開(nが2で割り切れる場合)

do i=１, n
do j=１, n, 2
do k=１, nC(i, j ) = C(i, j ) +A(i, k) * B(k, j )C(i, j+１) = C(i, j+１) +A(i, k) * B(k, j+１)
enddo
enddo
enddo
 A(i, k)をレジスタに置き、高速にアクセスできるようになる。
68

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Fortran言語）
 i-ループ2段展開(nが2で割り切れる場合)

do i=１, n, 2
do j=１, n
do k=１, n
C(i , j) = C(i , j) +A(i , k) * B(k , j)C(i+１, j) = C(i+１, j) +A(i+１, k) * B(k , j)
enddo
enddo
enddo
 B(i, j)をレジスタに置き、高速にアクセスできるようになる。
69

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Fortran言語）
 i-ループ、および j-ループ 2段展開

(nが２で割り切れる場合)

do i=１, n, 2
do j=１, n, 2
do k=１, n
C(i , j ) = C(i , j ) +A(i , k) *B(k, j )C(i , j+１) = C(i , j+１) +A(i , k) *B(k, j+１)
C(i+１, j ) = C(i+１, j ) +A(i+１, k) *B(k, j )C(i+１, j+１) =C(i+１, j+１) +A(i+１, k) *B(k, j +１)
enddo; enddo; enddo;
 A(i,j), A(i+１,k),B(k,j),B(k,j+１)をレジスタに置き、高速にアクセスできるようになる。
70

2017年度

計算科学技術特論A

ループアンローリングの例（行列-行列積、Fortran言語） コンパイラにわからせるため、以下のように書く方がよい

場合がある do i=１, n, 2
do j=１, n, 2

dc00 = C(i ,j ); dc01 = C(i ,j+１)
dc10 = C(i+１,j ); dc11 = C(i+１,j+１)
do k=１, nda0= A(i ,k); da1= A(i+１, k)
db0= B(k ,j ); db1= B(k, j+１)
dc00 = dc00+da0 *db0; dc01 = dc01+da0 *db1;
dc10 = dc10+da1 *db0; dc11 = dc11+da1 *db1;
enddo
C(i , j ) = dc00; C(i , j+１) = dc01
C(i+１, j ) = dc10; C(i+１, j+１) = dc11
enddo; enddo
71

2017年度

計算科学技術特論A

キャッシュライン衝突
とびとびアクセスは弱い

72

2017年度

計算科学技術特論A

不連続アクセスとは Ｃ言語の場合

 配列のデータ格納方式を考慮し連続アクセスすると速い（ループ内連続アクセス）NG

a［i］［j］
jfor (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

間隔４での不連続アクセス73

2017年度

計算科学技術特論A

キャッシュメモリの構成
キャッシュライン

ＣＰＵ

キャッシュメモリ（キャッシュ上のブロック）

演算器
演算
要求

10

演算
結果

0

6

2

14

レジスタ
データ供給

データ供給

メインメモリ

セット

データ蓄積

（ブロック
の並び）

ブロック

メインメモリ（記憶単位）

注）配列をアクセスすると、１要素分ではなくブロック単位のデータ（例えば32バイト（倍精度4変数分）
が同時にキャッシュに乗る（ブロックサイズと呼ぶ）
74

写像関数

データ蓄積

キャッシュメモリ

ブロックとキャッシュラインの対応

2017年度

8

9

10 11 12 13 14

0

1

2

計算科学技術特論A

3

4

6

7

キャッシュとキャッシュライン メインメモリ上とキャッシュ上のデータマッピング方式


読み出し： メインメモリ から キャッシュ へ ダイレクト・マッピング方式：

メモリバンクごとに直接的
 セット・アソシアティブ方式： ハッシュ関数で写像（間接的）

書き込み： キャッシュ から メインメモリ へ ストア・スルー方式：

キャッシュ書き込み時にメインメモリと中身を一致させる
 ストア・イン方式： 対象となるキャッシュラインが置き換え対象となったときに一致させる
キャッシュメモリ
キャッシュ
ライン
ライン０
ライン１
ライン２
ライン３
ライン４
ライン５

写像関数
メモリバンク…

75

メインメモリ…
2017年度

計算科学技術特論A

キャッシュライン衝突の例
直接メインメモリのアドレスをキャッシュに写像する、ダイレクト・マッピングを考える



物理結線は以下の通り

マッピング間隔を、ここでは４とする



メインメモリ上のデータは、間隔４ごとに、同じキャッシュラインに乗る

キャッシュラインは8バイト、メモリバンクも8バイトとする
配列aは 4×4の構成で、倍精度（8バイト）でメモリ確保されているとする

double a[4][4];

この前提で、格納方向と逆方向にアクセス（４とびのアクセス）する
メモリ（＝Ｃ言語の場合、i方向を連続アクセス）
メインメモリ連続方向
キャッシュ
キャッシュメモリ
１
２
３
４
ライン
ライン０
５
６
７
８
ライン１
９ １０ １１ １２
ライン２
１３ １４ １５ １６
ライン３
物理結線
配列アクセス方向

76

2017年度

計算科学技術特論A…

キャッシュライン衝突の例
 この前提の、＜実際の配列構成＞と＜メモリバンク＞の関係
実際は、以下のことがあるので、必ずしも、こうならないことに注意する

配列a[][]の物理メモリ上の配置はOSが動的に決定するので、ずれることがある

メモリバンクの容量は、8バイトより大きいダイレクト・マッピングではない Ｃ言語の場合

配列a［i］［j］

メインメモリ上のバンク構成

j
配列要素a[][] とメモリバンク構造と
が完全一致

１

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
77

2017年度

計算科学技術特論A…

キャッシュライン衝突の例
1.

a[0][0]があるバンク1がキャッシュライン0に乗る

2.

すぐに、a[1][0]があるバンク5がアクセスされる

3.

（物理結線先のキャッシュライン0に容量の空きがないので）キャッシュライン0のデータ(バンク1の内容)を追い出さないといけない

4.

バンク5のデータがキャッシュライン０に乗る

5.

すぐに、a[2][0]があるバンク９がアクセスされる

6.

キャッシュライン０のデータ（バンク5の内容）を追い出さないといけない…玉突きで、ライン1～3が空いていても、逐次的にキャッシュ上のデータが追い出される
キャッシュ
ライン

レジスタへ

78

メモリ連続

メインメモリキャッシュメモリ

１
５
９
１３

９
５
１
ライン０
ライン１
ライン２
ライン３

２
６
１０
１４

配列アクセス方向
2017年度

計算科学技術特論A

…

３
７
１１
１５

４
８
１２
１６

キャッシュライン衝突の例
 １～６の状態が連続して発生する。

メモリ→キャッシュの回線が常に稼働 ＜回線お話し中＞で、データが来るのが終わるまで、待たされる（回線レベルで並列にデータが持ってこれない） ストア・イン方式では、メモリにデータを書き戻すコストもかかる メモリからデータを逐次で読み出すのと同じ＜キャッシュがない＞のと同じ演算器にデータが届かないので計算を中断。
演算器の利用効率が悪くなる

以上の現象を＜キャッシュライン衝突＞と呼ぶ
79

2017年度

計算科学技術特論A

メモリ・インターリービング 物理的なメモリの格納方向に従いアクセスする時

データアクセス時、現在アクセス中のバンク上のデータは、
周辺バンク上のデータも一括して（同時に）、別のキャッシュライン上に乗せるハードウェア機能がある
キャッシュライン０のデータをアクセスしている最中に、
キャッシュライン１に近隣のバンク内データを（並列に）持ってくることが可能

メモリの＜インタリービング＞
演算機から見たデータアクセス時間が短縮
演算器が待つ時間が減少（＝演算効率が上がる）

物理的なデータ格納方向に連続アクセスするとよい80

2017年度

計算科学技術特論A

キャッシュライン衝突が起こる条件
 メモリバンクのキャッシュラインへの割り付けは２冪の間隔で行っていることが多い たとえば、３２、６４、１２８など
 特定サイズの問題（たとえば１０２４次元）で、性能が１／２～１／３、ときには１／１０になる
場合、キャッシュライン衝突が生じている可能性あり
double a[1024][1024];
NG

double precision a(1024, 1024)

実際は、OSやキャッシュ構成の影響で厳密な条件を見つけることは難しいが2冪サイズでの配列確保は避けるべき81

2017年度

計算科学技術特論A

キャッシュライン衝突への対応
 キャッシュライン衝突を防ぐ方法

パティング法： 配列に（２冪でない）余分な領域を確保
し確保配列の一部の領域を使う。

1.

余分な領域を確保して使う
 例： double A[1024][1025]; で1024のサイズをアクセス
 コンパイラのオプションを使う


2.

3.

82

データ圧縮法： 計算に必要なデータのみキャッシュライン衝突しないようにデータを確保し、かつ、必要なデータをコピーする。
予測計算法： キャッシュライン衝突が起こる回数を予測するルーチンを埋め込み、そのルーチンを配列確保時に呼ぶ。
2017年度

計算科学技術特論A

ブロック化
小さい範囲のデータ再利用

83

2017年度

計算科学技術特論A

ブロック化によるアクセス局所化
 キャッシュには大きさがあります。
 この大きさを超えると、たとえ連続アクセスしても、

キャッシュからデータは追い出されます。
 データが連続してキャッシュから追い出されると、メモリから転送するのと同じとなり、高速な
アクセス速度を誇るキャッシュの恩恵がなくなります。
 そこで、高速化のためには、以下が必要です
1. キャッシュサイズ限界までデータを詰め込む
2. 詰め込んだキャッシュ上のデータを、何度もアクセスして再利用する
84

2017年度

計算科学技術特論A

ブロック化によるキャッシュミスヒット
削減例 行列ー行列積
 行列サイズ：８×８
 double A[8][8]; キャッシュラインは４つ
 １つのキャッシュラインに４つの行列要素が載る
 キャッシュライン：4×8バイト(double)=32バイト

 配列の連続アクセスは行方向（C言語） キャッシュの追い出しアルゴリズム：

Least Recently Used (LRU)85

2017年度

計算科学技術特論A

配列とキャッシュライン構成の関係
 この前提の、＜配列構成＞と＜キャッシュライン＞の関係


ここでは、キャッシュライン衝突は考えません Ｃ言語の場合

配列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

 １×４の配列要素が、キャッシュラインに乗る
 どのキャッシュラインに乗るかは、<配列アクセス
パターン> と <置き換え
アルゴリズム>依存で決まる

キャッシュラインの
構成
１
3

格納方向

i
86

2017年度

計算科学技術特論A

24

行列-行列積の場合（ブロック化しない）
LRU:直近で最もアクセスされていないラインのデータを追い出すキャッシュミス①

キャッシュミス②
キャッシュミス③
キャッシュミス④＝

Ｃ＊キャッシュミス⑤

Ａ

Ｂ

キャッシュライン
ライン１
ライン２
ライン３
ライン４
87

2017年度※キャッシュライン４つ、置き換えアルゴリズム
LRUの場合計算科学技術特論A

行列-行列積の場合（ブロック化しない）

キャッシュミス１１

キャッシュミス⑥＝＊

キャッシュミス⑦
キャッシュミス⑧
キャッシュミス⑨

Ｃキャッシュミス⑩

Ａ

Ｂキャッシュライン
ライン１
ライン２
ライン３
ライン４
88

2017年度※キャッシュライン４つ、置き換えアルゴリズム
LRUの場合計算科学技術特論A

行列-行列積の場合（ブロック化しない）
ライン１
ライン２
ライン３
ライン４
キャッシュミスキャッシュミスキャッシュミスキャッシュミスキャッシュミスキャッシュミス＝キャッシュミスキャッシュミス＊キャッシュミスキャッシュミスキャッシュミスＣ

Ａ

Ｂ

キャッシュライン※２要素計算するのに、キャッシュミスヒット２２回
89

2017年度※キャッシュライン４つ、置き換えアルゴリズム
LRUの場合計算科学技術特論A

行列-行列積の場合（ブロック化する：2要素）
ライン１
ライン２
ライン３
ライン４

このブロック幅
単位で計算するキャッシュミス

１２

キャッシュミスキャッシュミス

１
①

②

＝

Ｃ

①

キャッシュミスキャッシュミス②＊

Ａ

キャッシュミスＢ

キャッシュライン※キャッシュライン４つ、置き換えアルゴリズム
LRUの場合90

2017年度

計算科学技術特論A

行列-行列積の場合（ブロック化する：2要素）
ライン１
ライン２
ライン３
ライン４

このブロック幅
単位で計算する

１２

キャッシュミスキャッシュミス１③＝④＊
キャッシュミス③

キャッシュミスキャッシュミス④

Ｃキャッシュライン

91

キャッシュミスＡ

Ｂ※２要素計算するのに、キャッシュミスヒット１０回
2017年度※キャッシュライン４つ、置き換えアルゴリズム
LRUの場合計算科学技術特論A

行列積コード（Ｃ言語）：キャッシュブロック化なし コード例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

B
k

i
92

j

2017年度

計算科学技術特論A

行列-行列積のブロック化のコード（Ｃ言語）
 nがブロック幅（ibl=１6）で割り切れるとき、以下のような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];} } } } } }93

2017年度

計算科学技術特論A

行列-行列積のブロック化のコード（Fortran言語）
 nがブロック幅（ibl=１6）で割り切れるとき、以下のような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;
94

2017年度

計算科学技術特論A

キャッシュブロック化時のデータ・アクセスパターン
ibl×iblの小行列単位で行列‐行列積
をする

ibl

ibl
ibl

ibl

ibl

＝

×

C
95

ibl

A
2017年度

計算科学技術特論AB

キャッシュブロック化時のデータ・アクセスパターン
ibl×iblの小行列単位で行列‐行列積
をする

ibl

ibl
ibl

ibl

ibl
×

＝

ibl

C
96

A
2017年度

計算科学技術特論AB

行列-行列積のブロック化のコードのアンローリング（Ｃ言語） 行列-行列積の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];} } } } } }97

2017年度

計算科学技術特論A

行列-行列積のブロック化のコードのアンローリング（Fortran言語） 行列-行列積の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;98

2017年度

計算科学技術特論A

その他の高速化技術

99

2017年度

計算科学技術特論A

共通部分式の削除（１） 以下のプログラムは、冗長な部分がある。

d = a + b + c;
f = d + a + b;
 コンパイラがやる場合もあるが、以下のように書く方が無難である。

temp = a + b;
d = temp + c;
f = d + temp;100

2017年度

計算科学技術特論A

共通部分式の削除（２） 配列のアクセスも、冗長な書き方をしないほうがよい。

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];}1012017年度

計算科学技術特論A

コードの移動 割り算は演算時間がかかる。ループ中に書かない。

for (i=0; i<n; i++) {a[i] = a[i] / sqrt(dnorm);} 上記の例では、掛け算化して書く。

dtemp = １.0d0 / sqrt(dnorm);for (i=0; i<n; i++) {a[i] = a[i] *dtemp;}
102

2017年度

計算科学技術特論A

ループ中のＩＦ文 なるべく、ループ中にＩＦ文を書かない。

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] = １.0d0;} } 以下のように書く。

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] = １.0d0;103

2017年度

計算科学技術特論A

ソフトウェア・パイプライニングの強化 基のコード（２段のアンローリング）
定義－参照の距離が近い→ソフトウェア的には何もできない ソフトウェアパイプライニングを強化したコード（２段のアンローリング）
定義－参照の距離が遠い→ソフトウェアパイプライニングが適用できる機会が増加！
104

2017年度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;}
計算科学技術特論A

数値計算ライブラリの利用

105

2017年度

計算科学技術特論A

数値計算ライブラリ 密行列用ライブラリ

行列の要素に0がない（というデータ構造を扱う）
 連立一次方程式の解法、固有値問題、FFT、その他
 直接解法（反復解法もある）
 BLAS、LAPACK、ScaLAPACK、SuperLU、MUMPS、FFTW、など
 疎行列用ライブラリ

行列の要素に0が多い
 連立一次方程式の解法、固有値問題、その他
 反復解法 PETSc、Xabclib、Lis、ARPACK、など
106

2017年度

計算科学技術特論A

疎行列用ライブラリの特徴
 疎行列を扱うアプリケーションはライブラリ化が難しい

疎行列データ形式の標準化が困難


COO、CRS(CCS)、ELL、JDS、BCSR、・・・

カーネルの演算が微妙に違う、かつ、カーネルは広い範囲に分散

陽解法（差分法）を基にしたソフトウェア

 数値ミドルウェアおよび領域特化型言語（Domain Specific Language, DSL）
解くべき方程式や離散化方法に特化させることで、処理（対象となるプロ
グラムの性質）を限定
 以上の限定から、高度な最適化ができる言語（処理系）の作成（DSL）や、ライブラリ化（数値ミドルウェア）ができる
 数値ミドルウェアの例




107

ppOpen-HPC(東大)、PETSc(Argonne National Laboratory, USA.) 、
Trilinos (Sandia National Laboratory, USA)、など
2017年度

計算科学技術特論A

BLAS
 BLAS（Basic Linear Algebra Subprograms、

基本線形代数副プログラム集）
 線形代数計算で用いられる、基本演算を標準化（API化）したもの。
 普通は、密行列用の線形代数計算用の基本演算の副プログラムを指す。
 疎行列の基本演算用の＜スパースBLAS＞というものあるが、まだ定着していない。
 スパースBLASはIntel MKL(Math Kernel Library)に入っているが、広く使われているとは言えない。
108

2017年度

計算科学技術特論A

BLAS
 BLASでは、以下のように分類わけをして、サブルーチンの命名規則を統一
1.
2.
3.
4.

演算対象のベクトルや行列の型（整数型、実数型、複素型）
行列形状（対称行列、三重対角行列）
データ格納形式（帯行列を二次元に圧縮）演算結果が何か（行列、ベクトル） 演算性能から、以下の３つに演算を分類
 レベル１ BLAS：

ベクトルとベクトルの演算 レベル２ BLAS： 行列とベクトルの演算 レベル３ BLAS： 行列と行列の演算
109

2017年度

計算科学技術特論A

レベル１ BLAS
 レベル１BLAS ベクトル内積、ベクトル定数倍の加算、など

y ← αx+ y データの読み出し回数、演算回数がほほ同じ
 データの再利用（キャッシュに乗ったデータの再利用によるデータアクセス時間の短縮）がほとんどできない 実装による性能向上が、あまり期待できない ほとんど、計算機ハードウエアの演算性能
 レベル１BLASのみで演算を実装すると、演算が本来持っているデータ再利用性がなくなる
 例：行列-ベクトル積を、レベル１BLASで実装
 例：

110

2017年度

計算科学技術特論A

レベル２ BLAS
BLAS
 行列-ベクトル積などの演算

 レベル２

 例： y ← α A x + β y 前進/後退代入演算、T x = y （Tは三角行列）をxについて解く演算、を含む
 レベル１BLASのみの実装よる、データ再利用性の喪失を回避する目的で提案
 行列とベクトルデータに対して、データの再利用性あり
 データアクセス時間を、実装法により短縮可能 （実装法により）性能向上がレベル１BLASに比べしやすい（が十分でない）
111

2017年度

計算科学技術特論A

レベル３ BLAS
BLAS
 行列-行列積などの演算

 レベル３C←αAB+βC
 共有記憶型の並列ベクトル計算機では、レベル２ BLASでも性能向上が達成できない。
 例：

 並列化により１PE当たりのデータ量が減少する。
 より大規模な演算をとり扱わないと、再利用の効果がない。

 行列-行列積では、行列データ

2

O(n ) に対して演算は O (n 3 ) なので、データ再利用性が原理的に高い。
 行列積は、アルゴリズムレベルでもブロック化できる。
さらにデータの局所性を高めることができる。
112

2017年度

計算科学技術特論A

典型的なBLASの性能
性能[FLOPS]

理論性能の限界

BLAS３

BLAS2
BLAS1
行列サイズ
113

2017年度

計算科学技術特論A

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 )
alpha
の値

Aが転置しているか
Bが転置しているか

Aの
アドレス

Mの大きさ

Aの１次元目
の要素数

Nの大きさ

Bの
アドレス

Bの１次元目
の要素数

Kの大きさ

114

2017年度

betaの値

計算科学技術特論A

Cのアドレス
Cの１次元目の要素数

BLASの機能詳細
http://www.netlib.org/blas/ 命名規則： 関数名：XYYYY 詳細はHP:
 X：

データ型
S:単精度、D：倍精度、C：複素、Z：倍精度複素 YYYY： 計算の種類
 レベル１：

例：AXPY：ベクトルをスカラー倍して加算 レベル２：
例：GEMV: 一般行列とベクトルの積 レベル３：
例：GEMM:一般行列どうしの積115

2017年度

計算科学技術特論A

GOTO BLASとは 後藤和茂 氏により開発された、ソースコードが無償入手可能な、高性能ＢＬＡＳの実装（ライブラリ） 特徴
マルチコア対応がなされている
 多くのコモディティハードウエア上の実装に特化 Intel Nehalem and Atom systems
 VIA Nanoprocessor
 AMD Shanghai and Istanbul等 テキサス大学先進計算センター（ＴＡＣＣ）で、GOTO BLAS2として、ソースコードを配布している
116

https://www.tacc.utexas.edu/research-development/tacc-software/gotoblas2
2017年度

計算科学技術特論A

LAPACK 密行列に対する、連立一次方程式の解法、および固有値の解法の“標準”アルゴリズムルーチンを無償で提供
 その道の大学の専門家が集結
カリフォルニア大バークレー校：
James Demmel教授
 テネシー大ノックスビル校：
Jack Dongarra教授
 HP

http://www.netlib.org/lapack/117

2017年度

計算科学技術特論A

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

118

2017年度

計算科学技術特論Aインタフェース例：DGESV （１/３）
 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
 線形方程式の数。行列Aの次元数。

119

2017年度N >= 0。

計算科学技術特論Aインタフェース例：DGESV （２/３）
NRHS (入力) – INTEGER 右辺ベクトルの数。行列Bの次元数。 NRHS >= 0。

A (入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
 入力時は、N×Nの行列Aの係数を入れる。

= P*L*Uを圧縮して出力する。
Lの対角要素は１であるので、収納されていない。

 出力時は、Aから分解された行列LとU

LDA (入力) – INTEGER 配列Aの最初の次元の大きさ。 LDA>= max(1,N)。

IPIVOT (出力) － DOUBLE PRECISION, DIMENSION(:)
 交換行列Aを構成する枢軸のインデックス。 行列のi行がIPIVOT(i)行と交換されている。

120

2017年度

計算科学技術特論Aインタフェース例：DGESV （３/３）
B (入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
N×NRHS 行列Bを入れる。
 出力時は、もし、INFO = 0 なら、N×NRHS行列である解行列Xが戻る。
 入力時は、右辺ベクトルのLDB (入力) －INTEGER 配列Bの最初の次元の大きさ。 LDB >= max(1,N)。



INFO (出力) ーINTEGER
 = 0: 正常終了
 < 0: もし INFO = -i

なら i-th 行の引数の値がおかしい。
 > 0: もし INFO = i なら U(i,i) が厳密に0である。 分解は終わるが、Uの分解は特異なため、解は計算されない。

121

2017年度

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

2017年度

計算科学技術特論AScaLAPACKのソフトウェア構成図
出典：http://www.netlib.org/scalapack/poster.html

分散メモリ用アルゴリズムのライブラリ

ScaLAPACK

分散メモリ用演算カーネル
ライブラリ

PBLAS
キャッシュ
最適化アルゴリズムのライブラリ

大域アドレス
局所アドレス
LAPACK
環境独立
環境依存
BLAS

ScaLAPACK用
通信ライブラリ

演算カーネル
ライブラリ

123

BLACS

Message Passing
Interface (MPI)
2017年度

計算科学技術特論A

汎用
通信ライブラリ

BLACSとPBLAS BLACS
ScaLAPACK中で使われる通信機能を関数化したもの。
 通信ライブラリは、MPI、PVM、各社が提供する通信ライブラリを想定し、ScaLAPACK内でコード修正せずに使うことを目的とする いわゆる、通信ライブラリのラッパー的役割でScaLAPACK内で利用

現在、MPIがデファクトになったため、MPIで構築された
BLACSのみ、現実的に利用されている。
 なので、ScaLAPACKはMPIでコンパイルし、起動して利用する

 PBLASBLACSを用いてBLASと同等な機能を提供する関数群
 並列版BLASといってよい。
124

2017年度

計算科学技術特論AScaLAPACKの命名規則
 原則：

LAPACKの関数名の頭に“P”を付けたもの
 そのほか、BLACS、PBLAS、データ分散を制御するためのScaLAPACK用関数がある。

125

2017年度

計算科学技術特論Aインタフェース例：PDGESV （１/４）
 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-１, JA:JA+N-１)の行列
 X と B は N×NRHS行列を分散したB(IB:IB+N-１, JB:JB+NRHS-１)の行列
 行交換の部分枢軸選択付きのLU分解 でsub(A) をsub(A) = P * L * U と分解する。ここで、P は交換行列、
L は下三角行列、Uは上三角行列である。
 分解されたsub(A) は、連立一次方程式sub(A) * X = sub(B)を解くのに使われる。
126

2017年度

計算科学技術特論Aインタフェース例：PDGESV （２/４）
N (大域入力) – INTEGER
 線形方程式の数。行列Aの次元数。

N >= 0。

NRHS (大域入力) – INTEGER 右辺ベクトルの数。行列Bの次元数。 NRHS >= 0。

A (局所入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
 入力時は、N×Nの行列Aの局所化された係数を配列A(LLD_A, LOCc( JA+N-１))を入れる。
 出力時は、Aから分解された行列LとU = P*L*Uを圧縮して出力する。
Lの対角要素は１であるので、収納されていない。

IA(大域入力) －INTEGER ：sub(A)の最初の行のインデックス
 JA(大域入力) －INTEGER ：sub(A)の最初の列のインデックス
 DESCA (大域かつ局所入力) – INTEGER 分散された配列Aの記述子。
127

2017年度

計算科学技術特論Aインタフェース例：PDGESV （３/４）


IPIVOT (局所出力) － DOUBLE PRECISION, DIMENSION(:)
 交換行列Aを構成する枢軸のインデックス。 行列のi行がIPIVOT(i)行と交換されている。分散された配列( LOCr(M_A)+MB_A )として戻る。
B (局所入力／出力) – DOUBLE PRECISION, DIMENSION(:,:)
N×NRHSの行列Bの分散されたものを(LLD_B, LOCc(JB+NRHS-１)) に入れる。
 出力時は、もし、INFO = 0 なら、N×NRHS行列である解行列Xが、行列Bと同様の分散された状態で戻る。
 入力時は、右辺ベクトルのIB(大域入力) －INTEGER sub(B)の最初の行のインデックスJB(大域入力) －INTEGER sub(B)の最初の列のインデックス



DESCB (大域かつ局所入力) – INTEGER
 分散された配列Bの記述子。

128

2017年度

計算科学技術特論Aインタフェース例：PDGESV （４/４）


INFO (大域出力) ーINTEGER
 = 0: 正常終了
 < 0:
 もし i番目の要素が配列で、 そのj要素の値がおかしいなら、

INFO = -(i*100+j)となる。
 もしi番目の要素がスカラーで、かつ、その値がおかしいなら、
INFO = -iとなる。
 > 0: もし INFO = Kのとき U(IA+K-１, JA+K-１) が厳密に0である。
分解は完了するが、分解されたＵは厳密に特異なので、解は計算できない。

129

2017年度

計算科学技術特論A

BLAS利用の注意
 C言語からの利用

BLASライブラリは（たいてい）Fortranで書かれている 行列を１次元で確保する

Fortranに対して転置行列になるので、BLASの引数で転置を指定

引数は全てポインタで引き渡す
 関数名の後に“_”をつける（BLASをコンパイルするコンパイラ依存） 例：dgemm_(...)
 小さい行列は性能的に注意キャッシュに載るようなサイズ（例えば100次元以下）の行列については、BLASが高速であるとは限らない


全体の行列サイズは大きくても、利用スレッド数が多くなると、
スレッド当たりの行列サイズが小さくなるので注意！


130

BLASは、大規模行列で高性能になるように設計されている

例） N=8000でも272スレッド並列だと、スレッドあたり約480x480 まで小さくなる
2017年度

計算科学技術特論A

その他のライブラリ（主に行列演算）
種類

問題

ライブラリ名

概要

密行列

BLAS

MAGMA

GPU、マルチコア、ヘテロジニ
アス環境対応疎行列

連立一次方程式

MUMPS

直接解法

SuperLU

直接解法

PETSc

反復解法、各種機能

Hypre

反復解法

連立一次方程式、 Lis
固有値ソルバ
Xabclib131

2017年度

反復解法（国産ライブラリ）
反復解法、自動チューニング
（AT）機能（国産ライブラリ）

計算科学技術特論A

その他のライブラリ（信号処理等）
種類
信号処理

問題
FFT

ライブラリ名
FFTW
FFTE

132

グラフ分割

離散フーリエ変換、
AT機能
離散フーリエ変換（国産ライブラリ）

Spiral
グラフ処理

概要

離散フーリエ変換、
AT機能
METIS、ParMETIS グラフ分割
SCOTCH、
グラフ分割
PT-SCOTCH

2017年度

計算科学技術特論A

その他のライブラリ（フレームワーク）
種類

問題

ライブラリ名概要

プログラミング
環境

マルチ
フィジックス、など

Trilinos

プログラミング
フレームワークと数値計算ライブラリ

ステンシル
演算

Phisisステンシル演算用
プログラミング
フレームワーク（国産ライブラリ）

FDM、FEM、DEM、BEM、FVM

ppOpen-HPC

5種の離散化手法に基づくシミュレーション
ソフトウェア、数値ライブラリ、AT機能（国産ライブラリ）

数値
ミドルウェア133

2017年度

計算科学技術特論A

レポート課題（その１） 問題レベルを以下に設定
問題のレベルに関する記述：•L00: きわめて簡単な問題。
•L10： ちょっと考えればわかる問題。
•L20： 標準的な問題。
•L30： 数時間程度必要とする問題。
•L40： 数週間程度必要とする問題。複雑な実装を必要とする。
•L50： 数か月程度必要とする問題。未解決問題を含む。
※L４０以上は、論文を出版するに値する問題。

 教科書のサンプルプログラムは以下が利用可能

Sample-fx.tar
 Mat-Mat-noopt-fx.tar
 Mat-Vec-fx.tar
 Mat-Mat-fx.tar
134

2017年度

計算科学技術特論A

レポート課題（その２）
1.

2.

3.

[L１０] 利用できる計算機で、行列-行列積について、
メモリ連続アクセスとなる場合と、不連続となる場合の性能を調査せよ。
[L１５] 行列-行列積のアンローリングを、i, j, k
ループについて施し、性能向上の度合いを調べよ。
どのアンローリング方式や段数が高速となるだろうか。
[L１０] FX10のＣＰＵであるSPARC64 IXfx、もしくはFX100のCPUであるSPARC64 XIfxの計算機アーキテ
クチャについて調べよ。特に、演算パイプラインの構成や、演算パイプラインに関連するマシン語命令について調べよ。
135

2017年度

計算科学技術特論A

レポート課題（その３）
4.

5.
6.
7.

[L１５] 利用できる計算機で、ブロック化を行った
行列-行列積のコードに対し、アンローリングを各ループについて施し性能を調査せよ。行列の大きさ（N）を変化させ、各Nに対して適切なアンローリング段数を調査せよ。
[L５] 身近にある計算機の、キャッシュサイズと、その構造を調べよ。
[L５] 身近にある計算機の、命令レベル並列性の実装の仕組みを調べよ。
[L５] 本講義で取り扱っていないチューニング手法
を調べよ。
136

2017年度

計算科学技術特論A