九州大学学術情報リポジトリ
Kyushu University Institutional Repository
Orphan指示文を用いた粗粒度OpenMP並列プログラム
鈴木, 厚
九州大学大学院数理学研究院
https://doi.org/10.15017/1467657
出版情報：九州大学情報基盤センター広報 : 全国共同利用版. 6 (1), pp.1-6, 2006-05. Computing and Communications Center Kyushu University
バージョン：
権利関係：
Orphan 指示文を用いた粗粒度 OpenMP 並列プログラム
鈴木 厚∗
1 はじめに
偏微分方程式の数値計算では,差分法や有限要素法による離散化の後,未知自由度に関する連立方 程式を解くことになる. 時間発展問題や非線形問題ではこの連立方程式を繰り返し解く必要があり, 連立方程式の求解が総計算時間の大部分を占めることが多い. このため連立方程を高速に解くことが 求められる. 差分法あるいは有限要素法で得られる連立方程式の係数行列は疎行列(成分の値に 0が 非常に多い)であるという特徴を持つ. このためその解法にはCG法, GMRES 法などの Krylov部 分空間法[1] が用いられる. Krylov 部分空間法では行列・ベクトル積,ベクトル内積が主要な演算で あり,これらの演算を高速に実行することが必要である. 本稿では並列計算機, 特に共有メモリー型 の構成において OpenMPライブラリーを用いて行列・ベクトル積演算を並列計算する際のプログラ ミング手法について提案を行う.
2 行列・ベクトル積の計算
CG 法, GMRES 法などの主要演算である行列・ベクトル積演算の並列化を考える. 行列を A ∈
RN×N,ベクトルをu, v∈RN とする. 既知データu に対して v:=A u ⇔ [v]i:= X
1≤j≤N
[A]i j[u]j
を計算する. A が疎行列である場合を考えるので各i行に関してj 列の行列成分が非零の場合のみ, ベクトル配列のうちj で示される添字の値との積演算をおこなう:
[v]i := X
j∈{j;[A]i j6=0}
[A]i j[u]j.
ここで,疎行列A の記憶にはCRS フォーマット [1]を用い,メモリーの効率的な使用をめざす.
2.1 行列・ベクトル積の演算ループ並列化
OpenMPは共有メモリー型の並列計算機での並列計算を記述する. ループ(外側ループ)を並列化
する指示子を用いることにより次のように実装することが可能である [2, pp.39-40], [3, p.25].
Implementation (OpenMP, ループ並列化)
#pragma omp parallel for private(j) for (i = 0; i < N; i++) {
[v]i := X
j∈{j;[A]i j6=0}
[A]i j[u]j; }
この演算で,作業用変数j はループ変数i に依存するためprivate であることを明示する.
∗九州大学 大学院数理学研究院
E-mail: [email protected]
2.2 行列のブロック分割
以下では疎行列の性質を利用してブロック分解し,共有メモリー型,分散メモリー型並列計算機で 共通に利用できるアルゴリズムを考える. ただし,実装においては共有メモリー型並列計算機を用い て,ループ並列実装との比較を行う.
ベクトルの添字 Λ :={1,2, . . . , N}を
Λ = Λ1⊕Λ2⊕ · · · ⊕ΛD−1⊕ΛD
に分割する. ⊕は直和であり,添字部分集合に重なりはないとする. up ∈RN は添字Λp でのみ値を 持ち他の添字では 0 となるベクトルとする. 行列Ap q ∈RN×N は 1≤p < q ≤D−1 なるp, q に 関しては[Ap q]i j= 0 (i∈Λp, j ∈Λq)を満たすものとする. この添字の分解により行列・ベクトル
積は
v1 v2 ... vD−1
vD
=
A1 1 A1D
A2 2 A2D
. .. ...
AD−1D−1 AD−1D AD1 AD2 · · · AD D−1 AD D
u1 u2 ... uD−1
uD
となる.
ΛD の添字に関するデータを更に分割する:
Λ(p)D :={j∈ΛD; [Ap D]i j6= 0 ∀i∈Λp} Λ(p)D は添字集合 ΛD の‘重なりのある’分割になる: ΛD = S
1≤p<DΛ(p)D . 行列 A(p)D D を AD D = P
1≤p<DA(p)D D を満たすAD D の部分添字 Λ(p)D への分割とする. これらの分割を用いると行列・ベクトル積は
A u= X
1≤p<D
"
Ap p Ap D AD p A(p)D D
# "
up u(p)D
#
となる.
2.3 領域分割による実現
係数行列 A を有限要素法の剛性行列とするとき,ベクトル添字の分解は‘重なりのない’領域分割 から構成することができる. 領域Ω をD−1 個の重なりのない部分領域Ωp(1≤p < D) に分割す る. 分割によって生じた内部境界を F と表す: F :=S
1≤p<q<D∂Ωp∩∂Ωq. 節点全体の集合は次の 分割に対応して D 個の部分に分割される:
Ω =¯ [
1≤p<D
( ¯Ωp\ F)∪ F.
剛性行列 A の分解は部分領域Ωp での領域積分で計算され,それぞれ次の節点の関係からなる: Ap p Ω¯p\ F の節点
Ap D Ω¯p\ F と∂Ωp∩ F の節点 A(p)D D ∂Ωp∩ F の節点
2.4 並列計算アルゴリズム
行列・ベクトル積の計算手順は次のようになる: Algorithm (ブロック並列化)
Step 1
プロセッサー p 毎に up, u(p)D を準備する.
( u(p)D はuD の部分ベクトルである. 共有メモリー型の場合,データ転送は必要ない. ) Step 2
プロセッサー p 毎に部分行列と部分ベクトルの積を計算する:
"
vp v(p)D
# :=
"
Ap p Ap D AD p A(p)D D
# "
up u(p)D
# .
Step 3
添字集合 ΛD に関する和をとる:
vD := X
1≤p<D
v(p)D .
ここで vp (1≤p < D) は他の値と足し合わせる必要はなく, Step 2が完了した時点で計算は完 了していることに注意する.
Step 3 は添字集合 ΛD の全ての成分で計算を実行する必要があるが,この部分の並列化のため添
字の分割を導入する.
ΛD = ΛD,1⊕ΛD,2⊕ · · · ⊕ΛD,D−1 なる D−1 個の添字の直和分解を用意する.
Step 30
プロセッサー p 毎に次を実行する:
[vD]j := X
1≤q<D
[vD(q)]j j∈ΛD, p.
Step 30 は全てのプロセッサーで Step 2が完了してから実行することになるため, Step 30 の直前で 同期をとる必要がある. 共有メモリー型の場合, Step 2とStep 30 の間には同期操作が必要なだけで, データ転送は必要ない.
内積演算はΛD の分割を用いると次のように計算できる: (u, v) = X
1≤p<D
(up, vp) + X
1≤q<D
(uD,q, vD,q).
2.5 OpenMP による実装
これらの計算法は分散メモリー型計算機にも実装することができるが,共有メモリー型でOpenMP を用いて記述すると容易に実装できる. 2.1節に述べたOpenMPのループ並列機能ではなく, orphaned 指示文 [4, pp.4,9] を用いる.
Implementation (OpenMP, orphaned 指示文)
#include <omp.h>
配列 uD,vD を確保; omp_set_dynamic(0);
omp_set_num_threads(D−1);
#pragma omp parallel default(shared) private (p) {
p=omp_get_thread_num();
配列 up,vp,u(p)D ,v(p)D を確保;
uD のデータを作業配列u(p)D へコピー;
"
vp vD(p)
# :=
"
Ap p Ap D AD p A(p)D D
# "
up u(p)D
#
;
#pragma omp barrier [vD]j :=P
1≤q<D[v(q)D ]j j∈ΛD, p; }
ここで, OpenMP の指示文 #pragma omp parallel { } で囲まれた並列リージョンが並列実行 される. プロセッサー番号p には OpenMPのスレッド番号を割り当てている.
配列 uD, vD は並列リージョンの外でメモリーを確保しているため,すべてのスレッドからアクセ スができる (shared) が, 配列 up, vp, u(p)D , v(p)D は並列リージョンの中でメモリーを確保しているた め, 他のスレッドからはアクセスできない (private). OpenMP ライブラリーの実装のうち多くの ものでは並列リージョンの中でメモリーを確保すると,プロセッサーpのローカルメモリーに確保さ れるため,高速な演算が期待される.
CG法などの計算手続きは並列リージョンの内部に記述する. スレッドの生成消滅は並列リージョ ン の中ではおこらないため,並列処理に関するオーバーヘッドは barrier文(orphaned指示文)のみ である. このように, OpenMP を用いても粗粒度の並列計算が実現できる.
プロセッサー p 毎での部分行列と部分ベクトルの積の演算においては添え字ΛD に ‘重なりのあ
る’分割 ∪1≤p<DΛ(p)D を行っているため,並列化前のアルゴリズムあるいはループ並列化のアルゴリ
ズムに比べて計算量が増えている.
2.6 数値実験結果
表1に示すサイズの有限要素法剛性行列を用いて,行列・ベクトル積の2種類のアルゴリズムの並 列計算効率を計測した. アルゴリズムの実行に用いた並列計算機は九州大学 情報基盤センター設備 の高性能演算サーバー IBM eSerever p5 モデル 595 (kyu-cc) である. この計算機は64 基の CPU からなる共有メモリー型計算機である. ここで領域分割数はD−1 = 48に固定し, プロセッサー数 はその約数に選んだ. 表2 に経過時間, MFLOPS値を示す. 図1 にMFLOPS値を両対数グラフを 用いて表し,並列効率を示す.
表1. 疎行列サイズ
ベクトル自由度 非零要素数 浮動小数点演算数(MFLOP)
orphaned 指示文 ループ並列
22,932,588 524,436,522 2,124.846 2,074.813
表2. 行列・ベクトル積演算の計算効率
orphaned 指示文 ループ並列
# CPU 経過時間(秒) MFLOPS 経過時間(秒) MFLOPS
1 4.27348 497.22 4.21108 492.70
2 2.11207 1,006.05 2.10430 985.99
3 1.45398 1,461.40 1.44805 1,432.83
4 1.06202 2,000.75 1.08669 1,909.29
6 0.706038 3,009.53 0.705740 2,939.92 8 0.535862 3,965.28 0.541420 3,832.20 12 0.367691 5,778.89 0.372019 5,577.17 16 0.278528 7,628.83 0.283959 7,306.78 24 0.196061 10,837.67 0.191501 10,835.97 48 0.115946 18,326.17 0.111125 18,670.98
1,000 10,000
1 2 3 4 6 8 12 16 24 48
MFLOPS
n = # CPU
orphaned directive loop parallelization O(n)
図1. 行列・ベクトル積演算の計算効率
図1 からorphaned指示文を用いた実装,ループ並列を用いた実装ともに並列演算効率が高いことが
わかる. orphaned指示文を用いた実装はループ並列を用いた実装に比較して総演算量が増えている
ため 1 CPUでの実行は遅くなっている. しかし 表 2より 2から 24 CPUでの実行ではFLOPS値 は大きく,実行時間はわずかではあるが短くなる場合が多いことがわかる.
3 まとめと今後の展望
疎行列とベクトルの積演算について,共有メモリー型並列計算機環境で外側ループを並列化するも のと領域分割を利用してブロック並列化するものの二種を比較した. 後者のアルゴリズムは粗粒度の
並列計算であり,分散メモリー環境向けであると考えられているが,共有メモリー環境においても有 効な実装法であることが数値実験により確認された.
共有メモリー型並列計算機を複数台結合した SMP クラスター並列計算環境では分散環境向けの MPI と共有環境向けの OpenMP を混在して用いることが多いが, 領域分割とループ並列を混在す るとプログラミングが非常に複雑になる. 本稿で述べた領域分割を利用した並列計算アルゴリズムは
MPI とOpenMP の双方に実装できるため両者の混在したプログラミング環境にも親和性が高いと
予想される. SMP クラスター環境への実装と性能評価が今後の課題である.
参考文献
[1] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der Vorst.Templates for the Solution of Linear Systems : Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.
[2] 南里 豪志, 渡部 善隆, OpenMP 入門(2), 九州大学情報基盤センター広報, Vol.2, No.1, pp.1 - 40, 2002.
[3] 南里 豪志,高性能演算サーバにおける並列処理 自動並列化, OpenMP, MPI,九州大学情報基盤 センター広報, Vol.5, No.1, pp.1 - 40, 2005.
[4] OpenMP C and C++ Application Program Interface Version 1.0 – October 1998, Document Number 004-2229-001,http://www.openmp.org.
