学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

jh190042-NAH

高性能・変動精度・高信頼性数値解析手法とその応用
中島研吾（東大情報基盤センター）
概要
本研究は，最先端のスパコン向けに開発された高性能数値アルゴリズムに対して，半精度から倍精度，倍々精度までの広範囲をカバーする変動精度演算を適用し，精度保証，そのための自動チューニング手法を開発する。開発された手法を様々なア
プリケーションに適用することで，低精度を中心とした変動精度演算の科学技術シ
ミュレーションへの有効性を検証する。開発したアルゴリズム，アプリケーションの消費電力の直接測定によって，各計算の特性と低精度演算の有効性を消費電力の観点から検討する。
 藤田航平（東大）高性能アルゴリズム・地震シ

1. 共同研究に関する情報(1) 共同研究を実施した拠点名

ミュレーション
 近藤正章（東大）精度保証・電力測定・AT 東 大 （ Oakforest-PACS ， Reedbush-U/H/L ， 奥田洋司（東大）工学シミュレーション・精度

Oakbridge-CX）
 東京工業大学（Tsubame-3）

保証(2) 共同研究分野

 森田直樹（東大）工学シミュレーション・精度■

超大規模数値計算系応用分野□

超大規模データ処理系応用分野

 荻田武史（東女大）精度保証・AT□

超大容量ネットワーク技術分野

 田中一成（早大）精度保証・AT□

超大規模情報システム関連研究分野

 尾崎克久（芝工大）精度保証・AT

保証(3) 参加研究者の役割分担

 片桐孝洋（名大）高性能アルゴリズム・精度保 中島研吾（東大）（代表）全体統括・高性能アル

証・AT

ゴリズム
 市村

 櫻井刀麻（名大）精度保証・AT

強（東大）（副代表）高性能アルゴリズム・

 八代

地震シミュレーション

尚（環境研）大気シミュレーション 河合直聡（理研⇒東大）高性能アルゴリズム・ 横田理央（東工大）（副代表）高性能アルゴリズ

量子科学シミュレーション

ム・精度保証 井上弘士（九大）精度保証・電力測定・自動チ 岩下武史（北大）高性能アルゴリズム

ューニング

 深谷

猛（北大）高性能アルゴリズム

 荒川隆（RIST）大気シミュレーション 塙

敏博（東大）性能最適化・電力測定 成瀬彰（NVIDIA）精度保証・AT

 伊田明弘（東大）高性能アルゴリズム

 堀越将司（インテル）精度保証・AT 星野哲也（東大）高性能アルゴリズム・性能最 大島聡史（名大）高性能アルゴリズム・性能最適化

適化

 坂本龍一（東大）電力測定・自動チューニング

 Gerhard（AT）

Wellein （ FAU

Erlangen-Nürnberg,

Germany)高性能アルゴリズム・量子科学シミュ 有間英志（東大）電力測定・AT

レーション・SELL-C- 古村孝志（東大）地震シミュレーション

 Achim Basermann（DLR, Germany）高性能アル1

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

演算精度を自動チューニング技術によって動

ゴリズム・量子科学シミュレーション的に制御する手法を確立する。

 Osni Marques （ Lawrence Berkeley National

 本研究によって開発された高性能・変動精度・Laboratory, USA）高性能アルゴリズム・AT

高信頼性数値解法を，自動チューニング機構を有するアプリケーション開発・実行環境

2. 研究の目的と意義

ppOpen-HPC （ JST-CREST

エクサスケールシステムにおける高性能数値ア2011-2018 ，

ルゴリズム実現には，メモリ・ネットワークの階https://github.com/Post-Peta-Crest/ppOpenHPC）及層の深化に対応した通信最適化（Serial，Parallel）

び（計算＋データ＋学習）融合を実現する革新とともに省電力・省エネルギー（以下「省電力」）

的ソフトウェア基盤 h3-Open-BDEC（科研費基に 向 け た 検 討 が 必 要 で あ る 。 Approximate

盤 S 2019-2023 ， http://nkl.cc.u-tokyo.ac.jp/h3-

Computing （ S. Mittal, A Survey of Techniques for

Open-BDEC/）に実装し，東大 Oakforest-PACS

Approximate Computing, ACM CSUR 48-4, 2016）は，（OFP），東大 Reedbush（RBH，RBL），東大

低精度演算の積極的活用により計算時間短縮，消Oakbridge-CX（OBCX），東工大 Tsubame-3（TSB3）

費電力削減を図る試みであり，従来は画像認識等で公開する。将来的には「富岳」も含む他のセの計算精度の要求されない分野を対象としていたンターのスパコンへの導入も視野に入れる。

が，昨今は数値計算において半精度から四倍精度まで演算精度を動的に変動させる変動精度

本研究は，最先端のスパコン向けに開発された高（Transprecision）の研究が進められている。数値計

性能数値アルゴリズムに対して，半精度から倍精算による近似解（数値解）は様々な計算誤差を含

度，倍々精度までの広範囲をカバーする変動精度

み，計算結果の信頼性の観点から，数値解の正し演算を適用し，精度保証，そのための自動チュー

さを数学的に保証する必要があり，低精度・変動

ニング手法を開発する試みとしては初めてのもの

精度使用時，悪条件問題には重要であるが，実問であり，開発された手法を様々なアプリケーショ

題で現れる大規模疎行列・H 行列への応用例はほンに適用することで，低精度を中心とした変動精

とんどない。本研究では，JHPCN システム群の中度演算の科学技術シミュレーションへの有効性をで消費電力当たり計算性能（GFLOPS/W 値）の高検証できる。開発したアルゴリズム，アプリケー

いシステムを主たるターゲットとして，以下を実ションの消費電力の直接測定によって，各計算の施する：特性と低精度演算の有効性を消費電力の観点から検討可能となる。

 疎行列演算，H 行列演算，ステンシル演算等の代表的数値アルゴリズム，各アプリケーション

3. 当拠点公募型研究として実施した意義（地震，大気科学，量子科学，構造力学）につJHPCN は多様な計算機環境を備え，実用的なシ

いて，Serial・Parallel 通信最適化に着目した高ステムとして最も GFLOPS/W 値の高い OFP，RBH，

性能最適化手法を各システムにおいて実装し，

RBL，OBCX，Tsubame-3 等の大規模システムを有低精度演算・変動精度演算について検討し，消し，本研究の目指す高性能・変動精度・高信頼性

費電力を測定する。

数値解法の研究には最適である。RBL，OBCX で

 疎行列演算や H 行列演算を対象として，特に悪は「ノード固定」における設定カスタマイズによ

条件問題における実用的な精度保証法を確立

り，個別ノードの消費電力測定が可能である。

する。更に①の各アルゴリズム，アプリケーシ

JHPCN は様々な分野の専門家を擁し，本研究のよョンについて所望の結果精度達成という条件

うな学際的研究を推進する体制を容易に構築でき，下で，計算時間や消費電力を最小化する最適な北大，東大，東工大，名大，九大各センターから2

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

様々な分野の研究者が参加している。JHPCN 各セの高精度化により，低精度・変動精度演算でもンターはオープンソースソフトウェア活用に積極安定した収束が得られる。

的であり，本研究の成果を公開，各センターのス

 従来の M 疎行列向け精度保証手法〔Ogita, T. et

パコンにデプロイし，講習会等の普及活動を協力

al., Computing, 2001〕は，相対誤差上限の見積もして行い，利用者拡大，ソフトウェアの更なる改りが厳しめである。

良が可能となる。
2018 年度研究では，図 1 に示すような不均質場

4. 前年度までに得られた研究成果の概要〔23〕における三次元熱伝導問題を有限体積法で離散化

本研究では，ステンシル計算（①構造格子，②して得られる連立一次方程式を ICCG 法で解くケ

半構造格子），疎行列演算（③一般行列，④悪条件ース（問題規模は 2,097,152=1283）を共通問題とし問題，⑤Adaptive CG（局所前処理による）），⑥Hて扱った。以下に，本アプリケーションを使用し行列，⑦精度保証，⑧消費電力測定，⑨自動チュ

た 2018 年度研究の研究成果のうち，1）性能最適

ーニング手法，の各項目についての研究開発を実化，（2）多様な計算精度，（3）精度保証，（4）実ア

施する。本研究は 2018 年度から，3 年計画としてプリケーションを対象とした消費電力測定に着目実施している。2018 年度は，上記①～⑩の項目にして紹介する。

  0@ z  zmax

ついて，（1）性能最適化，（2）多様な計算精度，（3）精度保証，（4）実アプリケーションを対象と

1した消費電力測定，を中心として予備的な検討をNZ実施し，後掲のように学会等で発表した他，ISC18，

2

SC18 等の国際会議の展示ブースでも紹介した。低1

精度・変動精度演算をアルゴリズム，最適化，精度保証，消費電力まで含めて扱った研究事例はほ      f  0

とんどなく，大きな反響があった。2018 年度に得図 1 不均質場における熱伝導問題

られた知見は以下の通りである：

1）性能最適化 従来，主として倍精度演算が適用されていた科
学技術計算に単精度・半精度演算を適用するこ

2.00E+00とによって，同等の精度の計算結果がより短時

CRS-ELL

1.885E+00

SELL-C-Sigma

1.50E+00

間で得られる場合がある。

9.360E-01

sec.

 低精度演算の適用により，演算密度が増加し時間当たり消費電力（Watt）が若干増加する場合

9.341E-01

1.00E+00

8.707E-01
7.483E-01
6.433E-01があるものの，①で述べた計算時間短縮に比例

5.00E-01

4.466E-01

して消費エネルギー（Joule）は減少する〕。
 悪条件問題では，低精度演算では正しい解が得0.00E+00

られない場合がある。特に半精度演算は変数の範囲が限定されるため注意が必要であり，精度

図 2 三次元熱伝導問題（図 1）の KNL での計算
時間（64 コア，128 スレッド）（ICCG 法）（-D：倍精度，-S：単精度，OPT：（星野他，2017）の最適化手法適用）のあまり要求されない反復法前処理等に適用するべきである。
 前処理付き反復法による疎行列演算では，反復改良法（Iterative Refinement），残差計算（r=b-Ax）
3

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

図 2 は図 1 に示す問題を対図 4 は，図 1 に示す問題を対象として，2 を変C化させた場合の Reedbush-U の 1 ノード（Intel Xeon

象として様々な行列格納手法に対して Oakforest-PACSs

Broadwell-EP（BDW），36 コア）を使用した場合の（ OFP ， Intel Xeon Phi

計算結果である。反復回数と ICCG 法の計算時間（Knights Landing，KNL））

について，倍精度演算の場合を 1 として，単精度

1 ノードを使用して実施し演算の値を示している。単精度演算の反復回数はた（64 コア，128 スレッド）若干増加するものの，計算時間は半分程度になっ場合の ICCG 法の計算時間ている。1/2 が増加し，係数行列の条件が悪化すである。CRS⇒ELL で 2 倍
以上，ELL⇒SELL-C-（図 3）

ると，単精度演算の反復回数は増加し，計算時間

図 3 SELL-C-形
式，C=2 and =8

も倍精度の場合を上回るようなケースがある。まで 15-25%程度の計算時間短た，1/2 が増加時は単精度演算による解の誤差も縮が得られている。更に倍精度⇒単精度で 30%程増加し，図 4 の1/2=106 のケースでは 10％にも達度の短縮が達成されている。C の最適値は倍精度している。倍精度⇒単精度⇒半精度によって扱う

で 8-16，単精度で 16 となっている。

ことのできる値の範囲，有効桁数は減少するため，低精度演算は，変数値の変動幅が大きい問題，大規模問題は不適当である。

2）多様な計算精度の適用
従来，科学技術計算には倍精度浮動小数点演算

図 5 は図 1 に示す問題の右辺，Z 方向メッシュ（FP64）を適用してきたが，2018 年度研究では，数などを変更した場合について，V100 上で単精度半精度（FP16）単精度（FP32），4 倍精度（FP128）

演算（前処理のみ半精度）の収束解が得られるパ等を上記②～⑥の各項目について適用し，様々なラメータの範囲を示している。倍精度，単精度（前CPU，GPU（Intel Xeon/Broadwell（BDW），Intel Xeon

処理含む）の場合は全て収束解が得られているが，Phi（KNL），NVIDIA Tesla P100（P100）
・V100（V100）） 前処理のみ半精度となった場合には解を得られるを使用して効果を検討した。単精度演算の適用に範囲が極端に少なくなっていることがわかる。

より，倍精度と比較して反復回数等の増加は見ら

4016
2048

れるものの，計算時間は 15～50%程度減少するこ

1024
512

とが確認された。計算時間の減少率はアルゴリズ

256

NZ

ム，アーキテクチャ，最適化の度合いによって異12864なっており，これらの関係について更なる検討が32必要であることが明らかとなった。

16

3.00

8

Iterations

sec.

4

2.50

2
1.E-04Ratio

1.E-02

1.E+00

1.E+02

1.E+04RHS

2.00

図 5 V100 上で単精度演算（前処理部分のみ半精度）を適用した場合の収束解が得られるパラメー
タの範囲（RHS：右辺，NZ：Z 方向メッシュ数，
1/2=100, 101, 102, 103）

1.50

1.00

0.50

0.00
1.E+003）精度保証1.E+01

1.E+02

1.E+03

1.E+04

1.E+05

1.E+062018 年度研究では，M 疎行列向けに提案されLambda-1/Lambda-2

た精度保証手法〔Ogita, T. et al., Computing,

図 4 単精度／倍精度演算時の反復回数比と1/2の関係（BDW）

2001〕を，図 1 の問題に適用した。求解と精度
4学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

ルギー（

保証を同時に実施するアルゴリズムを提案し，別々に計算する場合と比較して，

より

％程度計算

誤差が

の方が削減率が大きい。

の色数を増やしすぎると実行時間が長時間を削減できることが示された。この手法での場合，倍精度演算でも最大相対は， ）も大きく削減された。また，

くなり，エネルギー効率が悪化した。

のオーダーになってしまうことから，やや厳しすぎる判定であり，より精密に誤5. 今年度の研究成果の詳細

差の最大値を推定する手法の研究開発が必要で年度の重点実施項目は：ある（図 ）。
 反復改良法による低精度演算安定化（③）
 大規模分散並列問題で局所的・動的に演算精度を変動させる手法の検討（③，④）
 単精度より低い精度の演算の有効性の検討（⑤）
 疎行列・ 行列を係数行列とする実問題向けの精度保証手法の開発（⑦）
 様々な演算における消費電力の系統的測定（⑧）
 最適精度自動選択へ向けての検討（⑦，⑧，⑨）である。括弧内は， ．に示した①～⑨の研究開発

図  と最大相対誤差（実線）と相対残差ノルム（破線）の関係，図 の問題，倍精度項目である。予備的検討に留まった項目もあるが，当初予定の目標は順調に達成しており，研究成果，）消費電力・エネルギー測定

知見も多数得られている。以下①～⑨の各研究項

図目について順番に説明する。

の問題に対して，精度（倍精度・単精度），

データ配置方式（
法の色数（ ， ，，）等を変化させて測定した場

合の，計算時間と消費エネルギー（
図），①ステンシル演算：構造格子）の値を

変動精度を用いた効率的なステンシル計算の実に示す。

現のためには，同じステンシル計算の問題で，計算精度が異なる場合における，時空間タイリングの効果を明らかにする必要がある。そこで，今年度は，計算精度が異なるステンシル計算のプログ
ラムに，我々がこれまでに研究を行ってきた，異なるタイルレベルの並列性を持つ各種の時空間タ
イリングの手法を適用した，ベンチマークプログ
ラムの整備を実施した。ベンチマークプログラムの第一段階として，次元空間における熱伝導方

程式（拡散方程式）を差分法で離散化した場合に
図
倍精度
を単精度
消費エネルギーの比較化した場合の得られる， 点ステンシル型の計算を対象とし，倍精度浮動小数点（），単精度浮動小数点（

半精度の浮動小数点を想定した

低精度化することにより実行時間が大きく短縮（されることが確認でき，これに比例して消費エネ

），

整数型

），の三種類の演算精度（データ型）と，

時空間タイリングなし，および， 次元， 次元，
5

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

次元のタイルレベルの並列性を持つ時空間タイリングを適用した，計合わせて，合計期待できる。一方で，演算精度が異なる場合に対パターンの実装とを組みする，タイリングパラメータのチューニング方法

パターンの計算カーネルを実装が確立できていないので，今後の課題である。

し，問題サイズや時空間タイリングにおけるタイ

変動精度演算を想定した時空間タイリング付き

ルサイズ等のパラメータを自由に変化させて，統のステンシル計算のベンチマークプログラムの開一的に性能を評価することができるベンチマーク

発の第一段階は完了した。

プログラムを開発した。

示すように，本ベンチマークプログラムを利用す

節で提示した結果がることで，異なる演算精度間での性能評価を効率的に行うことが可能となるため，意義のある結果である。今後は，このプログラムを用いて，パラ
メータのチューニング方法の検討や，消費電力性
能を含む詳細な性能評価を行う予定である。また，
・年年度に対象としたが過去の研究で取り扱ったや，我々法等の，実ア

プリケーションを想定したベンチマークプログラ
ムの拡充も計画している。
②ステンシル演算：準構造格子のフルコード単精度実行を上での目指し，流体力学計算部分に引き続いて単精度化された諸物理過程の
を進めた。

による

最適化

には水の相変化と落下を解く雲

微物理過程，短波・長波放射による加熱・冷却を解く大気放射過程，格子スケールよりも小さい大図（ ノード）での結果（上：絶対性能，下：タイリングの効果）
東大の ノード（気運動による混合を解く乱流過程などの諸物理過程が導入されており，それぞれの過程について複数のスキームを実装している。このうち，最も頻コア）上で，

繁に利用する主要スキーム３つについて，先行し開発したベンチマークプログラムを用いて性能評価を実施した結果を図て単精度化が進められてきた。それぞれ，シングに示す。データサイズがルモーメント雲微物理スキーム

等しくなるように，各演算精度における問題サイズ （格子点数：

バンド大気放射スキーム）を設定し，性能（単位時間ーム

当たりに更新した格子点の数）と，各種時空間タ，乱流境界層スキである。これらのスキームの呼び出し部分の前後でドライバモジュール内の倍精度配列

イリング手法による性能向上の効果を検証した。

から単精度配列へのコピーを行い，膨大なソース

なお，タイルサイズ等のパラメータは，過去の研究を参考にして設定している。図，ブロード

コードの全体を変更することなく段階的な単精度から分かるよ化が進められた。

うに，演算精度に関わらず，時空間タイリングがの適用にはこの単 倍アダプタをそのまま利用し，デバイスとのコピー

効果的である。このことから，可能な範囲で低精イン，コピーアウトを実現した。対象とした物理度演算を利用し，同時に，適切な時空間タイリン過程の計算区間（コピー部分を除く）をそれぞれグ手法を組み合わせることで，大きな性能向上が，6，で測定した結学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

果を図に示す。３つのスキームはアルゴリズム

いては，メモリによって計算速度が律速されていの違いやコーディングの違いによって，それぞれる場合に高速化し，特にでは異なる結果の傾向を示した。

高速化した。一方で

（

）では単精

については，

倍での計算性能に劣化がみられた。これは主に度化による高速化がほとんどみられなかった。こ

巨大なループボディを持ち，多量のスカラー変数

れは前述のカーネル起動のオーバーヘッドと並列をループ内で利用する区間が原因であり，レジス

度の不足による性能頭打ちに起因すると考えられタスピルが頻発したことによるものであった。

る。加えて，⑧消費電力測定と連携し，ではレジスタ数が多いために起こらなかった。

については，の電力測定をで顕著な性能劣上で試行した。また，⑩

自動チューニングと連携して，を題材にサ化がみられる。これは巨大な配列の２次元目を切ブルーチン内の任意の変数の浮動小数点精度型をり出して，順にサブルーチンに渡して計算を行う

変更し，計算結果と計算速度を同時に評価する手区間が原因であった。１次元目の配列サイズは法についての具体的なテスト実験を進め，アプリで計算するには十分な大きさであり，データはケーション開発側の要望をフィードバックしを用いたステンシル計算の主要部分に上に保持されたまま計算が行われていたので，関する評価は上での計算は想定された速度が得られ年度の段階で一定の成果を得ていた。しかし，サブルーチンの外で２次元目のることが出来たが，大規模シミュレーションを実ループが回転するたびにのカーネル起現するためには，巨大な実アプリケーション全体動が頻発し，そのオーバヘッドが無視できないほに変動精度演算を適用する必要がある。今年度はど大きいことがわかった。

のカーネルは諸物理過程の最適化，混合精度化を進めた

サブルーチンをまたぐことが出来ない。上記２つが，その中でいくつかのプログラムコードの構造の性能劣化を解決するためには，コードリファク的問題が浮き彫りになった。今後は得られた知見

タリングを行う必要がある。

をもとにコードの改良を進め，大規模シミュレー
ションと消費電力性能評価を進めていく予定である。また，自動チューニング手法の研究と連携
し，ループ内に倍精度と単精度の計算を混在させ等での高速化を両立でざるを得ない時に，

きるコーディング等についても検討を進めていきたい。
③疎行列演算（混合精度））機械学習と混合精度・反復改良法を利用した連図
を用いた
物理過程の計算時間評価，それぞれ倍精度計算（左），単精度計算の結果を示す
ではスレッド数
でハイブリッド並列を用いており，
ではスレッド並列は用いていない

立一次方程式解法〔〕メニーコアプロセッサを用いたにおける近年の技術的課題のひとつは消費電力であり，省
電力かつ省エネルギーなアルゴリズムの開発，実用化が期待されている。混合精度演算については

また，各スキームのアダプタ部分でのデータコ
上で無視できないほど時間がか

既に多くの研究があるが，低精度利用時の精度保かっている。今後は諸物理過程ドライバ全体を証の研究は実施されているものの，実問題での大ピーは特に上で計算するよう，

規模疎行列への応用例はほとんどない。

適用領域を拡イブラリ

げていく必要がある。倍精度と単精度の比較にお
7

に

のラ
のサブ

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

ルーチンが追加されるなど，混合精度を用いた

り，便宜上，増加率の昇順にソートしている。1 反Iterative Refinement の手法も近年注目されている

復あたりの時間増加率は，行列でいかにばらつき

が，Iterative Refinement は必ずしもすべての行列にを生じるかがわかる。
2.50

法の収束判定値も適切な値を設定しないと効果が

2.00

Acceleration

有効ではなく，Iterative Refinement 内にある反復解得られない。また，機械学習への期待から，疎行列を対象として最適な数値計算ライブラリ選択に

1.50

1.00

関する研究が行われているが，Iterative Refinement

0.50への応用例はない。

0.00
0

本研究では解の精度を改善するための反復改良

20

40

50

80

100

120

Matrix ID

図 10 CG 法 1 反復あたりの時間増加率

法(Iterative Refinement)に着目し，有限要素解析で生成される大規模疎行列連立一次方程式の求解のための，高精度な混合精度演算手法を開発するこ次に，混合精度演算を用いた Iterative Refinementの数値実験を行った。混合精度演算を用いたとを目的とする。

Iterative Refinement のアルゴリズムを図 11 に示す。

Iterative Refinement のアルゴリズムではインナーループで単精度共役勾配法を実行するため，本研究の準備フェーズとして，単精度演算，倍精度
演算それぞれで共役勾配法の計算を実行し，収束判定条件を変化させて反復回数・計算時間を計測した。SuiteSparse Matrix Collection より正定値対象の疎行列を選び，右辺ベクトルは厳密解がすべて1 になるように定めた。また，プログラム中で疎行列は全て CRS 形式に変換して処理している。また，単精度演算を用いた共役勾配法の計算では，係数行列・右辺ベクトル共に単精度で計算している。
結果は次の 3 ケースの特徴に分けられることがわかった。

図 11 Iterative Refinement 最終的な収束判定値を10-8 に設定しインナーループの単精度 CG 法の収束判定値 η を変化させて計算時間を比較した（図12） <ケース 1> 反復回数は単精度も倍精度もほぼ等しく，計算時間が単精度の方が常に短い <ケース 2> 単精度の方が反復回数は多いが，Computvtion time of Itervtive Refinement400

収束判定値が緩いときは単精度の方が計算時

350
Computation time [s]間が短い <ケース 3> 倍精度が単精度よりも反復回数も少なく計算時間も短い

Computvtion time of CG solver by double precision

300
250

Best

200
150
Itervtive Refinement is slower100

また，共役勾配法の 1 反復にかかる単精度，倍50

精度それぞれの計算時間をさまざまな行列に対し

0
1.0E-01て測定し，単精度に対する倍精度の時間増加率を

1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E-07
Convergence criterion η of CG solver by single precision

図 12 インナーループ収束判定値 η の影響調査した（図 10）。図 10 の横軸は行列の番号であ8

1.0E-08

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

度（前処理・加法シュワルツ法部分のみ単精度）
図にの疎行列ソルバーを適用し，均質問題において，による計算結果の例を示す。このとき用いた疎行列は行列サイズが，非ゼロ要素数が図

領域倍精度の場合と同じ反復回数で同じ解を得た。

の行列である。

より，インナーループの単精度

判定値 ̀ の値を

全年度以降は不均質問題，動的負荷分散に継続して取り組む予定である。

法の収束に設定した時が最も計算が速く，̀ の設定値によっては倍精度法に比べて遅いケースがあることがわかる。
このようには適切なインナーループの単精度
ないと倍精度法の収束判定値 ̀ を設定し法に比べて良いパフォーマンス

図

並列有限要素法（ 分割）が得られないことがわかった。
今後は，機械学習を用いて行列の情報から ̀ の④疎行列演算：固有値分布を保存する行列縮退法

値を設定することを考える予定である。

について
有限要素法などから導かれる疎行列を係数行列）局所的な変動精度適用年とする連立一次方程式を反復法で解く場合，収束では，条件は係数行列の固有値分布に依存し，固有値分問題の条件に応じて局所的に演算精度を変動させ布を改善するための様々な前処理手法

る可能性について指摘された。

（

月に出席した国際会議

年度は，大規）が研究されている。近年，機械学

模分散並列問題で局所的・動的に演算精度を変動，習等の手法を用いて問題，すなわち係数に応じた

動的に負荷分散を適用する手法の開発を進めてい最適な前処理手法を選択する研究が盛んになってる。まず，領域分割に基づく並列有限要素法においるが，係数行列の固有値分布は最適前処理手法

いて，領域毎に異なる精度を適用した場合の検証

選択にあたって重要な指標となる。大規模な係数を，並列有限要素法による三次元構造解析におけ行列の固有値分布の算出には連立一次方程式を解

る加法シュワルツ法（

くのと同等のコストを要するため，効率的な算出

（

），図

）に基づく

法

手法が必要となる。本研究では，元の大規模行列ソルバーに対して実施した。

と同様の固有値分布を保ちながら，行列を縮退（）する方法について検討した。
本研究では，以下に示す種類の行列縮退手法

を提案している：




（

）

（

）


図
法シュワルツ法 による並列前処理安定化

加第 1 の手法（RDPM，ランダムな対角投影法）は固有値分布のスペクトルを，保持するものの最大分割した場合，固有値の計算精度が低い（図 14）。逆に第 2 の領域については倍精度， 領域については混合精方法（RCBPM，列単位投影法）は固有値分布の図に示すような規則形状を9

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

スペクトルは不正確であるが，最大固有値をより列群に対して，第 3 の手法を適用することによっ正確に求めることが出来る（図 15）。

て，100×100 程度の行列で大規模な疎行列の固有値分布を再現できることを確認した（図16）。
2019 年度は，固有値分布を保存する行列縮退について予備的な評価を実施した。今後は，縮退処理の並列化に取り組んで行く予定である。
⑤疎行列演算：図 14 第 1 の手法（RDPM，ランダムな対角投影法）による固有値分布計算（―：固有値を実際に求めた場合，―：RDPM による計算）（1138_bus（SuiteSparse Matrix Collection），Power
Network））精度変動演算用
〔

の

実装

〕
比較的条件の良い問題では，倍精度演算のかわりに単精度）を使うことによって計算時間を短縮できる場合があるが，機械学習などでは有効桁数が使われる半精度

桁程度のため，実問題に使用することは困難である。本研究
では

と

型である

の中間的な特性を持つデータを定義し，有限要素解析における前は符号処理において活用した。
数部

図 15 第 2 の手法（RCBPM，列単位投影法）による固有値分布計算（―：固有値を実際に求めた
場合，―：RCBPM による計算）（1138_bus（SuiteSparse Matrix Collection），Power
Network）

ビット・仮数部

からなり，
仮数部

符号

ビット・指

ビットの合計ビット

ビット・指数部

ビット・ビットの合計ビット と指数部のビッはと同じレンジを持つもののデータサイズがとなるため，メモト数が同一となる。

リバンド幅ネックのカーネルにおいては実行時間
がになると期待される。なお，現在の主要な計算機においてはの演算器はサポートされていないため，変数のメモリへの格納時のみにを使い，計算時には演算する。上記を
図 16 第 3 の手法（MP，混合投影法）による固有値分布計算（―：固有値を実際に求めた場合，―：MP による計算）（1138_bus（SuiteSparse
Matrix Collection），Power Network）

を

に変換して

にて実装し

で

計測した結果，想定通りの性能が得られ，の有効性が示された。本研究ではさらに上記をで実装し，ソースコードを公開することでより多くのアプリケーションで第 3 の方法は第 1 の手法（RDPM，ランダムなきるようにした〔

対角投影法）は第 2 の方法（RCBPM，列単位プを活用で〕。性能計測の結果，実際の地

震アプリケーションに組み込んだ際にロジェクション）を組み合わせる混合投影法によ

版においてって，最大固有値，固有値分布スペクトルの両者を確認した〔 〕。

を正確に模擬できるような行列縮退手法を提供するものである。SuiteSparse Matrix Collection の行
10版と同等の性能が得られること

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

FP16 演算を使うことができるようになる。また，レジスタ上でデータを使いまわす等の工夫をすることでデータアクセスコストを抑えた。開発手法
により Summit のほぼ全系にあたる 4544 ノード
図 17 低精度データ型 FP21 の模式図

(V100 GPU 27,264 枚)を使って 1.67×1012 自由度の2）AI 用低精度演算器の Equation-based Modeling

問題を解いたところ，ベースラインソルバー（3×3への活用〔4〕

ブロックヤコビ前処理による共役勾配法ソルバー）

2019 年 11 月現在で Top 500 リスト一位の米国比で 17.6 倍速，V100 GPU 用にチューニングを施Summit は主に AI や data analytics 用に設計されたした SC14 ゴードンベルソルバーに比べて 3.89 倍

Tensor Core を搭載する V100 GPU を多数搭載しての高速化を達成した。この際のピーク性能はソルいる。この Tensor Core では，FP16 において小規

バー全体で 416 PFLOPS, 行列ベクトル積部分で模な行列行列積を通常の FMA 演算器に比べて 4

1.10 Exa FLOPS となった〔4〕。このような，アル倍のピーク性能で実行できるため，物理シミュレ

ゴリズム的に equation-based modeling と AI を融合ーションに活用することで性能改善の可能性があさせる新しい HPC の方向性は他のアプリケーシる。その一方で，演算パターンが行列行列積に限ョンへも活用可能であると考えられ，今後の発展定されている点，演算精度が低い点，演算性能がが期待される。

高いためにデータアクセスがボトルネックになり⑥ H 行列〔5,6,17,18,27,30〕
本研究では，これまで H 行列，H2 行列，HSS 行

やすい点，の 3 点が制約となり equation-basedmodeling において使用が普及するに至っていない。

列，HOLDR 行列，BLR 行列など様々な階層的低

そこで本研究では，equation-based modeling の一例

ランク近似法（以下ではこれらの手法を簡単のたである非構造格子有限要素法を対象にこのようなめにまとめて H 行列と表記する）に対応できるラ演算器を活用できるアルゴリズムの構築例を示す。
通常の非構造格子有限要素法においては各要素

イブラリの開発を行ってきた。混合精度演算を行う場合，H 行列の中で最も工夫を要するのは，密の形状と節点の接続状況に応じた非均一な基底関行列を低ランク行列に圧縮する部分である。

数が使われるため，ランダムアクセスなどの複雑

Randomized SVD や Interpolatory Decomposition なな計算が生じ，結果として行列行列積に特化したど様々な圧縮方法が提案されているが，いずれのTensor Core の高い性能を発揮することが難しい。
そこで本研究では，AI の学習において用いられる手法も QR 分解を基礎としているため，2019 年度の JHPCN 課題では，多くの小さな QR 分解を GPU

計算パターンと類似の計算が生じるよう，ソルバ

上で並列に行う batched QR 分解を開発した。このー内で局所的・規則的な基底関数展開を使う。こ際 に ， Volta 世 代 の GPU に 搭 載 さ れ て い る

こでは，四面体二次要素メッシュを対象とした陰TensorCore を利用することで 10TFLOPS を超える解法ソルバーにおいて adaptive conjugate gradient

性能を得ることができた〔27〕。ただし，行列の大method を使い，前処理部に構造格子からなる共役

きさは 16x16 に限定されており，このままでは H

勾配法ソルバーを使う。このようにすることで，行列の圧縮に用いることは難しい。今年度はこれ

ソルバー内の大部分の計算を構造格子における演を TSQR に拡張し，任意の大きさの行列の低ラン

算に移すことができ，その内部の主要演算カーネク近似を高速に行うための QR 分解の枠組みを開ルとなる行列ベクトル積において，小規模な行列・発した〔6，30〕。Z 図 18 に行列サイズを変えたと行列積を用いることが可能となる。また，使用は

きの TSQR の V100 上での計算時間を示す。TSQR

前処理に限定されるため，低精度演算を用いてもでは列数 16 までの行列に対する QR 分解を行い，十分な精度の解を得ることができ，Tensor Core の11

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

行数が大きい場合にが開発している線形ができる。しかし再直交化を行うことで計算量はと比較して高速に代数ライブラリ増加するため，並列性の抽出・計算のオーバーラ分解を行うことができていることが確認された。
は積を計算する

つの行列を半精度として入力する必要があり，を用いて単ップを行うことで増加する計算時間を隠蔽することが重要となる。

精度行列積を計算する場合は半精度への変換を行う必要がある。これにより行列積の計算精度が単精度演算器で計算した場合と比較して劣化するが，半精度への変換により失われる仮数部を別途保持し，これを用いて行列積の計算精度を補正することができる。この精度補正計算を用いた場合，による計算量は 倍となるが，利用した場合計算時間は
り，

に

図

の計算時間

図

の相対残差

図

の直交性

倍程度にとどまっておを用いた場合と比較しても高速に分解が行われていることが確認された。一方で計算精度として相対残差（図）及び直交性（図）を見ると，精度補正を行うことでを用いた単精度分解により近い精度で計算を行えていることが確認された。

分解を行うのに必要な作業用メモリ量では，本研究でのを用いた分解と比較しは～の作業用メモリを削減できており，より大きな行列に対する分解が可能であることが確認され

た。
の実装では

変換を用いた分解を用いた。この分解では通常陽的に行列を作らず，行列ベクトル積とベク
トルの直積を用いて計算を行うことで理論計算量を削減する。

を用いてこの行う場合，陽的に分解を一方で，開発した行列を作りベクトル

クであることを前提としており，フルランクでなの直積と行列行列積を用いた方が高速に計算を行い行列が入力された場合には計算が破綻する。

うことができることが確認された。このようにアで行う

ルゴリズムを変更することが行列行列積を計算する回路であるでは列数がため，計算の破綻を防ぐ，または破綻を検知する
仕組みを導入する必要がある。また，今年度の研という制限がある

究によりため，この制限を解消し一般的な行列に対する分解を行うためにでは分解では，フルランク

でない行列が入力される場合が十分に想定されるをより効率的に利用する上で重要となる。
開発したは入力行列がフルランによる精度補正単精度行列積

を用いた

を用いる。

精度に

の直交性が劣化することが知られているが，再直交化を行うことで直交性の劣化を抑えること
12

では，
進数で

れている（図と比較して計算

桁程度の精度の劣化が確認さ，）。 こ の 精 度 劣 化 に よ り学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書の近似精度がどう影響を受ける

本研究では，有限体積法（７点差分メッシュ）に
よって離散化して得られた連立一次方程式に対すかを調査する必要がある。

る精度保証の改善に成功した。解析モデルを図 1,を用いた精度補正単精度行列積の既倍となる。 図 22 に示す。

による計算量が存手法では本研究ではこの精度補正計算の一部を省略し計算倍に削減することで高速に精度補正計算を量を行えており，

分解の精度の劣化も確認されていない。今後この精度補正計算の省略による計算精度への影響を注意深く調べるとともに，精度補正を効果的に行える行列の条件や逆に精度補正を行
いにくい行列の調査，行いにくい行列に対するア
プローチの調査を行っていく。
今後は行列の分散並列化における（境界は，すべてや図などのランタイムを用いた非同期実行
〔

〕や

行列の

分解〔

〕，

境界条件）
方程式）

解析モデル（三次元並列化における独自のランタイムを用いた非同期バッチ実係数行列が行〔 〕などへと証法〔による混合精度の精度行列性を持つ場合，高速な精度保〕を適用可能であるが，誤差限界か過大評価になっている可能保証を拡張していく予定である。

性がある。そこで，それを抑制するための方法
⑦精度保証

〔

）実問題向け精度保証手法の開発（疎行列）体的には，以下のようなアルゴリズムとなる。

〔〕を適用した。具

〕実問題に精度保証を適用するための基盤として， i.離散化して得られた連立一次方程式  を三次元解く。得られた近似解をとする。

方程式の数値解の精度保証を行った。数値シミュレーションでは，図ii.連立一次方程式   (はすべての要素が 1に示すようのベクトル)を解く。得られた近似解をとする。

に，微分方程式を離散化し，得られた連立一次方程式を解くことが多いが，離散化誤差については， iii.上記のを用いて係数行列の M 行列性の保証を行う。

独立に考えることができるため，本研究では，連iv.残差  

立一次方程式に対して精度保証を適用する。

を精度保証付きで計算する。残差の近似値を̀ ，誤差限界を௥ とする。

v.連立一次方程式  ̀ を解く。得られた近似解

PDE (3D Poisson equation)を ̀ とする。

Discretization by FDM, FEM, etc.

vi.以下の誤差評価式の右辺について，浮動小数点

Discretization error
Rounding error

演算の丸めモードの変更を適宜使用しながら，

Linear ssstem Ax n b
Iterative solver, such as
CG, GMRES, etc.
Truncation error
Rounding error

厳密な上限を求める。
At present,
we consider
only this part.

 ̀ ஶ ஶ    ̀ ஶ ஶ

（但し，右辺第 2 項の分母が 0 以下になった場

Approximate solution
図 21

ஶ

合は精度保証失敗）

今回の精度保証の対象
13

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書（ ノー

東京大学情報基盤センターのド）において数値実験を行った。連立一次方程式のソルバには法（マルチカラー前処理を用いて，による並列化（を行った。

計算時間（秒）表

精度保証

近似解の計算ଵ ଶ ଵ ଶ 

色）

スレッド）

଺

法の反復の停止条件は残差ノルム
図を用いて：
  ,   ,   ̀ , ̀

価が大幅に抑制できていることがわかる。また，ଶ ଶ  ଵ

表ஶ  ଶから，近似解と精度保証の計算時間は同程度であることがわかる。

 ̀ ଶ ̀ ଶ  ଷ

精度保証については，三次元方程式を離散化して得られた連立一次方程式に対しては，実のように与えた。但し，それぞれの残差ベクトルの計算についてはから，提案方式によって誤差限界の過大評

用的かつ高精度な精度保証法の適用に成功した。

法の反復中に計算される

今後は，単精度演算等の低精度演算での数値実験，

残差ベクトルを用いた。また，浮動小数点演算は，離散化誤差の評価方法の検討，階層型行列法におすべて倍精度を用いた。

ける精度保証法の検討等が課題である。

ここでは，以下のような問題設定とする。

）不連続係数ポアソン方程式に対する離散化誤差 ଵ  ଵ 

のଶ 評価〔に固定, ଶ を変化させる。

ିଵଶ, ଶ ିଶ, ଷ 

図 ，図ିଽ〕で説明した不均質場における温度分

布モデルに対する離散化誤差の評価を実現した。
本評価と前述の連立一次方程式に対する誤差評価

離散化して得られた連立一次方程式の近似解をを組み合わせることにより，対象モデルの真の解

計算し，精度保証法を用いて，最大相対誤差：とコンピュータによって得られた数値解のଶ 誤差を包括的に評価することができる。対象問題は同モデルを不連続係数ポアソン方程式－Ȃ Ȃ  Ȃ  

の上限及び相対残差ノルム̀$̀として定式化した楕円型偏微分方程式である。この上限を求めた。結果を図及び表こで Ȃ ଶ ̀は任意に与えられた関数とし， はに示す。

図 1 で説明した各部分領域で定数値をとる不連続係数である。本境界値問題とその有限要素法で離散化された方程式の解 ௛ に対して，まず以下の形

式での誤差評価を得た。

&－௛ &ఒ'Ă௅మここで，誤差評価式中のノルムȂఒ は拡散係数の図選び方に依存し*ଶఒ  Ȃ௞ ௞ Ȃ*௅ଶమ ሺஐೖሻ で与えら条件（ଵ ଶ ）毎の最大相対誤差（実線）と相対残差ノルム（破線）れる。また'Ăのとり方は離散化手法に依存する。

例えば立方体(正方形)メッシュ上で各方向での線
14

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

形近似を採用した場合は，メッシュの一辺の長さ質場における温度分布モデルに対する離散化誤差

をℎとしてのଶ 誤差評価を実現することができた。本成果を疎行列計算に対する精度保証法と組み合わせるこଵℎିℎ = max ௞ ଶ
 ௞とにより，対象問題の離散化誤差を含めた包括的誤差評価を高速に達成する理論基盤が整った。

と設定すれば良いことを明らかにした。このとき，

2020 年度には上記をベースとした手法を実装し，コンピュータ上で表現された実際の有限要素解のノルム⋅ఒ の意味での誤差評価はℎオーダーであすべての誤差（離散化誤差＋丸め誤差）を定量的り，またその誤差は領域中で最も小さい拡散係数に評価する。特に，離散化パラメータと演算精度

௞ に依存する。上記誤差評価式と Aubin-Niche ののバランス，例えば「倍精度など特定の演算精度技巧による式変形により，以下のଶ 誤差評価（ℎଶに対してどの程度メッシュサイズℎを小さくとる

オーダー）が得られる。

ことが妥当なのか」という問いに対してある種の － ௛ ௅మ ≤ (ℎ)ଶ  ௅మ

指標を与えることが今後の目標である。
また，2019 年度はଶ ノルムの意味での離散化誤

例 え ば ， 図 1 ， 図 22 の 設 定 下 で λଵ = 1, ଶ =ିଷ

差評価を得たが，2020 年度はこれをベースとしஶି଻

10 , ℎ = 10 とした場合ノルムの意味での誤差評価を目指す。これにより
真解 と数値解  ௛ の領域各点での差異を見積もるℎ ≤ 7.87 × 10ିଶことが可能となる。

ℎଶ ≤ 9.78 × 10ିଷのように評価できる。これに ( , , ) =+⑧消費電力測定：疎行列格納法，問題規模による+のଶ ノルムであるȂ2.5を乗じたものが，それぞれ低精度演算の計算時間，消費電力・エネルギーへ － ௛ ఒ と － ௛ ௅మ の上界となる。

の影響〔7,12,14,15,19,20,23,32,33,35〕以上の結果と前述の連立一次方程式に対する誤本研究では，図 1 のような形状に対して，有限差評価を組み合わせることにより，所望の真解体積法によるポアソン方程式ソルバーから導かれと数値解  ௛ (有限要素解のコンピュータによる数

る対称正定な疎行列を係数とする連立一次方程式

値的近似௛ )の誤差を包括的に評価することがでを不完全コレスキー分解前処理付き共役勾配法

きる。最も簡単には，連立一次方程式の ஶ 誤差評（ Preconditioned Conjugate Gradient Method by

価Incomplete Cholesky Factorization，ICCG 法）を対象として，倍精度・単精度演算の比較，消費電力，消 ௛ −  ௛ ஶ ≤ が得られ，かつ先と同様に費エネギーの測定を実施し，条件の比較的良い問௛,

題では，単精度演算により倍精度の場合と比較し ௛ がメッシュ内でて，50～60％程度の計算時間，消費エネルギーで各方向に線形補間されてる場合を考える。この場正しい計算を実施することができた〔23〕。

合，真解 と数値解  ௛ のଶ ノルムの意味での誤差

これまで実施してきた研究は，低精度演算の有は：効性を示すための予備的研究であるが，2019 年度
 －  ௛ ௅మ ≤  － ௛ ௅మ +  ௛ −  ௛ ௅మは，問題サイズ，行列格納手法，対象アーキテク≤ ℎଶ  ௅మ + |Ω|チャなどの様々な影響を考慮した検討を実施した。
本研究では以下の 5 種類の計算機環境を使用しで与えられる。

た：2019 年度の研究により，図 1，図 22 にある不均15

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

表3

マルチコア・メニィコア CPU

GPU 概要システム名
略称

 Reedbush-U（東京大学情報基盤センター）
 OFP-Mini（東京大学情報基盤センター）

GPU 名称 Oakbridge-CX（東京大学情報基盤センター）

Reedbush-L
P100
NVIDIA
Tesla P100

P9 Cluster
V100
NVIDIA
Tesla V100

5,300

7,500

16

32

557

855

250

300

GPU

理論演算性能（GFLOPS）
 P9-V100 クラスタ（東京大学情報基盤センター） 主記憶容量（GB）
メモリ性能（GB/sec）
STREAM Triad
本研究では，CPU については各 1 ノード，GPU に TDP (W) Reedbush-L（東京大学情報基盤センター）

ついては GPU1 基を使用して実施した。表 2 は CPUの各 1 ソケット，表 3 は各 GPU の概要である。

元のプログラムは Fortarn+OpenMP によって開

OFP-Mini は，最先端共同 HPC 基盤施設（JCAHPC）

発されているが，GPU 用にはこれに OpenACC をの運用する Oakforest-PACS と同じ CPU を使用し適用したバージョンを使用している。

た 8 ノードのクラスタである。また，P9-V100 ク

電力評価には RBU，OBCX については Intel がラスタは，IBM Power 9 と NVIDIA Tesla V100

提供する CPU の電力制御機能の RAPL を使用し（V100）から構成されるクラスタである。

た。また OFP については，同じく Intel の提供するIntel pcm-power.x を使用し，GPU（P100，V100）で表 2 計算機（CPU）各ソケット概要
システム名
略称

CPU 名称

ソケット当たりコア数
ノード当たりソケット数
理論演算性（GFLOPS）
主記憶容量
（GB）

Reedbush-U
RBUIntel Xeon
E5-2695
v4(Broadwell-EP)
18

OFPmini
OFPIntel
Xeon Phi
7250(Knights
Landing,
KNL)
68

は nvidia-smi を使用した。また，GPU 利用時は CPUの消費電力は測定していない。

OakbridgeCX
OBCX

本研究では，以下の 3 種類の行列格納形式について検討した（図 24）：

Intel Xeon
Platinum
8280(CascadeLake, CLX)

 Compressed Row Storage（CRS）
 Column-Wise な Sliced ELL（ELL）
 SELL-C-（C=8,  =1）（SCS）

28
C

2

1

2

604.8

3,046.4

2,419.2128

メモリ性能（GB/sec）
STREAM Triad

65.5

TDP (W)

120

MCDRA
M: 16
DDR4:
96
MCDRA
M: 490
DDR4:
84.5
215

s

96
CRS

ELL

Sliced ELL

SELL-C-s

図 24 様々な疎行列格納形式
101.0
図 25 は SELL-C-（SCS），OPT における前処理
205

計算・前進代入部分の実装方法である。
CRS⇒ELL/Sliced ELL⇒SELL-C-に従って，計算密度は増加する。また，RBU，OFP，OBCX につ
いては SCS を更に最適化してスレッド並列化によるオーバーヘッドを取り除いた実装を適用した。
16

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

本稿ではこれを「OPT」と呼ぶ〔星野他，OpenACC

 a:

Cache モードを用いた ICCG 法ソルバーの Pascal GPU における b:Flat モード，DDR4 のみ使用

性能評価，情報処理学会研究報告（2017-HPC-158-

 c:Flat モード，MCDRAM のみ使用

18），2017〕。

(b)

問題規模としては：

SCS: SELL-C-s!$omp parallml (…)
do ic= 2, NCOLORtot-1!$omp dodo ip= 1, PEsmpTOT
iq1= SMPindex((ip-1)*NCOLORtot + ic-1) + 1
iq2= SMPindex((ip-1)*NCOLORtot + ic)
ip0= (ip-1)*NCOLORtot + ic
iq0= (iq1-1)*8
do ib = iq1, iq2
ib0= (ib-iq1)*8!$omp simd
do is = 1, 8i= iq0 + ib0 + is
VAL= W(i,Z)
do k= 1, 3
VAL= VAL- AL(ib0+is,k,ip0)*W(itmmL(ib0+is,k,ip0),Z)
mnddo
W(i,Z)= VAL * W(i,DD)
mnddo
mnddo
mnddo
enddo
!$omp mnd parallml!$omp parallml (…)(…)
do ic= 2, NCOLORtot-1

 Tiny（t）:

32,768 DOF（=323）

 Small（s）:

262,144 DOF（=643）

 Medium（m）: 2,097,152 DOF（=1283）
16,777,216 DOF（=2563）

 Large（l）:の 4 種類を実施した。Large の場合でも，P100，
OFP（MCDRAM）の各メモリに充分収まっている。

OPT

(a)
iq1= SMPindex((ip-1)*NCOLORtot + ic-1) + 1
iq2= SMPindex((ip-1)*NCOLORtot + ic)
ip0= (ip-1)*NCOLORtot + ic
iq0= (iq1-1)*8
do ib = iq1, iq2
ib0= (ib-iq1)*8!$omp simd
do is = 1, 8i= iq0 + ib0 + is
VAL= W(i,Z)
do k= 1, 3
VAL= VAL- AL(ib0+is,k,ip0)*W(itmmL(ib0+is,k,ip0),Z)
mnddo
W(i,Z)= VAL * W(i,DD)
mnddo
mnddo
!$omp barrimr
enddo
(…)
!$omp mnd parallml

2.50

Ratio of Single/Double: Comp. Time

(a)

2.00

RBU-t
OFP/c-s
P100-tV100-sRBU-sOBCX-tP100-s

OFP/c-t
OBCX-s
V100-t

ELL

SCS

1.50

1.00

0.50

0.00
CRS

OPT

Matrix Storage(b)

RBU，OBCX については 1 ノード（2 ソケット）を使用した。RBU では合計 36 コア（36 スレッド），
OBCX では 48 コア（48 スレッド）を使用した。
また OFP では，先行研究における最適な設定である 64 コア（128 スレッド）を使用した。P100，V100
は GPU1 基を使用した。

5.00

Ratio of Performance: Computation Time

図 25 前処理計算・前進代入部分の実装方法（a）
SELL-C-（SCS），（b）OPT〔星野他，2017〕

4.00

RBC-D
OFP/c-D
OBCX-DP100-D
V100-D

RBC-S
OFP/c-S
OBCX-S
P100-S
V100-S

3.00

2.00

1.00

0.00
CRS

OFP は表 2 に示すように，通常の DDR4 メモリ

ELL

SCS

OPT

Matrix Storageに加えて，CPU パッケージに直接収められている

図 26 ICCG 法ソルバーの計算性能，「倍精度・CRS」の計算時間によって無次元化，（a）
Medium，（b）Large，1.00 より大きい場合は「倍精度・CRS」よりも ICCG 法の計算時間が短い

高速の MCDRAM を使用することが可能である。
Cache モードでは，MCDRAM を DDR4 のキャッシュとして使用ができる他，Flat モードでは，

本研究では，全計算を倍精度で実施した場合

DDR4 と MCDRAM を使いわけることが可能であ（Double）と全て単精度実施した場合（Single）のる。本研究では，以下の 3 ケースについて検討を

2 種類の計算を実施した。単精度のみの場合は，倍実施した：

精度の場合と比較して 20％程度反復回数が増加している。
17

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書
Ratio of Performance: Computation Time

5.00
RBC-t
OFP/c-s
P100-tV100-s

4.00

RBC-s
OBCX-tP100-sが比較的良い。また

少ないが，三者の中では

OFP/c-t
OBCX-s
V100-t

倍精度⇒単精度の速度向上は
図3.00％程度である。

は単精度を使用して実施した法部分の計算時間について，各計算機環境，各問題サイズ
における，
「倍精度・

2.00

」の計算時間によって無次元化したものである。値が1.00

は，「倍精度・より大きい場合」より単精度演算による法部分の計算時間が短いことを示している。

0.00
CRS

ELL

SCS

OPT

Matrix Storage
Ratio of Performance: Computation Time

5.00

4.00

RBC-m
OFP/c-m
OBCX-mP100-mV100-m

RBC-l
OFP/c-l
OBCX-l
P100-lV100-l

CRS

ELL

3.00

2.00

1.00

0.00
SCS

OPT

Matrix Storage
5.00

図

は（ ）

，
（ ）

の場合に

Ratio of Performance: Computation Time図法ソルバー（単精度）の計算性能，「倍精度・」の計算時間によって無次元化，
（ ）
，
，
（ ）
，
：
より大きい
場合は「倍精度・」よりも法の計算時間が短い法

部分の計算時間について，各計算機環境，各問題
サイズにおける，「倍精度・」の計算時間によって無次元化したものである。値が
い場合は，
「倍精度・

より大き
，

でやや遅くなり，，速くなっている。

，

SCS

OPT

図
法ソルバー（単精度）の計算性能，「倍精度」計算時間によって無次元化，（ ），，（ ），より大きい場合は「単精度」の計算時間が短い

は
でまた

図

（ともに

（ ）は

，

，図

（ ）は

，

については，

モード，のみを使用した，ケース （）のの場合である。

場合の結果である。全般的に問題規模が小さい場向け最適化であ合には，単精度適用による速度向上の効果は少なにおいては倍精度⇒単精度の速度向上

く，問題規模，データ移動量が大きくなるほど顕見られる（元々はELLと順調に性能向上が

⇒

るが）。

1.00

Matrix Storage

サイズにおいて顕著である。一方

では

2.00

CRS

プロセッサ）では倍精度⇒単精度の速度向上は特
に

3.00

0.00

で倍精度，単精度ともに計算性能が向上しているが，

4.00

RBC-S
OFP/c-S
OBCX-S
P100-S
V100-S

法部分の計」より算時間が短いことを示している。
⇒

RBC-D
OFP/c-D
OBCX-DP100-D
V100-D

⇒⇒，はと比較して少ない。特に十分なためと考えられる。

（

では単精度による速度向上が得られず，むしろ遅化が不速度向上が少ないのはベクトル化，，については，

著となる傾向がある。またでくなっているが，）の，，速度向上が得られており，図傾向は似ており，また行列格納方法による差異は
18

では

倍以上の

に示した場合と同学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

様，CRS⇒ELL⇒SCS・OPT と演算密度上昇とともためメモリ（DDR4）の消費電力は極めて少ない。

に順調に計算速度が増している。特に，「単精度・また，問題サイズが大きくなると OBCX ではメモ

OPT」の場合には問題規模が Medium，Large では，リ，CPU/Package ともに消費電力が増加している

AVX512 によるベクトル化の効果が顕著であり，ことがわかる。OFP/c では Medium の SCS，OPT で反復回数が 20%増加しているにもかかわらず，「倍は演算密度が CRS，ELL と比較して高いにもかか

精度・CRS」の場合の 4 倍以上の速度向上が得らわらず消費電力（W）が低下している。

20(a) 250.00

10

れている。OFP では全般的に問題規模増加ともに速度向の効果が顕著であるが，OBCX の場合は，MemoryPower Consumption (W)

問題規模増加とともにむしろ速度向上は低下している。
図 28 は各ケースにおいて単純に単精度と倍精度の計算時間を比較しているが図 27 と同様の傾向が見られる。

200.00CPU

150.00100.00

50.00

1.50
RBU-t
OFP/c-s
P100-tV100-sRBU-sOBCX-tP100-s

OFP/c-t
OBCX-s
V100-t

0.00

1.00

(b)

600.00
Memory500.00Power Consumption (W)Ratio of Power Consumption

(a)

0.50

0.00
CRS

ELL

SCS

OPT

Matrix Storage(b)
Ratio of Power Consumption

1.50
RBU-m

RBU-l

OFP/c-m

OFP/c-l

OBCX-m

OBCX-l

P100-mP100-lV100-mV100-lCPU

400.00300.00200.00100.000.00

1.00

図 30 ICCG 法ソルバーにおける消費電力（W）（単精度）（a）OFP/c，（b）OBCX（s：Small，m：Medium，l：Large）

0.50

図 31，図 32 は RBU，OFP/c，OBCX，P100，V100
各計算機環境の最適ケースにおける問題サイズ

0.00
CRS

ELL

SCS

OPT

「Medium，Large」の場合の ICCG 法の（a）消費Matrix Storage

図 29 ICCG 法ソルバーの消費電力，「倍精度・CRS」の消費電力（W）によって無次元化，（a）
Tiny，Small，（b）Medium，1.00 より大きい場合
は「倍精度・CRS」よりも消費電力（W）が大電力（W），（b）消費エネルギー（J）と計算時間である。消費エネルギーと計算時間はほぼ正比例していることがわかる。計算時間としては，OFP/c，

図 30 は（a）OFP/c（MCDRAM のみを使用した，

OBCX，P100，V100 はほぼ拮抗しているが，OBCX，

ケース c），（b）OBCX の各問題サイズ（Small，V100 がやや速い。OBCX は消費電力（W），消費

Medium，Large）における消費電力（W）の推移でエネルギ（J）が OFP/c，P100，V100 の 3 倍以上である。図 30 ではメモリ■，CPU/Package■の消費ある。OFP/c，P100，V100 は消費エネルギ（J）は電力をそれぞれ表示しているが，OFP/c では CPU

ほぼ等しいが，P100，V100 についてはホストの上に直接装着された MCDRAM のみ使用している

CPU の消費電力，消費エネルギー（アイドル時で19Power Consumption (W)

50

40

30

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書
600.00(a)

1.80⇒V100 の計算速度向上は Medium で 35-40%，Large

MemoryPackage400.001.60

で 45-60%である。全般的に消費電力（W）は

Time
1.40

400.00Medium⇒Large で増加しており，特に OFP（OFP/c）

1.20

Power Consumption (W)1.00

では 2 倍程度になっている。

sec.

300.000.80

200.000.60

250.00(a)4.50

Memory

CDU

Time

0.40

4.00

0.20

0.00

0.00
RBU/S

OFP/S

OBCX/S

P100/S

V100/S

RBU/D

OFP/D

OBCX/D

P100/D

V100/D

7001001180(b)

MemoryTime

Energy Consumption (J)600100

Package200.003.50

3.00
150.002.50

sec.

Power Consumption (W)100.002.00

1160
100.001.50

1140
5001001120

400100

1.00

50.00

0.50

sec.

1100
0.00

0180

3001000160

0.00SingleDouble2001000140(b)100100

OBCX/S

P100/S

V100/SRBU/DSingle

OFP/D

OBCX/D

P100/DTime

V100/D

Double

図 31 問題サイズ「Medium」の場合の ICCG 法の計算時間と（a）消費電力（W），（b）消費エネルギー（J），各計算機環境の最適ケース

3.50
500.003.00
400.002.50

sec.

Energy Consumption (J)

OFP/S

CDU

600.000100
RBU/S

Memory

4.00

0120

0100

4.50

700.002.00

300.001.50
200.001.00
100.000.50

600.0045.00

Memory(a)

Package500.00

0.00

Time
35.00

400.00Single25.00
300.00

Double

図 33 OFP：問題サイズ「Medium」の場合の ICCG
法の計算時間と（a）消費電力（W），（b）消費エネ
ルギー（J）

30.00

sec.

Power Consumption (W)

0.00

40.00

20.00

200.0015.00

図 33，図 34 は OFP において，

10.00
100.00
5.00

0.00
RBU/S

(b)

OFP/S

OBCX/S

P100/S

V100/S

RBU/D

OFP/D

OBCX/D

P100/D

15,000.00

Memory
Time

0.00

 a:

Cache モード

45.00

 b:

Flat モード，DDR4 のみ使用

 c:

Flat モード，MCDRAM のみ使用

V100/D

Package

40.00

10,000.00

30.00

を Medium（図 33），Large（図 34）の問題サイズ

25.00

sec.

Energy Consumption (J)35.00を適用し，ICCG 法の計算時間，消費電力（W），

20.00

5,000.00

消費エネルギー（J）を比較したものである。Cache

15.00

10.00

モード（OFP/a）と Flat/MCDRAM（OFP/c）はほぼ

5.00

0.00

同じ挙動である。これらのケースでは基本的に

0.00
RBU/S

OFP/S

OBCX/S

SingleP100/S

V100/S

RBU/D

OFP/D

OBCX/D

P100/D

V100/DCPU 部分に収納されている MCDRAM が主として

Double

使用されるため，メモリ部分の電力消費はほとん

図 32 問題サイズ「Large」の場合の ICCG 法の計算時間と（a）消費電力（W），（b）消費エネルギー（J），各計算機環境の最適ケース

どない。

も 1 ノード 120W 以上）は考慮されていない。P100
20

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書(a)250.0070.00

全般的に，問題規模が大きく，SIMD 化・ベクト

Memory60.00

Timeル化が進んでおり，演算密度が高い場合，単精度

200.0050.00

150.00

演算適用による計算速度の増加は顕著である。

40.00

OFP，GPU（P100，V100）は Intel Xeon（RBU，OBCX）

sec.

Power Consumption (W)

CDU

30.00

100.00と比較して消費電力，消費エネルギー共に小さい。

20.00

また，SELL-C-は，OFP を除くと必ずしも有効

50.00
10.00

ではなく，むしろ〔星野他 2017〕によるスレッド
0.00

0.00

並列化オーバーヘッドの除去（OPT）がより効果的
SingleDouble12,000.00

MemoryCDU

70.00

であることがわかった。今後は ELL への OPT の60.00実装を実施し，性能を評価する予定である。また，

50.00

今回は GPU を使用する場合は CPU の消費電力・

40.00

消費エネルギーを評価していないが，アプリケー

Time

10,000.00

8,000.00

sec.

Energy Consumption (J)(b)

6,000.00

30.00

ション全体としての効率を評価するためには，

20.00

CPU と合わせた消費電力・消費エネルギ評価が重

10.00

要である。

4,000.00

2,000.00

0.00

0.00

Single⑨自動チューニング：ppOpen-AT における演算精

Double度を考慮した実行時間最適化機能の設計〔13,22〕

図 34 OFP：問題サイズ「Large」の場合の ICCG
法の計算時間と（a）消費電力（W），（b）消費エネ
ルギー（J）

自動チューニング（AT）言語の ppOpen-AT では，
性能パラメータを設定し，プログラムをチューニングする際に性能パラメータを変化させて性能を問題サイズが大きくなると消費電力（W）はむ

測定することが AT 適用のシナリオとして定義さ

しろ OFP/b が他と比べて小さくなる。図 31，図 32

れている。そのため，プログラムの実行時に測定でも示したようにメモリ消費電力は RBU，OBCX

結果を参照できる。

と比較すると少ない。

ppOpen-AT により生成される最適化候補におい

表 2 に示すように MCDRAM と DDR のメモリて，どの候補を用いて実行するかを内部では性能

バンド幅の比は 6:1 程度であるため，本研究で扱パラメータにとる。対象箇所の実行時間の計測結っている疎行列ソルバーのような memory-bound

果から，実際にどの候補を利用するかを決定するなアプリケーションでは，その影響は特に顕著で仕組みがプログラムとして自動生成される。

あり，OFP/c-OPT と OFP/b-OPT の計算性能の比は

ここでは，演算対象の演算精度を変化させ，ユ倍精度の場合 5 倍程度である。

ーザによりあらかじめ提示された演算精度以下と本研究では，著者等が先行研究において開発しなる条件で最高速となる演算精度の組み合わせをた有限体積法における ICCG ソルバーをターゲッ探索する最適化を，新たな AT 機能として実現すトとし，実装方法，問題規模と低精度演算によるる。

性能改善の関係に注目し，5 種類のハードウェア

具体的には，以下である：

環境下での検討を実施した。実装方法（疎行列格
納形式），問題規模，ハードウェア環境によって，



低精度演算を使用することにより，計算時間，消AT ソフトウェア開発者による要求指定：基準となる演算に対する相対的な演算精費電力（W），消費エネルギー（J）への様々な効果

度劣化の許容値（許容相対誤差）を事前に与があることがわかった。

える。
21

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

!oat$ changeString definition region start

最適化対象：
基準となる演算精度に対して，より低い精real(DP) :: drhogqv(ijdim,kdim)

度演算（もしくは，より高い精度演算）によreal(DP) :: drhogqc(ijdim,kdim)

る混合精度演算により，(A)の要求を満たし，real(DP) :: drhogqi(ijdim,kdim)

最高速となる組み合わせを選ぶ。

…

以上の新らしい AT 機能を実現する，ppOpen-AT!oat$ changeString definition region end

における拡張機能（拡張した記述方式（ディレク!oat$ changeString definition region startreal(DP) :: wkティブ））を設計した。それを以下に示す。

(wk_nmax)real(DP) :: ml_Pconv(ijdim,kdim)
real(DP) :: ml_Pconw(ijdim,kdim)!oat$ static changeString region start(From 変換前 to 変換後)…

＜対象の文＞!oat$ changeString definition region end!oat$ static changeString region end!oat$ static changeString region end(From 変換前 to 変換後)

(From DP to SP)

加えて，上記ディレクティブで囲まれた範囲内で，

名 古 屋 大 学 情 報 基 盤 セ ン タ ー 設 置 の Fujitsu
PRIMEHPC FX100 を用いて予備評価を行った結果，同じグループとして扱う変数を以下で囲む。

ある変数宣言グループをまとまりとした単精度演!oat$ changeString definition region start

算（SP）で実行した場合，全て倍精度演算（DP）＜対象の文＞

で実行した場合に対して 0.7%速度向上が得られ!oat$ changeString definition region end

る場合があることを確認した。この場合，DP による演算精度を基準とした場合の相対誤差は8.50×10-13 であった。

2019 年度は，提案 AT 機能が有効かどうかを検証することを目的に予備実験を行なった。

以上の予備実験結果より，AT ソフトウェア開発

使用するプログラムは，全球雲解像モデル

者による許容相対誤差が 1.00×10-12 以下である場

NICAM (Non-hydrostatic ICosahedral Atmospheric

合は混合精度演算で高速化の可能性があることがModel)(非静力学 正二十面体 大気モデル)の，雲の判明した。その場合，0.7%程度の速度向上を得る

微 小 な 物 理 演 算 を 行 う カ ー ネ ル

可能性があることを明らかにした。

physicskernel_microphysics に含まれるサブルーチ

なお，この単精度演算との混合精度演算を用い

ン mp_snw6 を対象とした．mp_nsw6 の変数は変数ることによる速度向上の度合いは，DP に対する低を 139 個使用しており，もともとはすべて倍精度

精度演算の高速化率が大きく影響する。そのため，として宣言されている。

半精度演算を用いるとハードウェア的に DP に対本予備実験では，倍精度（double precision）から，

してより高速化される計算機環境を用いて再評価

一部に単精度（real）を用いた混合精度演算を想定をする予定である。

する。

6. 今年度の進捗状況と今後の展望

以下に mp_snw6 を対象とした場合の，ppOpen-

本研究は，計算科学・計算機科学・数値アルゴ

AT によるコード指定例を掲載する：

リズム分野の研究者の協力のもと，様々な数値ア!oat$ static changeString region start

ルゴリズムの最新アーキテクチャに向けた最適化，(From DP to SP)

低精度・変動精度導入による高速化，消費電力節22学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

約，および精度保証に基づく最適精度選択のため開発し，悪条件問題への有効性が示されたの自動チューニング選択手法の確立と実アプリケ〔10,11〕。

ーションでの検証を目指したものである。個別の 演算精度の影響は，問題規模・疎行列格納手法，

研究開発項目の進捗状況と今後の展望については

アーキテクチャによって多様である。
 局所的に演算精度を変更する手法の開発に着5.に示した通りである。
本研究では，ステンシル計算（①構造格子，②

手し，有効性が示された。

半構造格子），疎行列演算（③一般行列，④悪条件
問題，⑤Adaptive CG（局所前処理による）），⑥H

当初計画は，2018 年度：各アルゴリズム・アプ

行列，⑦精度保証，⑧消費電力測定，⑨自動チュ

リケーション最適化，多様な演算精度適用，消費ーニング手法，の各項目についての研究開発を実

電力測定，2019 年度：精度保証手法確立，2020

施する。本研究は 2018 年度から，３年計画として

年度：自動チューニング手法確立，であった。

実施している。2018 年度は，上記①～⑩の項目に2019 年度末時点で，ほぼ当初の予定通り目標

ついて，（1）性能最適化，（2）多様な計算精度，を達成する見込みだが，2020 年度は各アルゴリ（3）精度保証，（4）実アプリケーションを対象とズムの，演算精度，最適化，アーキテクチャのした消費電力測定，を中心として予備的な検討を他，問題規模も考慮して消費電力・エネルギーへ実施し，後掲のように学会等で発表した他，ISC18，の体系的な影響評価を継続して実施する他，低精SC18，ISC19，SC19 等の国際会議の展示ブースで度・混合／変動精度演算に関する研究開発を，こも紹介した。2018 年度の段階で低精度・混合／変れまでの研究成果〔7,10,11,35〕を元に継続して

動精度演算をアルゴリズム，最適化，精度保証，実施し，自動チューニング手法確立を目指す。精消費電力まで含めて扱った研究事例はほとんどな度保証手法については，疎行列演算に加えて，H

かったが，SC19（2019 年 11 月）では，低精度・行列，偏微分方程式解法向け手法の研究開発も実

混合／変動精度演算を数値計算に取り入れた研究施する（後者については 2019 年度途中から検討事例が多数見られた。2019 年度までに得られた知を開始）。各センターの保有する NVIDIA

見をまとめると以下の通りである：

V100，東大情報基盤センターに 2019 年度末に導入された富士通 FX700 クラスタ等を使用して，

 従来，倍精度演算が適用されていた科学技術計

新アーキテクチャ向け検討も実施する。特に富士算に単精度・半精度・混合演算を適用し，反復通 FX700 クラスタを使用した，半精度演算

改良法（Iterative Refinement）併用により，同等（FP16）のフィージビリティースタディを様々の精度の計算結果がより短時間で得られる場な手法，アプリケーションについて重点的に実施合がある〔7,12,14,15,19,20,32,33,35〕。一般に，する。①～⑨は密接に連携しているが，2020 年

計算時間短縮に比例して消費エネルギー（Joule） 度の各項目概要を以下に示す：は減少する〔7, 35〕。
 悪条件問題では，低精度演算では正解が得られ

 ステンシル計算（①，②）：差分法コード，地ない場合がある。特に半精度演算は変数の範囲

震波動伝播コード Seism3D（①），全球大気計が限定されるため注意が必要であり，精度の要算コード NICAM（②）を対象として，計算時求されない反復法前処理等に適用するべきで間・消費エネルギー最小化のための演算精度ある〔12,14,15,19,20,32,33〕。

選択手法の研究開発を継続して実施する。
 疎行列演算（③，④，⑤）：

 従来の M 疎行列向け精度保証手法〔Ogita, T. etal., Computing, 2001〕は，相対誤差上限の見積も

 一般行列，悪条件行列，Adaptive CG ともりが厳しめであったため，より現実的な手法をに，反復改良法適用による低精度・混合／
23

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

変動精度演算の安定化，反復法前処理に対率的な精度保証付き計算法を開発する。

する低精度演算適用の検討，様々な手法

 三次元不均質場における楕円型偏微分方程（Pipeline 法，Dynamic Loop Scheduling，

式の解の精度保証付き数値計算を実現す

hCGA 法，SELL-C-），計算条件，問題規る。有限要素法等の離散化する際に生じる模等の各演算精度による影響について，安誤差も考慮し，2019 年度迄に開発した疎行

定性・実行時間の評価，各 CPU，GPU での列計算に対する精度保証法と組み合わせる

消費電力測定を継続して実施する。

ことにより，対象問題の離散化誤差を含めた包括的誤差評価を高速に達成する。

 計算時間・消費エネルギーを最小にするた

 電力測定（⑧），自動チューニング（⑨）：

めの演算精度選択手法の研究開発を継続して実施する。

①～⑦の各項目と連携して，アルゴリズ

 2019 年度から新たに研究項目として追加しム，実装手法，問題規模，諸パラメー

た，大規模分散並列問題で局所的・動的にタ，演算精度，アーキテクチャと消費電

演算精度を変動，動的に負荷分散を適用す

力，計算時間，計算結果の関係を体系化る手法の研究開発を継続して実施する。

し，所望の計算精度を得るための最適パ： H 行列（⑥）


ラメータ設定を自動的に決定する手法の低精度・混合／変動精度演算に関する検検討を実施し，ppOpen-HPC，h3-Open-

討を継続して実施し，悪条件問題向けのBDEC への実装を行う。

前処理手法の開発，電磁気学，地震発生サイクル等シミュレーションによる性能

7. 研究成果（発表予定も含む）

検証，精度評価，電力測定を実施する。

(1) 学術論文

また，悪条件問題を中心とした精度保証なし

手法について検討を実施する（⑦参(2) 国際会議プロシーディングス照）。

[1] K. Fujita, M. Horikoshi, T. Ichimura, L. Meadows,

 精度保証（⑦）：

K. Nakajima, M. Hori and L. Maddegedara, 疎行列演算（反復解法）については，2019

Development of Element-by-Element Kernel

年度研究成果〔10,11〕を元に，悪条件問題

Algorithms in Unstructured Implicit Low-Orderに対して反復改良法と組み合わせた精度保Finite-Element Earthquake Simulation for Many-

証手法の研究開発を実施する。

Core Wide-SIMD CPUs, Proceedings of ICCS H 行列については，悪条件問題を中心とし2019, Lecture Notes in Computer Science 11536,

た精度保証の研究開発を実施する。階層的

267-280, 2019

低ランク近似法により，密行列における行[2] Nakajima, K., Gerofi, B., Ishikawa, Y., Horikoshi,

列計算の演算量を低減することができる。

M., Parallel Multigrid Methods on Manycore

低精度演算を用いる場合には H 行列のランClusters with IHK/McKernel, IEEE Proceedings

クを小さくすることができ，大幅に少ない

of 10th Workshop on Latest Advances in Scalable

演算量で同等の計算精度を実現できる。H

Algorithms

forLarge-Scale

Systems

行列のランクが特異値の減衰率に依存する

conjunction

with

SC19

International

ことに着目し，様々なアプリケーションでConference for High Performance Computing,

用いられる行列について特異値の観点か

Networking, Storage, and Analysis), 2019

ら、低精度演算と低ランク近似の適用範囲[3] Yamaguchi, T., Fujita, K., Ichimura, T., Naruse, A.,を明らかにし，密行列直接解法における効

Lalith, M., Hori, M., GPU implementation of a
24

(The

in

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

sophisticated implicit low-order finite element(3) 国際会議発表

solver with FP21-32-64 computation using[9] Hosokawa, H., Okuda, H., Mixed Precision

OpenACC, Proceedings of Sixth Workshop on

Conjugate

Accelerator

Directives

Hierarchical Memory Configurations, Proceedings(WACCPD) in conjunction with SC19 (Theof 15th U.S. National Congress on Computational

International Conference for High Performance

Mechanics (USNCCM15) (Austin, Texas, July 28-Computing, Networking, Storage, and Analysis),

August 1, Austin, 2019)

Programming

Using

Gradient

Iterations

Considering[10] Ogita, T., Nakajima, K., Accurate and verified

2019

solutions of large sparse linear systems arising[4] Ichimura, T., Fujita, K., Yamaguchi, T., Naruse, A.
et al., 416-PFLOPS fast scalable implicit solver on

from

low-ordered

Conference

unstructured

finiteelements

3D

Poissononequation,

MatrixInternational

Analysis

andits

accelerated by 1.10-ExaFLOPS kernel with

Applications (MAT TRIAD 2019), Liblice, Czech

reformulated AI-like algorithm: For equation-

Republic, September 10, 2019)

based earthquake modeling, Research Poster for[11] Ogita, T., Nakajima, K., Verified solutions of large

SC19 (The International Conference for Highsparse linear systems arising from 3D Poisson

Performance Computing, Networking, Storage,

equationand Analysis), 2019

Numerical

in

HPC

environments,

Mathematics

EuropeanandAdvanced[5] Ma, Q., Yokota, R., Runtime System for GPU-

Applications Conference 2019 (EnuMath 2019)

based Hierarchical LU factorization, Research(Egmond aan Zee, The Netherlands, October 2,

Poster for SC19 (The International Conference for

2019)

High

Performance

Computing,

Networking,[12] Nakajima, K., Innovative Methods for ScientificStorage, and Analysis), 2019

Computing in the Exascale Era by Integrations of[6] Ootomo, H., Yokota, R., TSQR on TensorCores,(Simulation+Data+

Learning)

Research Poster for SC19 (The International

Supercomputing“Society

Conference for High Performance Computing,

Universritaria

Networking, Storage, and Analysis), 2019 (BestSupercmoputo en la Tranformacion Digital)

Poster Candidate)(Guadalajara, Mexico) (Keynote Talk)

inde(S+D+L):
5.0”,

Supercomputo

Jornada2019(El[7] Sakamoto, R., Kondo, M., Fujita, K., Ichimura, T.,[13] Katagiri,

Nakajima, K., The Effectiveness of Low-Precision

Exascale

Floating Arithmetic on Numerical Codes: A Case

International Symposium on Computing andStudy on Power Consumption, ACM Proceedings

Networking), in 7th International Workshop on

of

High

Computer Systems and Architectures (CSA'19)

Performance Computing in Asia-Pacific Region,(Nagasaki, November 26-29, 2019) (Keynote(HPC Asia 2020), 2020

Talk)

the

International

Conferenceon[8] Nakajima, K., Parallel Multigrid Method on
Multicore/Manycore

Clusters,

IXPUG[14] Nakajima,T., Towards Auto-tuning Technology in
Era,

CANDAR'19

K., AnInnovative(The

SeventhMethodfor(IntelIntegration of Simulation /Data/Learning in the

Extreme Performance Users Group) HPC Asia

Exascale/Post-Moore Era, APCOM 2019: Asian2020, 2020

Pacific Congress on Computational Mechanics(December 18-21, 2020, Taipei, Taiwan) (SemiPlenary Talk)
25

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書[15] Nakajima, K., Parallel Multigrid Methods with[22] 片桐孝洋，櫻井刀麻，ポストムーア時代に向Adaptive Multilevel hCGA on Manycore Clusters,

けた自動チューニングに向けて～スレッド数

MS1603:

ParallelModels,の動的最適化，第24回日本計算工学会講演会

Algorithms

andExtreme（大宮，2019年5月30日）

Programming

Frameworks

forComputing & Big Data, APCOM 2019: Asian[23] 中島研吾，高性能・変動精度・高信頼性数値

Pacific Congress on Computational Mechanics

解析手法とその応用，第24回日本計算工学会(December 18-21, 2020, Taipei, Taiwan)

講演会（大宮，2019年5月30日）[16] Ozaki, K., Upper bound of maximum norm of[24] 伊田明弘，星野哲也，メニーコアクラスタに

inverse matrix in test set and its application,

おける階層型行列法の高速化に向けた性能評

MS1603:

ParallelModels,

価，第24回日本計算工学会講演会（大宮，2019

Algorithms

andExtreme

年5月30日）Programming

Frameworks

forComputing & Big Data, APCOM 2019: Asian[25] 森田直樹，橋本学，奥田洋司，深層学習によPacific Congress on Computational Mechanics

る反復法線形方程式ソルバの計算時間推定(December 18-21, 2020, Taipei, Taiwan)―学習モデルの提案と適用性評価―，第24回[17] Deshmukh, S., Yokota, R., Distributed Memory

日本計算工学会講演会（大宮，2019年5月31日）

Task-Based Block Low Rank Direct Solver, HPC[26] 深谷猛, グドール聖哉, 張臨傑, 岩下武史，倍Asia 2020 (poster) (Fukuoka, January 2020)

精度と単精度を用いた混合精度GMRES(m)[18] Apriansyah, M.R., Yokota, R., QR Decomposition

法の性能評価，第48回数値解析シンポジウム（福井，2019年6月12日）

of Block Low-Rank Matrices, HPC Asia 2020[27] 大 友 広 幸 , 横 田 理 央 ， Tensor コ ア を 用 い た(poster) (Fukuoka, January 2020)[19] Iwashita, T, Nakajima, K., Shimokawabe, T.,

TSQR の GPU 実 装 , 情 報 処 理 学 会 研 究 報 告Nagao, H., Ogita, T., Katagiri, T., Yashiro, H.,（2019-HPC-170-38）（北見，2019年7月24日）

Innovative[28] 堀越将司，中島研吾，Balazs Gerofi，石川裕，Software Platform for Scientific Computing in the

Parallel Preconditioned Itertative Solvers on

Exascale Era by Integrations of (Simulation + Data

Oakforest-PACS，2019年並列／分散／協調処+ Learning) , HPC Asia 2020 (poster) (Fukuoka,

理に関する『北見』サマー・ワークショップ

January 2020)（SWoPP北見2018），日本応用数理学会「行列・Matsuba,H.,

h3-Open-BDEC:

固有値問題の解法とその応用」研究部会[20] Nakajima, K., Innovative Methods for Scientific（MEPA）（北見，2019年7月25日）

Computing in the Exascale Era by Integrations of(Simulation+Data+ Learning), SIAM Conference[29] 深谷猛，グドール聖哉，張 臨傑，岩下武史，

on Parallel Processing for Scientific Computing

倍精度と単精度を用いた混合精度GMRES(m)(PP20), MS25/36: Progress and Challenges in法の収束性に関する実験的評価，日本応用数

Extreme Scale Computing and Big Data (Seattle,

理学会2019年度年会（東京，2019年9月3日）[30] 大 友 広 幸 , 横 田 理 央 , Tensor コ ア を 用 い た

WA, USA, February 12-15, 2020)

TSQR, 日本応用数理学会 2019 年度年会（東(4) 国内会議発表

京，2019 年 9 月 5 日）

[21] 細川洋輝, 奥田洋司, 階層的メモリ構造を考[31] 田中一成，楕円型方程式の弱解に対する正値

慮した反復法ソルバーの混合精度計算，第24性証明法, 日本応用数理学会 2019 年度年会回日本計算工学会講演会（大宮，2019年5月29（東京，2019 年 9 月 5 日）
[32] 中島研吾，塙 敏博，伊田明弘，下川辺隆史，日）
26

学際大規模情報基盤共同利用・共同研究拠点 2019 年度共同研究 最終報告書

三木洋平，星野哲也，有間英志，田浦健次朗，

日本応用数理学会 2019 年度連合発表会（東京，

工藤知宏，関谷勇司，中村 遼，Society 5.0 実

2020 年 3 月 5 日）

現に向けた（計算＋データ＋学習）融合，[35] 中島研吾，坂本龍一，星野哲也，有間英志，塙

AXIES 2019 福岡（福岡，2019 年 12 月 13 日）

敏博，近藤正章，低精度演算とアプリケーシ[33] 中島研吾，岩下武史，八代尚，下川辺隆史，長ョン性能，情報処理学会研究報告（2020-HPC-

尾大道，荻田武史，片桐孝洋，松葉浩也，（計

174-5）（第 174 回 HPC 研究会）（オンライ

算＋データ＋学習）融合によるエクサスケー

ン，2020 年 5 月 13 日）（in press）

ル時代の革新的シミュレーション手法，第11回 自動チューニング技術の現状と応用に関(5) その他（特許，プレス発表，著書等）

するシンポジウム（ATTA2019）（東京，2019年[36] FP21AXPY

12月22日）OpenACC

sourcehttps://github.com/y-mag-chi/fp21axpy[34] 田中一成, 中尾充宏, A priori error estimates for
Poisson's equation with discontinuous coefficients,

27code,