有限要素計算における全体剛性行列の作成法- ―

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 DO
44:
!$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 DO
62:
!$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)
!下限 <=0
24:
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
⋯
deg
2) この降順にて，各頂点に彩色可能な最小の色識別番
号を割り当てる．
このとき，必要色数
の上限は次式で与えられます．
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 ):: this
INTEGER,
INTENT(out):: vtx2col(:) !(nvtx)
INTEGER
::
INTEGER,ALLOCATABLE::
INTEGER,ALLOCATABLE::
INTEGER
::
INTEGER,ALLOCATABLE::
LOGICAL,ALLOCATABLE::
INTEGER
::
nvtx
!全頂点数
vtx (:) !(nvtx)頂点番号
deg1(:) !(nvtx)頂点次数+1
deg1_i, max_deg1
buf (:) !取出用バッファ
used(:) !(0:max_deg1)色用
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 DO
CALL 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 DO
END 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_class
4:
:
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) :: this
26:
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%a
58:
& 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
b
c
d
e
f
g
e
d
f
g
a
b
c
a
b
c
d
e
f
g
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1 1
1
1
1
1
1
1
1
1
e
d
f
g
a
b
c
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
(3)
(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) :: this
INTEGER,
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(:)
!次数 +1
LOGICAL, 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 = 1
n2o(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 DO
DEALLOCATE (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) :: this
INTEGER,
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 ) THEN
CALL insert_kernel( ptr%L, p_vtx )
ELSE
CALL insert_kernel( ptr%R, p_vtx )
END IF
ELSE
ptr => p_vtx
RETURN
END 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) 船曳信生ら，グラフ理論の基礎と応用，共立出版，
2012
2) Welsh, D.J.A., Powell, M.B., An upper bound for the
chromatic 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 Library
4) 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

Download Report