ｏｒｄｅｒｉｎｇ生成方法、プログラム及び共有メモリ型スカラ並列計算機
【課題】スパースな正値対称行列のコレスキーあるいは修正コレスキー分解を行う場合に、並列処理を高速化することができるorderingを生成する方法を提供する。
【解決手段】共有メモリ型並列計算機を用いてスパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行うにあたり、スパース行列が表す連立１次方程式が提示する問題における離散化された空間を、再帰的に２つの分割領域と、その間にある分割面とに分割する。分割を、分割面を構成するノードの数がスーパーノードの幅程度となったら止める。そして、再帰的に２分割されるごとに、分割領域内のノードに、分割面から遠いほうから順に番号付けを行う。分割面内のノードは、再帰的２分割のたびに、分割領域の番号付けの後に番号付けを行う。
【解決手段】共有メモリ型並列計算機を用いてスパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行うにあたり、スパース行列が表す連立１次方程式が提示する問題における離散化された空間を、再帰的に２つの分割領域と、その間にある分割面とに分割する。分割を、分割面を構成するノードの数がスーパーノードの幅程度となったら止める。そして、再帰的に２分割されるごとに、分割領域内のノードに、分割面から遠いほうから順に番号付けを行う。分割面内のノードは、再帰的２分割のたびに、分割領域の番号付けの後に番号付けを行う。
【発明の詳細な説明】
【技術分野】
【０００１】
スパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行なう場合における、orderingを生成する方法に関する。
【背景技術】
【０００２】
スパースな正値対称行列によって表される連立１次方程式（Ax=b、Aはスパースな正値対称行列、xは変数列ベクトル、bは定数列ベクトル）の解法において、変数列ベクトル内の変数の配置を並び替えることで、対応する行列のnonzero要素の位置も変わる。この並び替えをorderingと呼ぶ。適当なorderingを行うことでLDL＾T 分解のfill-in（元の行列ではzeroであったが、分解の結果nonzeroとなる行列要素） を削減することが知られている。LDL^T分解においては、スパース行列AをA=LDL^Tと分解して、連立１次方程式Ax=bを解く。ここで、Lは、下三角行列、^Tは、行列の転置を取ることを示す。
【０００３】
スパースな正値対称行列の対角要素を含む下三角行列部分L の列に関する、LDL^T分解での、行列要素間のデータ更新の依存関係はelimination treeと呼ばれる木構造で表される。elimination treeについては、非特許文献１を参照されたい。この依存性からleaf(elimination treeを構成する構成要素ノードのうち、子ノードを持たないノード）からroot（親ノードを持たないノード）にたどる方向に分解を行う。
【０００４】
elimination treeの互いに共通部分のない部分木は独立に分解の計算ができ、並列に処理することができる。
代表的なorderingにnested dissection がある。nested dissection とは、空間を２つの分割領域とこれらの分割領域を分割する分割面に分割して番号付けする方法である。nested dissectionについては、非特許文献１を参照されたい。さらに、スパース行列をグラフで表現して、グラフを構成ノード数がほぼ等しい２つの部分グラフに分割し、分割するために取り除くノード数が最小となるように分割をする。この2分割を、再帰的に繰り返し、もうこれ以上分割できないところまで分割して得られるorderingは一般化されたnested dissection によるorderingと呼ばれる。
【０００５】
分割は、分割面を形成するノードと分割された分割領域のノード集合に分かれる。あるレベルの分割に関して、分割面で分割された２つの集合をＡ，Ｂとしたとき、おのおのの領域をさらに分割しないようにＡに属する構成要素ノードに連続な番号を振ったあと、Ｂに属する構成要素ノードに連続な番号を振って、最後に分割面を構成する構成要素ノードに番号を連続に振る。この番号付けのルールは、再帰的に分割した部分にも適用する。
【０００６】
この結果生成された、nested dissection のorderingは、fill-in が少ないorderingであることが知られている。
上記では、LDL^T分解について述べたが、この分解は、修正コレスキー分解と呼ばれ、コレスキー分解（LL^T分解）とほぼ同じアルゴリズムが適用できる。したがって、以下の説明では、コレスキー分解について言及するが、同様のアルゴリズムが修正コレスキー分解にも適用できる。
【０００７】
スパースな正値対称行列の連立1次方程式をコレスキー分解で解くとき、行列の行あるいは列の並べ替えを行うことで、fill-in を少なくすることが一般的に考えられる。この並び替えをorderingと呼ぶ。並べ替えた行列をコレスキー分解する。
【０００８】
このコレスキー分解は、symbolic decomposition とnumeric decomposition の２つのフェーズから成る。symbolic decompositionで、各列の依存関係を解析し、elimination treeに表現する。この情報を利用して、スパース行列に存在するnonzero 要素の発生パターンを解析する。
【０００９】
コレスキー分解におけるnonzero要素の発生パターンは、elimination treeを生成して解析する。Ax=b（A=(aｉｊ））という連立方程式をコレスキー分解で解く場合、下三角行列をLとして、A=LL^Tと分解する。elimination treeは、以下の定理を用いて、最初の構成ノードから、子の構成ノードを見つけて連結することによって生成する。
・定理
コレスキー分解LL^T=Aにおいて、数値的なキャンセルを除くと、aki≠０、かつ、ｋ＞ｉであれば、ノードｉは、elimination treeにおけるノードｋの子ノードである。
（非特許文献１参照）
【００１０】
なお、ノードｋをノードｉの親ノードと呼び、ノードｉに対応するデータは、行列Aのｉ番目の列のデータである。また、elimination treeの構成要素ノードは、対応する番号の格子点に対応する。したがって、構成要素ノードｉは、離散化された空間の格子点ｉに対応し、構成要素ノードｉのデータは、行列Aの列ｉのデータである。
【００１１】
大規模なスパース正値対称行列のコレスキー分解を並列処理して性能を引き出す上で、elimination treeのお互いに交わりのない部分木は独立に計算できるので、独立な部分木については並列に計算を行なうことが考えられる。これ以外にも、elimination treeで示されるスパース行列の要素の更新の影響の伝搬関係に従い列を更新する処理を、順次パイプライン処理することで並列性能を引き出すことも考えられる。また、列を束ねてブロックを構成して計算量を大きくすることが並列効果を高める上でも重要である。一回のデータの読み込みで多くのデータを読み込んでおくと、速度の遅い、データの読み込み／書き出し処理の回数を減らすことができるので、並列処理効果が上がるのである。そのため、nonzero 要素発生パターンが同じ列をまとめてスーパーノードを構成する。スーパーノードを構成した場合には、スーパーノードを構成する列を結合したとき、nonzero な値を持つ行のデータのみを保持するようにする。したがって、zeroとなる行のデータを保持しなくて良くなるので、データを圧縮して格納して処理することができる。
【００１２】
さらに、コレスキー分解の性能を向上するために、条件を弱めて、分解の結果nonzeroパターンの似ている列も結合して、スーパーノードを構成する。この結合の結果、コレスキー分解の更新過程で、スーパーノードのデータ内に余分なzeroを含む行が発生し、演算量が増加する。しかし、単位演算量あたりのメモリアクセスを減らすことができる。このため、コレスキー分解の性能を向上することが期待できる。コレスキー分解の並列処理は、理論的には１列ごとの計算でも並列に行えるが、メモリアクセスの量や演算数が小さいため、並列処理のオーバーヘッドが無視できず、並列処理効果が引き出しにくい。
【００１３】
つまりメモリアクセスに対する演算数を上げること、および、並列処理のオーバーヘッドを削減するために粒度（granularity）を大きくすることが効率的な並列化を行うためには有効である。
【００１４】
理論的には、完全なnested dissection （空間の分割をこれ以上分割できないまでやって番号付けする方法）によるorderingを利用して行列を並び替え、コレスキー分解を、列を結合せずに、１列ごとに計算を行うとき、fill-in が最小になる。性能引き出しのために、１回のキャッシュへの読み込みで読み込まれるデータとして、列を結合してスーパーノードを作るが、elimination treeのleafに当る列から列の結合を行う処理を進める。nested dissection のorderingを利用してコレスキー分解で列の結合を行うとき、列を結合したブロックの幅を増やすためには、elimination treeの多くの分枝を有するノードを結合する必要がある。
【００１５】
分枝先ノードは分枝元ノードのnonzero要素とは共通のnonzero要素を多く持つが、分枝先のノード間では共通しないnonzero要素を多くもつ場合がある。このため多くのノードを結合すると、分枝ノードをかなり含むとき、zeros の数が急に増えることになる。そのため、ノードを順次結合して、zeroを少しずつ増やしながら、スーパーノードの適当なブロック幅を見つけることが難しいという問題点があった。
【００１６】
fill-in を削減する効果を期待してnested dissection orderingを使う場合、列の結合を行わないと、nested dissection によって発生したelimination treeの２分木構造のleafに相当する部分に、列一本に対して多くの計算が発生する。fill-in は少なくなるが、メモリアクセスの演算に対する割合や、並列処理の粒度が小さいことによるオーバーヘッドが顕在化してくるため、期待する性能が得られないことになる。
スパース行列の直接解法に使われる基本技術の詳細は、非特許文献１を参照されたい。
【００１７】
＜nested dissection から生じるorderingとそのelimination tree＞
２次元の正方形領域での問題を考える。この領域での問題を表現する偏微分方程式を離散化してスパースな正値対称行列を得る場合を考える。このとき、離散化された格子点は、領域の中で、隣接する格子点と結合されている。
【００１８】
図１のような、x-y 平面で考える。この領域を構成する格子点からｙ軸と平行な直線を構成する格子点の集合ａを取り除くことで、領域を構成する格子点の集合を２分割して、おのおのが構成する格子点の数がほぼ等しくなるようにする。
【００１９】
残った２つの領域のおのおのを、ｘ軸と平行な横方向の格子点の集合ｂを取り除くことで、同様に２分割して、おのおのが構成する格子点の数がほぼ等しくなるようにする。格子点の数がほぼ等しくなるように分割するには、今分割しようとする領域の全格子点数を２で割った数の格子点からなるように領域を分割すればよい。この手順を、分割によって生成された小さくなった領域が、これ以上分割できないようになるまで再帰的に分割を続ける。これ以上分割できないとは、分割して得られた小領域に含まれる格子点の数が３以下となって、次の４分割ができないような場合を示す。最終的に領域は2**nの小さな領域に分割される。
【００２０】
この処理は、領域をほぼ均等に分割して、取り除く境界線を構成する格子点数が少なくなるように選ぶことに相当している。
番号付けは、分割で生成した２つの領域の内の一つの領域を構成する格子点の集合に対して連続に番号を割り当てる。次に、もう一つの領域を構成する格子点の集合に対して引き続いて番号を連続に割り振る。最後に、２分割で取り除いた格子点に引き続き番号を割り振る。
【００２１】
このルールを再帰的に適用して、上記分割から生じた集合に番号付けを行う。この番号づけで並べた行列のelimination treeは２分木になる。
【００２２】
すなわち、例えば、離散化された空間上でラプラスの方程式（Δφ＝）を解くことを考える。一次元の問題であるとすると、ラプラシアンは、座標の２階偏微分になるので、φ（ｉ）を格子点ｉにおける関数値とすると、離散化されたラプラスの方程式は、
｛φ（ｉ＋１）＋φ（ｉ−１）−２φ（ｉ）｝／ｈ２＝０
となる。ここで、ｈは、隣り合う格子点間の距離である。これを行列で表すと、
【数１】
のようになる。これから明らかなように、隣り合った格子点が連続的に番号付けされている場合には、スパース行列内では、係数は隣接して発生する。しかし、格子点を上記のように分割すると、分割された領域内は、隣り合った格子点に連続した番号が現れるので、係数は固まって発生するが、分割して切り離され、別個に番号付けされた格子点に対応する係数は、他の係数とは離れた場所に発生する。したがって、このような行列のelimination treeを生成すると、係数が固まっている列同士に対応するツリーは互いに結合されるが、離れてしまった列同士は別のツリーとなる。したがって、領域を、ＡとＢ及びこれらの境界領域Ｃとに分割した場合、Ａ内の格子点に対応するelimination treeの構成要素ノードと、Ｂ内の格子点に対応するelimination treeの構成要素ノードとは、境界領域Ｃの格子点に対応するelimination treeの構成要素ノード以降において、２つに分枝することになる。
【００２３】
この２分木は、分枝するノードとノードが直線に並んだ枝から成る。
以下に、再帰呼び出しの手続きでの番号付けの方法を記述する。
以下のプログラムコードにおいて、関数partition は領域の分割を行う手続きで、領域を分割できるかを判定して、できるときは、seta,setb に２分割した領域を、cutsurfaceに分割面を返す。setは、今分割したい格子点の集合である。
call numbering(set)
recursive procedure numbering(set)
call partition(set,seta,setb,cutsurface,icon)
if(partition is done)then ! numbering を呼ばれた順に番号付けは連番で行う。
call numbering(seta) ! setaを番号付けする。
call numbering(setb) ! setbを引き続き番号付けする。
call numbering(cutsurface) ! cutsurfaceを番号付けする。
else
call numbering(set)
endif
return
end
なお、partitionのアルゴリズムについては、非特許文献２を参照されたい。
【００２４】
３次元問題に関しても同じように考えることができる。例えば、立方体ならy-z 平面に平行な面を構成する格子点の集合を取り除いて、領域をほぼ均等になるように分割する。引き続き、分割で生成された２つの領域を同様にx-z 平面に平行な面でおのおのを２分割する。これらを２次元の場合と同様に、これ以上分割できない状態まで繰り返す。
【００２５】
番号付けのルールも同じである。この番号付けで並び替えた行列のelimination treeはやはり２分木になる。
これらを、一般の連結する領域に対して適用するために、離散化した格子点の結合関係をグラフで表現し、グラフを分割する面を構成するノード数が最小で、分割した部分グラフを構成する構成要素ノードの数が均等になるように分割することを繰り返して、一般化されたnested dissection orderingを構成することが考えられている。この場合も同様に、elimination treeは２分木になる。
【００２６】
＜ノードの結合（スーパーノード）＞
elimination treeの構成要素ノードに対応する行列の列で、引き続くノード番号の構成要素ノードに関して、ノード番号の大きい方のnonzeroパターンがノード番号の小さい列と等しい列を結合する。また、枝を構成する構成要素ノードを調べてnonzero パターンが比較的に似ていて、結合してもzeroの割合が大きく増えない場合は結合する。この結果、コレスキー分解すると、結合した列ベクトルの行で、nonzero 要素を持つ行の数は増える。つまり、この行を圧縮して格納するのに必要な領域にzeroが含まれることになる。こうしてnonzeroパターンが似ていて列を結合してできたブロックをpanel と呼ぶ。
【００２７】
図２は、panelとindex listを示す図である。
panelには、index listが対応して設けられる。index listは、panelのデータがスパース行列の列の中の何行目に当るかを保持する。panelは、panelに保持される列に共通にzeroである要素以外のnonzero要素のデータしか保持しないので、index listは保持されているnonzero要素が列の中の何行目かを示すのである。なお、panelは、保持される列に共通にzeroとなる要素を保持しないものであるので、nonzero要素として保持されていても、実際の数値はzeroとなることがある。
【００２８】
panel ごとのコレスキー分解の処理は、ブロック幅が大きくなるほど演算量が増えるため、相対的にメモリのロードストアのコストが下がる。ただし、zero要素の部分は計算する必要のない無駄な計算が発生する。このため、結合によるブロック幅を、含まれるzeroの数を考慮しながら調整して、ブロック幅を決定する必要があり、このことが性能を引き出す上で重要である。
【００２９】
スーパーノードはelimination treeにおいて連続する構成ノードで、nonzero 要素のパターンが同じものを結合して、nonzero 要素がある行のみ集めて圧縮した形で格納する。すなわち、nonzero要素の行のみデータを保持し、zero要素の行のデータを保持しないことによって、保持すべきデータの量を圧縮するものである。
【００３０】
図３は、panelに不要なzeroを含む様子を示した図である。
elimination treeの親の構成要素ノードでnonzero 要素のパターンが一致するものを結合したものを、各列のnonzero 要素のみ圧縮して格納する。これはpanel となる。次の親ノードとはnonzeroパターンは一致しないが、包含関係にある。すなわち、次の親ノードにおいてnonzeroである要素に対応する行は、分解の処理における更新処理において更新され、nonzero要素となるので、次の親ノードのnonzeroパターンは、panel内に反映されることになる。このような関係を包含関係と呼んでいる。このような次の親ノードをpanelに結合すると、不要な0 をpanel に含むようになる。図３の例の場合、panelの下方に不要なzeroが生じており、この場合のインデクスは、panelと、結合された親ノードとで共通でない部分が大きい。
【００３１】
＜分枝するノードを結合すると大きくzeroが急激に増える。＞
分枝するノードから分かれる枝の先につながるノードは、nested dissection では分割された領域を構成するノードに対応する。
【００３２】
コレスキー分解の結果、ノードに隣接するノードは、ノードに対応する列のnonzero要素として現れる。コレスキー分解の過程で、elimination treeのrootへのパスをたどるとき、nonzero要素は、他の列への更新の際に影響を与える要素としてその影響が伝播することになる。つまり、ある面で分割された２つの領域の、この面以外で隣接する面のノードはnonzero要素として分枝ノードに伝わる。２分割されたノード群は分割面以外の境界面は共通要素を持たないので、これに隣接する領域外のノードに対応するnonzero要素は分枝ノードに受け継がれる。すなわち、分枝ノードにおいて、zeroである要素も、結合することによって、nonzero要素として扱う必要が生じる。つまり、分枝ノードを結合すると、分枝ノードの枝を構成するノードを結合したパネル幅に、zeroが急激に増加する。
【００３３】
nested dissectionでは、elimination treeの構成要素ノードをleafから結合していくと多くの分枝ノードに遭遇するためzeroの数とpanel の幅を調整することが難しい。
図４のような直方体がいくつかの分割の結果、生成されたとする。真ん中の部分は直方体を２分して取り除く面を構成するノードの集合を表す。分割前はこれら３つの部分はつながっている。また、分割前の表面は、他の領域と接していると仮定する。
【００３４】
分割領域内のノード（格子点）に番号付けするとき、領域を分割しないように適当に番号付けする。このときelimination treeでは、これらはY をさかさまにしたような木の部分として表現される。
【００３５】
分枝ノードで分かれる部分を考えると、elimination treeは、分枝ノードから下に２又に分かれる。コレスキー分解した結果の下三角行列L の構造は、親ノードのnonzero パターンの構造に包含される。そのため、ブロックA とC 以外で接する他の領域に存在するノードのnonzero は、出現した位置から伝播していく。
【００３６】
図５に示される列データを用いて説明すると、分枝ノードy に対し、分枝先の２つのノードz,w のnonzero パターンは図５の例のようになる。すなわち、ノードzは、対角要素のほかに、ブロックCの構成ノードから発生するnonzero要素ＣＣＮと、ブロックAのブロックC以外と接する外部ノードから発生するnonzero要素ＡＮからなる。ノードwは、対角要素のほかに、ブロックCの構成ノードから発生するnonzero要素ＣＣＮと、ブロックBのブロックC以外と接する外部ノードから発生するnonzero要素ＢＮとからなる。一方、分枝ノードyは、対角要素のほか、ブロックCの構成ノードから発生するnonzero要素ＣＣＮと、ブロックAのブロックC以外と接する外部ノードから発生するnonzero要素ＡＮと、ブロックBのブ、ロックC以外と接する外部ノードから発生するnonzero要素ＢＮと、ブロックCのブロックA,B以外と接する外部ノードから発生するnonzero要素ＣＮからなる。
【００３７】
分割は領域を、互いに均等な大きさになり、かつ、分割面が最小になるように分割しているとする。つまり、ＡＮはＢＮにほぼ等しく、外部に接する面が多ければ、ＡＮ〜ＢＮ＞ＣＣＮ＞ＣＮとなる。
【００３８】
この結果、図５の右側に示されるように、かなり多くのzeroが急激に増すことになる。
領域の分割をさらに進めてもう分割できないところまで続けたものを結合すると、この構造が繰り返される。
【先行技術文献】
【非特許文献】
【００３９】
【非特許文献１】T.DAVIS, "Direct Methods for Sparse Linear Systems", SIAM 2006
【非特許文献２】George Karypis and Vipin Kumar, "A fast and high quality multilevel scheme for partitioning irregular graphs", SIAM Journal on Scientific Computing, Vol. 20, No. 1, pp. 359 - 392, 1999
【発明の概要】
【発明が解決しようとする課題】
【００４０】
以上のように、スーパーノードの幅を広げて、計算の並列性を高め、高速に処理しようとしても、スーパーノードの幅が広いとスーパーノードに含まれるzeroの数が多くなり、無駄な計算が増えてしまう。これにより、スーパーノードの幅を広げたことによって並列性を高めた効果は少なくなってしまう。zeroは、elimination treeの分枝ノードを結合することによって多く発生する。したがって、elimination treeの形状が、分枝ノードが少ない形状であれば、このような問題は軽減される。ところで、elimination treeの形状は、行列のnonzero要素の配置によって決まり、nonzero要素の配置は、空間の格子点の番号付けの方法による。したがって、適切なorderingによって番号付けられた格子点を用いれば、zeroを急激に増やすことなく、スーパーノードの幅を広げることができる。
【００４１】
以下の実施形態は、スパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行なう場合に、並列処理を高速化することができるorderingを生成する方法を提供する。
【課題を解決するための手段】
【００４２】
以下に開示される実施形態の一側面におけるordering生成方法は、共有メモリ型スカラ並列計算機によってスパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法であって、共有メモリ型スカラ並列計算機は、スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行い、１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う。
【発明の効果】
【００４３】
以下の実施形態によれば、スパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行なう場合に、並列処理を高速化することができるorderingを生成する方法を提供することができる。
【図面の簡単な説明】
【００４４】
【図１】従来技術を説明する図（その１）である。
【図２】従来技術を説明する図（その２）である。
【図３】従来技術を説明する図（その３）である。
【図４】従来技術を説明する図（その４）である。
【図５】従来技術を説明する図（その５）である。
【図６】本実施形態で使用されるデータ構造を説明する図（その１）である。
【図７】本実施形態で使用されるデータ構造を説明する図（その２）である。
【図８】本実施形態で使用されるデータ構造を説明する図（その３）である。
【図９】本実施形態で使用されるデータ構造を説明する図（その４）である。
【図１０】本実施形態で使用されるデータ構造を説明する図（その５）である。
【図１１】本実施形態の全体処理のフローである。
【図１２】本実施形態の詳細な処理フロー（その１）である。
【図１３】本実施形態の詳細な処理フロー（その２）である。
【図１４】本実施形態の詳細な処理フロー（その３）である。
【図１５】本実施形態の詳細な処理フロー（その４）である。
【図１６】本実施形態の詳細な処理フロー（その５）である。
【図１７】本実施形態の詳細な処理フロー（その６）である。
【図１８】本実施形態の詳細な処理フロー（その７）である。
【図１９】本実施形態の詳細な処理フロー（その８）である。
【図２０】本実施形態の詳細な処理フロー（その９）である。
【図２１】本実施形態が適用される共有メモリ型スカラ並列計算機のハードウェア構成を説明する図である。
【発明を実施するための形態】
【００４５】
本実施形態においては、大きな粒度の並列性を引き出し、列を結合してスーパーノードの効果的なブロック幅を確保するために、ある深さでグラフの分割を止めて最後に分割面で２分割された２つのノードの集合内で番号を振る。このある深さでグラフの分割をやめるとは、取り除く分割面を構成するノード数が、ノードの結合を行うことを考えたときのスーパーノードのブロック幅の大きさ程度に縮小したかを判断して、その程度に達した場合に、その時点の深さで繰り返し行う再帰的分割を止めるように制御することを意味する。
【００４６】
そして、スーパーノードを構成する場合、elimination treeの枝を形成する連続する列を結合していく。このとき、あるノードを結合したときzeros の増加が急激に大きくならないような並びで番号付けすることを考える。すなわち、スーパーノードを列を結合して形成する場合、分枝ノードを結合するとzeroが大幅に増える。したがって、分枝ノードができるだけelimination treeに現れないようにするのが好ましい。ところで、elimination treeの形状は、行列のnonzero要素が行列内のどの位置に現れるかによって変化する。そして、行列のnonzero要素が行列のどの部分に現れるかは、行列要素の番号付けによって決まり、行列要素の番号付けとは、すなわち、空間を離散化した場合の格子点の番号付けの仕方によって決まる。格子点は、elimination treeの構成要素ノードに対応する。したがって、格子点（構成要素ノード）の番号付けを工夫することにより、zeroを大幅に増やすことなく、スーパーノードの幅を太くすることができる。スーパーノードの幅は、構成要素ノードの接続の仕方によって可変でよいが、上限値を大きく取ることができるようになる。スーパーノードの幅が太いほど並列性の効果を得やすいが、この幅は、当業者が使用する計算機の性能を考慮して、適宜設定すべきものである。例えば、スーパーノードの幅の上限値は、５０〜１００の列を結合した幅とすることが可能である。
【００４７】
このためには、分割を止めて残った領域の番号付けが、さらに領域を分断しないようにすることと、ノードの結合を行うとき、徐々にzeroが増加するような並びにすることとがpanel の幅とzeroの割合を制御し、性能を引き出す上で効果的なスーパーノードの構成を可能にする。空間を分割した結果のブロック内の番号付けでは境界のノードをどう順序付けするかという問題が存在する。
【００４８】
メモリアクセスのスピードに応じて２つの方法が考えられる。
１）ある程度の深さまででグラフの再帰的分割をやめる（取り除く分割面を構成するノード数が、ノードの結合を行うことを考えたときのスーパーノードのブロック幅の大きさ程度となったら分割をやめる）。最後の深さレベルを持つ２分割で生成された領域群のおのおのの領域を構成する構成要素ノードに番号を振るとき、２つに分割する面（分割面）の構成要素ノードからの距離が遠いものから順番に番号を振る方法。
２）もう一つは、生成された分割領域に接する分割面をすべて考えたとき、これらのnested dissection を生成するときに取り除く全ての分割面との距離が遠くなるような構成要素ノードから順番に番号を振る方法。
【００４９】
方式１および２とも、再帰的分割の最も深いレベルでの分割面は、帯型のような形状でnonzero 要素が並ぶ。
方式２の方は、多くの列を結合しなくても性能の引き出しが可能なメモリアクセスが高速である場合に適している。
【００５０】
方式１あるいは２を使用すると、ノードを結合してもzeroの増加はなだらかである。また、zeroが多く発生するのを抑えるためには、分枝ノードを結合することはやめて、elimination treeの直線を形成する枝の部分に着目して、スーパーノードの結合を行えばよい。
【００５１】
図６〜図１０は、本実施形態で使用されるデータ構造を説明する図である。
図６〜図１０は、ノードの隣接ノードを表すデータ構造を表す図である。このデータは、分割面を構成する構成要素ノードの集合、つまり最後に分割されたグラフの構成要素ノードの集合、分割面を表す集合から２つに枝分かれする構造を再帰的に表した木構造（離散化空間を分割した場合の各領域を示す領域ノードをtree構造で示したもの）、各集合に番号を振って識別するとき使う集合の構成ノード数からなる。
【００５２】
なお、以下では、elimination treeを構成するノードを構成要素ノード、領域の分岐の様子を示す木構造を構成するノードを領域ノードという。
図６は、各構成要素ノードに隣接する構成要素ノードを表す。
【００５３】
nadj(*)は、リストベクトルに各構成要素ノードの隣接ノードを格納する。nstart(n+1)は、ｉ番目の構成要素ノードの隣接ノードが１次元配列（nadj(*))の何番目から格納されるかを示す。
【００５４】
図７は、分割された領域を識別するデータを表す。
nodeinfは、分割面か分割された領域を識別する番号を、構成要素ノードのノード番号に対応させて格納する。構成要素ノードごとに、構成ノードポインター(チェイン構造で構成要素を表現する) 、領域識別番号、新番号、処理済flagを持つ。これらの要素をすべてのノード数分の配列の形で持つ。
【００５５】
図８は、分割領域の構造を示すデータ構造である。
領域を分割するときに、領域から分割面を取り除き、共通部分のない２つの領域を生成する。elimination treeにおいて、この２つの領域から分割面の構成要素ノードを経由せずに、他方の領域の構成要素ノードに向かうedgeはない。このような、分割を再帰的に繰り返し、ある深さで止める。ある深さで分割を止めるとは、再帰的分割において、取り除く分割面を構成するノード数が、ノードの結合を行うことを考えたときのスーパーノードのブロック幅の大きさ程度となった深さで分割をやめることを示す。
【００５６】
深さをS としたとき、分割されて残った領域は2**S個、分割面の数は2**S-1である。分割されて残った領域および分割面をノード（領域ノード）として木構造に表現する。最初の分割面をrootとして、２つに分かれた領域をさらに分割する分割面を木の枝として結合する。順次このように結合し、最後の分割でできた分割面は、その子として、２つの分割領域を表すノードを結合する。
【００５７】
図８（ａ）において、各ノード（分割領域を示すノード：領域ノード）のparentは一つレベルの小さな分割の分割面のノードを指す。child はレベルが増した分割面を示すか、最も深いレベルの分割領域を示すノード番号を指す。brother は配列で、child のbrother （深さが同レベルの領域ノード）の番号を指す。
【００５８】
領域ノードの表す分割面または分割領域の領域ノードを構成要素ノード（elimination treeの構成ノード）のチェイン構造の先頭からポインタで表現する。領域ノードが示す領域の構成要素ノードの総数を構成要素ノード総数に格納する。また、分割面および分割領域を識別するための番号をtreeノード番号に格納する。
【００５９】
図８（ｂ）において、図８（ａ）の構造のデータを、配列として保持する。深さレベルをSとすると、分割された領域と分割面のデータを、2**(S+1)-1個の要素を持つ配列として表現する。また、root領域ノードを指すポインタを別個に保持する。
【００６０】
図９は、領域を深さレベル３まで分割した場合の木構造を表す。
図９の番号付けは、treeをdepth first searchの順番で番号を振った。図９において、root領域ノードの番号は１５で、レベル０である。レベル３の領域ノード１、２に接続されるparent領域ノードは、レベル２の３である。レベル３の領域ノード４、５のparent領域ノードは、６である。レベル２の領域ノード３、６のparentノードは、レベル１の領域ノード７である。同様に、領域ノード８、９のparentが領域ノード１０、領域ノード１１，１２のparentが領域ノード１３、領域ノード１０、１３のparentが領域ノード１４である。root領域ノード１５は、領域ノード７、１４のparentである。逆に、領域ノード１、２は、領域ノード３のchildであり、領域ノード４、５は、領域ノード６のchildであり、領域ノード３、６は、領域ノード７のchildであり、領域ノード７、１４は、root領域ノード１５のchildである。また、領域ノード１に対する領域ノード２は、brotherである。以下同様である。
【００６１】
図１０は、番号付けで使うstackのデータ構造を示す。
stacksf およびstcknodeを用意する。スタックには、１次元配列を利用する。方式１では全処理について、それぞれstacksfとstacknodeを一つずつ設ける。方式２では、stacksfは、全処理について１つであるが、stacknode は、分割領域ごとに確保して利用する。
【００６２】
また、番号付けに際し、permutation vectorを設ける。permutation vectorは、ｉ番目の構成要素ノードにｊが振られた場合には、ｊという値に対し、ｉを返すものとする。すなわち、perm(j)=iとなるようにする。
【００６３】
図１１は、本実施形態の全体処理のフローである。
スパース行列をグラフ（elimination tree)に表現して、それを再帰的に分割するが、ノード（列）を結合してスーパーノードを作るときのブロック幅の上限から、これ以上分割するか制御する。これは、分割面を構成するノードを結合して適当な大きさのスーパーノードを構成できるようにするためである。
【００６４】
ステップＳ１０において、スパース行列をグラフ（elimination tree)に表現し、隣接ノードを表現するデータを作成する。ステップＳ１１において、グラフを分割して分割領域Ａ，Ｂおよび分割面Ｃを生成する。そしてＡ，Ｂを表すグラフを作成する。すべての分割面Ｃの構成ノード数の最低値が、スーパーノードのブロック数とほぼ等しいかを判断し、等しいなら、分割をやめて、ステップＳ１２の処理を行う。等しくないなら、すべての分割領域に関してステップＳ１１にもどって分割する。ステップＳ１２において、分割領域、分割面をノードとする２分木を作成し、構成要素をチェインとして格納する。また、２分木の領域ノードに属する構成要素ノードに分割領域、分割面を識別する番号を付与する。ここでは、stacksf、stacknodeを使用する。番号付けの結果は、permutation vectorに格納する。
【００６５】
以下の処理は、方式１と方式２で処理手順が異なる。
方式１の処理手順：
生成された、領域ノードからなるtreeの最大深さレベルをS としたとき、2**(S+1)-1個の領域ノードがある。
root領域ノードからdepth first searchで、順番に処理する。処理しようとするノードがleafノード（子ノードを持たないノード）のとき、その親ノードは分割面を表す。
【００６６】
１）depth first searchで取り出したノードがleaf node であるとき。
leafノードの親ノードである分割面を構成する構成要素ノードのチェインをたどり、構成要素ノードをスタック(stacksf) に積む。
このstacksf から構成要素ノードを一つずつ取り出して、取り出した構成要素ノードの隣接ノードをスキャンする。隣接ノードが現在対象としている分割領域に属する要素ノードであり、かつ、まだstacknode に積まれていないノードであれば、処理済のflagを設定して、stacknode に積む。
【００６７】
この後、stacknode の一番下から、ノードを一つずつスキャンして、スキャンした構成要素ノードの隣接ノードを順番に調べて、対象の分割領域に属し、かつ、まだstacknode に積まれていないなら、処理済みのflagを設定して、stacknode に積む。スキャンする構成要素ノードがなくなるまで、この処理を続ける。
【００６８】
その後、stacknodeのノードを上から下へと取り出して、相対番号を１から振る。全体での番号は、いままでに番号を振った総数＋相対番号になる。これを構成要素ノードごとに設定する。
【００６９】
２）depth first searchで取り出したノードがleafノードでないとき。
分割面を構成する構成要素ノードを表すチェインをたどり、相対番号を１から振る。全体での番号は、いままでに番号を振った総数＋相対番号になる。これをノードごとに設定する。
【００７０】
３）最後に、もとのノード番号の順番でノードの情報が格納されている１次元配列を元のノード番号順にスキャンする。ｉ番目の要素に新たに振られた番号ｊが格納されているとする。permutation を表すベクトルpermにperm(j)=i と格納してpermutation vectorを作成する。
【００７１】
方式２の処理手順：
生成された領域ノードからなるtreeの最大の深さをS とすれば、leafノードはm=2**S個ある。そこで2**S個のstacknode を各leafノードに対して確保し利用する。
そして、treeをdepth first searchでleafへ下るpathをたどる途中で分割面を通過するときに、つぎの処理を行う。
【００７２】
分割面の構成要素ノードを構成要素ノードのチェインをたどって、stacksf に積む。これが終わったら、stacksfから一つずつ取り出して、その隣接ノードをスキャンする。スキャンしたノードがleafノード、つまり最も深いノードで、分割領域の構成要素ノードであるとき、そのleafノードに対応する構成要素ノードをstacknode に積む。そして、そのノードは処理済みflagをonにする。
【００７３】
１）depth first searchで取り出したノードがleafノードであるとき。
depth-first searchで、leafノードに到達したら、stacknode の底から順番に構成要素ノードを取り出す。その構成要素ノードの隣接ノードをスキャンして、このleafノードの構成要素ノードであって、未処理のものについて、処理済のflagをonに設定して、stacknode に積む。この処理を繰り返してstacknode から取り出すノードがなくなったら、stacknode の一番上から下へと順番にノードを取り出して、構成ノードに相対番号を１から振る。
【００７４】
このleafノードまでにそれらの構成要素ノードに番号付けした総ノード数に相対番号を加えて番号を割り当てる。
【００７５】
２）depth first searchで取り出したノードがleafノードでないとき。
分割面を表す構成要素ノードを表すチェインをたどり、相対番号を１から振る。全体での番号は、いままでに番号を振った総数＋相対番号になる。これをノードごとに設定する。
【００７６】
３）最後に、もとのノード番号の順番でノードの情報が格納されている１次元配列を元のノード番号順にスキャンする。ｉ番目の構成要素ノードに新たに振られた番号ｊが格納されているとする。permutation を表すベクトルpermにperm(j)=i と格納してpermutation vectorを作成する。
【００７７】
図１２〜図２０は、本実施形態の詳細な処理フローである。
図１２においては、スパース行列をグラフ（elimination tree）で表現して、再帰的に２分割する。分割のレベルごとにすべての分割面を構成するノード総数の最低値が、スーパーノードのブロック数程度になったか否かを判断して分割を止める。その後、２分木の作成、分割面、分割領域の構成要素ノードのチェインの作成、どの分割面、分割領域に属するかの識別番号の設定を行う。
【００７８】
ステップＳ２０において、行列をグラフに表現する。ここでグラフといっているのは、elimination treeのことである。ステップＳ２１において、領域を表すグラフ（格子点とその接続関係を示す空間のグラフ）を、領域A、領域B、分割面Cに分割する。すなわち、離散化された空間を各分割領域と分割面に分割する。
【００７９】
ステップＳ２２において、直近の分割で生成された分割面の構成要素ノードの総数がスーパーノードのブロック幅（スーパーノードに結合される列数）の上限値程度か否かを判断する。ステップＳ２２の判断がＮｏの場合には、ステップＳ２３において、直近の分割で生じた領域ごとに、グラフ（elimination tree）を作成し、ステップＳ２１に戻る。ステップＳ２２の判断がＹｅｓの場合には、ステップＳ２４において、分割で生成された分割面を親ノードに、２つの分割領域を子ノードに持つ２分木（領域のtree、図８のデータ構造で、領域ノードを構成ノードとするtree）を作成する。各分割領域、分割面の構成要素ノードをチェインで関連付け、分割領域、分割面と対応付ける。この対応付けが識別できるように構成要素ノードに識別番号を対応付ける。すなわち、elimination treeの各構成要素ノードがどの領域に属するかを示す領域識別番号を、構成要素ノードごとに、図７のデータ形式の領域識別番号領域に入力する。
【００８０】
図１３〜図１６は、方式１の処理フローチャートである。
図１３及び図１４では、elimination treeの構成要素ノードの集合であって、分割面、分割領域を表す領域ノードのtreeをrootからdepth-first searchの順にたどって番号付けする。このtreeは２分木で、leaf ノードの深さレベルは、同じ。brother はあれば必ず１つとなる。
【００８１】
ステップＳ２５において、curparentという変数に、rootノードを設定する。また、curparentのchildをcurchildという変数に設定する。そして、curparentのもう１つのchildをbrohterに設定する。この場合、１回の分割で生じた、もう一つの領域をbrotherとして設定する。
【００８２】
ステップＳ２６において、curchildがnull、すなわち、存在しないか否かを判断する。ステップＳ２６の判断がＮｏの場合、ステップＳ２７において、curparentにcurchildを設定し、curparentのchildをcurchildを設定する。そして、curparentのchildに、brotherを設定し、ステップＳ２６に戻る。ステップＳ２６の判断がＹｅｓの場合には、ステップＳ２８において、grandpという変数に、curparentのparentを設定する。ステップＳ２９において、cutsurfaceにgrandpを設定し、nodesetにparentを設定し、分割領域の構成要素ノードの番号付けを行なうサブルーチンを呼び出す。
【００８３】
ステップＳ３０において、nodesetにgrandpのchildを設定し、grandpのchildに、brotherを設定する。ステップＳ３１において、設定されたcutsurface、nodesetで、分割領域の構成要素ノードの番号付けのサブルーチン（後述する）を呼び出す。ステップＳ３２において、curparentに、curparentのparentを設定する。
【００８４】
ステップＳ３３において、curparentのchildがnullか、すなわち、存在しないか否かを判断する。ステップＳ３３の判断がＮｏの場合には、ステップＳ３８において、curparentのchildをcurchildに設定し、curparentのもう１つのchildに、brotherを設定して、ステップＳ２６に戻る。
【００８５】
ステップＳ３３の判断がＹｅｓの場合には、ステップＳ３４において、nodesetに、curparentを設定する。ステップＳ３５において、cutsurfaceの構成要素ノードを番号付けするサブルーチンを呼び出す。ステップＳ３６において、curparentにcurparentのparentを設定する。この動作は、領域ノードからなる領域の分割を表すtreeをrootの方向に戻ることを意味する。ステップＳ３７において、curparentがnullか、すなわち、存在しないか否かを判断する。ステップＳ３７の判断がＹｅｓの場合には、全ての領域ノードについて処理し終わったので、処理を終了する。ステップＳ３７の判断がＮｏの場合には、ステップＳ３３に戻る。
【００８６】
図１５は、図１３のステップＳ２９、Ｓ３１のサブルーチンの処理フローである。
図１５においては、分割面を構成する構成要素ノードの集合から遠い構成要素ノードから、領域ノードの構成するtreeのleafに対応する分割領域の構成要素ノードに番号を振る。
【００８７】
ステップＳ４０において、図１０で示されたstacksfとstacknodeをクリアする。また、スタック内の位置を示すポインタstackptrを０に、stackptr2を１に初期化する。ステップＳ４１において、分割面を構成する構成要素ノードのチェインをたどって、たどった先の領域に属する構成要素ノードを取り出し、取り出したノードをstacksfに積む。ステップＳ４２において、stacksfから１つの構成要素ノードを取り出す。ステップＳ４３において、構成要素ノードの全ての隣接構成要素ノードを順番にスキャンし、領域番号を調べて、現在処理している分割領域に属しているが、stacknodeにまだ積まれていない構成要素ノードをstacknodeに積む。このとき、stackptr＝stackptr＋１として、構成要素ノードを積む位置を計算する。そして、新たに、構成要素ノードをスタックに積むとき、構成要素ノードの処理済flagをonにする（図７参照）。
【００８８】
ステップＳ４４において、stacksfが空か否かを判断する。ステップＳ４４の判断がＮｏの場合には、ステップＳ４２に戻る。ステップＳ４４の判断がＹｅｓの場合には、ステップＳ４５に進む。ステップＳ４５においては、stacknodeのstackptr2の指す位置から構成要素ノードを取り出す。取り出したらstackptr2＝stackptr2＋１とカウントアップする。ステップＳ４６において、スタックから取り出したノードの隣接ノードをスキャンして、領域番号が現在処理している分割領域に属し、また、スタックに積まれていない構成要素ノードがあれば、stackptr＝stackptr＋１として、stacknodeに積む。
【００８９】
ステップＳ４７において、stackptr2＞stackptrか否かを判断する。ステップＳ４７の判断がＮｏの場合には、ステップＳ４５に戻る。ステップＳ４７の判断がＹｅｓの場合には、ステップＳ４８において、counterという変数に、今までに既に番号付けされたノードの総数を設定する。ステップＳ４９において、stacknodeの上の方から構成要素ノードを一つずつ取り出し、counter＝counter＋１を演算後、構成要素ノードの新番号の領域（図７参照）に、この値を格納する。stacknodeの上の方から取り出すとは、スタックに最後に積まれた構成要素ノードから取り出すという意味である。ステップＳ５０において、stacknodeが空か否かを判断する。ステップＳ５０の判断がＮｏの場合には、ステップＳ４９に戻り、Ｙｅｓの場合には、サブルーチンを抜ける。
【００９０】
図１６は、図１４のステップＳ３５の詳細フローであり、分割面を構成する構成ノードに番号を振るサブルーチンの処理フローである。
ステップＳ５１において、counterに今まで番号付けした総ノード数を代入する。ステップＳ５２において、現在処理している分割面の構成要素ノードを、１つずつチェインをたどって取り出し、取り出された構成要素ノードの新番号として、counter＝counter＋１を格納する。ステップＳ５３において、stacknodeが空か否かを判断し、Ｎｏの場合には、ステップＳ５２に戻り、Ｙｅｓの場合には、サブルーチンを抜ける。
【００９１】
図１７〜図２０は、図１３〜図１６に対応する方式２のフローチャートである。
領域ノードの集合を表すtreeをrootからdepth-first searchの順にたどって番号付けする。depth first searchでleaf方向にtreeをchild 関係でたどるとき、分割面を表すノードを通過する。このとき、分割面の構成要素ノードの隣接ノードで分割領域に属する構成要素ノードを対応する分割領域のstacknode に積む。そして、分割領域の構成要素ノードをbreadth first searchで番号付けする出発集合とする。分割面の構成要素の番号付けは、方式１と同じ。この領域ノードからなるtreeは２分木で、leaf nodeの深さレベルは同じで、brother はあれば必ず１つとなる。
【００９２】
分割面の構成要素ノードの番号付けのルーチンは方式１と同じである。領域ノードからなる２分木の構成要素ノードの情報に通過したか否かを示すflagを持たせる。
図１７及び図１８において、まず、ステップＳ５５で、curparentに、領域ノードからなるtreeのrootノードを設定する。また、curparentのchildをcurchildに設定する。このとき、curparentの他のchildをbrotherに設定する。ステップＳ５６において、curchildがnullか、すなわち、存在しないか否かを判断する。ステップＳ５６の判断がＮｏの場合には、ステップＳ５７において、curparentをまだ通過していない、すなわち、まだ処理していない場合には、サブルーチン（後述）において、この領域ノードを表す分割面の構成要素ノードの隣接ノードが分割領域に属するとき、stacknodeに当該領域の構成要素ノードを積み、通過済みのflagを立てる。
【００９３】
ステップＳ５８において、curparent＝curchildとする。また、左記の様に更新されたcurparentのchildをcurchildに設定する。そして、curparentの他のchildにchildのbrohterを設定し、ステップＳ５６に戻る。
【００９４】
ステップＳ５６の判断がＹｅｓの場合には、ステップＳ５９において、grandpにcurparentのparentを設定する。ステップＳ６０において、cutsurface＝grandp、nodeset＝parentとして、分割領域の構成要素ノードの番号付けを行なうサブルーチンを呼び出す。ステップＳ６１において、nodesetにgrandpのchildを設定し、grandpの他のchildに、そのchildのbrotherを設定する。ステップＳ６２において、設定されたcutsurface、nodesetで、分割領域の構成要素ノードの番号付けを行なうサブルーチンを呼び出す。
【００９５】
ステップＳ６３において、curparentにcurparentのparentを設定する。ステップＳ６４において、curparentのchildがnullか、すなわち、存在しないか否かを判断する。ステップＳ６４の判断がＮｏの場合には、ステップＳ６９において、curparentのchildをcurchildに設定し、ステップＳ５６に戻る。
【００９６】
ステップＳ６４の判断がＹｅｓの場合には、ステップＳ６５において、nodesetにcurparentを設定し、ステップＳ６６において、cutsurfaceの構成要素ノードを番号付けするサブルーチンを呼び出す。ステップＳ６７において、curparentにcurparentのparentを設定する。すなわち、領域分割を表すtreeをroot方向に戻る。ステップＳ６８において、curparentがnullか、すなわち、存在しないか否かを判断する。ステップＳ６８の判断がＮｏの場合には、ステップＳ６４に戻り、Ｙｅｓの場合には、処理を終了する。
【００９７】
図１９のサブルーチンは、図１７のステップＳ５７のサブルーチンであり、depth first searchで、領域ノードからなるtreeをleafの方向にたどるとき、最初に、このノードを通過したとき呼び出される。
【００９８】
そのとき、分割面の構成要素ノード取り出して、その隣接ノードが分割領域に属するとき、分割領域に対応する、leafノードごとにstacknode を持つ。各スタックにノードを積む位置は、番号付け処理の最初にクリアしておく（stacknode ごとにstackptr=0）。
【００９９】
図１９において、ステップＳ７０で、分割面を構成する構成要素ノードのチェインをたどって、構成要素ノードをstacksfに積む。ステップＳ７１において、stacksfから１つのノードを取り出す。ステップＳ７２において、現在の構成要素ノードの全ての隣接ノードを順番に取り出し、同じ分割領域に属しているか否かを調べる。対応する分割領域のstacknodeに積まれていなければ、その分割領域のstacknodeに積む。対応するstacknodeに積む位置は、対応するポインタを、stackptr＝stackptr＋１と更新して、積む位置を計算する。新たに構成要素ノードを積むときは、構成要素ノードの処理済flagをonにする。ステップＳ７３において、stacksfが空か否かを判断し、Ｙｅｓならばサブルーチンを抜け、Ｎｏならば、ステップＳ７１に戻る。
【０１００】
図２０は、図１８のステップＳ６０、Ｓ６２、Ｓ６６で呼び出されるサブルーチンの処理フローであり、stacknode に積まれ、breadth first searchのための始点ノード群をスタックの底からたどりながら、その隣接ノードで同じ分割領域に属するノードで、まだスキャンされていないノードをstacknode に積む処理を行う。
【０１０１】
図２０において、ステップＳ８０においては、stackptr2＝１とポインタを初期化する。ステップＳ８１において、stacknodeのstackptr2の指す位置からノードを取り出す。取り出したら、stackptr2＝stackptr2＋１とカウントアップする。ステップＳ８２において、取り出したノードの隣接ノードをスキャンして、領域番号が同じ分割領域に属し、まだスタックに積まれていないノードがあれば、stackptr＝stackptr＋１として、stacknodeに積む。ステップＳ８３において、stackptr2＞stackptrであるか否かを判断する。ステップＳ８３の判断がＮｏの場合には、ステップＳ８１に戻る。ステップＳ８３の判断がＹｅｓの場合には、ステップＳ８４において、counterに今までに番号付けした総ノード数を設定する。ステップＳ８５において、stacknodeの上の方から構成要素ノードを一つずつ取り出し、構成要素ノードにcounter＝counter＋１をノードの番号として、新番号の領域に格納する（図７参照）。stacknodeの上の方から取り出すとは、スタックに最後に積まれた構成要素ノードから取り出すという意味である。ステップＳ８６において、stacknodeが空か否かを判断し、ＮｏであればステップＳ８５に戻り、Ｙｅｓのであれば、サブルーチンを抜ける。
【０１０２】
図２１は、本実施形態が適用される共有メモリ型スカラ並列計算機のハードウェア構成を説明する図である。
マルチコアＣＰＵ１０は、１つのＣＰＵノードの内部にＣＰＵコア（Ｌ１キャッシュも内臓）１２を複数封入したもので、各コアからL2キャッシュおよびバスインタフェース１１を共通に使う。論理的には、各ＣＰＵコアが１つのＣＰＵとして動作するように見える。
【０１０３】
図２１の例では、１つのＣＰＵノードに２つのＣＰＵコアが封入された場合を示している。１つのＣＰＵノードに４つの以上のＣＰＵコアが封入されるものもある。ＣＰＵノードが一つで、8 個のＣＰＵコアがそのノードに封入されていてもSMP（Shared Memory Palallel computer)と見なせる。ここでＣＰＵと呼ぶのは、ＣＰＵコアのことである。ＣＰＵには各スレッドが割り当てられて処理が行われるが、スレッドは、ＯＳが並列処理を行う処理単位のことである。スパース行列の要素データのうち、処理に必要なデータは、メモリモジュール１４に分散されて保持され、相互結合網（バス）１３を介して、L2キャッシュに格納される。そして、L2キャッシュから、スーパーノードを構成する、スパース行列の構成要素ノードのデータがL1キャッシュにロードされ、ＣＰＵコアによって演算される。スーパーノードの幅が大きいほどＣＰＵコアの演算量に比べ、L1キャッシュへのデータのロードストア回数を減らすことができるので、高速に演算を行うことができる。
【０１０４】
上記実施形態のアルゴリズムは、プログラムとして実装することが出来る。プログラムは、記録媒体１６に格納され、記録媒体読み取り装置１５によって読み取られる。記録媒体読み取り装置１５によって読み取られたプログラムは、メモリモジュール１４のいずれかに格納される。記録媒体１６は、ＣＤ−ＲＯＭ、フレキシブルディスク、ＤＶＤ、Ｂｌｕ−Ｒａｙなどの持ち運び可能な媒体などである。プログラムの実行は、メモリモジュール１４に格納された状態で実行されても良いし、記録媒体１６から読み取りながら実行しても良い。
【０１０５】
以上の実施形態は、シミュレーションや数理モデルから生じるスパースな正値対称行列で表される連立1次方程式の解を求めることで、数理モデルの分析を行う分野、有限要素法、偏微分方程式などの解析を行いシミュレーションを行う分野に適用可能である。
【０１０６】
上記実施形態のorderingによる行列の演算と、スパース行列を、グラフの分割に基づいた一般化されたnested dissection により生成されたorderingで並べ変えた行列の演算と比べる。
【０１０７】
並列化されたコレスキー分解において、本実施形態のorderingによるコレスキー分解を行う性能は、比較対象の方法で並び替えた行列に対してコレスキー分解を行う性能に対し、7 倍程度高速なものもある。
【０１０８】
本実施形態のorderingと、もう分割できないまで完全に分割を繰り返して生成されたorderingである一般化されたnested dissection ordering(METIS V4.0で生成)と比較を行った場合以下のような結果であった。
【０１０９】
＜１．orderingの説明＞
３次元の立方体でラプラス方程式を有限差分で離散化して得られる方程式を各次元を60分割した。つまり60^3個の変数が現れる場合を用いた。１つのレベルのnested dissectionで、各次元を２分割して境界面を取り除くと、8つの小さな立方体が生成される。これを3レベル行った。生成された8^3個の小さな（ほぼ）立方体について、境界面から遠い順に番号付けを行うことでorderingを得る。
【０１１０】
＜２．データ＞
スレッド数 1 2 4 8 16 32 64 128
３レベル 74.667 40.874 28.432 15.606 10.955 7.42 6.407 6.133
ND 359.599 198.51 179.318 105.55 78.111 61.029 45.925 40.959
【０１１１】
＜３．データの解説＞
計算機は共用メモリ型スカラ並列計算機M9000 (128cores, SPARC64VII quad core CPU 2.52GHz)を使用し、この計算機上でスパースな正値対称行列の連立１次方程式を以下のorderingでスパース行列のLDL^T分解によって解いた性能を示している。スパース行列の次数は60^3であった。単位は秒(sec)である。
【０１１２】
＜orderingの方法＞
“3レベル”は、本実施形態で示したものを使用し、“ND”は、一般化されたnested dissection ordering（METIS 4.0でサポートされているもの）を利用した。
【０１１３】
16〜64スレッドでの並列実行で本実施形態は、7倍以上高速である。
【符号の説明】
【０１１４】
１０ ＣＰＵノード
１１ L2キャッシュ及びバスインタフェース
１２ L1キャッシュ及びＣＰＵコア
１３ 相互結合網（バス）
１４ メモリモジュール
【技術分野】
【０００１】
スパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行なう場合における、orderingを生成する方法に関する。
【背景技術】
【０００２】
スパースな正値対称行列によって表される連立１次方程式（Ax=b、Aはスパースな正値対称行列、xは変数列ベクトル、bは定数列ベクトル）の解法において、変数列ベクトル内の変数の配置を並び替えることで、対応する行列のnonzero要素の位置も変わる。この並び替えをorderingと呼ぶ。適当なorderingを行うことでLDL＾T 分解のfill-in（元の行列ではzeroであったが、分解の結果nonzeroとなる行列要素） を削減することが知られている。LDL^T分解においては、スパース行列AをA=LDL^Tと分解して、連立１次方程式Ax=bを解く。ここで、Lは、下三角行列、^Tは、行列の転置を取ることを示す。
【０００３】
スパースな正値対称行列の対角要素を含む下三角行列部分L の列に関する、LDL^T分解での、行列要素間のデータ更新の依存関係はelimination treeと呼ばれる木構造で表される。elimination treeについては、非特許文献１を参照されたい。この依存性からleaf(elimination treeを構成する構成要素ノードのうち、子ノードを持たないノード）からroot（親ノードを持たないノード）にたどる方向に分解を行う。
【０００４】
elimination treeの互いに共通部分のない部分木は独立に分解の計算ができ、並列に処理することができる。
代表的なorderingにnested dissection がある。nested dissection とは、空間を２つの分割領域とこれらの分割領域を分割する分割面に分割して番号付けする方法である。nested dissectionについては、非特許文献１を参照されたい。さらに、スパース行列をグラフで表現して、グラフを構成ノード数がほぼ等しい２つの部分グラフに分割し、分割するために取り除くノード数が最小となるように分割をする。この2分割を、再帰的に繰り返し、もうこれ以上分割できないところまで分割して得られるorderingは一般化されたnested dissection によるorderingと呼ばれる。
【０００５】
分割は、分割面を形成するノードと分割された分割領域のノード集合に分かれる。あるレベルの分割に関して、分割面で分割された２つの集合をＡ，Ｂとしたとき、おのおのの領域をさらに分割しないようにＡに属する構成要素ノードに連続な番号を振ったあと、Ｂに属する構成要素ノードに連続な番号を振って、最後に分割面を構成する構成要素ノードに番号を連続に振る。この番号付けのルールは、再帰的に分割した部分にも適用する。
【０００６】
この結果生成された、nested dissection のorderingは、fill-in が少ないorderingであることが知られている。
上記では、LDL^T分解について述べたが、この分解は、修正コレスキー分解と呼ばれ、コレスキー分解（LL^T分解）とほぼ同じアルゴリズムが適用できる。したがって、以下の説明では、コレスキー分解について言及するが、同様のアルゴリズムが修正コレスキー分解にも適用できる。
【０００７】
スパースな正値対称行列の連立1次方程式をコレスキー分解で解くとき、行列の行あるいは列の並べ替えを行うことで、fill-in を少なくすることが一般的に考えられる。この並び替えをorderingと呼ぶ。並べ替えた行列をコレスキー分解する。
【０００８】
このコレスキー分解は、symbolic decomposition とnumeric decomposition の２つのフェーズから成る。symbolic decompositionで、各列の依存関係を解析し、elimination treeに表現する。この情報を利用して、スパース行列に存在するnonzero 要素の発生パターンを解析する。
【０００９】
コレスキー分解におけるnonzero要素の発生パターンは、elimination treeを生成して解析する。Ax=b（A=(aｉｊ））という連立方程式をコレスキー分解で解く場合、下三角行列をLとして、A=LL^Tと分解する。elimination treeは、以下の定理を用いて、最初の構成ノードから、子の構成ノードを見つけて連結することによって生成する。
・定理
コレスキー分解LL^T=Aにおいて、数値的なキャンセルを除くと、aki≠０、かつ、ｋ＞ｉであれば、ノードｉは、elimination treeにおけるノードｋの子ノードである。
（非特許文献１参照）
【００１０】
なお、ノードｋをノードｉの親ノードと呼び、ノードｉに対応するデータは、行列Aのｉ番目の列のデータである。また、elimination treeの構成要素ノードは、対応する番号の格子点に対応する。したがって、構成要素ノードｉは、離散化された空間の格子点ｉに対応し、構成要素ノードｉのデータは、行列Aの列ｉのデータである。
【００１１】
大規模なスパース正値対称行列のコレスキー分解を並列処理して性能を引き出す上で、elimination treeのお互いに交わりのない部分木は独立に計算できるので、独立な部分木については並列に計算を行なうことが考えられる。これ以外にも、elimination treeで示されるスパース行列の要素の更新の影響の伝搬関係に従い列を更新する処理を、順次パイプライン処理することで並列性能を引き出すことも考えられる。また、列を束ねてブロックを構成して計算量を大きくすることが並列効果を高める上でも重要である。一回のデータの読み込みで多くのデータを読み込んでおくと、速度の遅い、データの読み込み／書き出し処理の回数を減らすことができるので、並列処理効果が上がるのである。そのため、nonzero 要素発生パターンが同じ列をまとめてスーパーノードを構成する。スーパーノードを構成した場合には、スーパーノードを構成する列を結合したとき、nonzero な値を持つ行のデータのみを保持するようにする。したがって、zeroとなる行のデータを保持しなくて良くなるので、データを圧縮して格納して処理することができる。
【００１２】
さらに、コレスキー分解の性能を向上するために、条件を弱めて、分解の結果nonzeroパターンの似ている列も結合して、スーパーノードを構成する。この結合の結果、コレスキー分解の更新過程で、スーパーノードのデータ内に余分なzeroを含む行が発生し、演算量が増加する。しかし、単位演算量あたりのメモリアクセスを減らすことができる。このため、コレスキー分解の性能を向上することが期待できる。コレスキー分解の並列処理は、理論的には１列ごとの計算でも並列に行えるが、メモリアクセスの量や演算数が小さいため、並列処理のオーバーヘッドが無視できず、並列処理効果が引き出しにくい。
【００１３】
つまりメモリアクセスに対する演算数を上げること、および、並列処理のオーバーヘッドを削減するために粒度（granularity）を大きくすることが効率的な並列化を行うためには有効である。
【００１４】
理論的には、完全なnested dissection （空間の分割をこれ以上分割できないまでやって番号付けする方法）によるorderingを利用して行列を並び替え、コレスキー分解を、列を結合せずに、１列ごとに計算を行うとき、fill-in が最小になる。性能引き出しのために、１回のキャッシュへの読み込みで読み込まれるデータとして、列を結合してスーパーノードを作るが、elimination treeのleafに当る列から列の結合を行う処理を進める。nested dissection のorderingを利用してコレスキー分解で列の結合を行うとき、列を結合したブロックの幅を増やすためには、elimination treeの多くの分枝を有するノードを結合する必要がある。
【００１５】
分枝先ノードは分枝元ノードのnonzero要素とは共通のnonzero要素を多く持つが、分枝先のノード間では共通しないnonzero要素を多くもつ場合がある。このため多くのノードを結合すると、分枝ノードをかなり含むとき、zeros の数が急に増えることになる。そのため、ノードを順次結合して、zeroを少しずつ増やしながら、スーパーノードの適当なブロック幅を見つけることが難しいという問題点があった。
【００１６】
fill-in を削減する効果を期待してnested dissection orderingを使う場合、列の結合を行わないと、nested dissection によって発生したelimination treeの２分木構造のleafに相当する部分に、列一本に対して多くの計算が発生する。fill-in は少なくなるが、メモリアクセスの演算に対する割合や、並列処理の粒度が小さいことによるオーバーヘッドが顕在化してくるため、期待する性能が得られないことになる。
スパース行列の直接解法に使われる基本技術の詳細は、非特許文献１を参照されたい。
【００１７】
＜nested dissection から生じるorderingとそのelimination tree＞
２次元の正方形領域での問題を考える。この領域での問題を表現する偏微分方程式を離散化してスパースな正値対称行列を得る場合を考える。このとき、離散化された格子点は、領域の中で、隣接する格子点と結合されている。
【００１８】
図１のような、x-y 平面で考える。この領域を構成する格子点からｙ軸と平行な直線を構成する格子点の集合ａを取り除くことで、領域を構成する格子点の集合を２分割して、おのおのが構成する格子点の数がほぼ等しくなるようにする。
【００１９】
残った２つの領域のおのおのを、ｘ軸と平行な横方向の格子点の集合ｂを取り除くことで、同様に２分割して、おのおのが構成する格子点の数がほぼ等しくなるようにする。格子点の数がほぼ等しくなるように分割するには、今分割しようとする領域の全格子点数を２で割った数の格子点からなるように領域を分割すればよい。この手順を、分割によって生成された小さくなった領域が、これ以上分割できないようになるまで再帰的に分割を続ける。これ以上分割できないとは、分割して得られた小領域に含まれる格子点の数が３以下となって、次の４分割ができないような場合を示す。最終的に領域は2**nの小さな領域に分割される。
【００２０】
この処理は、領域をほぼ均等に分割して、取り除く境界線を構成する格子点数が少なくなるように選ぶことに相当している。
番号付けは、分割で生成した２つの領域の内の一つの領域を構成する格子点の集合に対して連続に番号を割り当てる。次に、もう一つの領域を構成する格子点の集合に対して引き続いて番号を連続に割り振る。最後に、２分割で取り除いた格子点に引き続き番号を割り振る。
【００２１】
このルールを再帰的に適用して、上記分割から生じた集合に番号付けを行う。この番号づけで並べた行列のelimination treeは２分木になる。
【００２２】
すなわち、例えば、離散化された空間上でラプラスの方程式（Δφ＝）を解くことを考える。一次元の問題であるとすると、ラプラシアンは、座標の２階偏微分になるので、φ（ｉ）を格子点ｉにおける関数値とすると、離散化されたラプラスの方程式は、
｛φ（ｉ＋１）＋φ（ｉ−１）−２φ（ｉ）｝／ｈ２＝０
となる。ここで、ｈは、隣り合う格子点間の距離である。これを行列で表すと、
【数１】
のようになる。これから明らかなように、隣り合った格子点が連続的に番号付けされている場合には、スパース行列内では、係数は隣接して発生する。しかし、格子点を上記のように分割すると、分割された領域内は、隣り合った格子点に連続した番号が現れるので、係数は固まって発生するが、分割して切り離され、別個に番号付けされた格子点に対応する係数は、他の係数とは離れた場所に発生する。したがって、このような行列のelimination treeを生成すると、係数が固まっている列同士に対応するツリーは互いに結合されるが、離れてしまった列同士は別のツリーとなる。したがって、領域を、ＡとＢ及びこれらの境界領域Ｃとに分割した場合、Ａ内の格子点に対応するelimination treeの構成要素ノードと、Ｂ内の格子点に対応するelimination treeの構成要素ノードとは、境界領域Ｃの格子点に対応するelimination treeの構成要素ノード以降において、２つに分枝することになる。
【００２３】
この２分木は、分枝するノードとノードが直線に並んだ枝から成る。
以下に、再帰呼び出しの手続きでの番号付けの方法を記述する。
以下のプログラムコードにおいて、関数partition は領域の分割を行う手続きで、領域を分割できるかを判定して、できるときは、seta,setb に２分割した領域を、cutsurfaceに分割面を返す。setは、今分割したい格子点の集合である。
call numbering(set)
recursive procedure numbering(set)
call partition(set,seta,setb,cutsurface,icon)
if(partition is done)then ! numbering を呼ばれた順に番号付けは連番で行う。
call numbering(seta) ! setaを番号付けする。
call numbering(setb) ! setbを引き続き番号付けする。
call numbering(cutsurface) ! cutsurfaceを番号付けする。
else
call numbering(set)
endif
return
end
なお、partitionのアルゴリズムについては、非特許文献２を参照されたい。
【００２４】
３次元問題に関しても同じように考えることができる。例えば、立方体ならy-z 平面に平行な面を構成する格子点の集合を取り除いて、領域をほぼ均等になるように分割する。引き続き、分割で生成された２つの領域を同様にx-z 平面に平行な面でおのおのを２分割する。これらを２次元の場合と同様に、これ以上分割できない状態まで繰り返す。
【００２５】
番号付けのルールも同じである。この番号付けで並び替えた行列のelimination treeはやはり２分木になる。
これらを、一般の連結する領域に対して適用するために、離散化した格子点の結合関係をグラフで表現し、グラフを分割する面を構成するノード数が最小で、分割した部分グラフを構成する構成要素ノードの数が均等になるように分割することを繰り返して、一般化されたnested dissection orderingを構成することが考えられている。この場合も同様に、elimination treeは２分木になる。
【００２６】
＜ノードの結合（スーパーノード）＞
elimination treeの構成要素ノードに対応する行列の列で、引き続くノード番号の構成要素ノードに関して、ノード番号の大きい方のnonzeroパターンがノード番号の小さい列と等しい列を結合する。また、枝を構成する構成要素ノードを調べてnonzero パターンが比較的に似ていて、結合してもzeroの割合が大きく増えない場合は結合する。この結果、コレスキー分解すると、結合した列ベクトルの行で、nonzero 要素を持つ行の数は増える。つまり、この行を圧縮して格納するのに必要な領域にzeroが含まれることになる。こうしてnonzeroパターンが似ていて列を結合してできたブロックをpanel と呼ぶ。
【００２７】
図２は、panelとindex listを示す図である。
panelには、index listが対応して設けられる。index listは、panelのデータがスパース行列の列の中の何行目に当るかを保持する。panelは、panelに保持される列に共通にzeroである要素以外のnonzero要素のデータしか保持しないので、index listは保持されているnonzero要素が列の中の何行目かを示すのである。なお、panelは、保持される列に共通にzeroとなる要素を保持しないものであるので、nonzero要素として保持されていても、実際の数値はzeroとなることがある。
【００２８】
panel ごとのコレスキー分解の処理は、ブロック幅が大きくなるほど演算量が増えるため、相対的にメモリのロードストアのコストが下がる。ただし、zero要素の部分は計算する必要のない無駄な計算が発生する。このため、結合によるブロック幅を、含まれるzeroの数を考慮しながら調整して、ブロック幅を決定する必要があり、このことが性能を引き出す上で重要である。
【００２９】
スーパーノードはelimination treeにおいて連続する構成ノードで、nonzero 要素のパターンが同じものを結合して、nonzero 要素がある行のみ集めて圧縮した形で格納する。すなわち、nonzero要素の行のみデータを保持し、zero要素の行のデータを保持しないことによって、保持すべきデータの量を圧縮するものである。
【００３０】
図３は、panelに不要なzeroを含む様子を示した図である。
elimination treeの親の構成要素ノードでnonzero 要素のパターンが一致するものを結合したものを、各列のnonzero 要素のみ圧縮して格納する。これはpanel となる。次の親ノードとはnonzeroパターンは一致しないが、包含関係にある。すなわち、次の親ノードにおいてnonzeroである要素に対応する行は、分解の処理における更新処理において更新され、nonzero要素となるので、次の親ノードのnonzeroパターンは、panel内に反映されることになる。このような関係を包含関係と呼んでいる。このような次の親ノードをpanelに結合すると、不要な0 をpanel に含むようになる。図３の例の場合、panelの下方に不要なzeroが生じており、この場合のインデクスは、panelと、結合された親ノードとで共通でない部分が大きい。
【００３１】
＜分枝するノードを結合すると大きくzeroが急激に増える。＞
分枝するノードから分かれる枝の先につながるノードは、nested dissection では分割された領域を構成するノードに対応する。
【００３２】
コレスキー分解の結果、ノードに隣接するノードは、ノードに対応する列のnonzero要素として現れる。コレスキー分解の過程で、elimination treeのrootへのパスをたどるとき、nonzero要素は、他の列への更新の際に影響を与える要素としてその影響が伝播することになる。つまり、ある面で分割された２つの領域の、この面以外で隣接する面のノードはnonzero要素として分枝ノードに伝わる。２分割されたノード群は分割面以外の境界面は共通要素を持たないので、これに隣接する領域外のノードに対応するnonzero要素は分枝ノードに受け継がれる。すなわち、分枝ノードにおいて、zeroである要素も、結合することによって、nonzero要素として扱う必要が生じる。つまり、分枝ノードを結合すると、分枝ノードの枝を構成するノードを結合したパネル幅に、zeroが急激に増加する。
【００３３】
nested dissectionでは、elimination treeの構成要素ノードをleafから結合していくと多くの分枝ノードに遭遇するためzeroの数とpanel の幅を調整することが難しい。
図４のような直方体がいくつかの分割の結果、生成されたとする。真ん中の部分は直方体を２分して取り除く面を構成するノードの集合を表す。分割前はこれら３つの部分はつながっている。また、分割前の表面は、他の領域と接していると仮定する。
【００３４】
分割領域内のノード（格子点）に番号付けするとき、領域を分割しないように適当に番号付けする。このときelimination treeでは、これらはY をさかさまにしたような木の部分として表現される。
【００３５】
分枝ノードで分かれる部分を考えると、elimination treeは、分枝ノードから下に２又に分かれる。コレスキー分解した結果の下三角行列L の構造は、親ノードのnonzero パターンの構造に包含される。そのため、ブロックA とC 以外で接する他の領域に存在するノードのnonzero は、出現した位置から伝播していく。
【００３６】
図５に示される列データを用いて説明すると、分枝ノードy に対し、分枝先の２つのノードz,w のnonzero パターンは図５の例のようになる。すなわち、ノードzは、対角要素のほかに、ブロックCの構成ノードから発生するnonzero要素ＣＣＮと、ブロックAのブロックC以外と接する外部ノードから発生するnonzero要素ＡＮからなる。ノードwは、対角要素のほかに、ブロックCの構成ノードから発生するnonzero要素ＣＣＮと、ブロックBのブロックC以外と接する外部ノードから発生するnonzero要素ＢＮとからなる。一方、分枝ノードyは、対角要素のほか、ブロックCの構成ノードから発生するnonzero要素ＣＣＮと、ブロックAのブロックC以外と接する外部ノードから発生するnonzero要素ＡＮと、ブロックBのブ、ロックC以外と接する外部ノードから発生するnonzero要素ＢＮと、ブロックCのブロックA,B以外と接する外部ノードから発生するnonzero要素ＣＮからなる。
【００３７】
分割は領域を、互いに均等な大きさになり、かつ、分割面が最小になるように分割しているとする。つまり、ＡＮはＢＮにほぼ等しく、外部に接する面が多ければ、ＡＮ〜ＢＮ＞ＣＣＮ＞ＣＮとなる。
【００３８】
この結果、図５の右側に示されるように、かなり多くのzeroが急激に増すことになる。
領域の分割をさらに進めてもう分割できないところまで続けたものを結合すると、この構造が繰り返される。
【先行技術文献】
【非特許文献】
【００３９】
【非特許文献１】T.DAVIS, "Direct Methods for Sparse Linear Systems", SIAM 2006
【非特許文献２】George Karypis and Vipin Kumar, "A fast and high quality multilevel scheme for partitioning irregular graphs", SIAM Journal on Scientific Computing, Vol. 20, No. 1, pp. 359 - 392, 1999
【発明の概要】
【発明が解決しようとする課題】
【００４０】
以上のように、スーパーノードの幅を広げて、計算の並列性を高め、高速に処理しようとしても、スーパーノードの幅が広いとスーパーノードに含まれるzeroの数が多くなり、無駄な計算が増えてしまう。これにより、スーパーノードの幅を広げたことによって並列性を高めた効果は少なくなってしまう。zeroは、elimination treeの分枝ノードを結合することによって多く発生する。したがって、elimination treeの形状が、分枝ノードが少ない形状であれば、このような問題は軽減される。ところで、elimination treeの形状は、行列のnonzero要素の配置によって決まり、nonzero要素の配置は、空間の格子点の番号付けの方法による。したがって、適切なorderingによって番号付けられた格子点を用いれば、zeroを急激に増やすことなく、スーパーノードの幅を広げることができる。
【００４１】
以下の実施形態は、スパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行なう場合に、並列処理を高速化することができるorderingを生成する方法を提供する。
【課題を解決するための手段】
【００４２】
以下に開示される実施形態の一側面におけるordering生成方法は、共有メモリ型スカラ並列計算機によってスパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法であって、共有メモリ型スカラ並列計算機は、スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行い、１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う。
【発明の効果】
【００４３】
以下の実施形態によれば、スパースな正値対称行列のコレスキー分解あるいは修正コレスキー分解を行なう場合に、並列処理を高速化することができるorderingを生成する方法を提供することができる。
【図面の簡単な説明】
【００４４】
【図１】従来技術を説明する図（その１）である。
【図２】従来技術を説明する図（その２）である。
【図３】従来技術を説明する図（その３）である。
【図４】従来技術を説明する図（その４）である。
【図５】従来技術を説明する図（その５）である。
【図６】本実施形態で使用されるデータ構造を説明する図（その１）である。
【図７】本実施形態で使用されるデータ構造を説明する図（その２）である。
【図８】本実施形態で使用されるデータ構造を説明する図（その３）である。
【図９】本実施形態で使用されるデータ構造を説明する図（その４）である。
【図１０】本実施形態で使用されるデータ構造を説明する図（その５）である。
【図１１】本実施形態の全体処理のフローである。
【図１２】本実施形態の詳細な処理フロー（その１）である。
【図１３】本実施形態の詳細な処理フロー（その２）である。
【図１４】本実施形態の詳細な処理フロー（その３）である。
【図１５】本実施形態の詳細な処理フロー（その４）である。
【図１６】本実施形態の詳細な処理フロー（その５）である。
【図１７】本実施形態の詳細な処理フロー（その６）である。
【図１８】本実施形態の詳細な処理フロー（その７）である。
【図１９】本実施形態の詳細な処理フロー（その８）である。
【図２０】本実施形態の詳細な処理フロー（その９）である。
【図２１】本実施形態が適用される共有メモリ型スカラ並列計算機のハードウェア構成を説明する図である。
【発明を実施するための形態】
【００４５】
本実施形態においては、大きな粒度の並列性を引き出し、列を結合してスーパーノードの効果的なブロック幅を確保するために、ある深さでグラフの分割を止めて最後に分割面で２分割された２つのノードの集合内で番号を振る。このある深さでグラフの分割をやめるとは、取り除く分割面を構成するノード数が、ノードの結合を行うことを考えたときのスーパーノードのブロック幅の大きさ程度に縮小したかを判断して、その程度に達した場合に、その時点の深さで繰り返し行う再帰的分割を止めるように制御することを意味する。
【００４６】
そして、スーパーノードを構成する場合、elimination treeの枝を形成する連続する列を結合していく。このとき、あるノードを結合したときzeros の増加が急激に大きくならないような並びで番号付けすることを考える。すなわち、スーパーノードを列を結合して形成する場合、分枝ノードを結合するとzeroが大幅に増える。したがって、分枝ノードができるだけelimination treeに現れないようにするのが好ましい。ところで、elimination treeの形状は、行列のnonzero要素が行列内のどの位置に現れるかによって変化する。そして、行列のnonzero要素が行列のどの部分に現れるかは、行列要素の番号付けによって決まり、行列要素の番号付けとは、すなわち、空間を離散化した場合の格子点の番号付けの仕方によって決まる。格子点は、elimination treeの構成要素ノードに対応する。したがって、格子点（構成要素ノード）の番号付けを工夫することにより、zeroを大幅に増やすことなく、スーパーノードの幅を太くすることができる。スーパーノードの幅は、構成要素ノードの接続の仕方によって可変でよいが、上限値を大きく取ることができるようになる。スーパーノードの幅が太いほど並列性の効果を得やすいが、この幅は、当業者が使用する計算機の性能を考慮して、適宜設定すべきものである。例えば、スーパーノードの幅の上限値は、５０〜１００の列を結合した幅とすることが可能である。
【００４７】
このためには、分割を止めて残った領域の番号付けが、さらに領域を分断しないようにすることと、ノードの結合を行うとき、徐々にzeroが増加するような並びにすることとがpanel の幅とzeroの割合を制御し、性能を引き出す上で効果的なスーパーノードの構成を可能にする。空間を分割した結果のブロック内の番号付けでは境界のノードをどう順序付けするかという問題が存在する。
【００４８】
メモリアクセスのスピードに応じて２つの方法が考えられる。
１）ある程度の深さまででグラフの再帰的分割をやめる（取り除く分割面を構成するノード数が、ノードの結合を行うことを考えたときのスーパーノードのブロック幅の大きさ程度となったら分割をやめる）。最後の深さレベルを持つ２分割で生成された領域群のおのおのの領域を構成する構成要素ノードに番号を振るとき、２つに分割する面（分割面）の構成要素ノードからの距離が遠いものから順番に番号を振る方法。
２）もう一つは、生成された分割領域に接する分割面をすべて考えたとき、これらのnested dissection を生成するときに取り除く全ての分割面との距離が遠くなるような構成要素ノードから順番に番号を振る方法。
【００４９】
方式１および２とも、再帰的分割の最も深いレベルでの分割面は、帯型のような形状でnonzero 要素が並ぶ。
方式２の方は、多くの列を結合しなくても性能の引き出しが可能なメモリアクセスが高速である場合に適している。
【００５０】
方式１あるいは２を使用すると、ノードを結合してもzeroの増加はなだらかである。また、zeroが多く発生するのを抑えるためには、分枝ノードを結合することはやめて、elimination treeの直線を形成する枝の部分に着目して、スーパーノードの結合を行えばよい。
【００５１】
図６〜図１０は、本実施形態で使用されるデータ構造を説明する図である。
図６〜図１０は、ノードの隣接ノードを表すデータ構造を表す図である。このデータは、分割面を構成する構成要素ノードの集合、つまり最後に分割されたグラフの構成要素ノードの集合、分割面を表す集合から２つに枝分かれする構造を再帰的に表した木構造（離散化空間を分割した場合の各領域を示す領域ノードをtree構造で示したもの）、各集合に番号を振って識別するとき使う集合の構成ノード数からなる。
【００５２】
なお、以下では、elimination treeを構成するノードを構成要素ノード、領域の分岐の様子を示す木構造を構成するノードを領域ノードという。
図６は、各構成要素ノードに隣接する構成要素ノードを表す。
【００５３】
nadj(*)は、リストベクトルに各構成要素ノードの隣接ノードを格納する。nstart(n+1)は、ｉ番目の構成要素ノードの隣接ノードが１次元配列（nadj(*))の何番目から格納されるかを示す。
【００５４】
図７は、分割された領域を識別するデータを表す。
nodeinfは、分割面か分割された領域を識別する番号を、構成要素ノードのノード番号に対応させて格納する。構成要素ノードごとに、構成ノードポインター(チェイン構造で構成要素を表現する) 、領域識別番号、新番号、処理済flagを持つ。これらの要素をすべてのノード数分の配列の形で持つ。
【００５５】
図８は、分割領域の構造を示すデータ構造である。
領域を分割するときに、領域から分割面を取り除き、共通部分のない２つの領域を生成する。elimination treeにおいて、この２つの領域から分割面の構成要素ノードを経由せずに、他方の領域の構成要素ノードに向かうedgeはない。このような、分割を再帰的に繰り返し、ある深さで止める。ある深さで分割を止めるとは、再帰的分割において、取り除く分割面を構成するノード数が、ノードの結合を行うことを考えたときのスーパーノードのブロック幅の大きさ程度となった深さで分割をやめることを示す。
【００５６】
深さをS としたとき、分割されて残った領域は2**S個、分割面の数は2**S-1である。分割されて残った領域および分割面をノード（領域ノード）として木構造に表現する。最初の分割面をrootとして、２つに分かれた領域をさらに分割する分割面を木の枝として結合する。順次このように結合し、最後の分割でできた分割面は、その子として、２つの分割領域を表すノードを結合する。
【００５７】
図８（ａ）において、各ノード（分割領域を示すノード：領域ノード）のparentは一つレベルの小さな分割の分割面のノードを指す。child はレベルが増した分割面を示すか、最も深いレベルの分割領域を示すノード番号を指す。brother は配列で、child のbrother （深さが同レベルの領域ノード）の番号を指す。
【００５８】
領域ノードの表す分割面または分割領域の領域ノードを構成要素ノード（elimination treeの構成ノード）のチェイン構造の先頭からポインタで表現する。領域ノードが示す領域の構成要素ノードの総数を構成要素ノード総数に格納する。また、分割面および分割領域を識別するための番号をtreeノード番号に格納する。
【００５９】
図８（ｂ）において、図８（ａ）の構造のデータを、配列として保持する。深さレベルをSとすると、分割された領域と分割面のデータを、2**(S+1)-1個の要素を持つ配列として表現する。また、root領域ノードを指すポインタを別個に保持する。
【００６０】
図９は、領域を深さレベル３まで分割した場合の木構造を表す。
図９の番号付けは、treeをdepth first searchの順番で番号を振った。図９において、root領域ノードの番号は１５で、レベル０である。レベル３の領域ノード１、２に接続されるparent領域ノードは、レベル２の３である。レベル３の領域ノード４、５のparent領域ノードは、６である。レベル２の領域ノード３、６のparentノードは、レベル１の領域ノード７である。同様に、領域ノード８、９のparentが領域ノード１０、領域ノード１１，１２のparentが領域ノード１３、領域ノード１０、１３のparentが領域ノード１４である。root領域ノード１５は、領域ノード７、１４のparentである。逆に、領域ノード１、２は、領域ノード３のchildであり、領域ノード４、５は、領域ノード６のchildであり、領域ノード３、６は、領域ノード７のchildであり、領域ノード７、１４は、root領域ノード１５のchildである。また、領域ノード１に対する領域ノード２は、brotherである。以下同様である。
【００６１】
図１０は、番号付けで使うstackのデータ構造を示す。
stacksf およびstcknodeを用意する。スタックには、１次元配列を利用する。方式１では全処理について、それぞれstacksfとstacknodeを一つずつ設ける。方式２では、stacksfは、全処理について１つであるが、stacknode は、分割領域ごとに確保して利用する。
【００６２】
また、番号付けに際し、permutation vectorを設ける。permutation vectorは、ｉ番目の構成要素ノードにｊが振られた場合には、ｊという値に対し、ｉを返すものとする。すなわち、perm(j)=iとなるようにする。
【００６３】
図１１は、本実施形態の全体処理のフローである。
スパース行列をグラフ（elimination tree)に表現して、それを再帰的に分割するが、ノード（列）を結合してスーパーノードを作るときのブロック幅の上限から、これ以上分割するか制御する。これは、分割面を構成するノードを結合して適当な大きさのスーパーノードを構成できるようにするためである。
【００６４】
ステップＳ１０において、スパース行列をグラフ（elimination tree)に表現し、隣接ノードを表現するデータを作成する。ステップＳ１１において、グラフを分割して分割領域Ａ，Ｂおよび分割面Ｃを生成する。そしてＡ，Ｂを表すグラフを作成する。すべての分割面Ｃの構成ノード数の最低値が、スーパーノードのブロック数とほぼ等しいかを判断し、等しいなら、分割をやめて、ステップＳ１２の処理を行う。等しくないなら、すべての分割領域に関してステップＳ１１にもどって分割する。ステップＳ１２において、分割領域、分割面をノードとする２分木を作成し、構成要素をチェインとして格納する。また、２分木の領域ノードに属する構成要素ノードに分割領域、分割面を識別する番号を付与する。ここでは、stacksf、stacknodeを使用する。番号付けの結果は、permutation vectorに格納する。
【００６５】
以下の処理は、方式１と方式２で処理手順が異なる。
方式１の処理手順：
生成された、領域ノードからなるtreeの最大深さレベルをS としたとき、2**(S+1)-1個の領域ノードがある。
root領域ノードからdepth first searchで、順番に処理する。処理しようとするノードがleafノード（子ノードを持たないノード）のとき、その親ノードは分割面を表す。
【００６６】
１）depth first searchで取り出したノードがleaf node であるとき。
leafノードの親ノードである分割面を構成する構成要素ノードのチェインをたどり、構成要素ノードをスタック(stacksf) に積む。
このstacksf から構成要素ノードを一つずつ取り出して、取り出した構成要素ノードの隣接ノードをスキャンする。隣接ノードが現在対象としている分割領域に属する要素ノードであり、かつ、まだstacknode に積まれていないノードであれば、処理済のflagを設定して、stacknode に積む。
【００６７】
この後、stacknode の一番下から、ノードを一つずつスキャンして、スキャンした構成要素ノードの隣接ノードを順番に調べて、対象の分割領域に属し、かつ、まだstacknode に積まれていないなら、処理済みのflagを設定して、stacknode に積む。スキャンする構成要素ノードがなくなるまで、この処理を続ける。
【００６８】
その後、stacknodeのノードを上から下へと取り出して、相対番号を１から振る。全体での番号は、いままでに番号を振った総数＋相対番号になる。これを構成要素ノードごとに設定する。
【００６９】
２）depth first searchで取り出したノードがleafノードでないとき。
分割面を構成する構成要素ノードを表すチェインをたどり、相対番号を１から振る。全体での番号は、いままでに番号を振った総数＋相対番号になる。これをノードごとに設定する。
【００７０】
３）最後に、もとのノード番号の順番でノードの情報が格納されている１次元配列を元のノード番号順にスキャンする。ｉ番目の要素に新たに振られた番号ｊが格納されているとする。permutation を表すベクトルpermにperm(j)=i と格納してpermutation vectorを作成する。
【００７１】
方式２の処理手順：
生成された領域ノードからなるtreeの最大の深さをS とすれば、leafノードはm=2**S個ある。そこで2**S個のstacknode を各leafノードに対して確保し利用する。
そして、treeをdepth first searchでleafへ下るpathをたどる途中で分割面を通過するときに、つぎの処理を行う。
【００７２】
分割面の構成要素ノードを構成要素ノードのチェインをたどって、stacksf に積む。これが終わったら、stacksfから一つずつ取り出して、その隣接ノードをスキャンする。スキャンしたノードがleafノード、つまり最も深いノードで、分割領域の構成要素ノードであるとき、そのleafノードに対応する構成要素ノードをstacknode に積む。そして、そのノードは処理済みflagをonにする。
【００７３】
１）depth first searchで取り出したノードがleafノードであるとき。
depth-first searchで、leafノードに到達したら、stacknode の底から順番に構成要素ノードを取り出す。その構成要素ノードの隣接ノードをスキャンして、このleafノードの構成要素ノードであって、未処理のものについて、処理済のflagをonに設定して、stacknode に積む。この処理を繰り返してstacknode から取り出すノードがなくなったら、stacknode の一番上から下へと順番にノードを取り出して、構成ノードに相対番号を１から振る。
【００７４】
このleafノードまでにそれらの構成要素ノードに番号付けした総ノード数に相対番号を加えて番号を割り当てる。
【００７５】
２）depth first searchで取り出したノードがleafノードでないとき。
分割面を表す構成要素ノードを表すチェインをたどり、相対番号を１から振る。全体での番号は、いままでに番号を振った総数＋相対番号になる。これをノードごとに設定する。
【００７６】
３）最後に、もとのノード番号の順番でノードの情報が格納されている１次元配列を元のノード番号順にスキャンする。ｉ番目の構成要素ノードに新たに振られた番号ｊが格納されているとする。permutation を表すベクトルpermにperm(j)=i と格納してpermutation vectorを作成する。
【００７７】
図１２〜図２０は、本実施形態の詳細な処理フローである。
図１２においては、スパース行列をグラフ（elimination tree）で表現して、再帰的に２分割する。分割のレベルごとにすべての分割面を構成するノード総数の最低値が、スーパーノードのブロック数程度になったか否かを判断して分割を止める。その後、２分木の作成、分割面、分割領域の構成要素ノードのチェインの作成、どの分割面、分割領域に属するかの識別番号の設定を行う。
【００７８】
ステップＳ２０において、行列をグラフに表現する。ここでグラフといっているのは、elimination treeのことである。ステップＳ２１において、領域を表すグラフ（格子点とその接続関係を示す空間のグラフ）を、領域A、領域B、分割面Cに分割する。すなわち、離散化された空間を各分割領域と分割面に分割する。
【００７９】
ステップＳ２２において、直近の分割で生成された分割面の構成要素ノードの総数がスーパーノードのブロック幅（スーパーノードに結合される列数）の上限値程度か否かを判断する。ステップＳ２２の判断がＮｏの場合には、ステップＳ２３において、直近の分割で生じた領域ごとに、グラフ（elimination tree）を作成し、ステップＳ２１に戻る。ステップＳ２２の判断がＹｅｓの場合には、ステップＳ２４において、分割で生成された分割面を親ノードに、２つの分割領域を子ノードに持つ２分木（領域のtree、図８のデータ構造で、領域ノードを構成ノードとするtree）を作成する。各分割領域、分割面の構成要素ノードをチェインで関連付け、分割領域、分割面と対応付ける。この対応付けが識別できるように構成要素ノードに識別番号を対応付ける。すなわち、elimination treeの各構成要素ノードがどの領域に属するかを示す領域識別番号を、構成要素ノードごとに、図７のデータ形式の領域識別番号領域に入力する。
【００８０】
図１３〜図１６は、方式１の処理フローチャートである。
図１３及び図１４では、elimination treeの構成要素ノードの集合であって、分割面、分割領域を表す領域ノードのtreeをrootからdepth-first searchの順にたどって番号付けする。このtreeは２分木で、leaf ノードの深さレベルは、同じ。brother はあれば必ず１つとなる。
【００８１】
ステップＳ２５において、curparentという変数に、rootノードを設定する。また、curparentのchildをcurchildという変数に設定する。そして、curparentのもう１つのchildをbrohterに設定する。この場合、１回の分割で生じた、もう一つの領域をbrotherとして設定する。
【００８２】
ステップＳ２６において、curchildがnull、すなわち、存在しないか否かを判断する。ステップＳ２６の判断がＮｏの場合、ステップＳ２７において、curparentにcurchildを設定し、curparentのchildをcurchildを設定する。そして、curparentのchildに、brotherを設定し、ステップＳ２６に戻る。ステップＳ２６の判断がＹｅｓの場合には、ステップＳ２８において、grandpという変数に、curparentのparentを設定する。ステップＳ２９において、cutsurfaceにgrandpを設定し、nodesetにparentを設定し、分割領域の構成要素ノードの番号付けを行なうサブルーチンを呼び出す。
【００８３】
ステップＳ３０において、nodesetにgrandpのchildを設定し、grandpのchildに、brotherを設定する。ステップＳ３１において、設定されたcutsurface、nodesetで、分割領域の構成要素ノードの番号付けのサブルーチン（後述する）を呼び出す。ステップＳ３２において、curparentに、curparentのparentを設定する。
【００８４】
ステップＳ３３において、curparentのchildがnullか、すなわち、存在しないか否かを判断する。ステップＳ３３の判断がＮｏの場合には、ステップＳ３８において、curparentのchildをcurchildに設定し、curparentのもう１つのchildに、brotherを設定して、ステップＳ２６に戻る。
【００８５】
ステップＳ３３の判断がＹｅｓの場合には、ステップＳ３４において、nodesetに、curparentを設定する。ステップＳ３５において、cutsurfaceの構成要素ノードを番号付けするサブルーチンを呼び出す。ステップＳ３６において、curparentにcurparentのparentを設定する。この動作は、領域ノードからなる領域の分割を表すtreeをrootの方向に戻ることを意味する。ステップＳ３７において、curparentがnullか、すなわち、存在しないか否かを判断する。ステップＳ３７の判断がＹｅｓの場合には、全ての領域ノードについて処理し終わったので、処理を終了する。ステップＳ３７の判断がＮｏの場合には、ステップＳ３３に戻る。
【００８６】
図１５は、図１３のステップＳ２９、Ｓ３１のサブルーチンの処理フローである。
図１５においては、分割面を構成する構成要素ノードの集合から遠い構成要素ノードから、領域ノードの構成するtreeのleafに対応する分割領域の構成要素ノードに番号を振る。
【００８７】
ステップＳ４０において、図１０で示されたstacksfとstacknodeをクリアする。また、スタック内の位置を示すポインタstackptrを０に、stackptr2を１に初期化する。ステップＳ４１において、分割面を構成する構成要素ノードのチェインをたどって、たどった先の領域に属する構成要素ノードを取り出し、取り出したノードをstacksfに積む。ステップＳ４２において、stacksfから１つの構成要素ノードを取り出す。ステップＳ４３において、構成要素ノードの全ての隣接構成要素ノードを順番にスキャンし、領域番号を調べて、現在処理している分割領域に属しているが、stacknodeにまだ積まれていない構成要素ノードをstacknodeに積む。このとき、stackptr＝stackptr＋１として、構成要素ノードを積む位置を計算する。そして、新たに、構成要素ノードをスタックに積むとき、構成要素ノードの処理済flagをonにする（図７参照）。
【００８８】
ステップＳ４４において、stacksfが空か否かを判断する。ステップＳ４４の判断がＮｏの場合には、ステップＳ４２に戻る。ステップＳ４４の判断がＹｅｓの場合には、ステップＳ４５に進む。ステップＳ４５においては、stacknodeのstackptr2の指す位置から構成要素ノードを取り出す。取り出したらstackptr2＝stackptr2＋１とカウントアップする。ステップＳ４６において、スタックから取り出したノードの隣接ノードをスキャンして、領域番号が現在処理している分割領域に属し、また、スタックに積まれていない構成要素ノードがあれば、stackptr＝stackptr＋１として、stacknodeに積む。
【００８９】
ステップＳ４７において、stackptr2＞stackptrか否かを判断する。ステップＳ４７の判断がＮｏの場合には、ステップＳ４５に戻る。ステップＳ４７の判断がＹｅｓの場合には、ステップＳ４８において、counterという変数に、今までに既に番号付けされたノードの総数を設定する。ステップＳ４９において、stacknodeの上の方から構成要素ノードを一つずつ取り出し、counter＝counter＋１を演算後、構成要素ノードの新番号の領域（図７参照）に、この値を格納する。stacknodeの上の方から取り出すとは、スタックに最後に積まれた構成要素ノードから取り出すという意味である。ステップＳ５０において、stacknodeが空か否かを判断する。ステップＳ５０の判断がＮｏの場合には、ステップＳ４９に戻り、Ｙｅｓの場合には、サブルーチンを抜ける。
【００９０】
図１６は、図１４のステップＳ３５の詳細フローであり、分割面を構成する構成ノードに番号を振るサブルーチンの処理フローである。
ステップＳ５１において、counterに今まで番号付けした総ノード数を代入する。ステップＳ５２において、現在処理している分割面の構成要素ノードを、１つずつチェインをたどって取り出し、取り出された構成要素ノードの新番号として、counter＝counter＋１を格納する。ステップＳ５３において、stacknodeが空か否かを判断し、Ｎｏの場合には、ステップＳ５２に戻り、Ｙｅｓの場合には、サブルーチンを抜ける。
【００９１】
図１７〜図２０は、図１３〜図１６に対応する方式２のフローチャートである。
領域ノードの集合を表すtreeをrootからdepth-first searchの順にたどって番号付けする。depth first searchでleaf方向にtreeをchild 関係でたどるとき、分割面を表すノードを通過する。このとき、分割面の構成要素ノードの隣接ノードで分割領域に属する構成要素ノードを対応する分割領域のstacknode に積む。そして、分割領域の構成要素ノードをbreadth first searchで番号付けする出発集合とする。分割面の構成要素の番号付けは、方式１と同じ。この領域ノードからなるtreeは２分木で、leaf nodeの深さレベルは同じで、brother はあれば必ず１つとなる。
【００９２】
分割面の構成要素ノードの番号付けのルーチンは方式１と同じである。領域ノードからなる２分木の構成要素ノードの情報に通過したか否かを示すflagを持たせる。
図１７及び図１８において、まず、ステップＳ５５で、curparentに、領域ノードからなるtreeのrootノードを設定する。また、curparentのchildをcurchildに設定する。このとき、curparentの他のchildをbrotherに設定する。ステップＳ５６において、curchildがnullか、すなわち、存在しないか否かを判断する。ステップＳ５６の判断がＮｏの場合には、ステップＳ５７において、curparentをまだ通過していない、すなわち、まだ処理していない場合には、サブルーチン（後述）において、この領域ノードを表す分割面の構成要素ノードの隣接ノードが分割領域に属するとき、stacknodeに当該領域の構成要素ノードを積み、通過済みのflagを立てる。
【００９３】
ステップＳ５８において、curparent＝curchildとする。また、左記の様に更新されたcurparentのchildをcurchildに設定する。そして、curparentの他のchildにchildのbrohterを設定し、ステップＳ５６に戻る。
【００９４】
ステップＳ５６の判断がＹｅｓの場合には、ステップＳ５９において、grandpにcurparentのparentを設定する。ステップＳ６０において、cutsurface＝grandp、nodeset＝parentとして、分割領域の構成要素ノードの番号付けを行なうサブルーチンを呼び出す。ステップＳ６１において、nodesetにgrandpのchildを設定し、grandpの他のchildに、そのchildのbrotherを設定する。ステップＳ６２において、設定されたcutsurface、nodesetで、分割領域の構成要素ノードの番号付けを行なうサブルーチンを呼び出す。
【００９５】
ステップＳ６３において、curparentにcurparentのparentを設定する。ステップＳ６４において、curparentのchildがnullか、すなわち、存在しないか否かを判断する。ステップＳ６４の判断がＮｏの場合には、ステップＳ６９において、curparentのchildをcurchildに設定し、ステップＳ５６に戻る。
【００９６】
ステップＳ６４の判断がＹｅｓの場合には、ステップＳ６５において、nodesetにcurparentを設定し、ステップＳ６６において、cutsurfaceの構成要素ノードを番号付けするサブルーチンを呼び出す。ステップＳ６７において、curparentにcurparentのparentを設定する。すなわち、領域分割を表すtreeをroot方向に戻る。ステップＳ６８において、curparentがnullか、すなわち、存在しないか否かを判断する。ステップＳ６８の判断がＮｏの場合には、ステップＳ６４に戻り、Ｙｅｓの場合には、処理を終了する。
【００９７】
図１９のサブルーチンは、図１７のステップＳ５７のサブルーチンであり、depth first searchで、領域ノードからなるtreeをleafの方向にたどるとき、最初に、このノードを通過したとき呼び出される。
【００９８】
そのとき、分割面の構成要素ノード取り出して、その隣接ノードが分割領域に属するとき、分割領域に対応する、leafノードごとにstacknode を持つ。各スタックにノードを積む位置は、番号付け処理の最初にクリアしておく（stacknode ごとにstackptr=0）。
【００９９】
図１９において、ステップＳ７０で、分割面を構成する構成要素ノードのチェインをたどって、構成要素ノードをstacksfに積む。ステップＳ７１において、stacksfから１つのノードを取り出す。ステップＳ７２において、現在の構成要素ノードの全ての隣接ノードを順番に取り出し、同じ分割領域に属しているか否かを調べる。対応する分割領域のstacknodeに積まれていなければ、その分割領域のstacknodeに積む。対応するstacknodeに積む位置は、対応するポインタを、stackptr＝stackptr＋１と更新して、積む位置を計算する。新たに構成要素ノードを積むときは、構成要素ノードの処理済flagをonにする。ステップＳ７３において、stacksfが空か否かを判断し、Ｙｅｓならばサブルーチンを抜け、Ｎｏならば、ステップＳ７１に戻る。
【０１００】
図２０は、図１８のステップＳ６０、Ｓ６２、Ｓ６６で呼び出されるサブルーチンの処理フローであり、stacknode に積まれ、breadth first searchのための始点ノード群をスタックの底からたどりながら、その隣接ノードで同じ分割領域に属するノードで、まだスキャンされていないノードをstacknode に積む処理を行う。
【０１０１】
図２０において、ステップＳ８０においては、stackptr2＝１とポインタを初期化する。ステップＳ８１において、stacknodeのstackptr2の指す位置からノードを取り出す。取り出したら、stackptr2＝stackptr2＋１とカウントアップする。ステップＳ８２において、取り出したノードの隣接ノードをスキャンして、領域番号が同じ分割領域に属し、まだスタックに積まれていないノードがあれば、stackptr＝stackptr＋１として、stacknodeに積む。ステップＳ８３において、stackptr2＞stackptrであるか否かを判断する。ステップＳ８３の判断がＮｏの場合には、ステップＳ８１に戻る。ステップＳ８３の判断がＹｅｓの場合には、ステップＳ８４において、counterに今までに番号付けした総ノード数を設定する。ステップＳ８５において、stacknodeの上の方から構成要素ノードを一つずつ取り出し、構成要素ノードにcounter＝counter＋１をノードの番号として、新番号の領域に格納する（図７参照）。stacknodeの上の方から取り出すとは、スタックに最後に積まれた構成要素ノードから取り出すという意味である。ステップＳ８６において、stacknodeが空か否かを判断し、ＮｏであればステップＳ８５に戻り、Ｙｅｓのであれば、サブルーチンを抜ける。
【０１０２】
図２１は、本実施形態が適用される共有メモリ型スカラ並列計算機のハードウェア構成を説明する図である。
マルチコアＣＰＵ１０は、１つのＣＰＵノードの内部にＣＰＵコア（Ｌ１キャッシュも内臓）１２を複数封入したもので、各コアからL2キャッシュおよびバスインタフェース１１を共通に使う。論理的には、各ＣＰＵコアが１つのＣＰＵとして動作するように見える。
【０１０３】
図２１の例では、１つのＣＰＵノードに２つのＣＰＵコアが封入された場合を示している。１つのＣＰＵノードに４つの以上のＣＰＵコアが封入されるものもある。ＣＰＵノードが一つで、8 個のＣＰＵコアがそのノードに封入されていてもSMP（Shared Memory Palallel computer)と見なせる。ここでＣＰＵと呼ぶのは、ＣＰＵコアのことである。ＣＰＵには各スレッドが割り当てられて処理が行われるが、スレッドは、ＯＳが並列処理を行う処理単位のことである。スパース行列の要素データのうち、処理に必要なデータは、メモリモジュール１４に分散されて保持され、相互結合網（バス）１３を介して、L2キャッシュに格納される。そして、L2キャッシュから、スーパーノードを構成する、スパース行列の構成要素ノードのデータがL1キャッシュにロードされ、ＣＰＵコアによって演算される。スーパーノードの幅が大きいほどＣＰＵコアの演算量に比べ、L1キャッシュへのデータのロードストア回数を減らすことができるので、高速に演算を行うことができる。
【０１０４】
上記実施形態のアルゴリズムは、プログラムとして実装することが出来る。プログラムは、記録媒体１６に格納され、記録媒体読み取り装置１５によって読み取られる。記録媒体読み取り装置１５によって読み取られたプログラムは、メモリモジュール１４のいずれかに格納される。記録媒体１６は、ＣＤ−ＲＯＭ、フレキシブルディスク、ＤＶＤ、Ｂｌｕ−Ｒａｙなどの持ち運び可能な媒体などである。プログラムの実行は、メモリモジュール１４に格納された状態で実行されても良いし、記録媒体１６から読み取りながら実行しても良い。
【０１０５】
以上の実施形態は、シミュレーションや数理モデルから生じるスパースな正値対称行列で表される連立1次方程式の解を求めることで、数理モデルの分析を行う分野、有限要素法、偏微分方程式などの解析を行いシミュレーションを行う分野に適用可能である。
【０１０６】
上記実施形態のorderingによる行列の演算と、スパース行列を、グラフの分割に基づいた一般化されたnested dissection により生成されたorderingで並べ変えた行列の演算と比べる。
【０１０７】
並列化されたコレスキー分解において、本実施形態のorderingによるコレスキー分解を行う性能は、比較対象の方法で並び替えた行列に対してコレスキー分解を行う性能に対し、7 倍程度高速なものもある。
【０１０８】
本実施形態のorderingと、もう分割できないまで完全に分割を繰り返して生成されたorderingである一般化されたnested dissection ordering(METIS V4.0で生成)と比較を行った場合以下のような結果であった。
【０１０９】
＜１．orderingの説明＞
３次元の立方体でラプラス方程式を有限差分で離散化して得られる方程式を各次元を60分割した。つまり60^3個の変数が現れる場合を用いた。１つのレベルのnested dissectionで、各次元を２分割して境界面を取り除くと、8つの小さな立方体が生成される。これを3レベル行った。生成された8^3個の小さな（ほぼ）立方体について、境界面から遠い順に番号付けを行うことでorderingを得る。
【０１１０】
＜２．データ＞
スレッド数 1 2 4 8 16 32 64 128
３レベル 74.667 40.874 28.432 15.606 10.955 7.42 6.407 6.133
ND 359.599 198.51 179.318 105.55 78.111 61.029 45.925 40.959
【０１１１】
＜３．データの解説＞
計算機は共用メモリ型スカラ並列計算機M9000 (128cores, SPARC64VII quad core CPU 2.52GHz)を使用し、この計算機上でスパースな正値対称行列の連立１次方程式を以下のorderingでスパース行列のLDL^T分解によって解いた性能を示している。スパース行列の次数は60^3であった。単位は秒(sec)である。
【０１１２】
＜orderingの方法＞
“3レベル”は、本実施形態で示したものを使用し、“ND”は、一般化されたnested dissection ordering（METIS 4.0でサポートされているもの）を利用した。
【０１１３】
16〜64スレッドでの並列実行で本実施形態は、7倍以上高速である。
【符号の説明】
【０１１４】
１０ ＣＰＵノード
１１ L2キャッシュ及びバスインタフェース
１２ L1キャッシュ及びＣＰＵコア
１３ 相互結合網（バス）
１４ メモリモジュール
【特許請求の範囲】
【請求項１】
共有メモリ型スカラ並列計算機によってスパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法であって、
共有メモリ型スカラ並列計算機は、
スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行い、
１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う
ことを特徴とするordering生成方法。
【請求項２】
再帰的分割を最後まで行なって得られた前記分割領域を構成する構成要素ノード群に対する番号付けは、生成された全ての前記分割面を構成する構成要素ノードから遠い構成要素ノードから順次番号付けすることを特徴とする請求項１に記載のordering生成方法。
【請求項３】
前記番号付けは、前記分割面を構成する構成要素ノードに接続された前記分割領域の構成要素ノードからの接続関係をたどりながら、構成要素ノードをスタックに積み上げ、全ての構成要素ノードを積み上げ終わった後、スタックの上から構成要素ノードを取り出して、順次番号付けすることによってなされることを特徴とすることを特徴とする請求項１または２に記載のordering生成方法。
【請求項４】
スパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法を共有メモリ型スカラ並列計算機に実現させるプログラムであって、
共有メモリ型スカラ並列計算機は、
スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行い、
１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う
ことを特徴とするプログラム。
【請求項５】
スパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法を実行する共有メモリ型スカラ並列計算機であって、
スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行う分割手段と、
１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う番号付け手段と、
を備えることを特徴とする共有メモリ型スカラ並列計算機。
【請求項１】
共有メモリ型スカラ並列計算機によってスパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法であって、
共有メモリ型スカラ並列計算機は、
スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行い、
１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う
ことを特徴とするordering生成方法。
【請求項２】
再帰的分割を最後まで行なって得られた前記分割領域を構成する構成要素ノード群に対する番号付けは、生成された全ての前記分割面を構成する構成要素ノードから遠い構成要素ノードから順次番号付けすることを特徴とする請求項１に記載のordering生成方法。
【請求項３】
前記番号付けは、前記分割面を構成する構成要素ノードに接続された前記分割領域の構成要素ノードからの接続関係をたどりながら、構成要素ノードをスタックに積み上げ、全ての構成要素ノードを積み上げ終わった後、スタックの上から構成要素ノードを取り出して、順次番号付けすることによってなされることを特徴とすることを特徴とする請求項１または２に記載のordering生成方法。
【請求項４】
スパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法を共有メモリ型スカラ並列計算機に実現させるプログラムであって、
共有メモリ型スカラ並列計算機は、
スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行い、
１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う
ことを特徴とするプログラム。
【請求項５】
スパースな正値対称行列をコレスキー分解あるいは修正コレスキー分解する際のordering生成方法を実行する共有メモリ型スカラ並列計算機であって、
スパース行列を表現したグラフを構成する構成要素ノードの集合を、２分割する分割面を取り除いて、２つの分割領域と分割面とに分割する処理を、該構成要素ノードのデータを結合してキャッシュに保持するために生成するスーパーノードの幅になるまで、再帰的に行う分割手段と、
１回の分割で構成要素ノードの集合が分割されて生成された分割領域の一方について分割面から遠い構成要素ノードから順次番号を付し、他方の分割領域についても分割面から遠い構成要素ノードから順次番号を付し、最後に分割面に属する構成要素ノードに番号付けする処理を再帰的に行う番号付け手段と、
を備えることを特徴とする共有メモリ型スカラ並列計算機。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【公開番号】特開２０１２−７３６８１（Ｐ２０１２−７３６８１Ａ）
【公開日】平成２４年４月１２日（２０１２．４．１２）
【国際特許分類】
【出願番号】特願２０１０−２１６１３６（Ｐ２０１０−２１６１３６）
【出願日】平成２２年９月２７日（２０１０．９．２７）
【出願人】（０００００５２２３）富士通株式会社 (25,993)
【Ｆターム（参考）】
【公開日】平成２４年４月１２日（２０１２．４．１２）
【国際特許分類】
【出願日】平成２２年９月２７日（２０１０．９．２７）
【出願人】（０００００５２２３）富士通株式会社 (25,993)
【Ｆターム（参考）】
[ Back to top ]