3189

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

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

1

はじめに

前回は、抽象データ型（Abstract Data Type, ADT）プロ
グラミングの考え方により、隣接行列 adj_mtx_c 型と疎行列 sp_mtx_c 型を定義し、それぞれに対する作用素としての総称名サブルーチン init（初期化）、add_
clique（組立て）などを、情報処理の視点から定義しました。今回の最終回は、1）要素のマルチカラー化に
よ る add_clique 周 り の Open-MP 並 列 化 と、2）連 立1次方程式の求解ライブラリ PARDISO を用いて、疎行列
sp_mtx_c 型に対する作用素 solve の実装を考えます。
また、3）隣接行列 adj_mtx_c 型に対して、全体自由度
番号を付替える作用素 renumber にも言及します。
先に、Open-MP 並列化を施して最終形とした主プロ
グラムを、次の program1に示しておきます。
（Program1：Open-MP 並列版の主プログラム）
1:PROGRAM main_openMP_ver
2: USE constants! 倍精度実数型の定数 DP の定義
3: USE adj_mtx_class ! adj_mtx_c 型に関する ADT 定義
4: USE sp_mtx_class! sp_mtx_c 型に関する ADT 定義
5: USE mesh_module! マルチカラー化など
6::
7: IMPLICIT none
8:
:

筆者紹介
ながい がくじ岐阜大学 機械工学科 准教授。1995年 東京理科大学
建築学科卒、2000年 東京工業大学 環境物理工学専
攻 博士課程修了。

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

Vol.19, No.4 2014（31）

9: TYPE(adj_mtx_c):: A_pp! 隣接行列10: TYPE( sp_mtx_c):: K_pp, M_pp
! 疎行列
11:
:
12: INTEGER, ALLOCATABLE :: elm2DOFg(:,:)! 接続行列
13: TYPE(jag),ALLOCATABLE :: col2elm (:) ! 色毎の要素
14: ! 第 ic 色の全要素は col2elm(ic)%c(:) ― ギザギザ配列
15::
16: ! 前回の具体例を再び用いて…
17: INTEGER, PARAMETER :: nelm = 9! 全要素数
18: INTEGER, PARAMETER :: nDOF_p= 7! 最大全体 DOF 番号
19: INTEGER, PARAMETER :: nDOF_n= 0! 最小〃
20: INTEGER, PARAMETER :: nDOFe = 2! 要素 DOF 数
21: INTEGER:: DOFe2g(nDOFe) !DOF の変換 TBL
22: REAL(DP):: Ke(nDOFe,nDOFe)! 要素剛性
23::
24: REAL(DP), ALLOCATABLE :: U(:) ! 変位 (nDOF_n:nDOF_p)
25: REAL(DP), ALLOCATABLE :: F(:) ! 外力 (1:nDOF_p)26::
27: !【第 0 段】要素のマルチカラー化 ― col2elm(ic)%c(:)
28: ALLOCATE( elm2DOFg(nDOFe,nelm) )29: DO k = 1, nelm
30:
elm2DOFg(:,k) = …
31: END DO
32: CALL mk_col2elm_from_elm2DOFg( col2elm, elm2DOFg )33: DEALLOCATE( elm2DOFg )34:
35: !【第 1 段：親】隣接行列 A_pp の作成 ―36: !$OMP PARALLEL
37: DO ic = 1, SIZE(col2elm)! 第 ic 番目色
38:!$OMP DO PRIVATE( p, k, DOFe2g )39:
DO p = 1, SIZE(col2elm(ic)%c)! 並列実行40:
k = col2elm(ic)%c(p)! 第 k 番目要素41:
DOFe2g(:) = …
42:
CALL add_clique( A_pp, DOFe2g )43:
END DO44:!$OMP END DO
45: END DO
46: !$OMP END PARALLEL
47::
48: CALL renumber( A_pp, node2DOFg )! 番号付替え
49::
50: !【第 2 段：子】疎行列 K_pp の作成 ―51: CALL init( K_pp, A_pp )! ダウンキャスト
52: !$OMP PARALLEL
53: DO ic = 1, SIZE(col2elm)! 第 ic 番目色
54:!$OMP DO PRIVATE( p, k, DOFe2g, Ke )55:
DO p = 1, SIZE(col2elm(ic)%c)! 並列実行56:
k = col2elm(ic)%c(p)! 第 k 番目要素
57:
DOFe2g(:) = …
58:
Ke(:,:)= …
59:
CALL add_clique( K_pp, Ke, DOFe2g )
60:
:
61:
END DO62:!$OMP END DO
63: END DO
64: !$OMP END PARALLEL
65::
66: CALL solve( K_pp, F, U(1: ) )
! 求解
67:
:
68:END PROGRAM main_openMP_verチュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)２

要素のマルチカラー化と組立て工程の並列化
隣接行列 とそれに対応する全体剛性行列 の組
立て工程 ⋃, ∑について，それぞれの Open-MP 並列化に先立ち，要素のマルチカラー化処理が必要です． すな わち ，最 終形 であ る Program1
中の第 35～46 行目および第 50～64 行目に先立ち，第27～33 行目が必要です．そこで，まず Open-MP 並列化に関して概説しておき，つぎに要素のマルチカラー化に入っていきたいと思います．
２.１ 組立て工程の Open-MP 並列化にあたって
一般に Open-MP 並列化は，Fortran コードへ !$OMP
ではじまる指示子 ― ディレクティブ ― を挿入することで実現します．様々な指示子がありますが，多くの場合は Program1 の DO ループに対するもののように，!$OMP PARALLEL ～ !$OMP END PARALLEL ブロックによる複数スレッドの起動指示と，その内部の !$OMP
DO PRIVATE(…) ～ !$OMP END DO ブロックによる DO
ループの各スレッドへの割当て指示で間に合います．
Open-MP 並列化はメモリ共有型環境に対するものゆえ，
デフォルトでは，変数は全スレッドで共有されます．
しかし，並列化の対象となる DO ループの制御変数と，
PRIVATE(…)内で列挙された変数だけは，スレッド毎に独立となります．本プログラムの場合，スレッド毎に独立な変数は，要素に関する制御変数 p，いま注目している要素番号 k，要素自由度番号から全体自由度番号への変換テーブル DOFe2g(:) ― from Element DOF
num. TO Global DOF num. ― ，要素剛性行列 Ke(:,:)
です．
Open-MP 並列化では，メモリ書込みに依存関係のな
いこと，すなわち複数スレッドが同じメモリ番地に同時に書込む可能性のないことが前提です．Open-MP 並列化された DO ループ内で関数・サブルーチンを呼出す
場合でも，その内部において同じです．メモリ書込みの依存関係を整理するのはコンパイラでなく開発者自身であり，データ構造とアルゴリズムに注意しておく必要があります．いま考えている組立て工程, ∑には総和計算があることか⋃
ら，メモリ書込みには依存関係があります．ところが，本プログラムで CALL している作用素としての総称名
サブルーチン add_clique は，これらの組立て工程の一部に過ぎません．したがって，メモリ書込みの依存関係を整理するには，作用素 add_clique 内で大域的変数への書込みがないことを前提として，以降で説明する
マルチカラー化の視点が必要です．
２.２ 要素のマルチカラー化の考え方
例として図 1(a) に示すような，4 節点四辺形要素による有限要素分割 ⋃を考えます．要素のマルチカラー化とは彩色問題 1) の一種であり，図 1(b) に示す
ように，節点 ― より厳密には自由度 ― を共有していない要素群を同色で塗り分けるものです．たとえば，
図中，赤で着色した要素群の各要素は，互いに自由度

計算工学(a) 4 節点四辺形要素による有限要素分割(b) 要素の彩色結果（計 6 色）(c) 要素を頂点，節点（自由度）を枝とみなしたグラフ
図1

要素のマルチカラー化の例

を共有していません．その他の色の要素群についても同じです．したがって組立て工程において，同色の要
素群に対する CALL add_clique(…)には，変数 A_pp やK_pp の同じメモリ番地に書込む可能性がないため，並列実行できます．色毎のこの並列実行を全色について順に繰り返せば，組立てが完了します．
数式上は，
を色数，
1, 2, ⋯ ,を色の識別番号， を 番目の色に属する要素集合として，

←

,
∈

←

(1)
∈

と記述できます．このようにメモリ書込みの依存関係を整理することで，依存関係がない内側の総和計算はOpen-MP 並列化できて，依存関係のある外側の総和計
算を逐次実行します．
番目の 色に属 する要 素集合 は， Program1 では
jagged 配列まがい ― 初回で説明 ― として実装し，ic 番 目 の 色 か ら 要 素 集 合 へ の 変 換 テ ー ブ ル
col2elm(ic)%c(:)としています．この変換テーブルは，組立て工程での使用に先立って作成します．すなわち，(31-2)

Vol.19, No.4 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)Program1 の第 27～33 行目がこれに該当し，要素から全体自由度番号への変換テーブル elm2DOFg(:,:)を一時的に作成したうえで，従来型サブルーチン
mk_col2elm_from_elm2DOFg を CALL します．なお，一時的であっても比較的大きな配列 elm2DOFg(:,:)を用意することは野暮です．しかし，連立 1 次方程式求解のために最大のメモリを割付ける前であるため，この野暮さは我慢することにします．
いよいよ要素の彩色法についてです．一般に彩色はグラフ ― 前回に説明：「頂点（vertex）」とそれを結ぶ「枝（branch）」で定義されるネットワーク ― に関する問題であり，グラフの頂点や枝をできるだけ少ない色数で塗り分けようとするものです．より好まし
い彩色結果を得るためには問題に応じたものが必要でしょう．しかし，本稿では簡便かつもっとも有名なア
ルゴリズムとして，グラフの頂点彩色に対する
Welsh-Powell のもの 1), 2) を紹介します．なお，本アル
ゴリズムは，次節にてグラフの隣接行列モジュール
adj_mtx_class 中に実装し，adj_mtx_c 型への作用素としての総称名サブルーチン do_coloring とします．
Welsh-Powell の彩色アルゴリズムはグラフの頂点に対するものなので，まずは図 1(c) に示すように，要素分割図を，要素を頂点，全体自由度を枝とするグラフに描き替えます．なお，前回のグラフは今回と違って，全体自由度を頂点，要素を枝としたものであり，互いに双対的となっていることに注意ください．今回のグラフに描き替えたうえで頂点に彩色するのが，先ほど言 及 し た Program1 中 の 従 来 型 サ ブ ル ー チ ン
mk_col2elm_from_elm2DOFg で す ． そ の 骨 子 を 次 のProgram2 に示します．本質的に要素と全体自由度を入れ替えただけのグラフなので，Program1 中の要素に関する組立て工程とほぼ同じです．
（Program2：mk_col2elm_from_elm2DOFg の実装）
1:MODULE mesh_module
2: IMPLICIT none
3: PRIVATE
4: PUBLIC :: mk_col2elm_from_elm2DOFg
5::
6:CONTAINS
7::
8: SUBROUTINE mk_col2elm_from_elm2DOFg &9:& ( col2elm, elm2DOFg )10:
USE adj_mtx_class!隣接行列(前回に説明)
11:
12:
TYPE(jag),ALLOCATABLE,INTENT(out)::col2elm(:)
13:
INTEGER,INTENT(in )::elm2DOFg(:,:)
14:
15:
TYPE(jag),ALLOCATABLE :: DOFg2elm(:)
16:
INTEGER, ALLOCATABLE :: elm2col (:)
17:
TYPE(adj_mtx_c) :: dual_A !要素が頂点のグラフ
18:
:
19:
! elm2DOFg(:,:)の逆引きDOFg2elm(:)%c(:)を作成
20:
CALL mk_Y2X_from_X2Y_2D( DOFg2elm, elm2DOFg )
21:
22:
nelm
= SIZE (elm2DOFg, dim=2)
23:
nDOF_n = LBOUND(DOFg2elm, dim=1)!下限 <=024:
nDOF_p = UBOUND(DOFg2elm, dim=1)!上限 >=0

計算工学

25:
26:
! 逆引きDOFg2elm(i)%c(:)に関する隣接行列を作成
27:
CALL init( dual_A, nelm )28:
DO i = nDOF_n, nDOF_p!全全体自由度を走査
29:
CALL add_clique( dual_A, DOFg2elm(i)%c )30:
END DO
31:
32:! グラフ頂点としての要素を彩色
33:
ALLOCATE( elm2col(nelm) )34:
CALL do_coloring( dual_A, elm2col )!Program3参照
35:
CALL final( dual_A )36:37:! elm2col(:)から逆引きcol2elm(ic)%c(:)を作成
38:
CALL mk_Y2X_from_X2Y_1D( col2elm, elm2col )39:
40: END SUBROUTINE mk_col2elm_from_elm2DOFg
41::
42:MODULE mesh_module

この Program2 では，まず第 19～20 行目にて，要素から全体自由度番号への変換テーブル elm2DOFg(:,:)
をもとに，その逆引きテーブル DOFg2elm(i)%(:)を作っています．つぎに第 26～30 行目にて，全体自由度番号について逐次 DO ループさせることで，要素をグラフ頂点とする隣接行列 dual_A を作成しています．その後に 第 32 ～ 35 行 目 で 頂 点 彩 色 し て ， そ の 結 果 を
elm2col(:)に返して，最後に第 37～38 行目でその逆引
き col2elm(ic)%c(:)を作成して，呼出し元に返しています．なお，順方向の整数変換テーブルから逆方向の整数変換テーブルを作成する従来型サブルーチン
mk_Y2X_from_X2Y_* 内は，一般に 2 回走査方式が素直
かと思います．すなわち，第 1 走査目で逆変換配列のための上下限を探索し，その動的メモリ割付をしてから，第 2 走査目で実際に逆変換テーブルを作成します― 詳説は略 ― ．
２.３ Welsh-Powell の頂点彩色アルゴリズム 1), 2)
最後に Welsh-Powell の頂点彩色アルゴリズム自体を記 す と と も に ， 隣 接 行 列 adj_mtx_c 型 へ の 作 用 素
do_coloring の実装までを記しておきます．
まずアルゴリズムは以下です．
1) グ ラ フ 各 頂 点 を ， 頂 点 に 接 続 し て い る 枝 の 数 ―「次数（degree）」と称し，隣接リスト長さから 1だけ引いたもの ― の降順に並び替える．すなわち，
頂点の 次 数 を degと し て ， deg⋯
deg2) この降順にて，各頂点に彩色可能な最小の色識別番
号を割り当てる．
このとき，必要色数の上限は次式で与えられます．
max min

, deg

1

deg

1

(2)

図 2 の例では，deg
10 ，
6 でした．
つ ぎ に そ の 実 装 で あ る do_coloring 本 体 を 次 のProgram3 に示します．

(31-3)

Vol.19, No.4 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)（Program3：adj_mtx_class へ do_coloring を追加）
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
23:
24:
22:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
30:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:

:
USE sort_and_search_module, ONLY: sort:
PUBLIC :: do_coloring
INTERFACE do_coloring
MODULE PROCEDURE do_coloring_by_Welsh_Powell
END INTERFACE:SUBROUTINE do_coloring_by_Welsh_Powell && ( this, vtx2col )TYPE(adj_mtx_c),INTENT(in ):: thisINTEGER,INTENT(out):: vtx2col(:) !(nvtx)
INTEGER:: nvtx!全頂点数
INTEGER,ALLOCATABLE:: vtx (:) !(nvtx)頂点番号
INTEGER,ALLOCATABLE:: deg1(:) !(nvtx)頂点次数+1
INTEGER:: deg1_i, max_deg1
INTEGER,ALLOCATABLE:: buf (:) !取出用バッファ
LOGICAL,ALLOCATABLE:: used(:) !(0:max_deg1)色用
INTEGER:: ncol, icol !色数, 色番号! 1) 頂点vtx(:)を(次数+1)のdeg1(:)順で昇順ソート
nvtx = SIZE(this%tree)
ALLOCATE( vtx(nvtx), deg1(nvtx) )DO i = 1, nvtx
vtx (i) = i
deg1(i) = get_ndat( this%tree(i) )END DOCALL sort( vtx, deg1 )max_deg1 = deg1(nvtx) !最大次数! 2) 頂点彩色は次数の降順かつ最小値の色でALLOCATE( buf ( max_deg1),&& used(0:max_deg1))!未着色を0にused(:) = .FALSE.
!全色未使用に
vtx2col(:) = 0!未着色=未訪問
ncol = 1
DO k = nvtx, 1, -1!次数の降順で
i = vtx(k)!注目頂点は i
CALL get_dat( this%tree(i), buf, deg1_i )! 頂点iの隣接頂点で使用されている色を調査
used(1:ncol) = .FALSE.
!全色未使用に
DO q = 1, deg1_i
j = buf(q)!隣接行列の(i,j)成分
col = vtx2col(j)!col=0の未着色(未訪問)はused(col) = .TRUE. ! used(0)に入れて無処理
END DO! 未使用色のうち，最小値の色を探索
DO col = 1, ncol
IF ( .NOT. used(col) ) EXIT
END DO
!完遂後はcol==ncol+1
vtx2col(i) = col!頂点の色が決定
IF ( ncol < col ) ncol = col !used(1:ncol)用END DOEND SUBROUTINE do_coloring_by_Welsh_Powell:３

連立１次方程式の求解
連立 1 次方程式を求解するために，疎行列モジュール sp_mtx_class 中に，sp_mtx_c 型に対する作計算工学

用素 solve を実装します．前回すでに sp_mtx_c 型の内部 デ ー タ 構 造 で あ る CRS 形 式 と ， 対 応 す る 作 用 素
add_clique までの実装を終えています．したがって，
もはや solve 内部で CRS 形式のデータを求解ライブラ
リに引渡すだけです．本チュートリアルの動機は，汎用の求解ライブラリを使って，並列かつメモリ効率よく直接求解したいというものでした．そこで，次に示
す Program4 にて，sp_mtx_c 型に対する作用素 solveの実装として，Intel MKL 版の PARDISO ライブラリを組込みます．なお，PARDISO 引数の詳細は MKL マニ
ュアル 3) に譲ります．
（Program4：sp_mtx_class に solve を追加）
1:INCLUDE 'mkl_pardiso.f90' !Intel MKLに同梱
2:
3:MODULE sp_mtx_class4::
5: USE mkl_pardiso!mkl_pardiso.f90内で定義
6::
7: TYPE, PUBLIC :: sp_mtx_c
8:
PRIVATE
9:
INTEGER, POINTER:: ind(:) => NULL() !区切り10:
INTEGER, POINTER:: DOF(:) => NULL() !列位置11:
REAL(DP),ALLOCATABLE:: a (:)!疎行列の成分12:!*** 上記変数は前回の再掲，下記変数は新規 ***13:
INTEGER:: state!状態遷移フラグ14:
TYPE(mkl_pardiso_handle)::pt(64)!MKLオマジナイ
15: END TYPE sp_mtx_c
16::
17: PUBLIC :: solve
18: INTERFACE solve
19:
MODULE PROCEDURE solve_mtx_CRS
20: END INTERFACE
21::
22:CONTAINS
23::
24: SUBROUTINE solve_mtx_CRS ( this, b, x )25:
TYPE(sp_mtx_c),INTENT(inout) :: this26:
REAL(DP),INTENT(inout) :: b(:) !(n)右辺
27:
REAL(DP),INTENT(inout) :: x(:) !(n)解28:!b(:)とx(:)のinout指定はmkl_pardiso由来
29:30:!以下はPARDISOに必要な定数・変数群の宣言
31:
INTEGER :: maxfct= 1!同一疎構造の行列の数
32:
INTEGER :: mnum =maxfct!行列の識別番号
33:
INTEGER :: mtype = 2!正定値行列（#ifdefで34:! = 11非対称行列 切り分け）
35:
INTEGER :: phase!ソルバの実行を制御
36:
INTEGER :: n!次元数
37:
INTEGER :: idummy(1)!ダミー配列38:
INTEGER :: nhrs = 1
!右辺の数
39:
INTEGER :: iparm(64)= 0 !各種パラメタ格納用
40:!( 1)= 0 で(2:64)はデフォルト値
41:
INTEGER :: msglvl= 0!1:情報出力,0:なし
42:
INTEGER :: error
43:
:
44:
SELECT CASE ( this%state ) !状態遷移に応じて，
45:
CASE ( not_yet_solved) ! その疎構造は初
46:this%pt(:)%dummy = 0! 初回だけ必須
47:
phase = 13
! シンボリック分解(1)→代入(3)48:CASE ( new_mtx_to_solve ) ! その係数は初49:phase = 23
! 数値分解(2)→代入(3)50:CASE ( solved_the_mtx) ! その係数は分解済(31-4)

Vol.19, No.4 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)51:phase = 33
! 代入(3)のみ
52:
END SELECT
53:
this%state = solved_the_mtx
54:
55:
CALL pardiso &56:& ( this%pt ,maxfct ,mnum
57:& phase,n,this%a58:& this%DOF ,idummy ,nrhs
59:& msglvl
,b
,x
60:
:
61: END SUBROUTINE solve_mtx_CRS
62::
63:END MODULE sp_mtx_class, mtype,&, this%ind ,&, iparm,&, error)この Program4 では，疎行列 sp_mtx_c 型の成分として，作用素 solve のために内部状態変数 state を組込
み，適宜遷移させています．内部状態は 3 つあり，「その疎構造の疎行列に対する求解は初」，「その係数の疎行列に対する求解は初」，「その係数の疎行列は分解済／前処理因子を作成済」です．PARDISO であれば，第44～53 行目の 1) 対象となる疎行列構造のシンボリッ
ク分解，2) 数値分解，3) 前進・後退代入から成る phase= 13, 23, 33 にそれぞれ対応します．
なお，汎用ライブラリによる直接求解でなく自作の反復求解を選ぶ場合には，内部状態変数 state に対応
させて 3 種類のルーチン： 1) 反復計算のための作業
配列の確保，2) 前処理因子の作成，3) 反復求解のルーチンに分けます．そのうえで，solve 内で効率的なFortran77 + Open-MP ルーチンへのバックドアを用意します．
RCM 法による全体自由度番号の付替え
汎用の直接求解ライブラリ PARDISO を使う分には，内部的に METIS 4) が使われているので，全体自由度番
号の付替え（re-ordering）は不要です．しかし，自作の反復求解では，メモリの空間的局在性を上げるため
４

f (2)

e (3)

括弧内は頂点の次数

g (4)

d (2)

a (3)

c (1)

b (3)

待ち行列を１次元配列で実現
1

0)

2

3

4

5

6

7

c
head

tail

b

2)

c

b

a

g

3)

c

b

a

g

f

a

g

f

図2

計算工学

…

7)

…

c

…

1)

c

b

d

e

に Reverse Cuthill-McKee（RCM）法程度があると好都合です．そこで，隣接行列 adj_mtx_c 型に対して全体
自由度番号を付替える作用素 renumber も実装しておきます．なお，この renumber は，adj_mtx_c 型の内部
データ構造が「計算モード」でなく「組立てモード」 ―前回で説明 ― にあることが前提です．
前回のグラフ（図 2 に再掲）を例にとり，RCM 法の元となる Cuthill-McKee（CM）法から述べます．CM
法は，隣接する頂点同士の番号ができるだけ近くなる
よう，地道にコツコツと各頂点に番号を付けていくものです．情報処理の視点 1), 5) からすると，グラフの全頂点を「幅優先探索（Breadth First Search, BFS）」により訪問する問題の一般化と考えることができます．幅
優先探索の本質は，「待ち行列（queue）」と「先入れ先出し（First-In-First-Out, FIFO）」です．
図 2 の図解は，CM 法により全頂点を訪問する手順です．図中のグラフで，頂点名 a, b, …, g 直後の括弧内は次数 ― 隣接頂点数 ― です．また，待ち行列を，
全頂点数分の長さの 1 次元配列により実現します．待
ち行列の先頭位置を head，末尾の直後位置を tail とします．このとき，CM 法による訪問手順を文章化す
ると，次の通りです．
0) 次数が最小の頂点を探す．頂点 c が該当するので，訪問．初訪問の頂点 c を待ち行列の末尾に追記．
1) 待ち行列の先頭にある頂点 c に注目し，これに隣接する全頂点を訪問．初訪問の頂点 b を，次数の昇順にて待ち行列の末尾に追記．
2) 待ち行列の先頭にある頂点 b に注目し，これに隣接する全頂点を訪問．初訪問の頂点 a, g を，次数の昇順にて待ち行列の末尾に追記．
3) 待ち行列の先頭にある頂点 a に注目し，これに隣接する全頂点を訪問．初訪問の頂点 f を，……．
::
7) 待ち行列の先頭にある頂点 e に注目し．これに隣接する全頂点を訪問．初訪問の頂点は皆無なので，先頭のみを１つ後ろへ繰る．head が tail に追い付いたので，訪問完了．
なお，待ち行列で並んでいる頂点群をグラフ上にプロットすると，前線を形作っており，少しずつ地道に陣地拡大していると解釈できます．
訪問完了後，１次元配列上に残されたデータ順 c, b,…, d, e が，CM 法による頂点番号の付替え結果です．
頂点名を辞書順として番号付けすると，この待ち行列を実現する１次元配列は，頂点について新番号から旧番号を引く逆引きとなっています．そこで，この配列
名を n2o(:) ― from New TO Old/Original ― としておきます．なお，この逆引きは１対１の対応にあるので，その順引き o2n(:)は in を反復の制御変数として，o2n(n2o(in))=in で作成できます ―．
RCM 法による頂点番号の付替えは，CM 法で得られたデータ順を逆順にして，e, d, …, b, c とするだけです．
付替え前後の隣接行列を示すと，

Cuthill-McKee 法による頂点番号の付替え例(31-5)

Vol.19, No.4 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)

a
a
b
c
d

b

c

f

d
f
g
a
b
c

d

f

g

1 1
1
1 1
1
1 1 1 1
1 1
1 1
1

1
1 1
e

e

e

1 1
1 1 1
1 1

e
g

d

f

g

a

b

(3)

c

1 1 1 1
1 1
1
1
1
1
1 1
1 1 1
1 1 1 11 1 1 11 1

(4)であり，付替え前に比べて付替え後は，対角近くに成分が集中しています．実装上は，この番号付替えを，2分探索木による隣接リスト表現に対して実施します．
次の Program5 に，隣接行列 adj_mtx_c 型に対する RCM
作用素 renumber の実装を示します．続く Program6 には，2 分探索木 tree_c 型に対する番号付替え作用素
renumber の実装を示します．
（Program5：adj_mtx_class に renumber を追加）
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:
33:
34:
35:
36:

:
PUBLIC :: renumber
INTERFACE :: renumber
MODULE PROCEDURE :: remunber_adj_mtx
END INTERFACE:SUBROUTINE renumber_adj_mtx ( this, node2DOFg )TYPE(adj_mtx_c),INTENT(inout) :: thisINTEGER,INTENT(inout) :: node2DOFg(:,:)
INTEGER,
ALLOCATABLE :: n2o(:), o2n(:)
TYPE(tree_c),ALLOCATABLE :: tree_tmp(:):!全体DOF番号を付替えるために順引きo2n(:)を作成
nDOF_p = SIZE(this%tree)!
ノイマンDOF数nDOF_n = MINVAL( node2DOFg ) !-ディレクレDOF数
ALLOCATE( o2n(nDOF_n:nDOF_p),&!順引き& n2o(
1:nDOF_p)) !逆引き
CALL mk_RCM_table( this%tree, n2o ) !39行目参照
FORALL (in=nDOF_n:0) o2n(in ) = in
FORALL (in=1:nDOF_p) o2n(n2o(in)) = in
DEALLOCATE( n2o ):! o2n(:)をthis%tree(:)とnode2DOFg(:,:)に反映
ALLOCATE( tree_tmp(nDOF_p) )!$OMP PARALLEL DO PRIVATE( in )DO io = 1, nDOF_p
in = o2n(io)tree_tmp(in) = this%tree(io)!行の移動
CALL renumber( tree_tmp(in), o2n(1:nDOF_p) )END DO! ↑Program6を参照!列の入替えCALL MOVE_ALLOC( tree_tmp, this%tree ) ! F2003
nDOFn = SIZE(node2DOFg,1)!節点あたりDOF数
nnode = SIZE(node2DOFg,2)!全節点数
FORALL (i=1:nDOFn, j=1:nnode) && node2DOFg(i,j) = o2n( node2DOFg(i,j) )
:

計算工学

37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:

END SUBROUTINE renumber_adj_mtx
SUBROUTINE mk_RCM_table ( tree, n2o )TYPE(tree_c),INTENT(in) :: tree(:)
INTEGER,INTENT(out ) :: n2o (:)
INTEGER, ALLOCATABLE :: buf (:)!バッファ
INTEGER, ALLOCATABLE :: deg1(:)!次数 +1LOGICAL, ALLOCATABLE :: visited(:) !訪問済フラグ:!CM法により，幅優先探索型でグラフの全頂点を訪問
max_deg1 = …
! 最大次数 +1
ALLOCATE( buf (max_deg1),&& deg1(max_deg1))
nvtx = SIZE(tree)
ALLOCATE( visited(nvtx) )visited(:) = .FALSE.
j= …
!最小次数の頂点から訪問
visited(j) = .TRUE.
tail = 1n2o(tail) = j!キュー末尾に追記
tail = tail +1
DO head = 1, nvtx
!幅優先型探索の本体
i = n2o(head)!キュー先頭の頂点iに注目
CALL get_dat( tree(i), buf, deg1_i )cnt = 0
!初訪問の頂点のみをカウント
DO q = 1, deg1_i
j = buf(q)IF ( .NOT. visited(j) ) THEN !頂点jは初訪問
visited(j) = .TRUE.
cnt = cnt +1
buf (cnt) = j! 初訪問の頂点を抜出
deg1(cnt) = get_ndat( tree(j) ) !頂点次数END IF
END DO!優先順位づけ：初訪問の隣接頂点buf(1:cnt)を，!
deg1(1:cnt)をキーとして昇順ソート
CALL sort( buf(1:cnt), deg1(1:cnt) )!キュー末尾に追記
n2o(tail:tail+cnt-1) = buf(1:cnt)
tail = tail +cnt
END DODEALLOCATE (visited, deg1, buf )!RCM法はn2o(:)を逆順 n2o(nvtx:1:-1)とする
CALL flip( n2o )!スタックオーバフローの回避:
END mk_RCM_table:（Program6：BS_tree_class に renumber を追加）
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

(31-6):PUBLIC :: renumber
INTERFACE :: renumber
MODULE PROCEDURE :: renumber_R
END INTERFACE:SUBROUTINE renumber_R ( this, o2n )TYPE(tree_c), INTENT(inout) :: thisINTEGER,INTENT(in) :: o2n(:)
TYPE(vtx), POINTER :: org_root
org_root => this%root !根をorg_rootに繋ぎかえ
this%root => NULL()!新しく木を茂らせる
CALL trvsl_post( org_root )CONTAINS
RECURSIVE SUBROUTINE trvsl_post( ptr ) !帰掛け
TYPE(vtx), POINTER, INTENT(inout) :: ptr
IF ( .NOT. ASSOCIATED(ptr) ) RETURN

Vol.19, No.4 2014

チュートリアル
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)CALL trvsl_post( ptr%L )CALL trvsl_post( ptr%R )
ptr%c = o2n(ptr%c)
ptr%L => NULL()
ptr%R => NULL()CALL insert_kernel( this%root, ptr )END SUBROUTINE trvsl_post
RECURSIVE SUBROUTINE insert_kernel( ptr, p_vtx )TYPE(vtx), POINTER, INTENT(inout) :: ptr
TYPE(vtx), POINTER, INTENT(inout) :: p_vtx
IF ( ASSOCIATED(ptr) ) THEN
IF ( p_vtx%c < ptr%c ) THENCALL insert_kernel( ptr%L, p_vtx )ELSE
CALL insert_kernel( ptr%R, p_vtx )END IF
ELSE
ptr => p_vtx
RETURNEND IF
END SUBROUTINE insert_kernel
END SUBROUTINE renumber_R:

この Program5 の第 13～21 行目では，RCM 法による全体自由度番号の付替え結果 n2o(:)を得て，その順引
き o2n(:)を作成しています．続く第 23～35 行目にて，この付替え結果を，隣接行列 adj_mtx_c 型の成分である隣接リスト表現 this%tree(:)と，節点番号から全体
自由 度 番 号 を 引 く 変 換 テ ー ブ ル node2DOFg(:,:)に反映させています．このとき，第 28 行目で 2 分探索木
tree_c 型同士の変数をコピーしていますが，構造型変数の成分に POINTER 属性がある場合には，Java や C な
どがそうであるように，ポインタ値のみのコピーです．
第 31 行目は Fortran2003 の組込みルーチンであり，ALLOCATABLE 属性の変数でも POINTER 属性の変数のよ
うに，割付メモリを指す先を変更するものです．この組込みルーチン MOVE_ALLOC により，メモリ再割付時のコピー操作が１回のみで済みます．この一連の操作は，隣接行列の行の入替えです．なお，列の入替えは第 29行目であり，続く Program6 が下請けします．
Program6 では，すでに茂っている 2 分探索木を，
n2o(:)による変換後の新しい番号にしたがってポイン
タを繋ぎ替え，新たな 2 分探索木に再構成しています．
具体的には，2 分探索木を「横断（traversal）」する帰り掛け（post-）に番号を付替え，新しい 2 分探索木のしかるべき位置に挿入していきます．これらは 2 つの内部再帰関数 trvsl_post と insert_kernel により実現
します．
Open-MP 並列化効率の数値実験
Open-MP 並列化効率のベンチマーク結果を示して，本チュートリアルを終えたいと思います．対象とする解析モデルは，初回に示した線形弾性問題であり，立方体形状の解析領域を 8 節点六面体要素で規則的に分割 し て い ま す ． 使 用 し た 計 算 環 境 は ， Intel Core
i7-5820K CPU (6 コア, 3.3GHz)，64GB メモリ，Windows
8.1，Intel Fortran Ver. 15.0，MKL Ver. 11.2 です．コン

５

計算工学

パ イ ル オ プ シ ョ ン は ， -fast -Qipo -Qopenmp-Qparallel -Qmkl です．すなわち，最適化を最大にして，Open-MP 並列化に加えて自動並列化も施しました．
図 3 に，マルチカラー化から求解までに要した実時
間を示します．PARDISO による直接求解では行列の正定値性を考慮しています．参考として，行列の対称性を反映せずに，対角スケーリング付き共役勾配法にて反復求解した場合も示します．全体として，6 コア時
で 1 コア時の 3 倍ほどの計算速度です．なお，立方体
領域を 50 3 要素で分割した場合（40 万自由度弱），MKL
PARDISO による直接求解で 5.5GB のメモリが必要でした．
時間の内訳に注目すると，まず PARDISO や自作反復ソルバの並列化効率がよくありません．メモリバン
ド幅が律速しているのかも知れません．しかし，今回の主題は Open-MP 並列化のための要素のマルチカラー化です．そこで，隣接行列 A_pp と係数行列 K_pp の組立てに注目すると，好ましい並列化効率を達成しています．特に，係数行列 K_pp の組立てはほぼ理想的な並列化効率です．CRS 形式における行列の成分位置(i,j)を，毎回 2 分探索するようにしたので処理時間が気に掛かっていましたが，杞憂でした．K_pp の組立て
は，非線形解析時に繰返し実行されること，反復求解自体の計算時間に比べて無視できるほどの時間でないことを考えると，マルチカラー化のメリットは大きいと考えます．
逆に，隣接行列 A_pp に基づき係数行列 K_pp を初期化 ― ダウンキャスト，内部的に「組立てモード」か
ら「計算モード」に移行 ― する場合のみ，幾分です
が並列化によりかえって時間を要するようになっています．この原因を調べてみると，Intel コンパイラのマ
ルチスレッドにおけるヒープメモリ管理に由来していました．すなわち，スレッド毎に独立なヒープメモリ
を設けているため，BS_tree_class で 2 分探索木を茂らせていく時には，メモリ割付の並列処理にほとんど問題ないものの，2 分探索木の終了処理時のメモリ解
放にしわ寄せがすべていくようです．2 分探索木のメ
モリ解放を止めると，このボトルネックは解消しました．しかし，ガーベージコレクションのことまで考えると，見掛けは不恰好でも，最大値をおさえて静的にメモリ確保しておくほうが良いのかもしれません．
６

おわりに計３回の本チュートリアルでは，有限要素計算における疎行列を圧縮行記憶（CRS）形式で作り上げる方法について，データ構造とアルゴリズムの視点から述べてきました．その核は，隣接行列という概念の切出
しと，その周辺を含めての実装でした．
なお，著者らの経験に基づき，コンパイラの最適化
にとって有利なプログラミング法程度まで言及しました．一方で情報処理の教科書 5) にありがちな形を示したいという意図から，抽象データ型プログラミング(31-7)

Vol.19, No.4 2014

チュートリアル

有限要素計算における全体剛性行列の作成法 ― 疎行列データ構造の視点から ― (3)（ADT）の考え方を壊すチューニングまでは述べませんでした．
本チュートリアルで示したプログラミング例について，より詳細なものを下記 URL に置いておきます．
http://www1.gifu-u.ac.jp/~gakuji/jsces2014/
ほぼすべてに詳細なコメントを入れてあります．ご笑覧いただければ幸いです．
参考文献
1) 船曳信生ら，グラフ理論の基礎と応用，共立出版，

20122) Welsh, D.J.A., Powell, M.B., An upper bound for thechromatic number of a graph and its application to
timetabling problems, The Computer Journal 10 (1):85-86, 1967
3) Reference Manual for Intel Math Kernel Library4) http://glaros.dtc.umn.edu/gkhome/views/metis
5) たとえば，石畑清, アルゴリズムとデータ構造, 岩
波書店, 1989

経過時間（秒）
0

20

40

60

80

100

120

要素分割数（コア数）32^3 (1)32^3 (6)50^3 (1)50^3 (6)50^3 (1)50^3 (6)
マルチカラー化
番号付替え
係数行列組立反復求解（対称性考慮せず）

隣接行列組立ダウンキャスト（隣接→係数行列）
直接求解（MKL PARDISO)

図 3 立方体領域の線形弾性問題の求解に要する実時間

計算工学(31-8)

Vol.19, No.4 2014