(57)【要約】【課題】ＳＭＰノード分散メモリ型並列計算機で高速に行列を処理することの出来る装置あるいは方法を提供する。【解決手段】ブロック化された行列のＬＵ分解法において、ネットワークで接続されたＳＭＰノードのそれぞれに、行列の内、更新べきブロックを縦にノード数分に分割し、それそれ分割された部分を各ノードに配置する。このような配置を更新すべきブロックについて、続けて行い、各ノードに、分割されたブロックの分割部分をサイクリックに割り当てるようにする。割り当てられた各ノードは、元のブロックの順番で分割されたブロックを更新する。元のブロックを順次更新することにより、各ノードの処理済みブロックの量が同じ分ずつ増えていき、負荷の分散をすることができる。【選択図】      図４
2004302928.tif 000002
【特許請求の範囲】
【請求項１】
複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理方法であって、
ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置ステップと、
該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離ステップと、
該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置ステップと、
該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解ステップと、
ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新ステップと、
を備えることを特徴とする行列処理方法を情報装置に実現させるプログラム。
【請求項２】
前記ＬＵ分解は、再帰的手続きにより、各ノードの各プロセッサで並列的に行われることを特徴とする請求項１に記載のプログラム。
【請求項３】
前記更新ステップにおいては、各ノードが、列ブロックを計算している間に、計算し終わった部分のデータであって、他のブロックの更新に必要なデータを該計算と平行して他のノードに転送することを特徴とする請求項１に記載のプログラム。
【請求項４】
前記並列計算機は、ＳＭＰ（ＳｙｍｍｅｔｒｉｃＭｕｌｔｉＰｒｏｃｅｓｓｏｒ）を各ノードとするＳＭＰノード分散メモリ型並列計算機であることを特徴とする請求項１に記載のプログラム。
【請求項５】
複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理装置であって、
ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置手段と、
該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離手段と、
該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置手段と、
該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解手段と、
ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新手段と、
を備えることを特徴とする行列処理装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、ＳＭＰ（Ｓｙｍｍｅｔｒｉｃ  ＭｕｌｔｉＰｒｏｓｅｓｓｏｒ）ノード分散メモリ型並列計算機における行列処理装置あるいは処理方法に関する。
【０００２】
【従来の技術】
ベクトルプロセッサをクロスバーで結合した並列計算機向けに開発した連立一次方程式の解法では、ブロックＬＵ分解の各ブロックを各ＰＥにサイクリックに配置してＬＵ分解を行っていた。ベクトルプロセッサではブロック幅を小さくしてもコストの高い行列積による更新部分の計算効率は非常に高かった。このためブロック幅１２程度でサイクリックな配置と見なして、まず、このブロックをＬＵ分解及び１つのＣＰＵで逐次的に計算してから、結果を部分的に分割して各プロセッサに転送して、行列積での更新を行っていた。
【０００３】
図２６は、スーパスカラ並列計算機用ＬＵ分解法のアルゴリズムを概略説明する図である。
配列Ａを外積形式のガウスの消去法をブロック化した方法でＬＵ分解する。ブロック幅ｄで分解する。
ｋ番目の処理で、更新部分Ａ（ｋ）を次の計算で更新する。
Ａ（ｋ）＝Ａ（ｋ）−Ｌ２（ｋ）×Ｕ２（ｋ）・・・・・（１）
ｋ＋１番目の処理では、Ａ（ｋ）を幅ｄで分解してｄだけ小さいマトリックスを同じ式で更新する。
Ｌ２（ｋ）、Ｕ２（ｋ）は以下の式で求める必要がある。
式（１）で更新を行う場合、
【０００４】
【数１】
2004302928.tif 000003
【０００５】
と分解し、Ｕ２（ｋ）＝Ｌ１（ｋ）−１Ｕ２（ｋ）と更新する。
上記のブロック化されたＬＵ分解の方法は、特許文献１に記載されている。
そのほか、並列計算機で行列を計算する技術として特許文献２には、連立１次方程式の係数行列を外部記憶装置に格納する方式が、特許文献３には、ベクトル計算機における方式が、特許文献４には、多枢軸同時消去を行う方式が、特許文献５には、スパース行列の各要素の構成を並び替えて、縁付きブロック対角行列にしてからＬＵ分解を行う方法が記載されている。
【０００６】
【特許文献１】
特開２００２−１６３２４６号公報
【特許文献２】
特開平９−１７９８５１号公報
【特許文献３】
特開平１１−６６０４１号公報
【特許文献４】
特開平５−２０３４９号公報
【特許文献５】
特開平３−２２９３６３号公報
【０００７】
【発明が解決しようとする課題】
上記スーパスカラ並列計算機用ＬＵ分解の方法を単純に一つのノードをＳＭＰとする並列計算機システムで行うと以下の問題が発生する。
【０００８】
ＳＭＰノードでの行列積を効率的に行うためにはベクトル計算機で１２と設定していたブロック幅を１０００程度に増やす必要がある。
（１）この結果、ブロック毎にそれが各プロセッサにサイクリックに配置されていると見なして処理を行うと、行列積での更新の計算量がプロセッサ間で不均一である割合が大きくなり並列効率が著しく低下する。
（２）また、１ノードで計算する幅１０００程度のブロックのＬＵ分解は、ノード内でのみ計算すると、他のノードはアイドル状態となる。幅の大きさに比例して、このアイドル時間が増えるため、並列化効率が著しく低下する。
（３）ＳＭＰノードを構成するＣＰＵ数を増やすと計算能力の増加に対して、転送スピードが相対的に劣化しているため、従来の方法は転送量が約０．５ｎ２×１．５要素（ここでの要素は、行列の要素である）であったが、相対的に増えて見える。このため効率がかなり落ちる。
（１）〜（３）までの劣化は全体で焼く２０〜２５％の性能ダウンを引き起こす。
【０００９】
本発明の課題は、ＳＭＰノード分散メモリ型並列計算機で高速に行列を処理することの出来る装置あるいは方法を提供することである。
【００１０】
【課題を解決するための手段】
本発明の行列処理方法は、複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理方法であって、ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置ステップと、該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離ステップと、該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置ステップと、該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解ステップと、ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新ステップとを備えることを特徴とする。
【００１１】
本発明によれば、各ノード間の計算負荷を分散し、並列化度を上げることが出来るので、より高速な行列処理が行える。また、演算と、データ転送を並列して行うことから、計算機の処理能力をデータ転送のスピードに制限されずに、向上することが出来る。
【００１２】
【発明の実施の形態】
本発明の実施形態においては、ブロック幅を大きくしても負荷バランスが完全に均一であり、１ＣＰＵで逐次計算していた部分をノード間で並列に処理する方式を提案する。
【００１３】
図１は、本発明の実施形態が適用されるＳＭＰノード分散メモリ型並列計算機の概略全体構成を示す図である。
図１（ａ）に示されるように、クロスバーネットワークにノード１〜ノードＮが接続され、相互に通信できるようになっている。各ノードは、図１（ｂ）に示されるように相互結合網１０によって、メモリモジュール１１−１〜１１−ｎ、及びプロセッサ１３−１〜１３−ｍとキャッシュ１２−１〜１２−ｍの組とが相互に結合され、通信可能となっている。データ通信用ハード（ＤＴＵ）１４は、図１（ａ）のクロスバーネットワークに接続され、他のノードと通信可能となっている。
【００１４】
まず、比較的ブロック幅の小さなコラムブロックをノードにサイクリックに配置する。各ノードに一巻き分（サイクリックにコラムブロックを配置した場合の１回でサイクリックに配置される分）あるブロックを一つに束ねたものを一つに行列と見なす。これは行列を２次元目を均等に分割し、各ノードに分散配置した状態と見なすことが出来る。これを１次元目を均等に分割した配置に並列転送を利用して動的に変更する。ここで、１次元目を分割、２次元目を分割とは、行列を長方形あるいは正方形とした場合、横方向を縦の線で分割することを１次元目を分割すると言い、縦方向を横の線で分割することを２次元目を分割するという。このとき一番上の正方形部分は各ノードが重複して持つようにする。
【００１５】
この分散配置の変更でクロスバーネットワークを利用した並列転送が使え、転送量はノード数分の１となる。１次元目を均等に分割した配置に変更した配置で、ノード間通信を使って、このブロックのＬＵ分解を並列に行う。このとき並列化効率があがり、かつＳＭＰの性能を引き出せるようにするために、更にブロックに分解して再帰的なＬＵ分解を行う。
【００１６】
このブロックＬＵ分解が終了した時点で各ノードには対角ブロック部分の情報と１次元目を均等に分割した部分の情報があるため、これを利用して行ブロック部分を更新して、保持している列ブロック部分とで更新できる部分を更新する。更新時に隣のノードにこの情報を転送して、次の更新の準備を行う。この転送は計算と同時に行える。これらの操作を繰り返して全ての更新部分の更新を行う。
【００１７】
図２は、本発明の実施形態に従った全体の処理フローチャートである。
まず、ステップＳ１０において、最後の一巻きか否かを判断する。ステップＳ１０の判断がＹＥＳの場合には、ステップＳ１５に進む。ステップＳ１０の判断がＮＯの場合には、ステップＳ１１において、対象となる一巻き分のブロックを結合したブロックを１次元目で分割した配置に並列転送を利用して変換する。このとき対角ブロックは全てのノードで共通に持つようにする。ステップＳ１２においては、１次元目を分割配置したブロックに関してＬＵ分解を行う。このときキャッシュの大きさを考慮したブロック幅までと、そのブロック幅より小さい部分の処理を再帰的な手続きで行う。ステップＳ１３では、ＬＵ分解した１次元目で分割配置されたブロックを並列転送をつかって元の２次元目を分割した配置に戻す。ステップＳ１４においては、この時点で各ノードには対角ブロックと残りをノード数に１次元目で分割した小ブロックが各ノードに割り付けられている。各ノードで共通に持っていた更新済みの対角ブロックを使ってブロック行を各ノードで更新する。このとき次の更新で必要となる列ブロックを隣のノードに計算と同時に転送する。ステップＳ１５では、最後の一巻きは各ノードに分割せずに冗長に配置して、同じ計算を行ってＬＵ分解を行う。各ノード部分に対応する部分をコピーバックする。そして、処理を終了する。
【００１８】
図３は、本発明の実施形態の一般概念図である。
図３に示されるように、行列を例えば、４等分して各ノードに分散配置する。各ノードは、列ブロックが割り当てられており、サイクリックな順序で処理する。このとき一巻き分を束ねて１つのブロックと見なす。これを対角ブロック部分を除き１次元目で分割し、通信を使って各ノードに再配置する。
【００１９】
図４及び図５は、比較的ブロック幅の小さなブロックをサイクリックに配置した状態を説明する図である。
図４及び図５に示すように、行列の一部の列ブロックを、更に小さい列ブロックに細分化し、各ノード（今の場合４つとしている）にサイクリックに割り当てる。このような配置の変更は、２次元目を分割されたブロックを１次元目を分割（対角ブロックは共通に保持）変更することになる。これはクロスバーネットワークの並列転送を利用して変更することが出来る。
【００２０】
これは、１巻きが結合されたブロックをメッシュに仮想的に分割したとき、対角線方向のブロックの並び（１１、２２、３３、４４）、（１２、２３、３４、４１）、（１３、２４、３１、４２）、（１４、２１、３２、４３）の各組のブロックを各ノードに（二次元目の示すプロセッサから１次元目の示すプロセッサに転送する）並列転送することで実現できる。このとき、対角ブロック部分も一緒に送ることで対角ブロック部分は各ノードが共通に持つことができる充分な大きさで、転送はプロセッサ数分の１になる。
【００２１】
このように分散配置を変更した列ブロックに対するＬＵ分解を、各ノードに対角ブロックと残りの部分を均等に分割したものを配置して、ノード間通信及びノード間で同期を取りながら処理を行う。また、ノード内でのＬＵ分解の処理はスレッド並列化を行う。
【００２２】
スレッド並列化でのＬＵ分解がキャッシュ上で効率的に行えるように、２重構造の再帰的手続で行う。つまり、あるブロック幅までの大きさで一次の再帰手続で行い、それより小さい部分に関しては、スレッド並列化のために、各スレッドで、そのブロックを対角部分と残りの部分を並列処理するスレッド数で均等に分割した部分を合わせて連続な作業域にコピーして処理を行う。このことでキャッシュ上のデータを有効に利用する。
【００２３】
また、ノード間で共有している対角ブロック部分の計算はノード間で冗長に計算されてノード間のＬＵ分解の並列化効率が劣化する。ＬＵ分解を２重の再帰的手続きで行うことで、各ノード内でスレッドで並列計算するときのオーバヘッドを減らすことが出来る。
【００２４】
図６は、図４及び図５で配置されたブロックの更新処理を説明する図である。図６の最も左のブロックは各ノードに対角ブロックを冗長に、かつ、残りのブロックを一次元目で均等に分割したものを作業域に配置したものである。あるノードでの状態と考える。最小ブロック幅まで１次の再帰手続きを行う。
【００２５】
最小ブロックのＬＵ分解が終わったら、この情報を使って、行ブロック及び更新部分の更新を更新する領域を均等に分割して、並列に更新する。
最小ブロック部分のＬＵ分解は、更に以下のように最小幅のブロックの対角部分を共通に、かつ、残り部分を均等に分割して、各スレッドの局所領域（キャッシュの大きさ程度）にコピーする。
【００２６】
この領域を使って、更に再帰的手続きでＬＵ分解を行う。ピボットを決めて、行の入れ替えを行うために各スレッドに、ピボットの相対的位置から、ノードでの相対位置、全体での位置に換算するための情報を保持しておく。
【００２７】
ピボットがスレッドの局所領域の対角部分内にあるときは、各スレッドで、独立に入れ替えを行える。
スレッドの対角ブロックを超えたときは、その位置が、以下の条件のときによって処理が異なる。
ａ）ピボットがノード間に分割配置したとき冗長に配置した対角ブロック内にあるとき。
【００２８】
このときは、ノード間で通信する必要はなく、各ノードで独立に処理できる。
ｂ）ピボットがノード間に分割配置した時冗長に配置した対角ブロックを超えたとき。
【００２９】
このときはスレッド間での最大値、つまりノードでの最大値を全ノードに通信して最大ピボットがどのノードに有るかを決定する。これが決まった後、最大ピボットを持つノードで行の入れ替えを行う。そのあと、入れ替えられた行（ピボット行）を他のノードに通信する。
【００３０】
このようなピボットの処理を行う。
２重構造を持つ再帰手続きでのＬＵ分解の二次のスレッド並列で行うＬＵ分解は、上記のピボット処理を行いながら、各スレッドの局所領域でＬＵ分解を並列に行うことができる。
【００３１】
ピボットの入れ替えの履歴は共用メモリに各ノードに冗長に保持する。
図７は、再帰的なＬＵ分解の手順を説明する図である。
再帰的なＬＵ分解の手順は以下のようになる。
【００３２】
図７（ｂ）のレイアウトを考える。図７（ｂ）の対角ブロック部分がＬＵ分解できると、ＵはＬ１を使って、Ｕ←Ｌ１−１Ｕ、Ｃ←Ｌ×Ｕと更新する。
再帰的手順は、ＬＵ分解する領域を前半と後半に分割し、分割した領域をＬＵ分解の対象と見なして、再帰的に行う方法である。ブロックの幅が、ある最小の幅より小さくなったとき、その幅に関しては従来通りのＬＵ分解を行う。
【００３３】
図６（ａ）は、領域を真ん中の太線で２分割し、その左側をＬＵ分解する過程で更に２分割したところである。太線で分割した左側は図６（ｂ）のレイアウトを当てはめられる。このレイアウトのＣの部分もＬＵ分解できたとき、太線から左側のＬＵ分解が終わる。
【００３４】
この左側の情報から、図６（ｂ）のレイアウトを全体にあてはめて、Ｃとなる右側の更新を行う。更新が終わったら、右側に図６（ｂ）のレイアウトを当てはめて同じようにＬＵ分解を行う。
・ブロックのＬＵ分解処理の後の行の入れ替えと行ブロックの更新及びｒａｎｋ  ｐ  ｕｐｄａｔｅでの更新
ノード間にブロックを再配置した状態でノード間通信及びスレッド並列を使ってＬＵ分解を並列に実行した後、各ノードには各ノードに共通に置かれた対角ブロックと残りの部分を均等に分割した部分のひとかけらがＬＵ分解された値を保持して残る。
【００３５】
各ノードでピボットの入れ替えの履歴の情報と対角ブロックの情報を使って、まず行の入れ替えを行う。その後、行ブロック部分の更新を行う。この後、対角ブロックの残り部分を分割した列ブロック部分と更新された行ブロック部分を利用して更新部分を更新する。この計算と同時に更新に使う分割された列ブロック部分を全ノードで隣のノードに転送する。
【００３６】
この転送は、次の更新で必要な情報を計算と同時に送り、次の計算の前までに準備を行うためであり、転送を計算と同時に行うことで計算を効率よく続けることができる。
【００３７】
また、部分的な行列積の更新をスレッド数が多くても効率的に行えるように各スレッドで計算する行列積の更新領域が正方形に近くなる用に分割する。各ノードで更新を受け持つ更新領域は、正方形である。この領域の更新を各スレッドに分担して、かつ、性能劣化を引き起こさないようにすることを考える。
【００３８】
このため、更新領域をできるだけ正方形に近い形に分割する。このことで更新部分の２次元目の大きさがかなり大きく取れ、行列積の計算で繰り返し参照される部分の参照をキャッシュ上に保持して有効利用することが比較的できるようになる。
【００３９】
このために、以下の手順で行列積の更新の各スレッドでの分担を決めて並列計算する。
１）スレッドの総数＃ＴＨＲＤの平方根を求める。
２）この値が整数でないとき、これを切り上げてｎｒｏｗとする。
３）２次元目の分割数をｎｒｏｗとする。
４）１次元目の分割数をｎｃｏｌを以下の条件を満たす最小の整数を見つける。
ｎｃｏｌ×ｎｒｏｗ＞＝＃ＴＨＲＤ
５）ｉｆ（ｎｃｏｌ＊ｎｒｏｗ＝＝＃ｔｈｒｄ）ｔｈｅｎ
１次元目をｎｃｏｌ等分、２次元目をｎｒｏｗ等分ｎｃｏｌ＊ｎｒｏｗに分割して各スレッドに更新を並列実行させる。
ｅｌｓｅ
１次元目をｎｃｏｌ等分、２次元目をｎｒｏｗ等分してｎｃｏｌ＊ｎｒｏｗに分割して（１、１）、（１、２）、（１、３）、・・・（２、１）、（２、２）、（２、３）・・・と＃ＴＨＲＤ個の部分を並列更新する。残りの領域は一般的に横に長い長方形となる。これを２次元目を均等に分割して全スレッドで負荷が均等になるように更新部分を分割して再度並列処理する。
ｅｎｄｉｆ
・ソルバー部分
図８は、対角部分以外の部分ブロックの更新について説明する図である。
【００４０】
ＬＵ分解された結果は、各ノードに分散配置された形で保存されている。各ノードには比較的ブロック幅の小さなブロックがＬＵ分解された状況で格納されている。
【００４１】
この幅の小さなブロックに関して前進代入、後退代入を行って次のブロックのある隣のノードに処理を渡す。このとき解を更新した部分を隣のノードに転送する。
【００４２】
実際の前進代入及び後退代入では細長いブロックで対角ブロック部分を除いた長方形部分を１次元目で均等にスレッド数で分割して並列更新を行う。
まず、一つスレッドでＬＤ×ＢＤ＝ＢＤを解く。
【００４３】
この情報を使って全スレッドで以下のようにＢを並列に更新する。
Ｂｉ＝Ｂｉ−Ｌｉ×ＢＤ
この１サイクルの更新で変更された部分を隣のノードに転送する。
【００４４】
前進代入が終わったら、今までの処理でノードに処理を渡してきたのとちょうど逆を辿るようにして後退代入を行う。
実際には、元の行列の各ノードに配置された部分をサイクリックに処理している。これは列ブロックを入れ替えて別の行列に変換していることに相当する。ＬＵ分解の過程でピボットをとる列は未分解部分のどの列を対象にしてもよいことに由来する。
ＡＰＰ−１ｘ＝ｂ→ｙ＝Ｐ−１ｘと置いてｙについて解くことに相当する。解いたｙを並び変えることでｘを求めることが出来る。
【００４５】
図９〜図１１は、行ブロックの更新処理を説明する図である。
列ブロックの計算が終わったら、今度計算された部分をもとの２次元目を分割した配置に戻す。ここで、２次元目を分割した形でのデータは各ノードに保持しておく。次に、行の入れ替え情報を元に、行の入れ替えを行ったあと、行ブロックを更新する。
【００４６】
各ノードに存在する列ブロックの部分を計算と同時に隣のノードにリング状に送ることで順次更新を進めていく。バッファをもう一つ持つことで可能となる。この領域には各ノードに対角ブロックを冗長に保持しているが、これも一緒に転送する。対角ブロック以外の部分のデータの量が多く、また、計算と同時に転送を行うので、転送時間は見えない。
【００４７】
図１０によれば、バッファＡからＢへのデータ転送を行う。次のタイミングではバッファＢからＡへのノードのリングに沿ってデータを送る。このようにしてスイッチしてデータ送る。更に、図１１において、更新が終わったら、列ブロックと行ブロックを除いた正方行列に対して大きさが縮小したもの対して同じ処理を繰り返す。
【００４８】
図１２〜図２５は、本発明の実施形態のフローチャートである。
図１２及び図１３はサブルーチンｐＬＵのフローである。このサブルーチンは、呼び出しプログラムであり、各ノードで１つのプロセスを生成してから呼び出すことで並列に処理を行う。
【００４９】
まず、解くべき問題の大きさを、単位ブロック数をｉｂｌｋｓｕｎｉｔ、ノード数をｎｕｍｎｏｒｄとして、ｎ＝ｉｂｌｋｓｕｎｉｔ×ｎｕｍｎｏｒｄ×ｍ（ｍは各ノードでの単位ブロック数）としたＬＵ分解を行う。各ノードに係数行列Ａの２次元目を均等に分割した共用メモリＡ（ｋ、ｎ／ｎｕｍｎｏｒｄ）（ｋ＞＝ｎ）及び行の入れ替えの履歴を格納するｉｐ（ｎ）を引数として受け取る。ステップＳ２０において、ｎｏｎｏｒｄにプロセス番号（１〜ノード数）を設定し、ｎｕｍｎｏｒｄにノード数（全プロセス数）を設定する。ステップＳ２１において、各ノードでスレッドを生成し、ｎｏｔｈｒｄにスレッド番号（１〜スレッド数）及びｎｕｍｔｈｒｄにスレッドの総数を設定する。ステップＳ２２において、ブロック幅の設定であるｉｂｌｋｓｍａｃｒｏ＝ｉｂｌｋｓｕｎｉｔ×ｎｕｍｎｏｒｄ、繰り返し回数であるｌｏｏｐ＝ｎ／（ｉｂｌｋｓｕｎｉｔ×ｎｕｍｔｈｒｄ）−１を計算し、更に、ｉ＝１、ｌｅｎｂｕｆｍａｘ＝（ｎ−ｉｂｌｋｓｍａｃｒｏ）／ｎｕｍｎｏｒｄ＋ｉｂｌｋｓｍａｃｒｏを設定する。
【００５０】
ステップＳ２３において、ｗｌｕ１（ｌｅｎｂｕｆｍａｘ，  ｉｂｌｋｓｍａｃｒｏ）、ｗｌｕ２（ｌｅｎｂｕｆｍａｘ，  ｉｂｌｋｓｍａｃｒｏ）、ｂｕｆｓ（ｌｅｎｂｕｆｍａｘ，  ｉｂｌｋｓｕｎｉｔ）、ｂｕｆｄ（ｌｅｎｂｕｆｍａｘ，  ｉｂｌｋｓｕｎｉｔ）の作業域を確保する。この領域をサブルーチンが実行の都度、実際の長さｌｅｎｂｕｆを計算して、必要な大きさだけ使う。
【００５１】
ステップＳ２４においては、ｉ＞＝ｌｏｏｐであるか否かを判断する。ステップＳ２４の判断がＹＥＳの場合には、ステップＳ３７に進む。ステップＳ２４の判断がＮＯの場合には、ステップＳ２５において、ノード間でバリア同期を取る。そして、ステップＳ２６において、ｌｅｎｂｌｋｓ＝（ｎ−ｉ×ｉｂｌｋｓｍａｃｒｏ）／ｎｕｍｎｏｒｄ＋ｉｂｌｋｓｍａｃｒｏを計算する。ステップＳ２７において、サブルーチンｃｔｏｂを呼び出し、各ノードにある幅ｉｂｌｋｓｕｎｉｔのｉ番目を対角ブロックと１次元目を均等分割した幅ｉｂｌｋｓｍａｃｒｏの部録を対角ブロックに結合し、ノードに持つ配置を変える。ステップＳ２８では、ノード間でバリア同期を取る。ステップＳ２９では、サブルーチンｉｎｔｅｒｌｕを呼び出して、配列ｗｌｕ１に格納され、分散再配置された、ブロックをＬＵ分解する。行の入れ替えの情報は、ｉｓ＝（ｉ−１）＊ｉｂｌｋｓｍａｃｒｏ＋１，ｉｅ＝ｉ＊ｉｂｌｋｓｍａｃｒｏとしてｉｐ（ｉｓ：ｉｅ）に格納されている。
【００５２】
ステップＳ３０において、ノード間でバリア同期を取り、ステップＳ３１において、サブルーチンｂｔｏｃを呼び出して、再配置されたブロックでＬＵ分解されたブロックを各ノードのもともと格納されていた場所に戻す。ステップＳ３２においてノード間でバリア同期を取り、ステップＳ３３において、サブルーチンｅｘｒｗを呼び出して、行の入れ替え及び行ブロックの更新を行う。ステップＳ３４においては、ノード間でバリア同期を取り、ステップＳ３５において、サブルーチンｍｍｃｂｔを呼び出して、各ノードにある列ブロックの部分（ｗｌｕ１に格納されている）と行ブロックの部部との行列積で更新する。計算と同時に列ブロック部分をプロセッサ間をリングに沿って転送し、次の更新の準備を行いながら更新する。ステップＳ３６においては、ｉ＝ｉ＋１として、ステップＳ２４に戻る。
【００５３】
ステップＳ３７では、ノード間でバリア同期を取り、ステップＳ３８において、生成したスレッドを消滅させる。ステップＳ３９において、サブルーチンｆｂｌｕを呼んで、最後のブロックのＬＵ分解を行いながら更新する。ステップＳ４０において、ノード間でバリア同期を取り、処理を終了する。
【００５４】
図１４及び図１５は、サブルーチンｃｔｏｂのフローである。
ステップＳ４５において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）、ｂｕｆｓ（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｕｎｉｔ）、ｂｕｆｄ（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｕｎｉｔ）を引数で受けて、各ノードのｉ番目の幅ｉｂｌｋｓｕｎｉｔのブロックをｎｕｍｎｏｒｄ個束ねたものの対角ブロック行列部分より下の部分をｎｕｍｎｏｒｄ個に分割したものと対角ブロックを加えたものとを各ノードに分散配置したものに転送を利用して配置換えする。
【００５５】
ステップＳ４６においては、ｎｂａｓｅ＝（ｉ−１）＊ｉｂｌｋｓｍａｃｒｏ（ｉは呼び出し元のメインループの繰り返し回数）、ｉｂｓ＝ｎｂａｓｅ＋１、ｉｂｅ＝ｎｂａｓｅ＋ｉｂｌｋｓｍａｃｒｏ、ｌｅｎ＝（ｎ−ｉｂｅ）／ｎｕｍｎｏｒｄ、ｎｂａｓｅ２ｄ＝（ｉ−１）＊ｉｂｌｋｓｕｎｉｔ、ｉｂｓ２ｄ＝ｎｂａｓｅ２ｄ＋１、ｉｂｅ２ｄ＝ｉｂｓ２ｄ＋ｉｂｌｋｓｕｎｉｔを計算する。ここで、送信データ数はｌｅｎｓｅｎｄ＝（ｌｅｎ＋ｉｂｌｋｓｍａｃｒｏ）＊ｉｂｌｋｓｕｎｉｔである。ステップＳ４７においては、ｉｙ＝１と設定し、ステップＳ４８において、ｉｙ＞ｎｕｍｎｏｒｄか否かを判断する。ステップＳ４８の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ４８の判断がＮＯの場合には、ステップＳ４９において、送信する部分、受信する部分を決める。すなわち、ｉｄｓｔ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｉｙ−１，ｎｕｍｎｏｒｄ）＋１（送信先ノード番号）、ｉｓｒｓ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｎｕｍｎｏｒｄ−ｉｙ＋１，ｎｕｍｎｏｒｄ）＋１（送信元ノード番号）を計算する。ステップＳ５０においては、各ノードで自分に割り付いている幅ｉｂｌｋｓｕｎｉｔの対角ブロック部分と、その下部分のブロックの１次元目をｎｕｍｎｏｒｄで分割した部分で、再配置した時保持する部分（転送先のノード数番目のもの）をバッファの下の部分に格納する。すなわち、ｂｕｆｄ（１：ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｕｎｉｔ）←Ａ（ｉｂｓ：ｉｂｅ，ｉｂｓ２ｄ：ｉｂｅ２ｄ）、ｉｃｐｓ＝ｉｂｅ＋（ｉｄｓｔ−１）＋ｌｅｎ＋１、ｉｃｐｅ＝ｉｓｐｓ＋ｌｅｎ−１、ｂｕｆｄ（ｉｂｌｋｓｍａｃｒｏ＋１：ｌｅｎ＋ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｕｎｉｔ）←Ａ（ｉｃｐｓ：ｉｃｐｅ，ｉｂｓ２ｄ：ｉｂｅ２ｄ）を演算する。このコピーは１次元目をスレッド数に分割して各スレッドで並列に処理する。
【００５６】
ステップＳ５１では、全ノードで送受信を行う。すなわち、ｂｕｆｄの内容おｉｄｓｔ番目のノードに送り、ｂｕｆｓに受信する。ステップＳ５２においては、送受信の完了を待つ。ステップＳ５３では、バリア同期を取り、ステップＳ５４において、ｗｌｕ１の対応位置に、ｉｓｒｓ番目のノードから受けたデータを格納する。すなわち、
ｉｃｐ２ｄｓ＝（ｉｓｒｓ−１）＊ｉｂｌｋｓｕｎｉｔ＋１，ｉｃｐ２ｄｅ＝ｉｃｐ２ｄｓ＋ｉｂｌｋｓｕｎｉｔ−１、ｗｌｕ１（１：ｌｅｎ＋ｉｂｌｋｓｍａｃｒ，，ｉｃｐ２ｄｓ：
ｉｃｐ２ｄｅ）←ｂｕｆｓ（１：ｌｅｎ＋ｉｂｌｋｓｕｎｉｔ，１：ｂｌｋｓｕｎｉｔ）を演算する。すなわち、１次元目をスレッド数で分割して各スレッドで並列コピーする。ステップＳ５５でｉｙ＝ｉｙ＋１とし、ステップＳ４８に戻る。
【００５７】
図１６及び図１７は、サブルーチンｉｎｔｅｒＬＵのフローである。
ステップＳ６０において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）、ｗｌｕｍｉｃｒｏ（ｎｃａｓｈ）を引数として受ける。ここで、ｗｌｕｍｉｃｒｏをＬ２キャッシュ（レベル２のキャッシュ）の大きさとし、各スレッドに確保されたものを受ける。ｗｌｕ１にＬＵブロック分解する幅ｉｂｌｋｓｍａｃｒｏのブロックで対角ブロックとその下位ブロックを１次元目でｎｕｍｎｏｒｄ個に分割した一つが各ノードの領域に格納されている。ピボットの検索と行の入れ替えに関してノード間転送を使いながら並列にＬＵ分解する。本サブルーチンは、再帰的に呼び出される。呼び出しが深くなるにつれてＬＵ分解したときのブロック幅は小さくなる。このブロックをスレッド並列してＬＵ分解したとき、各スレッドで計算する部分がキャッシュの大きさ以下になるところで、ＬＵ分解をスレッド並列化する別のサブルーチンを呼び出す。
【００５８】
スレッド並列は対象となる比較的幅の小さなブロックをこのブロックの対角行列部分を各スレッドで重複して持ち、対角ブロックより下位の部分を１次元目をスレッド数で均等分割して各スレッド（ＣＰＵ）にキャッシュの大きさより小さな領域ｗｌｕｍｉｃｒｏで処理できるようにコピーして処理を行う。ｉｓｔｍｉｃｒｏは小さなブロックの先頭位置であり、最初１に設定される。ｎｗｉｄｔｈｍｉｃｒｏは、小さなブロックの幅であり、最初は全体のブロック幅に設定される。ｉｂｌｋｓｍｉｃｒｏｍａｘは、小さなブロックの最大値であり、これ以上大きいときブロック幅を更に小さく（例えば、８０列に）する。ｎｏｔｈｒｄはスレッド番号、ｎｕｍｔｈｒｄはスレッド数、各ノードで重複して持つ１次元配列ｉｐ（ｎ）に行の入れ替え情報を入れる。
【００５９】
ステップＳ６１では、ｎｗｉｄｔｈｍｉｃｒｏ＜＝ｉｂｌｋｓｍｉｃｒｏｍａｘであるか否かを判断する。ステップＳ６１の判断がＹＥＳの場合には、ステップＳ６１において、ｉｂｌｋｓｍｉｃｒｏ＝ｎｗｉｄｔｈｍｉｃｒｏとし、各ノードに分担した領域にある対角ブロックと分割したブロックが格納されているｗｌｕ（ｌｅｎｍａｃｒｏ，ｉｂｌｋｓｍａｃｒｏ）のｗｌｕ（ｉｓｔｍｉｃｒｏ：ｌｅｎｍａｃｒｏ，ｉｓｔｍｉｃｒｏ：ｉｂｌｋｓｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏ−１）の部分に関して対角部分ｗｌｕ（ｉｓｔｍｉｃｒｏ：ｉｓｔｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏ−１，ｉｓｔｍｉｃｒｏ：ｉｓｔｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏ−１）を対角ブロックとする。また、ｉｒｅｓｔ＝ｉｓｔｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏとし、ｗｌｕ（ｉｒｅｓｔ：ｌｅｎｍａｃｒｏ，ｉｓｔｍｉｃｒｏ：ｉｓｔｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏ−１）を１次元目でスレッド数で均等分割したものを対角ブロックと結合して、各スレッド毎の領域ｗｌｕｍｉｃｒｏにコピーする。すなわち、ｌｅｎｍｉｃｒｏ＝（ｌｅｎｍａｒｏ−ｉｒｅｓｔ＋ｎｕｍｔｈｒｄ）／ｎｕｍｔｈｒｄとし、ｗｌｕｍｉｃｒｏ（ｌｅｎｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏ，ｉｂｌｋｓｍｉｃｒｏ）にコピーし、ｌｅｎｂｌｋｓｍｉｃｒｏ＝ｌｅｎｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏとする。そして、ステップＳ６３で、サブルーチンＬＵｍｉｃｒｏを呼び出す。これにおいては、ｗｌｕｍｉｃｒｏ（ｌｉｎｍｉｃｒｏ＋ｉｂｌｋｓｍｉｃｒｏ，ｉｂｌｋｓｍｉｃｒｏ）を受け渡す。ステップＳ６４では、ｗｌｕｍｉｃｒｏに分割していた部分を、対角部分は１つのスレッドから、他の部分は各スレッドのｗｌｕｍｉｃｒｏからｗｌｕに元々あった部分に戻す。そして、サブルーチンを抜ける。
【００６０】
ステップＳ６１の判断がＮＯの場合には、ステップＳ６５において、ｎｗｉｄｔｈｍｉｃｒｏ＞＝３＊ｉｂｌｋｓｍｉｃｒｏｍａｘまたは、ｎｗｉｄｔｈｍｉｃｒｏ＜＝２＊ｉｂｌｋｓｍｉｃｒｏｍａｘか否かを判断する。ステップＳ６５の判断がＹＥＳの場合には、ステップＳ６６において、ｎｗｉｄｔｈｍｉｃｒｏ２＝ｎｗｉｄｔｈｍｉｃｒｏ／２、ｉｓｔｍｉｃｒｏ２＝ｉｓｔｍｉｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ２、ｎｗｉｄｔｈｍｉｃｒｏ３＝ｎｗｉｄｔｈｍｉｃｒｏ−ｎｗｉｄｔｈｍｉｃｒｏ２とし、ステップＳ６８に進む。ステップＳ６５の判断がＮＯの場合には、ステップＳ６７において、ｎｗｉｄｔｈｍｉｃｒｏ２＝ｎｗｉｄｔｈｍｉｃｒｏ／３，  ｉｓｔｍｉｃｒｏ２＝ｉｓｔｍｉｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ２，  ｎｗｉｄｔｈｍｉｃｒｏ３＝ｎｗｉｄｔｈｍｉｃｒｏ−ｎｗｉｄｔｈｍｉｃｒｏ２とし、ステップＳ６８に進む。ステップＳ６８では、ｉｓｔｉｍｉｃｒｏは、そのまま、ｎｗｉｄｔｈｍｉｃｒｏとしてｎｗｉｄｔｈｍｉｃｒｏ２を渡してサブルーチンｉｎｔｅｒＬＵを呼び出す。
【００６１】
ステップＳ６９においては、ｗｌｕ（ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ−１）の部分を更新する。これは、一つのスレッドで更新すれば充分である。これにｗｌｕ（ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ２−１，ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ２−１）の下三角行列の逆行列を左から乗じたもので更新する。ステップＳ７０においては、ｗｌｕ（ｉｓｔｍｉｃｒｏ２：ｌｅｎｍａｃｒｏ，ｉｓｔｍｉｃｒｏ２：ｉｓｔｍｉｃｒｏ２＋ｎｗｉｄｔｈｍｉｃｒｏ３−１）をｗｌｕ（ｉｓｔｍｉｃｒｏ２：ｌｅｎｍａｃｒｏ，ｉｓｔｍｉｃｒｏ：ｉｓｔｍｉｃｒｏ２−１）×ｗｌｕ（ｉｓｔｍｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ２−１，ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ２：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ−１）を引いて更新する。このとき、１次元目をスレッド数で均等に分割して並列計算する。ステップＳ７１においては、ｉｓｔｍｉｃｒｏとして、ｉｓｔｍｉｃｒｏ２、ｎｗｉｄｔｈｍｉｃｒｏとしてｎｗｉｄｔｈｍｉｃｒｏ３を渡してサブルーチンｉｎｔｅｒＬＵを呼び出し、サブルーチンを終了する。
【００６２】
図１８及び図１９は、サブルーチンＬＵｍｉｃｒｏのフローである。
ステップＳ７５において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）、ｗｌｕｍｉｃｒｏ（ｌｅｎｉｂｌｋｓｍｉｃｒｏ，ｉｂｌｋｓｍｉｃｒｏ）を引数として受ける。ここで、ｗｌｕｍｉｃｒｏをＬ２キャッシュの大きさの各スレッドに確保されたものを受ける。本ルーチンでｗｌｕｍｉｃｒｏに格納された部分のＬＵ分解を行う。ｉｓｔは、ＬＵ分解するブロックの先頭位置で最初は、１とされる。ｎｗｉｄｔｈは、ブロック幅であり、最初は全体のブロック幅である。ｉｂｌｋｓｍａｘは、ブロック最大値（８程度）であり、これ以上小さくしない。ｗｌｕｍｉｃｒｏはスレッド毎に引数として渡される。
【００６３】
ステップＳ７６においては、ｎｗｉｄｔｈ＜＝ｉｂｌｋｓｍａｘか否かを判断する。ステップＳ７６の判断がＮＯの時は、ステップＳ８８に進む。ステップＳ７６の判断がＹＥＳの場合には、ステップＳ７７において、ｉ＝ｉｓｔとして、ステップＳ７８において、ｉ＜ｉｓｔ＋ｎｗｉｄｔｈか否かを判断する。ステップＳ７８の判断がＮＯの場合には、サブルーチンを抜ける。ステップＳ７８の判断がＹＥＳの場合には、ステップＳ７９において、各スレッドでｉ列目の絶対値最大の要素を見つけ、共用メモリ領域にスレッド番号順に格納する。ステップＳ８０においては、各ノードでのノード内の最大ピボットをこの中から見つけ、この要素とノード番号、位置をセットとして全ノードが各セットを持つように通信し、各ノードで全ノードでの最大ピボットを決定する。なお、各ノードで同じ方法で最大ピボットを決定する。
【００６４】
ステップＳ８１においては、このピボット位置が各ノードが持つ対角ブロックの中か判定する。ステップＳ８１の判断がＮＯの場合には、ステップＳ８５に進む。ステップＳ８１の判断がＹＥＳの場合には、ステップＳ８２において、最大ピボットの位置が各スレッドが重複して持つ対角ブロックの中かを判定する、ステップＳ８２の判断がＹＥＳの場合には、ステップＳ８３において、全ノードで保持する対角ブロック内での入れ替えで、かつ、全スレッドで重複して持つ対角部分内での入れ替えなので、スレッドで独立してピボットの入れ替えを行う。入れ替えた位置を配列ｉｐに格納し、ステップＳ８６に進む。ステップＳ８２における判断がＮＯの場合には、ステップＳ８４において、各ノードで独立にピボットとを交換する。交換すべきピボット行を共用域に格納して、各スレッドの対角ブロック部分と入れ替える。入れ替えた位置を配列ｉｐに格納し、ステップＳ８６に進む。
【００６５】
ステップＳ８５では、ノード間で通信して最大ピボットを有するノードから交換すべき行ベクトルをコピーする。その後ピボット行を入れ替える。ステップＳ８６においては、行を更新し、ステップＳ８７において、ｉ列と行で更新部分を更新し、ｉ＝ｉ＋１として、ステップ７８に戻る。
【００６６】
ステップＳ８８においては、ｎｗｉｄｔｈ＞＝３＊ｉｂｌｋｓｍａｘあるいは、ｎｗｉｄｔｈ＜＝２＊ｉｂｌｋｓｍａｘであるか否かを判断する。ステップＳ８８の判断がＹＥＳの場合には、ステップＳ８９において、ｎｗｉｄｔｈ＝ｎｗｉｄｔｈ／２、ｉｓｔ２＝ｉｓｔ＋ｎｗｉｄｔｈ２とし、ステップＳ９１に進む。ステップＳ８８の判断がＮＯの場合には、ステップＳ９０において、ｎｗｉｄｔｈ２＝ｎｗｉｄｔｈ／３、ｉｓｔ２＝ｉｓｔ＋ｎｗｉｄｔｈ２、ｎｗｉｄｔｈ３＝ｎｗｉｄｔｈ−ｎｗｉｄｔｈ２とし、ステップＳ９１に進む。ステップＳ９１では、ｉｓｔはそのまま、ｎｗｉｄｔｈとしてｎｗｉｄｔｈ２を引数として渡して、サブルーチンＬＵｍｉｃｒｏを呼び出す。ステップＳ９２では、ｗｌｕｍｉｃｒｏ（ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈ２−１，ｉｓｔｍｉｃｒｏ＋ｎｗｉｄｔｈ２：ｉｓｔｍｉｃｒｏ＋ｎｗｉｄｔｈｍｉｃｒｏ−１）の部分を更新する。ｗｌｕｍｉｃｒｏ（ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈ２−１，ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈ２−１）の下三角行列の逆行列を左から乗したもので更新する。ステップＳ９３では、ｗｌｕｍｉｃｒｏ（ｉｓｔｍｉｃｒｏ２：ｌｅｎｍａｃｒｏ，ｉｓｔｍｉｃｒｏ２：ｉｓｔｍｉｃｒｏ２＋ｎｗｉｄｔｈｍｉｃｒｏ３−１）をｗｌｕｍｉｃｒｏ（ｉｓｔｍｉｃｒｏ２：ｌｅｎｍａｃｒｏ，ｉｓｔｍｉｃｒｏ：ｉｓｔｍｉｃｒｏ２−１）×ｗｌｕｍｉｃｒｏ（ｉｓｔｍｉｃｒｏ：ｉｓｔｍａｃｒｏ＋ｎｗｉｄｔｈ２−１，ｉｓｔ＋ｎｗｉｄｔｈ２：ｉｓｔ＋ｎｗｉｄｔｈｍｉｃｒｏ−１）を引いて更新する。ステップＳ９４においてはｉｓｔとしてｉｓｔ２、ｎｗｉｄｔｈとしてｎｗｉｄｔｈ３を引数として受け渡して、サブルーチンＬＵｍｉｃｒｏを呼び出して、サブルーチンを抜ける。
【００６７】
図２０は、サブルーチンｂｔｏｃのフローである。
ステップＳ１００において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）、ｂｕｆｓ（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｕｎｉｔ）、ｂｕｆｄ（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｕｎｉｔ）を引数で受けて、各ノードのｉ番目の幅ｉｂｌｋｓｕｎｉｔのブロックをｎｕｍｎｏｒｄ個束ねたものの対角ブロック行列部分ｉｂｌｋｓｍａｃｒｏ×ｉｂｌｋｓｍａｃｒｏより下の部分をｎｕｍｎｏｒｄ個に分割したものと対角ブロックを加えたものを各ノードに分散配置したものに転送を利用して配置を変える。
【００６８】
ステップＳ１０１では、ｎｂａｓｅ＝（ｉ−１）＊ｉｂｌｋｓｍａｃｒｏ（ｉは呼び出しもとのメインループの繰り返し回数）、ｉｂｓ＝ｎｂａｓｅ＋１、ｉｂｅ＝ｎｂａｓｅ＋ｉｂｌｋｓｍａｃｒｏ、ｌｅｎ＝（ｎ−ｉｂｅ）／ｎｕｍｎｏｒｄ、ｎｂａｓｅ２ｄ＝（ｉ−１）＊ｉｂｌｋｓｕｎｉｔ、ｉｂｓ２ｄ＝ｎｂａｓｅ２ｄ＋１、ｉｂｅ２ｄ＝ｉｂｓ２ｄ＋ｉｂｌｋｓｕｎｉｔとし、送信データ数は、ｌｅｎｓｅｎｄ＝（ｌｅｎ＋ｉｂｌｋｓｍａｃｒｏ）＊ｉｂｌｋｓｕｎｉｔとする。
【００６９】
ステップＳ１０２において、ｉｙ＝１とし、ステップＳ１０３において、ｉｙ＞ｎｕｍｎｏｒｄか否かを判断する。ステップＳ１０３の判断がＹＥＳの場合、サブルーチンを抜ける。ステップＳ１０３の判断がＮＯの場合には、ステップＳ１０４において、送信する部分、受信する部分を決める。すなわち、ｉｄｓｔ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｉｙ−１，ｎｕｍｎｏｒｄ）＋１、ｉｓｒｓ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｎｕｍｎｏｒｄ−ｉｙ＋１，ｎｕｍｎｏｒｄ）＋１とする。ステップＳ１０５においては、計算結果が格納されているｗｌｕ１から元の位置に配置を戻すための送信のためにバッファに格納する。ｉｄｓｔ番目のノードに対応部分を送る。すなわち、ｉｃｐ２ｄｓ＝（ｉｄｓｔ−１）＊ｉｂｌｋｓｕｎｉｔ＋１、ｉｃｐ２ｄｅ＝ｉｃｐ２ｄｓ＋ｉｂｌｋｓｕｎｉｔ−１、ｂｕｆｄ（１：ｌｅｎ＋ｉｂｌｋｓｕｎｉｔ，１：ｉｂｌｋｓｕｎｉｔ）←ｗｌｕ１（１：ｌｅｎ＋ｉｂｌｋｓｍａｃｒｏ，ｉｃｐ２ｄｓ：ｉｃｐ２ｄｅ）とする。１次元目をスレッド数で分割して各スレッドで並列コピーする。
【００７０】
ステップＳ１０６では、全ノードで送受信する。ｂｕｆｄの内容をｉｄｓｔ番目のノードに送り、ｂｕｆｓに受信する。ステップＳ１０７で送受信の完了を待ち、ステップＳ１０８において、バリア同期を取る。ステップＳ１０９では、各ノードで自分に割り付いている幅ｉｂｌｋｓｕｎｉｔの対角ブロック部分と、その下の部分のブロックの１次元目をｎｕｍｎｏｒｄで分割した部分で再配置したときの部分（転送先のノード数番目のもの）を元々あった部分に格納する。Ａ（ｉｂｓ：ｉｂｅ，ｉｂｓ２ｄ：ｉｂｄ２ｄ）←ｂｕｆｓ（１：ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｕｎｉｔ）、ｉｃｐｓ＝ｉｂｅ＋（ｉｓｒｓ−１）＊ｌｅｎ＋１、ｉｃｐｅ＝ｉｓｐｓ＋ｌｅｎ−１、Ａ（ｉｃｐｓ：ｉｃｐｅ，ｉｂｓ２ｄ：ｉｂｅ２ｄ）←ｂｕｆｓ（ｉｂｌｋｓｍａｃｒｏ＋１：ｌｅｎ＋ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｕｎｉｔ）とする。このコピーは１次元目をスレッド数に分割して各スレッドで列毎に処理する。
【００７１】
ステップＳ１１０においては、ｉｙ＝ｉｙ＋１として、ステップＳ１０３に戻る。
図２１は、サブルーチンｅｘｒｗのフローである。
このサブルーチンは、行の入れ替え及び行ブロックの更新を行うものである。
【００７２】
ステップＳ１１５においては、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ、ｉｂｌｋｓｍａｃｒｏ）を引数として受ける。ｗｌｕ１（１：ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｍａｃｒｏ）には、ＬＵ分解された対角部分を全ノードが重複して持っている。ｎｂｄｉａｇ＝（ｉ−１）＊ｉｂｌｋｓｍａｃｒｏとする。ｉは呼び出し元のサブルーチンｐＬＵのメインループの繰り返し回数である。また、ピボットの入れ替えの情報が、ｉｐ（ｎｂｄｉａｇ＋１：ｎｂｄｉａｇ＋ｉｂｌｋｓｍａｃｒｏ）に格納されている。
【００７３】
ステップＳ１１６では、ｎｂａｓｅ＝ｉ＊ｉｂｌｋｓｕｎｉｔ（ｉは呼び出しもとのサブルーチンｐＬＵのメインループの繰り返し回数）、ｉｒｏｗｓ＝ｎｂａｓｅ＋１、ｉｒｏｗｅ＝ｎ／ｎｕｍｎｏｒｄ、ｌｅｎ＝（ｉｒｏｗｅ−ｉｒｏｗｓ＋１）／ｎｕｍｔｈｒｄ、ｉｓ＝ｎｂａｓｅ＋（ｎｏｔｈｒｄ−１）＊ｌｅｎ＋１、ｉｅ＝ｍｉｎ（ｉｒｏｗｅ，ｉｓ＋ｌｅｎ−１）とする。ステップＳ１１７では、ｉｘ＝ｉｓとする。
【００７４】
ステップＳ１１８では、ｉｓ＜＝ｉｅであるか否かを判断する。ステップＳ１１８の判断がＮＯの場合には、ステップＳ１２５に進む。ステップＳ１１８の判断がＹＥＳの場合には、ステップＳ１１９において、ｎｂｄｉａｇ＝（ｉ−１）＊ｉｂｌｋｓｍａｃｒｏ、ｊ＝ｎｂｄａｇ＋１として、ステップＳ１２０において、ｊ＜＝ｎｂｄｉａｇ＋ｉｂｌｋｓｍａｃｒｏであるか否かを判断する。ステップＳ１２０の判断がＮＯの場合には、ステップＳ１２４に進む。ステップＳ１２０の判断がＹＥＳの場合には、ステップＳ１２１において、ｉｐ（ｊ）＞ｊか否かを判断する。ステップＳ１２１の判断がＮＯの場合には、ステップＳ１２３に進む。ステップＳ１２１の判断がＹＥＳの場合には、ステップＳ１２２において、Ａ（ｊ、ｉｘ）とＡ（ｉｐ（ｊ），ｉｘ）を入れ替えて、ステップＳ１２３に進む。ステップＳ１２３においては、ｊ＝ｊ＋１として、ステップＳ１２０に戻る。
【００７５】
ステップＳ１２４においては、ｉｘ＝ｉｘ＋１とし、ステップＳ１１８に戻る。
ステップＳ１２５においては、バリア同期（全ノード、全スレッド）を取る。
ステップＳ１２６においては、Ａ（ｎｂｄｉａｇ＋１：ｎｂｄｉａｇ＋ｉｂｌｋｓｍａｃｒｏ，ｉｓ：ｉｅ）←ＴＲＬ（ｗｌｕ１（ｉ：ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｍａｃｒｏ））−１×Ａ（ｎｂｄｉａｇ＋１：ｎｂｄｉａｇ＋ｉｂｌｋｓｍａｃｒｏ，ｉｓ：ｉｅ）を全ノード、全スレッドで更新する。ここで、ＴＲＬ（Ｂ）は、行列Ｂの下三角部分を示す。ステップＳ１２７では、バリア同期（全ノード、全スレッド）を取って、サブルーチンを抜ける。
【００７６】
図２２及び図２３は、サブルーチンｍｍｃｂｔのフローである。
ステップＳ１３０において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ、ｉｂｌｋｓｍａｃｒｏ）、ｗｌｕ２（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）を引数として受ける。ｗｌｕ１に、ブロック幅ｉｂｌｋｓｍａｃｒｏのブロックをＬＵ分解した結果で、対角ブロックとその下位ブロックを１次元目でｎｕｍｎｏｒｄ個に分割した一つが格納されている。分割した順にノード番号に対応し、ノードに再配置される。これをノードのリングに沿って転送しながら（計算と同時に行う）行列積を行いながら更新する。計算の裏で性能に影響を与えないので計算に直接使用しない対角ブロック部分も一緒に送る。
【００７７】
ステップＳ１３１では、ｎｂａｓｅ＝（ｉ−１）＊ｉｂｌｋｓｍａｃｒｏ（ｉは呼び出しもとのサブルーチンｐＬＵのメインループの繰り返し回数）、ｉｂｓ＝ｎｂａｓｅ＋１、ｉｂｅ＝ｎｂａｓｅ＋ｉｂｌｋｓｍａｃｒｏ、ｌｅｎ＝（ｎ−ｉｂｅ）／ｎｕｍｎｏｒｄ、ｎｂａｓｅ２ｄ＝（ｉ−１）＊ｉｂｌｋｓｕｎｉｔ、ｉｂｓ２ｄ＝ｎｂａｓｅ２ｄ＋１、ｉｂｅ２ｄ＝ｉｂｓ２ｄ＋ｉｂｌｋｓｕｎｉｔ、ｎ２ｄ＝ｎ／ｎｕｍｎｏｒｄ、ｌｅｎｓｅｎｄ＝ｌｅｎ＋ｉｂｌｋｓｍａｃｒｏとし、送信データ数は、ｎｗｌｅｎ＝ｌｅｎｓｅｎｄ＊ｉｂｌｋｓｍａｃｒｏとする。
【００７８】
ステップＳ１３２において、ｉｙ＝１（初期値を設定）、ｉｄｓｔ＝ｍｏｄ（ｎｏｎｏｒｄ，ｎｕｍｎｏｒｄ）＋１（送り先ノード番号（隣ノード））、ｉｓｒｓ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｎｕｍｎｏｒｄ−１，ｎｕｍｎｏｒｄ）＋１（発信元ノード番号）、ｉｂｐ＝ｉｄｓｔとする。
【００７９】
ステップＳ１３３において、ｉｙ＞ｎｕｍｎｏｒｄであるか否かを判断する。ステップＳ１３３の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ１３３の判断がＮＯの場合には、ステップＳ１３４において、ｉｙ＝１か否かを判断する。ステップＳ１３４の判断がＹＥＳの場合には、ステップＳ１３６に進む。ステップＳ１３４の判断がＮＯの場合には、ステップＳ１３５において送受信の官僚を待つ。ステップＳ１３６では、ｉｙ＝ｎｕｍｎｏｒｄ（奇数の最後）であるか否かを判断する。ステップＳ１３６の判断がＹＥＳの場合には、ステップＳ１３８に進む。ステップＳ１３６の判断がＮＯの場合には、ステップＳ１３７において、送受信を行う。ｗｌｕ１の内容を（対角ブロックも含めて）隣のノード（ノード番号ｉｄｓｔ）に送る。かつ、ｗｌｕ２に（ノード番号ｉｓｒｓから）送られてくるデータを格納する。送受信データ長はｎｗｌｅｎとする。
【００８０】
ステップＳ１３８において、ｗｌｕ１のデータを使った更新のポジションを計算する。ｉｂｐ＝ｍｏｄ（ｉｂｐ−１＋ｎｕｍｎｏｒｄ−１，ｎｕｍｎｏｒｄ）＋１、ｎｃｐｔｒ＝ｎｂｅ＋（ｉｂｐ−１）＊ｌｅｎ＋１（１次元目の開始位置）とする。ステップＳ１３９では、行列積を計算するサブルーチンｐｍｍを呼び出す。このときｗｌｕ１を引き渡す。ステップＳ１４０において、ｉｙ＝ｎｕｍｎｏｒｄ（最後の処理が終わった）か否かを判断する。ステップＳ１４０の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ１４０の判断がＮＯの場合には、ステップＳ１４１において、行列積演算と同時に行っている送受信の完了を待つ。ステップＳ１４２において、ｉｙ＝ｎｕｍｎｏｒｄ−１（偶数の最後）であるか否かを判断する。ステップＳ１４２の判断がＮＯの場合には、ステップＳ１４４に進む。ステップＳ１４２の判断がＮＯの場合には、ステップＳ１４３において、送受信を行う。すなわち、ｗｌｕ２の内容を（対角ブロックも含めて）隣のノード（ノード番号ｉｄｓｔ）に送る。かつ、ｗｌｕ１に（ノード番号ｉｓｒｓから）送られてくるデータを格納する。送受信データ長はｎｗｌｅｎとする。
【００８１】
ステップＳ１４４では、ｗｌｕ２のデータを使った更新のポジションを計算する。すなわち、ｉｂｐ＝ｍｏｄ（ｉｂｐ−１＋ｎｕｍｎｏｒｄ−１，ｎｕｍｎｏｒｄ）＋１、ｎｃｐｔｒ＝ｎｂｅ＋（ｉｂｐ−１）＊ｌｅｎ＋１（１次元目の開始位置）とする。
【００８２】
ステップＳ１４５では、行列積を計算するサブルーチンｐｍｍを呼び出す。このとき、ｗｌｕ２を引き渡す。ステップＳ１４６において、ｉｙ＝ｉｙ＋２と、２を加えて、ステップＳ１３３に戻る。
【００８３】
図２４は、サブルーチンｐｍｍのフローである。
ステップＳ１５０において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｌｅｎｂｌｋｓ、ｉｂｌｋｓｍａｃｒｏ）、もしくは、ｗｌｕ２（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）をｗｌｕｘ（ｌｅｎｂｌｋｓ，ｉｂｌｋｓｍａｃｒｏ）に受ける。呼び出し元から渡された１次元目の開始位置ｎｃｐｔｒを使って正方形の領域を更新する。ｉｓ２ｄ＝ｉ＊ｉｂｌｋｓｕｎｉｔ＋１、ｉｅ２ｄ＝ｎ／ｎｕｍｎｏｒｄ、ｌｅｎ＝ｉｅ２ｄ−ｉｓ２ｄ＋１、ｉｓｌｄ＝ｎｃｐｔｒ、ｉｅｌｄ＝ｎｐｔｒ＋ｌｅｎ−１（ｉはサブルーチンｐＬＵの繰り返し数）、Ａ（ｉｓｌｄ：ｉｅｌｄ，ｉｓ２ｄ：ｉｅ２ｄ）＝Ａ（ｉｓｌｄ：ｉｅｌｄ，ｉｓ２ｄ：ｉｅ２ｄ）−ｗｌｕ（ｉｂｌｋｓｍａｃｒｏ＋１：ｉｂｌｋｓｍａｃｒｏ＋ｌｅｎ，１：ｉｂｌｋｓｍａｃｒｏ）×Ａ（ｉｓｌｄ−ｉｂｌｋｓｍａｃｒｏ：ｉｓｌｄ−１，ｉｓ２ｄ：ｉｅ２ｄ）（式１）とする。
【００８４】
ステップＳ１５１において、並列に処理するスレッド数の平方根を求めて切り上げる。ｎｕｍｒｏｏｔ＝ｉｎｔ（ｓｑｒｔ（ｎｕｍｔｈｒｄ））、もし、ｓｑｒｔ（ｎｕｍｔｈｒｄ）−ｎｕｍｒｏｏｔが０でないなら、ｎｕｍｒｏｏｔ＝ｎｕｍｒｏｏｔ＋１とする。ここで、ｉｎｔは小数点以下切り捨て、ｓｑｒｔは、平方根である。ステップＳ１５２において、ｍ１＝ｎｕｍｒｏｏｔ、ｍ２＝ｎｕｍｒｏｏｔ、ｍｘ＝ｍ１とする。ステップＳ１５３において、ｍ１＝ｍｘ、ｍｘ＝ｍｘ−１、ｍｍ＝ｍｘ×ｍ２とする。ステップＳ１５４において、ｍｍ＜ｎｕｍｔｈｒｄであるか否かを判断する。ステップＳ１５４の判断がＮＯの場合には、ステップＳ１５３に戻る。ステップＳ１５４の判断がＹＥＳの場合には、ステップＳ１５５において、更新する領域を１次元目をｍ１等分する。２次元目をｍ２等分して、ｍ１×ｍ２個の矩形にする。そのうち、ｎｕｍｔｈｒｄ個を各スレッドに割り当てて、（式１）の対応部分を並列に計算する。（１，１）、（１，２）、・・・（１，ｍ２）、（２，１）・・・・と２次元目の方向に順番にスレッドを対応付けていく。
【００８５】
ステップＳ１５６において、ｍ１＊ｍ２−ｎｕｍｔｈｒｄ＞０か否かを判断する。ステップＳ１５６の判断がＹＥＳの場合には、ステップＳ１５８に進む。ステップＳ１５６の判断がＮＯの場合には、ステップＳ１５７において、残りの矩形は最後の矩形の最後の行、１行の最後からｍ１＊ｍ２−ｎｕｍｔｈｒｄ個が更新されずに残っている。この矩形を結合して１つの矩形と考え、２次元目をスレッド数ｎｕｍｔｈｒｄで分割して（式１）の対応部分を並列に計算する。そして、ステップＳ１５８において、バリア同期（スレッド間）をとって、サブルーチンを抜ける。
【００８６】
図２５は、サブルーチンｆｂｌｕのフローである。
ステップＳ１６０において、Ａ（ｋ、ｎ／ｎｕｍｎｏｒｄ）、ｗｌｕ１（ｉｂｌｋｓｍａｃｒｏ、ｉｂｌｋｓｍａｃｒｏ）、ｂｕｆｓ（ｉｂｌｋｓｍａｃｒｏ、ｉｂｌｋｓｕｎｉｔ）、ｂｕｆｄ（ｉｂｌｋｓｍａｃｒｏ、ｉｂｌｋｓｕｎｉｔ）を引数で受けて、各ノードの幅ｉｂｌｋｓｕｎｉｔの最後のブロックをｎｕｍｎｏｒｄ個束ねたものを各ノードで重複して持つように利用不足部分を各ノードに送る。各ノードがｉｂｌｋｓｍａｃｒｏ×ｉｂｌｋｓｍａｃｒｏのブロックを重複して持った後、各ノードで同じ行列に対してＬＵ分解を行う。ＬＵ分解が完了したら、各ノードに配置されていた部分をコピーバックする。
【００８７】
ステップＳ１６１では、ｎｂａｓｅ＝ｎ−ｉｂｌｋｓｍａｃｒｏ、ｉｂｓ＝ｎｂａｓｅ＋１、ｉｂｅ＝ｎ、ｌｅｎ＝ｉｂｌｋｓｍａｃｒｏ、ｎｂａｓｅ２ｄ＝（ｉ−１）＊ｉｂｌｋｓｕｎｉｔ、ｉｂｓ２ｄ＝ｎ／ｎｕｍｎｏｒｄ−ｉｂｌｋｓｕｎｉｔ＋１、ｉｂｅ２ｄ＝ｎ／ｎｕｍｎｏｒｄとし、送信データ数はｌｅｎｓｅｎｄ＝ｉｂｌｋｓｍａｃｒｏ＊ｉｂｌｋｓｕｎｉｔとし、ｉｙ＝１とする。
【００８８】
ステップＳ１６２においては、バッファへのコピーを行う。すなわち、ｂｕｆｄ（１：ｉｂｌｋｓｍａｃｒｏ，１：ｉｂｌｋｓｕｎｉｔ）←Ａ（ｉｂｓ：ｉｂｅ，ｉｂｓ２ｄ：ｉｂｅ２ｄ）とする。ステップＳ１６３においては、ｉｙ＞ｎｕｍｎｏｒｄか否かを判断する。ステップＳ１６３の判断がＹＥＳの場合には、ステップＳ１７０に進む。ステップＳ１６３の判断がＮＯの場合には、ステップＳ１６４において、送信する部分、受信する部分を決定する。すなわち、ｉｄｓｔ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｉｙ−１，ｎｕｍｎｏｒｄ）＋１、ｉｓｒｓ＝ｍｏｄ（ｎｏｎｏｒｄ−１＋ｎｕｍｎｏｒｄ−ｉｙ＋１，ｎｕｍｎｏｒｄ）＋１とする。ステップＳ１６５では、全ノードで送受信する。ｂｕｆｄの内容をｉｄｓｔ番目のノードに送る。ステップＳ１６６においては、ｂｕｆｓにデータを受信し、送受信の完了を待つ。ステップＳ１６７において、バリア同期を取り、ステップＳ１６８において、ｗｌｕ１の対応位置にｉｓｒｓ番目のノードから来たデータを格納する。ｉｃｐ２ｄｓ＝（ｉｓｒｓ−１）＊ｉｂｌｋｓｕｎｉｔ＋１、ｉｃｐ２ｄｅ＝ｉｃｐ２ｄｓ＋ｉｂｌｋｓｕｎｉｔ−１、ｗｌｕ（１：ｉｂｌｋｓｍａｃｒｏ，ｉｃｐ２ｄｓ：ｉｃｐ２ｄｅ）←ｂｕｆｓ（１：ｉｂｌｋｓｕｎｉｔ，１：ｉｂｌｋｓｕｎｉｔ）とする。ステップＳ１６９において、ｉｙ＝ｉｙ＋１とし、ステップＳ１６３に戻る。
【００８９】
ステップＳ１７０では、バリア同期をとり、ステップＳ１７１では、ｗｌｕ１の上でｉｂｌｋｓｍａｃｒｏ×ｉｂｌｋｓｍａｃｒｏのＬＵ分解を各ノードで重複して行う。行交換の情報は、ｉｐに格納する。ＬＵ分解が終了したら、自ノード分を最後のブロックにコピーバックする。すなわち、ｉｓ＝（ｎｏｎｏｒｄ−１）＊ｉｂｌｋｓｕｎｉｔ＋１、ｉｅ＝ｉｓ＋ｉｂｌｋｓｕｎｉｔ−１、Ａ（ｉｂｓ：ｉｂｅ，ｉｂｓ２ｄ：ｉｂｅ２ｄ）←ｗｌｕ１（１：ｉｂｌｋｓｍａｃｒｏ，ｉｓ：ｉｅ）として、サブルーチンを抜ける。
（付記１）複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理方法であって、
ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置ステップと、
該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離ステップと、
該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置ステップと、
該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解ステップと、
ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新ステップと、
を備えることを特徴とする行列処理方法を情報装置に実現させるプログラム。
【００９０】
（付記２）前記ＬＵ分解は、再帰的手続きにより、各ノードの各プロセッサで並列的に行われることを特徴とする付記１に記載のプログラム。
（付記３）前記更新ステップにおいては、各ノードが、列ブロックを計算している間に、計算し終わった部分のデータであって、他のブロックの更新に必要なデータを該計算と平行して他のノードに転送することを特徴とする付記１に記載のプログラム。
【００９１】
（付記４）前記並列計算機は、ＳＭＰ（ＳｙｍｍｅｔｒｉｃＭｕｌｔｉＰｒｏｃｅｓｓｏｒ）を各ノードとするＳＭＰノード分散メモリ型並列計算機であることを特徴とする付記１に記載のプログラム。
【００９２】
（付記５）複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理装置であって、
ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置手段と、
該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離手段と、
該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置手段と、
該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解手段と、
ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新手段と、
を備えることを特徴とする行列処理装置。
【００９３】
（付記６）複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理方法であって、
ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置ステップと、
該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離ステップと、
該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置ステップと、
該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解ステップと、
ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新ステップと、
を備えることを特徴とする行列処理方法。
【００９４】
（付記７）複数のプロセッサとメモリを含む複数のノードをネットワークで接続した並列計算機における行列処理方法であって、
ノード毎にサイクリックに割り付けられた行列の部分の列ブロックの１巻き分を、該１巻き分をまとめたものを対象にして処理するために、各ノードに一つずつ分散して配置する第１の配置ステップと、
該１巻き分を結合したブロックに対して対角部分と該対角ブロックの下側にある列ブロックと他のブロックに分離する分離ステップと、
該対角ブロックを各ノードに冗長に配置すると共に、該列ブロックを１次元目で分割することによって得られるブロックを該複数のノードに、共に並列通信して一つずつ配置する第２の配置ステップと、
該対角ブロックと配置されたブロックを、各ノード間で通信しながら、各ノードで並列にＬＵ分解するＬＵ分解ステップと、
ＬＵ分解されたブロックを用いて、行列の他のブロックを更新する更新ステップと、
を備えることを特徴とする行列処理方法を情報装置に実現させるプログラムを格納する、情報装置読み取り可能な記録媒体。
【００９５】
【発明の効果】
ブロックを動的に１次元目の分割にして処理し、分解した後の各ノードの情報を使って更新し、転送は計算と同時に行える。このため更新部分は負荷はノード間で完全に均等になり、転送量はノード数分の１に削減できる。
【００９６】
ブロック幅を大きくすると負荷のバランスが崩れる従来の方法に対し負荷が均等になるため並列化効率が１０％程度向上する。また、転送量が減ることで３％程度の並列化率の向上に寄与でき、転送スピードがＳＭＰノードの計算性能に比べて遅くなっても影響は受けにくい。
【００９７】
ブロック部分のＬＵ分解をノード間で並列計算することによって、ブロック幅を大きくしたとき並列化出来ない部分の割合が増加するため並列化効率が落ちる部分をキャンセルできて約１０％の並列化効率の向上が見込める。また、ブロックＬＵ分解を、ミクロなブロックをベースにした再帰的プログラミングを使うことで、対角ブロックも含めてＳＭＰの並列化ができてＳＭＰでの並列処理での性能劣化を抑えることができる。
【図面の簡単な説明】
【図１】本発明の実施形態が適用されるＳＭＰノード分散メモリ型並列計算機の概略全体構成を示す図である。
【図２】本発明の実施形態に従った全体の処理フローチャートである。
【図３】本発明の実施形態の一般概念図である。
【図４】比較的ブロック幅の小さなブロックをサイクリックに配置した状態を説明する図（その１）である。
【図５】比較的ブロック幅の小さなブロックをサイクリックに配置した状態を説明する図（その２）である。
【図６】図４及び図５で配置されたブロックの更新処理を説明する図である。
【図７】再帰的なＬＵ分解の手順を説明する図である。
【図８】対角部分以外の部分ブロックの更新について説明する図である。
【図９】行ブロックの更新処理を説明する図（その１）である。
【図１０】行ブロックの更新処理を説明する図（その２）である。
【図１１】行ブロックの更新処理を説明する図（その３）である。
【図１２】本発明の実施形態のフローチャート（その１）である。
【図１３】本発明の実施形態のフローチャート（その２）である。
【図１４】本発明の実施形態のフローチャート（その３）である。
【図１５】本発明の実施形態のフローチャート（その４）である。
【図１６】本発明の実施形態のフローチャート（その５）である。
【図１７】本発明の実施形態のフローチャート（その６）である。
【図１８】本発明の実施形態のフローチャート（その７）である。
【図１９】本発明の実施形態のフローチャート（その８）である。
【図２０】本発明の実施形態のフローチャート（その９）である。
【図２１】本発明の実施形態のフローチャート（その１０）である。
【図２２】本発明の実施形態のフローチャート（その１１）である。
【図２３】本発明の実施形態のフローチャート（その１２）である。
【図２４】本発明の実施形態のフローチャート（その１３）である。
【図２５】本発明の実施形態のフローチャート（その１４）である。
【図２６】スーパスカラ並列計算機用ＬＵ分解法のアルゴリズムを概略説明する図である。
【符号の説明】
１０      相互結合網（バス）
１１−１〜１１−ｎ      メモリモジュール
１２−１〜１２−ｍ      キャッシュ
１３−１〜１３−ｍ      プロセッサ
１４      データ通信用ハード（ＤＴＵ）
【図１】
2004302928.tif 000004
【図２】
2004302928.tif 000005
【図３】
2004302928.tif 000006
【図４】
2004302928.tif 000007
【図５】
2004302928.tif 000008
【図６】
2004302928.tif 000009
【図７】
2004302928.tif 000010
【図８】
2004302928.tif 000011
【図９】
2004302928.tif 000012
【図１０】
2004302928.tif 000013
【図１１】
2004302928.tif 000014
【図１２】
2004302928.tif 000015
【図１３】
2004302928.tif 000016
【図１４】
2004302928.tif 000017
【図１５】
2004302928.tif 000018
【図１６】
2004302928.tif 000019
【図１７】
2004302928.tif 000020
【図１８】
2004302928.tif 000021
【図１９】
2004302928.tif 000022
【図２０】
2004302928.tif 000023
【図２１】
2004302928.tif 000024
【図２２】
2004302928.tif 000025
【図２３】
2004302928.tif 000026
【図２４】
2004302928.tif 000027
【図２５】
2004302928.tif 000028
【図２６】
2004302928.tif 000029