最新の気象予測～その理論と技術～
北海道大学地球惑星科学集中講義
札幌管区気象台 室井ちあし集中講義の内容• １日目– 気象予測の基礎– 気象とコンピュータ– 数値予報モデル• ２日目– 数値予報モデル（続）– データ同化– アンサンブル予報– （セミナー）北海道の気象災害リスクと防災気象
情報気象とコンピュータ

北海道大学地球惑星科学集中講義
札幌管区気象台 室井ちあしはじめに• 理論、実験を並び、シミュレーションは「第3の科学」• 近年は、データ解析は「第4の科学」と呼ばれる• 実験が困難な現象の解明や実験に時間や費用がかかりすぎる場合に、仮想的な実験手段として用いられる• 気象分野はとりわけ、シミュレーション、データ解析が盛んで、親和性が高い• 天気予報以外にも、気象研究に広く活用されている

気象分野でのコンピュータの利用• 予測モデル– 計算が複雑で量も多く、高速なコンピュータ、スー
パーコンピュータが利用されることが多い• データ解析– データ取扱量が多く、巨大なストレージ（記憶装置）が利用されることが多い• 可視化

予測モデル6

予測モデルを用いた研究パターン• 「できあい」の予測モデルを移植・利用して、結果を出す– 事例解析

• モデルの数値計算に関する研究• 物理プロセスに関する研究

日本の全球モデル• 気象庁GSM– 天気予報～気候予測ほか、様々な業務の基盤• 東大・JAMSTEC MIROC• 電脳倶楽部• 東大・JAMSTEC NICAM– 正二十面体型全球雲解像モデル• JAMSTEC MESG
– インヤン格子型全球・領域モデル

日本のメソモデル• 気象庁 NHM
– 気象研究所、気象庁で開発、研究でも広く利用されている• 名大 CReSS– 「地球シミュレータ」を契機に、雲解像モデルとして開発されている• 電脳倶楽部• 気象庁 asuca– NHMに代わる新モデル、今年春から実用化

世界のメソモデル• WRF– 米国のコミュニティモデル• ARPS– 米国オクラホマ大学で開発されたモデル• MM5– 米国NCARで開発されていたモデル• COSMO– ヨーロッパのコミュニティモデル

予測モデル研究の課題• 「できあい」のモデル– 効果的に結果が出せる一方、現代のモデルは巨大化していて、ブラックボックス化– 結果が正しいのかの検討が不十分になりがちスーパーコンピュータ

12スーパーコンピュータ• 普通のコンピュータと何が違うのか?
• 原理は同じ• クルマで言うところの「Ｆ１カー」• 「京」はスーパーコンピュータの中のスーパーコンピュータ → 事業仕分けの対象に

TOP500• スパコン間の性能比較ランキングで最も有名なもの• 年に２度、更新・発表される• Linpack （大規模密行列の直接法による解法プログラム）を性能指標とする– 長時間負荷をかけ耐久性が示せる、歴史的な価値がある一方、１つのプログラムで性能を示すと言えるのか？という問題もある• 消費電力といった環境を重視した“Green500” も別途存在し、年々競争が激しさを増している14スーパーコンピュータ性能ランキングの変遷○今回1位となった中国・国防科学技術大学（NUDT）のＴｉａｎｈｅ-２（天河２号）は、LINPACK実効性能33.86ペタFLOPS、実行効率61.6％。中国のスパコンが世界1位となるのは、2010年11月のTianhe-1A（天河1A号）以来、2回目。
○TOP10ランクイン状況としては、米国(5システム)が1位、中国及びドイツ(2システム)が2位、日本が1システム。
TOP100ランクインでは、米国が47システムと圧倒的優位な状況となっている。次いで、日本(10システム)が2位、英国(8システム)が3位、中国及びフランス(6システム)が4位、ドイツ及びインド(4システム)が6位と続いている。

平成24年6月
順位
システム名称
1 Sequoia
2 「京」 (K computer)
3 Mira4 SuperMUC
5 Tianhe-1A(天河1A号)
6 Jaguar
7 Fermi
8 JuQUEEN
9 Curie thin nodes
10 Nebulae(星雲)サイト
ベンダローレンスリバモア研
IBM
理研 計算科学研究機構(AICS) 富士通
アルゴンヌ研IBM
ﾗｲﾌﾟﾆｯﾂｽｰﾊﾟｰｺﾝﾋﾟｭｰﾃｨﾝｸﾞｾﾝﾀ（LRZ）IBM
天津スパコンセンタ
NUDT
オークリッジ研
Cray
Cinecaコンピューティングセンター IBM
ユーリヒ総合研究機構（FZJ）
IBM
フランス原子力庁
Bull SA
深圳スパコンセンタ
Dawning

12 Helios
14 TSUBAME2.0
18 Oakleaf-FX
36 BlueGene/Q
41 HA-PACS
70 Hitachi SR16000
73 Camphor
84 BX900

国際核融合エネルギー研究センタ
東工大学術国際情報センタ（GSIC)
東大情報基盤センタ
高エネルギー加速器研究機構
筑波大計算科学研究センタ
東北大学 金属材料研究所
京都大学
日本原子力研究開発機構(JAEA)

Bull SA
NEC/HP
富士通
IBM
Appro/Cray
日立
Cray
富士通

平成24年11月
Linpack
演算性能
国名 （テラFLOPS）

米
日
米
独
中
米
伊
独
仏
中

16,325
10,510
8,162
2,897
2,566
1,941
1,725
1,380
1,359
1,271

日
日
日
日
日
日
日
日

1,237
1,192
1,043
518
422
244
239
191

平成25年6月

順位
システム名称
1 Titan
2 Sequoia
3 「京」 (K computer)
4 Mira5 JuQUEEN
6 SuperMUC
7 Stampede
8 Tianhe-1A(天河1A号)
9 Fermi
10 DARPA Trial Subset

サイト
ベンダ
オークリッジ研
Cray
ローレンスリバモア研
IBM
理研 計算科学研究機構(AICS) 富士通
アルゴンヌ研IBM
ユーリヒ総合研究機構（FZJ）
IBM
ﾗｲﾌﾟﾆｯﾂｽｰﾊﾟｰｺﾝﾋﾟｭｰﾃｨﾝｸﾞｾﾝﾀ（LRZ）IBM
テキサス大学Dell
天津スパコンセンタ
NUDT
Cinecaコンピューティングセンター IBM
DOD国防高等研究計画局 IBM開発センタIBM

15 Helios
17 TSUBAME2.0
21 Oakleaf-FX
39 SGI Altix X
41 HIMAWARI
42 SAKURA
45 PRIMERGY CX400
51 HA-PACS
95 Hitachi SR16000
97 Camphor
100 Hitachi SR16000

国際核融合エネルギー研究センタ
東工大学術国際情報センタ（GSIC)
東大情報基盤センタ
電力中央研究所
高エネルギー加速器研究機構
高エネルギー加速器研究機構
九州大学
筑波大計算科学研究センタ
核融合科学研究所
京都大学
東北大学 金属材料研究所

Bull SA
NEC/HP
富士通
SGI
IBM
IBM
富士通
Appro/Cray
日立
Cray
日立

Linpack
演算性能
国名 （テラFLOPS）

米
米
日
米
独
独
米
中
伊
米

17,59016,325
10,510
8,162
4,141
2,897
2,660
2,566
1,725
1,515

日
日
日
日
日
日
日
日
日
日
日

1,237
1,192
1,043
582
518
518
460
422
253
251.7
243.9

順位
システム名称
1 Tianhe-2(天河2号)
2 Titan
3 Sequoia
4 「京」 (K computer)
5 Mira6 Stampede
7 JuQUEEN
8 Vulcan9 SuperMUC
10 Tianhe-1A(天河1A号)

サイト
ベンダ
国防科学技術大学
NUDT
オークリッジ研
Cray
ローレンスリバモア研
IBM
理研 計算科学研究機構(AICS) 富士通
アルゴンヌ研IBM
テキサス大学Dell
ユーリヒ総合研究機構(FZJ）
IBM
ローレンスリバモア研
IBM
ﾗｲﾌﾟﾆｯﾂｽｰﾊﾟｰｺﾝﾋﾟｭｰﾃｨﾝｸﾞｾﾝﾀ（LRZ）IBM
天津スパコンセンタ
NUDT

20 Helios
21 TSUBAME2.0
26 Oakleaf-FX
43 PRIMERGY CX400
47 SGI Altix X
50 SAKURA
51 HIMAWARI
62 HA-PACS
63 Cray XC30

国際核融合エネルギー研究センタ Bull SA
東工大学術国際情報センタ（GSIC) NEC/HP
東大情報基盤センタ
富士通
九州大学
富士通
電力中央研究所
SGI
高エネルギー加速器研究機構 IBM
高エネルギー加速器研究機構 IBM
筑波大計算科学研究センタ
Appro/Cray
国立天文台
Cray

Linpack
演算性能
国名 （テラFLOPS）

中
米
米
日
米
米
独
米
独
中

33,863
17,59017,173
10,510
8,587
5,168
5,009
4,293
2,897
2,566

日
日
日
日
日
日
日
日
日

1,237
1,192
1,043
621
582
537
537
422
420

（文部科学省資料）

TOP500（平成25年6月）のうち日本に設置されているスパコン
順位

サイト

ベンダ

システム名称

Linpack
演算性能（テラFLOPS）

4 理研 計算科学研究機構(AICS)

富士通「京」 (K computer)

10,51020 国際核融合エネルギー研究センタ
21 東工大学術国際情報センタ（GSIC)
26 東大情報基盤センタ
43 九州大学
47 電力中央研究所
50 高エネルギー加速器研究機構
51 高エネルギー加速器研究機構
62 筑波大計算科学研究センタ
63 国立天文台

Bull SA
NEC/HP
富士通
富士通
SGI
IBM
IBM
Cray
Cray

HeliosTSUBAME 2.0
Oakleaf-FX
PRIMERGY CX400
SGI Altix X
SAKURA
HIMAWARI
HA-PACS
Cray XC30

1,237
1,192
1,043
621
582
537
537
422
420

118 核融合科学研究所
121 京都大学
124 東北大学 金属材料研究所
129 分子科学研究所
138 サービスプロバイダ
161 日本原子力研究開発機構(JAEA)
177 サービスプロバイダ
207 九州大学
212 東大 物性研
273 エレクトロニクス関係
282 京都大学
318 サービスプロバイダ
334 金融関係
341 地球シミュレータセンタ
344 北海道大学情報基盤センタ
393 JAXA
432 北陸先端科学技術大学院大学
461 東大情報基盤センタ
467 東大ヒトゲノム解析センタ
480 理研情報基盤センタ(RIKEN)

日立
Cray
日立
富士通
HP
富士通
HP
富士通
SGI
IBM
Cray
HP
IBM
NEC
日立
富士通
Cray
日立
日立
富士通

Hitachi SR16000
Camphor
Hitachi SR16000
PRIMERGY CX250 & RX300
DL160 Gen8PRIMERGY BX900
BL460c Gen8PRIMEHPC FX10
SGI Altix ICE 8400EX
iDataPlex DX360M4
LaurelDL360p Gen
xSeries x3650M3
地球シミュレータ
Hitachi SR16000 Model M1
Fujitsu FX1Cray XC30
T2Kオープンスパコン
HA8000-tc/HT225
RICC

253
252
244
238
227
191
187
167
162
139
135
126
126
122
122
111
105
102
101
98

（文部科学省資料）

TOP500の各国1位の推移
世界で初めて10ﾍﾟﾀﾌﾛｯﾌﾟｽの壁を突破[FLOPS]
10ペタ

LINPACK性能

Tianhe-2
33.8PF

Tianhe-1A
RoadRunner 2.5PF
1.0PF

1ペタ

BlueGene/L
地球シミュレータ 70.7TF
35.8TF

100テラ｢京」
10.5PF

Titan
17.5PF

ヨーロッパ

アメリカ

日本
中国

10テラ

1テラ

100ギガ

（文部科学省資料）

TOP500

国別性能割合推移（文部科学省資料）

スーパーコンピュータ「京（けい）」の概要
・2011年6月と11月の二期連続で世界スパコン性能ランキング(TOP500)において1位を獲得
・ 「京」の利用研究が2年連続でゴードン・ベル賞（コンピュータシミュレーション分野での最高の賞）
を受賞
80 cm
○概要
206 cm◆平成23年11月にLINPACK性能※1 10ペタフロップス※2達成。
◆平成24年6月システム完成済（兵庫県神戸市の理化学研究所に設置）◆平成24年9月28日に共用開始

CPU（富士通製）
8万個以上を使用※1 スーパーコンピュータの性能を測るための世界的な指標（ベンチマークプログラム,※2 10ペタフロップス：一秒間に1京回（=10,000兆回＝1016回）の足し算，掛け算が可能な性能○プロジェクト経費 約1,110億円（H18~H24）

ラック

○特長
◆全ＣＰＵフル稼働時の連続実行時間は２９時間以上で世界最高水準の信頼性◆世界トップ10の実行効率(理論性能に対する実際の性能の比率）平均が78％のところ、「京」は93％◆アプリケーションプログラムの実行性能や使いやすさに関して高い性能◆水冷システムの導入により消費電力の削減や故障率の低減に寄与◆六次元メッシュ/トーラス結合の採用による高い利便性・耐故障性・運用性◆共用法に基づき、登録機関（高度情報科学技術研究機構）と理化学研究所が連携し、「京」を利用する体制を構築。
研究棟（文部科学省資料）

ＨＰＣＩの枠組み○ 「京」を中核とする国内のスパコンやストレージを高速ネットワークでつなぎ、ユーザー窓口の一元化などにより、利便性の高い利用環境を構築。
○ 「HPCIの整備・運営」として、各機関への委託事業により実施。
京大

阪大
スパコン運用北大
スパコン運用

情報学研究所・ユーザＩＤ管理システ
ム運用（シングルサイ
ンオン機能の提供）
・ＳＩＮＥＴ４運用

東北大
スパコン運用

スパコン運用
筑波大
スパコン運用「京」

東 大

九大
スパコン運用東工大

・スパコン運用
・共用ストレージ運用

スパコン運用

ＦＯＣＵＳ
アクセスポイントの設置・運用

ＲＩＳＴ
利用者選定

申
請

選
定

理 研
・「京」運用・共用ストレージ運用
・全体運営の企画調整

名大

ＲＩＳＴ

スパコン運用アクセスポイントの設置・運用

シングルサインオン
一つのアカウント
で全ての計算資源が利用可能

幅広いユーザ（文部科学省資料）

９大学情報基盤センターの計算リソースの概要
 大型計算機を運用管理するとともにその整備を図る
 学術研究等の共同利用に供する
 計算機の高度利用に関する研究および開発を行う

平成25年4月現在総理論演算性能 6,509Tflops
大阪大学：
SX-9 (16.4Tflops, 10TB)
SX-8R (5.3Tflops,3.3TB)
Express5800/120Rg-1 (6.1Tflops, 2TB)
Express5800/53Xh (16.6Tflops, 2.6TB)
京都大学：
Cray XE6 (300.8Tflops, 60TB)
APPRO GreenBlade8000 (242.5Tflops,
38TB)
APPRO 2548X (10.6Tflops, 24TB)
九州大学：
PRIMEHPC FX10 (181.6Tflops,24.6TB)
PRIMERGY CX400 S1 (811.9TF,185TB)
SR16000/L2 (25.3Tflops,5.5TB)
名古屋大学：FX1 (30.7Tflops,24TB)
HX600 (25.6Tflops,
10TB)
M9000 (3.84Tflops,3TB)

北海道大学：
SR16000/M1 (172.6Tflops,
22TB)
東北大学：
SX-9 (26.2Tflops, 16TB)
SX-9 (3.3Tflops, 2TB)
Express5800(1.7Tflops, 3TB)

筑波大学：
T2K-Tsukuba (95.4Tflops,21TB)
ﾌﾛﾝﾃｨｱ計算機ｼｽﾃﾑ (802Tflops,34TB)

東京大学：
T2K (140.1Tflops, 31TB)
SR16000/M1 (54.9Tflops,11TB)
PREMEHPC FX10 (1135.2Tflops,
150TB)
東京工業大学：TSUBAME2.0(2400Tflops,99TB)

2013年4月現在（文部科学省資料）

ＨＰＣＩ戦略プログラム戦略分野「京」を中核とするＨＰＣＩを最大限活用し、①画期的な成果創出、②高度な計算科学技術環境を使いこなせる人材の創出、③最先端コンピューティング研究教育拠点の形成を目指し、戦略機関を中心に戦略分野の「研究開発」及び「計算科学技術推進体制の構築」を推進する。

＜戦略分野＞

予測する生命科学・医療および創薬基盤分野１ ゲノム・タンパク質から細胞・臓器・全身にわたる生命現象を統合的に理解することにより、疾＜戦略機関＞
・理化学研究所

病メカニズムの解明と予測をおこなう。医療や創薬プロセスの高度化への寄与も期待される。

新物質・エネルギー創成

分野２ 物質を原子・電子レベルから総合的に理解することにより、新機能性分子や電子デバイス、更には各種電池やバイオマスなどの新規エネルギーの開発を目指す。

防災・減災に資する地球変動予測
分野３ 高精度の気候変動シミュレーションにより地球温暖化に伴う影響予測や集中豪雨の予測

・東大物性研（代表）・分子研・東北大金材研・海洋研究開発機構を行う。また、地震・津波について、これらが建造物に与える被害をも考慮した予測を行う。

分野４ 次世代ものづくり

先端的要素技術の創成～組み合わせ最適化～丸ごとあるがまま性能評価・寿命予測というプロ
セス全体を、シミュレーション主導でシームレスに行う、新しいものづくりプロセスの開発を行う。

物質と宇宙の起源と構造

分野５ 物質の究極的微細構造から星・銀河の誕生と進化の全プロセスの解明まで、極微の素粒子から宇宙全体に至る基礎科学を融合し、物質と宇宙の起源と構造を統合的に理解する。

・東大生産研（代表）・宇宙航空研究開発機構
・日本原子力研究開発機構

・筑波大（代表）
・高エネ研
・国立天文台※ スーパーコンピュータ「京」で、社会的・学術的に大きなブレークスルーが期待できる分野を「戦略分野」（５分（文部科学省資料）
野）とする。

スーパーコンピュータによって期待される成果の例（１）
心疾患のマルチスケール・マルチフィジックス
シミュレーション(研究代表者：東京大学・久田俊明)
細胞・組織・臓器を部分ではなく、心臓全体をありのままに 再現し、 心臓病の治療法の検討や薬の効果
の
評
価
に
貢献
心臓シミュレーション

地震・津波の予測精度の高度化に関する研究(研究代表者：東京大学・古村孝志、東北大学・今村文彦)

シミュレーションによる地震・津波の被害予測

50ｍ単位(ブロック単位)での予測から地盤沈下や液状化現象等の影響も加味した10ｍ単位(家単位)の詳細な予測を可能とし、都市整備計画への活用による災害に強い街
作 り や き め 細 か な
避難計画の策定等に貢献

創薬応用シミュレーション(研究代表者：東京大学・藤谷秀章)

新薬の候補物質を絞り込む期間を半減(約2年から約1年)して画期的な新薬
薬候補のタンパク質へのの開発に貢献
高精度結合シミュレーション

乱流の直接計算に基づく次世代流体設計システムの研究開発(研究代表者：東京大学・加藤千幸)

乱流の直接計算を工業製品の熱流体設計に適用することにより、従来行われていた風洞実験などを完
全
に
シ
ミ
ュ
レーションで代替し、設計車両挙動を解明する全乱流渦の 効 率 化 に 貢 献のシミュレーション（文部科学省資料）

スーパーコンピュータの課題• スケーリング• 消費電力の増大– 従来は、電力性能比を犠牲にしても、高速なマシンを求めてきた• 「京」・・15MW、一般家庭約3万～4万世帯分• 信頼性– 使用するパーツが増えるとともに、故障率の増大• プログラミング– 分割数が増えると、負荷分散が困難になり、メモリバッファが増大– 効率的なライブラリの整備• ＩＯ– 容量、消費電力、転送速度、信頼性・・

24高速なプログラミング

プログラミング言語• Fortran90– Fortran 95/2003, Co-array Fortran• C. C++• Ruby• Perl, javascript• GPUコンピューティグ
– CUDA– OpenACC

Fortran の特徴• 手続型コンパイル言語– 計算機言語の保守本流• 数値計算向け– 名前の由来は Formula Translation• 世界最古の高水準言語– 膨大なプログラム資産– ベクトル化、並列化などのコンパイラ技術の資産

Fortran 90： 現在• 規格では FORTRAN から Fortran になった
• 近代的な制御構造の導入（GOTO が不要）
• モジュールの導入（COMMON が不要）
• 配列機能の強化（並列プログラミング）• 内部副プログラムの導入• 構造型、ユーザ定義型、ポインタの導入• 自由形式（キーボードとディスプレイ向き）

配列演算： 添字の並べ方• Fortran の配列 a(i, j) はつぎのようにメモリ上に配置される (i = 1..10, j = 1..5)a(1, 1)a(2, 1)a(3, 1)a(4, 1)a(5, 1)a(6, 1)a(7, 1)a(8, 1)a(9, 1)a(10, 1)a(1, 2)a(2, 2)a(3, 2)a(4, 2)a(5, 2)a(6, 2)a(7, 2)a(8, 2)a(9, 2)a(10, 2)a(1, 3)a(2, 3)a(3, 3)a(4, 3)a(5, 3)a(6, 3)a(7, 3)a(8, 3)a(9, 3)a(10, 3)a(1, 4)a(2, 4)a(3, 4)a(4, 4)a(5, 4)a(6, 4)a(7, 4)a(8, 4)a(9, 4)a(10, 4)a(1, 5)a(2, 5)a(3, 5)a(4, 5)a(5, 5)a(6, 5)a(7, 5)a(8, 5)a(9, 5)a(10, 5)

– 二重の DO ループでは、変数 j の繰り返しを外側に、変数 iの繰り返しを内側にすれば、アクセスがメモリ上で連続になって高速（C 言語の a[j][i] とは i と j の位置が逆）高速なプログラム• プログラムは物理法則、数式通りに書くことが基本• しかし、速く計算する、ということも重要– 日々の天気予報はいうまでもないが– 研究目的でも迅速性は求められる• 高速なコンピュータを使うことのほかに、最適化や並列化といった工夫が行われる

コンピュータの性能を決める要素• CPU の演算速度• メモリバンド幅• ノード間（プロセス間）の通信速度• 近年では CPU の演算速度の向上やメモリ容量の増大に比べてメモリアクセスや通信の速度の向上が遅れている– 一回の通信あたりの演算量が多いほうが有利

最適化• アルゴリズムやデータ構造、ソースコードの書き方を工夫して高速化すること• アルゴリズムやデータ構造を工夫してメモリアクセスを工夫すること• 計算時間を多く消費する「ホットスポット」を見つけることがスタート– 実行時間の計測、プロファイラの利用

最適化実例
計算順序の入れ替え• 例： 3次式
a * x ** 3 + b * x ** 2 + c + x + d((a * x + b) * x + c) * x + d• 数学的には同じ• 四則演算回数は上の式は１０回、下の式は６回、よって下の式の方が高速

最適化レベルとおもな項目（その1）• レベル 0, 3, 4

– べき乗の乗算化
a = b ** 3 ⇒ a = b * b * b– 偶数の乗算にシフト命令を適用• レベル 3, 4
– 共通部分式の削除
x=a*b
y=a*b*c⇒y=x*c

最適化レベルとおもな項目（その2）• レベル 3, 4

– 不変式のループの外への移動
do i = 1, 100x=a+b
y = c + delta(i)end dox=a+b
do i = 1, 100
y = c + delta(i)end doa + b がループ中で不変なら
t=a+b
do i = 1, 100
x=t
y = c + delta(i)end dox への代入がループ中だけなら

最適化レベルとおもな項目（その3）• レベル 4

– 短いループの展開： 繰り返し回数がコンパイルのときに確定している場合のみ

do j = 1, 100
do i = 1, 3a(i, j) = b(i) * c(j)end doend dodo j = 1, 100a(1, j) = b(1) * c(j)a(2, j) = b(2) * c(j)a(3, j) = b(3) * c(j)end do最適化が逆効果になることがある• ループ内の不変式がめったに実行されないIF 文の条件下にある場合
do j = 1, 10000
do i = 1, 10
if (a(j) > 0.0) thena(j) = sin(i * 2.0)
end ifend doend do

最適化によって i のループにおける不変式
sin(i * 2.0) の計算が i のループの外に移動された場合、j のすべての繰り返し（10000回）でこの計算が必要になる。
もし IF 文の条件 a(j) > 0.0 が 1度も成り立たなければ、元のコードでsin(i * 2.0) の計算は不要だった。

メモリアクセスの局所性
行列ベクトル積 y = Ax, y_i = sum(A_i,j * x_j)
• オリジナルコード
do i = 1, n
do j = 1, ny(i) = y(i) + a(i,j)*x(j)end doend doa(I,j) のアクセスがとびとびになる改良版• 連続アクセスに改良
do j = 1, n
do i = 1, ny(i) = y(i) + a(i,j)*x(j)end doend doa(I,j) のアクセスが連続になったが、y(i) が毎回呼び出され、書き込まれてしまう

さらに改良• アンローリングに改良
do j = 1, n, 2
do i = 1, ny(i) = y(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
end doend doa(i,j)は連続アクセスで、y(i)の読み書きは半減！j が偶数か奇数かで場合分けが必要
40

並列化の概要• 計算を複数に分散させる– 1 つあたりの計算量を少なくすることにより処理を高速化– 大容量のメモリを利用できるようになる– 全体の計算量は変わらないか、オーバーヘッド
逆に増える• ２つの手法（負荷が効率よく分散するように選択される）– データパラレル– タスクパラレル並列計算機と多段構造
MemoryL3

Bus Controller
L1

CPU
L2

L1

CPU

MemoryBus Controller
MemoryCPU: 1.9 GHz/ 2.1 GHzPOWER5Bus Controller

MCM: Multi Chip
ModuleBus Controller
MCM

MemoryMemoryBus Controller
MemoryBus Controller
MemoryBus Controller
MemoryBus Controller
MCM

1 node: POWER5 processor x 16

並列化の種類• 共有メモリ並列化（スレッド並列）– コンパイラによる自動並列化（要素並列化）– コンパイラに指示行（ディレクティブ）を与えるOpen MP• 分散メモリ並列化（プロセス並列）– ノード間の通信を明示的に指示する Message
Passing Interface (MPI) が標準

並列化（１）ノード内• プログラム中の DO ループを複数のスレッドに分割し複数のプロセッサで並列に実行– 繰り返しの順序に依存関係がない場合のみ

このループ内の計算は繰り返しの順序に依存しない

do i = 1, 100
x(i) = a * y(i) + b
end do4 つのスレッドに分割する場合

do i = 1, 25
x(i) = a * y(i) + b
end dodo i = 26, 50
x(i) = a * y(i) + b
end dodo i = 51, 75
x(i) = a * y(i) + b
end dodo i = 76, 100
x(i) = a * y(i) + b
end doプロセッサ 1プロセッサ 2プロセッサ 3プロセッサ 4

データパラレル・分散メモリ並列化• モデルの計算領域を東西方向と南北方向に 2 次元分割

North– 分割した各領域 Rnn を各プロセスで独立に計算– 移流など各領域間のやり取りは「のりしろ」(interface) の変数を
MPI により通信

R02

R12

R22

R01

R11

R21

R00

R10

R20

South
West

computational domain

Eastinterface並列化（２）ノード外
program mpi_ex ! 超カンタンなプログラム例
include 'mpif.h'! おまじない
integer :: myrank, idata, ierrinteger :: status(MPI_STATUS_SIZE)

コンパイルの例$ mpif90 mpi_ex.f90 –o mpi_ex
CPU-0

MPI_Send で データ送信
call MPI_Init(ierr)! おまじない
call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ierr)(idata = 1)
if (myrank == 0) then
! CPU（ノード）ランクが0なら
MPI_Recv で データ受信
idata = 1
call MPI_Send(idata, 1, MPI_INTEGER, 1, 1234, &CPU-1& MPI_COMM_WORLD, ierr) ! データを送信
write(6,*) myrank, idata
else if (myrank == 1) then
! CPU（ノード）ランクが1なら受信送受信の手続きを
call MPI_Recv(idata, 1, MPI_INTEGER, 0, 1234, &ソースコードに書く!!
& MPI_COMM_WORLD, status, ierr)write(6,*) myrank, idata
実行方法：mpiexec –n 2 ./mpi_ex
endifcall MPI_Finalize(ierr)! おまじない
実行結果
end program（CPU-0の出力）
0
11
2010年6月4日

1
No. 46（CPU-1の出力）

大きなプログラムを高速で動かすには
要素並列化・OpenMPとMPI通信を上手に組み合わせる技術が必要

分散ﾒﾓﾘ並列MPI
通信

SMPノード
要素並列化
OpenMP

SMPノード
要素並列化
OpenMP

SMPノード
要素並列化
OpenMP

SMPノード
要素並列化
OpenMP2010年6月4日

No. 47

分散メモリ並列化： 出力専用プロセス
communication and output
by 0-th proc.
0

1

2

communication and output by 0-th
proc. (for I/O only)

output byeach proc.

3

0

1

2

3

0

1

2

3

時間
プロセス
computation

communication

output数値計算上の問題

計算機による実数の表現• 数学的に同じと、数値計算の結果が同じ、は意味が違う• コンピュータは０と１しか扱えない• 数学の実数は連続だが、計算機で扱うことができるのは有限桁の浮動小数点数• ほとんどの計算機が IEEE 754 という規格に定められた単精度と倍精度の表現を採用– 単精度（32 ビット）
• REAL(4)– 倍精度（64 ビット）
• REAL(8) または DOUBLE PRECISIONIEEE 754 が定める浮動小数点数• 規格は IEEE Std 754-1985 “IEEE Standard for
Binary Floating-Point Arithmetic”– 単精度（32 ビット）• 符号 1 ビット、指数部 8 ビット、仮数部 23 ビット– 倍精度（64 ビット）• 符号 1 ビット、指数部 11 ビット、仮数部 52 ビット– つぎの例外に対する割り込みの有無を指示できる• 無効な演算、ゼロによる除算、オーバーフロー、アンダー
フロー、不正確IEEE 754 が定めるゼロ、無限大、非数• ゼロ– 指数部と仮数部のビットがすべて 0• 無限大– 指数部のビットがすべて 1– 仮数部のビットがすべて 0• 非数 (NaN)– 指数部のビットがすべて 1– 仮数部のビットがすべて 0 ではない

丸め誤差
• 浮動小数点数は実数を有限の長さの 2進数で近似するため、一般に誤差を伴う– 丸め誤差 ＝ 浮動小数点値 － 真の値• 有効数字は有限– 10進数の 0.1 を 2進数で表すと循環小数0.000110011001100 ...
になる⇒ 10倍しても 1 にならないかも0.1 の 10.0 倍 == 1.0 ?（4バイト実数）
program eqzero4
implicit nonereal(4) :: val4
val4 = 0.1
write(*, *) 'val4 = ', val4
if (val4 == 0.1) thenwrite(*, *) 'val4 == 0.1'
else
write(*, *) 'val4 /= 0.1'
end ifif (val4 * 10.0 == 1.0) thenwrite(*, *) 'val4 * 10.0 == 1.0'
else
write(*, *) 'val4 * 10.0 /= 1.0'
end ifend program eqzero4

左のプログラム eqzero4.f90 を SR11000 の
日立最適化 FORTRAN でコンパイルして実行する$ f90 -o eqzero4 eqzero4.f90$ ./eqzero4
val4 = 0.100000001
val4 == 0.1
val4 * 10.0 == 1.0

0.1 の 10.0 倍 == 1.0 ?（8バイト実数）
program eqzero8
implicit nonereal(8) :: val8
val8 = 0.1
write(*, *) 'val8 = ', val8
if (val8 == 0.1) thenwrite(*, *) 'val8 == 0.1'
else
write(*, *) 'val8 /= 0.1'
end ifif (val8 * 10.0 == 1.0) thenwrite(*, *) 'val8 * 10.0 == 1.0'
else
write(*, *) 'val8 * 10.0 /= 1.0'
end ifend program eqzero8

左のプログラム eqzero8.f90 を SR11000 の日
立最適化 FORTRAN でコンパイルして実行する$ f90 -o eqzero8 eqzero8.f90$ ./eqzero8
val8 = 0.100000001490116119
val8 == 0.1
val8 * 10.0 /= 1.0

実定数 0.1 は単精度実定数 0.1e0 と同じ

0.1 の 10.0 倍 == 1.0 ?（8バイト実数）
program eqzero8d
implicit nonereal(8) :: val8
val8 = 0.1d0
write(*, *) 'val8 = ', val8
if (val8 == 0.1d0) thenwrite(*, *) 'val8 == 0.1'
else
write(*, *) 'val8 /= 0.1'
end ifif (val8 * 10.0d0 == 1.0d0) thenwrite(*, *) 'val8 * 10.0 == 1.0'
else
write(*, *) 'val8 * 10.0 /= 1.0'
end ifend program eqzero8d

左のプログラム eqzero8d.f90 を SR11000 の
日立最適化 FORTRAN でコンパイルして実行する$ f90 -o eqzero8d eqzero8d.f90$ ./eqzero8d
val8 = 0.100000000000000006
val8 == 0.1
val8 * 10.0 == 1.0

倍精度で扱うためには 0.1d0 とする

情報落ち（積み残し）• 絶対値が大きく異なる実数の和を計算すると、小さい数が無視される–例
1000000.0 + 0.0000001⇒ 浮動小数点数： 丸め誤差 ?????? を伴う1000000.0??????
+)
0.0000001??????
1000000.0??????

+ 0.0000001 は丸め誤差に埋没桁落ち• 絶対値がほぼ等しい実数の差を計算すると、結果の有効数字が短くなる–例
1000000.1 - 1000000.0⇒ 浮動小数点数： 丸め誤差 ?????? を伴う1000000.1??????
-)1000000.0??????
0.1??????

結果の有効数字が 1 桁に

桁落ちの実例• 2次方程式の解はax  bx  c  02

 b  b  4acx
2a2• もし b と b 2  4acの絶対値がほぼ等しいと、
複号 ± のうちどちらかで桁落ちが起こる桁落ちを避けるために？
• 桁落ちを避けるために、絶対値の大きい解をb  0 ならば b  b  4ac
x1 
2a
2

 b  b  4ac
x1 b  0 ならば
2acから計算し、もうひとつの解を x2 
ax1
によって計算する
2

気象計算と桁落ち• 微分方程式、微分で表現される物理量（渦度など）が物理法則にあらわれる気象計算では、桁落ちの問題がよくあらわれる• 平均場はあらかじめ引いておいて、あとから加える、といった処理が必要

計算順序の入れ替えの影響• 例： 3次式
a * x ** 3 + b * x ** 2 + c + x + d
において、べき乗を乗算に置き換え、乗算の回数を減らすため、数学的には等価な((a * x + b) * x + c) * x + dに変更すると結果が変わる可能性がある• 最適化によって計算順序が変わったときには注意（四則演算回数は上の式は１０回、下の式は６回、よって下の式の方が高速で、かつ下の式の方が情報おちが少ない）データ解析、可視化

気象分野の特徴• 非定常な問題を扱うことが多く、空間３次元・
時間１次元の膨大なデータを出力し解析する
ことから、ビッグ・データの問題が常につきまとう• データ圧縮技術、優れたユーザーインター
フェース、人間が理解しやすい応用処理（可視化など）が求められる

データセットのこと• 効率的な数値予報データ取り扱いのため、様々なデータ形式が提案されている– GRIB, GRIB2: 国際交換で標準的に用いられる– netcdf: 米国でメジャー– nusdas: 気象庁の標準

可視化• Grads• 他になにかないのか?
• 個人的にはもう少しなんとかしたい・・GPUコンピューティング

67

TSUBAME 2.0の4000 GPUを用いた次世代気象モデルの大規模高性能計算

下川辺 隆史東京工業大学
創造エネルギー専攻（学術国際情報センター）
第８回ASE研究会（東京大学
Dec. 6, 2010情報基盤センター）

What’s GPU ?
 Graphics Processing Unit

 もともと PC の3D描画専用の装置
 パソコンの部品として量産されてる。

= 非常に安価

3D GameComputer Graphics

GPU

http://www.nvidia.co.jp

GPGPU
http://gpgpu.org/

 General Purpose computation on GPU

汎用GPU計算、GPUコンピューティング
 数値流体力学（CFD）.

 N体問題 高速フーリエ変換（FFT）
 …

Graphics

CFD

 プログラムはGPUむけの開発言語・環境を用いる
 CUDA (NVIDIA) ATI Stream (AMD) OpenCL (Khronos Group)

GPUによるHigh Performance Computing が現実に。

Peak Performance [GFLOPS]

1500

CPU and GPUの演算性能

GPU

1250

1000
750
500
250
0

CPU
Date

71

http://developer.nvidia.com/object/cuda_3_2_downloads.html

Memory Bandwidth [GByte/s]CPU
andGPUのメモリバンド幅
200
GPU

160
120
80

40

CPU

0

Date
格子計算のアプリケーションでは多くの場合、演算性能よりもメモリバンド幅が重要
72

WRFとASUCAの高速化のアプローチ
の違い
 WRF
GPU
Acceleration
Accelerator Approach
Initial condition

Dynamics

Physics

OutputCPU
GPU

Full GPU Approach ASUCA GPU Computing
Initial condition

数十倍の高速化の実現

Dynamics

Physics

Output

CPU

GPU
73

マルチGPU計算：境界領域のデータ交換 MPIを用いたGPUとCPUによるデータ交換
境界領域

GPU

Node 1

Node 2(1) GPU→CPU(3) CPU→GPUCPU(2) CPU → CPUyx

GPUは他のGPU上のメモリへは直接アクセスできない

74

ASUCAによる台風計算の例
4792 x 4696 x 48 mesh (水平解像度 500 m)
437 GPUs on TSUBAME 2.0

75

ASUCA: Fortran から CUDA へ フル GPU アプリケーション
 ゼロから書き換え

Fortran

C/C++

気象庁における
配列の順序の交換
オリジナルコードCUDA
GPU コード

３次元配列の要素順序
z,x,y (k,i,j)-ordering

x,z,y (i,k,j)-ordering

x,z,y (i,k,j)-ordering

GPUコードでのメモリアクセスパフォーマンスを向上

76

TSUBAME 2.0 ベンチマークの日々
00:00

06:00

12:00

LINPACK by Endosan

10/14
10/15

18:00

HPL in HPCC
by Nukada-san

ASUCA

1.134PFHPL inHPCC
1.140PF

10/16

Full HPCC

10/17
1.023PFASUCALINPACK by Endosan

ASUCA

LINPACK by Endosan

Crash!!
ASUCA

10/18

LINPACK by Endosan

1.192PFWe can use all nodes only for four days.

24:00

自動GPU化する OpenACC が登場、しかし遅い

78

Fortran に指示行を入れて、自動CUDA化、OpenMP化するコンバータを共同開発中

79