Hierarchical Diagonal Blocking形式における疎行列ベクトル積の高速化

全文

(1)情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. Hierarchical Diagonal Blocking 形式における疎行列ベクトル積の高速化室燎†1 田中輝雄†1 藤井昭宏†1 概要：Guys.E らによって Hierarchical Diagonal Blocking(HDB)形式を用いた SpMV が提案されている．HDB 形式は疎行列に対しパーティショニング，リオーダリング，ブロッキングを再帰的に繰り返し，疎行列を階層的にブロック化する．本研究では HDB 形式を用いた SpMV について 3 つの課題に取り組む．第 1 に HDB 形式のリオーダリング，ブロッキング，三角行列保持によるメモリ帯域の削減の作用がどのように性能に影響を与えているか詳細に評価されていない点に着目し，詳細な性能評価を行う．第 2 に HDB 形式の持つ階層構造は計算する階層が浅くなるにつれて並列度が低下する．この課題に対して，ワークベクトルを用いて並列度の低下を防ぐ手法を提案する．第 3 にもともと対角部分に多く非ゼロ要素ある疎行列に対してリオーダリングを行うと，非ゼロ要素が散らばりキャッシュ効率が低下する．この課題に対して浅い階層でのリオーダリングを無効化する手法を提案する．結果 OpenMP による並列化と負荷分散を行なった CRS 形式と比較し，提案手法を含めた HDB 形式を用いることで最大 1.70 倍の性能向上に成功した．キーワード：大規模数値計算，疎行列ベクトル積. 1. はじめに対称性を持つ疎行列の疎行列ベクトル積計算(SpMV)は，. く評価されていない点に着目し，HDB 形式の詳細な性能評価を行う． ② HDB 形式の持つ階層構造は SpMV を行うとき階層間. 大規模数値計算における反復解法に使用される[1]．疎行列. で部分的にベクトル𝒚への書き込みに依存関係があるため，. とは行列のほとんどの成分がゼロである行列である．疎行. 計算する階層が浅くなるにつれて並列度が低下する課題が. 列ベクトル積の計算の際にはゼロ要素は計算する必要がな. ある．. いため，ゼロを省いて格納する．. この課題に対して本研究では浅い階層の計算にワークベ. 疎行列の格納には広く CRS(Compressed Row Storage)形式. クトルを用いて，階層間の依存関係を切断する手法を提案. が用いられている．CRS 形式を用いた SpMV(𝒚 = 𝐴𝒙)の課. し，その影響について評価を行う．. 題は，ベクトル部が非連続的なアクセスとなるため CPU キ. ③ HDB 形式が行うリオーダリングの課題として，もとも. ャッシュミスが多発し性能が律速されることと，対称疎行. と対角部分に要素が集まっている疎行列に対してリオーダ. 列を並列化し計算するときには上三角行列と下三角行列の. リングを行うと，対角部分にあった非ゼロ要素が非対角部. 両方を保持しなければならないため，必要なロード回数が. 分に移動し，キャッシュ効率が下がる課題があげられる．. 増え，メモリ帯域を圧迫することである．. この課題に関して，本研究では浅い階層でのリオーダリ. この課題に対して HDB(Hierarchical Diagonal Blocking)形. ングを無効化しブロッキングのみを行う手法を提案し，評. 式を用いることが Guys.E らによって提案されている[2]．. 価する．. HDB 形式は疎行列に対してパーティショニング，リオー. 本論文の 2 章では，一般手法である CRS 形式と，本研究. ダリング，ブロッキングを再帰的に繰り返し，疎行列を階. で扱う HDB 形式の疎行列格納と SpMV の方法について解. 層的にブロック化して格納する．. 説を行う．3 章では，HDB 形式の持つ課題に対する提案手. パーティショニング(#1)とリオーダリング(#2)の作用に. 法について解説を行う．4 章では，リオーダリング，ブロ. より非ゼロ要素を対角に集約し，ブロッキング(#3)の作用. ッキング，三角行列化が及ぼす性能への影響を調査し，さ. により対角に集約された非ゼロ要素を階層的にブロック化. らに提案手法についての考察を行う．. する．#1，#2 の操作によりメモリアクセスの空間的局所性. 5 章結論では数値実験の結果をもとに，HDB 形式を用い. を高め，#3 の操作により時間的局所性を高める．. た SpMV の性能について結論を述べる．. HDB 形式ではブロック単位に並列化を行い，ブロックの計算は逐次的に行うため，対角を含む三角行列のみ保持し. 2. 疎行列格納方式. 計算が行えるため必要なロード回数が減りメモリ帯域の要. 2.1 CRS 形式. 求を削減することができる．. CRS 形式の概要を以下の図 1 に示す．CRS(Compressed. 本研究では，HDB 形式を用いた SpMV について，以下の. Row Storage)形式は行列の非零要素を格納する val[]，非零. 3 つの課題に取り組む．. 要素の列番号を格納する index[]，val[]と index[]の配列の行. ① HDB 形式のリオーダリング，ブロッキング，三角行列. 先頭となる要素の配列番号を格納する row_ptr[]の 3 つの配. 化の処理が SpMV の性能にどのように影響しているか詳し. 列で構成する．. †1 工学院大学 . ⓒ 2017 Information Processing Society of Japan. 1.

(2) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. val[]. セスの時間的局所性を向上させることができる．. 10 1 2 3 4 5 6 7 8 9 10 …. !=2,. $=2. index[] 0 3 7 1 2 5 1 2 4 5 0 …. row_ptr[]. 1 2. 1 1. 0 3 6 10 14 17 22 26 30 2 1. 図 1 CRS 形式の格納方法 2.2 CRS 形式の SpMV 対称疎行列を CRS 形式で計算する際には，2 通りの計算. 2 2. 2 3. 2 4. 図 4 HDB 形式の概略図. 方法がある．. 2.4 パーティショニング. 第 1 の方法は，対称疎行列の対称性を用いずに疎行列全. 図 5 にパーティショニングの概略図を示す．パーティシ. 体を保持し計算を行う．第 2 の方法は対称疎行列を対角を. ョニングでは対称疎行列に対してグラフ分割アルゴリズム. 含む下三角(または上三角)行列のみ保持し対称性を利用し. を適用し，疎行列の持つグラフ構造の依存関係をなるべく. 計算を行う．便宜的に，第 1 の方法を正方 CRS 形式，第 2. 切断せずに，疎行列を幅𝒘個の領域に分割する．. の方法を対称 CRS 形式と呼ぶこととする．. 本研究ではグラフ分割に metis ライブラリ [3] の. 正方 CRS 形式の SpMV の擬似コードを図 2，対称 CRS. METIS_Part_GraphKway 関数を用いた．. 形式の SpMV の擬似コードを図 3 に示す． for(i=0;i<rowsize;i++){ for(j=index[i];j<index[i+1];j++){ y[i] += A.val[j] * x[index[j]]; }} 図 2 正方 CRS 形式の SpMV 擬似コード for(i=0;i<rowsize;i++){ for(j=index[i];j<index[i+1]-1;j++){ y[i] += A.val[j] * x[index[j]];//as Upper y[index[j]] += A.val[j] * x[j];//as Lower } y[i] += A.val[j] * x[index[j]];//as Diagonal } 図 3 対称 CRS 形式の SpMV 擬似コード. 疎行列の行番号(列番号)がグラフの頂点に，非ゼロ要素. 正方 CRS 形式では疎行列全体を保持しなければならないため，疎行列の格納に必要なメモリ容量が増え，非ゼロ要素のロード回数が増加する欠点があるが，行単位での並列化が容易なこと，また，ベクトル𝒚への書き込みが連続的に行えるという利点がある．一方，対称 CRS 形式は，疎行列を下三角行列しか保持し. がグラフの辺に相当する．図 5 は疎行列を𝒘 = 2として疎行列を 2 つの領域に分割した例である．灰色の要素がグラフのエッジカットされた辺に相当する．. 0 1 2 3 4 5 6 7. 0 1 2 3 4 5 6 7. ! = 1, % = 2. 7. 0 5. 3. 1 6. 2. 4. 図 5 パーティショニングの概略図 METIS_Part_GraphKway 関数は各頂点が属する領域番号を格納した part[]を出力する．. ないため，正方 CRS 形式に比べメモリ容量とロード回数が. 2.5 リオーダリング. 減る利点があるが，行単位での並列化を行う際にはベクト. リオーダリングの概略図を図 6 に示す．リオーダリング. ル𝒚に書き込み競合が発生するほか，ベクトル𝒚へのアクセ. では，パーティショニングにより得られた part[]配列を用い. スが非連続的になる欠点がある．. て疎行列の行番号，列番号を並び替える．. 2.3 HDB 形式概要. この処理によって疎行列の非ゼロ要素を対角上に集める. HDB 形式の概略図を図 4 に示す．HDB 形式は疎行列に. ことができる．これにより，SpMV を行う際にベクトル部. 対してパーティショニング，リオーダリング，ブロッキン. へのアクセスの空間的局所性が上がり，性能が向上するこ. グを再帰的に繰り返し，疎行列を階層的にブロック化する．. とが期待できる．. HDB 形式における階層数を深さ𝒅，各階層のブロック分割. 2.6 ブロッキング. 数を幅𝒘とする．. ブロッキングの概略図を図 7 に示す．ブロッキングでは，. パーティショニングとリオーダリングによって疎行列の. リオーダリングにより対角に集まった非ゼロ要素を，幅𝒘. 非零要素を疎行列の対角上に集約する．これによってメモ. 個のブロック(小行列)に分割し，階層構造(木構造)化する．. リアクセスの空間的局所性を高める．ブロッキングでは対. ブロッキングにより，階層構造の葉に近いブロックは. 角に集まった非零要素を小さな疎行列(ブロック)として切. SpMV を行う際にアクセスするベクトル部の範囲が縮小し，. り出し木構造に結びつける．ブロック化によりメモリアク. メモリアクセスの時間的局所性の向上が期待できる．. ⓒ 2017 Information Processing Society of Japan. 2.

(3) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. 0 3 6 7 1 2 4 5. 0 1 2 3 4 5 6 7 0 3 6 7 1 2 4 5. 0 1 2 3 4 5 6 7. 3. HDB 形式の課題 3.1 キャッシュ効率と並列度による幅パラメータの制約 2.7 節で述べたように，HDB 形式の各ブロックは自らの子ブロックの計算が完了してから自らのブロックの計算を行うため，計算する階層が浅くなるにつれて並列度が低下する．そのため，並列度の観点からは浅い階層では 1 階層あたりの幅𝒘を大きくすることで並列度が低下しないようにすることが求められる．. 図 6 リオーダリングの概略図. しかしながら，1 階層あたりの分割数を増やすとパーティショニングでエッジカットされる辺が増えるため，浅い階層のブロックに残る非ゼロ要素数が増加する[5]．図 9 に深さを𝒅 = 1に固定し，幅を大きくして行ったときに，どのように根ブロックの非ゼロ要素数が増えるかを示す．横軸は幅𝒘で，縦軸は𝒘 = 2の時の非ゼロ要素数と比較した時に，根に残る非ゼロ要素数が何倍になっているか示. 1 1. 1 2. している．図 9 から幅パラメータを大きくしたときに親ブロックに取り残される非ゼロ要素数が増加しているのがわかる．. 図 7 ブロッキング概略図ィショニング，リオーダリング，ブロッキングを階層の深さパラメータ𝒅に達するまで再帰的に変換を行うことで， HDB 形式の階層構造を作り上げていく． 2.7 HDB 形式の SpMV HDB 形式を用いた SpMV は，1 ブロックの計算を 1 タスクとしてスレッドに処理を割り当てるタスク並列で並列化を行う．タスク並列の実装には Cilk Plus[4]を用いた． HDB 形式の同じ階層のブロックは書き込みを行うベクトル𝒚に対して重複した領域に書き込み処理を行わないため，並列に処理することができる．異なる階層のブロックに関しては，各ブロックの持つ子ブロックはベクトル𝒚の同一の領域に書き込み処理を行うため，書き込み競合を排除するために子ブロックの計算の完了を待ってからと自らのブロックの計算を行う． HDB 形式の SpMV の擬似コードを図 8 に示す．タスク並列の実装は擬似コード中の cilk_for で実現する．cilk_for は再帰的に呼び出される関数であっても並列処理することができる．cilk_for はループ内の処理が全て完了してから継続処理を行う． DoSpMV(TREE *node,x,y){ if(子がいたら){ cilk_for(i=0;i<子の数;i++){ DoSpMV(node->children[i],x,y); } } CRS_SpMV(node->matrix,x,y); } 図 8 HDB 形式の SpMV の擬似コード. 25 根ブロック⾮ゼロ要素数増加倍率. ブロッキングで作成された子ブロックに対してもパーテ. audikw_1. 20 bone010. 15. G3_circuit. 10. inline_1. 5 0. pwtk 2. 6 10 14 18 22 26 30 幅パラメータ. 図 9 深さ 1 で幅を大きくした時の非ゼロ要素数の変化浅い階層のブロックは，リオーダリングによる空間的局所性の向上やブロッキングによる時間的局所性の向上が見込めない．さらに根ブロックは逐次に計算が行われるため，深さ 1 の幅𝒘を大きくした時の根ブロックの計算量増加は，大きく性能を律速させる要因となる．そのため，キャッシュの利用効率の観点からは，浅い階層で幅𝒘を小さくすることが望まれ，並列度との両立が難しい． 3.2 ワークベクトルを用いた幅パラメータの制約の緩和本研究では，3.1 節で述べた浅い階層での幅パラメータに関する相反する 2 つの制約に対して，浅い階層のブロックの計算にワークベクトルを用いて階層間の依存関係を切断する手法を提案する．本手法は浅い階層のブロックの計算結果をベクトル𝒚とは異なるベクトル(ワークベクトル)に計算結果を書き込むことで書き込み競合をなくし，浅い階層のブロックは子ブロックの計算を待つことなく独立したタスクとして計算できるようにする．これにより，1 階層あたりの分割数を無理に大きくすることなく並列度の低下を防ぐことができる．. ⓒ 2017 Information Processing Society of Japan. 3.

(4) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. 本研究では，ワークベクトルを 1 本用いて根のブロック. の依存関係の切断について，③浅い階層でのリオーダリン. を独立したタスクにした場合と，ワークベクトルを 2 本用. グの無効化について，それぞれが SpMV の性能に対してど. いて根と深さ 1 のブロックの計算を独立したタスクにした. のように影響を与えているかを比較し，HDB 形式を評価す. 場合の 2 種類について実装を行なった．ワークベクトルを. る．本実験では，下記の表 1 に示す 8 種類の手法を用いる．. 1 本用いる場合の擬似コードを図 10 に，2 本用いる場合の. CRS①と CRS②はともに逐次実行で SpMV を行う．CRS. 擬似コードを図 11 に示す． Calculate(TREE *root,x,y,work_vector){ cilk_spawn CRS_SpMV(root->matrix,x,work_vector); cilk_for(i=0;i<子の数;i++){ DoSpMV(root->children[i],x,y); } cilk_sync(); reduction(y,work_vector); } 図 10 ワークベクトルを 1 本用いた SpMV 擬似コード Calculate(TREE *root,x,y,work_vector1,work_vector2){ cilk_spawn SpMV(root->matrix,x,work_vector); cilk_for(i=0;i<根の子の数;i++){ TREE *node = root->children[i] cilk_spawn CRS_SpMV(node->matrix,x,work_vector2); cilk_for(j=0;j<子の子の数;j++){ DoSpMV(node->children[j],x,y); } } cilk_sync(); reduction(y,work_vector1,work_vector2); } 図 11 ワークベクトルを 2 本用いた SpMV 擬似コード. ①は対称疎行列を対称 CRS 形式で保持し，CRS②は正方 CRS 形式で保持する． CRS③と CRS④はともに行単位で並列化をした正方 CRS 形式で，CRS④は HDB 形式が行なっている再帰的なリオーダリング処理を疎行列に対して適用する．CRS③と CRS ④はともにスレッドごとに計算する非ゼロ要素がなるべく均一になるように負荷分散を行なっている． HDB①と HDB②はともに全階層でリオーダリングを行なっている HDB 形式で，HDB①は各ブロックの疎行列を正方 CRS 形式で保持し，HDB②は対称 CRS 形式で保持している． HDB③と HDB④はともに各ブロックの疎行列を対称 CRS 形式で保持する．HDB③は 3.3 節で述べた浅い階層でリオーダリングを無効化する手法で，HDB④は全ての階層でリオーダリングを行わない手法である．表 1 計測手法一覧. cilk_spawn キーワードで呼び出す関数は，現在実行している関数と並列して処理を実行することができ，これにより浅い階層のブロックを独立したタスクとして計算するこ. D. 0. 2. 0. H. とができる．cilk_spawn でタスク化した関数は cilk_sync() を用いることで同期を行うことができる．浅い階層のブロックを独立したタスクとして計算し，ワークベクトルに書き込まれた計算結果は cilk_sync()を用いて同期処理を行いベクトル𝒚に計算結果を集約する．. B. C. 3.3 浅い階層でのリオーダリングもともと非ゼロ要素が対角部分に集まっている疎行列に対してリオーダリングを行うと，非ゼロ要素が非対角部分. 4.2 探索を行うパラメータについて. に移動してしまう．. HDB 形式①〜HDB 形式④と，HDB 形式の再帰的なリオ. 浅い階層でのリオーダリングによって非対角部分に移動. ーダリングを適用した CRS④について，計測を行うパラメ. した非ゼロ要素は，下の階層のリオーダリングに影響され，. ータの組み合わせを表 2 に示す．. より非対角部分で無秩序に非ゼロ要素が散らばってしまう. 深さパラメータを 1 としたときには 16 パターン，深さ. 傾向がある．. を 2 としたときには 16 パターン，深さを 3 としたときに. そこで本研究では，深さ0から𝒅までの階層構造のうち，. は 48 パターンとなり，合計 80 パターンの組み合わせを，. 𝒅/2階層まではリオーダリングを無効化しブロッキングの. それぞれの手法に対して計測を行なった．表 2 深さ・幅パラメータの探索範囲. みを行う手法を提案する．この手法により，もともと非ゼロ要素が対角上に集まっている疎行列に対してリオーダリ. 深さ 1. ングを行なった際に，浅い階層でブロックの非ゼロ要素が，無秩序に散らばってしまう可能性を小さくすることが期待できる．. 4. 数値実験. 幅1. 2,4,6,…,32. 幅1. 2,4,8,16. 幅2. 2,4,8,16. 深さ 3. 幅1. 2,4,8. 幅2. 2,4,8,16. 幅3. 2,4,8,16. 深さ 2. 4.3 実験環境本実験は東京大学情報基盤センターのスーパーコンピュ. 4.1 実験内容. ータ Reedbush-U[6]の 1 ノード利用して計測を行った．. 本実験では，①リオーダリング・ブロッキング・三角行. 現状ではプログラムを NUMA 構成に最適化していない. 列化の有無について，②ワークベクトルの用いた浅い階層. ⓒ 2017 Information Processing Society of Japan. 4.

(5) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. ため，numactl コマンドを用いて 2CPU のうちの片 CPU の. 表 5 性能計測結果. みを使用し計測を行った．実験環境の詳細を表 3 に，実験に使用する疎行列の詳細を表 4 に示す．今回使用した疎行列は SuiteSparseMatrixCollection [7]より引用した．表 3 実験環境 3GMBE BHG O. .9 .9. eh. [. 12. V. * .HKB. .9. 5 / M . ACB. 4- .HKB. 5 3GLMKNAMDHG . ACB. 4- .HKB. 5 . ACB. 4- .HKB. 5 . ACB. 6- .9. * 1- .9. [. a. * 1- L .9. i V X ]. 3GMBE . .H IDEBK. V X U. P.8 8 HIBG I. c. GN R. U. AME R7. 表 4 実験に使用する疎行列 G N 5 8 6. O. G. //3 102. G. 4. (. 9. ). ,. ) ,( ( ,). . 47 7. ,. 94. , (, ) (,. ). , ,. ) ((. ). 4.4 逐次実行の対称 CRS 形式と正方 CRS 形式の性能比較表 5 に各手法の性能計測結果を示す． CRS①(対称 CRS 形式)と CRS②(正方 CRS 形式)では，5 種類中 4 つの疎行列で CRS②の性能がより高かった．CRS ①の性能は CRS②と比較し，性能が向上したのが G3_circuit の 1.21 倍で，性能が最も性能が低下したのが pwtk の 0.625 倍であった．対称疎行列を三角行列で保持し対称性を利用し計算する CRS①は疎行列の非零要素のロード回数が減ることがメリットであるが，その代わりベクトル𝒚へのアクセスが非連続的になるデメリットがある．ブロッキングを行わない CRS 形式ではアクセスするベクトル𝒚の範囲が大きくなるため，キャッシュミスが発生しやすくなり性能が低下する. 4.6 HDB 形式における正方行列保持と三角行列保持. 疎行列が多いのだと考えられる．. CRS①(逐次実行・対称 CRS 形式)と CRS②(逐次実行・正. 4.5 CRS 形式におけるリオーダリングの性能への影響. 方 CRS 形式)の比較では，多くの場合で CRS①(逐次実行・. CRS③(HDB 形式の再帰的リオーダリングなし)と CRS④. 対称 CRS 形式)で性能が低下していたが，HDB 形式では. (HDB 形式の再帰的リオーダリングあり)の比較では，5 種. HDB①(ブロックを正方 CRS 形式で保持)よりも HDB②(ブ. 類中 4 種類の疎行列でリオーダリングによる性能向上が確. ロックを対称 CRS 形式で保持)の方が総じて性能が高くな. 認できた．. った．. リオーダリングによって非零要素を対角上に集めること. 三角行列で対称性を利用する SpMV は，ベクトル𝒚への. で，ベクトル𝒙へのアクセスでキャッシュの空間的局所性. アクセスが非連続になりキャッシュヒットミスが発生しや. が向上したためだと考えられる．ただし，全ての疎行列で. すくなるが，HDB 形式ではブロッキングを行うため，ア. 性能が向上している訳ではないことから，疎行列によって. クセスするベクトル部の範囲が小さくなり，ベクトル𝒚の. リオーダリングが逆効果になるケースがあるということが. ランダムアクセスによるキャッシュミスの確率が下がる．. わかる．. そのためキャッシュミスペナルティが減り，ロード回数が少なくなる効果と相まって HDB②で性能が向上している. ⓒ 2017 Information Processing Society of Japan. 5.

(6) 情報処理学会研究報告 IPSJ SIG Technical Report のだと考えられる．. Vol.2017-HPC-162 No.3 2017/12/18. 表 6 ワークベクトルを用いない各手法の性能計測結果. 4.7 HDB 形式における対称 CRS 形式保持の効果 S. HDB③は浅い階層でリオーダリングを無効化する手法であり，HDB④は全階層でリオーダリングを無効化する手法である． HDB②(全階層でリオーダリング)と比較して，bone010 は. G. H9. HDB③で 1.08 倍，pwtk は HDB④で 1.11 倍性能が高くなっている．この 2 つの疎行列はもともと非ゼロ要素が対角部分に集まっている疎行列である．. CB. 3.3 節で述べたように，対角部分に非ゼロ要素が集まっている疎行列に対してリオーダリングを行うと，もともと対角部分にあった非ゼロ要素の一部が非対角部分に移動してしまい，空間的局所性が低下する．そのため，そのような. 0 9. GF. 疎行列では一部，ないしは全体でリオーダリングを無効化しブロッキングのみを行うことが有効だと考えられる．逆に，元から非ゼロ要素が非対角部分に多い audikw_1,G3_circuit,inline_1 は全階層でリオーダリングを. B B 9. 行う HDB②で性能が高く，リオーダリングを一部ないしは全階層で行わない HDB③,HDB④で性能が低下している．これは非対角ブロックの行列に残る非ゼロ要素が増加し，キャッシュ効率の悪い非対角ブロックの計算量が増加した. O. LRP. DHF. ためと考えられる．. -56 1. 1. 1. 1. -56 1. 1. 1. 1. -56 1. 1. 1. 1. -56 1. 1. 1. 1. -56 1. 1. 1. 1.. 70 23468 ) ( ( (. (. ( ) ). ( ( )(. (. ( ) ( ( ). (. 4.8 ワークベクトルが性能に及ぼす影響. 素数を抑えつつブロッキングを行うことで性能が向上する. ワークベクトルを用いない各 HDB 形式の手法の性能を. ことがわかった．. 以下の表 6 に示す．表 6 からわかるように，ワークベクト. ワークベクトルの有無に関しては，本実験環境ではいず. ルを用いない HDB 形式は多くのケースで性能が CRS 形式. れの状況でもワークベクトルを用いた手法で性能が高く，. よりも下回っていたが，ワークベクトルを用いることで，. HDB 形式が課題としていた，浅い階層での並列度の低下. 5 種類中 4 種類の疎行列で CRS 形式よりも HDB 形式で性. に対して，ワークベクトルを用いた浅い階層間の依存関係. 能が上回り，全体的に性能が向上していることがわかる．. の切断は有効な手段であることがわかった．. CRS③と比較しワークベクトルを用いない HDB 形式は. 全体を通して，一般的な CRS 形式に比べ，提案手法を. 最大で 1.14 倍の高速化なのに対して，ワークベクトルを. 含めた HDB 形式を用いることで 5 種類中 4 種類の疎行列. 用いた HDB 形式は最大で 1.70 倍の高速化となり，浅い階層での並列度の低下に関して，ワークベクトルの使用が効果的であることがわかった．. で性能が向上し，最大 1.70 倍の高速化に成功した．. 参考文献. た．. [1] Barrett, R., et al.: Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM (1994) [2] Guy E. Blelloch, et al. , Hierarchical Diagonal Blocking and Precision Applied to Combinatorial Multigrid, Super Computing 2010 (2010) [3] METIS, http://glaros.dtc.umn.edu/gkhome/lab [4] 菅原清文：Cilk がやってきた，カットシステム(2011) [5]花上直樹，佐々木信一，菱沼利彰，藤井昭宏，田中輝雄： Hierarchical Diaonal Blocking を用いた疎行列ベクトル積の特性評価，情報処理学会第 77 回全国大会(2015)． [6] 東京大学情報基盤センタースーパーコンピューティング部門：Reedbush スーパーコンユーターシステム https://www.cc.u-tokyo.ac.jp/system/reedbush/ [7] SuiteSparseMatrixCollection, https://sparse.tamu.edu. リオーダリングに関しては，疎行列の形状に合わせてリ. 謝辞本研究の一部は，JSPS 科研費 15K15998 の助成を. オーダリングを行うか否かを選択し非対角部分の非ゼロ要. 受けたものです．. 5. 結論本研究では HDB 形式のリオーダリング，ブロッキング，三角行列化，ワークベクトルがどのように性能に影響しているかを網羅的に調査した．ブロッキングを行わない CRS 形式では，三角行列で疎行列を保持し対称性を利用し計算すると，ベクトル𝒚へのアクセスが非連続になるペナルティが大きく，結果として性能が低下していたが，HDB 形式ではキャッシュブロッキングの効果によりこのペナルティが軽減され，さらにロード回数が減ることにより性能が向上することがわかっ. ⓒ 2017 Information Processing Society of Japan. 6.

(7) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. 付録以下の表は、表 5 の各手法が計算を行った疎行列を示す。 CRS①. CRS①. CRS②,. CRS②,. CRS③,. CRS③,. HDB④. HDB④. CRS④. CRS④. 深さ:1. 深さ:3. 幅:32. 幅:2-16-4. audikw_1. bone010 HDB① 深さ:1 幅:10. HDB② 深さ:3 幅:4-4-4. HDB③ 深さ:2 幅:2-4. ⓒ 2017 Information Processing Society of Japan. HDB① 深さ:1 幅:12. HDB② 深さ:2 幅:4-8. HDB③ 深さ:2 幅:8-2. 7.

(8) 情報処理学会研究報告 IPSJ SIG Technical Report. G3_circuit. ⓒ 2017 Information Processing Society of Japan. Vol.2017-HPC-162 No.3 2017/12/18. CRS①. CRS①. CRS②,. CRS②,. CRS③,. CRS③,. HDB④. HDB④. CRS④. CRS④. 深さ:3. 深さ:2. 幅:2-8-2. 幅:8-4. HDB①. inline_1. HDB①. 深さ:2. 深さ:1. 幅:4-2. 幅:14. HDB②. HDB②. 深さ:2. 深さ:3. 幅:2-8. 幅:8-4-8. HDB③. HDB③. 深さ:2. 深さ:2. 幅:2-8. 幅:2-8. 8.

(9) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2017-HPC-162 No.3 2017/12/18. CRS①. CRS②, CRS③, HDB④. CRS④ 深さ:3 幅:4-8-2. pwtk. HDB① 深さ:1 幅:16. HDB② 深さ:1 幅:16. HDB③ 深さ:3 幅 :16-164. ⓒ 2017 Information Processing Society of Japan. 9.

(10)