日本ソフトウェア科学会第 33 回大会 (2016 年度) 講演論文集

並列有限体積法のための領域特化言語の C++への埋め込み
伊藤 正勝

宮島 敬明

藤田 直行

有限体積法のような科学技術計算では、プログラムを高速化するために、プログラミングのしやすさが犠牲にされがちである。これに対して、我々は有限体積法のための領域特化言語を C++に埋め込むことで、プログラムの記述性と実効速度の両立を図っている。我々の領域特化言語は、そのソースコードを有限体積法の記法に近づけることで、プログラムの記述性を高め、そのソースコードを OpenMP で自動並列化することで、同時に実行速度も向上させる。この領域特化言語で熱伝導方程式ソルバをプログラミングして、その実行速度を計測した結果、手動並列化されたプログラムに劣らないことが確認された。

In the development of application programs with finite volume method, the expressiveness of source program is sacrificed for improving the speed of calculation, which makes the programming time-cosuming. To
resolve this problem, we are implementing an domain-specific language embedded in C++ for parallelizing
finite volume method. A prototype of our domain-specific language can automatically parallelize an simple
heat-equation solver program so that it can run as fast as another corresponding program which is manuallyparallelized.

は、行列演算をプログラミングするための領域特化言

1 はじめに語 (Domain-Specific Language(DSL)) としてソース

科学技術計算プログラミングには、ソースコードのコードの記述性を高めつつ、並列スケルトン [2] によ

記述性とプログラムの実行速度が両立しないというる自動並列化機能によってプログラムを高速化するこジレンマがある。そして、多くの場合、実行速度のほとができる。

うが優先されて、記述性が犠牲にされるため、プログこのように、数式テンプレート技法が拡張され、ラミングには手間がかかることになる。こうした犠牲

C++ライブラリに領域特化言語としての機能を持つは、有限体積法のように行列やベクトルの演算処理ものが現れたのは、数式テンプレートの変換をコンに多くの計算時間が必要となるようなアプリケーショパイラによる中間コードの変換のように利用できるン領域では、特に顕著である。

からである。しかも、最近では、このような数式テン

しかし、行列やベクトルのコーディングにおけるプレート技法の詳細を理解しなくとも、Boost.Proto

ジレンマが動機となって数式テンプレート技法 [14]ライブラリ [9] を使うことで、数式テンプレートを[13] が開発され、これまでに様々な数式テンプレー中間コードとする領域特化言語を作成できるようにトライブラリのなかで、その技法は拡張されてきたなった。実際、NT2 は Boost.Proto を利用して開発

2[6]。なかでも、Numerial Template Toolbox(NT ) [3]A domain specific language embedded in C++ for parallelizing finite volume method
Masakatsu Ito, Takaaki Miyajima and Naoyuki Fujita, 宇宙航空研究開発機構 航空技術部門 数値解析技術研究ユニット, Numerical Simulation Research Unit,
Japan Aerospace Exploration Agency(JAXA).

された。
このような領域特化言語は、そのコードが C++と同じ構文規則に従うが、異なる意味を持ちうる [8]。
この点で、埋め込み型の領域特化言語 (EmbeddedDomain-Specific Language(EDSL)) の一種であると言える。

このように、EDSL としての C++テンプレートラ分布は無視して、温度分布は T (x) で表現する。周囲イブラリが、様々なアプリケーション領域で開発さの空気の温度を T∞ とすると、この問題は、微分方れつつある。科学技術計算においては有限要素法の程式[11] [10] [7] [1] が存在しており、実行速度を犠牲にする

d2− s) T (x) = −s T∞(1)dx2によって記述される。左辺の第１項は熱拡散、残りことなく、ソースコードを有限要素法の記法に近づけの項は周囲の空気の対流による冷却あるいは加熱速ることで記述性を高めている。

度を表す。ここで、k は熱伝導係数、s は対流の熱輸分野で、feel++を始めとする複数の EDSL 実装 [12]ただ、このように記述性と実効速度を両立させるた(k送に関わる係数である。

めの EDSL の開発は、有限体積法の分野には及んでこの問題の温度分布は、円柱の端点における境界条いない。しかし、有限体積法は数値流体力学のような件によって決定される。例えば、円柱の左端は温度一分野で用いられている重要な手法であり、有限要素法定の高熱源に接していて、右端は断熱されているとすにおけるのと同様に EDSL が有用であると期待できる。高熱源の温度を TB 、左端から右端までの長さをる。そこで、我々は、有限体積法のための C++テンL とすると、これらの境界条件は、プレートライブラリを開発して、EDSL として提供

T (0)dT (L)
dxすることを目指している。

=TB

(2)

=0

(3)

本稿では、我々の EDSL、有限体積法テンプレートの開発 †1 と、その自動並列化機能について報告すと表現される。

る。以下、第 2 節で有限体積法の概略を、第 3 節で方程式で記述されるが、計算機では連続関数に含まれ有限体積法テンプレートの設計・実装を、第 4 節でる無限の情報を扱うことができない。そこで、連続なそれを利用したプログラミングについて述べる。第 5

データーを離散化して、データ量を有限にする必要が節はまとめである。

ある。この離散化を行う手法のひとつが、有限体積法

このように、物理現象は連続関数についての偏微分である。

2 有限体積法プログラミングの概略
有限体積法 [5] [15] とは偏微分方程式を数値的に解くための手法である。偏微分方程式は様々な物理現象を記述できるので、例えば、拡散、対流といった流体現象を、有限体積法を使ってシミュレートすることができる。

2. 1 有限体積法による微分方程式の解き方

図1

連続領域の離散化

この問題でも、温度分布を連続関数ではなく離散有限体積法による物理現象のシミュレーションの特的なデータの列として扱わねばならない。そのため徴を示すために、空気中に置かれた円柱での熱伝導問に、円柱に対応する長さ L の連続領域を、図 1 のよ

題を考える。

うに、グリッドに置き換えることで、連続関数 T (x)問題を単純化するために、円柱の断面積方向の温度†1 我々は 有 限 体 積 法 テ ン プ レ ー ト を 実 装 す る た
め に 、オ ー プ ソ ー ス ソ フ ト ウェア の コ ー ド を
ベ ー ス と し て 利 用 し た 。そ こ で 、有 限 体 積 法 テン プ レ ー ト を 、オ ー プ ン ソ ー ス・ソ フ ト ウェアの 派 生 物 と し て 、動 作 が 確 認 で き た 試 作 版 を 、
https://github.com/masa-ito/ProtoToPoisson で公開している。

をグリッド点 xi の上の温度 Ti として、離散化する。
換言すれば、連続関数 T (x) を、離散化することで、i番目の要素が Ti であるようなベクトル t に変換する。
これに加えて、微分方程式のなかの演算子d2(4)k 2 − s Iˆddxも離散化することで、行列 C に変換できる。

こうして、微分方程式 (1) を離散化によって、Ct=b(5)書きなおさねばならないコードが多いということでもある。のような行列方程式（連立方程式）に変換できる。あそして、手動でのチューニングと並列化が、さらにとは、この方程式を境界条件 (2)、(3) のもとで補正再利用なコードを減らし、アルゴリズムとデータ構造して、数値的に解けば、温度分布 t が得られる。の変更を難しくしている。

2. 2 抽象化しないプログラミングの煩雑さ

2. 3 領域特化言語の効用

ただし、上記に比べて、実際に行なわれているプロ

このプログラミングの煩雑さが、有限体積法のためグラミングは遥かに煩雑であることが多い。この煩の領域特化言語を開発する動機である。ソースコード

雑さは、演算子、関数、グリッド、行列、ベクトルとの記述性とプログラムの実行速度を両立できさえす

いった概念を抽象データ型として記述せずに、その内れば、これまでの困難はかなり軽減されると期待で部詳細をむき出しにしたままでコーディングすることきる。

から生じている。

まず、アプリケーションプログラムの開発者は、有このようなプログラミングは面倒であるが、高速化限体積法の詳細に煩わされることなく、問題のモデリのために行なわれている。というのは、行列やベクトングに集中できる。また、モデリングへの専念は、領ルのような抽象概念をオブジェクトとして、数式をオ

域特化言語による並列化コードの生成（自動並列化）

ブジェクトと演算子の組み合わせとしてコーディングによっても、容易になる。

すると、その数式コードが評価されるときに、余計な一時オブジェクトが作られ、実行速度が低下するから

3 有限体積法テンプレートの設計と実装である [14]。

我々は、次のような要件を満たすように、有限体積とはいえ、ソースコードの記述性を下げるという犠法テンプレートの設計を行った。

牲には、非常に大きな苦痛が伴う。抽象化をしないと• 領域特化言語の中間コードを数式テンプレートいうことは、有限体積法の実装詳細をカプセル化せとすることで、オブジェクトを含む数式を評価すず、再利用もしないということだからである。

る時にも、実行速度が低下しないようにする。

例えば、微分方程式から行列方程式への置き換え• 既存の C++コードと領域特化言語のコードをの詳細 [5] [15] は、次のようになっている。微分方程式(1) を、グリッド点まわりの区間 ∆x †2 ごとに、混在させることができるように、領域特化言語を

d2s T∞ dx
k 2 T (x) − s T (x) dx = −dx∆xi∆xi†3のように積分して 、
k2kk
Ti−1 −
Ti +Ti+1 + s∆x Ti = −s T∞ ∆x∆x
∆x∆xのように、連立方程式に変換している。（ここで、• コンパイル時に、有限体積法の意味論に沿って、∫()i∫

独立型でなく埋め込み型とする。

オブジェクトの数式は基本的な C++コードに変換される。

• 行列方程式の数値解法を並列化スケルトンを組み合わせることで自動並列化する。

∆x = xi+1 − xi としている。）
抽象化をしないプログラミングでは、この変換を異3. 1 Boost.Proto ライブラリの利用なった偏微分方程式を扱うたびにやり直して、プログ

Boot.Proto は、EDSL を作成するための EDSL でラミングをしなおさねばならない。これはソフトウェある。この C++ライブラリを使うことで、数式テンアを開発するときに、再利用可能なコードが少なく、プレートを中間コードとする領域特化言語を作成す†2 ((xi−1 + xi )/2, (xi + xi+1 )/2) の区間、有限体積法
では検査体積と呼ばれる。
†3 高次元の問題ではガウスの発散定理を使う。

ることができる。
作成される領域特化言語は、次のような構成を持つ。

• フロントエンド – 開発者が Boost.Proto の枠組

算子は• バックエンド – 領域特化言語のコードから C++
コードへの変換規則

..


.

1/∆xd2 1/∆x→
dx2


−2/∆x

1/∆x

1/∆x−2/∆x

みで定義したオブジェクト

• 意味論 – バックエンドで定義された変換規則

1/∆x• （構文規則 – C++と共通）
まず、領域特化言語のフロントエンドとしてオブ1/∆x ..

.
(6)のように、行列へと変換される。

ジェクトを定義する。これで、そのオブジェクトを含

この数式変換を、領域特化言語ではコード 1 のよむ数式はコンパイル時に、自動的に数式テンプレート

うに、グリッドオブジェクト grid に微分演算子オブに変換されようになる。そして、バックエンドで、こジェクト DiffOpr の積を渡して、行列オブジェクトの数式テンプレートが最終的にどのような C++ コー

matrix へと変換させるコードとして表現する。

ドに変換されるを定義する。
この変換規則は、領域特化言語のコードがどのよう

matrix = grid . d i s c r e t i z e O p e r a t o r (diffOpr * diffOpr ) ;な C++コードに変換されるかを定義しているので、領域特化言語の意味論でもある。こうして、領域特化

コード 1：二階微分演算子を離散化して行列へ変換

言語のコードには、C++コードとは別の意味を持たこの微分演算子オブジェクトは領域特化言語のフせることができる。しかし、EDSL のコードなので、ロントエンドで定義されているので、オブジェクトのその構文規則は C++と同じである。

数式である diffOpr * diffOpr はコンパイル時に、数式テンプレートに変換される。更に、この数式テン

3. 2 領域特化言語の構成と数理的構造

プレートはバックエンドで定義された変換規則にした

Boost.Proto による領域特化言語の構成に、有限がって、二階の「差分」演算子オブジェクトへと変換体積法の数理的構造を対応させると、以下のようになる。

される。
ここで、二階の差分演算子クラスも領域特化言語• フロントエンド – 有限体積法を構成する概念をオブジェクトとして定義– グリッド、微分演算子、関数、行列、ベクトル• バックエンド – 有限体積法の数式の変換規則を定義のフロントエンドでコード 2 のように、Boost.Protoが提供するクラス proto::callable の派生クラスとして定義されている。
struct S e c o n d D i f f e r e n c e O p r: proto :: callable{// ...

– 微分演算子 → 行列オブジェクト– 定数項、関数 → ベクトルオブジェクト– 行列・ベクトル演算 → 基本的な C++コード• 意味論 – 有限体積法の数式が、どのようにして初等的な式に分解・還元されるか、を定義

3. 3 有限体積法の意味論でコードを変換する仕組み

S e c o n d D i f f e r e n c e O p r () {}SecondDifferenceOpr (const S e c o n d D i f f e r e n c e O p r & expr ) {}// 行列の対角要素を計算するため
double operator () ( double spacing ,double s p a c i n g O r B o u n d a r y S p a c i n g )const {
return ( - 1.0 / spacing - 1.0 /spacingOrBoundarySpacing );}2. 1 節で述べたように、有限体積法では二階微分演// 行列の非対角要素を計算するため

double operator () ( double spacing ,FaceToWestTag )const {
return 1.0 / spacing ;}double operator () ( double spacing ,FaceToEastTag )const {
return 1.0 / spacing ;}};コード 2：二階差分演算子クラスの定義
このクラスでは、行列の対角要素と非対角要素を計算するためのメンバ関数が別々に定義されているの行列とベクトルの積のように縮約演算であるときには、reduce スケルトンを割り当てる。
このように、数値解法アルゴリズムを構成する演算処理の性質に応じて、並列スケルトンを使い分けるのは、実行速度が改善するためである。

4 有限体積法テンプレートを用いたプログラ
ミング
4. 1 一次元熱伝導方程式ソルバ
2. 1 節の熱伝導問題の解き方を、有限体積法テンプレートを使い、領域特化言語でプログラミングすると、コード 4 のようになる。

で、二階差分演算子オブジェクトがグリッドオブジェクトのなかで評価されると、式（6）のように対角要

int main ( int argc , char * argv []) {// ... 定数の定義 ...

素と非対角要素が別々の値をもつ行列が生成される。

3. 4 行列・ベクトル演算の自動並列化
同様に、領域特化言語でコーディングされた行列やベクトルの数式も、コンパイル時に基本的な C++コードへと変換される。
例 え ば 、行 列 オ ブ ジェク ト A と ベ ク ト ル
オ ブ ジェク ト r,b,x で 構 成 さ れ る 数 式 コ ー
ド r = b - A * x; は、OpenMP で並列化されたコード 3 へと変換される。
# paragma omp parallel shared ( r )for ( int i = 0; i < sz ; i ++) {r . data [ i ] = 0.0;
for ( int j = 0; j < sz ; j ++)r . data [ i ] +=b . data [ i ] A . data [ i ][ j ] * x . data [ j ];}コード 3：並列化された行列・ベクトル演算
並列化の鍵となるのは、コード 3 の先頭にあるOpenMP の指示行である。この指示行は、並列スケルトンを組みわせる途上で追加される。その組み合わせは、領域特化言語のバックエンドの数式テンプレート変換規則で決定される。この変換規則は、行列・ベクトルの演算が、行列同士あるいはベクトルどうしの線型結合のように、要素ごとの演算であるときには、map スケルトンを割り当てる。これに対して、// グリッドオブジェクト
FVM :: Grid1D < FVM :: C e n t D i f f S c h e m e T a g >grid ( NumCtrlVol , Cyl inderL ength ) ;// 境界条件 (2)grid . a d d D i r i c h l e t B o u n d a r y (-1 , 0 , H otTem peratu re ) ;// 境界条件 (3)grid . a d d N e u m a n n B o u n d a r y (NumCtrlVol , NumCtrlVol - 1 , 0.0) ;// 恒等演算子
const FVM :: I d e n t i t y O p e r a t o r T y p e
IdOpr = FVM :: I d e n t i t y O p e r a t o r T y p e () ;// 微分演算子
const FVM :: D i f f e r e n t i a l O p e r a t o r T y p e
const
DiffOpr = FVM ::D i f f e r e n t i a l O p e r a t o r T y p e () ;// 微分方程式 (1) の左辺の演算子
auto opr = proto :: deep_copy (T h e r m a l C o n d u c t i v i t y * Area * DiffOpr* DiffOpr- ConvectiveHeatTransCoeff *Circumference * IdOpr ) ;// 微分方程式 (1) の離散化
Matrix coeffMat ( NumCtrlVol ,NumCtrlVol ) ;coeffMat = grid . d i s c r e t i z e O p e r a t o r (opr ) ;Vector rhsVec ( NumCtrlVol ) ;rhsVec = grid . d i s c r e t i z e F u n c t i o n (- ConvectiveHeatTransCoeff *Circumference *A mb i en t Te mp e ra t ur e ) ;// 境界条件による行列とベクトルの補正
FVM :: Bo und ary Cor re cto r bCorrector ( grid, opr ) ;bCorrector . applyTo ( coeffMat ) ;bCorrector . applyTo ( rhsVec ) ;// 数値解法の前処理として、対角スケーリングを使う
SLA :: D i a g o n a l P r e c o n d i t i o n e r precond (coeffMat ) ;# include< P a r a l l e l i z a t i o n T y p e T a g / OpenMP . hpp >と書いて、OpenMP のオプションをつけてコンパイルすることで、マルチスレッド化された実行ファイルが得られる。
並列パフォーマンスを測定したときには、コンパイ
ラは GNU C++ 4.8.5 を用い、最適化オプションとして// 前処理付きの共役勾配法で行列方程式を解く
SLA :: ConjugateGradient < DLA :: Matrix ,SLA :: D i a g o n a l P r e c o n d i t i o n e r >cg ( coeffMat , precond ) ;const Vector tempGuess ( NumCtrlVol ,(100.0 + 20.0) / 2.0) ;Vector temperature ( NumCtrlVol ) ;temperature = cg . solve ( rhsVec ,tempGuess , c o n v e r g e n c e C r i t e r i o n ) ;を指定した。

// ... 温度分布の数値解を書き出す ...

された点ヤコビ前処理付き共役勾配法の実行速度であreturn 0;- std = c ++11 - O3 - Winline \-- param max - inline - recursive - depth =32\-- param max - inline - insns - single =2000

4. 2. 1 自動並列化と手動並列化の比較
図 2 は、この実行ファイルの並列パフォーマンスを計測した結果である。図の縦軸は、マルチスレッド化る。この値を得るために、有限体積法の計算はグリッ
ド点 1500 個で行なわせて、共役勾配法の経過時間}の測定を繰り返して、80 回の測定値を平均した。横コード 4：一次元熱伝導方程式のソルバ

軸は、スレッド数である。測定は２プロセッサ x 12

このように、我々の領域特化言語を使えば、有限体コアのマシン（Intel Xeon E5-2697 v2( 30M Cache,積法の内部詳細に煩わされずにすみ、モデルの構成要2.70 GHz)、メモリ 128Gb）で行ったので、最大ス素である、グリッド、微分方程式、境界条件、数値解レッド数は 24 である。

法の指定に専念できる。つまりは、モデリングをするように、プログラミングをすることができる。

4. 1. 1 数値解の検証
コード 4 から得られた数値解は、この熱伝導問題の厳密解cosh(n(L − x))(7)cosh(nL) √と比較することで検証した。（ここで、n = s/k でT (x) = T∞ + (TB − T∞ )ある。）
プログラムは、L = 1 (m)、TB = 100 （度）、T∞ = 20 （度）、k = 1.0 × 103 （熱伝導係数） 、
n = 25(/m2 ) 、10 個のグリッド点で走らせた。その結果、誤差は最大でも 2%であることを確認できた。

4. 2 並列パフォーマンスの計測

図2

並列パフォーマンスの比較

我々の領域特化言語をプログラミングをするときには、ソースファイルの先頭に、

図 2 の青く丸いドットは、領域特化言語によって自動並列化されたプログラムのパフォーマンスで、赤く三角のドッドが、手動並列化されたプログラムのものである。手動で並列化プログラミングをしたときには、あえて、オブジェクを使わずに 2. 2 節のスタイルで、共役勾配法をコーディングした。

Kuchen が行った比較によれば、並列スケルトンで生成された行列乗算プログラムの実行速度は、MPIで並列化されたプログラムの 80%から 95%であった[4]。そして、これに劣らず、我々の自動並列化は手動並列化に匹敵するパフォーマンスを出している。
ただし、並列化のために用いた手法は単純であり、手動並列化と自動並列化の両方で、行列・ベクトル演算が要素ごとの計算か、縮約を含む計算かに応じて、OpenMP 指示行の入れ方を変えただけである。このため、自動並列化、手動並列化のいずれも、緑の破線で示された理想的なスケーリングに並列パフォーマン
スは及ばず、スレッド数が 16 を超えると実行速度が改善されなくなっているものと思われる。

5 まとめと今後の展望
以上のように、領域特化言語を有限体積法のプログラミングに使うことで、ソースコードの記述性とプログラムの実行速度の両方を改善することができた。
今後は、領域特化言語が提供する離散化モデル、疎行列、数値解法アルゴリズムの種類を増やしていく予定である。これにより、有限体積法の様々な離散化モデルに応じて粗行列のデータ構造を切り替え、また、データ構造に応じて行列方程式の数値解法を切り替えられるようにすることを目指す。また、モデルから行列方程式が生成され、その方程式を数値解法で解くという処理は、有限体積法以外のアプリケーション分野でも行なわれることが多い †4 。そこで、そのような分野でも、プログラムの記述性と実行速度の両立に役立つように、領域特化言語を拡張していきたい。

参 考 文 献[ 1 ] Bart Janssens, Támas Bányai, K. L. and Bosschaerts, W.: Finite Element Assembly Using anEmbedded Domain Specific Language, Scientific†4 例えば、機械学習

Programming, Vol. 2015, Article ID 797325(2015),
pp. 22.
[ 2 ] Cole, M.: Algorithmic Skeletons: Structured
Management of Parallel Computation, MIT Press,Cambridge, MA, USA, 1991.
[ 3 ] Estérie, P., Falcou, J., Gaunard, M., Lapresté,J.-T., and Lacassagne, L.: The Numerical Template Toolbox, J. Parallel Distrib. Comput., Vol. 74,No. 12(2014), pp. 3240–3253.
[ 4 ] H, K.: Proc. Annual European Conference on
Parallel Processing (Euro-Par2002), Lecture Notes
in Computer Science 2400.
[ 5 ] H.K.Versteeg and Malalasekera, W.: An Introduction of Computational Fluid Dynamics – The
Finite Volume Method 2nd Ed., Pearson Educational Limited, 2007.
[ 6 ] Klaus Iglberger, Georg Hager, J. T. and Rude,U.: Expression Templates Revisited: A Performance Analysis of Current Methodologies, SIAM J.
Sci. Comput., Vol. 34, No. 2(2012), pp. C42C69.
[ 7 ] Lani, A., Villedieu, N., Bensassi, K., Kapa,L., Vymazal, M., Yalim, M., and Panesi, M.:
COOLFluiD: an open computational platform for
multi-physics simulation and research, 21th AIAACFD Conference, San Diego (CA), 2013.
[ 8 ] M.Fowler and R.Parsons: Domain-Specific Languages, Addison-Wesley, 2011.
[ 9 ] Niebler, E.: Proto: A Compiler Construction
Toolkit for DSELs, Proceedings of the 2007 Symposium on Library-Centric Software Design, LCSD
’07, New York, NY, USA, ACM, 2007, pp. 42–51.
[10] Pietro, D. A. D., marc Gratien, J., and Prudhomme, C.: A domain-specific embedded language
in C++ for lowest-order discretizations of diﬀusive
problems on general meshes, BIT Numerical Mathematics, (2013), pp. 152.
[11] Prud’Homme, C., Chabannes, V., Doyeux, V.,
Ismail, M., Samake, A., Pena, G., Daversin, C., andTrophime, C.: Advances in Feel++: a domain specific embedded language in C++ for partial diﬀerential equations, Eccomas’ 12-European Congress
on Computational Methods in Applied Sciences andEngineering, 2012.
[12] Prud homme, C., Chabannes, V., Pena, G.,
and Veys, S.: Feel++: finite element embedded language in C++, Free Software available at
http://www. feelpp. org. Contributions from A.
Samake, V. Doyeux, M. Ismail and S. Veys.
[13] Vandevoorde, D. and Josuttis, N. M.: C++
Templates: The Complete Guide, Addison-Wesley
Professional, 1 edition, November 2002.
[14] Veldhuizen, T.: Expression Templates, C++
Report, Vol. 7, No. 3(June 1995), pp. 26 – 31.
[15] H.K. Versteeg(原著) W. Malalasekera(原著) 松下洋介ら (共訳): 数値流体力学第２版, 森北出版株式会社,2011.