性能検証

この節の最後に，Bečka の動的順序を⽤いたときの性能測定結果を⽰す．ここで⽤いたプログラム では，重みの計算を，単純なフロベニウスノルムの近似ではなく，Bečka が⽤いている，列ブロック を正規化したものに対する近似を⾏っている．また，重みは近似による誤差によって⼀定以上0に近 づかないため，tolの値を，巡回順序を⽤いたときの値√𝑛uから 𝑛𝑙 uに変更している．図 5.20と 図 5.21は，前節と同様に SGI Altix ICE を⽤いて⽚側ブロックヤコビ法を実⾏したときの，実⾏時間 の内訳である．前節の結果と違い，動的順序を⽤いているため，重み計算や貪欲法によるマッチング，

ハンガリアン法による割り当ての時間が追加されている．そのため，𝑘 = 1のときにおける実⾏時間

1 2 4 6 8 12

timeinsec.

prepost col-exchange gram cholesky small-svd tri-solve matmul wc greedy-pair hungarian barrier others

図5.21 SGI Altix ICE の 144CPU を⽤いて Bečka の動的順序を⽤いた⽚側ブロックヤコビ法の実

⾏時間を測定したときの，実⾏時間の内訳．⾏列は⼤きさが , となっていること以 外，図 5.14と同じ． の結果は左軸，それ以外のものは右軸を⽤いている．

表5.6 Modified Modulus 順序と Bečka の動的順序の場合における縦⽅向の並列数 と反復回数

⾏列サイズ (𝑛) ⼿法 𝑘 = 1 𝑘 = 2 𝑘 = 4 𝑘 = 6 𝑘 = 8 𝑘 = 12 4, 320 MMO 2, 016 1, 008 432 288 180 120 Bečka 961 464 225 145 106 70 8, 640 MMO 2, 016 1, 008 432 288 180 120 Bečka 966 463 224 147 108 71

が増⼤したため，𝑘 = 1の場合のみ左軸を⽤い，他は右軸を⽤いている．

この結果から，𝑘が⼩さい場合は動的順序で追加された計算，特に重み計算とハンガリアン法によ る割り当ての時間が⼤きな時間を占めてしまっていることがわかる．貪欲法によるマッチングはそれ らよりも⼩さいが，これは近似解法を⽤いているためであり，ハンガリアン法と同じオーダーとなる Edmond の Blossom 法によって厳密解を求めた場合はハンガリアン法と同程度の時間となることが予 想される．しかしこれらの計算時間は𝑘が増⼤するに従い劇的に減少する．これは予測された結果で あり，理論計算量に従うと，重み計算の計算量は𝑘に反⽐例し，貪欲法によるマッチングは𝑘 の反⽐

例よりも速く減少，ハンガリアン法は𝑘 に反⽐例するためである．

表 5.6は動的順序を⽤いたときの並列ステップ数である．ただし⽐較のため MMO の並列ステップ数 についても再掲している．MMO のものと⽐べると，約半分のステップ数となっており，動的順序の効 果がわかる．しかしながら，全体の実⾏時間⾃体はあまり変化しておらず，𝑛 = 4, 320のとき，MMO では1.76秒であったが，動的順序の場合，1.45秒，𝑛 = 8, 640のとき，MMO では5.75秒であったが，

動的順序の場合，4.84秒であった．これは，前処理・後処理の実⾏時間は両者で変化しないことや，重 み計算などの追加計算があること，そして，列ブロック交換における実⾏時間が増加したことによる．

表5.7 Modified Modulus 順序と Bečka の動的順序における 1 並列ステップあたりの列ブロック交 換の実⾏時間 (msec)

⾏列サイズ (𝑛) ⼿法 𝑘 = 1 𝑘 = 2 𝑘 = 4 𝑘 = 6 𝑘 = 8 𝑘 = 12 4, 320 MMO 0.357 0.445 0.353 0.331 0.584 0.362 Bečka 0.622 0.650 0.754 0.909 0.900 0.865 8, 640 MMO 2.39 2.60 1.80 2.40 1.99 2.03 Bečka 3.08 3.27 2.68 3.66 3.48 3.40

表 5.7は MMO と Bečka の動的順序における，1 並列ステップ当たりの列ブロック交換にかかった時 間である．実⾏時間が短いため，単位はミリ秒としている．この表からわかる通り，動的順序によっ て列ブロック交換にかかる時間は増⼤し，最⼤で 3 倍弱の時間となっていることがわかる．

以上の通り，Bečka の動的順序について，収束性の理論的解析を⾏い，また，Level-3 BLAS を活⽤

した重み計算の実装や，⼆次元ブロック分割に適した計算⼿法と，列ブロックのノードへの割り当て

⼿法を開発し，実⾏時間の⽐較とその内訳を検証した．その結果，縦⽅向の並列数𝑘が⼩さい場合，動 的順序に必要な新たな計算の負荷によってむしろ遅くなるが，最適な𝑘を⽤いた場合における実⾏時 間では MMO を上回る性能となった．将来的にノード数が増⼤した場合には，重み計算の計算速度の 改善や，ハンガリアン法の並列化，近似解法の適⽤などの改良が必要になると推測される．

第 6 章

部品の性能評価と実装技術

この章では，将来の⾼性能計算機へ向けた実装を考えるため，⽚側ブロックヤコビ法に⽤いられる 部品の内，Gram ⾏列の計算を取り出して議論する．Gram ⾏列の計算は前章の通り，⽚側ブロックヤ コビ法におけるオーダーの意味で⽀配的な 2 つの計算の内 1 つであるため，⾏列サイズ𝑛やノード数 𝑝がより⼤きくなった場合には重要となるものと考えられる．Gram ⾏列の計算は BLAS で定義され た⾏列積のパターンの 1 つである xSYRK であり，通常の⾏列積 xGEMM と同等の性能が期待される が，xSYRK の持つ対称構造のために性能効率が低く，前章の結果の中では，2 倍の演算量である列ブ ロックの更新と同程度の実⾏時間となってしまっている．これは，将来の計算機に頻繁に利⽤される と考えられるメニーコア CPU においてはより顕著となる．

そこでここでは xSYRK について，メニーコア CPU の 1 つである Xeon Phi (Knights Corner, KNC) への実装⼿法と，性能向上⼿法について議論する．KNC は，単体では遅い速度のコアを多数集積する ことで，全体としては⾼い性能を実現するものであり，1 つのコア当たりの速度は倍精度において約 17GFlop/s であるが，これを 60 個程度，集積することで，1 プロセッサでは約 1TFlop/s の⾼い理論ピ ーク性能を持つ．このように多数のコアを持つため，KNC において⾼い性能を実現するためには，⾼

い並列性を持つプログラムが必要であり，xSYRK の実装においてもこのことを考慮しなければなら ない．

KNC に対する⾏列積 (DGEMM) の実装⼿法については Smith ら [127, 128] が研究しており，BLIS と いう名前を持つ KNC などのメニーコア向け BLAS 実装も公開しているが，xSYRK については⼗分な 研究がされていない．そこで，Smith らの研究や，他にマルチコアプロセッサ向けの Level-3 BLAS 実 装についての後藤ら [129] の研究をベースに，KNC への DSYRK の並列化⼿法の実装と，性能解析を

⾏う．またここでは倍精度向けの DSYRK について議論するが，並列化⼿法やキャッシュ利⽤⼿法な どは他の精度や複素数においても役⽴つものと考えられる．

この章の内容は関連論⽂ [d] に基づくものであり，DSYRK において 3 つの並列化軸すべてを⽤いた 並列化⼿法を⽰し，実際に KNC 上で DSYRK の⾼速化を実現している点が新規である．また，将来的 に DSYRK を⽚側ブロックヤコビ法に⽤いる際に役⽴つものだと考えられる．