3144

チュートリアルFortran（FORmula TRANslation）は半世紀以上の歴史を持ち、一部では時代遅れと言われながらも、今なお数値計算に利用する研究者が多いプログラミング言語です。Fortran90/95の機能や使用例を改めて理解したいという研究者のために、Fortran90/95による近年の有限要素法プログラムを題材にして、岐阜大学の永井学志先生と橋本一輝氏に解説をお願いいたしま
した。今回のチュートリアルでは、有限要素法における全体剛性行列の作成方法について解説していただきます。なお、チュートリアル記事は1ページ目のみを本誌に掲載し、続きは日本計算工学会 HP 上で公開していますので、そちらも併せてご参照ください。

有限要素計算における全体剛性行列の作成法－ 疎行列データ構造の視点から －永井 学志 橋本 一輝

1

はじめに

前回は、Fortran90/95と一部にFortran2003を用いて、抽象データ型（Abstract Data Type，ADT）プログラミングの考え方と実例を示しました。これは、オブジェクト指向プ
ログラミングから動的多態性と継承を除いたものです。
データを流していくプログラミングのスタイルに慣れていらっしゃる方には、取っ付きづらい概念かもしれません。しかし、内部のデータを隠蔽 ― カプセル化 ― し、
公開インターフェイスのみを介して外部とやり取りするというこのスタイルは、一度慣れてしまうと便利なものです。なぜならば、開発するプログラムには「ブツ（名詞）に、〜する（動詞）」という言語化がダイレクトに反映されるためです。したがって、対象とする問題を分割統治し損なわないかぎり、数年前の自分 ― 赤の他人 ―が書いたプログラムも読み返し易い（？）ように感じています。
今回は、このADTプログラミングの考え方に基づき、FEM計算で必須の疎行列、すなわち全体剛性行列Kと全体整合質量行列Mの組立てまでを述べます。行列の疎性を活かした可逆圧縮のデータ構造には色々とありますが、ここでは圧縮行記憶（Compressed Row Storage, CRS）
形式を対象とします。この圧縮データ構造は、非ゼロ成
分のみを記憶するもののなかで、基本的なものの1つで
しょう。
筆者紹介
ながい がくじ岐阜大学 機械工学科 准教授。1995年 東京理科大学
建築学科卒、2000年 東京工業大学 環境物理工学専

ところで、FEM創成期に諸先輩方が著された教科書には、少し不満があります。振り返ってみて、それは仕方がないと納得できることです。しかし仮に「疎行列やその組立て過程は、グラフ理論なる学問体系に他ならない
ので、そちらを学習しなさい」と強調しておいてくだ
さっていれば…、という思いです。小国力編著の本1）にも、グラフ理論はあまり知られていない旨の記述があります。全体自由度番号の付替えにはじまり、ウェーブフ
ロント法やマルチフロンタル法、さらには消去木とシン
ボリックLU分解2）など、いずれもこの理論の枠組みのようです。これらは我流の泥臭いアルゴリズムでも対応できるでしょうが、教科書的なアルゴリズムがあれば、そちらを学習するほうがきっと効率的です。
グラフ理論は図1に示すようなネットワーク ― 単にグ
ラフという ― を分析する道具であり、対応するアルゴ
リズムが色々と提案されています。この応用は多岐に渡っており、カーナビなどの経路探索や、コンパイラによるレジスタ割付け解析などに使われているようです。
グラフは、データ構造とアルゴリズムに関する教科書3）の後ろの方に載っています。配列とリストにはじまり、
ソート、2分探索、スタックとキュー、木、などを経てグラフに至ります。なお、筆者らがグラフというキー
ワードを知ったのは、電磁場解析で辺要素4）を使わざるを得なかったためです。
本稿では、第2章で疎行列の組立て過程をグラフ視点から俯瞰したのち、第3章で2分探索木のADTプログラミ
ングを、第4章でCRS形式までのADTプログラミングを述べます。なお、今回の逐次実行プログラムソースは、実用版を含め本学会誌HPに置いています。

攻 博士課程修了。

はしもと かずき
名古屋大学 複雑系科学専攻 修士課程在学中。2013年 岐阜大学 数理デザイン工学科卒。

図1

（26）

頂点（vertex）と枝（edge）で定義されるグラフ（ネットワーク、路線、管網、フレームなど）

計算工学

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)２全体剛性行列の組立ての概要と全体設計
まず，基本事項の確認からはじめさせてください．
骨組みのマトリクス構造解析や電気回路の節点解析な
どでは，対象とする系の変位や電位などの未知量を求めるために，「要素」と「節点」という概念を導入し，
図 1 のようなネットワーク図を描きます．そのうえで，系を要素にバラして ― よって節点の概念を要素節点と全体節点に分化させて ― 考えます．すなわち，
1) 節点内力の概念を導入し，各要素が自由体として満たすべき平衡方程式を立て，
2) こちらの都合で一度バラしてみたものの，本来は要素節点と全体節点で変位や電位などが一致している，という適合条件を課し，
3) 全体節点において，力や電流などがつりあっている，という平衡条件を課すことで，系全体が満たすべき平衡方程式を立て，未知変数について解くことになります．
線型問題であれば，系全体の係数行列 𝐊 は 𝑒 番目要素に関する係数行列を 𝐊 𝑒 として，数式上𝑛elm

𝐊 ← � 𝐋T𝑒 𝐊 𝑒 𝐋𝑒

(1)

𝑒=1

と記述できます．ここで，𝑛elm は全要素数， 𝐋𝑒 とその転 置 𝐋T𝑒 は 𝑒 番 目 要素 に 関す る 論理 型 の 長方 形行 列 であり，それぞれ上記手順 2), 3) に起因して出てくるものです．改めて，𝐊 は一般に大規模な疎行列，𝐊 𝑒 はせ
いぜい数十次元の密行列であることに注意ください．
非線型問題であれば，𝐊 はニュートン法のためのヤコビ行列であり，反復求解過程と増分過程で繰返し作成
されます．
ここから具体例として，図 1 を何らかの物理問題として，かつ節点がスカラ自由度である場合，すなわち
節点と自由度を同一視できる場合について，話を進めます．まず単純化のため，行列 𝐊 𝑒 と𝐊 について，それらの成分がゼロか否かについてのみ注目してみます．
成分がゼロか否かを 0 と 1 に対応づけて，それぞれ論
理型行列 𝐀𝑒 と𝐀 とします．𝐊 𝑒 は一般に密行列なので，𝐀𝑒 = �

1 1
�
1 1

(2)

とすべて 1 にしておくのが最も安全です．𝐀 について
は，式(1)の総和 Σ を論理和 ∪ で置き換えて，𝑛elm𝐀 ← � 𝐋T𝑒 𝐀𝑒 𝐋𝑒

(3)

𝑒=1

と表記できます．この例での結果は，
11b(2) ⎡ 1c(3) ⎢⎢d(4)𝐀=⎢
e(5) ⎢
f(6) ⎢ 1
g(7) ⎣ 1
a(1)

計算工学

2

1
1
1
1

3

1
1

4

5

1
1

1
1 1
1 1
1

1

6

1

7

1
1 ⎤
⎥
⎥
1 ⎥
1 ⎥
⎥
1 ⎦

となります．ここで，行列 𝐀 の空欄成分は 0 であり，
行と列それぞれの並び順は全体節点の名称の辞書順― a, b, c, …, g ― です．
式(3)の組立ての実際は後回しにして，その結果である式(4)を解釈してみます．この式は，全体節点間のつながり関係を真偽，1 と 0 として，リーグ表にしただ
けです．たとえば，図 1 のグラフでは，全体節点 a(1)は要素を介して，自分自身を含め全体節点 a(1), b(2),f(6), g(7) につながっています．この意味で，式(2)は𝑒 番
目要素の要素節点同士がすべて相互につながっている
ことを，式(3)はこの関係の全体節点への変換と統合を表しています．
ここで今後の説明を容易にするため，グラフ理論 3)における用語をいくつか紹介させてください．論理型
行列 𝐀 を「隣接行列（adjacency matrix）」と称します．
また，通常は自由度 ― ここでは節点 ― を「頂点（vertex）」，2 頂点間のつながり ― ここでは要素 ―を「枝（branch）」と称します．さらに，全成分が 1 である 𝐀𝑒 に対応するグラフを「完全グラフ」，その頂点群を「クリーク（clique：小集団）」と称します．なお，
この例の場合，図 1 のグラフから直接的に隣接行列を構成するほうが容易です．しかし，一般の有限要素計
算への応用においては，一度メッシュ絵を離れて隣接行列から新たにグラフを思い描くほうが好都合です．
本題に入ります．本稿の主題は式(1)の疎行列 𝐊 の組
𝑛elm
立て過程 ∑𝑒=1
です．主題の多くは，いま式(3)の隣接
𝑛elm
行列 𝐀 の組立て過程 ⋃𝑒=1に要約されました．また，非線型計算では，式(3)の組立ては基本的に一度で済みますが，式(1)の組立ては繰返し必要です．そこで，これら 2 つの組立て過程をそのまま親子 2 段構えで継承
し，親である隣接行列（ADJacency MaTriX）に関するモジュール adj_mtx_class と，子である疎行列（SParse
MaTriX）に関するモジュール sp_mtx_class に切り分けます．
両モジュールは，前回で述べた ADT プログラミン
グの考え方で記述します．すなわち，データ構造とア
ルゴリズムの抽象化により，抽象データ型とそのイン
ターフェイスとして公開する総称名サブルーチン ―抽象データ「ブツ」に寄り掛かる作用素「～する」 ―と，カプセル化すべき実装本体とに分離できます．そこで，実装本体は次章以降で述べるとして，ここでは使用者の視点から使い勝手を考え ― 実際には作成者の視点から計算効率も考え ―，両モジュール内の抽象
データ型 adj_mtx_c と sp_mtx_c と，それぞれに対する作用素の用法を先に定義します．具体的に，次のprogram1 のような用法にて定義します．Program1 では
すべてのサブルーチン CALL が作用素に該当し，それぞれ第 1 引数が作用を受ける抽象データです．

(4)（Program1：逐次実行版の主プログラム）
1:PROGRAM main_serial_ver
2: USE constants! 倍精度実数型の定数DPの定義
3: USE adj_mtx_class ! adj_mtx_c型に関するADT定義
4: USE sp_mtx_class! sp_mtx_c 型に関するADT定義(26-2)

Vol.19, No.3 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)5::
6: IMPLICIT none
7::
8: TYPE(adj_mtx_c):: A_pp!隣接行列
9: TYPE( sp_mtx_c):: K_pp, M_pp
!疎行列
10:
:
11: INTEGER, PARAMETER :: nelm = 9!図1の要素数12: INTEGER, PARAMETER :: nDOF_p= 7!正の全体DOF数
13: INTEGER, PARAMETER :: nDOF_n= 0!非正の 〃
14: INTEGER, PARAMETER :: nDOFe = 2!要素DOF数
15: INTEGER:: DOFe2g(nDOFe) !DOFの変換TBL
16: REAL(DP):: Ke(nDOFe,nDOFe),& !要素剛性17:& Me(nDOFe,nDOFe)
!要素質量
18:
:
19: !【第1段：親】隣接行列A_ppの作成
20: CALL init( A_pp, nDOF_p )𝑛elm
21: DO k = 1, nelm
! ⋃𝑒=1
22:
DOFe2g(:) = …
23:
CALL add_clique( A_pp, DOFe2g )24: END DO
25: ! 確認用途：完成したA_ppのファイル出力
26: CALL wrt( A_pp, 'A_pp.dat' )27::
28: !【第2段：子】疎行列K_ppとM_ppの作成
29: CALL init( K_pp, A_pp )30: CALL init( M_pp, A_pp )𝑛elm
31: DO k = 1, nelm
! ∑𝑒=1
32:
DOFe2g(:) = …
33:
Ke(:,:)= …
34:
Me(:,:)= …
35:
CALL add_clique( K_pp, Ke, DOFe2g )36:
CALL add_clique( M_pp, Me, DOFe2g )37::
38: END DO
39: ! 確認用途：完成したK_ppのファイル出力40: CALL wrt( K_pp, 'K_pp.dat' )41::
42:END PROGRAM main_serial_ver

こ こ で ， プ ロ グ ラ ム 中 の 自 由 度 （ DOF: Degree Of
Freedoms）に関して，補足させてください．全体自由度番号は正負に渡る連番として，ディレクレ条件に関して非正の番号（対応して接尾語_n）を，それ以外に正の番号（接尾語_p）を付けます．計算ではディレク
レ条件以外の自由度についてのみ疎行列が必須なので，数式上𝐀=

デ ィ レク レ ≤0 非ディ レ クレ 0<

ディ レク レ ≤0

非ディ レク レ 0<

�

𝐀nn
𝐀Tnp

𝐀np
𝐀pp

�

(5)

と分割して，𝐀pp に関してのみ考えます．また，配列に よ る 変 換 テ ー ブ ル DOFe2g(:) ― 要 素 自 由 度 番 号（Element DOF num.）から全体自由度番号（Global DOFnum.）を引く ― の命名法に関しても補足させてくだ
さい．変換は函数𝒚 = 𝒇(𝒙) の言い換えに過ぎないので，𝒇 に対応する配列名を x2y（ from x to y）とすることで，可読性を上げたつもりです．
話を戻して，Program1 を見ていきましょう．まず，3，4 行目でモジュール adj_mtx_class と sp_mtx_classの使用を宣言します．8，9 行目でそれぞれの抽象デー
タ型 adj_mtx_c と sp_mtx_c で，隣接行列 A_pp と疎行

計算工学

列 K_pp，M_pp を宣言します．そのうえで，20 行目で隣接行列 A_pp を初期化（CALL init）し，21～24 行目
で小集団（完全グラフ）を加算（CALL add_clique）します．26 行目は必須でないですが，完成した隣接行列
A_pp のファイル出力（CALL wrt）です．疎行列 K_pp，
M_pp についても同様に，29, 30 行目で初期化（CALL
init）し，31～38 行目で要素剛性行列 Ke と整合質量
行列 Me の小集団を加算（CALL add_clique）します．
より詳細に見ると，29, 30 行目の CALL init では，
隣接行列 A_pp を元に疎行列 K_pp，M_pp を初期化しています．元々，疎行列(1)から隣接行列(3)を切り出して考えたので，元に戻す処理が必要です ― オブジェク
ト指向でいうところの「継承」と「ダウンキャスト」 ―．
ここで，抽象データ型 adj_mtx_c と sp_mtx_c の内訳を先読みしておきます．疎行列 K_pp，M_pp は同じ隣接行
列 A_pp からの派生であるため，あえてデータのカプセ
ル化を POINTER で破り，両疎行列から A_pp の同じデー
タ実体を指すようにします ― 詳細は次章 ―．このことのメリットは，メモリ効率の良さ以上に，複数の疎行列が対象としている問題の同一性が言えることです．
たとえば，時刻歴解析のニューマークβ法では，スカラ値を c として K_pp+c*M_pp を計算する作用素 add が必要ですが，両疎行列は同一の隣接行列 A_pp のデータ実体を指していますので，対象問題の同一性は既に保障されています．以上が使用者向けの話です．
３ 結局，抽象データ型としての隣接行列の組立て
３.１ 隣接行列と疎行列のデータ関係の切り分け
ADT プログラミングに向けて，隣接行列と疎行列のデータ関係の切り分けと，疎行列 sp_mtx_c 型の内部データ構造の設計からはじめます ― ここからは使用者
でなく作成者の視点 ―．出発点は，式(4)の隣接行列 𝐀に等価な圧縮表現の 1 つ，
{ 1, 2, 6, 7 },⎧{1, 2, 3, 7 },⎪c(3) ⎪ { 2, 3 },d(4) { 4, 5, 7 },e(5) ⎨ { 4, 5, 6, 7 },f(6) ⎪⎪ { 1, 5, 6 },g(7) ⎩ { 1, 2, 4, 5, 7 }
a(1)

b(2)

(6)

です．この表現は，隣接行列 𝐀 の成分が 1 の列位置を行毎にまとめて集合 {… } とし，さらにそれらを縦に並べたものです．グラフ理論 3) では各集合 {… } を「隣接リスト（adjacency list）」と称します．この名称は，図1 のグラフと式(6)を見比べると，実に自然です．全体
節点に関して，a(1)とつながっているものをすべてリ
ストアップすると，自分自身を含めて a(1), b(2), f(6),g(7)です．
隣接行列 𝐀pp を元に初期化される疎行列 𝐊 pp などでは，計算効率の視点から圧縮データ構造として CRS 形式が好まれます．CRS 形式は，表現(6)を区切りのための補助配列とともに 1 次元配列に展開・拡張したものといえます．具体的に sp_mtx_c 型の内訳となる CRS(26-3)

Vol.19, No.3 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)形式は，次の program2 のようになります．
親：隣接行列 adj_mtx_c 型（adj_mtx_class モジュール内）（Program2：sp_mtx_c と adj_mtx_c の関係）
1:MODULE sp_mtx_class
2: USE constants! 倍精度実数型の定数DPの定義
3: USE adj_mtx_class ! adj_mtx_c型に関するADT定義
4::
5: IMPLICIT none
6: PRIVATE!基本は非公開
7::
8: TYPE, PUBLIC :: sp_mtx_c!抽象データ型
9:
PRIVATE
10:
INTEGER, POINTER :: ind(:) => NULL() !区切り11:
INTEGER, POINTER :: DOF(:) => NULL() !列位置
12:
REAL(DP), ALLOCATABLE :: a(:)!疎行列の成分13::
14: END TYPE sp_mtx_c
15::
16: PUBLIC init
!総称名の公開
17: INTERFACE init
!作用素として
18:
MODULE PROCEDURE init_mtx_CRS!(源氏名の正体)
19: END INTERFACE
20::
21:CONTAINS
22::
23: SUBROUTINE init_mtx_CRS ( this, adj )24:
TYPE( sp_mtx_c), INTENT(inout) :: this25:
TYPE(adj_mtx_c), INTENT(inout) :: adj !隣接行列
26:
:
27:
CALL associate( adj, this%ind, this%DOF )28:!この自作ダウンキャスト作用素はadj_mtx_class29::
30: END SUBROUTINE int_mtx_CRS
31::
32:END MODULE sp_mtx_class

前章の最後で先読みしたように，16～30 行目の初期化
init 中，27 行目 CALL associate により adj_mtx_c 型
データ実体を指すようにしています．たとえば表現(6)の 場 合 ， 列 位 置 の 配 列 this%DOF(:) に は (/1,2,6,7,
1,2,3,7, 2,3, 4,5,7, 4,5,6,7, 1,5,6, 1,2,4,5,7/)，
区切りのための補助配列 this%ind(:)には(/1, 5, 9,
11, 14, 18, 21, 26/)が入っています．
しかし，上述のように組立て終えたのちに CRS 形式とするのは問題ないですが，隣接行列 𝐀 pp の組立て過
𝑛elm
程 ⋃𝑒=1
を CRS 形式で行うのは無理があり，これを抽象データ型 adj_mtx_c の内訳とはできません．この問
題は組立て過程𝑒 = fe(1)𝑒 = fe(1) ~ ed(4)∪ { },∪ { 6, 1 },⎧⎧b(2)∪ { },∪ { },⎪⎪
c(3) ⎪ ∪ { },⎪ ∪ { },
⇒⋯⇒
⇒⋯
d(4)
∪ { },∪ { 5, 4 },⎨⎨ ∪ { 6, 5; 5, 7; 5, 4 },e(5)∪ { 6, 5 },⎪ {f(6) ⎪⎪ ∪ { 6, 5 },⎪ ∪ 6, 5; 6, 1 },g(7) ⎩ ∪ { }⎩ ∪ { 5, 7 }
a(1)

(7)

を考えると自明です．ここで，要素番号 𝑒 の番号並びは，要素名を頂点名 a, b, c,…を 2 つ書き連ねて，図 1の上から右向きに fe(1), fa(2), eg(3), ed(4), …, bc(9) としています ― 辞書順にすると次節の説明が最悪ケー
計算工学!組立てモード用! ()!計算モード用(CRS 形式)
INTEGER, POINTER:: ind(:)
INTEGER, POINTER:: DOF(:)

CALL associate(…)

REAL(DP), ALLOCATABLE :: a(:)!!
:
子：疎行列 sp_mtx_c 型（sp_mtx_class モジュール内）

図2

隣接行列 adj_mtx_c 型と疎行列 sp_mtx_c 型の設計

スに（笑） ―．式(7)で注意頂きたいことは，組立て
𝑛elm
過程⋃𝑒=1を完了しないかぎり各隣接リスト {… } の成
分数は未知なので，CRS 形式の 1 次元配列 this%DOF(:)に前もって展開できないという点です．各隣接リスト{… } は集合なので，並び順と重複に意味がないことに
𝑛elm
も注意く ださい ．な お，組 立て過程 ⋃𝑒=1
（program1の 19～24 行目）の前に全要素をループで走査すれば，各 隣 接 リ ス ト {… } の 成 分 数 の 上 限 だ け は お さ え ら れて，配列を割付けることができますが，泥臭いデータ
構造とアルゴリズムになってしまいます．
そこで，ADT プログラミングの考え方では使用者が抽象データ型の内訳を知り得ないことを逆手にとります．すなわち，隣接行列 adj_mtx_c 型の内訳として，
図 2 に模式的に示すように，データ構造を 2 つ，組立てモード用の何らかのデータ構造と，計算モード用のCRS 形式とを用意します．そのうえで，隣接行列 A_ppを元に初めて疎行列をダウンキャスト初期化するタイ
ミ ン グ ― program1 の 29 行 目 で CALL init(K_pp,
A_pp)，これから program2 の 16～30 行目が呼出され，さらに同 27 行目で CALL associate ― にて，内部データ構造を組立てモード用から計算モード用に，不可
逆で移行させます．
３.２ 抽象データ型としての隣接行列の組立て
ここでは，モジュール adj_mtx_class における抽象
データ型 adj_mtx_c と，それに対する作用素 init，
add_clique，associate について，実装本体を考えます．データ構造は，前節で 2 つのモードに切り分けたので，組立てモードに特化しましょう．式(7)から分かるように，隣接行列を可逆圧縮した状態のまま組立てることから，工夫が必要です．式(6),(7)では縦方向と横方向があるため，2 層構造として考えるのが自然で𝑛elmに要する時間は後の非線
す．また，組立て過程 ⋃𝑒=1
型計算を考えると無視できるので，行毎に伸びていく
メモリ割付けの空間効率を重視しましょう．
𝑛elm
組立て過程 ⋃𝑒=1の例(7)などをメモリ効率よく実現するには，配列表現でなくリスト表現の類が適してい
ま す ． リ ス ト は 構 造 型 宣 言 TYPE と ポ イ ン タ 宣 言
POINTER を組み合わせて，再帰的に数珠つなぎしていくものです．ここで，配列表現に比べてリスト表現による初期化がどれほど遅くなるのか，連結リスト（linked list） ― 線形リスト ― の場合で実験してみ(26-4)

Vol.19, No.3 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)ら一巡して根に戻ってくることを「横断（traversal）」と称します．このとき，図中の各成分周りに付した領
域 a, b, c のうち，領域 b を通り掛け（in-order）する
際にその成分を順次ピックアップしていくと，(/ 1, 2,
4, 5, 7 /) と昇順整列済みのデータ列を得ることができます．したがって，組立てモード用のデータ構造を
2 分探索木としておき，初めて associate が呼出された際に，CRS 形式に移行させるのは容易です．これで図 2 中の空欄（ ）が「2 分探索木の森」で埋まり，模式図が完成しました．
そこで，隣接行列モジュール adj_mtx_class の本体実装では，2 分探索木（Binary Search Tree）モジュー
ル BS_tree_class の抽象データ型 tree_c を下請けとして切り分けて，次の program3 のようにします．

root
a

5

c

b
a

1

a

c

c

b

b
NULL

7

4

a

c

NULL

NULL

b
a

2

c

NULL

b
NULL

図3

NULL

頂点 g(7)の隣接リスト組立てに対する 2 分探索木

ま し た ． 約 50 倍 遅 く な り ま し た ． 原 則 的 に 逐 一
ALLOCATE 命令でシステム側 にメモリ割当てを要求す
るので仕方ないことです．しかし，4byte データを 25M
回（約 100Mbyte）繰返し割当てるために要する時間はわずか約 2 秒でしたので，無視します．
いま，隣接リストの集合 ∪ {… } 毎に考えます．組立て過程(7)と完成版(6)を比べると，常に昇順を保ちつつ，
clique（小集団）の加算毎に重複分については無視し，新規分のみをしかるべき位置に挿入していくのが得策でしょう．この作業に最適なデータ構造は，「2 分探索木（binary search tree）」3) です．これも TYPE と POINTERの組み合わせで実現されるものです．なお，Java 言語
では TreeSet という木のクラスがあらかじめ備わっています．
2 分探索木を完成させた一例を，図 3 に示します．
この例は，式(7)の最終行 g(7)を完成させた隣接リスト∪ { 5, 7; 1, 7; 7, 4; 7, 2 } に対するものです．この木の作り方は，図中に示すように最上部を根（root）として，再帰的に 2 分しつつ下に木を茂らせていきます．このとき，右下に直上よりも大きな値を，左下に直上よりも小さな値を入れます．この例で最初から描くと，1)
値 5：根が指す直下のメモリ域に挿入，2) 値 7：根から辿り値 5 の右側ポインタが指すメモリ域に挿入，3)
値 1：根から辿り値 5 の左側ポインタが指すメモリ域に挿入，4) 値 7：根から辿り値 5 の右側を辿ると値 7は挿入済みなので探索終了，…となります．描画上は2 分ですが，処理上は，小さい場合，等しい場合，大きい場合の 3 分であることに注意ください．
式(7)全体では，このような木の「森（forest）」とな
ります．2 分探索木のメリットは，整列済みの各リス
ト の 成分 数 を 𝑚 ― 疎 行 列 では 数 十程 度 ― と す れば，探索に平均 𝑂(log 2 𝑚)，その後の挿入には 𝑂(1) の軽い処理回数となることです．なお，配列では工夫す
れば探索を𝑂(log 2 𝑚) とできるものの，その後のしかるべき位置への挿入には再割付けの手間を除いても𝑂(𝑚)，逆に連結リストではしかるべき位置への挿入には 𝑂(1) なものの，探索に𝑂(𝑚) を要します．
図 3 のような 2 分探索木があるとき，図中に示すよ
うに，根の左側をスタートして木を常に左手に見なが計算工学（Program3：adj_mtx_class の実装）
1:MODULE adj_mtx_class
2: USE BS_tree_class ! tree_c型に関するADT定義
3::
4: IMPLICIT none
5: PRIVATE!基本は非公開
6::
7: TYPE, PUBLIC :: adj_mtx_c!抽象データ型
8:
PRIVATE
9:! 組立てモード用（2分探索木の森）
10:
TYPE(tree_c), ALLOCATABLE :: tree(:)
11:! 計算モード用（CRS形式）
12:
INTEGER, POINTER :: ind(:) => NULL() !区切り13:
INTEGER, POINTER :: DOF(:) => NULL() !列位置
14: END TYPE adj_mtx_c
15::
16: PUBLIC init
!総称名の公開
17: INTERFACE init
!作用素として
18:
MODULE PROCEDURE init_adj_mtx
19: END INTERFACE
20:
21: PUBLIC add_clique!総称名の公開
22: INTERFACE add_clique!作用素として
23:
MODULE PROCEDURE add_clique_adj_mtx
24: END INTERFACE
25:
26: PUBLIC associate!総称名の公開
27: INTERFACE associate!作用素として
28:
MODULE PROCEDURE associate_adj_mtx
29: END INTERFACE30::
31:CONTAINS
32::
33: SUBROUTINE init_adj_mtx ( this, nDOF_p )34:
TYPE(adj_mtx_c), INTENT(inout) :: this35:
INTEGER,INTENT(in) :: nDOP_p
36:
ALLOCATE( this%tree(nDOF_p) )!木の森を準備
37: END SUBROUTINE init_adj_mtx
38:
39: SUBROUTINE add_clique_adj_mtx ( this, DOFe2g )40:
TYPE(adj_mtx_c), INTENT(inout) :: this41:
INTEGER,INTENT(in) :: DOFe2g(:)
42:
INTEGER :: buf(SIZE(DOFe2g))!自動割付け配列
43::44:
cnt = 045:
DO ie = 1, SIZE(DOFe2g)
46:
i = DOFe2g(ie)(26-5)

Vol.19, No.3 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)47:
IF ( i > 0 ) THEN48:
cnt = cnt +1
49:
buf(cnt) = i!正の全体DOF番号のみを抽出
50:
END IF
51:
END DO
52:
DO ie = 1, cnt!木々に複数データを一括挿入
53:
i = buf(ie)54:CALL insert( this%tree(i), buf(1:cnt) )55:
END DO
56: END SUBROUTINE add_clique_adj_mtx
57:
58: SUBROUTINE associate_adj_mtx ( this, ind, DOF )59:
TYPE(adj_mtx_c), INTENT(inout) :: this60:
INTEGER, POINTER, INTENT(out ) :: ind(:)
61:
INTEGER, POINTER, INTENT(out ) :: DOF(:)
62:
:
63:
IF ( ALLOCATED(this%tree) ) THEN !初CALL時のみ
64:
nnz = 0! モード移行65:
max_len = 0
66:
nDOF_p = SIZE(this%tree)
67:
DO i = 1, nDOF_p!CRS作成の準備
68:
length = get_ndat( this%tree(i) )69:nnz = nnz + length!Num. of Non-Zeros70:
IF ( max_len < length ) max_len = length71:
END DO
72:
ALLOCATE( this%ind(nDOF_p+1) )73:
ALLOCATE( this%DOF(nnz))
74:
ALLOCATE( adj_list(max_len ) ) !取出用バッファ
75:
p = 1
76:
this%ind(1) = p
77:
DO i = 1, nDOF_p!CRS作成78:
CALL get_dat( this%tree(i), adj_list, length)79:!対称行列はプリプロセッサ命令 #ifdef で
80:
DO m = 1, length
81:
this%DOF(p) = adj_list(m) !列位置82:
p = p +1
83:
END DO84:
this%ind(i+1) = p!区切り
85:
CALL final( this%tree(i) )!木の解放(必須)86:
END DO87:
DEALLOCATE( adj_list )88:
DEALLOCATE( this%tree )89:
END IF90:
ind => this%ind(:)!区切りデータの共有
91:
DOF => this%DOF(:)!列位置 〃〃
92: END SUBROUTINE associate_adj_mtx
93::
94:END MODULE adj_mtx_class

ここで，抽象データ型 adj_mtx_c の内訳は 7～14 行目
であり，それに対する作用素 init の実装は 16～19，
33～37 行目，add_clique は 21～24，39～56 行目，
associate は 26～29，58～92 行目です．なお，疎行列の対称性を考慮する場合には，associate 内にプリプ
ロセッサ命令#ifdef ～ #else ～ #endif を，SJIS の5C 文字問題に注意しつつ書くのも一案です．
また，2 分探索木の実装は ADT の考え方によりほぼ隠されたことに注意ください．2 行目でそのモジュー
ル BS_tree_class の使用を宣言したうえで，10 行目で
adj_mtx_c の 1 成分として，その抽象データ型 tree_cの動的割付け配列である森 tree(:)を宣言しています．
それに対する作用素の用法はそれぞれ，54 行目：木の探索および挿入の insert，68 行目：木の成分数を返す
計算工学

get_ndat，78 行目：木の全成分をバッファ配列に取出
す get_dat，85 行目：木を解放する final です．
蛇足ですが，コンパイラとシステムに負荷を掛けがち な 処 理 と し て ， 42 行 目 の 自 動 割 付 け 配 列
buf(SIZE(DOFe2g))と，第 54 行目の部分配列の実引数
buf(1:cnt)があります．前者について，自動割付け配
列はスタックメモリに積まれることが多いですが，小さな配列なので特に問題ありません．後者について，
部分配列を実引数とするとき配列記述子 ― 上下限やストライドなど ― が対応する仮引数に渡されますが，コンパイラがストライド 1 に限定した最適化を掛ける
とき ― Intel fortran では手続き間最適化の –ipo オプ
ション ―，速度低下の問題は起こりません．
３.３ 最下請けとしての 2 分探索木
プログラムを完結させるため，2 分探索木モジュー
ル BS_tree_class についても，実装本体の一例を次のprogram4 に示しておきます．抽象データ型 tree_c に対する作用素，insert と get_dat については，内部サ
ブルーチンによる再帰版を示します．スタックメモリのオーバフローを恐れる場合には，煩雑ですが非再帰で書き換えられます．
（Program4：BS_tree_class の実装）
1:MODULE BS_tree_class
2: IMPLICIT none
3: PRIVATE!基本は非公開
4: TYPE, PRIVATE :: vtx!2分探索木の基本要素
5:
INTEGER :: c
6:
TYPE(vtx), POINTER :: L => NULL()
7:
TYPE(vtx), POINTER :: R => NULL()
8: END TYPE vtx
9:
10: TYPE, PUBLIC :: tree_c!抽象データ型
11:
PRIVATE12:
TYPE(vtx), POINTER :: root => NULL()
13:
INTEGER :: nvtx = 0!頂点数(for効率化)
14: END TYPE tree_c
15:
16: PUBLIC :: insert!総称名の公開
17: INTERFACE insert!作用素
18:
MODULE PROCEDURE insert_data_R
19: END INTERFACE
20:
21: PUBLIC :: get_ndat!作用素
22:
23: PUBLICT get_dat!総称名の公開
24: INTERFACE get_dat!作用素
25:
MODULE PROCEDURE traversal_asend_R
26: END INTERFACE
27::
28:CONTAINS
29::
30: SUBROUTINE insert_data_R ( this, a )31:
TYPE(tree_c), INTENT(inout) :: this32:
INTEGER,INTENT(in) :: a(:)!挿入データ
33:
:
34:
DO i = 1, SIZE(a)35:
CALL insert_kernel( this%root ) !内部Sub.
36:
END DO(26-6)

Vol.19, No.3 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)37: CONTAINS!内部Sub.（上位変数は有効）
38:RECURSIVE SUBROUTINE inset_kernel ( ptr )39:
TYPE(vtx), POINTER, INTENT(inout) :: ptr
40:
IF ( ASSOCIATED(ptr) ) THEN
!終端でない41:IF( a(i) == ptr%c ) THEN
42:
RETURN!探索終了
43:
ELSE IF ( a(i) < ptr%c ) THEN
44:
CALL insert_kernel( ptr%L ) !左下へ探索
45:
ELSE
46:
CALL insert_kernel( ptr%R ) !右下へ探索
47:
END IF
48:
ELSE!終端に到達49:
ALLOCATE( ptr )!メモリ割付50:
ptr%c = a(i)!挿入
51:
this%nvtx = this%nvtx +1
52:
RETURN53:
END IF
54:
END SUBROUTINE insert_kernel
55: END SUBROUTINE insert_data_R
56:
57: FUNCTION get_ndat ( this ) RESULT( ans )58:
TYPE(tree_c), INTENT(in) :: this59:
INTEGER:: ans60:
ans = this%nvtx
61: END FUNCTION get_ndat
62:
63: SUBROUTINE traversal_asend_R ( this, buf, n )64:
TYPE(tree_c), INTENT(in) :: this65:
INTEGER,INTENT(out ) :: buf(:)
66:
INTEGER,INTENT(out ) :: n
67:
:
68:
cnt = 069:
CALL trvsl_inorder( this%root ) !内部sub
70: CONTAINS!内部Sub.（上位変数は有効）
71:
RECURSIVE SUBROUINE trvsl_inorder( ptr )72:
TYPE(vtx), POINTER, INTENT(in) :: ptr
73:
IF ( .NOT. ASSOCIATED(ptr) ) THEN !終端に到達
74:
RETURN75:
ELSE!終端でない
76:
CALL trvsl_inorder( ptr%L )! a)左下へ
77:
cnt = cnt +1! b)中央78:
buf(cnt) = ptr%c!
出力79:CALL trvsl_inorder( ptr%R )! c)右下へ80:
END IF81:
END SUBROUTINE trvsl_inorder82: END SUBROUTINE traversal_asend_R
83::
84:END MODULE BS_tree_class４

抽象データ型としての疎行列の組立て
疎行列 sp_mtx_c 型に対する作用素 add_clique の実
装を次の program5 のようにして，program2 に追加し
ま す ． CRS 形 式 の 区 切 り this%ind(:) と 列 位 置
this%DOF(:)を元に疎行列を組立てるには，this%a(:)
において疎行列 𝐊 pp の (𝑖, 𝑗) 成分が格納される位置 p を探索する必要があります．この探索には，整列済み配列データに対する「2 分探索（binary search）」 3) を用います．
2 分探索の一例として図 4 に，g(7)の隣接リストから
値 5 の格納位置 4 を探索する様子を示します．再帰的に 2 等分を繰返して追い込んでいくことから，この探索は𝑂(log 2 𝑚) であり，いま𝑚 は数十程度なので実質的に 10 回未満です ― それでも重い処理ですが… ―．
計算工学

g(7)の隣接リスト（昇順整列済み）
1 回目:

2 回目:

(1)

(2)

(3)

(4)

(5)

1

2

4

5

7

L

R
m=(1+5)/2
L
m=(4+5)/2

3 回目:(探索終了条件 L>R)
図4

R

R

L

昇順整列済みデータに対する 2 分探索の例

22～26 行目にあるように，疎行列の第 i 行目のすべての 列 位 置 デ ー タ は 昇 順 に て this%DOF(this%ind(i):this%ind(i+1)-1)に 格 納 され て お り ， こ こ から CALL
binary_search して第 j 列目のオフセット格納位置 loc
を探索します．よって，𝐊 pp の (𝑖, 𝑗) 成分の格納位置は
p = this%ind(i) -1 +loc と計算できます．
（Program5：sp_mtx_class への実装追加分）
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:

:
USE sort_and_search_module, ONLY: binary_search:
PUBLIC add_clique!総称名の公開INTERFACE add_clique
MODULE PROCEDURE add_clique_CRS_mtx
END INTERFACE:SUBROUTINE add_clique_CRS_mtx ( this, Ke, DOFe2g )TYPE(adj_mtx_c), INTENT(inout) :: thisREAL(DP),INTENT(in) :: Ke (:,:)
INTEGER,INTENT(in) :: DOFe2g(:):
nDOFe = SIZE(DOFe2g)
DO ie = 1, nDOFe
i = DOFe2g(ie)IF ( i > 0 ) THEN
DO je = 1, nDOFe
j = DOFe2g(je)!対称行列はプリプロセッサ命令 #ifdef で
IF ( j > 0 ) THEN
loc = binary_search( j,&& this%DOF(this%ind(i):&&this%ind(i+1)-1)) )p = this%ind(i) -1 + loc
this%a(p) = this%a(p) + Ke(ie,je)
END IF
END DOEND IF
END DOEND SUBROUTINE add_clique_CRS_mtx:

2 分探索の関数 binary_search の本体は，program6に示すように，ADT プログラミングの考え方でなく従来型として，とりあえず関数・サブルーチンの単なる寄せ集めモジュール sort_and_search_module に入れておきます．等号の存在ゆえに，名称にある 2 分でなく実質的に 3 分となるためか，短いながらもややこしいプログラムの代表例 3) です．

(26-7)

Vol.19, No.3 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (2)（Program6：整列済み配列データに対する 2 分探索）
1:MODULE sort_and_search_module
2: IMPLICIT none
3: PRIVATE
4: PUBLIC binary_search
5::
6:CONTAINS
7::
8: FUNCTION binary_search ( j, a ) RESULT( loc )9:
INTEGER, INTENT(in) :: j!探索すべき値
10:
INTEGER, INTENT(in) :: a(:) !昇順整列済み
11:
:
12:
left = 1
13:
right = SIZE(a)14:
DO WHILE ( left <= right )15:
middle = ( left + right ) / 2
16:
IF ( j <= a(middle) ) right = middle -1
17:
IF ( a(middle) <= j ) left = middle +1
18:
END DO19:
it_exists = ( left - right == 2 ) !存在の有無
20:
IF ( it_exists ) THEN21:
loc = middle!格納された位置
22:
ELSE
23:
loc = left
!挿入すべき位置24:
END IF25::
26: END FUNCTION binary_search
27::
28:END MODULE sort_and_search_module

計算工学

５

おわりに
思っていた以上に，隣接行列の話が長くなりました．
今回は，前回示した ADT プログラミングの考え方を復習しつつ，グラフや 2 分探索木を新たに導入しました．前回に続き，掲載が遅くなり申し訳ございませんでした．次回の最終回では，要素のマルチカラー化 ―ふたたびグラフの話 ― による作用素 add_clique 周
りの Open-MP 並列化や，疎行列 sp_mtx_c 型に対する
作用素 solve 内への PARDISO ソルバの組込み，さらに，いくつかの落穂 ― 連立 1 次方程式の反復求解と相性がよい，Reverse Cuthill-Mckee（RCM）法による全体自由度番号の並び替えなど ― を拾う予定です．
参考文献
1) 小国力編著, 行列計算ソフトウェア : WS、スーパ
ーコン、並列計算機, 丸善, 1991
2) たとえば，G.W. Stewart, Building an Old-Fashioned
Sparse Solver, 2003, http://hdl.handle.net/1903/1312(2014 年 6 月現在)
3) た と え ば ， 石 畑 清 , ア ル ゴ リ ズ ム と デ ー タ 構 造 ,岩波書店, 1989
4) 五十嵐一ら, 新しい計算電磁気学, 培風館, 2003

(26-8)

Vol.19, No.3 2014