学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

jh150006-NA05

コデザインアプローチによる高性能電磁場解析基盤の確立
岩下武史（北海道大学）

概要 本研究では，計算機科学研究者と電磁場解析分野の研究者の協働により，次世代の電磁場解析の基盤となる技術の確立を目指した研究活動を行った．具体的には，課題
１：電動モータを対象とした電磁場解析，課題２：マイクロマグネティクス計算，課題３：複数電磁場解析実行フレームワークの 3 つを主な研究対象とした．課題１については，多相化した TP-EEC 法と並列処理を組み合わせた高速化技術を確立し，その有効性を示した．また，有限要素解析において必要となる線形ソルバの高速化について，電磁場解析向けの新しい前処理手法の理論的解明を進めた．課題２については，H 行列法と呼ばれる行列近似手法をマイクロマグネティクス計算に導入することを行った．課題３については，フレームワークの基盤となるジョブ並列スクリプト言語 Xcyrpt を北海道大学情報基盤センターの Hitachi SR16000 に移植し，その概念設計を実施した．

1.

共同研究に関する情報

析に対する性能向上の要求も極めて大きい状況にある．また解析の精度面においても，機器の効率(1) 共同研究を実施した拠点名に大きな影響を与える鉄心材料の高精度なモデリ

北海道大学，京都大学

ング手法が求められている．その他にも実応用設計の現場では，設計の初期段階において，多数の(2) 共同研究分野■

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



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



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



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

小・中規模モデル（数万～百万自由度程度）をなるべく短時間に解きたいという要望がある．こうした技術的要請に応えるために，本研究では主に以下の 3 つの研究課題を掲げ，その実現を研究目的とする．
研究課題１ 電動モータを対象とした非線形電(3) 参加研究者の役割分担

磁場解析の高速化

岩下武史 北海道大学情報基盤センター

研究課題２ マイクロマグネティクスに基づく－研究統括，高性能計算に関する助言

高精度な鉄損評価法の確立

美舩健 京都大学大学院工学研究科

研究課題３ 並列計算機上における複数電磁場－アルゴリズム開発

解析の効率的な実行法の確立

高橋康人 同志社大学理工学部

これらの研究課題に対し，計算機科学研究者と－電磁場解析プログラム開発

電磁場解析分野の研究者の協働により，高性能計

平石拓 京都大学学術情報メディアセンター

算（HPC）技術と解法技術の両面からその解決に－複数電磁場解析システム開発取り組んでいく．

2.

研究の目的と意義
3. 当拠点公募型共同研究として実施した意義

電磁場解析は様々な電気機器，電子デバイスの本研究の申請段階において，我々は次世代の高設計や評価に広く用いられ，実応用上の設計作業性能電磁場解析基盤の確立を目的とした場合，単では欠くことのできないものとなっている．近年，に電磁場解析プログラムの開発者が既存のライブ省エネルギー性能等の機器性能に関する要求は高ラリや解法を組み合わせ，さらに並列化するといまる一方であり，機器設計の基盤である電磁場解
1学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

った従来型のアプローチでは不十分であるとの考えを持っていた．即ち，次世代の電磁場解析手法の確立には，最新の解法や実装方法と応用分野固有の知識の相互活用が必須であり，計算機科学者と計算科学者の協働が不可欠との認識を持っていた．
本研究を当拠点公募型共同研究として実施した(a) Mesh

結果，上記の申請時に想定していた計算機科学者と計算科学者の協働が円滑に行われ，研究推進上

1200
Line voltage [V]の大きな意義があった．例えば，電動モータの解析では，TP-EEC 法という新しい電磁場解析手法と並列処理に関する知見を組み合わせることにより，電磁場解析の高性能化に寄与することができ

0-1200

た．また，マイクロマグネティクス解析では，高0性能計算分野において開発された分散 H 行列法の本研究の成果の多くは電磁場解析に関する応用研

TP-EEC1500

U-phase current [A]れたものであり，当拠点公募型共同研究として実施した意義は大きいと言える．

Steady-state (TDPFEM)
Parallel TP-EEC (256 processes)

1000
500
0-500-1000
Parallel TP-EEC (16 processes)-1500

4. 前年度までに得られた研究成果の概要

120 180 240 300 360
Electric angle [deg](b) Voltage waveform

電磁場解析分野への導入に成功した．このように，究者と計算機科学者の共働の結果として産み出さ

600

該当なし

180540720360Electric angle [deg]9001080(c) U-phase current waveform
図15. 今年度の研究成果の詳細

PWM インバータ駆動 IPM モータの解析結果

図 1 に示す IPM モータを対象に非線形渦電流解① 電動モータ解析の高速化（研究課題 1）

析を行い，実機解析における有効性を検証した．

本研究では，時間分割型の並列有限要素法とし

図 1(c)に，U 相電流波形を示す．各プロセスで独て並列化 TP-EEC 法を開発した．あらゆる時間周立に過渡解析を実行しているため初期ステップで期波形は，時間ステップの並べ替え操作により任は誤差が大きいが，並列化 TP-EEC 法による補正

意の多相交流時間周期性を有する問題に変換できを繰り返すことで最終的に定常状態が得られている（上記操作を相変換と呼ぶ）．分担者の高橋らのることが確認できる．

提案手法では，各プロセス（またはスレッド）に表 1 に，並列化 TP-EEC 法，従来提案されていそれぞれ初期値を与え，過渡解析を独立に実行すた時間分割型並列有限要素法である並列化

る．相変換により多相交流時間周期問題と解釈す

TPFEM，および定常解高速求解法である（逐次）ることで，各プロセスで独立に得られた過渡解にTP-EEC 法を用いた過渡解析において 1 周期の定多相交流 TP-EEC 法を適用し，収束の遅い誤差成分やプロセス間に生じる誤差を抽出し，補正する．

常解を得るまでに要した計算時間を示す．なお，妥当な計算時間の比較を行うために，まずこの補正を繰り返すことで，最終的に定常解を算

TPFEM について定められた相対残差以下の精度出することができる．
2が得られるまで解析を行い，得られた解析結果か学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

表1セス数を 256 とした解析において，誤差の補正にPWM インバータ駆動 IPM モータの計算時間
Calculation time [s]Number oftime steps
in oneperiod

1024

2048
4096

Number ofprocesses

Parallel
TPEEC

Parallel
TPFEM

163264128
256
TP-EEC(sequential)
256
256

2492.31871.21523.71714.22446.911736.1
5937.72970.91500.5753.4

関する計算が全体の計算時間において支配的となっていた．そのため，1 周期あたりの時間ステップ数を 2048 や 4096 に増加させた場合，過渡解析に要する計算時間は増大するものの総計算時間はそれほど大きく増加しない結果となった．従って，
1 周期あたりの時間ステップ数を 4096 とした解19915.7
2802.92998.3

析のようにプロセス数に対して時間ステップ数が2484.38273.2

相対的に大きい場合，並列化 TP-EEC 法の使用がより有効となると考えられる．

らトルクを算出した．その後，TP-EEC 法（逐次／並列）による計算では，トルクをあるタイムステップ間隔で算出し，この値が TPFEM における② マイクロマグネティックス計算の高性能化（研究課題 2）

トルク値と十分な精度で一致するまで時間発展のマイクロマグネティックス計算は，強磁性体の計算を行った．
使用計算機は北海道大学のスーパーコンピュー

磁区構造解析手法として，これまで主に磁気記録

タ Hitachi SR16000 であり，Flat-MPI 方式で 1の分野で物理現象の解釈やデバイスの設計に用いコアに 1 プロセスを割り当てて数値実験を行った． られてきたが，高精度な磁性材料モデリング手法としてもモータ内の鉄損評価等の様々な応用上の表 1 に示す通り，並列化 TPFEM は高並列計算環展開が期待されている手法である．

境下において優れた台数効果を示し，256 プロセス使用時に最も高速な計算が実行可能である．一マイクロマグネティックスでは，外部磁界によ方，並列化 TP-EEC 法の場合，64 プロセス使用るエネルギー，磁気モーメント同士の相互作用で時が最速であり，並列化 TPFEM において 128 プある反磁界によるエネルギー，隣り合うスピン間ロセスを使用した場合と同等の性能を示している． の相互作用である交換磁気エネルギー，磁気異方性に起因する磁気異方性エネルギーの総和が極小
また比較のために，1 周期を 2048 分割および4096 分割した場合の過渡解析も実行した．表 1 に，256 プロセス使用時の計算時間を示す．分割数がになるよう磁化の分布を決定する．この時，磁化M の歳差運動を表す方程式が，以下の増加するにつれて，並列化 TPFEM の計算時間はLandau-Lifshitz-Gilbert（LLG）方程式である．

急激に増加するのに対し，並列化 TP-EEC 法の計算時間は緩やかに増加しており，1 周期あたり(1 + α ) ∂∂Mt = −γ (M × H ) − αMγ M × (M × H )4096 ステップの場合，並列化 TP-EEC 法の方がM s は飽和磁化の大きさ，H eff は実効磁界，
ここで，並列化 TPFEM と比べ 2 倍以上高速な解析を実現

α は減衰定数，γ はジャイロ磁気定数である．実2eff(1)

効磁界 H eff は，外部磁界 H E，反磁界 H D，磁気異している．本現象は以下のように説明できる．

方性エネルギーに起因した磁界 H K，および交換

並列 TP-EEC 法は，各プロセスで独立に実行す磁気エネルギーに起因した磁界 H Aの和として表る過渡解析と多相交流 TP-EEC 法によるプロセスされる．ここで，反磁界 H Dは間の誤差の補正の 2 つの計算部分から成り立っている．前者の計算量は各プロセスに割り当てられN 
1 ⌠(r − r ' )dV 'H D (r ) = ∑  ( − ∇'⋅M (r ' ) )| r − r ' |3l =1  4πµ 0 ⌡ Vl(2)1(r − r ' )+d S ' M (r ' ) ⋅ n
4πµ 0 ∫ ∂Vl| r − r ' |3

た時間ステップ数に比例するのに対し，後者の計算量は 1 周期あたりの時間ステップ数には無関係で，使用プロセスの総数に比例する．本解析の場
合，1 周期あたりの時間ステップ数を 1024，プロeffs3

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月の式により求められる．ここで，V と∂V は磁化の析分野では，連立一次方程式の求解法（線形ソル定義された要素とその境界を表す．磁化ベクトルバ）として，前処理付き反復法を用いることが通を要素内で一定とすれば，(2)式右辺第一項は 0 と例となっており，特に ICCG（不完全コレスキーなる．

分解前処理付き共役勾配）法が事実上の標準解法静磁界相互作用（反磁界）計算ではすべての要となっている．しかしながら，ICCG 法を用いた素からの寄与を考慮する必要があるため，演算量場合，問題サイズの拡大に伴う反復回数の増加はは要素数 N の 2 乗に比例し，全計算時間の大部分避けられないため，マルチグリッド法等の大規模を占める．そこで本研究では，反磁界計算に課題問題に有効な線形ソルバに関する研究が行われて代表者が別のプロジェクト（JST CREST）で開発いる．このような背景の下で，美舩らが提案したしている分散 H 行列ライブラリ HACApK を導入

幾何的ブロック対角前処理は，有限要素法を用いし，演算量の大幅な削減を試みた．H 行列法は，た大規模電磁界計算に伴う計算負荷を削減するこ(2)式の右辺項に見られるような，距離の 2 乗や 3とを目的として開発された計算手法である．以下乗に反比例する相互作用を表現する行列に対してに，幾何的ブロック対角前処理に関する研究につ適用可能である．このような場合，互いに遠方にいて述べる．

ある二つの要素間の相互作用は近接する要素間の電動モータの解析において標準解法となってい相互作用と比べて相対的に小さいため，相互作用

る有限要素法によって 3 次元電磁界計算を行う場の表現行列の一部をよりデータ量の少ない低ラン合，電磁界分布を表現するために用いる形状関数

ク行列で近似しても，例えば行列ノルムにおいてとして，辺要素と呼ばれるベクトル関数を使用す大きな変化はない．H 行列法では本性質を用いるるのが一般的となっている．辺要素を使用する場ことで，行列の一部を低ランク化し，演算量，メ合，有限要素法における未知変数は格子分割の各モリ使用量を低減させる．本研究では，マイクロ辺に割り当てられる．有限要素法によって電磁界マグネティックス計算において，分散 H 行列ライの基礎方程式を離散化することにより，これらのブラリ HACApK を用いることにより，H 行列法未知数が満たすべき連立一次方程式が導出される．

を導入することに成功した．
（なお，本導入におい以下では，有限要素法で使用する計算格子として最も問題となったのはマイクロマグネティクスて，直方体のみを用いた規則的な分割を想定する．

計算のプログラムが C 言語ベースであったのに対格子分割を規則的なものに限定する場合，解析対して HACApK が Fortran ベースであった点であ象の形状表現に関する柔軟性が失われるが，一方った．HACApK は通常の数値計算ライブラリとはで，格子作成にかかるコストを劇的に削減できる異なり，H 行列生成のために内部でユーザ関数をだけでなく，その規則性を積極的に利用した高速呼び出す．この場合，メインのマイクロマグネテ計算手法を使用できる可能性が生じる．近年ではィクス計算からみて C → Fortran → C のよう規則的な格子分割を用いる有限要素解析の有用性な形で関数を順次呼びだす必要がある．アプリケに期待が高まっており，形状表現に関する柔軟性ーション固有の調整に加えて，このような特性への欠落についてその克服を目的とした研究も広くの対応に時間を要した．
）

為されている．
計算格子の各辺が x, y, z 方向のいずれかに平行である場合，有限要素法により導出される連立一③ 電磁場解析の高速化（新しい前処理手法に関する研究，研究課題 1 に関連）

次方程式は次の形に表すことができる．

電動モータの解析では，有限要素解析が用いらK
 xx
 K yx K zx

れることが多いが，その場合，計算時間の大部分は連立一次方程式の求解に費やされる．電磁場解
4K xy
K yy
K zy

K xz  ux   bx    K yz  u y  =  by K zz  uz   bz (3)学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

ここで u x , u y , u z はそれぞれ x, y, z 方向に平行な[E ]ij = ∫ N i ⋅ [∇(∇ ⋅ N j )]dV

辺に割り当てられた未知数を成分とするベクトル(8)であり，b x , b y , b z は右辺に関する同様のベクトルと与えられ，行列 D, E と式(3)の係数行列の間にである．また K xx などは未知ベクトル，右辺ベクは，トルを式(3)の形に表したときの係数行列の部分行列を表している．
式(3)のような連立一次方程式を前処理付き反復法により解く場合，前処理行列 M が
 K xxM ≈  K yx

 K zx

K xy
K yy
K zy

K xz 
K yz 
K zz (4)OO   K xx
 O ≈ O
M zz   OO O K zz (9) O
E ≈  K yx

 K zx

K xy

K xz 
K yz 
O (10)

OK zy

ここで，式(3)に対して前処理行列を
M =D

これに対して，幾何的ブロック対角前処理ではOM yyO
K yyOの関係が成り立つ．

を満たすように M を設計するのが一般的である．

 M xx
M = O O

 K xxD= O OO
K yyOO O K zz (11)としたとき，境界条件に関する適当な仮定の下に式(7), (8)を用いて，行列 D -1E の条件数µが少なく(5)ともµ = O(log n )を満たすように M を設計する．式(5)中の O は零(12)を満たすことを示せることが本研究により明らか行列を表す．
現在までに，3 次元電磁界解析の例題について，となった．ここで n は係数行列の次元である．こ式(5)に基づいて IC (Incomplete Cholesky)分解あれは，幾何的ブロック対角前処理を用いる反復解

るいは代数マルチグリッド法により構成した前処法が収束までに要する反復回数が，少なくとも特理が，通常の前処理と比較して優れた性能を示す

別なケースについては問題サイズの増大に対してことを確認している．しかしながら式(5)を用いる緩やかにしか増えないことを示している．

ことの妥当性に関する理論については，不明な部④ 複数電磁場解析の高効率化（研究課題 3）分が多く検討が十分になされていない．そこで本
研究では，幾何的ブロック対角前処理によって前並列計算機上における複数電磁場解析の効率的処理された係数行列の固有値等について検討を進な実行法の確立について，これまでに「多数の電め，その理論的妥当性を明らかにすることを行っ磁場解析シミュレーション実行を効率化するフレた．

ームワーク」の具体的な要求要件の定義と基本設

透磁率が解析領域中で均一な場合，静磁界解析計およびプロトタイプ実装を行った．また同フレにおける式(3)の係数行列の成分は辺要素形状関ームワークの実装に用いるジョブ並列スクリプト数 N i を用いて

言 語 Xcrypt の 北 海 道 大 学 情 報 基 盤 セ ン タ ー[K ]ij = ∫ N i ⋅ [∇ × (∇ × N j )]dV
= [D ]ij + [E ]ij

Hitachi SR16000 環境への移植を行った．
(6)

これまでの設計に基づくフレームワークの概念
図を図 2 に示す．実行する電磁場解析プログラムと表される．ただし境界条件に由来する項についは，任意のプロセス数・スレッド数で実行可能なては省略している．ここで行列 D, E の成分は超関ようにハイブリッド並列化されており，また，い数の意味でくつかのチューニング可能な性能パラメータが存[D]ij = − ∫ N i ⋅ [∇ 2 N j ]dV在するものとする．そのようなプログラムをスパ(7)コンや並列計算サーバ環境で，それぞれ異なる入5

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

FDTD カーネルに適用し，上記のパラメータ値のXcryptライブラリ （提案システム）
解析実行順序，並列数，性能パラメータ等の制御入力データ最適化に利用できることを確認した．

・・・
起動

ユーザスクリプト

電磁場解析プログラム

入力集合や実行条件の定義

また，Xcrypt は複数のスパコン環境を同時に利ジョブ投入・実行

job

プロファイル結果反映

job
job

job

job

job

用するジョブ並列実行を 1 つのスクリプトで容易job
jobに記述できる仕組みを提供している．そこで，開・・・

性能プロファイル

・・・

発するフレームワークも図 2 のように複数のスパ解析結果・・・

試行データ追加

コンの同時利用に対応させることも視野に入れて

図 2 複数電磁場解析実行フレームワークいる．
なお，Xcrypt の SR16000 環境への移植につい力（解析空間サイズや解析タイムステップ数等）ては，Xcrypt 自体がジョブ投入コマンドやジョブを与えられて多数実行する状況を考え，そのようスクリプトファイルの書式などの環境依存の部分な実行をジョブ並列実行によりなるべく早く終わを定義した設定ファイルを作成することにより新らせることを要求要件とした．

たなスパコン環境での実行が可能になる仕組みを

本要求要件を実現するためには，それぞれのジ提供していたため，基本的にはその方式に従ってョブ実行がどれだけの計算資源を割り当てた場合実装することができた．にどの程度の時間で完了するかを予測し，それに6.

基づいて実行順序や各実行の並列数を調整する必要がある．たとえば，l 台の計算ノードが利用可能

今年度の進捗状況と今後の展望① 電動モータ解析の高速化（研究課題 1）

な環境で強スケール性に限界がある 2 つのシミュ本研究にて開発した並列 TP-EEC 法をさまざまレーション実行を行う場合，両者の負荷がほぼ等な電気機器解析に適用し，他の時間分割型並列有しければそれぞれに l / 2 ノードを割り当てて同時限要素法との得失比較や複雑な物理現象を有するに実行開始したほうがよいが，負荷に偏りがあれ実機解析における有用性検証を進めていく予定でば l ノードを割り当てた実行を順に行ったほうがある．また，現在の磁化状態が過去の遷移状態によいと考えられる．さらに，プログラムの性能パ依存する磁気ヒステリシス特性を考慮した磁界解ラメータ値や MPI プロセス数，OpenMP スレッ析の高速化を目指し，並列 TP-EEC 法のさらなるド数等を実行環境や入力サイズごとに最適になる高度化を検討する．

ように調整する必要もある．
このような調整を適切に行うため，ジョブの実② マイクロマグネティックス計算の高性能化

行ごとに性能プロファイルを取得・保存し，その（研究課題 2）履歴に基づいて後続の実行の所要時間の予測およ

HACApK を導入したマイクロマグネティック

び実行順序や性能パラメータの調整を行うシステ

スシミュレーションに対しては，高並列計算環境ムを Xcrypt を用いて開発した．プロトタイプ実装での性能評価の実施に至らなかったため，数 10～

では，与えられたプログラムを様々な入力サイズ，性能パラメータ値，MPI プロセス数，OpenMP ス

100 ノード程度の並列実験による性能評価を行うとともに，実規模問題への適用とその有用性検証

レッド数で実行し，それぞれの実行で得られた性を実施する予定である．

能値をデータベースに保存する機能を実現した．
また，その結果を取り出して利用することでプロ③ 電磁場解析の高速化（新しい前処理手法にグラムの性能チューニングを行えるようにした．

関する研究）

このプロトタイプのシステムを，ハイブリッド並

幾何的ブロック対角前処理の妥当性について考

列実行に対応した時空間タイリング付き三次元6

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

察し，収束までの反復回数に関する理論的結果をNakashima, “Dynamic Load Balancing for導出したが，解析領域中の透磁率を均一と仮定す

Parallel Computation of Hierarchical Matrices,”

る等その適用範囲は限定的である．今後はさらにISC 2015 HPC in Asia Posters, (2015.7)

一般化された理論的結果の導出に取り組む予定で
A. Ida, T. Iwashita, “HACApK: Library forある．

hierarchicalmatrices

withadaptivecross④ 複数電磁場解析の高効率化（研究課題 3）

approximation,”

研究課題 3 については，プロトタイプ実装としWorkshops on Advances in Computational

the3rdInternational

Mechanics (IWACOM-III) (2015.10)て，プロファイル情報の取得・保存・取り出し機
能や，それを用いた単一ジョブ実行レベルの最適化は実現したが，ジョブの実行順序や各ジョブに(4) 国内会議発表

割り当てる計算資源量を制御することによる複数「非線形磁気特性
髙橋康人 ，美舩健 ，藤原耕二：

ジョブ実行全体の最適化の実現は今後の課題としを考慮した時間調和渦電流解析へのニュートン・て残されている．また，本研究における他課題のラフソン法の適用に関する一考察」，電気学会静止成果として得られた解析プログラムを用いて，よ器・回転機合同研究会資料，SA-15-97/RM-15-135

り実践的な複数電磁場解析実行最適化のために本（2015.9）

システムを利用し，その有効性の評価を行いたい「時間分割型並列有限要素法を用いた髙橋康人 ：と考えている．

回転機の磁界解析」，電力エネルギーフォーラム7.

第 25 回電磁界数値解析に関するセミナー，pp.

研究成果リスト

65-79（2015.12）(1) 学術論文
Y. Takahashi, K. Fujiwara, T. Iwashita, H.
Nakashima, “Parallel Finite-Element Analysis

髙橋康人 ，北尾純士，藤原耕二，阿波根明，松尾

ofDomain

哲司：「並列化 TP-EEC 法を用いたヒステリシス磁Nonconforming

界解析の高速化」電気学会静止器・回転機合同研

Rotating

Machines

Decomposition

Based

Consideringon

究会資料，SA-16-16/RM-16-16（2016.1）

Mesh Connection,” IEEE Trans. Magn., vol. 52,
no. 3, 7401604 (2016.3)

A. Ida, T. Iwashita, T. Mifune, Y. Takahashi,(2) 国際会議プロシーディングス“Variable Preconditioning of Krylov Subspace
Methodsfor(3) 国際会議発表

AdaptiveCross

HierarchicalMatrices

Approximation,”

AnnualMeeting on Advanced Computing System andR. Plasser, Y. Takahashi, G. Koczka, O. Biro ,Infrastructure (ACSI2016) (2016.1)“Comparison of Various Methods for the Finite
Element Analysis of Nonlinear 3D periodic(5) その他（特許，プレス発表，著書等）

Eddy-current Problems,” 20th InternationalConference

ontheComputationof

Electromagnetic Fields (COMPUMAG 2015),
PB1-17 (2015.7)
T. Hiraishi, K. Munakata, A. Ida, T. Iwashita, H.

with7