4. 各会合での情報提供資料
本 WG の各会合にて、以下の情報提供を実施した。次ページ以降に資料を掲載する。
No

実施日

タイトル

1

2011/08/23

2

2011/10/28

3
4

2011/10/28
2012/03/01

5
6

2012/03/01
2012/07/23

7

2012/07/23

8
9

2012/07/23
2013/02/22

Japan's Grand Challenge Activities inLife Science
多倍長計算の紹介と K-computer 向け固有値
ソルバ「Eigen-K」について
高精度線形代数演算ライブラリMPACKの紹介
宇宙プラズマでのペタスケール数値計算(惑
星磁気圏シミュレーション)
ペタスケールのアプリケーションへ向けて
JAXA における数値計算ライブラリの利用状
況と OPL 適用評価
並列擬似乱数発生アルゴリズムとペタスケー
ル時代のモンテカルロ法の展望
FFTE について
フラグメント分子軌道（FMO）法と並列 FMO
プログラム OpenFMO について

情報提供者 (敬称略,所属は当時)
姫野 龍太郎 (理化学研究所)
今村 俊幸 (理化学研究所計算科学研究機構)
中田 真秀 (理化学研究所)
深沢 圭一郎 (九州大学)
町田 昌彦 (日本原子力研究開発機構)
松尾 裕一 (宇宙航空研究開発機構)
三浦 謙一 (国立情報学研究所)
高橋 大介 (筑波大学)
稲富 雄一 (九州大学)

4.1 Japan’s Grand Challenge Activities in Life Science

Japan’s Grand Challenge Activities in Life Science
日本における生命科学分野でのグランドチャレンジ
Computational Science Research Program, RIKEN
理化学研究所、次世代計算科学研究開発プログラム、副プログラムディレクタ
Ryutaro Himeno
姫野龍太郎

概要
次世代スーパーコンピュータ京を使って、その性能がフルに発揮できるソフト
ウェアを開発し、それ以前では不可能だった計算を可能にするための生命科学分野でのグランドチャレンジを 2006 年 10 月から理化学研究所を拠点に開始した。このプロジェクトには全国の１３機関から約２００名の研究者が参加して2013 年３月末まで行われた。この間、分子スケール、細胞スケール、臓器全身スケール、大規模生命化学系データ解析、脳神経系など多岐にわたるアプリケ
ーションプログラムと、大規模可視化ツールや並列化ツールなど、合計３４本を開発していた。これらの開発状況と並列性能を紹介した。これらのソフトのうち、この時点では１１本が京の上で実行できるようになり、十分な並列性能を達成していた。さらに５本は京の上で非常に良い性能を出していた。特に分子動力学プログラム cppmd は京で１PFLOPS を超える性能を達成した。

Japan’s'Grand'Challenge'
Ac2vi2es'in'Life'Science
Ryutaro'Himeno'
Next%Genera*on,Integrated,Simula*on,of,Living,Ma9er,
Computa*onal,Science,Research,Program,
RIKEN

Schedule of the Next-generation supercomputer
R&D project
We are here.
FY2006

Buildings

Applications

SystemNext-Generation
Integrated
Nanoscience
Simulation

Next-Generation
Integrated
Life Simulation

Computerbuilding
Researchbuilding

FY2007

Conceptual
design

FY2008FY2009FY2010FY2011

Prototype,
Prototype,evaluation
evaluation

Production,
Production, installation,
installation,andand adjustment
adjustment

Development,
Development, production,
production, andand evaluation
evaluation

Verification
Verification

Detailed
Detailed design
designDevelopment,
Development, production,
production, andand evaluation
evaluation

DesignDesignDesignDesignFY2012TuningTuning andandimprovement
improvement

Verification
Verification

Construction
Construction
Construction
ConstructionTotal Simulation of Living Matter

Goal
44Understand the life phenomena: reproducing life program on computer
Develop new medicine and medical equipment

Approachcombination of analytical approach and Data-driven approach on Peta-scalecomputer
Data-Based•Understand theMoleculeapproach tolife phenomena:
scale
Analyticalanalysis
reproducing life
approachprogram on
Cell
Peta
scale
simulation
computerscale•Develop new
New
medicine andmethodology
Organ/Whole
Fusion of Data
medicalbody scale
based approach andequipment
Ana5ytical approach
Basic principal(Equations)

Experimental data

List'of'developing'so?ware'1/2

9

1
6



List'of'developing'so?ware'2/2
9

5

4



Current'status'of'developmentPhase'I'<'256

Phase'II'<'8192

Phase'III'on'K

SPHERE'LSV'
cppmd'IOSSIM'
NeuroMorphoKit'
VSM'
CMDN'
NEST'
MEGADOCK'LiSDAS'SiGNHSSM'
SiGNHL1'
SiGNHBN'ExRAT'
NGS'analyzer'
ParaHaplo'
UTHeart'ZZHLUNG'HIFU'ZZHDOSE'ZZHRBC'ZZHEFSI'
RICS'
PlatypusHCGM/MM'
PlatypusHQM/MM'
PlatypusHQM'
ProteinDF'
CafeMol'
MARBLE'
PlatypusHREIN'
PlatypusHMM/CG'

H1''

cppmd'>'1'Peta

Goal:'>'1'Pata'FLOPS
1'

11/31'are'running'on'K

2''

5''



8''

1''
SPHERE'LSV'

cppmd'IOSSIM'

NeuroMorphoKit'

VSM'

CMDN'

NEST'

MEGADOCK'LiSDAS'SiGNHSSM'

SiGNHL1'

SiGNHBN'ExRAT'

NGS'analyzer'

1,000,000''

ParaHaplo'

UTHeart'ZZHLUNG'HIFU'ZZHDOSE'ZZHRBC'ZZHEFSI'

RICS'

PlatypusHCGM/MM'

PlatypusHQM/MM'

PlatypusHQM'

ProteinDF'

CafeMol'

MARBLE'

PlatypusHREIN'

PlatypusHMM/CG'

Maximum'No.'of'parallel'execu2on'
16/31'>'8000'parallelized

100,000''

10,000''1,000''

100''

10''



Parallelized'ra2o'(1)

100.000''

11/31'>'99.75%

99.950''99.900''99.850''99.800''99.750''99.700''99.650''99.600''99.550''99.500''Parallelized'ra2o'(2)
100.000''

3/31'>'99.99%

99.990''99.980''99.970''99.960''99.950''99.940''99.930''99.920''LSV'SPHERE'

cppmd'IOSSIM'

VSM'

NeuroMorphoKit'

NEST'

CMDN'

LiSDAS'MEGADOCK'SiGNHL1'

SiGNHSSM'

ExRAT'SiGNHBN'

NGS'analyzer'

UTHeart'

ParaHaplo'

HIFU'

ZZHLUNG'ZZHRBC'ZZHDOSE'RICS'

ZZHEFSI'

PlatypusHCGM/MM'

PlatypusHQM'

PlatypusHQM/MM'

CafeMol'

ProteinDF'

MARBLE'

PlatypusHREIN'

99.900''

PlatypusHMM/CG'

99.910''

Strong'scaling'shown'by'cppmd

+!"( 
0)"-+#/+).   
0/+##-!&0.2 2)-$&*0)"-+#,&-.,"-/+)   &*)-$&*
 +0*/.#+- 
./",.    
+,+(+$1 +# +),0//&+*3 3
*+!"0)"-+#*+!".

0)"-+# +-".
 
%"+-"/& (,"',"-#+-)* "  ( 0(/&+* /&)" #+-  ." ./",0./&*"!,"-#+-)* "
##& &"* 1 Measured'At'the'end'of'March,'2011 Library'used'in'our'developing'
so?ware'package
FFTW
Marbleo

Platypuso

NetCDF

HDF

Lapack,

ScaLapack GSL

o

o

o'

o

ProteinDF

others

o

RICS

o

BoostHIFU

SHEREZZHEFSI

SHERE

Megadock

o

CMDN

o

IOSSIM

Sundail'InterView

EuroMorpho'
Kit

o

o

VSM

o

o

GD,'zlib

LSV

Libxml2,'zlib

SHERE

Libxml2,'zlib

11'

Molecular'Dynamics'Simula2on'Code:'cppmd
• Designed'for'K'computer'• Achieved'about'50%'of'theore2cal'peak'performance'

of'K

12'

Molecular simulation of multi-drug
efflux transporter AcrB

Platelet'and'Blood'clohng'simula2onHeart'Simula2on

Neuron'simula2on'(silkworm)

Summary• We are developing 34 application softwarefor Life Science researchers• 31 applications that are highly tuned for K will
be offered to use on K by the end of FY2012.
• 16 out of 31 have been checked its parallel
performance on clusters.
• 5 applications shows good parallel
performance• cppmd: MD simulation core program has
already achieved more than 1 Peta FLOPsContributors'

Workshop'in'2009'

Workshop'in'2008'18'

Organiza2on'13'(RIKEN,'U.'of'Tokyo,'Osaka'U.,'Kyoto'U.'Tohoku'U.,'
Tokai'U.,'JAIST,'Chiba'U.,'Keio'U.'Yokohama'City'U.,'
Kobe'U.,'TiTch,'ISM)''
About'200'person,'including'60'postDoc'
19'

4.2 多倍長計算の紹介と K-computer 向け固有値ソルバ「Eigen-K」について

多倍長計算の紹介と K-computer 向け固有値ソルバ「Eigen-K」について
報告者：今村俊幸（電気通信大学大学院 情報理工学研究科，
現理化学研究所計算科学研究機構）
共同研究：山田進, 町田昌彦 (JAEA)
１．

多倍長のソフトウェア開発

ペタスケールの計算機を数週間稼働させたとき、その演算回数は 10^21 に達する。理論上、浮動小数点計算の蓄積誤差はその有効精度を超えてしまう。実際には丸め誤差は相当数がキャンセルするため、多くの場合はシビアではない。しかし、ペタスケール計算機で解くべき問題規模においては、問題そのものの表現が倍精度で十分かの議論も登場しつつある。
そのような背景の中で、我々は double-double データフォーマットに基づく多倍長数値計算
ソフトウェア(multiple-precision numerical software)の研究開発を進めている。Knuth 以
来、多倍長計算の基本アルゴリズムは確立しているが、本研究で採用するのは Bailey, Hidaらによる double-double である。通常の多倍長計算が integer を多数並べて実現するところを, double によって実現するものである。基本演算が double の演算で構成できるため, マ
イクロプロセッサの持つ高い浮動週数点演算性能を活かして高速計算が可能という利点を持つ。Double-double の基本演算自身は Bailey, Hida らによって qd ライブラリの一部という形で公開されている。また、数値計算ライブラリも中田により MPACK の一部に収録されている（MPACK 自身は double-double 以外にも多くの多倍長フォーマットをサポートしている）。
我々は、DD(double-double)に基づく基本線形計算ルーチンの整備から始め、実機向けに最適化を施し実用レベルまで性能を引き出すことを考えている。次に、MPI 通信ライブラ
リ、更に高度なソフトウェア群（LAPACK, ScaLAPACK）へと開発を進め、大規模なスパ
コン環境下での実用的な多倍長数値計算ライブラリの展開を進める計画でいる。
この報告書を作成している現在（平成２５年４月）、我々は BLAS として知られている線
形 基 本 演 算 サ ブ プ ロ グ ラ ム 群 の double-double 版 QPBLAS の 公 開 を 開 始 し て い る（http://ccse.jaea.go.jp/ja/download/qpblas.html）。さらに、東大奥田洋司教授の協力を得て GPU 版も公開をしている(http://ccse.jaea.go.jp/ja/download/qpblas_gpu.html)。また、次章で紹介する「京」向け固有値ソルバ EigenK についても double-double による実装
QPEigenK が完成し、その速報が SC2012(Salt Lake City, USA)にてポスター報告されている(T.Imamura 他, Preliminary Report for a High Precision Distributed Memory Parallel

Eigenvalue Solver, SC12, Nov. 2012)。

２．

京コンピュータ向け固有値ソルバ「Eigen-K」

我々が地球シミュレータ時代から開発を進めてきた固有値ソルバーを、2008 年ごろから「京」コンピュター向け（開発初期は次世代計算機と呼んでいた）に改良をした Eigen-K が昨年
６月に JAEA より公開されている(http://ccse.jaea.go.jp/ja/download/eigenk.html)。報告資料にもあるように、MPI/OpenMP ハイブリッドモデルによる２階層の並列処理が可能であり、ベクトル化もしくは SIMD による性能チューニングも実施されている。ベンダー提供の高性能 BLAS を最大限に利用することで、高い性能を発揮できるように設計されている。
開発当時は「京」コンピュータ実機が利用できなかったため、東大情報基盤センターの T2K
スパコンや FX10 を使った開発体制をとっていたが、平成２４年９月の一般共用開始により「京」上でのチューニングや機能強化が進んでいる（「京」上での機能強化版の公開は近日
を予定）
。
現在、
「京」から「ポストペタ」に向けた固有値ソルバへの開発を指向して、JST CREST「ポストペタスケール高性能計算に資するシステムソフトウェア技術の創出」領域の中の「ポストペタスケールに対応した階層モデルによる超並列固有値解析エンジンの開発」課
題が進行している。我々は当該課題において Eigen-K を基にした Eigen-Exa の開発を進めている。Eigen-Exa は「京」以上に高並列化が進むと予想される次世代計算環境を想定して、メニイコア向けの階層的なアルゴリズム構成、省通信最適化、自動チューニングなど
次世代計算機環境向けの研究技術が盛り込まれた数値計算ライブラリになる予定である。

‫ᧈ̿ٶ‬ᚘምỉኰʼể
K-computerӼẬ‫׍‬ஊ͌ἏἽἢ
ẐEigen-KẑỆếẟề
ᩓൢᡫ̮‫ܖٻܖٻ‬ᨈᴾ ऴ‫ܖ߻ྸإ‬ᄂᆮᅹʻ஭̢࠳σӷᄂᆮ :ᴾ ‫ޛ‬ဋᡶ, ထဋଡࢠᴾ (JAEA)
2011/10/28

SSᄂ

1

ᾀᢿᾉᴾ ྸᛯὉᚘምඥỉᛇኬỊɶဋ൞Ệấ˓ẶẲộẴẇ

‫ᧈ̿ٶ‬ỉἏἧἚỸỹỹỴ᧏ႆ
2011/10/28SSᄂ2‫ᙹٻ‬೉҄ỆợỦችࡇồỉࢨ᪪ ̿ችࡇᚘምỉᄊዧờᡈẟὉὉὉ இ‫ٻ‬Ⴛ‫ݣ‬ᛚࠀᾉFrank Matrix κ(A)=O(N^2)ẆᾁỴἽἆἼἌἲỂ

2011/10/28

SSᄂ

310ɢឬἅỴồỉ‫ޒ‬ஓ ࣱᏡˌ‫ٳ‬ᾉችࡇ N=105ửឬảỦểẆவˑૠầ O(N㻞)ỉբ᫆ỂችࡇɧឱầႆဃẴỦẇẮủỊ̿ችࡇỂỉᘙྵỆؕỀẪஜឋႎễբ᫆ίᡦɶỉᚘምඥોծỂ˴ểẦểẟạἾἫἽỂỊễẟὸ
̊ᾉFrank matrix ỺỿἇἋἃὊἽỂỊN=10^6ầἑὊἄἕἚ

Æ᭗ᡮ҄ờ᣻ᙲẻầችࡇỉբ᫆ờvisibleỆ ἅỴỉஊјМဇ ί̊ảịὸ࠹ếẦỉἅỴỆ‫ᧈ̿ٶ‬ᚘም QD,DDỆợỦ‫̖ܤ‬ễွ˩8or4̿ችࡇᚘም ̊ᾉᾧᾑᾛᾐᾢễỄỉέᘍᄂᆮẝụ DDửỶὅἻỶὅ‫᧏ޒ‬ỂẨỦợạܱᘺẴủịẆἅἋἚỊ20̿ᆉࡇỉ‫ف‬ьỂOK, ẰỤỆᙲ൭Byte/flopỊ˯ẟỉỂ˄ьႎễᚘምểẲềɼỉ̿ችࡇᚘምử‫ן‬ᡐẲễẟểʖेẇ

2011/10/28

SSᄂ

4৆ဦႎễ‫ᧈ̿ٶ‬ἏἧἚૢͳỉ᣻ᙲࣱ ር‫׊‬Ị࠼ẟỉỂẆἑὊἄἕἚửᘍЗᚘምỆ዁Ếềẇ ᾑᾛᾐᾢỉdouble-double༿ỉ‫ܦ‬μܱᘺ ߻ૠỊᾎҥኝỆ10ʴஉỂỂẨỦẦᾎ இᢘ҄ỊỴὅἿὊἼὅἂᆉࡇ ‫ݣ‬ᝋἽὊἓὅầ‫ٶ‬ૠÆ‫ٳ‬දἫὊἋỉ˺ಅẇ ἅὊἙỵὅἂਃ࢘ᙐૠẇૠ͌ᚘምỉჷᜤầễẪềờӧᏡễር‫׊‬Ểẇ ᭗ࡇễἨἿἕἁ҄ỊᇹʚἋἘἕἩ

 DD-BLAS [2011] 2011࠰12உửႸ೅Ệπ᧏ửแͳ LAPACKễỄỉɥˮૠ͌ᚘምỴἽἆἼἌἲồỉ‫᧏ޒ‬Ị,ɶဋ൞ỉMPACKểᡲઃửẇ
2011/10/28

5SSᄂ

Double-double (Knuth,Dekker,Bailey,Hida) ွ˩ᾃ̿ችࡇᾉᴾ a=(a.hi/a.lo) ̿ችࡇ๫ምᾁ̾Ểನ঺Ẇ(̊ᾉDD.MUL)̿ችࡇỉьም (a+b → sum+err)
a t b ửˎ‫ܭ‬Ẳễẟ

a t b ửˎ‫ܭ‬
VXEURXWLQHTXLFNBWZRBVXP DEVXPHUU UHDO DEVXPHUUVXP DEHUU E VD 
HQG

VXEURXWLQHWZRBVXP DEVXPHUU UHDO DEVXPHUUYVXP DE
 Y VXPDHUU D VXPY  EY 
HQG

4̿ችࡇỉьም (c=a+b)ᵿᵛᵿᵌᶆᶇᵉᵿᵌᶊᶍᵊᴾᶀᵛᶀᵌᶆᶇᵉᶀᵌᶊᶍᵊᴾᶁᵛᶁᵌᶆᶇᵉᶁᵌᶊᶍᴾ ểẴỦᴾ
ᴾᴾᴾᴾᴾᴾᵆᵿᵆᵏᵇᵛᵿᵌᶆᶇᵊᴾᵿᵆᵐᵇᵛᵿᵌᶊᶍᵊᴾᶀᵆᵏᵇᵛᶀᵌᶆᶇᵊᴾᶀᵆᵐᵇᵛᶀᵌᶊᶍᵊᴾᶁᵆᵏᵇᵛᶁᵌᶆᶇᵊᴾᶁᵆᵐᵇᵛᶁᵌᶊᶍᵇᴾ

VXEURXWLQHDGG DEF UHDO D  E  F  VXPHUU
FDOOWZRBVXP D  E  VXPHUU HUU HUUD  E  FDOOTXLFNBWZRBVXP VXPHUUF  F  
HQG
2011/10/28

SSᄂ
S
Sᄂ

6

Double-double (Knuth,Dekker,Bailey,Hida)
 (̊ᾉDD.MUL, IntelỉSSEἾἊἋἑử̅ạ‫ئ‬ӳ)̿ችࡇỉʈም (a*b → prod+err)
VXEURXWLQHVSOLW DKLOR 
UHDO DKLOR
KL ,$1' D[IIIIIIIIIF 
OR DದKL
HQG

VXEURXWLQHWZRBSURG DESURGHUU UHDO DESURGHUUUHDO DBKLDBOREBKLEBORSURG D EFDOOVSOLW DDBKLDBOR FDOOVSOLW EEBKLEBOR 
HUU  DBKL EBKLದSURG DBKL EBOR DBOR EBKL DBOR EBOR 
HQG

4̿ችࡇỉʈም (c=a*b)
ᵿᵛᵿᵌᶆᶇᵉᵿᵌᶊᶍᵊᴾᶀᵛᶀᵌᶆᶇᵉᶀᵌᶊᶍᵊᴾᶁᵛᶁᵌᶆᶇᵉᶁᵌᶊᶍᴾ ểẴỦᴾ
ᴾᴾᴾᴾᴾᴾᵆᵿᵆᵏᵇᵛᵿᵌᶆᶇᵊᴾᵿᵆᵐᵇᵛᵿᵌᶊᶍᵊᴾᶀᵆᵏᵇᵛᶀᵌᶆᶇᵊᴾᶀᵆᵐᵇᵛᶀᵌᶊᶍᵊᴾᶁᵆᵏᵇᵛᶁᵌᶆᶇᵊᴾᶁᵆᵐᵇᵛᶁᵌᶊᶍᵇᴾ

VXEURXWLQHPXO DEF UHDO D  E  F  SURGHUU
FDOOWZRBSURG D  E  SURGHUU 
HUU HUUD  E  D  E  FDOOTXLFNBWZRBVXP SURGHUUF  F   
HQG
SSᄂ
S
Sᄂ

2011/10/28

7

4.6, DD-BLAS dd_axpyỉࣱᏡ ʻ஭ầችɟ஽ἓἷὊἝὅἂẲẺ‫ئ‬ӳỂ
Xeon X3330 (2.67MHz) single core / 10.67GFLOPS
900
800

700
600
500

1

400

2

300

3

200

4

Мဇ
ἅỴૠ

100
8
10
100
200
400
600
800
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
1
11000
1
12000
13000
14000
15000

0

2011/10/28SSᄂ8৆ဦႎễ‫ᧈ̿ٶ‬ἏἧἚૢͳỉ᣻ᙲࣱ৆ဦႎỆỊ
DD-BLASễỄỉؕஜἽὊἓὅ፭ửૢͳẰỤỆ᭗ࡇễૠ‫ܖ‬ἏἧἚỸỹỴ᧙ૠ፭LAPACKởFFTỉૢͳ
Doubl-DoubleửМဇẲẺἩἿἂἻἱὅἂỉɟᑍ҄ίHida-Bailey
ỉqdἻỶἨἻἼὸ

1.
2.
3.

 C++: qdἁἻἋửἫὊἋỆẲềẆợẟࣱᏡầ஖ࢳẰủỦ Fortran90: moduleửဇẟỦẮểờӧᏡẻầࣱᏡỊợỨẲẪễẟẇἅὅἣỶἻẦ᭗ࡇễἩἼἩἿἍἕἇẇᐯѣἓἷὊἝὅἂỉ২ᘐầ࣏ᙲẇ4.

ɳЗ҄ ᾜᾟᾘ‫ࣖݣ‬ỊܾତẇOpenMPỊreductionầ἟ἕἁẇ GPGPUồỉ‫᧏ޒ‬

 DoubleầἇἯὊἚẰủỦˌЭỊdouble-floatờẝẾẺẇ
 Double-doubleỊܾତỆܱᘺӧᏡẦếẆଏ‫܍‬ᄂᆮẝụ౟Ὁ᭗೛ẆɶဋẆɶ᣺Ụ5.

2011/10/28

ỴἩἼἃὊἉἹὅἸὊἈἾἫἽỂỉ࣏ᙲࣱỉᛐᜤSSᄂ

9৆ဦႎễ‫ᧈ̿ٶ‬ἏἧἚૢͳỉ᣻ᙲࣱᢿЎắểỆỂẝ ẺụᢿЎắểỆỂẝẾẺụẆC++ἸὊἈẝẺụίᚘምᅹ‫ܖ‬ỂễẟἈẝẺụίᚘምᅹ‫ܖ‬Ểễễẟ᧏ႆᎍኺဌὸỂỉ᧏ႆỊʻࢸẝụạỦẻỨạẇ

ᾗᾟᾒỉᙻໜầỪẦỦᾢᾢᄂỆợỦྙέẲẺૢͳ˳СầஓộủỦẇẐʮẑởἯἋἚἬἑồỉ‫᧏ޒ‬

2011/10/28

ᾠᾉFujitsuỊ‫ࣖݣ‬ẴỦỉẦᾎ

SSᄂ

10

ᾁᢿᾉ

݅ᘍЗ
З༿‫׍‬ஊ͌ἏἽἢ
EIGENἉἼὊἌ
㻞㻜㻝㻝㻛㻝㻜㻛㻞㻤㻌

㻿㻿ᄂ

㻝㻝㻌

EigenἉἼὊἌᾌЎ૝ɳЗᚘምೞӼẬ݅‫ݣ‬ᆅᘍЗ‫׍‬ஊ͌ἏἽἢὊἣἕἃὊἊ :: ᧏ႆἤἋἚἼὊ


ᾁ὿὿ᾅ


SuperComputing2006, GBP Finalists
 Quantum Physics Simulation codeồỉࣖဇ ‫ྶע‬ἉἱἷἾὊἑểSGI AltixӼẬ


ᾁ὿὿ᾇ




ᾣᾁᾚᢃဇ᧏‫ڼ‬ႺࢸẦỤẐʮẑӼẬ‫׍‬ஊ͌ἏἽἢὊểẲề᧏ႆ᧏‫ڼ‬ᾁ὿ᾀ὿

T2KỆấẟề256ἠὊἛ(4096ἅỴ)ѣ˺ᄩᛐ ˂ẆBX2, FX1, RICCễỄ‫׎‬ϋɼᙲễ‫ᙹٻ‬೉ἁἻἋἑỂỉѣ˺౨ᚰ GPGPUồỉ‫ࣖݣ‬ờ᧏‫(ڼ‬eigen_g)


ᾁ὿ᾀᾀ


T2KỆấẟề(8192ἅỴ)Ểỉѣ˺ᄩᛐฎỚ(10TFLOPSឬẆࢸᡓὸ K᧏ႆ᧙ᡲἂἽὊἩỆ੩̓ί̊ᾉRSDFTὸ ỼὊἩὅἏὊἋểẲềỉπ᧏ʖ‫(ܭ‬eigen_K)
 Open Petascale Library ồỉӋьʖ‫ܭ‬
 CRESTẐἯἋἚἬἑἋἃὊἽ᭗ࣱᏡᚘምỆ᝻ẴỦἉἋἘἲἏἧἚỸỹỴ২ᘐỉоЈẑỉờểỂẆἯἋἚἬἑἋἃὊἽᚘም࿢‫ؾ‬ӼẬ‫׍‬ஊ͌ἏἽἢὊ᧏
ႆ᧏‫(ڼ‬eigen-Exa)
2011/10/28

SSᄂ

12

Quick Review of the algorithms,
Diagonalization, for a dense real symmetric matrix
 Householder-D&C
Method:Householder
O(N3)Dense matrix

Householder:

Divide&ConqureO(N2)~O(N3)

Tridiagonal

Back Transform.
O(N3)Eigenpairs

Eigenpairs

Self-development

Divide & Conqure:Porting ScaLAPACK

Back Transform:

Self-development

2011/10/28

Block Algorithm
Loop unrolling2-D distribution

13

SSᄂ

ᾔᾢỂỉ‫׍‬ஊ͌ἏἽἢὊỉࣱᏡ,
㼀㼛㼠㼍㼘
㻾㼑㼐
㻱㼕㼓
㻮㼍㼏㼗㼠㼞㼍㼒㼛3hours Earth Simulator:
4,096 VPU’s

Memory:ᴾ
2~3 T Bytes

[sec]

㻝㻜㻜㻜㻜
㻥㻜㻜㻜
㻤㻜㻜㻜
㻣㻜㻜㻜
㻢㻜㻜㻜
㻡㻜㻜㻜
㻠㻜㻜㻜
㻟㻜㻜㻜
㻞㻜㻜㻜
㻝㻜㻜㻜
㻜

Matrix dim.:
280,000~375,000㻺㻩㻞㻤㻜㻷 㻺㻩㻟㻜㻜㻷 㻺㻩㻟㻡㻜㻷 㻺㻩㻟㻣㻡㻷[Dimension of Matrix]
Reported at SC|06, Tampa.

Wee confirmed stability of the solver up to 375,000 X 375,000
375
5,䯲ᴾ Accuracy is pretty excellent up to 300K X 300K (confirmed).
㻞㻜㻝㻝㻛㻝㻜㻛㻞㻤㻌

㻿㻿ᄂ㻌
ᄂ

㻝㻠㻌

Eigen-s/Eigen-sx ഏɭഏɭˊᚘምೞӼẬỆ᧏ႆɶỉ‫׍‬ஊ͌ἏἽἢὊ
ഏɭˊˊᚘምೞӼẬ
ˊᚘ
ˊᚘም
ẬỆ
Ệ᧏ႆɶ
᧏ႆɶỉ
᧏ႆ


Eigen_s (‫ݣ‬ᆅᘍЗӼẬἏἽἢὊ)ࢼஹ᧏ႆỉ‫׍‬ஊ͌ἏἽἢὊ(SGI Altix, ES)ửἰἽἓἅỴὉἰἽἓἩἿἍἕἇӼ
ẬỆોᑣ(T2KỖẦ)ẆHouseholderᾂ᣻‫ݣ‬ᚌ҄ὺЎлወ඙ඥὺᡞ‫٭‬੭

Eigen_sx࠘ᘍЗửኺဌẴỦỴἽἆἼἌἲử੔ဇẇ‫ᙹٻ‬೉բ᫆Ể᭗ᡮẇ Ⴘ೅ೞᏡMPI/OpenMPἡỶἨἼἕἛἴἙἽỂỉѣ˺μềỉἴἙἽỂỉ᭗ࣱᏡᢋ঺(ScaLAPACKợụỊ᭗ᡮ) ᭗ẟἋἃὊἻἥἼἘỵ 104ἅỴộỂỊỖỗἼἝỴỂẆ105ˌᨀἦὊἁỉ1%ửႸ೅ἓἷὊἝὅἂἣἻἳὊἑ
ScaLAPACKểỊီễụẆࣱᏡἣἻἳἑỊЎ૝૾ࡸể཯ᇌѣႎỆСࣂӧᏡ BLASờᆢಊႎỆМဇẇˬẲẆἋἾἕἛԠỎЈẲ૾ࡸἰἽἓἅỴӼẬỴἽἆἼἌἲỉ੔ဇ


Narrow-band-reductionỴἽἆἼἌἲ+࠘ᘍЗӼẬЎлወ඙ඥ

2011/10/28

15

ᄂ
SSᄂ

Strong scaling and weak scaling for eigen_s/T2K Tokyo
3HUI
3
3HUIRUPDQFHVFDODELOLW\RIRXUHLJHQVROYHULQVWURQJVFDOLQJUIIRUPDQFH VFDODELOLW\ RI RXU HLJHQVROYH
YH
HU LQ VWURQJV
VFDOLQJ

1 

1 
1 

OLQHDU
SHDN

Performance scalability of our eigensolver in weak-scaling












10000




ካዙኦ㟿

1000
GFLOPS*)/236



100

Tridiagonalization(ideal)Backtransform(ideal)101㻞㻜㻝㻝㻛㻝㻜㻛㻞㻤㻌

㻿㻿ᄂ㻌

2

43

84

165

32 6

䝜䞊䝗ᩘ䠄⾜ิ䛾ḟඖ䛿䝜䞊䝗ᩘ䠴9600䠅 㻝㻢㻌

64 7ࣱᏡஇᢘ҄ὉἓἷὊἝὅἂỆếẟềẐբ᫆ໜẑᴾ᫆ໜ ᡈᡈˊႎễἰἽἓἅỴἩἿἍἕἇỊᚘምᏡщỆ‫ݣ‬ẲềẆႎễἩᚘምᏡ Ệ‫ݣ‬ẲềἩἿἍἕἇὊἳἴἼ᧓ỉἢὅἛࠢếộụẆἙὊἑ̓ዅᏡщầ˯ẟ ᭗ᡮ҄ỉਦᤆᾉ1.
2.

ἅὅἣỶἻỉஇᢘ҄৖ѣỂỉἏὊἋἅὊἛἓἷὊἝὅἂ ἽὊἩನᡯởἙὊἑỉ‫୼٭‬ ἽὊἩỴὅἿὊἼὅἂễỄίᾗᾟᾒỉἘỿἋἚỆẝỦؕஜ২ᘐὸ

3.

᭗ᡮễܱᘺửỂẨỦẻẬМဇẴỦ ᾑᾛᾐᾢᴾ (ATLAS, GotoBLAS, MKL, CUBLAS, MAGMA) ૠ͌ᚘምἻỶἨἻἼ

4.

ɥᚡỉẐբ᫆ໜẑửᚐൿẴỦợạễỴἽἆἼἌἲỉોᑣ ཎỆἨἿἕἁ҄ỆợỦLevel3 BLASỉ෇ဇ5.

μẪီễỦᚘምҾྸỉ‫ݰ‬λ2011/10/28

SSᄂ

17

Eigen_sxỉஇ‫ࣱٻ‬Ꮱᴾ (ᙹ೉: 512node,N=300K,
HA8000 ἁἻἋἑἉἋἘἲ 512 ἠὊἛἇὊἥἋᴾ Feb.2011)
N=
300000 NM=
9392
NV=
9392
NUM.OF.PROCESS=
2048 (
32
NUM.OF.THREADS=
4calc (u,beta)
569.725378274918
mat-vec (Au)2382.50363969803
7555.07765028280
2update (A-uv-vu) 372.795704603195
48283.8181281065
calc v
327.551038742065
v=v-(UV+VU)u149.752300024033
UV post reduction 1.38964557647705
COMM_STATBCAST :: 110.417062759399
REDUCE :: 685.837212562561
REDIST :: 618.934033155441
GATHER :: 0.000000000000000E+000
TRD-BLK
300000 3831.22618889809
9396.46949175666
GFLOPS
TRD-BLK-INFO
30000048
before PDSTEDC 0.361311912536621
PDSTEDC 401.591748952866
COMM_STAT

64 )BCAST :: 0.125138998031616
REDUCE :: 0.000000000000000E+000
REDIST :: 0.000000000000000E+000
GATHER :: 0.000000000000000E+000
D&C 408.536870002747
ERRCODE=TRBAK= 1203.30433797836
COMM= 221.230289697647
44876.4275966327
GFLOPS
54519.0569849592
GFLOPS
56566.4621065163
GFLOPS
COMM_STATBCAST :: 88.3142938613892
REDUCE :: 126.187926292419
REDIST :: 0.000000000000000E+000
GATHER :: 0.000000000000000E+000
TRBAKWY 1203.30767083168
TRDBAK
300000 1203.34389090538
44874.9512761035
GFLOPS

0

Total 5443.10710787773
16534.6736099997
GFLOPS
max|w(i)-w(i).true|/|w|= 3.233048489790741E-005
දॖᾉЎлወ඙ඥỉἅἋἚử໯ᙻẲềẟỦỉỂẆ
ܱᨥỉFLOPS͌Ị୼Ệ᭗ẟ

T2K 512ἠὊἛ(8192cores),ἛྸᛯἦὊἁ75.4TFLOPS㻞㻜㻝㻝㻛㻝㻜㻛㻞㻤㻌

㻿㻿ᄂ㻌

㻝㻤㻌

Eigen-KửМဇʖ‫ܭ‬ỉἏἧἚ ᇹʚЎ᣼ ᾡᾢᾓᾕᾣᾉᴾ ਀‫ޛ‬έဃẆ‫ޥ‬ဋẰỮίி‫ٻ‬ὸẆ ʻဋέဃẆɤඑẰỮίி‫ٻ‬ὸ Ẹỉ˂ẆཋࣱᄂỉἸὊἈễỄ ᇹׄЎ᣼ ᾟᾗᾐᾢᾔ: ᱅ဋẰỮίྸᄂὸ МဇἸὊἈਘ‫ٻ‬ỉẺỜỆ ỺἽἱὊἚπ᧏༿ỉૢͳίᾗᾁᾂ࠰ࡇὸ ᧏ႆ࿢‫ؾ‬ỉૢͳ

ᴾ ᴾ ÅẮủỊ‫׍‬ஊ͌ἏἽἢẻẬỉᛅỂỊễẟỉỂẆẮẼỤờ৆ဦႎỆỴἳἼỽỉἛὅỾἻ૙੉ỉợạỆởỤễẟểẟẬễẟὲ
2011/10/28

SSᄂ

21

SSᄂ

22

ộểỜ
2011/10/28

ஜଐỉấᛅỉộểỜ ‫ᧈ̿ٶ‬ᚘምἩἿἂἻἲỉ᣻ᙲࣱ ἬἑἋἃὊἽ଺ˊỆễụẆἙὊἑᘙྵἾἫἽỂችࡇɧឱầ๼‫נ‬ႎỆẝụẆẸủầvisibleỆễẾềẨềẟỦẇ ৆ဦႎỆ‫ᧈ̿ٶ‬ἏἧἚỸỹỴửૢͳẲềẆɭမửἼὊἛὶᚘ

ምᅹ‫ܖ‬ἉἱἷἾὊἉἹὅỉឋႎ‫҄٭‬ửࣱᏡểểờỆẆችࡇờԃỜềᎋảợạẇᾚӼẬ‫׍‬ஊ͌ἏἽἢ


2011/10/28

ᾔᾸᾶᾴ᾽-Kᴾ PFLOPSầӧᏡễἏἽἢπ᧏ʖ‫ܭ‬
Мဇᎍਘ‫ٻ‬ỉẺỜỉԗᡀἏἧἚỸỹỴૢͳSSᄂ

23

4.3 高精度線形代数演算ライブラリ MPACK の紹介

高精度線形代数演算ライブラリ MPACK の紹介

中田真秀

線形代数は有史以来、人類の知的好奇心を満足させるだけでなく、具体的にも応用面で役立ってきた、本質的に重要な分野である。特に計算を行うということに限れば、コンピュータの発展はより大規模な線形代数の問題を、より高速に解くこととパラレルになっているといって良い。
今までのアプローチでは、大規模な線形代数の問題を、高速に解くことにほぼ限られていた。つまり正確に解くということに関してはほとんど光があたって来なかった。私は線形代数の問題をコンピュータ上で正確に解くというパラダイムを提案したい。
コンピュータでの演算には誤差が入るが、これまで演算誤差、結果の正確性についてはあまり注目されて来なかった。IEEE 754 倍精度は 10 進 16 桁で十分だったからだ。
しかし、そもそも倍精度では足りない問題は多くある。例としては、反復法での連立一次方程式を解くことや半正定値計画問題を解くことだろう。ただ、これは氷山の一角であると考えられる。研究者たちは実際に倍精度では足りない問題などは、うまくゆかない問題、運が悪い問題などと諦めることが多く表に出ることが少ないからだ。さらに CPU の処理速度が爆発的に伸び、超巨大な問題を長時間かけて解く現在、精度が不足することは容易に考えられる。
こ れ ら を 踏 ま え 、 高 精 度 線 形 代 数 演 算 パ ッ ケ ー ジ で あ る MPACK(http://mplapack.sourceforge.net/) を開発してきた。線形代数演算パッケージ
は、BLAS, LAPACK が有名であるが、これをベースにし、高精度な演算ができるようにした。以下に特徴を挙げる。コンピュータ資源の許す限り高精度な計算や、四倍精度程度だが高速な計算も可能で、様々なニーズに応えられるようになっている。BLAS, LAPACK のように、線形代数演算のビルディングブロッ
ク、参照実装、プログラムインターフェースを提供する。C++で書かれており、5 つの高精度演算のための型をサポートしている (GMP, MPFR, QD, DD, 四倍精度)、さまざまなプラットフォーム (Linux, FreeBSD, Windows, Intel Xeon
Phi, NVIDIA GPU) をサポートし、オープンソースのライセンス(2-clause BSD
ライセンス)で配布している。

4.4 宇宙プラズマでのペタスケール数値計算（惑星磁気圏シミュレーション）

宇宙プラズマでのペタスケール数値計算（惑星磁気圏シミュレーション）
九州大学情報基盤研究開発センター 深沢 圭一郎
宇宙プラズマ研究において、我々は主に太陽から吹いてくる磁場を伴ったプラズマの風（太陽風）と地球の磁場が相互作用して起こる様々な現象を研究ターゲットにしている。
これらは宇宙空間で起きる現象であるため探査機を打ち上げて観測を行うが、基本的に”その場”の観測しか行えない。そのため、宇宙プラズマ計算機シミュレーションがこの分野の理論の発展、また観測結果の理解の促進に非常に重要な役割を果たしてきている。宇宙プ
ラズマの振る舞いはブラソフ方程式によって正確に記述されるが、特に我々が注目している太陽風と惑星の磁場の相互作用によって形成される磁気圏というグローバルな領域では、ブラソフ方程式を近似した電磁流体（MHD）近似というものがよく成り立ち、MHD 方程式を用いて惑星時磁気圏のシミュレーションをおこなっている。
近年計算機の発達により高精度なシミュレーションが行えるようになってきた。例として、土星磁気圏では今までのシミュレーション結果からは見ることができなかったスモー
ルスケールのプラズマ対流の渦構造が磁気圏内で見えてきている。この渦構造は電流を生成するため、惑星のオーロラ発光と関連づけられることがある。今回の結果は現在土星を観測しているカッシーニ探査機により撮像された土星オーロラとの相関性が見え、土星オーロラの発光に渦構造が関係しているということが強く示唆された。
一方でこのような高精度の計算も、実際に我々が見たい物理現象のスケールからはまだまだ粗く、今後の計算機の発展が期待されている。その際は超大規模並列計算を行うことになるが、そのような状況で効率良く計算できるか把握しておくことは重要である。そこで九州大学の計算機システムの一つである、富士通 PRIMERGY RX200 S6 を用いて MHD
シミュレーションコードのベンチマークを行った。5000 並列近くまで利用でき、weakscaling で性能を測定した。その結果実行効率は 30％を超え、非常に良い性能であったが、スケーラビリティが悪く、1000 並列時に比べ、5000 並列近くでは～10%の性能劣化が見えた。今後より多くの CPU を利用し、計算が行われることを考慮すると、この並列性能劣化は非常に大きな問題である。

Ᏹᐂ䝥䝷䝈䝬䛷䛾䝨䝍䝇䜿䞊䝹ᩘ್ィ⟬䠄ᝨᫍ☢Ẽᅪ䝅䝭䝳䝺䞊䝅䝵䞁䠅῝ἑᆂ୍㑻1, 2, 3
1. ஑ᕞ኱Ꮫ᝟ሗᇶ┙◊✲㛤Ⓨ䝉䞁䝍䞊
2. ஑ᕞ኱Ꮫᐂ✵⎔ቃ◊✲䝉䞁䝍䞊
3. CREST, JST

March 1, 2012

Context
1. Introduction
9 ኴ㝧ᆅ⌫ᝨᫍ⣔⛉Ꮫ䛾⤂௓
9 ᝨᫍ䛸≉ᚩ䛸☢Ẽᅪ䛾ᵓ㐀
9 MHD᪉⛬ᘧ䛸Vlasov᪉⛬ᘧ

2. ᅵᫍ☢Ẽᅪ䛾䝅䝭䝳䝺䞊䝅䝵䞁
9 䝥䝷䝈䝬 ᵓ㐀䛸䜸䞊䝻䝷䛾㛵㐃

3. MHD䝁䞊䝗䛾ᛶ⬟ ᐃ
9 ᵝ䚻䛺ィ⟬ᶵ䝅䝇䝔䝮䛷䛾ᛶ⬟

4. 䜎䛸䜑

1

Introduction -12ኴ㝧ᆅ⌫ᝨᫍ⣔⛉Ꮫ䛸䛿

Introduction -2

3

Introduction -3

4

Introduction -45ᆅᆅ⌫䚸ᮌᫍ䚸ᅵᫍ䛻䛚䛡䜛㟁☢Ẽⓗ≉ᚩJupiter

SaturnEarth

Magnetic field [nT]420,000

21,00031,000Magnetic polarity

N pole is north N pole is northN pole is south

Rotation period [hr]

10

10.65

24

Main plasma source

Io, ionosphere

Enceladus,
ionosphereionosphere

Equatorial Radius [km]71,49260,2686378

From Sun [A.U.]

5.2

9.55

1

ᮌᫍ䛿ᕧ኱䛺☢ሙ䛸㇏ᐩ䛺䝥䝷䝈䝬䜢ᣢ䛳䛶㧗㏿⮬㌿䛩䜛.
ᅵᫍ䛿㇏ᐩ䛺䝥䝷䝈䝬䜢ᣢ䛳䛶㧗㏿⮬㌿䛩䜛䚹Introduction -56ᆅᆅ⌫䚸ᅵᫍ☢Ẽᅪ䛾ᙧ

Enceladus
Plasma disk
Fig. 1. A schematic of Terrestrial magnetosphere

Fig. 2. Schematic of Saturn’s magnetosphere[Kivelson, 2006]

ᅵᫍ䛿㧗㏿⮬㌿䛧䛶䛔䜛䛯䜑䚸☢Ẽᅪෆ䛾䝥䝷䝈䝬ᑐὶ䛜」㞧䚹

Introduction -6
Spacecrafts

 8ᶵ䛾᥈ᰝᶵ䛜ᮌᫍ䜢ほ Pioneer10 (1973), Pioneer11 (1974), Voyager1 (1979),
Voyager 2 (1979), Ulysses (1992), Galileo (1995-2003),
Cassini (2000), and New Horizon (2007) 4ᶵ䛾᥈ᰝᶵ䛜ᅵᫍ䜢ほ Pioneer11 (1979), Voyager1 (1980),
Voyager 2 (1981), Cassini (2004 - now), Future missions
 JUNO (launch at 2011, arrive at Jupiter in 2016)
 EJSM (Jupiter mission: launch after 2020)
 TANDEM (Saturn mission: launch after 2020)7

Motivation
ᕧᕧ኱䛷」㞧䛺ᝨᫍ☢Ẽᅪ䜢ィ⟬ᶵ䝅䝭䝳䝺䞊䝅䝵䞁䛷ㄪ䜉䛯䛔䛺䛬ᩘ್䝅䝭䝳䝺䞊䝅䝵䞁䛛 ᝨᫍ䛾100ಸ௨ୖ䛻ᗈ䛜䜛☢Ẽᅪ඲య䜢ぢ䜛䛻䛿ほ䛷䛿㞴䛧䛟䚸䝅䝭䝳䝺䞊䝅䝵䞁䛧䛛↓䛔䚹 ⾨ᫍほ 䛷䛿䛒䜛᫬㛫䛾✵㛫୍Ⅼ䛧䛛䜟䛛䜙䛪䚸᫬✵㛫䛾ศ㞳䜒㞴䛧䛟䚸≀⌮⌧㇟䛾⌮ゎ䛜㞴䛧䛔䚹 ☢Ẽᅪ䛿10m䡚109m䝇䜿䞊䝹䜎䛷䛾⌧㇟䜢ྵ䜐䝬䝹䝏䝇䜿䞊䝹⎔ቃ䛷䛒䜚䚸䝅䝭䝳䝺䞊䝅䝵䞁䛻᭱㐺䚹 ᭱⤊ⓗ䛻䛿Ᏹᐂ⎔ቃ䜢⌮ゎ䛧䚸ணሗ䛩䜛䛣䛸䛷ே㢮䛜Ᏹᐂ✵㛫䛻ฟ䛶⾜䛟㝿䛻Ᏻ඲䛺᝟ሗ䜢ᥦ౪䛧䛯䛔䚹㟁☢ὶయ䝁䞊䝗 -1ᏱᏱᐂ䝥䝷䝈䝬䜢ྲྀ䜚ᢅ䛖᪉⛬ᘧ(1)Vlasov᪉⛬ᘧ ↓⾪✺Boltzmann᪉⛬ᘧ䛸Maxwell᪉⛬ᘧ䛛䜙ᡂ䜛䝥䝷䝈䝬䛾᣺䜛⯙䛔䜢᭱䜒ṇ☜䛻⾲⌧䛷䛝䜛᪉⛬ᘧ⣔㏿ᗘศᕸ㛵ᩘf(x, v, t)䜢⪃䛘䜛䛸䚸
wf
wf
wf q
 v   (E  v u B ) 
wt
wx m
wv

0

䛣䜜䛸Maxwell᪉⛬ᘧ䜢㐃❧䛧䛶ゎ䛟䚹 䛧䛛䛧䚸఩⨨(x, y, z)3ḟඖ䚸㏿ᗘ(vx, vy, vz)3ḟඖ䛸᫬㛫䛛䜙ᡂ䜛㠀⥺ᙧ᪉⛬ᘧ⣔䛷䚸ゎ䛟䛣䛸䛜ᅔ㞴䚹ព࿡䛾䛒䜛f(x, y, z, vx, vy, vz, t)䜢ィ⟬䛩䜛䛻䛿8PB⛬ᗘ䛾䝯䝰䝸䛜ᚲせ䛷䛒䜚䚸⌧ᅾ䛾ィ⟬ᶵ䛷䛿䝯䝰䝸୙㊊䚹㟁☢ὶయ䝁䞊䝗 -2ᏱᏱᐂ䝥䝷䝈䝬䜢ྲྀ䜚ᢅ䛖᪉⛬ᘧ(2)MHD (Magnetohydrodynamics)᪉⛬ᘧ Vlasov᪉⛬ᘧ䛾nḟ䝰䞊䝯䞁䝖ྲྀ䜛䛣䛸䛷䚸ồ䜑䜙䜜䜛䚹0ḟ(㏿ᗘ✵㛫䛷✚ศ)䚸1ḟ(v䛛䛡䛶✚ศ)䚸2ḟ䠄v2䛛䛡䛶✚ศ䠅䜘䜚䚸wU  ( vU )
wt
wv
1
1
 ( v   ) v  p  J u B
wt
U
U
wp
( v  ) p  Jp  v
wt䜢ᚓ䜛䚹䛣䜜䜙䛸☢ሙ䛾ㄏᑟ᪉⛬ᘧwBwt

 u ( v u B)䜢䜎䛸䜑䛶MHD᪉⛬ᘧ䛸䛔䛖䚹䛣䛾ሙྜḟඖ䛜f(x, y, z, 8)䛻䛺䜚䚸౑⏝䝯䝰䝸䛜⃭ῶ䛩䜛䚹ᅵᫍ☢Ẽᅪ䝅䝭䝳䝺䞊䝅䝵䞁111ᅵᅵᫍ☢Ẽᅪ኱つᶍ䝅䝭䝳䝺䞊䝅䝵䞁 ᵓ㐀䛸䜸䞊䝻䝷䛾㛵ಀ䜢ㄪ䜉䛶䛔䜛䚹ᕥഃ䛻 ᵓ㐀䛜ぢ䛘䛯ᚋ䚸ྑഃ䛻䜒 ᵓ㐀䛜ぢ䛘䜛䚹
ほ 䛷䜒ྠᵝ䛾ᵓ㐀䛜ሗ࿌䛥䜜ጞ䜑䛶䛔䜛䚹

Fig. 3. ᅵᫍ☢Ẽᅪ㉥㐨㠃䛻䛚䛡䜛☢ሙᙉᗘ[Fukazawa et al.,
JGR, 2012]ᅵᫍ䛻䜒䜸䞊䝻䝷䛜䛒䜛12ᅵᫍ☢Ẽᅪ䝅䝭䝳䝺䞊䝅䝵䞁213

 ἢ☢ຊ⥺㟁ὶ䛸ᵓ㐀䛾㛵ಀἢ☢ຊ⥺㟁ὶ䛾ᙉ䛔⟠ᡤ䛛䜙☢ຊ⥺䜢ఙ䜀䛩䛸䠊䠊䠊኱つᶍィ⟬䛻䜘䜛㧗⢭ᗘ䛺䝅䝭䝳䝺䞊䝅䝵䞁䛻䜘䜚䚸ᅵᫍ ᵓ㐀䛸ᩬⅬ䜸䞊䝻䝷ᵓ㐀䛾㛵㐃ᛶ䜢ึ䜑䛶♧၀䚹ほ ⤖ᯝ䛻ఝ䛯ᵓ㐀䜢෌⌧䚹

Fig. 4. ᅵᫍ☢Ẽᅪ䛻䛚䛡䜛[Fukazawa et al., JGR, 2012]ᵓ㐀䛸ἢ☢ຊ⥺㟁ὶ䛾㛵ಀ

Fig. 5. ᅵᫍ༡ᴟ䛻䛚䛡䜛䜸䞊
䝻䝷䛾㍤ᗘ[Grodent et al., 2011]䛥䜙䛻㧗⢭ᗘ䛾ィ⟬䜈䜸䜸䞊䝻䝷䛾䜘䛖䛺⌧㇟䜢⢭ᗘⰋ䛟ゎ䛟ᚲせ䛜䛒䜛 ௒䜎䛷䛾ィ⟬䛷䛿0.1 RS䠄RS䛿ᅵ༙ᫍᚄ = 60,300 km䠅䛾᱁Ꮚᖜ䜢฼⏝䚹 䝬䜽䝻䛸䝭䜽䝻䛾㑄⛣㡿ᇦ䠄MHD㏆ఝ䛾㝈⏺㡿ᇦ䠅䜢ィ⟬䛩䜛䛯䜑䛻䛿䚸0.01 RS䛾᱁Ꮚᖜ15,0003⛬ᗘ䛾䜾䝸䝑䝗䛜ᚲせ䚹 CPU䛾࿘Ἴᩘ䛜䛭䜜䜋䛹ୖ䛜䜙䛺䛔䛸௬ᐃ䛩䜛䛸
15,0003䜾䝸䝑䝗䜢⌧ᐇⓗ䛺᫬㛫䛷ィ⟬䛩䜛䛯䜑䛻䛿200୓䝁䜰䛜ᚲせ䚹 ୍᪉䛷䚸HA8000@ᮾ኱䜢฼⏝䛧䛯8192୪ิ䛾ィ⟬䛷䛿䚸୪ิ໬ຠ⋡䛜1024୪ิ䛸ẚ䜉⣙5䠂పୗ䛧䛶䛔䜛䚹㉸኱つᶍ୪ิ䛿኱୔ኵ䛛?

15

MHD䝅䝭䝳䝺䞊䝅䝵䞁䛻䜘䜛ᛶ⬟ホ౯MHD䝅䝭䝳䝺䞊䝅䝵䞁䛾ᛶ⬟ホ౯᮲௳16䝧䝧䞁䝏䝬䞊䜽䝉䝑䝔䜱䞁䜾ゎἲ䛿Leap-frog+Lax Wendroffἲ䛾ΰྜᡭἲ ☢ሙ䜢ᢅ䛖᪉⛬ᘧ䛷䛿div䞉B=0䜢ಖド䛩䜛ᡭἲ䜢⏝䛔䜙↓䛡䜜䜀ᡂ䜙䛪䚸㐺⏝ྍ⬟ᡭἲ䛻㝈⏺䛜䛒䜛䚹 ௒䜎䛷䛾䝧䞁䝏䝬䞊䜽䛸䛾ẚ㍑䚸䜎䛯䚸simple䛺ᡭἲ䜢⏝䛔䜛䛣䛸䛷䚸ィ⟬ᶵ䛾ᛶ⬟䛜ฟ䜔䛩䛔䛸⪃䛘䜙䜜䜛䚹ᇶᮏ䝃䜲䝈

 1core䛒䛯䜚f(x, y, z, 8) = 691200grid 䜢タᐃ䚹 ィ⟬ୖ䛣䛾7ಸ䛾㓄ิᩘ䜢౑⏝䛩䜛䚹

MPI୪ิ 䛾䜚䛧䜝䛿๓ᚋ୍䛴䚹OpenMP䛿౑⏝䛫䛪䚹䜻䝱䝑䝅䝳ホ౯ 䝇䜹䝷ᶵ䛷䛿䜻䝱䝑䝅䝳䝠䝑䝖ྥୖ䜢ぢ㎸䜣䛷䚸㓄ิ䛾୪䜃᭰䛘䜢⾜䛖䛸ᛶ⬟䛜ୖ䛜䜛ᶵ✀䛜䛒䜛䠄㠀x86⣔䠅䛾䛷䚸䛭䛾ຠᯝ䜒ホ౯䛩䜛 Type A䠖f(i, j, k, m)䚸Type B䠖f(m, i, j, k)

MHD䝁䞊䝗䛻䛚䛡䜛ᛶ⬟ホ౯ᡭἲ173✀㢮䛾㡿ᇦศ๭ἲ1ḟඖィ⟬᫬㛫䠄TS䠅TS 1

k 2 n 2 ( p  1)㏻ಙ᫬㛫䠄TC䠅

TC1

k 2 n ( p  1)
2

2ḟඖ

TS 2

k 2 n 2 ( p  1)

TC 2

2k 2 n ( p  1)
2

1
2

3ḟඖ

TS 3

k 2 n 2 ( p  1)

TC 3

3k 2 n ( p  1)2

1
3

Estimated Time18

PRIMERGY RX200 S6

19

2011ᖺ5᭶䛻✌ാ䛧䛯஑ᕞ኱Ꮫ䛾ィ⟬ᶵ䝅䝇䝔䝮
PRIMERGY RX200 S6CPUᆺᘧ

6 core Xeon X5670

࿘Ἴᩘ

2.93 GHz (70.32GFlops)

Cache

L1: 64KB/core
L2: 256MB/core
L3: 12MB/CPU

Memory Bandᖜ
B/F್

32/70.32

0.45

Node

Coreᩘ

12

䝯䝰䝸

48GB

Nodeᩘ

392 (4704core)⌮ㄽᛶ⬟

55.13TFlops

ᐩኈ㏻〇PC䜽䝷䝇䝍ᆺィ⟬ᶵWestmere䝁䜰

32GB/s /CPU

SystemNode㛫㏻ಙ InfiniBand 4x QDR (4GB/s)MHD䝁䞊䝗䛾ᛶ⬟ホ౯(2)20sugoka䠄PRIMERGY RX200 S6䠅䛷䛾

Scalability

Performance
20,000

GFlops

0.96
0.95

0.94
0.92
0.91

3DdecompositionTypeB

16,0000.99

1 0.97

17,627

3DdecompositionTypeA

ᐃ

0.93
0.9 0.92

12,0000.940.94

0.91

9,023

0.89

12,031

0.87

8,000
0.8

4,737
3,506

4,000

3DdecompositionTypeA

6,113

3DdecompositionTypeB

1,495
767

3,158
2,265
0
1,009
0 501
1,000

0.7
2,000

3,000

4,000

5,000

0

1,000

Number of cores

2,000

3,000

4,000

5,000

Number of cores• ᐇຠᛶ⬟䛿TypeA䛷4700䝁䜰᫬䛻᭱኱17TFlops䠄ᐇ⾜ຠ⋡33%䠅ฟ䛶䛔䜛䚹• 䝇䜿䞊䝷䝡䝸䝔䜱䛿ప୪ิ᫬䛻୙Ᏻᐃ䛷䚸⦆䜔䛛䛻ୗ䛜䛳䛶䛚䜚䚸TypeB䛾䜋䛖䛜10%㏆䛟㧗䛔䚹

MHD䝁䞊䝗䛾ᛶ⬟ホ౯(3)21⌧≧䛾ྛ✀䝇䝟䝁䞁䛸䛾ẚ㍑
Coreᩘ
ᩘ
/CPUᩘ

ᐇຠᛶ⬟
[GFlops]⌮ㄽᛶ⬟[GFlops]

ᐇ⾜ຠ⋡
[%]

㡿ᇦศ๭

CPU✀㢮

SX-9

64/64

2188

6553

33

2ḟඖ

䝧䜽䝖䝹

SX-8R

8/8

80

282

28

1ḟඖ

䝧䜽䝖䝹

HA8000 8192/1024

10038

75366

13

3ḟඖA

Opteron

HX600

1024/256

2166

10240

21

3ḟඖA

Opteron

FX1

1024/256

2081

10240

21

3ḟඖB

SPARC64VII

SR16000 1344/672

5375

25267

21

3ḟඖB

POWER6PG S64704/144

18267

55130

33

3ḟඖA

Xeon

PG S3

1536/768

2536

18432

14

3ḟඖA

Xeon

MHD䝁䞊䝗䛾ᛶ⬟ホ౯(4)22SX-9䜢1䛸䛧䛯䛸䛝䛾ྛ✀䝇䝟䝁䞁CPU䚸core䛾ᛶ⬟
@core

@CPU

✀㢮
CPU✀

SX-9

1.00001.0000

䝧䜽䝖䝹

SX-8R

0.29250.2925䝧䜽䝖䝹

HA80000.03580.1432

Opteron

HX600

0.06180.2474

Opteron

FX1

0.05940.2377SPARC64VII

SR16000

0.11700.2340POWER6PG S60.11870.7122

Xeon

PG S3

0.0483

0.0966Xeon

ᐇ䜰䝥䝸䛷䛿FX1䚸
SR16000䜔HX600䛜4CPU
䛷SX-9䛾1CPU䛸ྠ➼䛾ᛶ⬟䛻䛺䛳䛶䛔䜛䚹
PRIMEGY RX200 S6䛷䛿
1CPU䛷SX-9䛾CPU䛾7๭
㏆䛟䛾ᛶ⬟

Summary23Ᏹᐂ䝥䝷䝈䝬䛷䛾䝨䝍䝇䜿䞊䝹ᩘ್ィ⟬ ᝨᫍ☢Ẽᅪ䝅䝭䝳䝺䞊䝅䝵䞁 ほ 䛷䛿᝟ሗ䛜㊊䜚䛺䛔䛯䜑䚸ᩘ್䝅䝭䝳䝺䞊䝅䝵䞁䛜ຊ䜢Ⓨ᥹䛷䛝䜛ศ㔝䚹 ᅵᫍ☢Ẽᅪ䛷䛿 ᵓ㐀䛸䜸䞊䝻䝷䛾㛵ಀ䛜ほ 䛛䜙♧၀䛥䜜䛶䛔䜛䛜䚸䝅䝭䝳䝺䞊䝅䝵䞁䛷䜒☜ㄆฟ᮶䛶䛔䜛䚹 䜎䛰䜎䛰ィ⟬㈨※䛜ᚲせ䛷⌮ㄽⓗ䛻䛿䜶䜽䝃䜎䛷ᑐᛂྍ⬟䚹 MHD䝁䞊䝗䛾ᛶ⬟ᐃ sugoka䛷䛿㧗䛔ᐇ⾜ຠ⋡䛜ฟ䛶䛔䜛୍᪉䛷䚸䝇䜿䞊䝷䝡䝸䝔䜱䛾పୗ䛜ၥ㢟䚹 ィ⟬ᶵ䝅䝇䝔䝮䛤䛸䛻᭱㐺䛺୪ิ໬ᡭἲ䜒㐪䛖䛯䜑䚸䝷䜲䝤䝷䝸䛷ᑐᛂྍ⬟䛰䛸䛖䜜䛧䛔䚹

4.5 ペタスケールのアプリケーションへ向けて

ペタスケールのアプリケーションへ向けて：日本原子力研究開発機構・システム計算科学センター・シミュレーション技術開発室の取り組み
日本原子力研究開発機構
システム計算科学センター
シミュレーション技術開発室
山田進、佐々成正、中村博樹、奥村雅彦、町田昌彦
理化学研究所今村俊幸
日本原子力機構・システム計算科学センター・シミュレーション技術開発室では、現在、原子力研究開発における喫緊の課題に対し、最先端の大規模シミュレーションを通して、
その課題解決を図るため、ペタスケール規模のシミュレーション技術の研究開発を理化学研究所等と連携し進めている。本公開資料では、まず、原子力機構全体の主たる大規模シ
ミュレーションによる取組みを紹介し、原子力研究開発における大規模シミュレーションの役割とその成果について概観した後、当該研究室が取り組んでいるペタスケールライブ
ラリ開発について、その現状と課題について述べる。また、現在、原子力分野において極めて大きな課題となっている福島原発事故後の環境修復の研究課題については、計算科学による研究の現状を紹介することとした。広範囲に拡散した放射性核種による汚染からの環境修復のための研究は、原子力分野がこれまでに経験してきた問題のレベルを遥かに超え、国家的課題となっており、分野を超えて日本（世界も）の英知を結集して取り組むべ
き課題であることから、本紹介を通して様々な研究者のご意見を伺い、また、可能であれば研究連携を実施して行きたいと考えている。
ペタスケールライブラリ開発において、当該研究室が特に力を入れているのは、行列の対角化に係るソルバー開発である。当該研究室では、原子炉構造材料、原子炉核燃料、そして原子力研究開発に関係する機能材料の第一原理計算を様々なレベル（超並列大規模シ
ミュレーションによる核燃料の高温挙動や、構造材料スクリーニングのためのパラメータ
サーベイ・シミュレーション）で進めている他、次世代デバイス材料と考えられている強相関電子系の電子構造の高精度計算を実施するため、大規模行列の高速対角化に焦点を当て研究開発を進めてきた。その結果として、当該研究室では Eigen-K と名付けた「京」上
にて高速に動作する行列対角化ソルバーの開発に成功し、研究室 HP より現在、公開している。尚、対象とする行列は密行列であり、1 万コア程度までは、ほぼリニアにスケールする並列性能を示し、1 万要素程度以上の行列になると、ScaLAPACK より遥かに十分に高い性能を示すことが分かっている。また、エルミート行列用のライブラリーも準備されており、実対称及びエルミートの大規模密行列を高速に対角化するコードが公開されている。
福島原発事故後、当該研究室では放射性物質セシウムの土壌表面での吸着機構を明らかにするための計算科学研究課題に取り組んできた。土壌には、粘土粒子と称される直径２μｍ以下の細粒成分が多く含まれているが、その原子分子レベルでの吸着化学機構を明らかにするため、第一原理計算手法や分子動力学手法を駆使して研究を行っている。本発表
ではその一端を紹介し、物質科学における計算科学研究がどれほど、現在立ちはだかる大きな課題（吸着土壌から脱離させ除染による廃棄土壌を減容化する）の解決に貢献できるかを議論した。物質科学から見ると、土壌のような極めて複雑な系に対し、効果的な議論を進めるためには、できる限り現象の普遍性を捉えるモデル化を進めて、正確な計算手法
を用いて解くことと、大規模計算を実施し、非一様性の中での普遍性を逆に確認する必要がある。超並列計算が今後は重要な役割を果たすはずとして、今後の展開をできれば先導
したいものと考えている。

䈽䉺䉴䉬䊷䊦䈱䉝䊒䊥䉬䊷䉲䊢䊮䈻ะ䈔䈩㩷䋺䉲䊚䊠䊧䊷䉲䊢䊮ᛛⴚ㐿⊒ቶ䈱ข䉍⚵䉂㩷ᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴䉲䉴䊁䊛⸘▚⑼ቇ䉶䊮䉺䊷䉲䊚䊠䊧䊷䉲䊢䊮ᛛⴚ㐿⊒ቶ↸↰᣽ᒾਥ䈢䉎⎇ⓥදജ⠪䊤䉟䊑䊤䊥䋭㐿⊒䋺ጊ↰ㅴ䇮૒䇱ᚑᱜ䇮੹᧛ବᐘ䋨㔚ㅢᄢ䇮⃻䊶ℂ⎇䋩䉲䊚䊠䊧䊷䉲䊢䊮⎇ⓥ䋺ᅏ᧛㓷ᒾ䇮ਛ᧛ඳ᮸ን჻ㅢSS⎇䇮2012ᐕ3᦬1ᣣ
1⊒⴫ౝኈ䋱䇮ේሶജᯏ᭴䈱ਥ䈢䉎䉝䊒䊥䉬䊷䉲䊢䊮
2䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂
3䇮⑔ፉ੐᡿ᓟ䈱⸘▚⑼ቇ䈱ข䉍⚵䉂䋱䇮ේሶജᯏ᭴䈱ਥ䈢䉎䉝䊒䊥䉬䊷䉲䊢䊮㩷㩷ේሶജ⎇ⓥ㐿⊒䋴ᄢ䊒䊨䉳䉢䉪䊃

3䋱䇮ේሶജᯏ᭴䈱ਥ䈢䉎䉝䊒䊥䉬䊷䉲䊢䊮㩷㩷ᾀᴾఋᗡӳἩἻἌἰʏ්ἉἱἷἾὊἉἹὅᴾ

ᾁᴾҾ‫܇‬щೞ֥ϋඖᬸʚႻ්ᚐௌᴾ

ᾂᴾἜἚἼỸἲϬҲ᭗ᡮ໐ỆấẬỦ᩼‫්ࠝܭ‬ѣྵᝋᴾ

ᾃᴾᇹɟҾྸᚘምỆợỦ஬૰э҄ἳỽἝἌἲᴾ ίቩမᏤ҄ὸỉᄂᆮᴾ

ᾄᴾẸỉ˂ᚘምᅹ‫ܖ‬ᄂᆮʙ̊ᴾ

4

ᾀᴾ

ᩭⲢว䊒䊤䉵䊙ੂᵹ䉲䊚䊠䊧䊷䉲䊢䊮㩷ᚘምἅὊἛӸᾉᵥᵲᵓᵢᴾ
Ӹᾉᵥᵲ
ᵲᵓᵢ
ᵢ

‫׎‬ᨥ༏ఋᗡӳ᬴ܱ໐ᴾᵧᵲᵣᵰᴾᵓഏΨᢃѣᛯἴἙἽỆợỦఋᗡӳ໐࣎μ˳ỉᴾἩἻἌἰʏ්᠞ᡛἉἱἷἾὊἉἹὅᴾᐮမἩἻἌἰᚾ᬴ᘺፗᴾᵨᵲᵔᵎᵳᴾ

‫׹ݱ‬ᘺፗἴἙἽᴾ
ᵆྵཞἇỶἌᵇᴾ‫׹ݱ‬ᘺፗᴾ

ᵨᵲᵔᵎᵳᴾ

ᵧᵲᵣᵰᴾ

໐࣎˳ᆢ‒

6m3

60m3

800m3

఍‫܇‬ૠ‒

5.24 x 109

1.38 x 1011

1.23 x 1012

∇∈∐᣽‒

800GB

14TB

147TB

ᚘም᣽‚ܱјࣱᏡ∝ଐ‛‒

30TFlopsὉday

750TFlopsὉday

7PFlopsὉday
5

ᾃᴾ╙৻ේℂ⸘▚䈮䉋䉎㩷╙㩷᧚ᢱഠൻ䊜䉦䊆䉵䊛䋨☸⇇⣀ൻ䋩䈱⎇ⓥ㩷

ᚘምἅὊἛӸᾉӸᾉᵴᵟᵱᵮ
ᵮᴾ

э҄ἳỽἝἌἲử৭੮ẲễẬủịẆᅹ‫ܖ‬ႎὉӳྸႎఌਗỆؕỀẪэ҄ʖยầỂẨễẟᴾ‫ٳ‬਴ỆợỦᴾ‫ٳ‬਴ỆợỦ
சஹỉэ҄ʖยᴾ

Ҿ‫܇‬໐ನᡯ஬૰ỉэ҄ἳỽἝἌἲᚐଢửᡫẲềẆҾ‫܇‬໐ỉ᭗ኺ࠰҄‫ݣ‬ሊỆᅹ‫ܖ‬ႎὉӳྸႎఌਗửᴾ
੩̓ẴỦẮểẇᴾ
ஜᄂᆮỂỊẐቩမᏤ҄ẑểẟạэ҄ỆếẟềẆҾ‫܇‬Ὁᩓ‫܇‬ἾἫἽẦỤẆẸỉἳỽἝἌἲửଢỤẦỆẴỦẇᴾ

஬૰э҄ᴾᚘምỉႸႎẆಒᙲᴾ
ወᚘϼྸᴾ࣯ນỆэ҄ẴỦᴾӧᏡࣱỊᾎᴾэ҄ἙὊἑᴾ
ܱย͌ᴾ˴࠰έộỂᴾ‫ٳ‬਴ỂẨỦẦᾎᴾ
ᢅӊᴾ

ྵ‫נ‬ᴾ

சஹᴾ

଺᧓ᴾ

ᚘᴾ ምᴾ ϋᴾ ܾᴾኽ୒ቩမỉҾ‫܇‬ἴἙἽử˺ụẆᇹɟҾྸᚘምỆợẾềɧኝཋΨእỉቩမồỉᨼộụởẴẰẆቩမࢍࡇỉ‫҄٭‬ửᛦỔỦẇᴾ
6䋭ਥⷐ䈭䉮䊷䊄䈪䈱䊤䉟䊑䊤䊥䊷㩷㩷ἅὊἛӸᴾἅὊἛಒᙲᴾɼᙲἏἽἢὊᴾ

ᵥᵲᵓᵢᴾ

ἧỻἕỽὊὉἩἻὅἁ૾ᆉࡸὺἯỴἏὅ૾ᆉࡸᴾᡲᇌɟഏ૾ᆉࡸᴾ

ᵲᵮᵤᵧᵲᴾᵬᵿᶔᶇᶃᶐᵋᵱᶒᶍᶉᶃᶑ૾ᆉࡸểἯὁἏὅ૾ᆉࡸᴾᡲᇌɟഏ૾ᆉࡸᴾᵟᵯᵳᵟᴾҥႻᾂഏΨ൮ဇ༏්ѣᚐௌἩἿἂἻἲᴾ
ʏ්‫ئ‬ᶉᵋ᷁ἴἙἽᴾἣἕἃὊἊᴾίᾎὸᴾᵣᶖᵿᶁᶒᵽᶂᶇᵿᶅᴾ
ᵢᵫᵰᵥᴾ

ἡἱἽἚἝỴὅᘍЗỉؕࡁཞ७ửᚘምᴾ

ίျᘍЗὺ݅ᘍЗὸ‫ݣ‬ᚌ҄ᴾᵬᵭᵮᵧᵡᴾᵫᵿᶖᶕᶃᶊᶊ૾ᆉࡸểЎಊ૾ᆉࡸᴾᡲᇌɟഏ૾ᆉࡸᾎᴾ

ᵲᵢᵢᵤᵲᴾܱᆰ᧓ܱ଺᧓ႆ‫ޒ‬ᇹɟҾྸᚘምᴾ
ᴾျᶍᶐ݅ᘍЗ‫ݣ‬ᚌ҄ᾎᵴᵟᵱᵮᴾ݅ࡇ൮᧙ૠඥ‫˳׍‬ᩓ‫܇‬ನᡯᚘምἩἿἂἻἲᴾἣἕἃὊἊᴾίᵤᵤᵲὺ‫ݣ‬ᚌ҄ὸᴾ

8䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷ኒⴕ೉ኻⷺൻ䝮䉲䊚䊠䊧䊷䉲䊢䊮䈪㜞䈇㔛ⷐⶄᢙ࿕᦭⁁ᘒ䉕ⷐ᳞䈜䉎䉲䊚䊠䊧䊷䉲䊢䊮䈲ዋ䈭䈒䈭䈇㩷 㩷 㩷 㩷 㩷 㩷 ㊂ሶജቇ䈮ၮ䈨䈒‛ℂ䊶ൻቇ䊶↢‛䈠䈱ઁ䈱ಽ㊁㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷 䋪ᚢ⇛ಽ㊁䋲䈱䋳⺖㗴䈻ᣢ䈮ឭଏᷣ䉂㜞ᕈ⢻䋨㜞ਗ೉䋩࿕᦭୯䉸䊦䊋䊷㐿⊒㩷㩷㩷 㩷 㩷 㩷 㩷 㩷 㜞ᐲ䈭䋨ᢙቇ⊛⢛᥊䋫⸘▚ᯏ䋩⍮⼂䈏ᔅⷐ㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷 㩷 ⸘▚⑼ቇജ䈱⸽᣿

9䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷Eigen_K 䋨ታኻ⒓䋩 & Eigenh_K 䋨䉣䊦䊚䊷䊃ⴕ೉ะ䈔䉸䊦䊋䊷䋩㩷㩷䈖䉏䉁䈪䈱㐿⊒ታ❣䋨Altix/SGI, ES/NEC䋩*ES਄䈪ᦨᄢ䋷䋰䋦એ਄䈱ᕈ⢻㩷 䋨SC䋰䋶䊶GBP Finalist䋩
MPI/OpenMP䊊䉟䊑䊥䉾䊄䊝䊂䊦䈪䈱േ૞㩷 㩷 㩷 㜞䈇䉴䉬䊷䊤䊎䊥䊁䉞㩷㩷 㩷 㩷 䋪䋱ਁ䉮䉝䉁䈪䈾䈿䊥䊆䉝㩷 㩷 㩷 䋪䋱䋰ਁ䉮䉝એ㒠䊶䊏䊷䉪䈱䋱䋰䋦䉕⋡ᜰ䈜㩷 㩷 㩷 GPU䋯Cell䋯䉝䉪䉶䊤䊧䊷䉺䈻䈱ኻᔕ䉅⸘↹ਛ㩷 㩷 㩷 ో䈩䈱䊝䊂䊦䈪㜞ᕈ⢻㆐ᚑ䋨ScaLAPACK䉋䉍㜞ㅦ䋩㩷㩷㩷㩷10䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷䂹ᕈ⢻䊌䊤䊜
䊜䉺䈲ScaLAPACK䈫⇣䈭䉍ಽᢔᣇᑼ䈫⁛┙䈮േ⊛䈮೙ᓮ䂹BLAS䉅Ⓧᭂ⊛䈮೑↪䋨ૉ䈚䇮䉴䊧䉾䊄๭䈶಴䈚ᣇᑼ䋩䂹䊙䊙䊦䉼䉮䉝ะ䈔䉝䊦䉯䊥䉵䊛䈱ណ↪
㩷 㩷 㩷 Narrow-band-reduction䉝䊦䉯䊥䉵䊛+Ꮺⴕ೉ะ䈔ಽഀ⛔ᴦᴺ

11䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷

Eigen_S: ᮡḰ⊛䈭䋳Ბ㓏䉝䊦䉯䊥䉵䊛䋺㩷 䊊䉡䉴䊖䊦䉻䊷䋳㊀ኻᩰൻᴺ 㩷 Cuppen䈱ಽഀ⛔ᴦᴺ 㩷 䊊䉡䉴䊖䊦䉻䊷ㅒᄌ឵ᮡḰ⊛䈭ᣇᴺ䈪䈅䉍䇮ታⵝ䈪䈱㜞ㅦൻ䈮ᵈജ䇯
ScaLAPACK䈱pdsyevd䈫ᯏ⢻⊛䈮ห╬ᾢᾲᾷᾴᾼᾴᾉ

Householder
O(N3)Divide&ConqureO(N2)~O(N3)SC2006䊶GBP Finalistᦨᄢല₸70%/ ES

Back-Transformation
O(N3)

12

䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷

Eigen_K: Eigen_S䉕ᕈ⢻਄㒢䊜䊝䊥䊋䊮䊄᏷䈮ᡷ⦟

5㊀ኻⷺൻᴺ䊊䉡䉴䊖䊦䉻䊷



Ꮺኻ⒓ⴕ೉ะ䈔ಽഀ⛔ᴦᴺ䊊䉡䉴䊖䊦䉻䊷ㅒᄌ឵䋲䈧䈱䉝䊦䉯䊥䉵䊛䈱ਗ೉ 䉕ᣂⷙ䈮㐿⊒㩷 Eigen_s䈫ห᭽䈱ታⵝ䈫ว䉒䈞䈩㜞ㅦൻ䉕ታᣉ䇯㩷 ᦝ䈭䉎㜞ㅦൻ䈏ᔅⷐ䈭႐㕙䈪ផᅑ䈘䉏䉎

Scheme:Householder Divide&Conqure Back-Transformation
O(N3)
O(N3)
O(N2)~O(N3)13䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷⿥ਗ೉ൻ䋺㩷
䇸੩䇹䈪ផᅑ䈘䉏䉎ో䈩䈱ਗ೉ൻ䊝䊂䊦䉕ណ↪






MPI: ಽᢔ䊜䊝䊥ਗ೉ಣℂ
OpenMP: ౒᦭䊜䊝䊥ဳ䉴䊧䉾䊄ਗ೉ಣℂ
SIMD: ᯏ᪾⺆䊧䊔䊦ਗ೉ಣℂ

Tofu䉟䊮䉺䊷䉮䊈䉪䊃䉕᦭ല೑↪




䊉䊷䊄䉕䋲
䋲ᰴర䉫䊦䊷䊒䈮ಽ䈔䇮ⴕ䋭೉䈠䉏䈡䉏䈪㐽䈛䈢ㅢା䈮㒢ቯ䋨䋱ኻ䋱ㅢା䇮㓸࿅ㅢା䋩
ㅢା㊂䈭䉌䈶䈮⽶⩄ဋⴧ䉕ᦨዊ䈮䈜䉎䇮䋲ᰴర䉰䉟䉪䊥䉾䉪䊂䊷䉺ಽ
ഀ䉕ណ↪

14

䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷

േ૞ታ❣㩷


T䋲K䉴䊌䉮䊮䈪䈱ታ᷹ᦨᄢ䋴䋰䋹䋶䉮䉝䉁䈪MPI+OpenMP䈱䊊䉟䊑䊥䉾䊄ਗ೉ⅣႺേ૞䋨䋩䋨Flat MPI䈪䉅䋸䋱䋹䋲䉮䉝䉁䈪േ૞ᬌ⸽ᷣ䉂䋩

䋨1024䉮䉝䋩
15䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷
T䋲K䉴䊌䉮䊮䈪䈱ታ᷹ᣢሽ䉸䊐䊃䉡䉢䉝ScaLAPACK䉋䉍䉅㜞ㅦ䈮േ૞䈜䉎΄ΡΖΖΕΦΡͫ͑ΖΝΒΡΤΖΕ͑ΥΚΞΖ͑ΣΒΥΚΠ͙͑΄ΔΒͽͲ΁Ͳʹͼ͠ͶΚΘΖΟΐΤ͚͑

䊉䊷䊄ᢙ

16

䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷䂓㩷 ᕈ⢻੍᷹䋺㪮㪼㪸㫂㩷㪪㪺㪸㫃㫀㫅㪾㩷㩿㪜㫀㪾㪼㫅㪶㪪㪀㩷N

T2KἋἣἅὅỂỉܱยẦỤ๫ም᣽Ẇᡫ̮᣽ί‫ׅ‬ૠὸửምЈẐʮẑἅὅἦἷὊἑỉे‫ܭ‬ẰủỦࣱᏡẦỤਖ਼‫͌ܭ‬ửᚘምታⴕᤨ㑆ផቯ㩿㪪㪚㪘㪄㫋㫉㫀㪄㪉㪛㪄㪩㪋㪀㪉㪈㪇㪇㫇㪻㪾㪼㫄㫄䋨ォㅍ䋩㫋㫉㪼㪻㪈㪶㫍䋨ォㅍ䋩ᤨ㑆㪲⑽㪴

㫋㫉㪼㪻㪈㪶㪸㫌㩿ォㅍ㪀㪈㪉㪇㪇㫋㫉㪼㪻㪈㪶㫌㩿ォㅍ㪀
䈠䈱ઁṶ▚

㪐㪇㪇

㫋㫉㪹㫂㪉䋨Ṷ▚䋩㫋㫉㪹㫂㪈䋨Ṷ▚䋩㪍㪇㪇㫇㪻㫃㪸㪼㪻㪈䋨Ṷ▚䋩㫋㫉㪼㪻㪈㪶㪉㫌㫇㪻㪸㫋㪼㪶㪹㫆㪻㫐䋨Ṷ▚䋩㪊㪇㪇㫋㫉㪼㪻㪈㪶㪸㫌㪶㪹㫆㪻㫐㪈䋨Ṷ▚䋩㪇⚂㪈ਁ ⚂㪉ਁ ⚂㪊ਁ ⚂㪋ਁ ⚂㪍ਁ ⚂㪏ਁ
䊒䊨䉶䉴ᢙ

ታലല₸

㫋㫉㪼㪻㪈㪶㫍㫆㫉㪿㫋㫆䋨ォㅍ㪀

㪈㪌㪇㪇

⚂1ਁ
⚂2ਁ
⚂3ਁ
⚂4ਁ
⚂6ਁ
⚂8ਁ

ታലല₸䋨㪪㪚㪘㪄㫋㫉㫀㪄㪉㪛㪄㪩㪋䋩㫋㫉㪹㪸㫂㫎㫐䋨ォㅍ䋩

㪈㪏㪇㪇

䊒䊨䉶䉴ᢙ

437,149
549,908
693,452
794,145908,238
1,000,0
00

㪋㪌㩼
㪋㪇㩼
㪊㪌㩼
㪊㪇㩼
㪉㪌㩼
㪉㪇㩼
㪈㪌㩼
㪈㪇㩼
㪌㩼
㪇㩼
⚂㪈ਁ

⚂㪉ਁ

⚂㪊ਁ

⚂㪋ਁ

⚂㪍ਁ

⚂㪏ਁ

䊒䊨䉶䉴ᢙ
17

䋲䇮䋨䈽䉺䉴䉬䊷䊦䋩䊤䉟䊑䊤䊥䊷㐿⊒䈻䈱ข䉍⚵䉂㩷䂓㩷 ᕈ⢻੍᷹䋺㪮㪼㪸㫂㩷㫊㪺㪸㫃㫀㫅㪾㩷㩿㪜㫀㪾㪼㫅㪶㪢㪀㩷N

T2KἋἣἅὅỂỉܱยẦỤ๫ም᣽Ẇᡫ̮᣽ί‫ׅ‬ૠὸửምЈẐʮẑἅὅἦἷὊἑỉे‫ܭ‬ẰủỦࣱᏡẦỤਖ਼‫͌ܭ‬ửᚘምታⴕᤨ㑆ផቯ㩿㪪㪚㪘㪄㫋㫉㫀㪄㪉㪛㪄㪙㪉㪄㪩㪋㪀㪉㪈㪇㪇㫇㪻㪾㪼㫄㫄䋨ォㅍ䋩㫋㫉㪼㪻㪈㪶㫍䋨ォㅍ䋩㫋㫉㪼㪻㪈㪶㫍㫆㫉㪿㫋㫆䋨ォㅍ䋩㫋㫉㪼㪻㪈㪶㪸㫌䋨ォㅍ䋩

㪈㪉㪇㪇

ታലല₸

ᤨ㑆㪲⑽㪴

㪈㪌㪇㪇

㫋㫉㪼㪻㪈㪶㫌䋨ォㅍ䋩
䈠䈱ઁṶ▚

㪐㪇㪇

㫋㫉㪹㫂㪉䋨Ṷ▚䋩㫋㫉㪹㫂㪈䋨Ṷ▚䋩㪍㪇㪇㫇㪻㫃㪸㪼㪻㪈䋨Ṷ▚䋩㪊㪇㪇㫋㫉㪼㪻㪈㪶㪉㫌㫇㪻㪸㫋㪼㪶㪹㫆㪻㫐䋨Ṷ▚䋩㫋㫉㪼㪻㪈㪶㪸㫌㪶㪹㫆㪻㫐㪈䋨Ṷ▚䋩㪇⚂㪈ਁ ⚂㪉ਁ ⚂㪊ਁ ⚂㪋ਁ ⚂㪍ਁ ⚂㪏ਁ

⚂1ਁ
⚂2ਁ
⚂3ਁ
⚂4ਁ
⚂6ਁ
⚂8ਁ

ታലല₸䋨㪪㪚㪘㪄㫋㫉㫀㪄㪉㪛㪄㪙㪉㪄㪩㪋䋩㫋㫉㪹㪸㫂㫎㫐䋨ォㅍ䋩

㪈㪏㪇㪇

䊒䊨䉶䉴ᢙ

437,149
549,908
693,452
794,145908,238
1,000,0
00

㪋㪌㩼
㪋㪇㩼
㪊㪌㩼
㪊㪇㩼
㪉㪌㩼
㪉㪇㩼
㪈㪌㩼
㪈㪇㩼
㪌㩼
㪇㩼
⚂㪈ਁ

⚂㪉ਁ

⚂㪊ਁ

⚂㪋ਁ

⚂㪍ਁ

⚂㪏ਁ

䊒䊨䉶䉴ᢙ

䊒䊨䉶䉴ᢙ

18

ᔕ↪଀

SC2012 Gordon Bell Prize㪩㪪㪛㪝㪫㩿㪩㪼㪸㫃㩷㪪㫇㪸㪺㪼㩷㪛㪼㫅㫊㫀㫋㫐㩷㪝㫌㫅㪺㫋㫀㫆㫅㪸㫃㩷㪫㪿㪼㫆㫉㫐㩷㪺㫆㪻㪼㪀㩷㩷• First principle electronic structure simulation• Implemented real space method (not use FFT)
• Main procedure are conjugate-gradient method, ortho-normalization, andsubspace diagonalization• Use MPI + OpenMP• Parallelized in grid space and orbitals• Use DGEMM, and eigensolver library (ScaLAPACK ,EIGEN)

10nm

10nm

Silicon nanowire
39,696 atoms

10648-atom cell of Si crystal and its electron density
19Performance of challenging simulations- Si nw more than 100,000 atoms Si nw, 107,292 atoms, Grids:576x576x192, Orbitals:229,824
Procedure block

Execution time Computation(sec.)
time (sec.)

Communication time (sec.)
Adjacent/spaceGlobal/spaceGlobal/orbitalWait/orbital

Performance(PFLOPS/%)

SCF

5,456.21

4,417.15

83.18

899.05

15.87

40.93

3.08/ 43.63

DIAG3,710.01

3,218.72

27.70

458.874.70-

2.72/ 38.52

MatE

1,084.70

717.4527.70

337.85

4.70-

3.09/ 43.72

EigenSolve

1,322.16

879.61-442.39--0.04/ 0.61

RotV

1,298.16

1,177.15-

121.01--

5.18/ 73.25

CG

209.29

57.66

55.48

96.14

0.01

-

0.05/ 0.74

GS

1,536.90

1,140.76-

344.04

11.16

40.93

4.37/ 61.87
20㪧㪼㫉㪽㫆㫉㫄㪸㫅㪺㪼㩷㫆㪽㩷㪜㪠㪞㪜㪥㩷㩿㪻㪼㫋㪸㫀㫃㪀㩷Si nw, 107,292 atoms , Grids:576x576x192, Orbitals:229,824
Procedure block

TRD1(Householder)DC(Divide&Conqure)
TRKBAKWY(Back Transform)

Total (sec.)

Computation(sec.)

Communication(sec.)

853.443

141.123712.32072.00872.0060.000235.28113.99221.360Communication time of TRD1 is most largest.

Is this performance too low?
••

Matrix element is too small compare to total parallel number.
It is not suitable to map to Tofu network topology.
9 EIGEN: 2 dimension V.S OTHERS: 4 dimension
21㪙㪼㫅㪺㪿㫄㪸㫉㫂㩷㫆㪽㩷㪜㪠㪞㪜㪥㩷
Nearly Weak scaling
Element ofmatrixNumbe Executio
r ofnprocess time(sec.)

151,448288(16x18)1,085.6

694.8

191,232576(24x24)

1,167.1

240,768

1,152(32x36)

302,592

382,464

Performance / core

Breakdown

TRD(sec.)

DC(sec.)

TRKBA
K(sec.)GFLOPS

Efficiency
(%)

84.3

306.3

37.0

29.0

790.1

66.3

311.1

34.6

27.1

1,287.1

846.9

116.6

323.6

31.3

24.5

2,304
(48x48)1,355.1

949.5

76.3

329.2

29.5

23.1

4,608
(64x72)

1,713.7

1,164.
3

184.5

364.8

23.6

18.4

22

3䇮⑔ፉ੐᡿ᓟ䈱⸘▚⑼ቇ䈱ข䉍⚵䉂⸘▚‛ᕈ‛ℂ䈮䉋䉎ේሶജἴኂ䈻䈱ข䉍⚵䉂㩷 㩷 㩷 㩷 䋺࿯ფᳪᨴ䈱ℂ⸃෸䈶䈠䈱ઁᣣᧄේሶജ⎇ⓥ㐿⊒ᯏ᭴䊶䉲䉴䊁䊛⸘▚⑼ቇ䉶䊮䉺䊷䊶䉲䊚䊠䊧䊷䉲䊢䊮ᛛⴚ㐿⊒ቶ㩷ᅏ᧛㓷ᒾ*㩷 ਛ᧛ඳ᮸㩷 ↸↰᣽ᒾ*⑔ፉᛛⴚᧄㇱ㩷 ⑔ፉⅣႺ቟ో䉶䊮䉺䊷⻢ㄉ䋺ᴡ᧛㓶ⴕ䋨ጟጊᄢ䋩䇮ศ↰ༀⴕ䊶ᄢ⽾ᢅᒾ䊶ጊጯഞ䊶ᢪ⮮౏᣿䋨ේሶജᯏ᭴䋩23᡼኿ᕈ‛⾰䈱࿯ფ䈻䈱ๆ⌕䈱‛ℂൻቇᒻᘒ㩷㩷᡼኿ᕈ‛⾰䋨䉶䉲䉡䊛䋩䈱࿯ფᳪᨴ䈱ታᘒ໧㗴䋺 ࿯ფ䈱⴫ጀ䈮䈱䉂ᒝ䈒࿕ቯ䋨⴫ጀ䈎䉌䈱⒖േ䈏ዋ䈭䈇䋩䊶᡼኿⢻䈱ᷫ⴮䈏ㆃ䈇㜞䈇᡼኿⢻ᳪᨴ࿾ၞ䈏ᱷሽ䊶࿾ਅ᳓䈻䈱⒖ⴕ䉇䈠䈱ઁ䈱⒖ⴕ䈮䉋䉎ⵍኂ䈏ዋ䈭䈇ᳪᨴ䈱᜛ᄢ䈏ᛥ೙䈘䉏䉎24

᡼኿ᕈ‛⾰䈱࿯ფ䈻䈱ๆ⌕䈱‛ℂൻቇᒻᘒ㩷㩷᡼኿ᕈ‛⾰䋨䉶䉲䉡䊛䋩䈱࿯ფᳪᨴ䈱ታᘒ⴫ጀ䈮⇐䉁䉎ᳪᨴ╳ᵄᄢ䊶ᕲ↰ᢎ᝼䈮䉋䉎䉴䊤䉟䊄䈎䉌25

᡼኿ᕈ‛⾰䈱࿯ფ䈻䈱ๆ⌕䈱‛ℂൻቇᒻᘒ㩷㩷᡼኿ᕈ‛⾰䋨䉶䉲䉡䊛䋩䈱࿯ფᳪᨴ䈱ታᘒ⴫ጀ䈮⇐䉁䉎ᳪᨴ࿕⌕䈚䈢䉶䉲䉡䊛⣕㔌䈲࿎㔍䋣╳ᵄᄢ䊶ᕲ↰ᢎ᝼䈮䉋䉎䉴䊤䉟䊄䈎䉌26

᡼኿ᕈ‛⾰䈱࿯ფ䈻䈱ๆ⌕䈱‛ℂൻቇᒻᘒ㩷㩷䉶䉲䉡䊛䈱࿯ფᳪᨴ䋨ේሶ䊶ಽሶ䊧䊔䊦䋩䈫䈲䋿䈲䋿࿯ფ䈱䊙䊦䉼䉴䉬䊷䊦ⅣႺ

H2O
Cs

K

ጀ
ጀ⁁᭴ㅧ
䊄䊷䊏䊮䉫
27

᡼኿ᕈ䉶䉲䉡䊛䈱ๆ⌕䉲䊅䊥䉥

http://www.iai.ga.a.u-tokyo.ac.jp/mizo/seminar/cliffslide.pdf

28᡼኿ᕈ‛⾰䈱࿯ფ䈻䈱ๆ⌕䈱‛ℂൻቇᒻᘒ

㩷㩷
㘑ൻ࿯ფ

Ყセ⊛ᣂ䈚䈇࿯ფ⑔ፉ䈱࿯ფ䈮ㄭ䈇䋿

K2+

Mg2+

䋨Ca2+䋩

᳓ಽ䈏᜽䉁䉍⤘Ả䈭☼࿯29᡼኿ᕈ䉶䉲䉡䊛䈱ๆ⌕䉲䊅䊥䉥

30᡼኿ᕈ䉶䉲䉡䊛䈱ๆ⌕䉲䊅䊥䉥 䋨ᣂ䈚䈇࿯ფ☼࿯㋶‛䋩

NHK䊶䉰䉟䉣䊮䉴䉷䊨䋨2/17)

http://www.iai.ga.a.u-tokyo.ac.jp/mizo/seminar/cliffslide.pdf

31᡼኿ᕈ䉶䉲䉡䊛䈱ๆ⌕䉲䊅䊥䉥 䋨ฎ䈇࿯ფ☼࿯㋶‛䋩㘑ൻ䈮䉋䉎䈾䈧䉏

32

᡼኿ᕈ‛⾰䈱࿯ფ䈻䈱ๆ⌕䈱‛ℂൻቇᒻᘒ㩷㩷⸘▚⑼ቇ䈱⎇ⓥኻ⽎䋱䇮࿯ფ☼࿯㋶‛䈱╙৻ේℂ⸘▚䋨ၮᧄ⊛⸘▚䋩䉶䉲䉡䊛䈱቟ቯᕈ䉕ታ㓙䈮⏕⹺䈜䉎䋲䇮䉲䊅䊥䉥䉕⏕䈎䉄䉎䋨╙৻ේℂMD or QM/MM)ๆ⌕೨ᓟ䈱䉣䊈䊦䉩䊷ᄌൻ䈱⹏ଔ䋳䇮䉶䉲䉡䊛⣕㔌䈱䉻䉟䊅䊚䉪䉴‛ℂ⊛䉝䊒䊨䊷䉼䋨䊙䉟䉪䊨ᵄ䉇䊧䊷䉱䊷ᾖ኿䇮ടᾲ䈮䉋䉎ើ⊒䋩ൻቇ⊛䉝䊒䊨䊷䉼䋨䉟䉥䊮੤឵䈠䈱ઁ䋩
33ᓧ䉌䉏䈢⎇ⓥᚑᨐ䋨ᅏ᧛䊶ਛ᧛䊶↸↰䋩᳓ਛ䈱㓁䉟䉥䊮䈱᳓๺䉣䊈䊦䉩䊷Cs䈱ൻቇ⊛
㩷 ቟ቯᕈ

H2O

㓁䉟䉥䊮䈱⟎឵䈮䉋䉎⚿᥏ᱡ䉂䉣䊈䊦䉩䊷

Exchangeg pattern
p

Li

Na

A.KKKK->KKKM

72.08

B.KKKCs->KKMCs

K

Rb

Cs

NH4

30.61

3.87

22.82

-187.74

86.65

39.37

-2.11

7.62

-189.1C.KKKCs->KMKCs

82.7

36.18

-0.45

11.48

-192.23

D.KKCsCs->KMCsCs

90.33

39.75

-4.53

1.61

-196.91

E.KCsCsCs->MCsCsCs

100.36

47.04

-8.28

-8.05

-194.58

F.KKKCs->KKKM

49.26

7.79

-22.82-18.95-210.56

G.KKCsCs->KKMCs

79.03

31.75

-7.62

-9.73

-196.72

H.KCsKCs->KMKCs

71.22

24.7

-11.48-11.93-203.71

I.KCsCsCs->KMCsCs

88.72

38.14

-1.61

-6.14

-198.52

J.CsCsCsCs->MCsCsCs

108.41

55.09

8.05

-0.23

Cs
ฦ⒳㔕Უ᭽
☼࿯㋶‛

䉟䉥䊮੤឵
䉣䊈䊦䉩䊷⴫
䋨5⒳㋶‛䋩

ૉ䈚䊋䊦䉪⁁ᘒ䊋㩷 㩷 䈪䈱੤឵㩷

-186.53䉝䊮䊝䊆䉝䈏䉟䉥䊮੤឵૕䈫䈚䈩䈅䉌䉉䉎㔕Უ᭽☼࿯㋶‛䈮ኻ䈚᦭ല
34䉁䈫䉄䈫੹ᓟ䈱⺖㗴㩷ේሶജ⎇ⓥ㐿⊒䈮䈍䈔䉎⸘▚⑼ቇ䈪䈱䊤䉟䊑䊤䊥䊷㩷







ㅪ┙৻ᰴᣇ⒟ᑼ䇮ⴕ೉࿕᦭୯䋨ኒ䋧⇹䋩䇮㪝㪝㪫╬䈏૶䉒䉏䈩䈇䉎䇯㩷ฦ䉸䊦䊋䊷䈱⒖ᬀ䈫ᕈ⢻⹏ଔ䈲น⢻䇯㩷⃻࿷䈲ኒⴕ೉࿕᦭୯䉸䊦䊋䊷㩿䇸੩䇹䋩㐿⊒䇯㩷ᩭⲢว෻ᓳᴺ䉸䊦䊋䊷㐿⊒䋨䇸੩䇹䋩䋫㪋୚♖ᐲ㪙㪣㪘㪪╬䈱㐿⊒㩷⑔ፉ੐᡿એᓟ䈮ᄢ䈐䈒ᄌ䉒䈦䈢⸘▚⑼ቇ㩷࿯ფᳪᨴ䈱ේሶ䊶ಽሶ䊧䊔䊦䈪䈱ℂ⸃㹢ല₸⊛㒰ᨴᚻᴺ䈱㐿⊒㩷᡼኿ᕈ‛⾰䈱↢૕ౝャㅍ䈱ේሶ䊶ಽሶ䊧䊔䊦䈪䈱ℂ⸃㩷ቢᚑᐲ䈱㜞䈇⸘▚䉮䊷䊄䋨╙৻ේℂ⸘▚䇮ಽሶ゠㆏⸘▚䋩䈏ᔅⷐ䋣㩷

35

4.6 JAXA における数値計算ライブラリの利用状況と OPL 適用評価

JAXA における数値計算ライブラリの利用状況と OPL 適用評価
松尾裕一（JAXA)
本報告では、宇宙航空研究開発機構（JAXA)における数値計算ライブラリ（以下、ライブラリ）の利用動向調査と幾つかのアプリケーション（以下、アプリ）に対する Open Petascale Library(OPL)の適用評価の結果、さらにそれらを踏まえた分析と考察、最後に今後の当該分野におけるライブラリ利用に関する私見を述べる。
JAXA のスパコンシステムは、富士通製 FX1 の 3,008 ノード（SPARC64 VII [2.5GHz、4 コア]、
120TFLOPS、94TB メモリ、DDR IB によるファットツリーIC）を中核とするメインシステム（「JSS-M」と呼ばれる）と、大規模共有メモリやベクトルから成るサブシステムから成る。アプリの 9 割は計算流体力学（CFD)
であり、計算コードの大半は内製1、学術系の大規模ジョブより工学系のパラメータスタディ用ジョブ（並列度は 256 以下）が多い特徴がある。この JSS-M システムにおけるライブラリの利用動向を 2012 年の 1 月
から 2 月にかけて調査2したところ、延べ 350 人程度のシステム利用者のうち、ライブラリ利用者は、（コンパイラが自動で呼び出している場合も含め）5 割以上あるものの、直接ライブラリを呼びだしている利用者はそれほど多くはなく（述べ 50 人程度）、用途も画一的限定的であり、プロセス並列版の利用といった高度利用は行われていないということがわかった。（調査期間がやや短いが）このような結果は、ベクトルの時代から「メモリを節約するためにコードを全部自分で書く」という分野特有の慣習や、コードの長寿命性（一旦作ったコードは 10 年程度は使う）ゆえ、最初からある程度予想されていたとはいえ、そういった流れとは異なる新しい利用も特には進んでいなかった。調査結果を分析する傍ら、ライブラリを利用している 2 件の実アプリを抽出し、FFT と複素連立一次方程式につき OPL の適用を試みた。もともとは SSLII のライブラリを呼び出している部分（dcftm,dlcx）に、OPL の FFTE（zfft1d）と PLASMA（PLASMA_zgesv_incpiv）を適用し、OPL の方が SSLII より若干速いという結果が得られた。適用性は SSLII とほぼ変わらないが、初期化等の用法が全く同じではないので置き換えに伴うデバッグが難しいとか適当なマニュアルがないといった課題も見出された。また、連立一次等の利用者が自分で書いている部分を OPL に置き替えようと試みたが、複雑すぎて短時間では無理と判断した。このあたりにもライブラリを使う難しさの一端が垣間見られた。少ない事例から何ともいえない面もあるが、SSLII も相当に成熟しているライブラリであることを勘案すると、OPL を使うメリットはまずはその高度利用にあると思われる。次に、ライブラリが使われない（使わない）理由について実際に利用者へ無作為にヒアリングしてみたところ、多い順に、①移植性がなくなるから、②使うのが結構面倒だから、③呼び出しの時間的なオーバーヘッドが大きいから、④使い方・メリットがわからないから（「労多くして益少なし」の印象）、⑤余分なメモリを食うから、という回答が得られた。（ただし、これは「利用者がそう思っている」ということであって、実際にそうかどうかは別であることに注意する。）工学
利用は納期が重要なので、様々なプラットフォームへの対応は重要であり、移植性がトップに来たのは意外ではない。面倒さとかオーバーヘッドについては、「オープンソースや市販コード、もっと単純なアルゴリ
ズムを使います」という追加コメントが返って来そうであり、一理ある考え方である。また、流体解析では行
列をきっちり解かなくても何とかなることもあり、ライブラリを使わなければならない理由・局面に乏しいことも利用が進まない原因の一つと考えられる。
このような JAXA におけるライブラリの利用状況は、航空宇宙や流体アプリ特有の事情が少なからず関与していることは確かだが、今後、システムが複雑化し、並列度が高まって行く中では、ライブラリ利用の重要性や必要性は増して行くのではないかと予想される。理由の一つは、ハイブリッド並列等の並列処理の複雑化、GPU 等によるシステム利用の複雑化に対して、利用者がこれら全てに自分でコードを書くという対応はどこかで破綻するのではないかと思われるからである。また、数値解析法の成熟とともに、従来の連立一次、固有値、フーリエ、乱数、補間といった処理に加えて、マルチグリッド、解適合（AMR）、クリロフ
反復、領域分割といった複雑な処理が定番化していることも、これらのライブラリ化がうまく進めばライブラ
1
2

FORTRAN77/90 ベースで書かれ、並列化は MPI または XPF による。最近は市販コードの利用が増えつつある
2012/1/4-2/23 間に JSS-M でコンパイルを行ったユーザを対象に、コンパイル時のログ情報を解析リ利用を加速させる（もっと言えば必須となる）可能性がある。OPL においてこのような複雑処理のライブラリ化が行われることを期待したい。ただし、ライブラリだけ作ってどこかに置いておくだけではだめで、利用者とライブラリを近づける部分を担う人材の育成や、ライブラリの有効利用（特に高度利用）に関するノウ
ハウの蓄積、啓蒙活動が同時に必要であろう。また、単にライブラリ利用ということではなく、もう少し広い意味でのフレームワークやミドルウェアといった考え方を取り入れた方が高度利用にはしっくり来るのではないかとも思われ、ちょっと大袈裟かもしれないが、旧来の数値解析プロセスの見直しや欧米のような分業を取り入れる時に来ているような気もするのだが如何。

SS⎇䊕䉺䉴䉬䊷䊦䊤䉟䊑䊤䊥WG

䎭䎤䎻䎤䬺䬙䬠䭚ᢙ୯⸘▚䮰䭫䮝䮰䮱䬽೑↪⁁ᴫ䬷䎲䎳䎯ㆡ↪⹏ଔ䎃᧻የ⵨৻ቝቮ⥶ⓨ⎇ⓥ㐿⊒ᯏ᭴⎇ⓥ㐿⊒ᧄㇱᢙ୯⸃ᨆ䉫䊦䊷䊒
2012ᐕ7᦬23ᣣ

ႎ๔ౝኈz JAXA䮀䮘䭺䮺䯴JSS䯵䬽᭴ᚑz JSS-M䬺䬙䬠䭚䮰䭫䮝䮰䮱䬽೑↪േะz ታ䭩䮞䮱䬺䬙䬠䭚ㆡ↪⹏ଔ੐଀z ಽᨆ䫺⠨ኤ䫺੹ᓟ䬽ኻᔕ

2

JAXA䮀䮘䭺䮺䯴JSS䯵䬽䭾䮀䮍䮧᭴ᚑ⸘▚䭯䮺䭿䮺ㇱ䎃ᄢⷙᮨਗ೉⸘▚ᯏ䭾䮀䮍䮧䯴䎰䎏䎳䯵䎃౒᦭䮨䮩䮱⸘▚ᯏ䭾䮀䮍䮧䯴䎤䎏䎹䯵䎃䮀䮏䮳䯃䭿ㇱ䎃䮀䮏䮳䯃䭿䭾䮀䮍䮧䎃

ಽᢔⅣႺ⛔วㇱ䎃㆙㓒೑↪䭾䮀䮍䮧䎃ಽᢔ䮎䯃䮆౒᦭䭾䮀䮍䮧䎃㜞ㅦ䮔䮊䮏䮶䯃䭶䎃3⸘▚䭯䮺䭿䮺ㇱ᭴ᚑ䭾䮀䮍䮧ฬ⒓䎃

䎰䭾䮀䮍䮧䎃

䎳䭾䮀䮍䮧䎃

䎦䎳䎸䮆䭫䮞䎃

䎤䭾䮀䮍䮧䎃

䎹䭾䮀䮍䮧䎃

䮀䭲䮰䯃䎃䭾䮀䮍䮧䮆䭫䮞䎃

䮠䭶䮏䮲䎃

䎰䎳䎳䎃

䎶䎰䎳䎃

䮕䯃䮐ᢙ䎃

䎖䎏䎓䎓䎛䎃

䎖䎛䎗䎃

䎔䎃

䎖䎃

䎦䎳䎸ᢙ䎒䮕䯃䮐䎃

䎔䎃

䎔䎃

䎖䎕䎃

䎔䎙䎃

䭺䭩ᢙ䎒䎦䎳䎸䎃
䎋ో䭺䭩ᢙ䎌䎃

䎗䎃䯴䎔䎕䎏䎓䎖䎕䯵䎃䎗䎃䯴䎔䎏䎘䎖䎙䯵䎃

䎗䎃䯴䎔䎕䎛䯵䎃䎔䎃䯴䎗䎛䯵䎃䮛䯃䭶ᕈ⢻䎾䎷䎩䎯䎲䎳䎶䏀䎃䯴䮕䯃䮐䬑䬮䭙䎾䎪䎩䎯䎲䎳䎶䏀䯵䎃

䎔䎕䎓䎃
䯴䎗䎓䯵䎃

䎔䎘䎃
䯴䎗䎓䯵䎃

䎔䎑䎕䎃
䯴䎔䎏䎕䎓䎓䯵䎃

䎗䎑䎛䎃䯴䎔䎏䎙䎓䎓䯵䎃

䮨䮩䮱ኈ㊂䎾䎷䎥䏀䎃䎋䮕䯃䮐䬑䬮䭙䎃䎾䎪䎥䏀䎌䎃

䎜䎗䎃
䯴䎖䎕䯵䎃

䎙䎃
䯴䎔䎙䯵䎃

䎔䎃
䯴䎔䎏䎓䎓䎓䯵䎃

䎖䎃䯴䎔䎏䎓䎓䎓䯵䎃

ን჻ㅢ䎶䎨䎰䎜䎓䎓䎓䎃

䎱䎨䎦䎃䎶䎻䎐䎜䎃⵾ຠฬ䎃ን჻ㅢ䎃䎩䎻䎔䎃

4

JSS-M䬺䬙䬠䭚䮰䭫䮝䮰䮱䬽೑↪േะz ⺞ᩏ
¾ ኻ⽎▸࿐䰆2012/1/4-2012/2/23䬽㑆䬺M䭾䮀䮍䮧䬶䭺䮺䮘䭫䮲䭡䎃ⴕ䬲䬮৻⥸䮭䯃䭽䯴▤ℂ䮭䯃䭽䬙䭗䭂ታⴕ䬯䬠䬽৻⥸䮭䯃䭽䬾ኻ⽎ᄖ䯵
¾ 㓸⸘ᣇᴺ䰆䭾䮀䮍䮧䬽䭺䮺䮘䭫䮲ᤨ䮴䭷ᖱႎ䭡⸃ᨆ䬦䫺᧦ઙ䬺䎃⹥ᒰ䬨䭚䮭䯃䭽䭡ੱᢙන૏䬶䭲䭭䮺䮏䯴䭩䮞䮱䬽9ഀ䬾CFD䯵z 㓸⸘⚿ᨐ᧦ઙੱᢙഀว䭺䮺䮘䭫䮲೑↪⠪ో૕

348

100.0%

ㅙᰴ೑↪

31

8.9%

ਗ೉೑↪

317

91.1%

䮞䮴䮂䮀ਗ೉䬽䭎೑↪

10

2.9%

䮀䮳䮊䮐ਗ೉䬽䭎೑↪

74

21.3%

䮖䭫䮝䮱䮊䮐ਗ೉䬶೑↪

233

67.0%

SSLII, BLAS╬䭡೑↪

192

55.2%

ᐔဋਗ೉ᢙ䋺50䌾200䉮䊮䊌䉟䊤䈏⥄േ䈪
BLAS䊦䊷䉼䊮䈮ᄌ឵䈜䉎䉋䈉䈭႐ว䉅฽䉃䋨଀䈋䈳DAXPY䈭䈬਄૏5䊦䊷䉼䊮䉕ᑧ䈼50ੱ䈏೑↪䋩

5

JSS-M䬺䬙䬠䭚䮰䭫䮝䮰䮱䬽೑↪േะ䯴2䯵z 㓸⸘⚿ᨐ䯴䬳䬴䬜䯵¾ SSLII㑐ᢙᲤ䬽೑↪ੱᢙ䬚䭘਄૏10୘䭡᛽಴䬦䬮⚿ᨐ䎃 䎃 䎃 䯴䮭䯃䭽䬛᣿␜⊛䬺๭䭂಴䬦䬵䬓䭚䭑䬽䯵
㑐ᢙฬ

ੱᢙ

1

SPLV(R8)

11

3.2% 3ᰴspline⵬㑆ᑼ䬺䭗䭚⵬㑆䯂ᢙ୯ᓸಽ2

RANU2(R4)82.3% ৻᭽ੂᢙ䯴0,1䯵䬽↢ᚑ

3

LCX(C16)

5

1.4% ⶄ⚛ⴕ೉䬽ㅪ┙৻ᰴᣇ⒟ᑼ䯴䭶䮰䭭䮏ᴺ䯵

3

ALU(R8)5

1.4% ታⴕ೉䬽LUಽ⸃

5

CEIG2(C16)41.1% ⶄ⚛ⴕ೉䬽࿕᦭୯࿕᦭䮠䭶䮏䮲䯴QRᴺ䯵

5

CFTM(C16)

41.1% ᄙᰴర㔌ᢔဳⶄ⚛䮜䯃䮱䭯ᄌ឵䯴ᷙวၮᐩFFT䯵

5

EIG1(R8)

4

1.1% ታⴕ೉䬽࿕᦭୯࿕᦭䮠䭶䮏䮲䯴2ᲑQRᴺ䯵

5

LUIV(R8)

4

1.1% LUಽ⸃䬤䭛䬮ታⴕ೉䬽ㅒⴕ೉

5

SMLE2(R8)41.1% ᦨዊੑਸ਼ㄭૃᄙ㗄ᑼ䬺䭗䭚ᐔṖൻ䯴ਇ╬㑆㓒㔌ᢔὐ䯵

10

VCFM1(C16)3

0.9% 1ᰴర㔌ᢔဳⶄ⚛䮜䯃䮱䭯ᄌ឵䯴2,3,5䬙䭗䭂7䬽ᷙวၮᐩ䯵

Rank

ഀว

ห䈛ੱ䈎䉅䈚䉏䈭䈇ᯏ⢻

6

JSS-M䬺䬙䬠䭚䮰䭫䮝䮰䮱䬽೑↪േะ䯴3䯵z ⋥ធ䮰䭫䮝䮰䮱䭡๭಴䬦䬵䬓䭚䮭䯃䭽䬾䬑䭍䭙ᄙ䬞䬹䬓䫻䎃 䎃↪ㅜ䭑㒢䭘䭛䫺㜞ᐲ೑↪䬾⷗䭘䭛䬹䬓䫻¾ ㅪ┙৻ᰴ䫺࿕᦭୯䫺䮜䯃䮱䭯䫺ੂᢙ䯂ㄭૃ䯂⵬㑆䫺╬䬺㒢ቯ
¾ 䮞䮴䮂䮀ਗ೉䬽೑↪䬾䬹䬚䬲䬮䯴䮀䮳䮊䮐orㅙᰴ 䬽䭎䯵z 䮭䯃䭽䬽䮰䭫䮝䮰䮱೑↪䬺㑐䬨䭚േะ䬽ಽᨆ䫺⠨ኤ䬾ᓟㅀz 䮭䯃䭽䬛⥄ಽ䬶ᦠ䬓䬵䬓䭚ㇱಽ䯴ㅪ┙৻ᰴ䫺╬䯵䭡䮰䭫䮝䮰䮱䬺⟎䬜឵䬗䭗䬕䬷䬦䬮䬛䫺ⶄ㔀䬨䬝䬵ήℂ䬷್ᢿ&¾ ଀䰆 Q D  R  1  F Q Q   F Q    Q S  forward sweep
2
*

Qi

i

1

i

*



i

Q*i  DijjL i 

11 F Q j
2 jU i 

j

*

j

A

j

ij

Q j   F Q j    A Q j Sij

 backward sweepz 䮰䭫䮝䮰䮱䭡૶↪䬦䬵䬓䬵䫺ㆇ↪஥䬶ታⴕⅣႺ䭡଻᦭䬦䬵

䬓䭚2ઙ䬽ታ䭩䮞䮱䬺㑐䬦䬵䫺OPL䬽ㆡ↪⹏ଔ䭡ታᣉ䎃
7

OPL䬽ㆡ↪⹏ଔ䭸䯃䮀⟎឵䬗ኻ⽎SSLII㑐ᢙኻᔕ䬨䭚䭱䯃䮞䮺䮡䮆䮀䭸䯃䮲䮰䭫䮝䮰䮱

䭩䮞䮱ฬ

ฬ೨

ᯏ⢻

ฬ೨

CHANL

dcftm
ᄙᰴర㔌ᢔဳⶄ⚛䮜䯃䮱䭯 zfft1d(FFTE)(SSL2ၮᧄ) ᄌ឵䯴ᷙวၮᐩ䯵

Crux

dlcx
ⶄ⚛ⴕ೉䬽ㅪ┙䯽ᰴᣇ⒟(SSL2ၮᧄ) ᑼ䯴䭶䮰䭭䮏ᴺ䯵PLASMA_
zgesv_incpiv(PLASMA)ᯏ⢻䯽ᰴరⶄ⚛䮜䯃䮱䭯ᄌ឵䯴ᷙวၮᐩ䯵

ⶄ⚛ⴕ೉䬽ㅪ┙䯽ᰴᣇ⒟ᑼ

ᵈ䋩 Cruxਛ䈪䈲࿕᦭୯䉅๭䉖䈪䈇䉎䋨DCEIG2䋩䈏䇮OPL஥䈱੐ᖱ䈪੹࿁䈱ㆡ↪⹏ଔ䈮䈲㑆䈮ว䉒䈝

8

CHANL䬽႐วz CHANL䭺䯃䮐᭎ⷐ
¾ 3D㕖࿶❗ᕈ䮑䮚䭯䮀䮏䯃䭶䮀ᣇ⒟ᑼ䬺䭗䭚⋥ធ䭾䮦䮬䮳䯃䭾䮮䮺¾ ᤨ㑆Ⓧಽ䰆2ᰴ䭶䮰䮺䭶䮒䭺䮲䮄䮺+3ᰴ䮲䮺䭹䭶䮊䮆¾ ⓨ㑆㔌ᢔൻ䰆4ᤨ♖ᐲਛᔃᏅಽ+2ᤨ♖ᐲਛᔃᏅಽ
¾ ࿶ജ䮤䭩䮄䮺䮄䮲䮗䰆FFT¾ 䮞䮴䮂䮀ਗ೉䰆XPFortran䯸䮀䮳䮊䮐ਗ೉䰆VISIMPACT+OMP¾ 䮞䮴䭷䮰䮧㊂䰆20,000ⴕ⒟ᐲ
z FFT䬺OPL䯴FFTE䯵䭡ㆡ↪
z OPL FFTE䬽᭎ⷐ¾ V4.1䯴⃻࿷䬾V5.0䯵䫺Fortran77¾ OMP䫺OMP+MPI¾ ⶄ⚛ᢙ/ⶄ⚛ᢙ䫺ၮᐩ=2,3,5䫺ᰴర=1,2,3䫺୚♖ᐲ

9

CHANL䬽႐ว䯴2䯵
z SSLII䬷FFTE䬽Ყセ
¾ ᯏ⢻

¾ ᒁᢙ

SSLII䯴DCFTM䯵

FFTE䯴zfft1d䯵

ೋᦼൻ

ਇⷐ

ᔅⷐ

ᱜⷙൻ

䬦䬹䬓

䬨䭚

䮗䮊䮜䭨

ਇⷐ

ᔅⷐ

䭯䮰䯃䭺䯃䮐䬑䭙䬹䬦

SSLII䯴DCFTM䯵

FFTE䯴zfft1d䯵

╙1ᒁᢙ౉ജ䰆㈩೉䬽ታㇱ䎃಴ജ䰆ᄌ឵ᓟ㈩೉䬽ታㇱ䯴਄ᦠ䬜䯵౉ജ䰆ⶄ⚛㈩೉䎃಴ജ䰆ᄌ឵ᓟ䬽ⶄ⚛㈩೉䯴਄ᦠ䬜䯵

╙2ᒁᢙ౉ജ䰆㈩೉䬽⯯ㇱ䎃಴ജ䰆ᄌ឵ᓟ㈩೉䬽⯯ㇱ䯴਄ᦠ䬜䯵㈩೉䬽ⷐ⚛ᢙ

╙3ᒁᢙ㈩೉䬽ⷐ⚛ᢙ

0䰆ೋᦼൻ䫺-1䰇ᱜᄌ឵䫺1䰆ㅒᄌ឵╙4ᒁᢙᰴర䯴ᧄ੐଀䬶䬾1䯵䮗䮊䮜䭨↪ⶄ⚛㈩೉
ᄌ឵䬨䭚㈩೉䬽2୚䬽ⷐ⚛ᢙ䭡ᜰቯ

╙5ᒁᢙ

1䰆ᱜᄌ឵䫺-1䰆ㅒᄌ឵╙6ᒁᢙ䭯䮰䯃䭺䯃䮐

10

CHANL䬽႐ว䯴3䯵z ᦠ䬜឵䬗଀䯴ᱜᄌ឵䯵ㆡ↪೨
:1%.524'#&&1+5'&1,,)
  ࡮࡮࡮
ECNNIGVVQF VU QORFQ
&1++)
&1--)9#< - 9$ +- 9#< - &%106+07'%%#..&%(6/ 9#<9#<0'+%10 %&1--) 9$ +- 9#< - 9$ +- 9#< - %106+07'
%106+07'
QORGPFFQECNNIGVVQF VU   ࡮࡮࡮:1%.'0&524'#&ㆡ↪ᓟ
:1%.524'#&&1+5'&1,,)
  ࡮࡮࡮
ECNNIGVVQF VU QORFQ
&1++)
&1--)9#< - 9$ +- %106+07'%%#..\HHVF 9#<0'  9#< ೋᦼൻ
%#..\HHVF 9#<0'  9#<  &A((6%&1--) 9$ +- FTGCN 9#< - 9$ +- FKOCI 9#< - %106+07'
%106+07'
QORGPFFQECNNIGVVQF VU   ࡮࡮࡮
:1%.'0&524'#&

11

CHANL䬽႐ว䯴4䯵z ᦠ䬜឵䬗଀䯴ㅒᄌ឵䯵ㆡ↪೨
:1%.524'#&&1+5'&1,,)
  ࡮࡮࡮
ECNNIGVVQF VU QORFQ
&1++)
&1--) 9#< - 9$ +- 9#< - 9$ +- %106+07'
&1--)
9#< -) - 9$ +- 9#< -) - 9$ +- %106+07'%%#..&%(6/ 9#<9#<0'+%10 %&1--)9$ +- 9#< - %106+07'
%106+07'
QORGPFFQECNNIGVVQF VU   ࡮࡮࡮:1%.'0&524'#&ㆡ↪ᓟ
:1%.524'#&&1+5'&1,,)
  ࡮࡮࡮
ECNNIGVVQF VU QORFQ
&1++)
&1--) 9#< - FEORNZ 9$ +- 9$ +%106+07'&1--)
9#< -) - FEORNZ 9$ +- 9$ +- %106+07'%%#..\HHVF 9#<0'  9#<  &A((6%&1--)9$ +- 0'  FTGCN 9#< - %106+07'
%106+07'
QORGPFFQECNNIGVVQF VU   ࡮࡮࡮
:1%.'0&524'#&

12

CHANL䬽႐ว䯴5䯵z ታⴕ᧦ઙ
¾ ⸘▚ⷙᮨ
¾ 䮞䮴䮂䮀ᢙ¾ 䮀䮳䮊䮐ᢙ

512䯴䮀䮳䮊䮐䯵䐙320䯴XPF䯵䐙256䯴FFT䯵
32
4z ㆡ↪⚿ᨐ
䭼䮝䮲䯃䮈䮺䎃

ㆡ↪೨䎃

ㆡ↪ᓟ䎃

䎦䎲䎶䎩䎩䎷䎕䎧䎩䯴ᱜᄌ឵䯵䎃䎓䎑䎓䎖䎖䎓䎘䎚⑽䎃

䎓䎑䎓䎕䎛䎜䎘䎗⑽䎃

䎦䎲䎶䎩䎩䎷䎕䎧䎵䯴ㅒᄌ឵䯵䎃䎓䎑䎓䎖䎔䎜䎗䎘⑽䎃

䎓䎑䎓䎓䎙䎚䎕䎖⑽䎃

¾ ㅒᄌ឵䬛ㅦ䬓䬽䬾䯸ೋᦼൻ䬛ᔅⷐ䬹䬓䬮䭐

¾ 䮀䮳䮊䮐ਗ೉ൻ䬾䬤䭛䬵䬓䭚

13

Crux䬽႐ว
z Crux䭺䯃䮐䬽᭎ⷐ¾ ⿥㖸ㅦ䬽౞㍙䭍䭞䭙䬽ゲኻ⽎ᵹ䫺⿥㖸ㅦ䬽ਥ⠢਄㕙䬽ᵹ䭛¾ NSᣇ⒟ᑼ䭡ਗⴕᵹ䬷✢ᒻᕈ䬽ㄭૃ䭡ᣉ䬦䫺eNᴺ䬶⸃䬞¾ ࿕᦭୯⸘▚䰆䭺䮴䭸䯃䭾䮮䮺ᴺ¾ 䮞䮴䮂䮀ਗ೉䰆䬹䬦䯸䮀䮳䮊䮐ਗ೉䰆VISIMPACT+OMP¾ 䮞䮴䭷䮰䮧㊂䰆1,000ⴕ⒟ᐲ
z ⶄ⚛ㅪ┙৻ᰴᣇ⒟ᑼ䬺OPL䯴PLASMA䯵䭡ㆡ↪
z OPL PLASMA䬽᭎ⷐ
¾ V2.4.2䯴⃻࿷䬾V2.4.5䯵䫺Fortran77/C¾ pthread¾ ⶄ⚛ᢙ/ⶄ⚛ᢙ䫺୚♖ᐲ
14

Crux䬽႐ว䯴2䯵
z PLASMA䬷SSLII䬽Ყセ㪪㪪㪣㪠㪠㩿㪛㪣㪚㪯㪀
¾ ᯏ⢻

¾ ᒁᢙ

㪧㪣㪘㪪㪤㪘㩿㪧㪣㪘㪪㪤㪘㪶㪱㪞㪜㪪㪭㪶㪠㪥㪚㪧㪠㪭㪀䍏䍷䍘䍼䍶䍛䍼䍯

䉪䊤䉡䊃ᴺ

㪣㪬ಽ⸃

ೋᦼൻ

ਇⷐ

ᔅⷐ

ฝㄝ

䊔䉪䊃䊦䋨䋱ᰴర䋩㈩೉䋨䋲ᰴర䋩䊋䉾䊐䉜ᔅⷐ䋨໧㗴ⷙᮨ䈮ว䉒䈞䈢㈩೉䉕ᜰቯ䋩ᔅⷐ䋨ኾ↪䊤䉟䊑䊤䊥䈪⏕଻䈚䈢㈩೉䉕ᜰቯ䋩ᒁᢙ㩷㪪㪪㪣㪠㪠㩿㪛㪣㪚㪯㪀㩷㪧㪣㪘㪪㪤㪘㩿㪧㪣㪘㪪㪤㪘㪶㪱㪞㪜㪪㪭㪶㪠㪥㪚㪧㪠㪭㪀㩷╙䋱ᒁᢙ㩷ଥᢙⴕ೉㪘䋨ⶄ⚛ᢙ䋩㩷ଥᢙⴕ೉㪘䈱ᰴᢙ㩷╙䋲ᒁᢙ㩷ଥᢙⴕ೉㪘䈱ᢛวኸᴺ䋨㻢㪥䋩㩷ฝㄝ䈱㈩೉䌢䈱ᰴᢙ䋨㻢㪇䋩㩷╙䋳ᒁᢙ㩷ଥᢙⴕ೉㪘䈱ᰴᢙ㫅㩷౉ജ䋺ଥᢙⴕ೉㪘䋨ⶄ⚛ᢙ䋩㩷಴ജ䋺㪣㪬ಽ⸃⚿ᨐ㩷╙䋴ᒁᢙ㩷౉ജ䋺ቯᢙ䊔䉪䊃䊦䌢䋨ⶄ⚛ᢙ䋩㩷಴ജ䋺⸃䊔䉪䊃䊦㫏㩷ଥᢙⴕ೉㪘䈱ᢛวኸᴺ䋨㻢㪥䋩㩷╙䋵ᒁᢙ㩷䊏䊗䉾䊃䈱⋧ኻ㔖್ቯ୯䋨㻢㪇㪅㪇䋩㩷಴ജ䋺ⴕ೉ᑼ䉕᳞䉄䉎䈢䉄䈱ᖱႎ㩷╙䋶ᒁᢙ㩷䋱䋺䋱⚵⋡䈱ᣇ⒟ᑼ䉕⸃䈒䇯㩷䋲䋺䋲⚵⋡એ㒠䈱ᣇ⒟ᑼ䉕⸃䈒䇯㩷಴ജ䋺䊏䊗䉾䊃ᖱႎ㩷╙䋷ᒁᢙ㩷಴ജ䋺ⴕ೉㪘䈱ⴕ೉ᑼ䉕᳞䉄䉎䈢䉄䈱ᖱႎ㩷ฝㄝ䈱㈩೉䌢㩷╙䋸ᒁᢙ㩷૞ᬺ㗔ၞ䋨ᄢ䈐䈘㫅䈱ⶄ⚛ᢙဳ䋱ᰴర㈩೉䋩㩷ฝㄝ䈱㈩೉䌢䈱ᢛวኸᴺ㩷╙䋹ᒁᢙ㩷૞ᬺ㗔ၞ䋨ᄢ䈐䈘㫅䈱৻ᰴర㈩೉䋩㩷㵪㩷╙䋱䋰ᒁᢙ㩷

䉣䊤䊷䉮䊷䊄㩷

㵪㩷

15

Crux䬽႐ว䯴3䯵z ᦠ䬜឵䬗଀
ㆡ↪೨

ㆡ↪ᓟ


KPENWFGRNCUOCHJKPVGIGTEQTGKPHQ
RCTCOGVGT EQTG 
KPVGIGT JN  JRKX  GZVGTPCN2.#5/#AKPKV2.#5/#A#..1%A914-52#%'A<)'58A+0%2+8
2.#5/#A<)'58A+0%2+8

GZVGTPCN2.#5/#A&'#..1%A*#0&.'2.#5/#A(+0#.+<'
KUY

FQLPPECNN2.#5/#A+0+6 EQTGKPHQ FQKPP
\G K D KL 
ECNN2.#5/#A#..1%A914-52#%'A<)'58A+0%2+8 PPJNJRKXKPHQ GPFFQKECNNVKOGTAUVC  
ECNNVKOGTAUVC  
ECNNFNEZ CPPPP\GGRU\KUYKU\XYKXYKEQP  ECNN2.#5/#A\IGUXAKPERKX PPPPCPPJNJRKXDPPKPHQ ECNNVKOGTAGPF  
ECNNVKOGTAGPF  
FQLPPFQKPPFQKPP
\E KL \G K 
GPFFQK\E KL D KL 
KUY
GPFFQ
GPFFQLGPFFQ࡮࡮࡮
ECNN2.#5/#A&'#..1%A*#0&.' JNKPHQ ECNN2.#5/#A&'#..1%A*#0&.' JRKXKPHQ ECNN2.#5/#AHKPCNK\G KPHQ 16Crux䬽႐ว䯴4䯵z ታⴕ᧦ઙ
¾ ⸘▚ⷙᮨ
¾ 䮞䮴䮂䮀ᢙ¾ 䮀䮳䮊䮐ᢙ

118䐙61䐙84䬹䬦
4z ㆡ↪⚿ᨐ
ㆡ↪೨䎃ㆡ↪ᓟ䎃䎔䎛䎗䎑䎓䎚⑽䎃

䎔䎓䎗䎑䎗䎓⑽䎃

¾ 㜞ㅦൻ䬾䫺๭䭂಴䬦䬛৻࿁䬺ᷫ䬲䬮䬮䭐¾ ኻ⽎ㇱಽ䬾䮀䮳䮊䮐ਗ೉ൻ䬤䭛䬵䬓䬹䬓¾ ᧦ઙ䬺䭗䭙⺋Ꮕ⊒↢¾ PLASMAౝㇱ䬶䭯䮰䯃ᬌ಴䬦䬮႐ว䬽䮎䮗䮊䭷䬛࿎㔍¾ 㘃ૃ䬽㑐ᢙ䬛ⶄᢙ䬑䭙䫺䬸䭛䬛ㆡಾ䬚್䭙䬺䬞䬓17䬢䬢䭍䬶䬽䭍䬷䭐
z OPL䭡ㆡ↪䬨䭚䬺䬑䬮䭙䫺䮭䯃䭽䭩䮞䮱䬺䬙䬠䭚䮰䭫䮝䮰䮱䬽೑↪⁁ᴫ䭡⺞ᩏ¾ 䮰䭫䮝䮰䮱䬛䬑䭍䭙૶䭞䭛䬵䬓䬹䬓䬢䬷䬛್᣿䯴䬬䭢䬹䭑䬽?䯵
¾ 䮞䮴䮂䮀ਗ೉䭡฽䭏䭗䬕䬹ᄢၞ⊛䬹䮰䭫䮝䮰䮱䬽೑↪䬾䬹䬓z ૶䭞䭛䬵䬓䭚䭑䬽䬚䭘䫺FFT䬷ⶄ⚛ㅪ┙৻ᰴᣇ⒟ᑼ䬺䬳䬜OPL䭡ㆡ↪⹏ଔ¾ ㆡ↪ᕈ䯴૶䬓ൎᚻ╬䯵䬺䬳䬜䫺ᄢ䬜䬹໧㗴䬹䬦• ⚿ᨐ䬛㆑䬕႐ว䬽ኻᔕ䬺࿎┆䯸䮥䮒䮬䭩䮲䭑䬹䬓¾ 㜞ㅦൻ䬾䭫䮥䭫䮈䯴䭍䬑䫺䬢䭢䬹䭑䬽?䯵

z એਅ䬶䫺䮰䭫䮝䮰䮱䬛૶䭞䭛䬹䬓ℂ↱䭓੹ᓟ䬽ኻᔕ䬺䬳䬜⠨ኤ
18૶䭞䭛䬹䬓ℂ↱䬺䬳䬓䬵⠨䬗䬵䭎䬮z 䮨䮩䮱䬛ዋ䬹䬚䬲䬮䮠䭶䮏䮲ᤨઍ䬚䭘䬽ฬᱷ䬶䫺ዋ䬦䬶䭑䮨䮩䮱䭡ല₸⦟䬞૶䬕䬮䭐䬺⥄ಽ䬶ᦠ䬞≻䯴⠌ᘠ䯵䬛ᱷ䬲䬵䬓䭚z ᧂ䬯䬺૛ಽ䬹䮨䮩䮱䭡㘩䬕䯴䬷ᕁ䬲䬵䬓䭚䯵䬽䬶૶䬓䬮䬞䬹䬓¾ 䮶䯃䭶㈩೉䬛䮨䮩䮱䭡㘩䬕z ๭䭂಴䬦䭡䬨䭚ᤨ㑆⊛䭱䯃䮗䯃䮟䮊䮐䬛ᄢ䬜䬓䯴䬷ᕁ䬲䬵䬓䭚䯵¾ ೋᦼൻ䫺㑐ᢙ䬽ᒁᢙ㈩೉䭇䬽䭺䮛䯃䫺㑐ᢙ๭䭂಴䬦䫺⚿ᨐ䬽䭺䮛䯃䬹䬸z ૶䬕䬽䬛⚿᭴㕙ୟ
¾ 䮥䮒䮬䭩䮲䯴ᣣᧄ⺆䯵䬛䬹䬓䫺ᒁᢙ䬽ᗧ๧䬽ℂ⸃䬺ᤨ㑆䬛䬚䬚䭚¾ ૶䭞䬹䬠䭛䬿䬹䭘䬹䬓䭊䬸ⶄ㔀䬹ಣℂ䯴or䮢䮊䮏䮀䮤䮊䮏䯵䬛䬹䬓z ૶䬓ᣇ䯂䮨䮱䮊䮏䭡⍮䭘䬹䬓䯴䬄ഭᄙ䬞䬦⋉ዋ䬹䬦䬅⊛ශ⽎䯵
¾ ᧄᒰ䬺ㅦ䬓䬽䬚䰋䎃 ᧄᒰ䬺૶䬓䭓䬨䬓䬽䬚䰋z ૶䬕䬷⒖ᬀᕈ䬛䬹䬞䬹䭚z 㐿⊒஥䫺䱨䲇䱪䱚஥䬺૶䬲䬵䭑䭘䬕䬮䭐䬽໪⫥ᵴേ䬛ਇ⿷z 䮝䮰䮊䭶䮣䮊䭶䮀⊛䬺䬹䭚䬽䭡ෂᗋ䫺╵䬗䬛ᱜ䬦䬓䬚ਇ቟19૶䭞䭛䬹䬓ℂ↱䬺䬳䬓䬵⡞䬓䬵䭎䬮z 䮨䮩䮱䬛ዋ䬹䬚䬲䬮䮠䭶䮏䮲ᤨઍ䬚䭘䬽ฬᱷ䬶䫺ዋ䬦䬶䭑䮨䮩䮱䭡ല₸(0)䎃

⦟䬞૶䬕䬮䭐䬺⥄ಽ䬶ᦠ䬞≻䯴⠌ᘠ䯵䬛ᱷ䬲䬵䬓䭚z ᧂ䬯䬺૛ಽ䬹䮨䮩䮱䭡㘩䬕䯴䬷ᕁ䬲䬵䬓䭚䯵䬽䬶૶䬓䬮䬞䬹䬓(1)¾ 䮶䯃䭶㈩೉䬛䮨䮩䮱䭡㘩䬕¾ ⸃䬜䬮䬓ᄌᢙ㈩೉䬺ኻ䬦䬵⋥ធ Ꮕಽᑼ䭡ᦠ䬜ਅ䬦䬵✭๺Ṷ▚䭡䬚䬠䬵䬦䭍䬲䬵䬓䭍䬨䫻 䬢䭛䬯䬷೨ಣℂ䭡䬚䬠䬛䬮䬓䬽䬶䫺ᧄ᧪䫺෼᧤䬛ᖡ䬞䬹䬲䬵䬸䬕䬚䫺 䬷䬓䬕䬢䬷䭑䬑䭚䬽䬶䬨䬛䫺ᐘ䬓೨ಣℂ䬦䬹䬞䬷䭑䫺 ᢙච࿁⒟ᐲ䬶ᔅⷐ䬹♖ᐲ䭍䬶෼᧤䬨䭚䬽䬶䫺䭗䬓䫺䬷䬓䬕䬢䬷䬺䬦䬵䬓䭍䬨䫻z ๭䭂಴䬦䭡䬨䭚ᤨ㑆⊛䭱䯃䮗䯃䮟䮊䮐䬛ᄢ䬜䬓䯴䬷ᕁ䬲䬵䬓䭚䯵(2)¾ ೋᦼൻ䫺㑐ᢙ䬽ᒁᢙ㈩೉䭇䬽䭺䮛䯃䫺㑐ᢙ๭䭂಴䬦䫺⚿ᨐ䬽䭺䮛䯃䬹䬸z ૶䬕䬽䬛⚿᭴㕙ୟ

(3)¾ 䮥䮒䮬䭩䮲䯴ᣣᧄ⺆䯵䬛䬹䬓䫺ᒁᢙ䬽ᗧ๧䬽ℂ⸃䬺ᤨ㑆䬛䬚䬚䭚¾ ૶䭞䬹䬠䭛䬿䬹䭘䬹䬓䭊䬸ⶄ㔀䬹ಣℂ䯴or䮢䮊䮏䮀䮤䮊䮏䯵䬛䬹䬓¾ 㘩䭞䬩ህ䬓䬯䬠䬚䭑䬦䭛䭍䬪䭢20૶䭞䭛䬹䬓ℂ↱䬺䬳䬓䬵⡞䬓䬵䭎䬮䯴2䯵z ૶䬓ᣇ䯂䮨䮱䮊䮏䭡⍮䭘䬹䬓䯴䬄ഭᄙ䬞䬦⋉ዋ䬹䬦䬅⊛ශ⽎䯵

(2)¾ ᧄᒰ䬺ㅦ䬓䬽䬚䰋䎃 ᧄᒰ䬺૶䬓䭓䬨䬓䬽䬚䰋z ૶䬕䬷⒖ᬀᕈ䬛䬹䬞䬹䭚(4)¾ SSLII䬽ઍ䭞䭙䬺䮜䮱䯃䬽䮋䯃䮲䭡૶䬲䬵䬓䭚䫻SSLII䬽ᣣᧄ⺆䬽䮥䮒䮬䭩䮲䬛ℂ⸃ਇ⢻䬶䫺 ৻ᣇ䫺 䮜䮱䯃䬽䮋䯃䮲䬽⧷⺆䬽䮐䭴䮬䮨䮺䮏䬾⸃⺑䬛ৼካ䬶್䭙䭓䬨䬚䬲䬮¾ ੹䬽ᦨᄢ䬽ℂ↱䬾⒖ᬀᕈ䬶䬨䫻Fortran䬽႐ว䬾ᚻర䬽PC䬯䬷 Gnu 䬚 Intel 䬶䫺Linux↪ን჻ㅢFortran 䬾32bit 䬾⚳䭞䭙䫺64bit 䭑v1䬽䭍䭍㐳䭘䬞ᣂ䬦䬞䬹䬲䬵䭍䬪䭢䬚䭘䬼䫻ᤄ䫺NWT䬛䬶䬜䬮䬢䭜䬺䫺ਗ೉䮰䭫䮝䮰䮱䬛䬿䬲䬰䭙ᢛ஻䬤䭛䬵䬓䬮䭘䫺䬶䬜䭚䬯䬠䬬䬲䬰䭡૶䬕䫺䬷䬓䬕䬢䬷䬺䬹䬲䬵䬓䬮䬚䭑⍮䭛䭍䬪䭢䬛䫺 ⥄ಽ䬽䮄䯃䮀䮞䮴䭷䮰䮧䭡⒖ᬀ䬨䭚䬢䬷䬶Ά䬗ዧ䬜䭍䬦䬮䫻䬢䭛䬚䭘䬯䬷ోㇱC++䭓JAVA䬶ᦠ䬞䬷ታലᕈ⢻䬛಴䬹䬓ಽ䭡䮰䭫䮝䮰䮱䬶 䭲䮗䯃䬶䬜䭚䫺䬹䭢䬵㘑䬺䬹䭚䬷䬓䬓䬽䬶䬦䭖䬕䬼¾ ੹䬾䫺䮘䮄䭺䮺䬶䬶䬜䭚⸘▚䭑䬑䭙䭍䬨䫻䬤䭘䬺PC䭶䮰䮀䮆䫺䮂䮺䮆䯃䮥䭾䮺䬷ⶄᢙ䬽䮞䮰䮊䮏䮢䯃䮧䭡૶䬕䬢䬷䬾䬑䬮䭙䭍䬗䬽ᤨઍ䬹䬽䬶䫺⒖ᬀᕈ䬛଻⸽䬤䭛䬵䬓䬹䬓䬷૶䬓㔍䬓䬷ᕁ䬓䭍䬨䫻ઍ⴫⊛䬹䮞䮰䮊䮏䮢䯃䮧䭡✂⟜䬦䬮䮰䭫䮝䮰䮱䭑䬑䭙䭍䬨䬛䫺ታ䬾․ቯ䬽䮞䮰䮊䮏䮢䯃䮧ะ䬠䬺䬦䬚䮈䮬䯃䮒䮺䭷䬤䭛䬵䬓䬹䬓䫺䬷䬓䬕႐ว䭑ᄙ䬓䭗䬕䬶䬨䫻
¾ 䭫䮺䮆䯃䮜䭮䯃䮀䬾਎䬽ਛ䬺ᄙ䬞಴࿁䬲䬵䬓䭚䭑䬽䬷ห䬧䬶䫺 ታⵝ䬾⸘▚ᯏ࿕᦭䬽ᯏ⢻䭡㚟૶䬦䬵ᩰᲑ䬺ㅦ䬓䭑䬽䭡૞䭚䫺 䬷䬓䬕ᣇะ䬶ദജ䬦䬵㗂䬠䭛䬿䬷ᕁ䬓䭍䬨䫻
21૶䭞䭛䬹䬓ℂ↱䬺䬳䬓䬵⡞䬓䬵䭎䬮䯴3䯵
z 䬬䬽ઁ

(2)

¾ 䮰䭫䮝䮰䮱䬶৻⇟ലᨐ䬾䬑䭚䬽䬾䫺ⴕ೉䭡⸃䬓䬮䭙䫺 ࿕᦭୯䭡᳞䭐䬮䭙䬷䬓䬕ㇱಽ䬶䬨䬛䫺 ᵹ૕⸘▚䬾ၮᧄ⊛䬺㑆ធᴺ䬹䬽䬶䫺⋥ធᴺ䬽䮰䭫䮝䮰䮱䭡 䬬䬽䭍䭍૶䬕䬽䬾ዋ䬹䬓䬷ᕁ䬓䭍䬨䫻㑆ធᴺ䬾䬬䬢䭍䬶 䭺䯃䮎䭪䮺䭷䭑㔍䬦䬞䬹䬓䬶䬨䬦䫻 ࿕᦭୯䭑ᵹ૕䬶䬾䬑䭍䭙૶䭞䬹䬓䬷ᕁ䬓䭍䬨䫻 ᭴ㅧ⸃ᨆ䬽FEM䭺䯃䮐䬯䬲䬮䭘䫺⛘ኻ䮰䭫䮝䮰䮱䭡૶䬕䬷ᕁ䬓䭍䬨䫻
¾ CFD䬽႐ว䫺✢ᒻઍᢙ䬺ઍ⴫䬤䭛䭚䮞䮱䮦䮍䭪䮝䬹䮰䭫䮝䮰䮱䬶䬾䬹䬞䫺ઍ⴫⊛䬹
⒖ᵹ䮀䭴䯃䮧䭓ੂᵹ䮩䮎䮲䬽䮰䭫䮝䮰䮱䬛䬑䭛䬿䫺ᗧᄖ䬷㔛ⷐ䬛䬑䭚䭗䬕䬹᳇䬛䬦䭍䬨䫻 ✢ᒻઍᢙ䯴╬䯵䬽䮰䭫䮝䮰䮱䬽႐ว䫺Navier-Stokesᣇ⒟ᑼ䬛5ⴕ5೉䬽䮝䮴䮊䭶ⴕ೉䭡ⷐ⚛䬷䬨䭚䬢䬷䬚䭘䫺ᄢⷙᮨ䬹⇹ⴕ೉䭡ᛒ䬕䮰䭫䮝䮰䮱䬶䬾䬹䬞䫺ዊⷙᮨ䬹ⴕ೉䭡৻ᐲ䬺ᄢ㊂䬺ല₸䭗䬞ᛒ䬗䭚䮰䭫䮝䮰䮱䬽ᣇ䬛䫺㔛ⷐ䬛䬑䭚䭗䬕䬹᳇䬛䬦
䭍䬨䫻․䬺䫺䮨䮒䭫䭺䭩䭓GPGPU䭑ⷞ㊁䬺౉䭛䭚䬷䫺଀䬗䬿ᄢⷙᮨ⇹ⴕ೉䬽ㅒⴕ೉䭡⋥ធᴺ䬺䭗䭙᳞䭐䭚䬢䬷䬾⃻ታ⊛䬶䬾䬹䬞䬵䭑䫺ዊⷙᮨⴕ೉䬽ㅒⴕ೉䭡ᄢ㊂䬺᳞䭐䭚䬢䬷䬾ᗧᄖ䬷น⢻䬺䬹䭙䬬䬕䬹᳇䬛䬦䭍䬨䫻䯴䭺䭩ᢙn䬷⸘▚䬨䭈䬜ዊⴕ೉䬽ᢙN䬛ห⒟ᐲ䬽䭱䯃䮇䯃䬶䬑䭛䬿໧㗴䬹䬞⸘▚น⢻䬶䬨䫻੹ᓟ䬾䫺䭺䭩䬑䬮䭙䬽⸘▚☸ᐲ䬛䬸䭢䬸䭢ዊ䬤䬞䬹䭚䬷⠨䬗䭘䭛䭚䬽䬶䫺ᬌ⸛䬺୯䬨䭚ᤨᦼ䬺䬹䭙䬳䬳䬑䭚䬽䬶䬾䬹䬓䬚䬷ᕁ䬓䭍䬨䫻䯵22ಽᨆ䫺⠨ኤz ૶䭞䭛䬹䬓ℂ↱䭡䭍䬷䭐䭚䬷¾ ⒖ᬀᕈ䰊㕙ୟ䰊䭱䯃䮗䯃䮟䮊䮐䰊䮨䮱䮊䮏ਇ᣿¾ 䬮䬯䬦䫺⧯䬓ੱ䬽෻ᔕ䬾㆑䬕䬚䭑䬦䭛䬹䬓z ಽᨆ¾ Ꮏቇ䬾⚊ᦼ䬛㊀ⷐЈ᭽䫽䬹䮞䮰䮊䮏䮜䭰䯃䮧䭇䬽ኻᔕЈ⒖ᬀᕈ
¾ 䭱䯃䮗䯃䮟䮊䮐ᄢ䫺㕙ୟЈ৻ℂ䬑䭚• 䮰䭫䮝䮰䮱䭡૶䬕೨䬺䭱䯃䮞䮺䮄䯃䮀䫺Ꮢ⽼䭺䯃䮐䫺න⚐䬹䭩䮲䭻䮱䮁䮧ൻ¾ ૶䭞䬹䬠䭛䬿䬹䭘䬹䬓ℂ↱䯂႐㕙䬺ਲ䬦䬓Јಽ㊁․᦭?

• 䮥䮏䮱䭶䮀䭡䬜䬲䬰䭙⸃䬚䬹䬞䬵䭑⦟䬓
¾ 䮨䮱䮊䮏䬛ਇ᣿⏕䫺⺑᣿ਇ⿷¾ 䭩䮲䭻䮱䮁䮧䭗䭙B/FᲧ䬽਎⇇䫺䮠䭶䮏䮲ᯏ䬽ሽ࿷
¾ CFDኾ↪䬽䮰䭫䮝䮰䮱?
23ಽᨆ䫺⠨ኤ䯴2䯵z 䬦䬚䬦䫺䬬䭛䬶⦟䬓䬽䬚¾ ୘ੱ⊛䬺䬾…䮰䭫䮝䮰䮱䬾㊀ⷐ䫺੹ᓟ䬽HPC䬺䬾ᔅ㗇• ਗ೉ಣℂ䯴䱧䲄䱙䱮䲈䫺䱶䲉䲅䱨䱧䫺䱴䱜䱶䲈䲂䱙䱮䲈䯵 䫺䮨䮩䮱㓏ጀ䬽ⶄ㔀ൻ 䎃 䎃 䎃 䎃 䎃 䎃䎃 䎃 䎃 䎃 Ј䎃 ᕈ⢻䭡಴䬨䬽䬛㔍䬦䬞䬹䭚• 䮥䮲䮈䭷䮱䮊䮐䫺AMR䫺䭶䮱䮴䮜෻ᓳ䬷䬓䬲䬮䯴ⶄ㔀䯵ቯ⇟ಣℂ䬽Ⴧട¾ ⥶ⓨቝቮ䬾䫺એ೨䬾HPC䬽‧ᒁᓎ䬯䬲䬮䬛੹䬾?

• ᚑഞ૕㛎䬽⟂䫺䭫䮕䮠䯃䭾䮮䮺䬛ᔅⷐ• ⧯ᚻ䬛⢒䬲䬵䬓䬹䬓z ੹ᓟ䬽ኻᔕ䬷䬚¾ ૶䭞䭛䬹䬓ℂ↱䬽⺞ᩏ䯂ಽᨆ䭡⛮⛯¾ 䭩䮞䮴䯃䮈䭡ᄌ䬗䬵䭎䭚• ᅤ૗䬺ᄢᄌ䬚/ᄢᄌ䬶䬹䬓䬚䭡ታ㓙䬺䭓䬲䬵䭎䭚• ᚑഞ䯂ᄬᢌ੐଀䬽⫾Ⓧ䫺⚻㛎⠪䬽Ⓧ• ໪⫥ᵴേ

24䯴ෳ䎃 ⠨䯵

25

JAXA䭩䮞䮱䯴CFD䯵䬽㘃ဳൻ
Capability♽᭴ㅧൻ䮞䮴䭷䮰䮦䮺䭷♽

Capacity♽

ቇⴚ೑↪

Ꮏቇ೑↪

New Faceቇⴚ䯸Ꮏቇ․ᱶ

ᩰሶὐᴺ

☸ሶᴺ

䮲䯃䮞ਗ೉

㗔ၞಽഀਗ೉

㗔ၞಽഀਗ೉

න৻㗔ၞ

ⶄᢙ㗔ၞ

ⶄᢙ㗔ၞ

᭴ㅧᩰሶ

᭴ㅧ/㕖᭴ㅧ/⋥੤ᩰሶ

F77, F90

F77, F90

F90, C, C++䮤䭫䮺䮆䯸᭴ㅧ૕೑↪

ⷙೣ⊛ᄙ㊂ㅢା

ਇⷙೣዋ㊂ㅢା

ਇⷙೣዋ㊂ㅢା

䯴ォ⟎䯸AlltoAll, 䮗䮱䭩䯵 䯴⴫㕙䬽䭎䯸ᄢၞㅢା䯵 䯴⴫㕙䬽䭎䯸ᄢၞㅢା䯵

XPF, 䯴MPI䯵

MPI

MPI

F90
䮴䯃䮐䭫䮺䮗䮰䮺䮀
MPI⸘▚䬷ㅢା䬽ಽ㔌
ᄢਗ೉

ਛ䱊ዊਗ೉

ਛ䱊ዊਗ೉

ᄢਗ೉

ዊ㗫ᐲᄙ㊂I/Oᄙ㗫ᐲዋ㊂I/Oᄙ㗫ᐲዋ㊂I/Oᄙ㗫ᐲዋ㊂I/O

ዋᢙ䭿䮮䮝

ᄙᢙ䭿䮮䮝

ᄙᢙ䭿䮮䮝

ዋᢙ䭿䮮䮝
26

CFD䬽ᢙ୯⸃ᴺz ၮ␆ᑼ䰆䮑䮚䭯䮀䮏䯃䭶䮀z ⸃ᴺ䯴ౖဳ䯵䎃
p = n䰆㓁⸃ᴺ, 䎃 p=n+1䰆㒶⸃ᴺ䎃ㅪ┙৻ᰴ䋨⇹ⴕ೉䋩

MG䇮AMR
27

CFD䬽ᢙ୯⸃ᴺ䯴2䯵
z ⸃ᴺ䯴䬳䬴䬜䯵¾ 㕖࿶❗ᕈ䫺MHD Ј䮤䭩䮄䮺Ј䭶䮱䮴䮜෻ᓳ䫺FFT䫺AMG¾ ቟ቯᕈ໧㗴
Ј࿕᦭୯
¾ ☸ሶ♽
Јੂᢙ

z 䬬䬽ઁ䬺૶䭞䭛䬵䬓䭚䭑䬽¾ METIS,Zoltan,SCOTCH 㕖᭴ㅧᩰሶ↪㗔ၞಽഀ䮋䯃䮲¾ OpenFOAM
䭱䯃䮞䮺䮄䯃䮀CFD䮋䯃䮲䯴C++䯵
¾ Paraview ,Visitนⷞൻ䮋䯃䮲
¾ CHEMKINൻቇ෻ᔕ䮍䯃䮝䮲

28

4.7 並列擬似乱数発生アルゴリズムとペタスケール時代のモンテカルロ法の展望

並列擬似乱数発生アルゴリズムとペタスケール時代のモンテカルロ法の展望
発表者：国立情報学研究所

特任教授 三浦謙一

ＨＰＣシステムのコア数は増大する傾向にあり、京コンピュータで６５Ｋコア、ローレ
ンスリバモア研究所の Sequoia(IBM BlueGene/Q)ではすでに 1 Ｍコアを超えている。このように並列度が高くなるとこれまでの数値解法アルゴリズムの中には十分にスケールせず性能が飽和あるいは低下するものもあり、数値計算アルゴリズムそのものの再評価が必要となってきている。そこでモンテカルロ法に期待が寄せられる。モンテカルロ法の利点
としては① Embarrasingly Parallel でコア間の通信量が少なく、スケーラビリティが良い。
② コアあたりの所要メモリー量が少ない。
③ 統計的な計算のため、ある程度の Fault Resilience が期待できる。
したがってミリオンコア時代のアルゴリズムとしてモンテカルロ法をこれまで以上にいろいろな応用分野に適用することを考えてみる価値がある。実際コンピュータの黎明期には線形代数、楕円型偏微分方程式の境界値問題（ディリクレ問題）等に対してモンテカルロ法の適用が研究されたこともあり、後者についてはその後も研究が続けられており、新しいアルゴリズムが提案されている。
モンテカルロ法による計算を数 10 万から 100 万規模のコア数に対して行うためには質の良い擬似乱数の発生アルゴリズムが不可欠である。具体的な要件としては以下が挙げられる。
① 周期が十分長い。
② 統計的な質が良く、TESTU01 のような検定プログラムをパスする。
③ 擬似乱数の発生が高速である。
④ 多数の CPU コアに対して擬似乱数の初期値を高速に計算し配布するアルゴリズムが存在する。
（これが特に重要である。
）
本報告者が提唱する、Multiple Recursive Algorithm(MRG)に基いた発生法である
MRG8 はこれらの要件をすべて満足するものであり、今後の展開が期待される。
なお本アルゴリズムの実装（逐次版）は SS 研ユーティにも登録されている。

SS◊䝨䝍䝇䜿䞊䝹ᩘ್ィ⟬䝷䜲䝤䝷䝸WG

୪ิᨃఝ஘ᩘⓎ⏕䜰䝹䝂䝸䝈䝮䛸䝨䝍䝇䜿䞊䝹᫬௦䛾䝰䞁䝔䜹䝹䝻ἲ䛾ᒎᮃ
2012ᖺ7᭶23᪥ᅜ❧᝟ሗᏛ◊✲ᡤ≉௵ᩍᤵ㻌୕ᾆ㻌 ㅬ୍

䝰䞁䝔䜹䝹䝻䟿

䜰䜴䝖䝷䜲䞁
䞉䝇䜿䞊䝷䝤䝹䜰䝹䝂䝸䝈䝮䛸䛧䛶䛾䝰䞁䝔䜹䝹䝻ἲ䞉஘ᩘ䛻ồ䜑䜙䜜䜛ᛶ㉁䞉୪ิ໬䛾ᐜ᫆ᛶ
䞉䝰䞁䝔䜹䝹䝻ἲ䛾฼⏝ศ㔝䛾ᣑ኱䞉 Open Petascale Libraries Project

Million Core᫬௦䛾୪ิ䜰䝹䝂䝸䝈䝮䡚䝰䞁䝔䜹䝹䝻ἲ䛾䝯䝸䝑䝖䡚-㻌 ᩘከ䛟䛾ఏ⤫ⓗ䛺ᩘ್ィ⟬䜰䝹䝂䝸䝈䝮䛜䚸୪ิᛶ䛾ほⅬ䛛䜙㻌 㻌 ぢ┤䛧䛜ᚲせ䚹-㻌 䝰䞁䝔䜹䝹䝻ἲ䛿CPU䝁䜰㛫䛾㏻ಙ䛜ᑡ䛺䛔䚹㻌㻌ĺ“embarrassingly parallel” 䛷䇾 highly scalable䇿䛷䛒䜛䚹㻌 㻌- CPU䝁䜰ᙜ䛯䜚䛾ᡤせ䝯䝰䝸㔞䛜ᑡ䛺䛟䛶῭䜐䚹-㻌 ⤫ィⓗฎ⌮䛾䛯䜑䚸䛒䜛⛬ᗘ”Fault Resilient” 䛷䛒䜛䚹䝁䞁䝢䝳䞊䝍䛾ึᮇ䠄1950ᖺ௦䠅䛻䛿䛔䜝䛔䜝⪃䛘䜙䜜䛯䛜䚸㻌 䛭䛾ᚋ䛒䜎䜚ὀ┠䛥䜜䛶䛔䛺䛛䛳䛯ศ㔝䛷䛾ᛂ⏝䛜ᮇᚅ䛥䜜䜛䚹- ೫ᚤศ᪉⛬ᘧ䠄ቃ⏺್ၥ㢟䠅䛾ゎἲ䕧u =0 (Laplace Equation)
䕧u =f (Poisson’s Equation)
䕧u + u=0  - 㐃❧᪉⛬ᘧ䛾ゎἲ- ᅛ᭷್ၥ㢟䛾ゎἲ

with u=g on the boundary(Dirichlet ⬒䟛㧰ↅ.)

4ᨃఝ஘ᩘ䜰䝹䝂䝸䝈䝮䛸䛧䛶ᚲせ䛺᮲௳- ࿘ᮇ䛜㛗䛔㻌 ྛ䝁䜰ෆ䡔CPUᩘ䛾⊂❧䛺䝅䞊䜿䞁䝇-Ⰻ䛔⤫ィⓗ䛺ᛶ㉁䠄㏲ḟ䛸୪ิ䠅㻌 㻌 㻌 TestU01䛾䜘䛖䛺᳨ᐃ䝥䝻䜾䝷䝮䛜᭷⏝- 㧗㏿䛻Ⓨ⏕䛷䛝䜛㻌 㻌 㻌 㝖⟬䚸ศᒱ䜢ྵ䜎䛺䛔- ㉸୪ิ䛻㐺䛧䛯ึᮇ್ศ㓄䜰䝹䝂䝸䝈䝮䜢ᣢ䛴㻌 㻌 㻌 Jump-ahead䛾ᐜ᫆ᛶ
5Random Number Generation Algorithms(1) Linear Congruential Method: X n = (a X n-1 + c) mod M
Multiplicative 䠖 c=0.  Period㻌 = 2j-2
Mixed 䠖㻌 㻌 㻌 c .ne. 0. 㻌 Period㻌 = 2jwhere M=2j 䠄usually Machine Word Size)(2) Binary M-sequence with Primitive Trinomial:
X n = (Xn-m .eor. Xn-k ) mod 2㻌 㻌 Period = 2k-1 (m<k)(3) Generalized Fibonacci Method with Primitive Trinomial:
X n = (Xn-m op. Xn-k ) mod M㻌 㻌 Period =(2k-1)2j-1~ 2 k+j-1
㻌 㻌 㻌 㻌 where op. is {+, -, *}.
(4) Generalized Recurrence Method With Large Prime Modulus(Multiple Recursive Generator or MRG)㻌 㻌 㻌 㻌 Period =pk-1~ 2 j*k㻌 ,where p ~ 2 j is a Prime Number (e.g., 2 31 – 1)

Characteristics (and Advantages) of MRG
with 8th-order Full Primitive Polynomials
From here on, consider p=231 – 1 and k = 8.
f(x) = (x8- a1x7 – a2x6 – a3x5 – a4x4 – a5x3 – a6x2 – a7x – a8) mod ( p )(1) Better chance of finding a generator which possesses good Lattice Structurewith non-constrained full coefficients㻌 㻌 㻌 dt > 1/sqrt(1 + a12 + a22 + a32 + a42 + a52 + a62 + a72 + a82 )㻌(Ref. P.L’Ecuyer , INFORMS Journal on Computing, Vol.9, No.1,pp.57-60, 1997)(2) Reasonably long period: (231-1)8 – 1 ~ 4.5*1074(3) Ease of applying to vector/parallel processing, due to small
dimension of the states(4) Many primitive polynomials to choose from:I(pk – 1)/k/(pk – 1) = 2.2%
(5) Simple and straight-forward implementation is possible ( and hopefully,
good performance too!)

A Sample Polynomial with
Good Lattice Structure (LatMRG)
Ref: Prof. Pierre L’Ecuyer, Univ. Montreal
xn = a1xn-1 + a2xn-2 + a3xn-3 + a4xn-4 + a5xn-5 + a6xn-6 + a7xn-7 + a8xn-8 mod( p )where㻌 a1 = 1089656042
a5 = 189748160
a2 = 1906537547
a6 = 1984088114
a7 = 626062218
a3 = 1764115693
a4 = 1304127872
a8 = 1927846343
㻌 p = 231 – 1=2147483647
Figure of Merit: M32=.62729

Out of 345597 tested polynomials with a8 as primitive roots,
31843 are primitive polynomials䠄9.2%䠅.

Other Reported Full Polynomials
with Good Lattice Structure
3rd Order:

xn = a1xn-1 + a2xn-2 + a3xn-3 mod( p ), where p=231 -1 =2147483647

Grube-Dieter
㻌 a1 = 518175991
a2 = 510332243
a3 = 71324449
4th Order:

Dietera1 = 388425559
a2 = 227651891
a3 = 5412951

L’Ecuyera1 = 2021422057
a2 = 1826992351
a3 = 1977753457

xn = a1xn-1 + a2xn-2 + a3xn-3 + a4xn-4 mod( p ), where p=231 -1 =2147483647

L’Ecuyer㻌 a1 = 2001982722
a2 = 1412284257
a3 = 1155380217
a4 = 1668339922
Ref. http://crypto.mat.sbg.ac.at/results/karl/server/node7.html

As of 8/10/06

Performance Measurement of MRG8 Rev. 06/14/11
on Machines with Different Architectures
SystemArchitecture
㻌㻌㻌㻌

Clock

Peak
Rate of
Performance
Generation
䠄GH䡖䠅㻌 㻌 㻌 㻌 (Gflops)(106 dp rng/sec)

Fujitsu VPP5000
Vector䠄Proprietary)㻌 㻌 㻌 .3
.5
NEC SX-6
Vector䠄Proprietary)
Fujitsu PrimePower Scalar (Sparc64V)
1.3
㻌 HPC2500
1.3
IBM p-series 690
Scalar (Power4)
SGI Altix3700 㻌 㻌 㻌 㻌 Scalar (Itanium 2) 㻌
1.3
AMD Opteron
Scalar
2.0
Intel Woodcrest Xeon Scalar (1 core)
Fujitsu Primergy㻌 㻌 㻌 Scalar䠄Xeon EM64T)
RX200S2
Fujitsu PrimeQuest480 Scalar (Itanium 2)
Fujitsu PrimeQuest Scalar (Itanium2)
㻌 RXI300
RIKEN/Fujitsu “K”
Scalar (SPARC64VIIIfx)(1 core)9.6 㻌 㻌
8.0
5.2

201.8
224.3
23.2

5.2
5.2
4.0

11.4
39.7
34.0

2.66
3.6

10.6
7.2

108.1
79.9

1.5
1.6

6.0
6.4

80.3
83.5

2.0

16.0

34.5
(2 Integer Ops./clock)Order of Recurrence vs Rate of Generation
MRG 2

Platform:Fujitsu PrimePower HPC2500

MRG 4Rate of Generation (Million rng/sec)

30

MRG 8
20

MRG 16
10

MRG 32
MRG 64

10100㻌 㻌 10200 㻌 㻌 10300 㻌

10400

10500

10600

Period

1D and 2D Uniformity Tests of MRG8
File Name(95% Chisq.)
RNG5a
RNG5b
RNG5c
RNG51
RNG55
RNG56
RNG57
RNG58
RNG59

n
106
106
106
106
106
106
106
106
106

Mean S.D. 1D(100) 1D(1000) 2D(10^2).50000 .28868 (123.3) (1073.6 )(123.3 ).50007 .28871 83.82
1008.0799.89.49977 .28858 96.65
1024.52111.40.49953 .28883 114.35
1042.9895.19.50015 .28869 83.99
963.66106.71.49988 .28851 79.81
1074.5983.61.49929 .28885 117.17
1013.80
101.32.50033 .28865 95.40
1021.6797.12.50010 .28876 99.61
979.15102.52.49989 .28842 84.101
999.03119.95

2D(32^2)(1098.6)
1026.831033.65
943.461027.241032.411022.111001.07965.891044.58

A Sample Autocorrelation Function of
MRG8.002139 @ 11
0.002

1.96/sqrt(N)

0.001

20

40

60

80

100

W(Lag)-0.001- 1.96/sqrt(N)-0.002-.002521 @ 34-.002457 @ 67

N=106, File: RNG59

Parallelization 䠄1䠅~㻌 Random Initialization䡚
Same as “Birthday Problem” (e.g., Feller)• Period of Random Numbers䠖N
Usage of Random Numbers in each process/thread䠖n
Total number of chunks which should not overlap: N/n = m
Total Number of Cores䠖p• Probability of no collision = 1䞉(1-1/m)䞉(1-2/m)䞉䞉䞉䞉䠄1-(p-1)/m)~= exp[ -(p-1)p/(2m)] ~= 1 – (p-1)p/(2 m)
Probability that at least one collision occurs㻌 ~=㻌 䠄p-1)p/(2m)Example䠖N=10^74, n=10^15, m=10^59, p=10^6 ~ 5䞉䞉10^(-48)㻌

Parallelization 䠄2䠅
~ Jump-Ahead (First Order Recurrence) ~Xn+1 = a Xn mod(M)

X16
a16
a8

X0

a2

a
X1

X8

a4
X4

X2

Parallelization of First Order Recurrence
by Recursive Doublingx䡅 = ai xi-1x䡅 = ai ai-1 …… a3 a2 a1 x0

*
*
*
a1

a2

*
*

a3

*

a4

*

*

*

*

*
a5

a6

*
a7

a8

Parallelization 䠄2䠅~㻌 Jump-Ahead (the 8th Order Recurrence) ~xn = a1xn-1 + a2xn-2 + a3xn-3 + a4xn-4 + a5xn-5 + a6xn-6 + a7xn-7 + a8xn-8 mod( p )Define a Transfer Matrix A:A =a 1 a 2 a 3 a 4 a 5 a6 a 7 a 8
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0
0 0 0 0 0 0 1 0

Then

xn-1
xn-2
xn-3
x = xn-4
xn-5
xn-6
xn-7
xn-8 ,

,

xn
xn-1
xn-2
x’ = xn-3
xn-4
xn-5
xn-6
xn-7

x’ = A X Mod(p).

Jumping Ahead the Sequence of MRG
X23X’㻌 = A X mod(M)A: Transfer Matrix

X19
X’’’’
A16

X’

X
X0 X

1

X2

X4

X7

X8

X17
X16

Jump-Ahead for Arbitrary DistanceIn order to compute xn = An x0 mod(p) for an arbitrary n:j2(1) Compute and store Aj = A mod(p) (j=0,1,2,3,4,……).
(2) Represent “n” in the binary form,e.g., (bm-1,…..,b2,b1,b0).
(3) Multiply Aj mod(p)’s together only when bj=1
(j=0,….,m-1), to obtain An.
Note: The same strategy works with the polynomial formulation with fewerarithmetic operations (Knuth).

Random Number Testing on Quantum
Diffusion Monte Carlo Application

Ref. K.Hongo, R.Maezono and K.Miura, J. Comp. Chem. 31 ,pp.2186-2194 2010

Architecture/Software Issues• Language- Fortran-C
Parallel Programming Model- shared memory model (Open MP)- Distributed memory model (MPI)- Hybrid model• Portability• Multicore CPUs
• Fast Integer Arithmetic (64 bits mandatory)䝰䞁䝔䜹䝹䝻ἲ䛾෌ㄆ㆑䝁䞁䝢䝳䞊䝍䛾ึᮇ䠄1950ᖺ௦䠅䛻䛿䛔䜝䛔䜝⪃䛘䜙䜜䛯䛜䚸䛭䛾ᚋ䛒䜎䜚ὀ┠䛥䜜䛶䛔䛺䛛䛳䛯ศ㔝䛷䛾ᛂ⏝䛜ᮇᚅ䛥䜜䜛䚹౛䛘䜀㻌 㻌 -೫ᚤศ᪉⛬ᘧ䠄ቃ⏺್ၥ㢟䠅䛾ゎἲ䕧u =0 (Laplace Equation)
䕧u =f (Poisson’s Equation)
䕧u + u=0   with u=g on the boundary (Dirichlet B.C.) -㐃❧᪉⛬ᘧ䛾ゎἲ㻌 㻌 -ᅛ᭷್ၥ㢟䛾ゎἲ౛䠖䝷䝥䝷䝇᪉⛬ᘧ䛾ゎἲ㟁☢ሙ䚸ຊሙ䛾ィ⟬䛻୙ྍḞ'< 0 in *<  g on J䠍䠊ᕪศἲ/᭷㝈せ⣲ἲ䠎䠊ቃ⏺せ⣲ἲ䠏䠊䝰䞁䝔䜹䝹䝻ἲ- Walk on Spheres Method- Walk on Rectangles Method- Walk on Boundary Methodゎἲ䛾ẚ㍑䞉㻌 ᭷㝈せ⣲ἲ䛩䜉䛶䛾᱁ᏊⅬ䛷䛾್䜢ィ⟬䛩䜛㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 ᚲせ䛒䜚• ቃ⏺せ⣲ἲቃ⏺᮲௳䛜ኚ䜟䜛䛤䛸䛻ቃ⏺᪉⛬ᘧ㻌㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 䜢ồ䜑┤䛩ᚲせ䛒䜚• 䝰䞁䝔䜹䝹䝻ἲ㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 -㻌 ᚲせ䛺᱁ᏊⅬ䛷䛾ィ⟬䛰䛡䛷Ⰻ䛔㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 -㻌 ቃ⏺䛾್䛜᫬㛫ⓗ䛻ኚ䜟䛳䛶䜒㻌㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 㻌 ༶᫬ồゎ䛜ྍ⬟㟁ሙෆ䛷䛾Ⲵ㟁⢏Ꮚ䝅䝭䝳䝺䞊䝅䝵䞁䛾䜘䛖䛻䚸⢏Ꮚ䛜Ꮡᅾ䛩䜛ሙᡤ䛾್䛾䜏䛜ᚲせ䛺ሙྜ䛻᭷ຠ䛛䚹౛䠖䝷䝥䝷䝇᪉⛬ᘧ䛾ゎἲྛ᱁ᏊⅬ䛷䛾䝷䞁䝎䝮䜴䜷䞊䜽䜢⾜䛖ᚲせ䛿↓䛔䚹    Method0'< 0 in *<  g on Jbg(x,y)

䐟୍␒㏆䛟䛾ቃ⏺䛻᥋䛩䜛䜘䛖䛺㻌 ෇䠄䠏ḟඖ䛺䜙⌫䠅䜢ᥥ䛟䚹䐠෇࿘ୖ䛷䝷䞁䝎䝮䛻Ⅼ䜢䛸䜛䚹䐡䛭䛾Ⅼ䜢୰ᚰ䛻᭦䛻෇䜢ᥥ䛟䚹䐢Ⅼ䛜ቃ⏺䛛䜙H௨ෆ䛻䛺䛳䛯䜙㻌 㻌 ቃ⏺Ⅼ䛷䛾䡃䠄x,y)䜢⣼⟬䛩䜛䚹

*

1

-a

0

ィ⟬౛
䠍0,000N=䠍

1P
No. of walks = 10000 Sum of hits = 5098
Solution at (x,y) =0.5098
Last radius of sphere = 0.000013661
Generated random numbers = 121292

N=40,000
No. of walks = 40000 Sum of hits = 20106
Solution at (x,y) =0.50265
Last radius of sphere = 0.0000106012
Generated random numbers = 484677

N=1,000,000
No. of walks = 1000000 Sum of hits = 499734
Solution at (x,y) =0.499734
Last radius of sphere = 1.71 x 10^-6
Generated random numbers = 12083215

ὀ䠖ཝᐦゎ䛿0.50000䛷䛒䜛䚹

J

-b
a

WOS䛾⪃ᐹ䛸ㄢ㢟䞉㻌 㻌㻩㻌㻝㻜㻙㻠䛾ሙྜ෇䛜ቃ⏺䛻฿㐩䛩䜛䜎䛷ᖹᆒ㻌㻌 㻝㻞ᅇ䛾㼃㼍㼘㼗䛜ᚲせ㻌 䠄⌮ㄽୖ䛿㻌㼘㼛㼓㻔 㻕㻌䛻ẚ౛䠅㻌㻌• ௚䛾ィ⟬ἲ䛸䛾ᐃ㔞ⓗ䛺ẚ㍑䚸≉䛻ቃ⏺せ⣲ἲ䛸䛾⢭ᗘ㻌 㼢㼟㻌ィ⟬᫬㛫䛾ẚ㍑䛜ᚲせ㻌㻌• ௚䛾㻼㻰㻱䜈䛾ᣑᙇ㻌䜎䛸䜑• 䝬䝹䝏䝨䝍䝇䜿䞊䝹᫬௦䛻ධ䜚䚸䝁䜰ᩘ䜒
O(10^6)䛻✺ධ㻌 㻌 䜶䜽䝃䛷䛿䛹䛖䛺䜛䠛• 䛔䛟䛴䛛䛾ィ⟬䜰䝹䝂䝸䝈䝮䛾ぢ┤䛧䛜ᚲせ• 䝰䞁䝔䜹䝹䝻ἲ䛿䜎䛩䜎䛩᭷ᮃ䛻㻌 㻌 㻌 㻌 㐺⏝ศ㔝䛾ᣑ኱䠄౛䠖PDE)• 䛭䛾䛯䜑䛻䛿Ⰻ䛔୪ิ䠄ᨃఝ䠅஘ᩘ䛜ᚲせ㻌 㻌 㻌 㻌 MRG䛾ᥦ᱌• ᪤Ꮡ䛾䝷䜲䝤䝷䝸䛻䛒䜛஘ᩘ䝹䞊䝔䜱䞁䛿せὀព• Open Petascale Libraries䠄OPL) Project䛷䜒㻌 㻌 ஘ᩘ䛿㔜せ䛺఩⨨௜䛡䛷䛒䜛

Acknowledgement
The author would like to thank• Prof. P.L’Ecuyer : Sample Primitive Polynomials• Dr. C.Chen at Fujitsu Computer System:Measurement on VPP5000, Primergy, PrimeQuest• Dr. M.Kurokawa, RIKEN:Measurement on SX-6, Woodcrest• Mr. Hayakawa, AMD Japan:Measurement on AMD Opteron• Mr. Hotta, Fujitsu Limited (Makuhari):Measurement on Sparc64VIIIfx

4.8 FFTE について

FFTE について
筑波大学システム情報系（計算科学研究センター）
高橋大介
高速フーリエ変換（Fast Fourier Transform，以下 FFT）は，科学技術計算において今日広く用いられているアルゴリズムである．
本発表では FFT を計算するためのライブラリである FFTE ライブラリの概要および二次
元分割を用いた並列三次元 FFT アルゴリズムについて述べている．
FFTE は 2002 年 5 月にバージョン 1.0 を公開し，最新のバージョンは 2011 年 11 月に公開したバージョン 5.0 である．シングルプロセッサ，共有メモリ型並列（OpenMP），分散メモリ型並列（MPI），ハイブリッドプログラミングモデル（OpenMP+MPI）をそれぞれサポートしており，http://www.ffte.jp/よりソースコードを入手可能である．
FFTE ライブラリの設計方針としては，キャッシュブロッキングを行うことでキャッシュ
ミスをできるだけ少なくすると共に，ノード内およびノード間のスケーラビリティの確保を目指している．
2012 年 6 月の Top500 において，2 システムが 10PFlops の大台を突破しているが，今後出現すると予想される 100PFlops 級マシンでは，ほぼ確実に 100 万コアを超える規模のものになると予想される．
この場合，MPI と OpenMP のハイブリッド実行を行ったとしても，
MPI プロセス数が 10 万個以上になる可能性がある．
これまで，並列三次元 FFT においては三次元（x，y，z 方向）のうちの一次元（例えばz 方向）のみを分割して配列を格納することが多かったが，MPI プロセスが 1 万個の場合，
z 方向のデータ数が 1 万個以上でなければならず，三次元 FFT の問題サイズに制約があるという問題があった．この問題を解決するために，FFTE の並列三次元 FFT では一次元分割に加えて y，z 方向の二次元で分割する，二次元分割もサポートしている．
二次元分割では一次元分割に比べて，全体の通信量は 2 倍となるが，MPI プロセス数が大きく，かつ全対全通信のレイテンシが大きい場合には，二次元分割の方が一次元分割に比べて通信時間が短くなる場合がある．
T2K 筑波システムの 4096 コア（256 ノード）までを用いた性能評価の結果，並列三次元
FFT において，二次元分割により通信時間を削減することで，MPI プロセス数が多い場合に一次元分割に比べて性能が改善されることを示した．
二次元分割を用いた並列三次元 FFT は，MPI プロセス数が 1 万個を超えるようなペタス
ケール環境において，より効果的であると考えられる．

FFTEについて
高橋大介
筑波大学システム情報系
計算科学研究センター

2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

1

発表内容
• FFTEライブラリの概要• 二次元分割を用いた並列三次元FFT
アルゴリズム2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

2

FFTEライブラリの概要• FFTEは高速フーリエ変換（Fast Fourier Transform，
以下FFT）を計算するためのライブラリである．
• 2002年5月にバージョン1.0を公開し，現在のバージョ
ンは2011年11月に公開したバージョン5.0．
• 逐次，共有メモリ（OpenMP），分散メモリ（MPI），ハイブリッド（OpenMP+MPI）をそれぞれサポート．
– 二次元および三次元の実数FFT– 一次元，二次元および三次元の複素数FFT– それぞれ，2，3，5のべき乗からなるデータ点数に対応している．
– 三次元の並列複素数FFTでは二次元分割にも対応．

• http://www.ffte.jp/ よりソースコードを入手可能．
SS研 ペタスケール数値計算ライブラリWG
第4回会合

2012/7/23

3

FFTEライブラリの設計方針• 性能– キャッシュミスをできるだけ少なくする．
• 明示的なキャッシュブロッキングを行っている．

– スケーラビリティの確保（ノード内，ノード間）• API– 引数を少なくして使いやすさを重視．
– SGI SCSLやAMD ACMLに類似している.

• ポータビリティ– ノード間通信：MPI– 演算部分：Fortran + OpenMP
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

4

背景
• 2012年6月のTop500において，2システムが10PFlopsの大台を突破している．
– Sequoia: 20.133 PFlops (1,572,864 Cores)– K computer: 10.510 PFlops (705,024 Cores)• 今後出現すると予想される，100PFlops級マシンは，ほぼ確実に100万コアを超える規模のものになると予想される．
– MPIとOpenMPのハイブリッド実行を行ったとしても，
MPIプロセス数が10万個以上になる可能性がある．
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

5

目的
• 並列三次元FFTにおける典型的な配列の分散方法– 三次元（x，y，z方向）のうちの一次元のみ（例えばz方向）のみを分割して配列を格納．
– MPIプロセスが1万個の場合，z方向のデータ数が1万点
以上でなければならず，三次元FFTの問題サイズに制約．

• x，y，z方向に三次元分割する方法が提案されている[Eleftheriou et al. ’05, Fang et al. ’07]．
– 各方向のFFTを行う都度，全対全通信が必要．

• 二次元分割を行うことで全対全通信の回数を減らしつ
つ，比較的少ないデータ数でも高いスケーラビリティを得ることを目的とする．
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

6

z方向に一次元ブロック分割した
場合の並列三次元FFT
1. x方向FFT

2. y方向FFT

z

z

z
x

y

3. z方向FFT

x

y

y

x

各プロセッサでslab形状に分割
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

7

三次元FFTの超並列化• 並列アプリケーションプログラムのいくつかにおいては，三次元FFTが律速になっている．
• x，y，zのうちz方向のみに一次元分割した場合，超並列化は不可能．
– 1,024×1,024×1,024点FFTを2,048プロセスで分割できない（1,024プロセスまでは分割可能）• y，zの二次元分割で対応する．
– 1,024×1,024×1,024点FFTが1,048,576（=1,024×1,024）プロセスまで分割可能になる．
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

8

y，z方向に二次元ブロック分割した場合の並列三次元FFT
1. x方向FFT

2. y方向FFT

z

z

z

x

y

3. z方向FFT

x

y

y

x

各プロセッサで直方体形状に分割
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

9

二次元分割による並列三次元FFT• 二次元分割した場合，P × Q 個のプロセッサにおいて，– P 個のプロセッサ間で全対全通信を Q 組– Q 個のプロセッサ間で全対全通信を P 組行う必要がある．
• MPI_Comm_Split()を用いてMPI_COMM_WORLDを
y方向（ P プロセッサ）とz方向（ Q プロセッサ）でコミュニ
ケータを分割した．
– 各コミュニケータ内でMPI_Alltoall()を行う．

• 入力データがy，z方向に，出力データはx，y方向に二次元ブロック分割されている．
– 全対全通信はy方向で1回，z方向で1回の合計2回で済む．
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

10

一次元分割の場合の通信時間• 全データ数を N ，プロセッサ数を P × Q ，プロセッサ間
通信性能を W （Byte/s），通信レイテンシを L （sec）とする．
• 各プロセッサは N /(PQ) 2 個の倍精度複素数データを
自分以外の PQ − 1 個のプロセッサに送ることになる．
• 一次元分割の場合の通信時間は⎛16 N ⎞⎟⎟
T1dim = ( PQ − 1)⎜⎜ L +2( PQ) ⋅ W ⎠⎝16 N≈ PQ ⋅ L +（sec）
PQ ⋅ W
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

11

二次元分割の場合の通信時間• y方向の P 個のプロセッサ間で全対全通信を Q 組行う．
– y方向の各プロセッサは N /( P 2Q ) 個の倍精度複素数データを，y方向の P − 1 個のプロセッサに送る．
• z方向の Q 個のプロセッサ間で全対全通信を P 組行う．
– z方向の各プロセッサは N /( PQ 2 )個の倍精度複素数データを，ｚ方向の Q − 1 個のプロセッサに送る．
• 二次元分割の場合の通信時間は⎛⎛16 N ⎞
16 N ⎞⎟⎟ + (Q − 1)⎜⎜ L +⎟⎟T2 dim = ( P − 1)⎜⎜ L + 2
2
P Q ⋅W ⎠PQ ⋅W ⎠⎝⎝≈ ( P + Q) ⋅ L +2012/7/23

32 N
（sec）
PQ ⋅ W

SS研 ペタスケール数値計算ライブラリWG
第4回会合

12

一次元分割と二次元分割の場合の通信時間の比較（1/2）• 一次元分割の通信時間

16 N
T1dim ≈ PQ ⋅ L +PQ ⋅ W• 二次元分割の通信時間

32 N
T2 dim ≈ ( P + Q) ⋅ L +PQ ⋅W• 二つの式を比較すると，全プロセッサ数 P × Q が大きく，
かつレイテンシ L が大きい場合には，二次元分割の方が通信時間が短くなることが分かる．
SS研 ペタスケール数値計算ライブラリWG
第4回会合

2012/7/23

13

一次元分割と二次元分割の場合の通信時間の比較（2/2）• 二次元分割の通信時間が一次元分割の通信時間よりも少なくなる条件を求める．

32 N
16 N
( P + Q) ⋅ L +< PQ ⋅ L +PQ ⋅ W
PQ ⋅ W

を解くと，

N<( LW ⋅ PQ)( PQ − P − Q )16• 例えば， L = 10 −5 (sec)， W = 10 9 (Byte/s)，P = Q = 64を上の式に代入すると， N < 1010 の範囲では二次元分割の通信時間が一次元分割に比べて少なくなる．
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

14

性能評価
• 性能評価にあたっては，二次元分割を行った並列三次
元FFTと，一次元分割を行った並列三次元FFTの性能比較を行った．
• Strong Scalingとして N = 323 , 643 , 1283 , 2563 点の順方向
FFTを1～4,096MPIプロセスで連続10回実行し，その平均の経過時間を測定した．
• 評価環境
– T2K筑波システムの256ノード（4,096コア）を使用– flat MPI（1core当たり1MPIプロセス）– MPIライブラリ：MVAPICH 1.2.0– Intel Fortran Compiler 10.1
– コンパイルオプション：”ifort –O3 –xO”（SSE3ベクトル命令）
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

15

T2K筑波システムの概要• T2K筑波システムの諸元– ノード台数：648（Appro Xtreme-X3 Server）– ノード構成：4ソケット／ノード– CPU：Quad-Core AMD Opteron 8356
（Barcelona 2.3GHz）– 総メモリ容量：20TB– ネットワークインターフェース：DDR InfiniBand
Mellanox ConnectX HCA×4– ネットワークトポロジ：Fat Tree
– Full-bisectionバンド幅：5.18TB/s
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

16

二次元分割を用いた
並列三次元FFTの性能

1000

N=32^3
100
GFlopsN=64^3
10
N=128^3
1

N=256^3

96
40

24

10

6
25

64

16

4

1

0.1

Number of cores2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

17

考察（1/2）• N = 32 3 点FFTでは良好なスケーラビリティが得られていない．
• これは問題サイズが小さい（データサイズ：1MB）ことから，全対全通信が全実行時間のほとんどを占めているからであると考えられる．
• それに対して，N = 256 3 点FFT（データサイズ：512MB）
では4,096コアまで性能が向上していることが分かる．
– 4,096コアにおける性能は約401.3 GFlops（理論ピーク性能の約1.1%）– 全対全通信を除いたカーネル部分の性能は約10.07 TFlops（理論ピーク性能の約26.7%）
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

18

1000

256^3点FFTにおける一次元分割と二次元分割の性能比較

一次元
分割
二次元
分割

GFlops

100
10
1

96
40

24
10

6
25

64

16

4

1

0.1

Number of cores2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

19

並列三次元FFTの実行時間の内訳(256cores, 256^3点FFT)
0.06

Time (sec)

0.05

演算時間
通信時間

0.04
0.03
0.02
0.01
0
一次元分割

2012/7/23

二次元分割

SS研 ペタスケール数値計算ライブラリWG
第4回会合

20

考察（2/2）• 64コア以下の場合には，通信量の少ない一次元分割
が二次元分割よりも性能が高くなっている．
• 128コア以上では通信時間を少なくできる二次元分割
が一次元分割よりも性能が高くなっていることが分かる．
• 二次元分割を行った場合でも，4,096コアにおいては96%以上が通信時間に費やされている．
– 全対全通信において各プロセッサが一度に送る通信量がわずか1KBとなるため，通信時間においてレイテンシが支配的になるためであると考えられる．

• 全対全通信にMPI_Alltoall関数を使わずに，より低レベ
ルな通信関数を用いて，レイテンシを削減する工夫が必要．
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

21

二次元分割の場合のFFTの限界性能（1/2）• 以下のようなシステムを仮定し，N = 10243 点FFTを計算– プロセッサ数： P × Q = 256 × 256 = 65,536– プロセッサ当たりの理論ピーク演算性能：「∞」 Flops– プロセッサ間通信性能： W = 「∞」 Byte/s– 通信レイテンシ: L = 1μsec• N 点FFTの演算量は 5 N log 2 N とする．
• N点のデータ（倍精度複素数）の通信量： 16 N バイト• 全対全通信に要する時間は，
32 N= 5.12 ×10 − 4 sec
T2 dim ≈ ( P + Q) ⋅ L +PQ ⋅W
3• N = 1024 点FFTの限界性能は，(5 N log 2 N ) / T2 dim ≈ 314.57 TFlops
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

22

二次元分割の場合のFFTの限界性能（2/2）• 全対全通信にshuffle-exchangeアルゴリズムを使った場合の通信時間は，⎛⎛
16 N ⎞
16 N ⎞⎟⎟ + (log 2 Q)⎜⎜ L +⎟⎟T2 dim ≈ (log 2 P)⎜⎜ L + 2
2
P Q ⋅W ⎠PQ ⋅ W ⎠⎝⎝= 1.6 ×10 −5 sec• N = 10243 点FFTの限界性能は，(5 N log 2 N ) / T2 dim ≈ 10.07 PFlops
2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

23

まとめ
• FFTEライブラリの概要および二次元分割を用いた並列三次元FFTアルゴリズムについて述べた．
• 並列三次元FFTにおいて，二次元分割により通信時間を削減することで，MPIプロセス数が多い場合に性能が改善されることを示した．
• 二次元分割を用いた並列三次元FFTは，プロセッサ数
が1万個を超えるようなペタスケール環境において，
より効果的であると考えられる．

2012/7/23

SS研 ペタスケール数値計算ライブラリWG
第4回会合

24

4.9 フラグメント分子軌道（FMO）法と並列 FMO プログラム OpenFMO について

フラグメント分子軌道（FMO）法と並列 FMO プログラム OpenFMO について
九州大学情報基盤研究開発センター

稲富 雄一

分子軌道（Molecular Orbital, MO）法は、原子や分子の中にある電子の振る舞い（分布）を求めるための量子力学に基づいた計算手法である。MO 法では、分子の安定構造や化学反応における反応エネルギーや活性化エネルギーなどの有用な情報が得られるため、理論化学の研究者だけでなく、実験化学分野の研究者にも幅広く利用されている。しかしながら、経験的なパラメタを一切用いない非経験的な MO 法では、最も計算量が小さ
い Hartree-Fock（HF）法と呼ばれる近似計算手法を用いた場合でも、O(N4)の計算量であり、MO 法の応用が期待される分野の一つであるタンパクや DNA などの巨大生体分子に対する計算を行うことができない。
北浦和夫現神戸大学特命教授によって開発されたフラグメント分子軌道（FMO）法は、従来困難であった大規模分子に対する MO 計算を行うために開発された計算手法である。FMO 法では、計算対象となる巨大な分子を 20～40 原子程度の小さなフラグメントに分割して、各フラグメント（モノマー）、および、フラグメントペア（ダイマー）に対する小規模な MO 計算を行うことで、巨大分子全体の状態を近似する。複数のモノマ
ーやダイマーに対する MO 計算を並列に処理できること、さらに、各小規模 MO 計算
自身も並列処理が可能であることなどから、大規模並列処理向きの計算手法である。そこで我々は、京コンピュータをはじめとした超並列計算機を用いて、効率よく並列 FMO計算を行うための最適化を並列 FMO プログラム OpenFMO に対して行った。OpenFMOは九州大学と九州先端科学技術研究所でスクラッチから開発された C 言語で記述されている並列 FMO プログラムであり、HF 法に基づいた FMO 計算に特化しているためソースコードが 5 万行程度と小規模であるため、最適化が行い易い。
我々がこれまでに行った主な最適化は、MPI/OpenMP ハイブリッド並列化、グローバル
カウンタを用いた動的負荷分散による負荷均等化、および、FMO 計算で用いる中間データの保存方法の改良、の３つである。これらの改良を行った OpenFMO プログラムを用いて２万並列までのベンチマークを行った結果、非常に高い並列性能を示すことが分かった。また、２万並列時に１万６千原子を超える分子に対する FMO 計算が 30 分弱で終わることが分かった。このことから、より大規模な並列計算を行うことで、大規模分子に対する電子状態計算を 10 分未満で行うことが近い将来に可能であり、創薬等の分野での応用が、より簡単になる可能性が見えてきた。
一方で、小規模 MO 計算部分においては、99％以上の計算量である分子積分計算は、並列化が容易なために、並列処理により計算時間が大幅に短縮されたが、その代わりに、一般化固有値問題求値（対角化）の時間がノード間並列化が困難な部分として、並列化
効率を低下させる要因になっていることが分かった。今後は、対角化などのノード間並列が困難な部分の処理時間を隠蔽するなどの工夫をすることで、更なる並列処理速度の改善をする必要がある。

フラグメント分子軌道（FMO）法と並列
FMOプログラムOpenFMOについて
九州大学情報基盤研究開発センター

稲富 雄一
2013年2月22日

1

目次
1.
2.
3.
4.

目的
量子化学計算フラグメント分子軌道（FMO）法OpenFMOの最適化と性能評価
FMO計算における行列計算

5. まとめ2

目標
大規模生体分子に対する第一原理（量子力学）に基づいた電子
状態計算を高速に行うことができるようにしたい！
大規模分子軌道計算の応用が期待される分野¾創薬•薬の候補分子と生体分子（タンパク質、核酸、糖鎖、ウイルスなど）との結び付きやすさが分かるため、生理活性の推定や候補物質の絞り込み（スクリーニング）に使える
右図：H5N1型鳥インフルエンザウイルスのノイラミニダーゼ（PDB ID=2HU4, 23,856原子）
感染した細胞からウイルスを分離させるために必要。このたんぱくの働きが阻害されると、細胞内で合成されたウイルスが他の細胞に移動できなくなるので、増殖できない（抗インフルエンザ薬のタミフル、リレンザはノイラミニダー
ゼ阻害剤）

計算機シミュレーションによるスクリーニングを行う
ためには、多数回の計算を短時間で行う必要がある！

3

ついでにもうひとつ（リボゾーム）
60Sユニット（PDB ID=1S1I）

40Sユニット（PDB ID=1S1H）

タンパク質合成（ペプチド結合形成）を担う

m-RNAの暗号解読を担う¾細菌とヒトではリボゾームの構造が違うので、細菌のリボゾームの機能を選択的に阻害する薬剤（リボゾーム阻害剤）は細菌性感染症の薬になる（ストレプトマイシン、テトラサイクリンなどの抗生物質はリボゾーム阻害剤）¾タンパクの折りたたみの機構の解明につながる（かも）4

主に、Hartree-Fock（HF）法について

量子化学計算（≒分子軌道法）

5

量子化学計算とは？
量子力学に基づいた第一原理計算
分子中の電子に対するSchrödinger方程式を解いて、分子内での電子の運動（分布）を求める
Ĥ elec Ψ = Eelec ΨnZA ⎞ n n
1⎛ 1 ⎞ n N ⎛Hˆ elec = ∑ ⎜ − ∇i2 ⎟ + ∑∑ ⎜⎜ −⎟⎟ + ∑∑
2 ⎠ i A ⎝ ri − R A ⎠ i j <i ri − r ji ⎝n = 電子数
N = 原子数
Z A = 原子Aの原子番号得られた波動関数 Ψ が電子状態を表している

一般には３粒子以上が関係するので、近似的に解く必要がある6

分子軌道（Molecular Orbital, MO）法とは？
¾ 電子のSchrödinger方程式の近似解法の１つ¾ 電子が分子軌道（MO）と呼ばれる軌道（のようなもの＝1電子
波動関数）を運動している、という化学者に分かりやすい描像を与える¾ 実験値を利用する（半）経験的方法から、実験値を用いない非経験的な方法、近似の粗い方法や高精度な方法など、たくさんの方法が存在量子化学計算≒分子軌道計算

7分子軌道法のなかま（これでも、まだ一部・・・）¾ 半経験的方法– （拡張）Huckel法、PPP法、CNDO法、PM5、AM3, DFTB, … etc.

¾ 非経験的（ab initio）分子軌道法– Hartree-Fock（HF）法• RHF法、UHF法、ROHF法– Post Hartree-Fock法• （Møller-Plesset,MP）法（摂動法）– MP2、MP3、MP4 …
• 配置間相互作用（Configuration Interaction, CI）法– CID、CISD、CISDT、full CI …
• クラスター展開（Coupled Cluster, CC）法– CCSD, CCSD(T), CCSDT、…
高精度（計算量大）• 多配置参照（Multi Reference, MR）法– MCSCF、MR-MP、 MR-CI、MR-CC …
非経験的分子軌道法は、ほぼすべて、HF法の解を0次近似として用いている

8

MO法で分かること （１） 分子のエネルギー
化学反応の起こりやすさ、結合の強さ、吸収波長などが分かる
遷移状態

正反応
⎯⎯⎯
→[ AB]A + B ←⎯⎯
⎯
逆反応

エネルギー

正反応の活性化エネルギー
逆反応の活性化エネルギー

A＋B
（始状態）
結合エネルギー（反応エネルギー）

[AB]
（終状態）

反応座標

9

MO法で分かること （２） 分子軌道
化学反応が起こる場所（活性部位）が分かる
Lysozyme分子のフロンティア軌道 (1,961原子、 10,969関数)
赤-青＝Highest Occupied Molecular Orbital (HOMO)
黄色-水色＝Lowest Unoccupied Molecular Orbital (LUMO)

10

Hartree-Fock法で解くべき基本方程式
Schrödinger方程式からHartree-Fock-Roothaan方程式まで

Ĥ elec Ψ = Eelec Ψエネルギー期待値を最小化する

n N
n
n
ZA
1 n
1
Hˆ elec = − ∑ ∇i2 − ∑∑
+ ∑∑
2 i
i
A ri − R A
i j <i ri − r j

Eelec =

Ψ Hˆ elec ΨΨ|ΨHartree-Fock近似
制限付きスピン軌道

χ 2 k (x) = ψ k (r )α (σ )k = 1, 2,...
χ 2 k +1 (x) = ψ k (r ) β (σ )% HF =Ψ≈Ψ

χ1 (x1 ) χ 2 (x1 ) L χ n (x1 )1 χ1 (x 2 ) χ 2 (x 2 ) L χ n (x 2 )
M

n!

M

O

M

χ1 (x n ) χ 2 (x n ) L χ n (x n )基底関数展開 ψ i (r) =NAO∑cμϕ (r; n , R )
μ=1

i

μ

μ

μ

Hartree-Fock-Roothaan方程式

FC = SCε11

Hartree-Fock-Roothaan（HFR）方程式
Schrödinger方程式に各種近似を適用して得られる代数方程式

FC = SCεF = Fock行列（対称行列）
S = 重なり行列（正定値対称行列）{}
ε = MOエネルギー（固有値） ≡ diag {ε , ε ,K , ε }C = MO係数行列（固有ベクトル） ≡ c1 , c 2 ,K , c NAO
1

2

N AO

エネルギー期待値を極小化する（スピン）軌道を求める問題が、分子軌道を基底関数で展開する線形結合係数（MO係数）を求める問題に置き換えられているNAO

ψ i (r) = ∑ cμi ϕμ (r; nμ , Rμ )
μ=1

12

HFR方程式の特徴

FC = SCε

Fock行列F = Hcore + G

１電子ハミルトン行列

Hcore = T + V

運動エネルギー積分
核引力積分
２電子ハミルトン行列⎛ 1 ⎞Tμν = ∫ dr1 ϕ μ (r1 ) ⎜ − ∇12 ⎟ ϕν (r1 )⎝ 2 ⎠⎛ N −Z A ⎞Vμν = ∫ dr1 ϕ μ (r1 ) ⎜⎜ ∑⎟⎟ ϕν (r1 )⎝ A r1 − R A ⎠
n

n

σ

λ

Gμν = ∑∑ Dσλ {2 ( μν | σλ ) − ( μλ | σν )}２電子積分( μν | σλ ) = ∫ dr1 ∫ dr2 ϕμ (r1 )ϕν (r1 )密度行列

Dμν = ∑ Cμi Cν i

1
ϕσ (r2 )ϕλ (r2 )r1 − r2n /2i•Fock行列Fと重なり行列Sを入力とした一般化固有値問題である•

与えるべきFock行列Fの計算には、この方程式の解であるMO係数行列Cが必要であるため、この方程式は非線型方程式であり、反復解法（Self
Consistent Field, SCF法）によって解かなければならない•

計算量はO(N4)（２電子積分の計算量に依存）

13

HFR方程式を解く手順
重なり積分S, １電子ハミルトン行列
Hcore、および、２電子積分を計算
重なり積分をCholesky分解（S=UTU）
密度行列の初期値（D0）を設定
密度行列と２電子積分から２電子
ハミルトン行列Gを計算

n

n

σ

λ

Gμν = ∑∑ Dσλ {2 ( μν | σλ ) − ( μλ | σν )}( μν | σλ ) = ∫ dr1 ∫ dr2 ϕμ (r1 )ϕν (r1 )F=Hcore+Gを計算

密度行列を更新

FC=SCεを解く
Cから新たな密度行列Dを計算

no

密度行列が収束した？

n /2

Dμν = ∑ Cμi Cν iiyes
収束したF, C, Dなどを用いて、望みの期待値（エネルギーなど）を計算

1
ϕσ (r2 )ϕλ (r2 )r1 − r2

14

そのままでは大規模計算ができない・・・
計算量がO(N4)かそれ以上である– コンピュータが1000倍速くなっても、入力データの大きさは約5.6倍（=1000の4乗根）にしかできない

参考：SC|05, T. Ikegami et al.,technical paper
20,581原子77,754電子
164,442関数
1,398フラグメント

従来分子軌道法(Hartree-Fock法）で計算した場合の予想計算時間＝およそ２０万年！
FMOでの計算時間＝72.5時間（Opteron (model 246), 600プロセッサ使用時）

15

大規模分子に対する電子状態計算を行うための計算手法

フラグメント分子軌道（FMO）法16

フラグメント分子軌道（FMO）法とは？
• 大規模生体分子に対するMO計算を行うために北浦和夫神戸大学特命教授により開発された並列処理向きの計算手法• 分子を20～40原子程度のフラグメントに分割して、各フラグメ
ント（モノマー）とフラグメントペア（ダイマー）に対する電子状態
で、分子全体の電子状態を近似• フラグメント間ではクーロン相互作用のみを考慮• フラグメント間距離に応じて、いくつかの近似を適用することで計算時間を大幅に削減

17

FMO計算のフローチャート
初期化
Self-Consistent Charge loop
as Many-task Computation
モノマー電子状態計算

No

モノマー電子状態計算

・・・

モノマー電子状態計算

収束？
Fragment-Pair Correction
as Many-task ComputationYes
ダイマー電子状態計算

ダイマー電子状態計算

・・・

ダイマー電子状態計算

全エネルギー計算

Schematic Flow of FMO

FMO計算のほとんどはモノマーやダイマーに対する小規模電子状態計算（フラグメント計算）

18

並列FMOプログラムの構造
初期化

Self-Consistent Charge loop
as Many-task Computation

master process
モノマー計算Job assignment,collection of results

worker processes

No

group 2group Ngrank 0rank 0rank 0・・・

ダイマー計算全エネルギー計算

・・・

モノマー計算Dimer Correction
as Many-task Computation
ダイマー計算・・・

ダイマー計算Schematic Flow of FMO

rank 1
・・・

・・・

・・・

rank 1

Converged?
Yesgroup 1

rank 1

モノマー計算data transfer of monomer density matrices¾ いくつかのgroupに分割して、各groupに小規模電子状態計算を割り当てる（粗粒度並列処理）¾ group内プロセスで小規模電子状態計算を並列処理（細粒度並列処理）する¾ Master-Worker型（ただしワーカー（グループ）間通信が存在する・・・）

19

超並列化に向けた最適化

OpenFMOの最適化と性能評価

20

OpenFMOとは？
• 九州大学、九州先端科学技術研究所で開発された並列FMO
プログラム• Hartree-Fock法に基づいた計算に特化しており、コードがコン
パクト（C言語、約5万行）
さまざまな最適化の適用が容易！

• 中間データ（フラグメント密度行列データ）の保存方法を工夫することで、大規模入力にも対応• 並列化における標準API（MPIやOpenMP）を用いる• フリーのツール（MPICH2, OpenMPI, GCCなど）でも使えるよ
うにする• ターゲットはSMPクラスタ型並列計算機

21

これまでに行ってきた最適化¾ MPI/OpenMPハイブリッド並列化¾ 中間データの分散保存、アクセス方法の改良¾ グローバルカウンタを用いた動的負荷分散による負荷均等化22

MPI/OpenMP hybrid並列化「京」のノード構成¾SMPクラスタ型並列計算機を構成する各ノードには、プロセッサ（コア）数よりも少ないNICしか搭載されていないケースが多いので、Flat MPIプログラム
では、通信資源の競合による性能低下を引き起こす恐れがあるSPARC64TMVIIIfx
core

core

core

core

core

core

core

core

通信プロセッサ¾MPIライブラリが使用するメモリ量や通信量の観点から、起動MPIプロセス数が少ない方が望ましい

通信エ
ンジン

通信エ
ンジン

通信エ
ンジン

通信エ
ンジン

OpenMPによるノード内スレッド並列とMPIによるノード間プロセス並列を組み合わせたhybrid並列化を適用した
ただし、分子積分部分にのみ手を加えており、行列計算に関連する部分は、ライブラリ任せ23

モノマー密度行列の保存・アクセス方法
数万～10万フラグメントの超巨大分子に対するFMO計算を行う場合を想定して、2種類のモノマー密度行列の保存・アクセス方法を検討した
方法１：全体で１（複数）セットのデータを持つ
group 1
rank 0rank 1group 2
rank 0rank 1group 3rank 0rank 1group 4rank 0rank 1

方法２：データ保存専用プロセスの利用
storage nodes
rank 0rank 1

rank 2

rank 2

rank 2

rank 2

rank 2

rank 3rank 3rank 3rank 3rank 3

参照処理

worker nodes

データの
やり取り

group 1group 2group 3rank 0rank 0rank 0rank 1

rank 1

rank 1

rank 2

rank 2

rank 2

rank 3rank 3rank 3

ターゲットプロセスからGet+計算グループ内
Bcast

参照処理

更新処理

ターゲットプロセスへのPut

更新処理

ターゲットプロセスへのSend

同期処理

記憶グループ間の対応するプロセス間でのBcast（または、Allreduce）

同期処理

記憶グループ間の対応するプロセス間でのBcast（または、Allreduce）• （MPI-2の）片側通信が必要となる• データ保存のための特別なプロセスが必要ない（コードがシンプル）

ターゲットプロセスからのRecv＋計算グループ
内Bcast• MPI-1規格でも記述できる• 記憶グループ用と計算グループ用の少なくとも2種類のコードが必要（コードが複雑）24

アクセス方法による性能の違い
モノマー電子状態計算部分のプロファイル
入力＝アクアポリン（PDB ID=2F2B，492フラグメント、基底関数＝STO-3G）
計算機＝PowerEdge M610(quad core Xeon×2)×5, interconnect=10GbE, MPICH2-1.3.1, intel C++ compiler (12.0.0)

中間データ保存専用プロセスを利用した場合

MPI-2の片側通信を用いた場合•MPI-2の片側通信の実装によっては、密度行列データへのアクセス性能が著しく低下•データ保存専用プロセス用いた場合には密度行列データへの効率的なアクセスができた

25

方法２の大規模実行時の通信性能
データ保存プロセス数、計算プロセス（グループ）数とモノマー密度行列データ取得時間
測定環境: Riken Integrated Cluster of Clusters (RICC), 63 groups, 32 processes/group, 4threads/process (128threads/group), OpenMPI(version 1.4.3), Intel C++ compiler (version 11.0)
入力データ:アデノウィルスKNOBドメイン(PDB ID=1nob, 576 fragments), 基底関数 6-31G*

データ保存プロセス

ケース１

rank 0rank 1

rank 0rank 1 ・・・

group 1

rank 0rank 1 ・・・

group 2

rank 0rank 1 ・・・

group 63 rank 0

rank 1 ・・・

group 1

rank 0rank 1 ・・・

group 2

rank 0rank 1 ・・・

group 63 rank 0

rank 1 ・・・

・
・
・

ケース３

rank 0group 1

・
・
・

ケース２

rank 0

計算プロセス

データ取得時間*

総実行時間

ケース１

37.8 (37.8)

10773

ケース２

44.2 (0.70)

177.6

ケース３

37.3 (0.59)

177.8* モノマー密度行列データを取得するために要した時間の全グループでの合計（かっこ内はグループあたり）

大規模並列時でも効率よく密度行列
データにアクセスできる

26

粗粒度並列部分の並列化効率
Machine:FUJITSU PRIMERGY(CX400) in Kyushu University(16 processes x 8 threads/group (=128 threads/group),# of groups=10-160) Compiler: FUJISTU MPI + FUJITSU C compiler（version 1.2.0)
Input : Adenovirus KNOB domain (PDB ID=1nob, 576 fragments, 16,764 atoms), basis set= 6-31G* (142,974 functions)16
linear scale

speedup ratio

8

30分で
計算終了

total FMO calculation

4

2

11280 (10)2560 (20)5120 (40)10240 (80) 20480 (160)# of worker threads (# of groups)

２万並列時でも高い並列化効率を保っている！

27

小規模電子状態計算部分の模擬コード
calculation_SCF_energy(int ifrag) {for (id=0; id<(Nifc4c+1); id++ ) make_cutoff_table(id);// カットオフテーブル計算
calc_twoint ( myrank, nprocs, ERI_buffer);// 2電子積分
V=0.0;for ( i=0; i<Nifc4c; i++ ) V+=calc_ifc4c( myrank, nprocs, i); // 4中心クーロン積分
for ( i=0; i<Nifc3c; i++ ) V+=calc_ifc3c( myrank, nprocs, i); // 3中心クーロン積分
for ( i=0; i<Nfrag; i++ ) V+=calc_ifc2c( myrank, nprocs, i); // 2中心クーロン積分
if ( myrank == 0 ) {calc_oneint( S, H );// 1電子積分
calc_projection( P );// 射影演算子}
MPI_Allreduce( V, comm );SCF_procedure( S, H, P, V, ERI_bufferi, Difrag ,comm );// SCF計算}
計算内容

1回あたりの計算時間

小規模電子状態1回あたりの呼び出し回数

カットオフテーブル計算

0.01～0.2秒

10～30

2電子積分

数10～数100秒

1*

4中心クーロン積分

数10～数100秒

10～30

3中心クーロン積分

数秒

10～20

1電子積分、射影演算子

0.05～0.5秒

1

2中心クーロン積分

0.05～0.5秒

フラグメント数（～10万）*; 大量の2電子積分（106～109個）が、すべてがメモリに保存できた場合

計算時間の多くを占める各種分子積分計算の負荷均等化が並列性能向上の上で重要！

28グローバルカウンタを用いた負荷均等化の効果
小規模電子状態計算実行プロファイル
入力＝アデノウィルスのKNOBドメイン（PDB ID=1NOB，576フラグメント、基底関数＝STO-3G）
計算機＝PowerEdge M610(quad core Xeon×2)×5, interconnect=10GbE, MPICH2-1.3.1, intel C++ compiler (12.0.0)

グローバルカウンタによる負荷均等化なし

グローバルカウンタによる負荷均等化あり

分子積分部分に負荷不均衡が見られる

分子積分計算の負荷均等化が図られている

29

フラグメント計算部分の速度向上比

速度向上比

Input＝Adenovirus KNOB domain（PDB ID=1NOB, 576 monomers, basis set＝6-31G*）
Machine＝FUJITSU PRIERGY CX400 (8-core Xeon×2/node, 1476 nodes, Interconnect=InfiniBand FDR)
Compiler=Fujisu C compiler, MPI=Fujitsu MPI library, SSL II

16

linear scaling

14

モノマー積分12

モノマー全体
ダイマー積分10

ダイマー全体

8
6
4
2
0
0

50

100

150

200

250

ワーカースレッド数

分子積分計算部分だけなら、256並列時で約97%の高い並列化効率を達成しているが、全体では並列化効率が悪い・・・（特に、ダイマー計算）

30

ダイマー計算における各部分の計算時間
Input＝Adenovirus KNOB domain（PDB ID=1NOB, 576 monomers, basis set＝6-31G*） Machine＝FUJITSUPRIERGY CX400 (8-core Xeon×2/node, 1476 nodes, Interconnect=InfiniBand FDR)
Compiler=Fujisu C compiler, MPI=Fujitsu MPI library, SSL II, 8 threads/process

160,000

計算時間（秒）

分子積分
対角化

120,000

初期化
その他

80,00040,000

0
0

50

100

150

200

250

ワーカースレッド数• 分子積分計算時間はスレッド数の増加とともに減少• 対角化部分の計算時間は、ノード間並列していないため、ワーカー
スレッド数に依存せず一定

31

SCF計算で用いている行列演算関数
初期密度行列（D0）生成重なり行列（S）、１電子ハミルトン行列（H）の計算
重なり行列をCholesky分解

操作内容

S = Ut U

BLAS, LAPACK関数名
dpotrf２電子積分を計算して、バッファに保存

バッファに保存された２電子積分と密度行列を用いてFock行列を計算
5~10収束加速のためのFock行列の更新

密度行列
更新

Fdiis = ∑ ci Fii

一般化固有値問題を解いてMO係数行列（C）を求める
新たな密度行列、エネルギーの計算
no

エネルギーが収束した？

yes
終了（次の処理）daxpy（ dsygst, dsysv）

FC = SCε⎛ F ' = U − t FU −1 ⎞ dsygst⎜⎟ dsytrd, dorgtr, dsteqr⎜ F ' C ' = C 'ε ⎟⎜ C = U −1C ' ⎟ dtrsm⎝⎠
E = 2tr {D(H + F)} ddot
N elec /2
678
D = C nC (n = diag(1,1,...,1, 0,...0))t

32

フラグメント計算部分で扱う行列のサイズ
700
600

モノマー

500

ダイマー

頻度

400
300
200
100
0

行列サイズ

一般化固有値問題の特徴• 密行列である• 1/2～1/3の固有値と対応する固有ベクトルが必要

33

まとめ¾ FMO計算は、解くべき問題を分割して解くことで、これまで困難であった大規模生体分子に対する電子状態計算を可能にした¾ FMO計算の粗粒度並列部分の並列化効率は非常に高い¾ 細粒度並列部分（フラグメント計算の並列処理）では、計算量の最も大きな分子積分計算は高並列が可能であるが、非並列部分の存在などによる並列性能低下が大きい¾ FMO計算で扱う行列サイズは小さいため、行列計算部分の並列化による性能向上は困難かも
更なる並列性能向上のためには、対角化と他の処理（他のフラグメ
ントの積分計算など）をオーバーラップさせるなどの対策が必要！

34

謝辞
以下のみなさま（敬称略）の協力を得ています– 眞木淳（九州先端研）– 本田宏明（九州大学）– 高見利也（九州大学）– 小林泰三（九州大学）– 青柳睦（九州大学）
本研究の一部は理化学研究所と九州大学との共同研究で行っている。また、科学研究費補助金基盤研究（C）「超並列フラグメント分子軌道法プログラムライブラリの開発」（課題番号22550015）、および、JST,CRESTの研究領域「ポストペタス
ケール高性能計算に資するシステムソフトウェア技術の創出」の研究課題「省メモリ技術と動的最適化技術によるスケーラブル通信ライブラリの開発」の支援を受けている。さらに、学際大規模情報基盤共同利用・共同研究拠点の公募型共同研究平成２３年度採択課題「超並列フラグメント分子軌道法プログラム
OpenFMOの 性能評価と高性能化」による計算機利用を行っている。

35

ご静聴ありがとうございました