消去順序

Rijk の順序と列ソーティング

Rijk の消去順序 [110, Section 2.4] は部分的に動的なペアの選択を⾏うことで収束を⾼めることを⽬

的としている．Rijk の消去順序は⾏巡回順序を基にしている．⾏巡回順序は次のように考えることが できる．始めに，⾏を選択し，その⾏の要素を選択する．例えば，第⼀に1⾏⽬を選択し，その中の要 素(1, 2) (1, 3) (1, 4) ⋯を選択していく．次に2⾏⽬を選択し，その中の要素(2, 3) (2, 4) (2, 5) ⋯を 選択していく．ヤコビ法における “消去する要素” は⽚側ヤコビ法においては “直交化する列ベクトル のペア” となっている．そこで，⾏巡回順序は，1 つ列ベクトル，ピボット列を選択し，もう⼀⽅の列 ベクトルを次々の選択していく順序だと⾔える．Rijk の消去順序では，ピボット列の選択時にノルム 最⼤の列の置き換えを⾏う．つまり，⾏巡回順序で 1 列⽬を選択するときに，1, … , 𝑛列⽬の中から最 もノルムが⼤きいものを選び，1 列⽬とノルム最⼤の列を置き換えて，残りは同じように進⾏する．2 列⽬を選択するときは同様に2, 3, … , 𝑛列⽬の中からノルム最⼤のものと置き換えて，残りは同じよう

に進⾏する．これは選択ソートを消去順序に組み込んだものだと⾔える．選択ソートははじめに，1か ら𝑛要素⽬までの要素の中らか最⼤のものを探索し，第1要素と置き換え，次は2から𝑛要素⽬まで の要素の中から最⼤のものを探索し，第2要素と置き換える，といった動作を繰り返す．そこで Rijk の消去順序は，Givens 回転の⾓度𝜃が常に0であれば，列ベクトルをノルム⾮昇順に並び替えること ができる．𝜃が0でなければ列のノルムが変化しながら進⾏するため，完全にソートすることは不可 能だが，反復が進⾏し，列のノルムの変化が無視できる程度になれば，選択ソートと同等のことが発

⽣する．

Rijk は彼の消去順序の利点を，古典的ヤコビ法との関係で解説している．つまり，ヤコビ法にお ける⾮対⾓要素は𝑐_, = 𝑎 𝑎 = 𝛽_, ‖𝑎 ‖ 𝑎 と書ける．ただし𝛽_, は𝑎 と𝑎 の余弦⾓度だとす る．そこで𝑎 がノルム最⼤の列ベクトルだとすると，⾮対⾓要素の値も⼤きい傾向にあり，𝑎 との ペアから順に消去していくことは，絶対値の⼤きな⾮対⾓要素を優先的に消去することにつながる．

Mascarenhas [111] は 2.41 次収束のための消去順序で有名であるが，ソーティング⾃体が収束に与える 影響についても調べており，事前に⾏列のソートを⾏うだけでも収束性を改善する結果を⽰している．

Rijk の消去順序は動的な選択が⾏われるためそのままの形では並列化が難しい．そこで，共有メモ リ並列化においては，巡回毎に消去を始める前に列ベクトルをノルム⾮昇順に並び替える⽅法を採⽤

する．これは Rijk の消去順序において消去の中にソーティングが組み込まれているという巧妙さと⽐

べると，ソーティングを消去では別に⽤意しなければならず，冗⻑となるが，並列化には適している．

ヤコビ法の収束を考えれば，巡回が進⾏すると列のノルムの変化が⼩さくなり，列の⼊れ替えは少な くて済むだろうと考えられる．そこで，ほとんど整列済みの列に対して効率的なソーティングアルゴ リズムを⽤いることが適切である．本研究では，Python や Java のソートアルゴリズムとして使われ ている TimSort を⽤いている．TimSort については Auger らの論⽂ [112] が詳しいが，ほぼ整列した 数列に対して改良された Merge Sort の⼀種である．また，ノルムのソーティングと列ベクトルの並び 替えを別に⾏うことで，列ベクトルの移動を最⼩化することが適切だと考えられる．

並列化とキャッシュ利⽤

xGESVJ においては Rijk の順序に Block Oriented 順序（ブロック要素単位の⾏巡回順序とブロック 内部の⾏巡回順序の組み合わせ）を⽤いることでキャッシュ再利⽤性を⾼めている [3, Algorithm 1]．

本研究の実装では Rijk の動的選択の効果をソーティングによって実現できるものとし，並列化のため BOMO を⽤い，キャッシュ再利⽤性のため，ブロック内部ではブロック化列巡回順序を⽤いる．そこ で 2 種類のブロック幅が必要となる．1 つは，BOMO ⽤のブロック幅であり，これは BOMO の並列数 がスレッド数と⼀致するように設定する．また，もう 1 つのブロック幅はブロック内部のブロック化 列巡回順序のブロック幅であり，⾏列サイズとキャッシュサイズに依存して決定する．

三⾓構造の利⽤

xGESVJ は⼊⼒⾏列が三⾓⾏列である場合，ブロック三⾓要素単位の⾮零構造を利⽤して演算量を 削減する．また，⾮零構造のため演算量が少ない部分⾏列に対して複数回，消去を⾏うことで，演算 量を抑えつつ収束をより進める．ここでは Hari の V2 ⼿法を⽤いたときに出現する上三⾓⾏列の場合 における xGESVJ の動きを説明する．

xGESVJ ではまず⼊⼒⾏列𝐴を4つの列ブロックに分割する．𝐴は上三⾓であるため，

𝐴 =

⎡

⎢

⎢

⎣

𝐴 𝐴

𝐴 𝐴

𝟘 𝟘

𝟘

⎤

⎥

⎥

⎦

(5.1)

という形を持つ．すべての列ブロックが同じ幅を持つならば，𝐴 は𝐴の 4 分の 1 の⾏数となるため，

𝐴 に対する⽚側ヤコビ法は他の部分⾏列よりも⼩さなコストで⾏うことができる．そこで，𝐴 に対 しては 2 巡回の⽚側ヤコビ法を適⽤し，次に𝐴 ，そして𝐴 と𝐴 のペア，最後に𝐴 に対して⽚側 ヤコビ法を適⽤する．これを消去順序として書くと(1, 1)(1, 1)(2, 2)(1, 2)(3, 3)となる．これは Block Oriented 順序のブロック要素単位の順序であるため，対⾓ブロックを選択するものも含まれる．また，

すべてのブロック要素を選択するものではないため巡回順序でも擬似巡回順序でもない．これはあく まですぐに破壊される三⾓構造を利⽤して⼩さな演算量で消去を⾏うことを⽬的としているためであ り，収束するまでの反復は後の巡回に任せるためである．

三⾓構造の利⽤における並列化は，それぞれの𝐴 , 𝐴 , 𝐴 に対する⽚側ヤコビ法を BOMO を⽤いて

⾏うことで実現できる．𝐴 , 𝐴 のペアに対する並列化では⾮対⾓ブロックの消去であるため，⻑⽅領 域が出現するが，Modulus 順序から𝐴 のみや𝐴 のみに対するペアを除いた順序を⽤いて並列化する．

また，全体の並列順序と同じように，再帰的にブロック化列巡回順序と組み合わせることでキャッシ ュ再利⽤性を⾼める．

プログラムの概要

以上の機能を組み込んだ擬似プログラムが図 5.1である．この擬似プログラムでは簡単のため，擬 似プログラム中に現れる割り算がすべて割り切れるものとしている．また，冗⻑となるため，三⾓構 造を利⽤する部分については省いているが，この擬似プログラムを再帰的に呼び出すことで実現可能 である．擬似プログラム中のnb1は内側のブロック化列巡回順序のためのブロック幅であり，定数 CASHSIZEによって決める．最外側のdoループは巡回毎のループであり，収束判定と𝐴の列ベクトル のノルム⾮昇順によるソートを⾏う．ただし列ベクトルのソートでは𝐴と𝑉を同じように並び替えな ければ𝐴と𝑉が対応しなくなるため，sort-columnsは𝐴と𝑉の両⽅を引数とし，𝐴の並び替え順序 によって𝑉を並び替える．CASHSIZEは CPU のコアごとのキャッシュサイズを基に設定する値であ るが，実装においては⼿作業で探索した値を⽤いており，実装環境において⾏列サイズ𝑚 = 𝑛 = 3, 000 付近で最も⾼速となったCASHSIZE= 24, 000に設定している．この擬似プログラムは OpenMP によ って並列化することを前提とし，omp-parallel-doによって，CPU のコア数と同じ数だけスレッド を⽴ち上げる．そしてomp_num_threadsとomp_thread_numによってスレッド数とスレッド番号 を取得する．その内側のrに関するループは並列ステップ数であり，modulus-pairによって対応す るペアを⽣成する．またmodulus-pairが “あまりのペア” を⽣成した場合はflagがtrueを返し，

そうでない場合はfalseを返すため，flagの値によって，対⾓ブロックを選ぶか⾮対⾓ブロックを 選ぶかを選択する．if⽂の内側は 3 重ループとなっており，プログラムの最初に決定したnb1の幅 のブロック化列巡回順序を⾏うようになっている．

1 subroutine MTOSJ( [ ⋯ ], tol):

2 nb1 = CASHSIZE / n

3 [ ⋯ ] ,

4 do

5 maxt = 0

6 sort-columns( , )

7 omp-parellel-do

8 id = omp_thread_num()

9 nt = omp_num_threads()

10 nb = n/nt

11 for r=1:2*nt

12 (pp,qq,flag) = modulus-pair(2*nt, id, r)

13 if flag == false:

14 for p2=nb*(pp-1)+1:nb*pp:nb1

15 for =nb*(qq-1):nb*qq:

16 for =p2:p2+nb1-1

17 t = orthcols( , , , )

18 maxt = max(maxt, t)

19 end

20 end

21 end

22 else if flag == true:

23 for p2=nb*(pp-1)+1:nb*pp:nb1

24 for =p2+1:nb*pp:

25 for =p2:q

26 t = orthcols( , , , )

27 maxt = max(maxt, t)

28 end

29 end

30 end

31 for p2=nb*(qq-1)+1:nb*qq:nb1

32 for =p2+1:nb*qq:

33 for =p2:q

34 t = orthcols( , , , )

35 maxt = max(maxt, t)

36 end

37 end

38 end

39 end

40 end

41 while maxt > tol

42 for = 1,

43 ‖ ‖

44 end

45 ( , , … , )

46

47 return ( , , , )

図5.1 共有メモリ並列⽤の⽚側ヤコビ法の擬似プログラム