close

ログイン

パスワードの再発行?

OpenIDでログイン

大規模疎行列演算と 3次元有限要素並列計算 - 九州大学

エンベッド
大規模疎行列演算と3次元有限要素並列計算
大規模疎行列演算と3次元有限要素並列計算
鈴木厚*
1 はじめに
有限要素法[2, 9】は偏微分方程式の離散化手法の一つであり,領域の形状に対する汎舟性があ
ること,数学理論が整備され数値解に対する誤差評価を持つスキームを用いることができること
などの長所がある.筆者らは地球マントル対流問題を記述する無限プラントル数を持つブシネ
スク流体の熱対流の有限要素スキームと3次元計算コードの開発を行ってきた匝, 6, 7, 8].
3次元問題に対する有限要素法では,債域を四面体に分割し離散化を行うことで,節点での未
知量に関する連立方程式を得る.六面体を用いた要素分割も可能であるが,要素分割が容易であ
ること,演算量が少ないことを考慮して四面体分割を採用する.この連立方程式を構成する行列
は疎行列である.有限要素法では,行列データに規則性の無い非構造データを扱うことになる.
クリロフ部分空間法[1, 3]などの反復解法を用いて連立方程式を解く.この解法では行列と
ベクトルの積計算と,ベクトルの内積計算を繰り返し行う.本稿では,ワークステーション型の
RISCプロセッサーの大規模疎行列とベクトルの積演算の性能を調べることと,キャッシュメモ
リーの利用と並列化を行い演算の高速化を図ることを目的とする.
2 流れ問題に現れる有限要素行列
有限要素法では支配方程式にテスト関数を乗じ,対象とする領域で積分を行い,弱形式を得る.
さらにこの弱形式を離散化することで,剛性行列からなる有限要素方程式を得る.本稿では, 3
次元の流れ問題を考える際に必要な次の剛性行列を取り扱う・ Wih≦i<Nを流速の有限要素基
底回とする.流速の自由度Ⅳは節点自由度の3倍である. Ⅳ×Ⅳ行列Aを次のものする.
Wij蝣蝣-2
∑ Dki(ipj)Dkl((pi)dx (0 ≦ ij < N). (1)
/n 1≦k,l≦3
ここにD(u)は変形速度テンソルDuiu) - (∂uk/∂xi4-∂ul/∂Ek)/2 (1 ≦ k,I ≦ 3)である.こ
の行列は,基底関数が局所的であることから疎行列となり,要素分割に自由度があることから一
般に零成分の分布に規則性のない非構造データとなる.また, Aは対称行列であることに注意す
る.
以下,行列と数ベクトルのインデックスはC言語の慣例に従い0から数えるものとする.ベ
クトル百のインデックスiで示される値をbli,行列Aのi行j列の成分を[A]i]と書く・
3 疎行列の記憶
CRS (Compressed Row Storage)フォーマット国を用いて,疎行列を記憶する.このフォー
マットは,各行の零成分を取り除いて圧縮し,行列の非零成分のみを記憶するため,省メモリー
*九州大学大学院数理学研究院E-mail: asuzuki馳ath.kyushu-u.ac.jp
135-
九州大学情報基盤センター広報
Vol. 1 No. 3 2∝)1
研 究 開 発
である. nxを行列Aの行の数,nAを行列Aの非零要素の数とする. CRSフォーマットでは
行列Aの値を格納する長さnAの浮動小数点ベクトルV-A ∈RnA,列のインデックスを記憶する
長さnxの整数ベクトルC-A ∈ Znxと,各行での列の開始位置を記憶する長さnx+lの整数ベ
クトルr-A∈%nx+lを用いる・ Aのk行l列にある非零要素[A]klを配列屯のインデックス
βの場所fedβに格納する・この時C-Aのインデックスβでの値3dβはlを記憶している蝣rA
は各行で列を開始するインデックスを記憶するため> [^]; ≦β< Rdk+iである.ただし, rAの
n*+l番目の値Rd"xはnAとする.
行列Aが対称行列の場合,下三角部分と対称部分のみを記憶すれば良いことに注意する.
4 行列ベクトル積演算と正射影演算
図1 (左)にCRSフォーマットによるデータ構造(C言語の構造体による記述)を示す. SIZE
はベクトル長さ,すなわち行列の行数nxを表しROWSlはr-Aの長さを表す.またNONZEORS
は行列の非零要素の数nAを表す.説明の簡単化のため,これらのサイズを静的に決定するが,
実際の有限要素プログラムコードでは,プログラム中で動的に決定する.これは次の理由による.
ベクトル長さは自由度数,非零要素数は自由度数と要素分割に依存する.要素分割に関するデー
タを入力データとして与えることで,領域形状に対するプログラムの汎用性を実環することがで
きる.
for (i=0; iくSIZE; i++) i
#define S工ZE /****/ v_tmp = 0.0;
#define ROVSl (SIZE + 1) j_begin = aa.row[i] ;
#define NONZEROS /****/
j_end= aa.rowLi + 1] - 1;
for (j = j.begin; jくj_end; j++) {
typedef struct {
j_tmp = aa.index[j] ;
int row[ROWSl] ;
v_tmp += aa.value[j] * u[j_tmp] ;
int index[NONZEROS] ;
v[j_tmp] += aa.value[j] * u[i];
I
double value [NONZEROS] ;
> sparse_matrix;
j_tmp = aa.index[j_end] ;
v_tmp += aa.value[j_end] * u[j_tmp] ;
v[i] += v_tmp;
sparsejnatrix aa;
)
図1:疎行列データ構造(左)と行列ベクトル積(右)
図1・ (右)に対称行列Aとベクトルu-の積演算を行い,ベクトル百を得る手続き
v:=Au
を示す・この演算ではベクトルtf(uD), tT(vD)の要素のアクセスに,間接参照が必要である.
これは,各行の零成分を取り除いて圧縮して成分を記憶することから生じている.行列が対称で
あることを利用し,間接参照のインデックスjに関するforループで行列の下三角と上三角部
九州大学情報基盤センター広報
Vol. 1 No.3 2001
-136-
大規模疎行列演算と3次元有限要素並列計算
S-tmp = 0.0;
for (i=0; iくSIZE; i++) -C
s_tmp += x[i] * yLi];
I
for (i=0; iくSIZE; i++) {
Ⅹ[i】 -= s_tmp * y[i];
)
図2:正射影演算
分の積を計算した後,対角部分の積の計算を行っている.ここで,配列Ⅴ[]は0に初期化されて
いるものとする.
図2.にベクトルの正射影演算を示す.ベクトルE ∈Rnxから長さが1に規格化されたベク
トルy-∈ョーnx ((y,y) - i)の成分を除去する演算
x:-x-(x,y)y
である・ここで>%y)はベクトル3;-とy-の内積である
(x,y):- ∑ [x]iffli蝣
0≦i<nx
正射影演算は,非圧縮流れ問題において,圧力の自由度を除去する時などに必要な演算である.
表1, 2に数値実験に使用した計算機とCコンパイラオプションをそれぞれ示す.
表3に図1 (右)に示される行列ベクトル積の実行結果,表4に図2に示される正射影演
算の実行結果を示す.両者とも1 CPUにより実行した.表3の結果は,図7に示される球殻
領域の有限要素分割と離散化によって得られた流速の剛性行列(1)に対するものである.ベク
トルの長さnx (SIZE)は973,956であり,行列Aの下三角部分と対角部分の非零要素の数nA
(NONZEROS)は21,909,498である.行列Aの一行あたりの平均的な非零要素の数は23 × 2程
度であり,行列Aの非零成分は要素全体(nx x nx個)に対して(0.0023× 2)%程度である.
表3,表4とも演算量(Flop)はSR8000/MPPのハードウェア機能により計測したものであ
り, MFLOPS値は(演算量)/(実行時間)により算出した.
SR8000/MPPの1 CPUでの理論性能は1.8GFLOPSであるが,行列ベクトル積の演算では
10%程度,正射影演算は22%程度の性能になっている. SR8000/MPPのベクトル化機能では
図1 (右)のコードは十分にベクトル化されないため,スカラー型のCPUとして利用している.
ワークステーション型の計算機GP7000F, GS320では主記憶-のアクセス速度が演算の制約に
なっていると考えられる.
CRSフォーマットはデータ構造が単純であり,省メモリーの観点から優れた記憶方式である
が,計算効率は高くない.ベクトル化に向いたJDSフォーマット国があるが,本稿ではスカ
ラー型のRISCプロセッサーでの速度向上を目指す.
-137-
九州大学情報基盤センター広報
Vol. 1 No. 3 2001
研 究 開 発
表1:数値実験で使用した計算機
計算機名 CPU仕様 設置場所
GP7000F SPARC 64 九州大学情報基盤センター
GS320 Alpha 21264 (731MHz)九州大学情報基盤センター
SR8000/MPP RISCプロセッサー 東京大学情報基盤センター
(1CPUあたりの
理論性能1.8GFLOPS)
表2:コンパイラオプション
GP7000F fee -Kfast -Kfast_GP=2 -Kprefetch -KV9FMADD
GS320 cc -fast -04 -arch host -tune host -inline speed -non_shared
SR8000/MPP cc -Os -64 -restrict -noparallel -prefetch -nopreload
表3:行列ベクトル積{nx - 973,596, nA - 21,909,498,演算量88.612M Flop)
計算機名 実行時間(秒) MFIJOPS
GP7000F 2.03
GS320 0.446
SR8000/MPP 0.490
表4:正射影(n* - 973,596,演算量3.894M Flop)
計算機名 実行時間(秒) MFLOPS
GP7000F 0.0936
GS320 0.0330
SR8000/MPP 0.00962
5 高速化
5.1 合同な部分領域への分割
計算領域の特徴を活用し,合同な部分領域-の分割を行う.領域は〟個の重りなのない部分
領域Q<m> (0 ≦m < M)に分割されているとし,それぞれ,ある直交変換Rim)で参照領域n<->
を写したものであるとする.また,部分領域での要素分割も参照領域での要素分割をR(m)で写
したものであるとする.
成分が-i, O, iのみからなる直交変換RSm>に対して,次の性質が成り立つ.
九州大学情報基盤センター広報
Vol. 1 No.3 2001
- 138-
大規模疎行列演算と3次元有限要素並列計算
定理1 A をQMでの流速の剛性行列とする.直交変換の符号より定まる,対角要素が-1
と1からなる対角行列β(m)が存在して,
J^rn) - g{m)j^(O)g(m)
(2)
が成り立つ(証明は[7]を参照).
この定理より,流速の剛性行列は領域全体で作成すること無く,参照領域fi<->でのみ作成し記憶
すれば良いことになり,大幅なメモリー節約が可能になる(詳細は[7]を参照).領域分割を並列
計算に利用することにより計算の高速化を実現することができる.
図7に球殻領域の48部分領域-の領域分割を示す.ここで,領域分割と要素分割が見易いよ
うに,参照領域を取り除き,各部分領域をシフトしている.
5.2 キャッシュメモリーの活用とOpenMPによる並列化
式(2)を利用して行列ベクトル積の演壷を行う.領域の分割によって生じた内部境界上のデー
タを足し合わせる操作が必要であるが,部分領域毎に独立して次の計算を行うことができる.
g(m).-S{m)^m) (q≦m<M),
〆m) :-A(0)」(m) (o≦m<M),
〆m).-g(m)〆m) (o≦m<M).
M ( SUBDOMAINS )個の領域に対してP ( PROCS )個のCPUを使用して計算を行う・ PをM
の約数になるように選ぶと,各CPUはM/P (SIZE2)個の部分領域を扱うことなり,負荷分散
の良い並列計算が可能である.図3, 4に(4)の疎行列ベクトル積に対応する演算を示す・図3
は,領域分割のループを行列ベクトル積の外側に置いたもの,図4は,再内側に領域分割のルー
プを置くようにループの入れ替えと配列の記憶方法の変更を行ったものである. idはCPUの
インデックスを表す( 0 ≦id<PROCS ).また, S工ZElは部分領域での自由度を表す・
並列化はOpenMP国を用いて行った. loop 1, loop 2ともにOpenMPの並列実行領域の中
に記述され,インデックスidに対する並列化が行われる. OpenMPでは変数と配列に対して,
スレッド間で共有するかしないかの指示が必要であるが,ここではOpenMPによる並列化コー
ドの詳細は省略する.
図5と図6に部分領域数を〟-48としたときの,全体領域の行列とベクトルの積の演算
(内部境界での足し合わ漣処理を含む)に対するGP7000FとGS320での並列化効率を示す.
loop2では,行列の値aa.value[j]がインデックスqに関するループで再利用されるため,
キャッシュメモリーを有効に利用していると推測される. 1 CPUでの実行時, Ioop2はIoopl
に比較し, GP7000Fでは1.73倍, GS320では1.42倍の速度向上が得られている.
また,図4からGP7000Fでは,単体CPUあたりの演算能力はGS320に比べ低いが, OpenMP
による並列化が効率良く行われることがわかる. Ioop2で6 CPUまでは,ほぼCPU台数分の
速度向上が得られている. 24 CPU使用時の速度向上は13.2倍である.
6 おわりに
ワークステーション型のRISC計算機による疎行列ベクトル演算は,主記憶-のランダムな
間接アドレス参照を含むことなどから,ピーク性能を実現できない.しかし,行列データにある
-139-
九州大学情報基盤センター広報
Vol. 1 No. 3 2001
研 究 開 発
SIZE2 = SUBDOMAINS / PROCS;
for (q=0; qくSIZE2; q++) {
p=
id*S工ZE2+q;
for (iー=0; iくSIZEl; i++) i
v_tmp = 0.0;
j_begin - aa.row[ij ;
j_end= aa.row[i + 1】 - 1;
for (j = j.begin; jくj_end; j++) {
j_tmp = aa.index[j] ;
v_tmp += aa.value【j】 * u[p] [j_tmp] ;
Ⅴ[pHj_tmp] += aa.value[j] * u[p] [i];
I
j_tmp = aa.index[j_end] ;
v_tmp += aa.value[j_end] * u[p] [j_tmp] ;
v[p] [i] += v_tmp;
)
I
図3:部分領域ループを外側に持つIoop 1
SIZE2 = SUBDOMAINS / PROCS;
for (i=0; iくSIZEl; i++) i
for (q=0; qくSIZE2; q++) {
p=id*SIZE2+q;
v_tmp[p】 0.0;
I
j_begin = aa.row[i] ;
j_end = aa.row[i + 1] - 1;
for (j = j_begin; jくj_end; j++) 1
j_tmp = aa.index[j] ;
for (q=0; q<SIZE2; q++) {
p=id*S工ZE2+q;
v_tmp[p] += aa.value[j] * u[j_tmp] [p] ;
Ⅴ【j_tmp] [pコ+= aa.value[j] * u【i] Cp] ;
I
I
j_tmp = aa.index【j_endコ;
for (q=0; q<SIZE2; q++) {
p=id*SIZE2+q;
v_tmp【p] += aa.value【j_end] * u【j-tmp] 【p] ;
Ⅴ[iコ[pコ+= v_tmp【pコ;
)
)
図4:部分領域ループを再内側に持つIoop 2
九州大学情報基盤センター広報
Vol. 1 No.3 2001
-140-
大規模疎行列演算と3次元有限要素並列計算
ni
('08S)9∈!1
3 4 6 8 12 16 24
P : numberofCPUs
図5: GP7000Fでの領域分割計算の効率
■l
("08S)9∈!1
3 4 6 8 12 16 24
P : numberofCPUs
図6: GS320での領域分割計算の効率
-141-
九州大学情報基盤センター広報
Vol. 1 No.3 2001
研 究 開 発
種の構造があり,データの再利用を用いてキャッシュメモリーを利用できる場合には,演算速度
向上が得られることを確認した.また, OpenMPを用いた並列計算によっても高速化が実現で
きた.
共有メモリー型計算機とOpenMPを利用する並列計算は分散メモリー型の計算機を用いる
場合と比較し,並列コード-の書き換えとデバッグが比較的容易である. OpenMPによる並列
コードの実行性能の高い計算機が望まれる.
本研究は,九州大学大学院数理学研究院田端正久教授との共同研究の成果の一部であり,文
部科学省科学技術振興調整費「全地球ダイナミクス」の援助を受けた.
図7:球殻領域の合同な部分領域-の分割
参考文献
[1] R. Barett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout,
R. Pozo, C. Romine, and H.Van der Vorst, Templates for the Solution of Linear Systems:
Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.長谷川里美,長谷川秀
彦,藤野清次訳,応用数値計算ライブラリ反復法Templates,朝倉書店1996.
2 P. G. Ciarlet, The Finite ElemenまMethod for Elliptic Problems, North-Holland, Amster-
dam, 1978.
九州大学情報基盤センター広報
Vol. 1 No.3 2∝)1
- 142
__」
大規模疎行列演算と3次元有限要素並列計算
[3] G. H. Golub and C. F. Van Loan, Matrix Computations, 3rd ed. The John Hopkins
University Press, Baltimore, 1996.
[4] OpenMP C and C++ Application Program Interface Version 1.0 - Octorber 1998,
http : //www. openmp. org.
囲A. Suzuki, M. Tabata, and S. Honda, Numerical solution of an unsteady Earth's mantle
convection by a finite elementーmethod, Theoretical and Applied Mechanics, 48 (1999),
371-378.
[6] M. Tabata and A. Suzuki, A stabilized finite element method for the Rayleigh-Benard
equations with infinite Prandtl number in a spherical shell, Comp. Meth. Appl. Mech.
Engrg., 190 (2000), 387-402.
円鈴木厚,田端正久,合同な部分領域-の分割による3次元球殻内熱対流問題の並列計算,計
算工学講演会論文集, 5 (2000年5月), 165-168.
[8] M. Tabata and A. Suzuki, Mathematical modeling and numerical simulation of Earth's
mantle convection, to appear in Proceedings of the International Symposium on Mathematical Modeling and Numerical Simulation in Continuum Mechanics, Lecture Notes in
Computational Science and Engineering, P. G. Ciarlet I. Babuska and T. Miyoshi (eds),
Springer.
[9] 0. C. Zienkiewicz and F. C. Scott, On the principle of repeatability and its application
in analysis of turbine and pump impellers, Int. J. Num. Meth. Eng., 4 (1972) 445-452.
-143-
九州大学情報基盤センター広報
Vol. 1 No. 3 2001
Document
カテゴリー
SOFTWARE AND INTERNET
本日のアクセス数
21
ファイルサイズ
1 621 KB
タグ
1/--ページ
報告