(57)【特許請求の範囲】
【請求項１】
非ゼロ要素の配置が対角要素を挟んで対称である行列を、下三角行列と上三角行列とを含む複数の行列に分解する計算機において、
前記行列の分解結果を記憶する記憶手段と、
連続する対角要素を共有する列群と行群とを分解の処理単位として、非ゼロ要素の配置に基づいて前記行列を複数の処理単位に分け、前記複数の処理単位それぞれについて、含まれる列及び行分のデータ量に、所定数の列及び行分のデータ量を加算した合計データ量に相当する記憶容量の第１記憶領域を、前記記憶手段に確保する第１記憶領域確保手段と、
前記複数の処理単位それぞれを所定の順番で対象処理単位とし、前記対象処理単位の列及び行と、対角要素の値が所定値以下であるために分解処理実施済みの処理単位から前記対象処理単位に移動された列及び行とに対して、前記第１記憶領域を使用して分解処理を実施する第１分解手段と、
前記対象処理単位に移動された列及び行、並びに前記対象処理単位の分解処理の結果、対角要素の値が所定値以下であるために移動することになった列及び行のうち、前記第１記憶領域に入りきらない列及び行がある場合、前記第１記憶領域に入りきらない列及び行のデータ量に相当する記憶容量の第２記憶領域を、前記記憶手段に確保する第２記憶領域確保手段と、
前記第１記憶領域に入りきらない列及び行に対して、前記第２記憶領域を使用して分解処理を実施する第２分解手段と、
を有する計算機。
【請求項２】
前記第１分解手段と前記第２分解手段とは、移動対象がｉ（ｉは１以上の整数）番目の行とｉ番目の列の場合、ｉ列内の行番号がｉ以上の要素を、前記対象処理単位に属する列の次に移動し、ｉ行内の列番号がｉより大きい要素を、前記対象処理単位に属する行の次に移動する、
請求項１記載の計算機。
【請求項３】
所定の規則に従って、前記複数の処理単位間の親子関係を決定する親子関係決定手段をさらに有し、
前記第１分解手段は、分解処理実施済みの処理単位のうち、前記対象処理単位の子に該当する処理単位の分解処理において、対角要素の値が所定値以下であるために処理できなかった列及び行を、前記対象処理単位に属する列及び行の次に移動して、分解処理を実施する、
請求項１または２記載の計算機。
【請求項４】
非ゼロ要素の配置が対角要素を挟んで対称である行列を、下三角行列と上三角行列とを含む複数の行列に分解する行列分解方法において、
コンピュータが、
連続する対角要素を共有する列群と行群とを分解の処理単位として、非ゼロ要素の配置に基づいて前記行列を複数の処理単位に分け、前記複数の処理単位それぞれについて、含まれる列及び行分のデータ量に、所定数の列及び行分のデータ量を加算した合計データ量に相当する記憶容量の第１記憶領域を、記憶手段に確保し、
前記複数の処理単位それぞれを所定の順番で対象処理単位とし、前記対象処理単位の列及び行と、対角要素の値が所定値以下であるために分解処理実施済みの処理単位から前記対象処理単位に移動された列及び行とに対して、前記第１記憶領域を使用して分解処理を実施し、
前記対象処理単位に移動された列及び行、並びに前記対象処理単位の分解処理の結果、対角要素の値が所定値以下であるために移動することになった列及び行のうち、前記第１記憶領域に入りきらない列及び行がある場合、前記第１記憶領域に入りきらない列及び行のデータ量に相当する記憶容量の第２記憶領域を、前記記憶手段に確保し、
前記第１記憶領域に入りきらない列及び行に対して、前記第２記憶領域を使用して分解処理を実施する、
行列分解方法。
【請求項５】
非ゼロ要素の配置が対角要素を挟んで対称である行列を、下三角行列と上三角行列とを含む複数の行列に分解する処理をコンピュータに実行させる行列分解プログラムにおいて、
前記コンピュータに、
連続する対角要素を共有する列群と行群とを分解の処理単位として、非ゼロ要素の配置に基づいて前記行列を複数の処理単位に分け、前記複数の処理単位それぞれについて、含まれる列及び行分のデータ量に、所定数の列及び行分のデータ量を加算した合計データ量に相当する記憶容量の第１記憶領域を、記憶手段に確保し、
前記複数の処理単位それぞれを所定の順番で対象処理単位とし、前記対象処理単位の列及び行と、対角要素の値が所定値以下であるために分解処理実施済みの処理単位から前記対象処理単位に移動された列及び行とに対して、前記第１記憶領域を使用して分解処理を実施し、
前記対象処理単位に移動された列及び行、並びに前記対象処理単位の分解処理の結果、対角要素の値が所定値以下であるために移動することになった列及び行のうち、前記第１記憶領域に入りきらない列及び行がある場合、前記第１記憶領域に入りきらない列及び行のデータ量に相当する記憶容量の第２記憶領域を、前記記憶手段に確保し、
前記第１記憶領域に入りきらない列及び行に対して、前記第２記憶領域を使用して分解処理を実施する、
処理を実行させる行列分解プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、計算機、行列分解方法、及び行列分解プログラムに関する。
【背景技術】
【０００２】
偏微分方程式などで表される数理モデルに基づきシミュレーション（流体解析、回路解析、心臓などの生体解析など）や数理計画問題を解く際、スパース行列（疎行列）で表される連立１次方程式の解を計算することがある。スパース行列とは、成分のほとんどが０の行列である。スパース行列の連立１次方程式は、スパース行列をＬＵ分解またはＬＤＬT分解して解くことができる。ＬＵ分解とは、正方行列を下三角行列（Ｌ：lower triangular matrix）と上三角行列（Ｕ：upper triangular matrix）との積に分解することである。ＬＤＬT分解は、正方行列を下三角行列（Ｌ）、対角行列（Ｄ：diagonal matrix）、及び下三角行列の転置行列（ＬT）の積に分解することである。
【０００３】
スパース行列の連立１次方程式を、ＬＵ分解やＬＤＬT分解を用いて解く方法として、例えばスーパーノーダル法がある。スーパーノーダル法では、非ゼロ要素のパターンが近い列を集めてスーパーノードが構成され、スーパーノードに対応する分解結果を格納できるだけの記憶領域が確保され、その記憶領域を使い数値的な分解が行われる。
【０００４】
また、ＬＵ分解、ＬＤＬT分解などの行列の分解では、ピボットと呼ばれる対角要素を含む行（ピボット行）の各要素の定数倍を、ピボット行より下の行の各要素から減算し、ピボットの下の要素を「０」にする処理が行われる。この処理はピボットの値が十分に大きいことが前提であり、ピボットとなる対角要素が小さな値やゼロとなると、分解が続けられなくなる。これを避けるための手法に、delayed pivotsと呼ばれる手法がある。delayed pivotsでは、列及び行の移動操作により、ピボットとするのに不適切な値を有する対角要素が後ろ（行番号と列番番号が大きくなる方向）に移動される。該当対角要素が後ろに移動すると、それより前の行に対する分解処理の過程で、該当対角要素の値が更新され、ピボットとしての条件を満たすような大きな値になることが期待できる。
【０００５】
スパース行列を用いたシミュレーションとしては、例えば電磁界解析などの数値シミュレーションにおいて最終的に解くべき方程式として現れる大型疎行列方程式の計算を高速に行うことができる行列方程式計算方法がある。また、疎行列を係数行列とする連立一次方程式の解を求める連立一次方程式求解方法の一例として、高速なシンボリック分解手法を持つ連立一次方程式求解方法も考えられている。さらに、ＳＯＲ（Successive Over-Relaxation）法（逐次的過剰緩和法）もしくはＳＳＯＲ（Symmetirc SOR）法（対称逐次的過剰緩和法）を用いた連立一次方程式の解法において、問題に応じて緩和パラメタの適切な値を高速かつ自動的に決定する技術も考えられている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１０−１２２８５０号公報
【特許文献２】特開２００９−２５９６２号公報
【特許文献３】特開２０１５−４９７２４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
スパース行列で表される連立１次方程式の、スーパーノーダル法を用いた直接解法に、delayed pivotsを適用すると、列及び行の並べ替えにより、スーパーノードを構成する行列要素が動的に変わる。そして列及び行を後ろに移動するときに、移動する列及び行の要素を格納するための領域を設けることとなり、使用するメモリ量も動的に変化する。しかし、使用するメモリ量が増加する度にメモリの記憶領域を動的に確保していたのでは、処理が煩雑になりすぎる。すなわち、スーパーノーダル法とdelayed pivotsとの両方を適用する場合、メモリ管理負荷が過大となるという問題がある。
【０００８】
１つの側面では、本件は、行列分解時のメモリ管理の効率化を図ることを目的とする。
【課題を解決するための手段】
【０００９】
１つの案では、非ゼロ要素の配置が対角要素を挟んで対称である行列を、下三角行列と上三角行列とを含む複数の行列に分解する計算機が提供される。当該計算機は、記憶手段、第１記憶領域確保手段、第１分解手段、第２記憶領域確保手段、及び第２分解手段を有数する。
【００１０】
記憶手段は、行列の分解結果を記憶する。第１記憶領域確保手段は、連続する対角要素を共有する列群と行群とを分解の処理単位として、非ゼロ要素の配置に基づいて行列を複数の処理単位に分け、複数の処理単位それぞれについて、含まれる列及び行分のデータ量に、所定数の列及び行分のデータ量を加算した合計データ量に相当する記憶容量の第１記憶領域を、記憶手段に確保する。第１分解手段は、複数の処理単位それぞれを所定の順番で対象処理単位とし、対象処理単位の列及び行と、対角要素の値が所定値以下であるために分解処理実施済みの処理単位から対象処理単位に移動された列及び行とに対して、第１記憶領域を使用して分解処理を実施する。第２記憶領域確保手段は、対象処理単位に移動された列及び行、並びに対象処理単位の分解処理の結果、対角要素の値が所定値以下であるために移動することになった列及び行のうち、第１記憶領域に入りきらない列及び行がある場合、第１記憶領域に入りきらない列及び行のデータ量に相当する記憶容量の第２記憶領域を、記憶手段に確保する。第２分解手段は、第１記憶領域に入りきらない列及び行に対して、第２記憶領域を使用して分解処理を実施する。
【発明の効果】
【００１１】
１態様によれば、行列分解時のメモリ管理の効率化を図ることができる。
【図面の簡単な説明】
【００１２】
【図１】第１の実施の形態に係る計算機の構成例を示す図である。
【図２】第２の実施の形態に用いるコンピュータのハードウェアの一構成例を示す図である。
【図３】delayed pivotsを説明する図である。
【図４】スーパーノードを示す図である。
【図５】ＬＬT分解の一例を示す図である。
【図６】生成されるelimination treeとrow subtreeの一例を示す図である。
【図７】コンピュータの機能を示すブロック図である。
【図８】スパースでありかつ構造的に対称な実行列のＬＵ分解の処理手順を示すフローチャートである。
【図９】スパースでありかつ構造的に対称な実行列のＬＵ分解におけるパネル構造の一例を示す図である。
【図１０】領域管理情報の一例を示す図である。
【図１１】行列分解作業領域のデータ構造例を示す図である。
【図１２】スーパーノード間の行と列の移動を示す図である。
【図１３】delayed pivotsの移動例を示す図である。
【図１４】移動させるdelayed pivotsの一例を示す図である。
【図１５】delayed pivots移動先のスーパーノードを示す図である。
【図１６】プライマリスーパーノードの分解に応じて発生したdelayed pivotsの移動例を示す図である。
【図１７】ＬＵ分解処理の手順を示すフローチャートである。
【図１８】サブルーチンrsupdateの処理を示すフローチャートである。
【図１９】サブルーチンdpcountの処理を示すフローチャートである。
【図２０】サブルーチンmvtonporderの処理を示すフローチャートである。
【図２１】サブルーチンcpupdateの処理を示すフローチャートである。
【図２２】サブルーチンrpupdateの処理を示すフローチャートである。
【図２３】サブルーチンcreatemvの処理を示すフローチャートである。
【図２４】サブルーチンcpupdatenewの処理を示すフローチャートである。
【図２５】サブルーチンrpupdatenewの処理を示すフローチャートである。
【図２６】スパースな不定値対称行列のＬＤＬT分解の処理手順を示すフローチャートである。
【図２７】スパースな不定値対称行列のＬＤＬT分解におけるパネル構造の一例を示す図である。
【図２８】不定値対称行列のＬＤＬT分解におけるdelayed pivotsの移動例を示す図である。
【図２９】ＬＤＬT分解処理の手順を示すフローチャートである。
【図３０】サブルーチンsymrsupdateの処理を示すフローチャートである。
【図３１】サブルーチンsymdpcountの処理を示すフローチャートである。
【図３２】サブルーチンsymmvtonporderの処理を示すフローチャートである。
【図３３】サブルーチンsymcpupdateの処理を示すフローチャートである。
【図３４】サブルーチンsymcreatemvの処理を示すフローチャートである。
【図３５】サブルーチンsymcpupdatenewの処理を示すフローチャートである。
【発明を実施するための形態】
【００１３】
以下、本実施の形態について図面を参照して説明する。なお各実施の形態は、矛盾のない範囲で複数の実施の形態を組み合わせて実施することができる。
〔第１の実施の形態〕
図１は、第１の実施の形態に係る計算機の構成例を示す図である。第１の実施の形態に係る計算機１０は、非ゼロ要素の配置が対角要素を挟んで対称である行列１を、下三角行列と上三角行列とを含む複数の行列に分解するものである。例えば行列１は、要素のほとんどが０のスパース行列である。計算機１０は、行列１のＬＵ分解、またはＬＤＬT分解を行う。計算機１０は、行列１を分解するために、記憶手段１１、親子関係決定手段１２、第１記憶領域確保手段１３、第１分解手段１４、第２記憶領域確保手段１５、及び第２分解手段１６を有する。
【００１４】
記憶手段１１は、行列１の分解結果を記憶する。行列１は、複数の処理単位２〜６に分けられる。処理単位２〜６は、非ゼロ要素の配置が同じ列群、及びその列群と対角要素を共有する行群の集合である。処理単位２〜６ごとに、分割処理が実施される。
【００１５】
親子関係決定手段１２は、所定の規則に従って、複数の処理単位２〜６間の親子関係を決定する。例えば処理単位３の分割処理の結果、処理単位４内の要素の値が更新される場合、２つの処理単位３，４の関係では、処理単位３が子、処理単位４が親となる。
【００１６】
第１記憶領域確保手段１３は、連続する対角要素を共有する列群と行群とを分解の処理単位として、非ゼロ要素の配置に基づいて行列１を複数の処理単位２，３に分ける。例えば第１記憶領域確保手段１３は、非ゼロ要素の配列が同一または近似する行群及び列群を、分解の処理単位とする。このような分解の処理単位は、スーパーノーダル法におけるスーパーノードである。そして第１記憶領域確保手段１３は、複数の処理単位それぞれについて、処理単位に含まれる列及び行分のデータ量に、所定数の列及び行分のデータ量（図１中の網掛け部分）を加算した合計データ量に相当する記憶容量の第１記憶領域１１ａを、記憶手段１１に確保する。
【００１７】
第１分解手段１４は、複数の処理単位２〜６それぞれを所定の順番で対象処理単位とする。そして第１分解手段１４は、対象処理単位の列及び行と、対角要素の値が所定値以下であるために分解処理実施済みの処理単位から移動された列及び行とに対して、第１記憶領域１１ａを使用して分解処理を実施する。例えば第１分解手段１４は、分解処理実施済みの処理単位のうち、対象処理単位の子に該当する処理単位の分解処理において、対角要素の値が所定値以下であるために処理できなかった列及び行を、対象処理単位に属する列及び行の次に移動する。このような列及び行の移動は、delayed pivotsと呼ばれる処理である。delayed pivotsで行の入れ替えだけでなく列も入れ替えることで、行列１の非ゼロ要素の配置の対称性を崩さずに済み、使用するメモリ容量が増加することが抑止できる。
【００１８】
なお、delayed pivotsで第１記憶領域１１ａに移動されるのは、第１記憶領域１１ａに入りきる分の列及び行のみである。delayed pivotsで対象処理単位へ移動する列及び行の数が過多となり、対象処理単位の第１記憶領域１１ａから溢れる列及び行がある場合、溢れた分については、第２記憶領域１１ｂの確保を待ってから移動される。従って、第１分解手段１４は、第１記憶領域１１ａに格納可能な分の列及び行に対してのみ分解処理を実施することとなる。
【００１９】
第２記憶領域確保手段１５は、分解処理実施済みの処理単位から移動された列及び行、並びに対象処理単位の分解処理の結果、対角要素の値が所定値以下であるために移動することになった列及び行が、すべて第１記憶領域１１ａに収まったか否かを判断する。そして第２記憶領域確保手段１５は、第１記憶領域１１ａに入りきらない列及び行がある場合、第１記憶領域１１ａに入りきらない列及び行のデータ量に相当する記憶容量の第２記憶領域１１ｂを、記憶手段１１に確保する。
【００２０】
第２分解手段１６は、第１記憶領域に入りきらない列及び行に対して、第２記憶領域１１ｂを使用して分解処理を実施する。
このような計算機１０によれば、行列１を分解する際に、処理単位２〜６それぞれについて、余分なスペースを含む第１記憶領域１１ａが確保される。そして、処理単位２〜６の分割処理において、第１記憶領域１１ａを超える量のdelayed pivotsによる列及び行の移動が発生した場合にのみ、第２記憶領域１１ｂが確保される。これにより、delayed pivotsによる列及び行の移動が、第１記憶領域１１ａに設けた余分なスペースに収まる場合、第２記憶領域１１ｂを確保せずに済み、メモリ管理が効率的である。
【００２１】
しかも、第２記憶領域１１ｂを確保する際には、delayed pivotsによって移動される多数の列及び行のうち、第１記憶領域１１ａに入りきらない部分に応じた記憶容量の第２記憶領域１１ｂが纏めて確保される。すなわち、追加で行われる記憶領域確保処理は、１つの処理単位に対して１回で済む。そのためdelayed pivotsによる列及び行の移動が発生するごとに記憶領域を確保する場合に比べ、メモリ管理の処理が効率化されている。
【００２２】
なお、第１分解手段１４と第２分解手段１６とは、列及び行の移動において、移動対象がｉ（ｉは１以上の整数）番目の行とｉ番目の列の場合、ｉ列内の行番号がｉ以上の要素を、対象処理単位に属する列の次に移動する。また第１分解手段１４と第２分解手段１６とは、ｉ行内の列番号がｉより大きい要素を、対象処理単位に属する行の次に移動する。すなわち、行列１における、行番号と列番号とのいずれもがｉより大きい要素のみが移動対象となる。これにより、移動先に設ける記憶容量を削減することができる。
【００２３】
なお、親子関係決定手段１２、第１記憶領域確保手段１３、第１分解手段１４、第２記憶領域確保手段１５，及び第２分解手段１６は、例えば計算機１０が有するプロセッサにより実現することができる。また、記憶手段１１は、例えば計算機１０が有するメモリまたはストレージ装置により実現することができる。
【００２４】
また、図１に示した各要素間を接続する線は通信経路の一部を示すものであり、図示した通信経路以外の通信経路も設定可能である。
〔第２の実施の形態〕
次に第２の実施の形態について説明する。第２の実施の形態は、偏微分方程式などで表される数理モデルに基づきシミュレーションや数理計画問題解く際に、スパース行列で表される連立１次方程式の解を効率的に計算するコンピュータである。
【００２５】
＜ハードウェア構成＞
まず第２の実施の形態に係るコンピュータのハードウェア構成について説明する。
図２は、第２の実施の形態に用いるコンピュータのハードウェアの一構成例を示す図である。コンピュータ１００は、プロセッサ１０１によって装置全体が制御されている。プロセッサ１０１には、バス１０９を介してメモリ１０２と複数の周辺機器が接続されている。プロセッサ１０１は、マルチプロセッサであってもよい。プロセッサ１０１は、例えばＣＰＵ（Central Processing Unit）、ＭＰＵ（Micro Processing Unit）、またはＤＳＰ（Digital Signal Processor）である。プロセッサ１０１がプログラムを実行することで実現する機能の少なくとも一部を、ＡＳＩＣ（Application Specific Integrated Circuit）、ＰＬＤ（Programmable Logic Device）などの電子回路で実現してもよい。
【００２６】
メモリ１０２は、コンピュータ１００の主記憶装置として使用される。メモリ１０２には、プロセッサ１０１に実行させるＯＳ（Operating System）のプログラムやアプリケーションプログラムの少なくとも一部が一時的に格納される。また、メモリ１０２には、プロセッサ１０１による処理に利用する各種データが格納される。メモリ１０２としては、例えばＲＡＭ（Random Access Memory）などの揮発性の半導体記憶装置が使用される。
【００２７】
バス１０９に接続されている周辺機器としては、ストレージ装置１０３、グラフィック処理装置１０４、入力インタフェース１０５、光学ドライブ装置１０６、機器接続インタフェース１０７及びネットワークインタフェース１０８がある。
【００２８】
ストレージ装置１０３は、内蔵した記憶媒体に対して、電気的または磁気的にデータの書き込み及び読み出しを行う。ストレージ装置１０３は、コンピュータの補助記憶装置として使用される。ストレージ装置１０３には、ＯＳのプログラム、アプリケーションプログラム、及び各種データが格納される。なお、ストレージ装置１０３としては、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）を使用することができる。
【００２９】
グラフィック処理装置１０４には、モニタ２１が接続されている。グラフィック処理装置１０４は、プロセッサ１０１からの命令に従って、画像をモニタ２１の画面に表示させる。モニタ２１としては、ＣＲＴ（Cathode Ray Tube）を用いた表示装置や液晶表示装置などがある。
【００３０】
入力インタフェース１０５には、キーボード２２とマウス２３とが接続されている。入力インタフェース１０５は、キーボード２２やマウス２３から送られてくる信号をプロセッサ１０１に送信する。なお、マウス２３は、ポインティングデバイスの一例であり、他のポインティングデバイスを使用することもできる。他のポインティングデバイスとしては、タッチパネル、タブレット、タッチパッド、トラックボールなどがある。
【００３１】
光学ドライブ装置１０６は、レーザ光などを利用して、光ディスク２４に記録されたデータの読み取りを行う。光ディスク２４は、光の反射によって読み取り可能なようにデータが記録された可搬型の記録媒体である。光ディスク２４には、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）などがある。
【００３２】
機器接続インタフェース１０７は、コンピュータ１００に周辺機器を接続するための通信インタフェースである。例えば機器接続インタフェース１０７には、メモリ装置２５やメモリリーダライタ２６を接続することができる。メモリ装置２５は、機器接続インタフェース１０７との通信機能を搭載した記録媒体である。メモリリーダライタ２６は、メモリカード２７へのデータの書き込み、またはメモリカード２７からのデータの読み出しを行う装置である。メモリカード２７は、カード型の記録媒体である。
【００３３】
ネットワークインタフェース１０８は、ネットワーク２０に接続されている。ネットワークインタフェース１０８は、ネットワーク２０を介して、他のコンピュータまたは通信機器との間でデータの送受信を行う。
【００３４】
以上のようなハードウェア構成によって、第２の実施の形態の処理機能を実現することができる。なお、第１の実施の形態に示した装置も、図２に示したコンピュータ１００と同様のハードウェアにより実現することができる。
【００３５】
コンピュータ１００は、例えばコンピュータ読み取り可能な記録媒体に記録されたプログラムを実行することにより、第２の実施の形態の処理機能を実現する。コンピュータ１００に実行させる処理内容を記述したプログラムは、様々な記録媒体に記録しておくことができる。例えば、コンピュータ１００に実行させるプログラムをストレージ装置１０３に格納しておくことができる。プロセッサ１０１は、ストレージ装置１０３内のプログラムの少なくとも一部をメモリ１０２にロードし、プログラムを実行する。またコンピュータ１００に実行させるプログラムを、光ディスク２４、メモリ装置２５、メモリカード２７などの可搬型記録媒体に記録しておくこともできる。可搬型記録媒体に格納されたプログラムは、例えばプロセッサ１０１からの制御により、ストレージ装置１０３にインストールされた後、実行可能となる。またプロセッサ１０１が、可搬型記録媒体から直接プログラムを読み出して実行することもできる。
【００３６】
＜行列分割処理の概要＞
図２に示したコンピュータ１００は、delayed pivotsを適用したスーパーノーダル法を用い、スパース行列の連立１次方程式を、ＬＵ分解、またはＬＤＬT分解する。
【００３７】
まず、コンピュータ１００が実施する行列分割処理の概要について説明する。スーパーノーダル法により行列を分解する場合、コンピュータ１００は、数値的な分解に先立ち、シンボリック（symbolic）分解を行う。コンピュータ１００は、シンボリック分解の結果に基づいて、非ゼロ要素のfill-in（非ゼロを入れること）や列番号（行番号）に対応したノードの依存関係などの情報を獲得する。ここでノードとは、対角要素を共有する行と列との組である。ノードの依存関係とは、一方のノードの分割処理により、他方のノード内の要素の値が更新されるような関係である。この場合、分割処理が行われたノードが子ノード、子ノードの分割処理により要素の値が更新されるノードが親ノードと呼ばれる。シンボリック分解の結果から得た情報に基づいて、コンピュータ１００は、非ゼロ要素のパターンが近い連続するノードを集めてスーパーノードとし、スーパーノード間の依存関係を示すsupernodal elimination tree, supernodal row subtreeを構成する。スーパーノードには、１または複数のノードが含まれる。
【００３８】
コンピュータ１００は、スーパーノードに対応する分解結果を格納するパネルの大きさを計算する。そしてコンピュータ１００は、パネルの大きさで、スーパーノードに割り付けたメモリ領域を使い、数値的な分解を行う。
【００３９】
コンピュータ１００は、分解の安定性を高め精度を向上するために、ピボットをスーパーノードに対応する対角ブロック（１または複数の対角要素を含む領域）内で取るが、条件を満たすものが取れなかったとき、delayed pivotsを実施する。
【００４０】
コンピュータ１００は、スーパーノーダル法においてdelayed pivotsを実施する場合、スーパーノードの分解結果を格納するパネルとして、列・行の移動を考慮したスペースを加えた記憶領域を設定する。追加する空きスペースの大きさを表す追加の本数は、シンボリック分解の前に、コンピュータ１００に指定されている。コンピュータ１００は、スペースを追加したパネルの大きさをシンボリック分解で得た情報から計算し、数値的分解前に、スーパーノードに割り付ける。この時点でパネルの割り付けを受けるスーパーノードを、プライマリスーパーノードとする。
【００４１】
コンピュータ１００は、追加したスペースでは記憶容量が不足した場合、移動先のスーパーノードにセカンダリスーパーノードを作成し、格納用のパネルを割り付けてdelayed pivotsを移動する。コンピュータ１００は、移動先のスーパーノードのＬＵ分解やＬＤＬT分解で生じた新たなdelayed pivotsの候補にスペースに入らなかったdelayed pivotsの残りに加えた分が格納できるような大きさのセカンダリスーパーノードの大きさを割り付ける。
【００４２】
セカンダリスーパーノード用のパネルを割り当てる、指定した大きさのメモリプールをスーパーノードのパネルを割り当てるときに同じく割り当て、そこからセカンダリスーパーノードに対するパネルを動的に確保して割り当てる。
【００４３】
このように、スーパーノーダル法に、delayed pivotsを組み込むために、スーパーノードのパネルにスペースを追加し、delayed pivotsを行うためのスペースが不足したときセカンダリスーパーノードを作成してdelayed pivotsを行う。
【００４４】
＜delayed pivots＞
以下、delayed pivotsについて詳細に説明する。
連立１次方程式の直接解法で行列をＬＵ分解するとき、ピボットとなる対角要素が小さな値やゼロとなって、分解が続けられなくなる場合がある。これを避けるためにdelayed pivotsと呼ばれる手法が用いられる。
【００４５】
図３は、delayed pivotsを説明する図である。delayed pivotsでは、ピボットとして不適切な値の対角要素を対称置換（Symmetric permutation）で後ろに移動することを考える。図３の例では、ｉ行ｊ列（ｉ，ｊは、１以上の整数）の対角要素が後ろに移動する場合を想定している。この場合、コンピュータ１００は、まず、ｉ行目とｉ＋１行目を対称置換する。その結果、２つの列及び行が入れ替わる。
【００４６】
例えば、行列Ｐを置換行列として、行列Ｐのａ行ｂ列（ａ，ｂは１以上の整数）の要素をｐa,bとする。この行列Ｐは、ｐi,i＝ｐi+1,i+1＝0，ｐi+1,i＝ｐi,i+1＝１でｐj,j＝１（ｊはｉ，ｉ＋１以外）であり、他の要素は「０」の直交行列であるものとする。分割対象の行列Ａに行列Ｐを左から掛けると、ｉ行とｉ＋１行が入れ替わる。また分割対象の行列Ａに行列Ｐを右から掛けると、ｉ列とｉ＋１列が入れ替わる。入れ替え後の行列をＢとすると、入れ替え操作は、以下の式で表される。
Ｂ＝ＰＡＰT
このような操作は、ｉ行ｉ列をｉ＋１行ｉ＋１列に対称置換することに相当する。同様の入れ替え操作を（ｉ，ｉ＋１），（ｉ＋１，ｉ＋２），…，（ｊ−１，ｊ）まで引き続いて行うと、ｉ＋１，・・・，ｊまでをｉ，・・・，ｊ−１に移動して、ｉ列、ｉ行をｊ列、ｊ行に移すことになる。
【００４７】
これはＰi，Ｐi+1，・・・，Ｐj-1を両側から順次掛けることに相当する。ただしＰiは左から掛けてｉ行とｉ＋１行を入れ替える直交行列とする。
このような移動を行って分解を続けると、引き続き行われる各列、行の分解で（ｊ，ｊ）の要素に更新がかかり、この要素の絶対値が大きくなることで、ピボットとしての条件を満たすようになることが期待できる。移動後に分解を続けても該当要素の値が十分大きくならなかった場合、さらに後ろに移動させることが繰り返される。
【００４８】
このように、ピボットを含む行と列とを後ろに対象置換していく処理が、delayed pivotsである。
ｉ列、ｉ行とｉ＋１列、ｉ＋１行を入れ替えた場合、ｉ行ｉ列の要素ａi,iは、ｉ＋１行ｉ＋１列の要素ａi+1,i+1と入れ替わる。その結果、ｉ列、ｉ行のＬＵ分解と外積型の更新を行うと、以下のように更新される。
ａi+1,i+1＝ａi+1,i+1−ａi+1,i×（ａi,i+1／ａi,i）
このため、移動して更新することでａi+1,i+1の絶対値が大きくなりピボットとして受け入れられるようになることが期待できる。
【００４９】
＜delayed pivotsのスーパーノーダル法への適用＞
スパース行列の連立１次方程式の直接解法であるスーパーノーダル法における行列分解中にdelayed pivotsが実施する場合を考える。スパース行列の連立１次方程式の直接解法であるスーパーノーダル法では、スーパーノードの対角ブロックの中で各種ピボットを取る。そのため、依存関係における親のスーパーノードの最後に対象の列・行を移動すれば、適当な入れ替えの後、妥当なピボットが選ばれるようになる。delayed pivotsの対象となった列・行をそのスーパーノードの親のスーパーノードの最後に移動すると、親のスーパーノードのＬＵ分解の結果を使った更新が移動した部分にかかることになる。
【００５０】
スパース行列の直接解法では、分解でのfill-inや依存関係を分解前に決定する。これをシンボリック分解と呼ぶ。対角成分より行番号が大きな非ゼロ要素からなる列はノードと呼ばれる。そしてこれらの情報を使って非ゼロパターンの近い列をまとめてスーパーノードを構成する。
【００５１】
図４は、スーパーノードを示す図である。スーパーノードは、上部三角領域の要素が非ゼロであり、各列が同じ非ゼロ構造を有する列の集合である。
スーパーノードの構成ノードを対角要素にもつ行で列番号が対角要素より大きなものをスーパーノードに対応する行と呼ぶ。そして、同じスーパーノードに対応する行をまとめることで、スーパーノードに対応する列と行が決まる。
【００５２】
そしてコンピュータ１００は、ＬＵ分解の結果はスーパーノードごとに、列をまとめたpanelと行をまとめたpanelに格納する。これらpanelは２次元配列とみなせる。またコンピュータ１００は、対角ブロックは、列のpanelに格納する。
【００５３】
コンピュータ１００は、シンボリック分解を行って各スーパーノードをＬＵ分解した分解結果を格納する領域の大きさを計算する。数値的分解を行う前に、計算した領域の大きさで格納領域を割り付ける。その後数値的な分解を行う。そして数値的分解の安定化のため、コンピュータ１００は、delayed pivotsを実施する。
【００５４】
＜delayed pivotsのスーパーノーダル法への適用時の課題＞
一般行列に関しては、対角要素に大きな要素を並べるようにする列の入れ替え（permutation）を求めるアルゴリズムが開発されている。それを適用すると、対角優位に近い行列に変換してＬＵ分解ができる。対角優位であれば、絶対値が閾値（threshold）を超えるものをピボットとして受け入れる閾値ピボティング（threshold pivoting）を行うと近くの要素でピボットサーチが完了することが多い。しかし構造的に対称な行列や不定値対称行列にこの前処理を適用すると元の行列が持っていた対称性が壊れてしまい、元の対称性を利用できなくなる。そのため、対称性を保って処理する上では、安定性に対してdelayed pivotsなどを組み入れることが重要になる。なお、構造的に対称であるとは、非ゼロ要素の位置が、対角要素を挟んで対称の位置に存在することである。
【００５５】
一般実行列のスーパーノーダル法では、スーパーノード内でピボットサーチを行っている。見つけられなかったときは、pivotを適当な大きさの数（倍精度実数の場合1.0D-8〜1.0D-10くらいの大きさ（おおよそ半分の精度くらい））に置き替えるstatic pivotを行って近似的に分解を行う。このとき分解の精度は劣化するため、精度を上げて（倍精度なら４倍精度を使って）反復改良を行うことで解の精度を回復する処理が行われる。
【００５６】
一般実行列Ａに対しては、それを包含するＡ＋ＡTを考え、そのelimination treeなどを使いＬＵ分解する方法がある。この方法にも対称性を壊さないdelayed pivotsを組み込むことで効果が見込める。つまり、分解の精度劣化を起こさず反復改良の回数を減らすことで性能向上が見込める。
【００５７】
スパース行列の連立１次方程式の直接解法に、delayed pivotsの処理を追加しようとすると、スーパーノードを構成する行列要素が動的に変わり、列及び行を後ろに移動するときに、移動する列及び行の要素を格納するための領域が増大する。さらにdelayed pivotsを実際行う時点になってからでないと正確な増分は分からない。そのため、格納域の大きさを動的に計算して、判明した時点で、記憶領域の確保・解放などを行うこととなる。
【００５８】
実際multifrontal法では、分解の中間結果を動的に確保した領域に一時的に保存し、依存関係を調べながらその都度分解の対象となる列・行を収集して処理を行う。このためメモリが多く使用され、処理も複雑である。
【００５９】
また、スーパーノーダル法では、数値分解の前にシンボリック分解を行い、各ノードの依存関係を解析し木構造で表現する。この木構造の情報を使って数値的な分解を行う。そして、スーパーノーダル法では、シンボリック分解で分解結果を格納する領域の大きさを計算し、分解結果を格納する領域をnumericな分解を行う前に割り付ける。しかし、delayed pivotsを動的に行うと利用する記憶領域の大きさが変化するため、ほぼすべてのメモリを動的に管理することとなり、処理が煩雑になりすぎる。
【００６０】
＜スーパーノード間の依存関係＞
第２の実施の形態では、シミュレーションや数理計画問題を解く際に、なるべくシンボリック分解の情報を利用して、事前に領域を割り付け、その結果依存関係を表す木構造を解析して行う処理方式を継承できるようにする。
【００６１】
ここで、第２の実施の形態の機能を説明する前に、スーパーノード間の依存関係の判別に用いるelimination tree, row subtreeの概要を説明する。
スパースな正値対称行列のＬＬT分解では、elimination tree, row subtreeが使われる。このelimination tree, row subtreeは、構造的対称な実行列や不定値対称行列にも適用できる。そこで、スパースな正値対称行列のＬＬT分解で使われるelimination tree, row subtreeについて説明する。
【００６２】
スパースな正値対称行列Ｐの非ゼロパターンから分解結果の行列Ｌの各列の間の依存関係を表すelimination treeが得られる。ＰはＰ＝ＬＬTのようにコレスキー分解できる。行列Ｌのｉ行ｊ列の要素をＬijとしたとき、ｍｉｎ｛ｉ｜ｉ＞ｊかつＬij≠０｝がｊの親となる。
【００６３】
親を辿り、これ以上親が辿れないノードをelimination treeのrootノードと呼ぶ。あるノードｑの親をｐとしたとき、ｑはｐの子ノードという。elimination treeのrootノードからdepth first searchして得られた番号をそのノードのpostorderと呼ぶ。elimination treeのrootノードからdepth first searchしてこれ以上深くならないノードをleafノードと呼ぶ。各leafノードから親ノードを辿りelimination treeのrootノードまで遡ったときに辿った各ノードに、このleafノードが対応付けられる。あるノードにleafノードが複数対応付けられたとき、これらのleafノードの中でpostorderのもっとも小さいものを、そのノードのfirst descendantと呼ぶ。
【００６４】
もとの正値対称行列の下三角行列のｊ番目の列をｂjとしたときＬの非ゼロパターンは、ｂjとｊの子ノードのＬの列ｌkのunion になる。このことより、Ｌのｉ番目の行の非ゼロ要素はelimination treeの部分木として表現できる。この部分木は、ｉ番目のノードをrootとするrow subtreeになる。
【００６５】
以下、図５、図６を参照して、正値対称行列に関するＬＬT分解の結果として得られる、elimination treeとrow subtreeとの具体例を説明する。
図５は、ＬＬT分解の一例を示す図である。図５では、上段に分割前の行列を示し、下段に分割後の下三角行列Ｌを示している。各行列の対角要素には、その要素の行番号が示されている。また各行列の対角要素は、非ゼロ要素であるものとする。
【００６６】
分割前の行列では、非ゼロ要素を黒丸で示している。この行列をＬＬT分解した場合、分解後の下三角行列Ｌには、新たに非ゼロ要素が追加されている。図５において、分割により非ゼロとなった要素を、網掛けの丸印で示している。ＬＬT分解によって生成された下三角行列（行列Ｌ）に基づいて、elimination treeとrow subtreeが生成できる。
【００６７】
図６は、生成されるelimination treeとrow subtreeの一例を示す図である。例えば図５に示す行列Ｌにおける第１列（ｊ＝１）では、ｍｉｎ｛ｉ｜ｉ＞ｊかつＬij≠０｝は６である。この場合、ノード「１」の親ノードはノード「６」となる。同様に、ノード「２」の親ノードはノード「３」となる。このような親子の関係をすべての列について判定し、親子関係を線で接続すると、図６に示すelimination tree９１が生成できる。
【００６８】
elimination tree９１のdepth first searchの順番は、以下のようになる。
「２→３→１→６→７→４→５→８→９→１０→１１」
そして、depth first searchの順番に従って、各ノードに１から１１のpostorderが振られる。
【００６９】
ここで図５を参照すると、分解前の行列の７行目は｛１，３，７｝に非ゼロ要素がある。図６に示すelimination tree９１では、ノード「１」はleafでありノード「７」まで遡れる。ノード「３」はノード「７」まで遡れる。そこで、elimination tree９１から部分木｛１，６，３，７｝が抽出され、ノード「７」をrootノードとするrow subtree９２となる。row subtree９２では、ノード「１」とノード「３」がleafである。
【００７０】
同様に、ノード「１１」のrow subtree９３も生成できる。
分解前の下三角行列（行列Ｌ）の行番号「１１」にある非ゼロ要素は対角要素を除き｛１，５，８｝である。これをpostorder順に取り出すと、順序は同じ「１→５→８」となる。ノード「１，５，８」それぞれのfirst descendantは「１，４，２」となる。
【００７１】
ノード「１」は、最初なのでleafである。ノード「５」のfirst descendantであるノード「４」のpostorderは、ノード「１」のpostorderより大きいので分枝ノードで枝分かれしている（depth firstなので分枝ノードは後回しにしてsearchする）。ノード「８」のfirst descendantであるノード「２」のpostorderは、ノード「５」のpostorderより小さいので、もとの枝に戻ったことになり、ノード「５」とノード「８」の間で枝分かれはない。
【００７２】
このように分解前の行の非ゼロ要素をpostorder順に取り出し、一つ前に取り出したノードのpostorderと今取り出したノードのfirst descendantのpostorderを比較する。今取り出したノードのfirst descendantの方が大きければ、今取り出したノードがこのrow subtreeのleafである。すなわち、depth first searchでpostorderを振っているので、この２つのノードのcommon ancestorで分枝しているときには、今取り出したノードのfirst descendantの方が大きくなる。
【００７３】
行の非ゼロ要素をpostorder順に取り出してひとつ前のノードを覚えておく代わりに、ひとつ前のleafノードを覚えておいて、新たなleafノードが見つかれば更新するようにしても結果は同じである。
【００７４】
row subtreeは各行の非ゼロ要素を表している。非ゼロ要素が構造的に対称な実行列の場合は、行列Ｌの行と対称な位置にあるＵの列の非ゼロ要素も表している。postorder順にノードに関する行列Ｌの列を更新するとき使われる行列Ｌの列及び行列Ｕの行を更新するとき使われる行列Ｕの行は、そのノードのrow subtreeにあらわれるノードということが知られている。
【００７５】
これらの関係は、行列Ｌの列の非ゼロ要素のパターンが近いノードを集めてひと塊にしたスーパーノードに対してもあてはまる。すなわち、複数のノードの集まりであるスーパーノードを１つのノードと捉えれば、図５、図６を参照して説明したのと同様の処理により、スーパーノード間の依存関係を示す木構造を生成できる。elimination tree及びrow subtreeと同様の手順で作成した、スーパーノード間の依存関係を示す木構造が、supernodal elimination tree及びsupernodal row subtreeとなる。
【００７６】
＜コンピュータの機能＞
図７は、コンピュータの機能を示すブロック図である。コンピュータ１００は、記憶部１１０と解析部１２０とを有する。
【００７７】
記憶部１１０は、領域管理情報１１１を記憶すると共に、行列を分解する際の演算に使用する行列分解作業領域１１２を有する。領域管理情報１１１は、行列分解作業領域１１２のサイズなどを示す管理情報である。行列分解作業領域１１２は、スーパーノードの分割時の要素を格納する記憶領域である。
【００７８】
解析部１２０は、連立１次方程式の求解を伴う各種解析処理を行う。例えば解析部１２０は、流体解析などのシミュレーションを行う。解析部１２０は、解析処理の過程で、連立１次方程式を解くために、スパース行列のＬＵ分解またはＬＤＬT分解を行う。
【００７９】
なお、図７に示した各要素の機能は、例えば、その要素に対応するプログラムモジュールをコンピュータに実行させることで実現することができる。図７に示す各要素のうち、解析部１２０は、図１に示した親子関係決定手段１２、第１記憶領域確保手段１３、第１分解手段１４、第２記憶領域確保手段１５、及び第２分解手段１６を包含する機能の一例である。また図７に示す記憶部１１０は、図１に示した記憶手段１１の一例である。
【００８０】
＜スパースでありかつ構造的に対称な実行列のＬＵ分解＞
まず、スパースでありかつ構造的に対称な実行列のＬＵ分解について説明する。
図８は、スパースでありかつ構造的に対称な実行列のＬＵ分解の処理手順を示すフローチャートである。以下、図８に示す処理をステップ番号に沿って説明する。
【００８１】
［ステップＳ１０１］解析部１２０は、column panelとrow panelに追加するスペースの大きさの指定入力を受け付ける。例えばユーザがキーボード２２などの入力デバイスを介して入力した、スペースの大きさを示す値を、解析部１２０が取得する。
【００８２】
［ステップＳ１０２］解析部１２０は、シンボリック分解を行う。このシンボリック分解により、elimination treeの生成、スーパーノードの検出、supernodal row subtreeの生成などが行われる。
【００８３】
［ステップＳ１０３］解析部１２０は、column panel, row panelそれぞれの、スペースを含んだ大きさを計算する。
［ステップＳ１０４］解析部１２０は、column panel用の領域とrow panel用の領域、及びセカンダリスーパーノード用の領域を割り当てるためのメモリプールをメモリ内に確保する。
【００８４】
［ステップＳ１０５］解析部１２０は、ＬＵ分解を行う。ＬＵ分解では、解析部１２０は、delayed pivotsの移動用に設けたスペースを利用する。また移動するdelayed pivotsの本数がスペースの大きさを超えた場合、解析部１２０は、セカンダリスーパーノードを作成して、delayed pivotsに対応する列及び行を格納する。
【００８５】
このような手順でＬＵ分解が進められる。以下、ＬＵ分解を実施する際の、メモリ１０２の有効活用方法について詳細に説明する。
スパースな構造的に対称な実行列のＬＵ分解を実施する場合、スーパーノードの要素をパネルと呼ばれる記憶領域に格納する。
【００８６】
図９は、スパースでありかつ構造的に対称な実行列のＬＵ分解におけるパネル構造の一例を示す図である。図９には、プライマリスーパーノードの分解結果を格納する領域であるプライマリスーパーノード領域３１のパネル構造と、プライマリスーパーノードの分解結果を格納する領域であるセカンダリスーパーノード領域３２のパネル構造とが示されている。
【００８７】
解析部１２０は、シンボリック分解で求めたスーパーノードの大きさにdelayed pivotsのためのスペース３３を設ける。解析部１２０は、スペース３３を、列をまとめて格納するcolumn panelと行をまとめて格納するrow panelの両方に設ける。各パネルは２次元であり、１次元目と２次元目との両方にスペース３３が設けられる。
【００８８】
ここで、プライマリスーパーノード領域３１の格納パネルに移動できるdelayed pivotsの最大数をnsp1とする。また、nsp1の分のスペースではスペースが不足したとき、動的に確保するセカンダリスーパーノード領域３２の格納panelの２次元目の大きさの最大数をnsp2とする。動的に確保する大きさは、プライマリスーパーノード領域３１の後ろに移動するdelayed pivotsの本数からnsp1を引いた数（nmpssp）である。つまりnmpsspがnsp2より小さい場合、処理を続けることができる。
【００８９】
なお、解析部１２０は、セカンダリスーパーノードにも１次元目にnsp3＝nsp1＋nsp2のスペースを確保しておく。
図９に示したようなパネル構造を示す情報が、領域管理情報１１１として記憶部１１０に格納される。そして、パネルと、パネルに付随する情報が行列分解作業領域１１２に格納される。パネルに付随する情報としては、非ゼロ要素の指標（index vector）とスーパーノード内で行ったピボットの入れ替えに関する情報（exchange history）がある。
【００９０】
以下、図１０と図１１とを参照し、記憶部１１０の領域管理情報１１１と行列分解作業領域１１２とについて詳細に説明する。
図１０は、領域管理情報の一例を示す図である。領域管理情報１１１には、共通情報１１１ａとスーパーノード別管理情報１１１ｂ，１１１ｃ，・・・とが含まれる。なお、図１０において、各データの名称に含まれるpspはprimary supernodeを表しており、sspはsecondary supernodeを表している。
【００９１】
共通情報１１１ａには、nsp1、nsp2、及びnsp3が含まれる。
スーパーノード別管理情報１１１ｂ，１１１ｃ，・・・には、プライマリスーパーノード情報１１１−１とセカンダリスーパーノード情報１１１−２とが含まれる。
【００９２】
プライマリスーパーノード情報１１１−１には、ipscp、ipsrp、ipslpindx、ipsupindx、ipsrex、ipscex、ndb、nbboff、nmppsp、npppspが含まれる。ipscpは、プライマリスーパーノードのcolumn panelを指し示すインデックスである。ipsrpは、プライマリスーパーノードのrow panelを指し示すインデックスである。ipslpindxは、プライマリスーパーノードのcolumn panelのインデックスリストを指し示すインデックスである。ipsupindxは、プライマリスーパーノードのrow panelのインデックスリストを指し示すインデックスである。ipsrexは、プライマリスーパーノードのrow exchangeのhistoryである。ipscexは、プライマリスーパーノードのcolumn exchangeのhistoryである。ndbは、シンボリック分解で求めたスーパーノードの大きさである。nbboffは、シンボリック分解で求めたスーパーノードの対角要素以外の大きさである。nmppspは、プライマリスーパーノードに移動されたdelayed pivotsの本数である。npppspは、プライマリスーパーノードにピボットとして残ったノード数である。
【００９３】
セカンダリスーパーノード情報１１１−２には、isscp、issrp、isslpindx、issupindx、issrex、isscex、nmpssp、nppsspが含まれる。isscpは、セカンダリスーパーノードのcolumn panelを指し示すインデックスである。issrpは、セカンダリスーパーノードのrow panelを指し示すインデックスである。isslpindxは、セカンダリスーパーノードのcolumn panelのインデックスリストを指し示すインデックスである。issupindxは、セカンダリスーパーノードのrow panelのインデックスリストを指し示すインデックスである。issrexは、セカンダリスーパーノードのrow exchangeのhistoryである。isscexは、セカンダリスーパーノードのcolumn exchangeのhistoryである。nmpsspは、セカンダリスーパーノードに移動されたdelayed pivotsの本数である。nppsspは、セカンダリスーパーノードにピボットとして残ったノード数である。
【００９４】
図１１は、行列分解作業領域のデータ構造例を示す図である。行列分解作業領域１１２には、すべてのプライマリスーパーノードについて、column panel４１、 row panel４２、column index４３、row index４４、row exchange history４５，column exchange history４６の各領域が確保される。各領域には、スペース３３の分の領域も含まれる。
【００９５】
スーパーノードごとの６つの領域の位置を示すインデックスは、それぞれに対応する１次元配列nofstcp, nofstrp, noftstcindx, nofstrindx, nofstrcx, nofstrexに格納される。ここで、スーパーノードの総数をnumspとしたとき、各１次元配列の大きさはnumsp＋１となる。numsp＋１番目の１次元配列には、領域の大きさ＋１の値が設定される。
【００９６】
あるスーパーノードに移動されるdelayed pivotsがプライマリスーパーノードに入りきらなかった場合、セカンダリスーパーノードが生成される。
生成されたセカンダリスーパーノードについても、図１１に示すような領域が確保される。セカンダリスーパーノードの大きさは、プライマリスーパーノードのパネルを分解（ＬＵ分解またはＬＤＬT分解）したとき生じたdelayed pivotsの数を、プライマリスーパーノードには入りきらず残ったdelayed pivotsに加えた本数の大きさになる。
【００９７】
プライマリスーパーノードとは別に、セカンダリスーパーノードに関するデータはメモリプール（メモリ１０２内の予約された記憶領域）に確保される。メモリプールの全体のサイズは、ユーザが予め指定したサイズである。セカンダリスーパーノードについての記憶領域は、予め確保されたメモリプール内から割り当てられる。図１１に示した各領域は、メモリプールの先頭からのオフセットを示すインデックスで特定される。これらインデックスに、例えばプライマリスーパーノードのpostorderによる番号を付与することで、識別できるようにすることができる。
【００９８】
例えば、２次元配列として２次元目でスーパーノードを指定し、１次元目に各種情報が入っているようにすることで、プライマリスーパーノード及びセカンダリスーパーノードに関するデータの格納域の場所などの情報を保持することができる。
【００９９】
次に、依存関係にあるスーパーノードについて説明する。第２の実施の形態では、依存関係にあるスーパーノードにおける子のスーパーノードから、delayed pivotsに対応する行と列が、そのスーパーノードの親のスーパーノードに移動する。
【０１００】
図１２は、スーパーノード間の行と列の移動を示す図である。子スーパーノード領域５１においてdelayed pivotsとされた対角要素を含む行と列が、親スーパーノード領域５２のスペース５２ａに移動する。例えば子スーパーノード内の最後尾（列番号が最も大きい）のノードがdelayed pivotsの移動対象となった場合、その子スーパーノード内では後ろに移動することはできないため、親スーパーノードに移動することとなる。この場合、親スーパーノード領域５２の非ゼロ要素が増加する。なお、子スーパーノードの兄弟のスーパーノードに関しては、依存関係がないので、非ゼロ要素は変化しない。
【０１０１】
ここで、シンボリック分解で決定されたスーパーノードへ、その子のスーパーノードから移動できる最大のdelayed pivotsの本数はnsp1+nsp2である。また、子のスーパーノードから移動する本数がnsp1を超えたときは、セカンダリスーパーノードが生成され、delayed pivotsに対応する行と列がセカンダリスーパーノードに移動する。このときのセカンダリスーパーノードの大きさは、プライマリスーパーノードで生じたdelayed pivotsの本数をその子のスーパーノードからの移動で残ったdelayed pivotsに加えた大きさである。
【０１０２】
このようにスペースを確保すればよいことは、以下のことより分かる。
前述の図３の例では、ｉ列ｉ行をｊ列ｊ行の後ろに移動している。このとき、Ａ（ｉ−ｎ，ｉ−ｎ）の部分行列を考えて、ｉ列ｉ行をｊの後ろに移動することを考える。この移動は、ｉ＋１〜ｊが親スーパーノードであって、１つ前の子スーパーノードの最後の列、行（ｉ列ｉ行）を、ｉ＋１〜ｊの親スーパーノードの後ろに移動する操作である。このとき移動先の親スーパーノードの列panelと行panelのスペースに、移動される列と行とが格納される。移動するdelayed pivotsの本数が増えると、使用されるスペースも増える。従って、親のスーパーノードの最後に移動すると移動本数分のスペースを確保すればよいことがわかる。移動してできたスーパーノードを通過して次のスーパーノードに移動するときも移動分のスペースを確保すればよい。
【０１０３】
次に、delayed pivotsの移動について詳細に説明する。
図１３は、delayed pivotsの移動例を示す図である。対称置換で移動される部分は、図１３の網掛けで示した部分になる。この移動される部分は、移動元のスーパーノードにおける、縦方向の下側、横方向の右側の縁に沿った部分である。すなわち、delayed pivotsの移動は、元の行列Ａ0のうちの、ｉ行以降及びｉ列以降の部分行列Ａk内に制限されている。このように、部分行列内に制限してdelayed pivotsに対応する行と列とを移動させることで、記憶容量の削減や処理の効率化が図れる。
【０１０４】
以下、複数のブロックからなるdelayed pivotsの動きを示す。
図１４は、移動させるdelayed pivotsの一例を示す図である。図１４では、最初のスーパーノードの分解で対角ブロックａ１，ｂ２，ｃ３で示される３本のdelayed pivotsが生じている。
【０１０５】
これらのうち対角ブロックａ１，ｂ２で示されるdelayed pivotsを対角ブロックｄ４で示されるスーパーノードにおけるプライマリスーパーノードのスペースに移動するものとする。対角ブロックｃ３で示される部分も合わせて移動され、セカンダリスーパーノードに格納されるものとする。
【０１０６】
図１５は、delayed pivots移動先のスーパーノードを示す図である。図１５の例では、対角ブロックがｄ４であるプライマリスーパーノードに対角ブロックがａ１，ｂ２のdelayed pivotsが移動されている。対角ブロックｃ３に対応するdelayed pivotsはプライマリスーパーノードに入りきらなかった分であり、セカンダリスーパーノードに格納されている。
【０１０７】
図１５に示すように、シンボリック分解で求めたスーパーノードの大きさndbに、nsp1の分のスペース、及びセカンダリスーパーノードの大きさであるnsp2に相当する領域内に、移動対象の要素が収まっている。
【０１０８】
ここで、プライマリスーパーノードの分解により、対角ブロックｂ２の部分がdelayed pivotsであると判断された場合を想定する。この場合、対角ブロックｂ２の部分は、対角ブロックｃ３で示されるdelayed pivotsの後ろに移動する。
【０１０９】
図１６は、プライマリスーパーノードの分解に応じて発生したdelayed pivotsの移動例を示す図である。図１６の例では、余ったdelayed pivots（対角ブロックｃ３）とプライマリスーパーノードの分解で生じたdelayed pivots（対角ブロックｂ２）とを格納する分のセカンダリスーパーノードが生成されている。そして生成されたセカンダリスーパーノードに、対角ブロックｃ３，ｂ２で示されるdelayed pivotsに対応する行と列が格納されている。
【０１１０】
このようにして、スーパーノーダル法に対してdelayed pivotsを適用して行列分割を実行する際に、効率的なメモリ１０２の使用が可能となる。
＜具体的なＬＵ分解処理手順＞
構造的に対称な実行列のスパース行列をＬＵ分解することを考える。分解対象の行列の非ゼロ要素は対称な位置にある。そのため、各ノードの分解での依存関係はelimination treeで表現される。
【０１１１】
非対称な実行列に対してもＡ⊆Ａ＋ＡTなる対称行列を考えて分解での依存関係を考えることができる。そのため構造的に対称な実行列で考えることにする。
解析部１２０は、親子関係にあるノードで非ゼロパターンが近いものをまとめてスーパーノードとする。解析部１２０は、スーパーノードに関して、スーパーノードの依存関係を表すelimination treeを生成し、supernodal elimination treeとする。さらに解析部１２０は、elimination treeと同じように、スーパーノードにもdepth first searchを順にpostorderを振る。ここで、分解された行列Ｌについて考える。このとき、スーパーノードに対応する行の集まりで、非ゼロ要素が存在する部分に対応するスーパーノードは、row subtreeに対応したsupernodal row subtreeを形成する。
【０１１２】
解析部１２０は、スーパーノードを構成しているノード及びそれらのノードより大きなノード番号を持つ非ゼロ要素を格納するために、column panelとrow panelを利用する。row panelは転置した形で行を格納する領域である。格納する非ゼロ要素の番号を知るために、インデックスリストとしての１次元配列が、column panelとrow panelとのそれぞれに設けられている。
【０１１３】
delayed pivotsを考慮せずに各ノードを分解する場合、ノードの分解は以下の順で行われる。
［１］解析部１２０は、postorder順にスーパーノードを取り出し［２］と［３］を繰り返す。
［２］解析部１２０は、選ばれたスーパーノード内の要素の値を、supernodal row subtreeからの寄与に応じて更新する。
［３］解析部１２０は、更新が終わったら、このスーパーノードのcolumn panelをＬＵ分解する。
【０１１４】
なお、行列Ｌのcolumn panelには行列Ｕの対角ブロック部分も含まれている。解析部１２０は、この行列Ｌのcolumn panelのＬＵ分解の対角ブロック部分の分解結果を使って、更新した行列Ｕのrow panelを更新する。
【０１１５】
delayed pivotsを組み込んだときの処理は以下のようになる。
解析部１２０は、スーパーノードのＬＵ分解に際して、ピボットを決定する。なお解析部１２０は、ピボットを、スーパーノードの対角ブロック内の要素からのみ決定する。例えば、対角pivotingを取ることを考える。対角pivotingは、スーパーノードの対角要素の最も大きなものを選んでpivotとする。このようなpivotの絶対値が所定値より小さいかまたはゼロになったとき、pivotが見つけられなかったことになる。
【０１１６】
このようなとき、解析部１２０は、ＬＵ分解ができたところまでの結果を使って、ピボットを見つけられなかった残りの部分の列・行の更新を行う。その後、解析部１２０は、対称置換を行って、現在分解中のスーパーノードの親のスーパーノードの後ろに、ピボットを見つけられなかった部分の列・行を移動する。なお、解析部１２０は、親スーパーノードに、移動のためのスペースを予め設けておく。
【０１１７】
このようにスーパーノードの対角ブロックの中でピボットを取ると、入れ替えの結果、対角ブロックに対応する部分の行インデックス、列インデックスも入れ替わる。そこで解析部１２０は、行と列との入れ替え結果を、インデックスに反映する。すなわち解析部１２０は、delayed pivotsを移動するとき、対応するインデックスも同じく移動する。移動した時点では、親のスーパーノードへ移動したdelayed pivotを含むcolumn panelの対角ブロック以下の部分の、子のスーパーノードの非ゼロ要素のパターンを表すインデックスリストは、親のスーパーノードの非ゼロ要素のパターンと共通となる。
【０１１８】
解析部１２０は、postorder順に行うスーパーノードのＬＵ分解を行う。スーパーノードの分解では、解析部１２０は、まずcolumn panelとrow panelの更新を行う。このとき解析部１２０は、このスーパーノードのsupernodal row subtreeの要素のスーパーノードからの寄与を計算し、要素の値を更新する。なお、supernodal row subtreeを構成するスーパーノードはdelayed pivotsの移動でブロックの幅（列数・行数）が変わるが、親のスーパーノードのpanelに更新をかけるノード番号の範囲に関する変更はない。
【０１１９】
解析部１２０は、構造的に対称な実行列の場合の対象となるスーパーノードの更新の計算を、以下の手順で行う。
［１］解析部１２０は、分解の対象となるスーパーノードのsupernodal row subtreeから構成するスーパーノードがなくなるまで取り出し、以下の［１．１］、［１．２］の処理を繰り返す。このとき解析部１２０は、supernodal row subtreeの構成要素のセカンダリスーパーノードはsupernodal row subtreeの構成要素として処理する。
［１．１］解析部１２０は、column panelを更新する。この処理の詳細は、後述のcolumn panelの更新手順（ａ．１〜ａ．４）に従う。
［１．２］解析部１２０は、１．１と同様にrow panelを更新する。
【０１２０】
［２］解析部１２０は、delayed pivotsの移動とＬＵ分解処理を行う。
［２．１］解析部１２０は、ＬＵ分解の対象となるスーパーノードの子からのdelayed pivotsの総数numdpを計算する。そして解析部１２０は、numdpの中でスペースに入る分のdelayed pivotsをスペースに入れる。
［２．２］解析部１２０は、ＬＵ分解をする。
［２．３］解析部１２０は、分解の対象となるスーパーノードのＬＵ分解でdelayed pivotsが生じたら、分解結果からdelayed pivotsに更新をかける。例えばdelayed pivotsとなったノードの列及び行の移動などが行われる。
［２．４．１］解析部１２０は、numdpのスペースに入りきらなかった分のrestdpのdelayed pivotsがあれば、それと［２．３］で生じたdelayed pivotsを合わせたものが入るセカンダリスーパーノードを確保する。そして解析部１２０は、該当するdelayed pivotsをセカンダリスーパーノードに移動する。その後、解析部１２０は、分解の対象となるスーパーノードのＬＵ分解結果を使ってrestdp部分を更新する。
［２．４．２］解析部１２０は、セカンダリスーパーノードをＬＵ分解する。この結果delayed pivotsが生じたら、解析部１２０は、分解結果からdelayed pivotsに更新をかける。例えばdelayed pivotsとなったノードの列及び行の移動などが行われる。
【０１２１】
以上が、スーパーノードの更新の計算の手順である。以下column panelの更新手順（ａ．１〜ａ．４）について詳細に説明する。
［ａ．１］解析部１２０は、スーパーノードに属するノードの行から対角ブロック部分を除いたものを転置してrow panelとする。解析部１２０は、row panelを２次元配列rpanel（nr1，nr2）に格納する。
［ａ．２］解析部１２０は、row panelに格納されている列で、更新するスーパーノードの列の番号を持つものを作業行列Ｂ（nr1，len）にコピーする。lenはこの列の本数である。解析部１２０は、column panelで更新するノードの番号を持つ先頭ntopを見つける。column panelは、２次元配列cpanel（nc1，nc2）で表される。
［ａ．３］解析部１２０は、行列Ｃ（nc1-ntop+1,len）←cpanel（ntop−nc1，nc2）×Ｂ（nr1，len）を計算する。ここでスーパーノードの大きさnc2とnr1は等しい。
［ａ．４］解析部１２０は、行列Ｃの各列を更新する列番号と同じ番号を持つスーパーノードの列で同じ行番号の要素から差し引いて更新をかける。
【０１２２】
ここで、構造的に対称な実行列のrow panelの更新は、row panelとcolumn panelを入れ替えて計算できる。
以下、余分なスペースを割り付けてdelayed pivotsを実施してＬＵ分解を行う処理の詳細を、図１７〜図２５のフローチャートを参照して説明する。
【０１２３】
図１７は、ＬＵ分解処理の手順を示すフローチャートである。以下、図１７に示す処理をステップ番号に沿って説明する。
［ステップＳ１１１］解析部１２０は、nporderに１を設定する。nporderは、処理対象のスーパーノードを示す番号である。すなわち分解対象の行列に含まれるスーパーノードには、１から昇順の識別番号が付与されている。そして、nporderで示される識別番号を有するスーパーノードが、処理の対象となる。
【０１２４】
［ステップＳ１１２］解析部１２０は、サブルーチンrsupdateを呼び出す。そして解析部１２０は、サブルーチンrsupdateを実行する。
図１８は、サブルーチンrsupdateの処理を示すフローチャートである。解析部１２０は、nporderのスーパーノードのsupernodal row subtreeにおけるnporder以外のスーパーノードからの寄与によりcolumn panelとrow panelを更新する。解析部１２０は、createされたノードがあればその分も含めて、column panelとrow panelを更新する（ステップＳ１１２ａ）。その後、ＬＵ分解処理（図１７）に戻る。
【０１２５】
［ステップＳ１１３］解析部１２０は、サブルーチンdpcountを呼び出す。そして解析部１２０は、サブルーチンdpcountを実行する。
図１９は、サブルーチンdpcountの処理を示すフローチャートである。解析部１２０は、nporderの子のスーパーノードからnporderのスーパーノードへ移動するdelayed pivotsの本数の合計を計算する（ステップＳ１１３ａ）。その後、ＬＵ分解処理（図１７）に戻る。
【０１２６】
［ステップＳ１１４］解析部１２０は、サブルーチンmvtonporderを呼び出す。そして解析部１２０は、サブルーチンmvtonporderを実行する。
図２０は、サブルーチンmvtonporderの処理を示すフローチャートである。解析部１２０は、nporderの子のスーパーノードから移動するdelayed pivotsでnporderに入る分を、nporderに移動する（ステップＳ１１４ａ）。その後、ＬＵ分解処理（図１７）に戻る。
【０１２７】
［ステップＳ１１５］解析部１２０は、サブルーチンcpupdateを呼び出す。そして解析部１２０は、サブルーチンcpupdateを実行する。サブルーチンcpupdateは、nporderのスーパーノードのcolumn panelを、delayed pivotsの分も含めて更新する処理である。
【０１２８】
図２１は、サブルーチンcpupdateの処理を示すフローチャートである。以下、図２１に示す処理をステップ番号に沿って説明する。
［ステップＳ１１５ａ］解析部１２０は、nporderのスーパーノードのcolumn panelを、スペースに移動したdelayed pivotsも含めてＬＵ分解する。
【０１２９】
［ステップＳ１１５ｂ］解析部１２０は、ＬＵ分解により、delayed pivotsとなるノードが発生したか否かを判断する。delayed pivotsとなるノードが発生した場合、処理がステップＳ１１５ｃに進められる。delayed pivotsとなるノードが発生していない場合、サブルーチンcpupdateの処理が終了する。
【０１３０】
［ステップＳ１１５ｃ］解析部１２０は、ＬＵ分解の結果に基づいて、delayed pivotsとなるノードを更新する。その後、ＬＵ分解処理（図１７）に戻る。
［ステップＳ１１６］解析部１２０は、サブルーチンrpupdateを呼び出す。そして解析部１２０は、サブルーチンrpupdateを実行する。サブルーチンrpupdateは、nporderのスーパーノードのrow panelを、delayed pivotsの分も含めて更新する処理である。
【０１３１】
図２２は、サブルーチンrpupdateの処理を示すフローチャートである。以下、図２２に示す処理をステップ番号に沿って説明する。
［ステップＳ１１６ａ］解析部１２０は、nporderのスーパーノードのrow panelを、column panelのＬＵ分解結果を使って更新する。
【０１３２】
［ステップＳ１１６ｂ］解析部１２０は、ＬＵ分解により、delayed pivotsとなるノードが発生したか否かを判断する。delayed pivotsとなるノードが発生した場合、処理がステップＳ１１６ｃに進められる。delayed pivotsとなるノードが発生していない場合、サブルーチンrpupdateの処理が終了する。
【０１３３】
［ステップＳ１１６ｃ］解析部１２０は、ＬＵ分解結果に基づいて、delayed pivotsとなるノードを更新する。その後、ＬＵ分解処理（図１７）に戻る。
［ステップＳ１１７］解析部１２０は、nporderのスーパーノードのスペースに入りきらずに残ったdelayed pivotsがあるか否かを判断する。残ったdelayed pivotsがある場合、処理がステップＳ１１８に進められる。残ったdelayed pivotsがなければ、処理がステップＳ１２１に進められる。
【０１３４】
［ステップＳ１１８］解析部１２０は、サブルーチンcreatemvを呼び出す。そして解析部１２０は、サブルーチンcreatemvを実行する。サブルーチンcreatemvは、nporderのスーパーノードに入りきらなかった残りのdelayed pivotsに、nporderのＬＵ分解で生じたdelayed pivotsを加えた分のスーパーノードを、作成する処理である。
【０１３５】
図２３は、サブルーチンcreatemvの処理を示すフローチャートである。以下、図２３に示す処理をステップ番号に沿って説明する。
［ステップＳ１１８ａ］解析部１２０は、セカンダリスーパーノードを作成する。
【０１３６】
［ステップＳ１１８ｂ］解析部１２０は、nporderのスーパーノードへ移動するべきdelayed pivotsの残りを、セカンダリスーパーノードのcolumn panel, row panelの前半に移動する。また解析部１２０は、nporderのスーパーノードに含まれるdelayed pivotsを、セカンダリスーパーノードのcolumn panel, row panelの後半に移動する。その後、ＬＵ分解処理（図１７）に戻る。
【０１３７】
［ステップＳ１１９］解析部１２０は、サブルーチンcpupdatenewを呼び出す。そして、解析部１２０は、サブルーチンcpupdatenewを実行する。サブルーチンcpupdatenewは、nporderのスーパーノードで分解された結果で、column panelの残りのdelayed pivots部分を更新したあと、全体をＬＵ分解する処理である。さらにdelayed pivotsが生じたらその部分の更新も行われる。
【０１３８】
図２４は、サブルーチンcpupdatenewの処理を示すフローチャートである。以下、図２４に示す処理をステップ番号に沿って説明する。
［ステップＳ１１９ａ］解析部１２０は、nporderのスーパーノードのＬＵ分解結果を使って、column panelの前半を更新する。
【０１３９】
［ステップＳ１１９ｂ］解析部１２０は、セカンダリスーパーノードのcolumn panelを、nporderのスーパーノードから移動したdelayed pivotsも含めてＬＵ分解する。
［ステップＳ１１９ｃ］解析部１２０は、ＬＵ分解により、delayed pivotsとなるノードが発生したか否かを判断する。delayed pivotsとなるノードが発生した場合、処理がステップＳ１１９ｄに進められる。delayed pivotsとなるノードが発生していない場合、cpupdatenewの処理が終了する。
【０１４０】
［ステップＳ１１９ｄ］解析部１２０は、セカンダリスーパーノードのＬＵ分解の結果に基づいて、column panelのdelayed pivotsとなるノードの部分を更新する。その後、ＬＵ分解処理（図１７）に戻る。
【０１４１】
［ステップＳ１２０］解析部１２０は、サブルーチンrpupdatenewを呼び出す。そして解析部１２０は、サブルーチンrpupdatenewを実行する。サブルーチンrpupdatenewは、nporderで分解された結果でrow panelを更新し、column panelの分解結果で更新する処理である。さらにdelayed pivotsが生じたらその部分の更新も行われる。
【０１４２】
図２５は、サブルーチンrpupdatenewの処理を示すフローチャートである。以下、図２４に示す処理をステップ番号に沿って説明する。
［ステップＳ１２０ａ］解析部１２０は、セカンダリスーパーノードのrow panelを、nporderのスーパーノードのＬＵ分解結果を使って更新する。
【０１４３】
［ステップＳ１２０ｂ］解析部１２０は、セカンダリスーパーノードのcolumn panelのＬＵ分解結果を使って、セカンダリスーパーノードのrow panelを更新する。
［ステップＳ１２０ｃ］解析部１２０は、ＬＵ分解により、delayed pivotsとなるノードが発生したか否かを判断する。delayed pivotsとなるノードが発生した場合、処理がステップＳ１２０ｄに進められる。delayed pivotsとなるノードが発生していない場合、rpupdatenewの処理が終了する。
【０１４４】
［ステップＳ１２０ｄ］解析部１２０は、セカンダリスーパーノードのＬＵ分解の結果に基づいて、セカンダリスーパーノードのrow panelのdelayed pivotsとなるノードの部分を更新する。その後、ＬＵ分解処理（図１７）に戻る。
【０１４５】
［ステップＳ１２１］解析部１２０は、nporderに１を加算する。
［ステップＳ１２２］解析部１２０は、nporderの値が、全スーパーノード数より大きいか否かを判断する。nporderの値が全スーパーノード数を超えた場合、ＬＵ分割処理が終了する。nporderの値が全スーパーノード数以下であれば、処理がステップＳ１１２に進められる。
【０１４６】
以上のような手順で繰り返しサブルーチンの呼び出しを行っていくことで、構造的に対象な行列を、delayed pivotsを用いて効率的にＬＵ分解することができる。
＜スパースな不定値対称行列のＬＤＬT分解＞
次に、スパースな不定値対称行列のＬＤＬT分解について説明する。不定値対称行列に対しては、対称性からcolumn panelについてだけ考えることになる。
【０１４７】
図２６は、スパースな不定値対称行列のＬＤＬT分解の処理手順を示すフローチャートである。以下、図２６に示す処理をステップ番号に沿って説明する。
［ステップＳ２０１］解析部１２０は、column panelに追加するスペースの大きさの指定入力を受け付ける。例えばユーザがキーボード２２などの入力デバイスを介して入力した、スペースの大きさを示す値を、解析部１２０が取得する。
【０１４８】
［ステップＳ２０２］解析部１２０は、シンボリック分解を行う。このシンボリック分解により、elimination treeの生成、スーパーノードの検出、supernodal row subtreeの生成などが行われる。
【０１４９】
［ステップＳ２０３］解析部１２０は、column panelのスペースを含んだ大きさを計算する。
［ステップＳ２０４］解析部１２０は、column panel用の領域、及びセカンダリスーパーノード用の領域を割り当てるためのメモリプールをメモリ内に確保する。
【０１５０】
［ステップＳ２０５］解析部１２０は、ＬＤＬT分解を行う。ＬＤＬT分解では、解析部１２０は、delayed pivotsの移動用に設けたスペースを利用する。また移動するdelayed pivotsの本数がスペースの大きさを超えた場合、解析部１２０は、セカンダリスーパーノードを作成して、delayed pivotsに対応する行を格納する。
【０１５１】
このような手順でＬＤＬT分解が進められる。以下、不定値対称行列のＬＤＬT分解を実施する際の、メモリ１０２の有効活用方法について詳細に説明する。
図２７は、スパースな不定値対称行列のＬＤＬT分解におけるパネル構造の一例を示す図である。図２７には、プライマリスーパーノード領域６１のパネル構造とセカンダリスーパーノード領域６２のパネル構造とが示されている。解析部１２０は、スペース６３を、列をまとめて格納するcolumn panelに設ける。図２７に示すように、スパースな不定値対称行列に関するスーパーノードのパネルについては、下三角行列Ｌに関する部分のみ格納すればよい。
【０１５２】
delayed pivotsの移動についても、下三角行列Ｌに関する部分のみ移動すればよい。
図２８は、不定値対称行列のＬＤＬT分解におけるdelayed pivotsの移動例を示す図である。移動される部分は、図２８の網掛けで示した部分になる。すなわち、delayed pivotsの移動は、元の行列Ａ0のうちの、ｉ行以降及びｉ列以降の部分行列Ａk内に制限されている。
【０１５３】
このように、不定値対称行列に関してはcolumn panelのみ考えればよい。
＜具体的なＬＤＬT分解処理手順＞
不定値対称行列においても非ゼロ要素は対称な位置にあるため、各ノードの分解での依存関係はelimination treeで表現される。解析部１２０は、elimination treeにおいて親子関係にあるノードで非ゼロパターンが近いものをまとめてスーパーノードとする。このようなスーパーノードに関してもスーパーノードの依存関係を表すelimination treeが構成できる。supernodal elimination treeと呼ぶ。解析部１２０は、elimination treeと同じようにスーパーノードにもpostorderを振る。ここで、分解された下三角行列Ｌについて考える。スーパーノードに対応する行の集まりで、非ゼロ要素が存在する部分に対応するスーパーノードは、row subtreeに対応したsupernodal row subtreeを形成する。
【０１５４】
解析部１２０は、各ノードの分解を、以下の手順で行う。
［１］解析部１２０は、postorder順にスーパーノードを取り出し［２］と［３］を繰り返す。
［２］解析部１２０は、選ばれたスーパーノードのsupernodal row subtreeからの寄与を更新する。
［３］解析部１２０は、更新が終わったら、このスーパーノードのcolumn panelをＬＤＬT分解する。
【０１５５】
解析部１２０は、不定値対称行列の場合の対象となるスーパーノードの更新の計算は以下の手順で行う。
［１］解析部１２０は、分解の対象となるスーパーノードのsupernodal row subtreeから構成するスーパーノードがなくなるまでスーパーノードを取り出し、以下［１．１］、［１．２］の処理を繰り返す。このとき解析部１２０は、supernodal row subtreeの構成要素のセカンダリスーパーノードは、supernodal row subtreeの構成要素として処理する。
［１．１］解析部１２０は、column panelを更新する。この処理の詳細は、後述のcolumn panelの更新手順（ａ．１〜ａ．４）に従う。
［１．２］解析部１２０は、［１．１］と同様にrow panelを更新する。
【０１５６】
［２］解析部１２０は、delayed pivotsの移動とＬＤＬT分解処理を行う。
［２．１］解析部１２０は、ＬＤＬT分解の対象となるスーパーノードの子からのdelayed pivotsの総数numdpを計算する。そして解析部１２０は、numdpの中でスペースに入る分のdelayed pivotsをスペースに入れる。
［２．２］解析部１２０は、ＬＤＬT分解をする。
［２．３］解析部１２０は、分解の対象となるスーパーノードのＬＤＬT分解でdelayed pivotsが生じたら、分解結果からdelayed pivotsに更新をかける。
［２．４．１］解析部１２０は、numdpのスペースに入りきらなかった分のrestdpがあれば、それと［２．３］で生じたdelayed pivotsを合わせたものが入るセカンダリスーパーノードを確保する。そして解析部１２０は、該当するdelayed pivotsをセカンダリスーパーノードに移動する。その後、解析部１２０は、分解の対象となるスーパーノードのＬＤＬT分解結果を使ってrestdp部分を更新する。
［２．４．２］解析部１２０は、セカンダリスーパーノードをＬＤＬT分解する。この結果delayed pivotsが生じたら、解析部１２０は、分解結果からdelayed pivotsに更新をかける。
【０１５７】
以上が、スーパーノードの更新の計算の手順である。以下column panelの更新手順（ａ．１〜ａ．４）について詳細に説明する。
column panelの更新の計算は以下の手順で行う。
［ａ．１］解析部１２０は、supernodal row subtreeからスーパーノードが終わるまで、スーパーノードを取り出し、以下の処理を繰り返す。
［ａ．２］対称性より分解行列ＬTに相当するものはcolumn panel（nc1，nc2）を転置したものであり、これがrow panelに相当する。rowを束ねたものを転置したものが、column panelに格納されている。解析部１２０は、転置されたcolumn panelに格納されている列で、更新するスーパーノードの列の番号を持つものを作業行列Ｂ（nc2，len）にコピーする。lenはこの列の本数である。解析部１２０は、column panelで更新するノードの番号を持つ先頭ntopを見つける。column panelは、２次配列cpanel（nc1，nc2）で表される。
［ａ．３］解析部１２０は、行列Ｃ（nc1−ntop＋1，len）←cpanel（ntop−nc1，nc2）×Ｄ×Ｂ（nc2，len）を計算する。
［ａ．４］解析部１２０は、行列Ｃの各列を更新する列番号と同じ番号を持つスーパーノードの列で同じ行番号の要素から差し引いて更新をかける。
【０１５８】
ここで、delayed pivotsを移動してから１×１、１×２の対角ブロックを持つＬＤＬT分解を行う。子のスーパーノードからのdelayed pivotsの移動で、スペースが不足するときは、解析部１２０は、セカンダリスーパーノードを作成する。１×１、２×２の対角ブロックをＤと表す。解析部１２０は、この対角ブロックを、column panelの対角及び副対角に格納する。
【０１５９】
不定値対称行列の分解では、解析部１２０は、１×１または２×２の小行列になるdiagonal pivotを利用する。ここで行列Ｐを１×１または２×２の小行列として、行列Ａを以下のように分解する。これを再帰的に続ける。
【０１６０】
【数１】
0006610350.tif 000002
【０１６１】
これをcolumn panelに適用すると、column panelはＬＤＬT分解される。行列Ｄは１×１、２×２の対称な対角ブロック行列になる。対角要素が１でＤが２×２のブロックとなるところは、副対角要素の位置が０となる。
【０１６２】
図２９は、ＬＤＬT分解処理の手順を示すフローチャートである。以下、図２９に示す処理をステップ番号に沿って説明する。
［ステップＳ２１１］解析部１２０は、nporderに１を設定する。
【０１６３】
［ステップＳ２１２］解析部１２０は、サブルーチンsymrsupdateを呼び出す。そして解析部１２０は、サブルーチンsymrsupdateを実行する。
図３０は、サブルーチンsymrsupdateの処理を示すフローチャートである。解析部１２０は、nporderのスーパーノードのsupernodal row subtreeにおけるnporder以外のスーパーノードからの寄与によりcolumn panelを更新する。解析部１２０は、createされたノードがあればその分も含めて、column panelを更新する（ステップＳ２１２ａ）。その後、ＬＤＬT分解処理（図２９）に戻る。
【０１６４】
［ステップＳ２１３］解析部１２０は、サブルーチンsymdpcountを呼び出す。そして解析部１２０は、サブルーチンsymdpcountを実行する。
図３１は、サブルーチンsymdpcountの処理を示すフローチャートである。解析部１２０は、nporderの子のスーパーノードからnporderのスーパーノードへ移動するdelayed pivotsの本数の合計を計算する（ステップＳ２１３ａ）。その後、ＬＤＬT分解処理（図２９）に戻る。
【０１６５】
［ステップＳ２１４］解析部１２０は、サブルーチンsymmvtonporderを呼び出す。そして解析部１２０は、サブルーチンsymmvtonporderを実行する。
図３２は、サブルーチンsymmvtonporderの処理を示すフローチャートである。解析部１２０は、nporderの子のスーパーノードから移動するdelayed pivotsでnporderに入る分を、nporderに移動する（ステップＳ２１４ａ）。その後、ＬＤＬT分解処理（図２９）に戻る。
【０１６６】
［ステップＳ２１５］解析部１２０は、サブルーチンsymcpupdateを呼び出す。そして解析部１２０は、サブルーチンsymcpupdateを実行する。サブルーチンsymcpupdateは、nporderのスーパーノードのcolumn panelを、delayed pivotsの分も含めて更新する処理である。
【０１６７】
図３３は、サブルーチンsymcpupdateの処理を示すフローチャートである。以下、図３３に示す処理をステップ番号に沿って説明する。
［ステップＳ２１５ａ］解析部１２０は、nporderのスーパーノードのcolumn panelを、スペースに移動したdelayed pivotsも含めてＬＤＬT分解する。
【０１６８】
［ステップＳ２１５ｂ］解析部１２０は、ＬＤＬT分解により、delayed pivotsとなるノードが発生したか否かを判断する。delayed pivotsとなるノードが発生した場合、処理がステップＳ２１５ｃに進められる。delayed pivotsとなるノードが発生していない場合、サブルーチンsymcpupdateの処理が終了する。
【０１６９】
［ステップＳ２１５ｃ］解析部１２０は、ＬＤＬT分解の結果に基づいて、delayed pivotsとなるノードを更新する。その後、ＬＤＬT分解処理（図２９）に戻る。
［ステップＳ２１６］解析部１２０は、nporderのスーパーノードのスペースに入りきらずに残ったdelayed pivotsがあるか否かを判断する。残ったdelayed pivotsがある場合、処理がステップＳ２１７に進められる。残ったdelayed pivotsがなければ、処理がステップＳ２１９に進められる。
【０１７０】
［ステップＳ２１７］解析部１２０は、サブルーチンsymcreatemvを呼び出す。そして解析部１２０は、サブルーチンsymcreatemvを実行する。サブルーチンsymcreatemvは、nporderのスーパーノードに入りきらなかった残りのdelayed pivotsに、nporderのＬＤＬT分解で生じたdelayed pivotsを加えた分のスーパーノードを、作成する処理である。
【０１７１】
図３４は、サブルーチンsymcreatemvの処理を示すフローチャートである。以下、図３４に示す処理をステップ番号に沿って説明する。
［ステップＳ２１７ａ］解析部１２０は、セカンダリスーパーノードを作成する。
【０１７２】
［ステップＳ２１７ｂ］解析部１２０は、nporderのスーパーノードへ移動するべきdelayed pivotsの残りを、セカンダリスーパーノードのcolumn panelの前半に移動する。また解析部１２０は、nporderのスーパーノードに含まれるdelayed pivotsを、セカンダリスーパーノードのcolumn panelの後半に移動する。その後、ＬＤＬT分解処理（図２９）に戻る。
【０１７３】
［ステップＳ２１８］解析部１２０は、サブルーチンsymcpupdatenewを呼び出す。そして、解析部１２０は、サブルーチンsymcpupdatenewを実行する。サブルーチンsymcpupdatenewは、nporderのスーパーノードで分解された結果で、column panelの残りのdelayed pivots部分を更新したあと、全体をＬＤＬT分解する処理である。さらにdelayed pivotsが生じたらその部分の更新も行われる。
【０１７４】
図３５は、サブルーチンsymcpupdatenewの処理を示すフローチャートである。以下、図３５に示す処理をステップ番号に沿って説明する。
［ステップＳ２１８ａ］解析部１２０は、nporderのスーパーノードのＬＤＬT分解結果を使って、column panelの前半を更新する。
【０１７５】
［ステップＳ２１８ｂ］解析部１２０は、セカンダリスーパーノードのcolumn panelを、nporderのスーパーノードから移動したdelayed pivotsも含めてＬＤＬT分解する。
［ステップＳ２１８ｃ］解析部１２０は、ＬＤＬT分解により、delayed pivotsとなるノードが発生したか否かを判断する。delayed pivotsとなるノードが発生した場合、処理がステップＳ２１８ｄに進められる。delayed pivotsとなるノードが発生していない場合、symcpupdatenewの処理が終了する。
【０１７６】
［ステップＳ２１８ｄ］解析部１２０は、セカンダリスーパーノードのＬＤＬT分解の結果に基づいて、column panelのdelayed pivotsとなるノードの部分を更新する。その後、ＬＤＬT分解処理（図２９）に戻る。
【０１７７】
［ステップＳ２１９］解析部１２０は、nporderに１を加算する。
［ステップＳ２２０］解析部１２０は、nporderの値が、全スーパーノード数より大きいか否かを判断する。nporderの値が全スーパーノード数を超えた場合、ＬＤＬT分割処理が終了する。nporderの値が全スーパーノード数以下であれば、処理がステップＳ２１２に進められる。
【０１７８】
以上のような手順で繰り返しサブルーチンの呼び出しを行っていくことで、構造的に対象な行列を、delayed pivotsを用いて効率的にＬＤＬT分解することができる。
［その他の実施の形態］
以上、実施の形態を例示したが、実施の形態で示した各部の構成は同様の機能を有する他のものに置換することができる。また、他の任意の構成物や工程が付加されてもよい。さらに、前述した実施の形態のうちの任意の２以上の構成（特徴）を組み合わせたものであってもよい。
【符号の説明】
【０１７９】
１  行列
２〜６  処理単位
１０  計算機
１１  記憶手段
１１ａ  第１記憶領域
１１ｂ  第２記憶領域
１２  親子関係決定手段
１３  第１記憶領域確保手段
１４  第１分解手段
１５  第２記憶領域確保手段
１６  第２分解手段
【図１】
0006610350.tif 000003
【図２】
0006610350.tif 000004
【図３】
0006610350.tif 000005
【図４】
0006610350.tif 000006
【図５】
0006610350.tif 000007
【図６】
0006610350.tif 000008
【図７】
0006610350.tif 000009
【図８】
0006610350.tif 000010
【図９】
0006610350.tif 000011
【図１０】
0006610350.tif 000012
【図１１】
0006610350.tif 000013
【図１２】
0006610350.tif 000014
【図１３】
0006610350.tif 000015
【図１４】
0006610350.tif 000016
【図１５】
0006610350.tif 000017
【図１６】
0006610350.tif 000018
【図１７】
0006610350.tif 000019
【図１８】
0006610350.tif 000020
【図１９】
0006610350.tif 000021
【図２０】
0006610350.tif 000022
【図２１】
0006610350.tif 000023
【図２２】
0006610350.tif 000024
【図２３】
0006610350.tif 000025
【図２４】
0006610350.tif 000026
【図２５】
0006610350.tif 000027
【図２６】
0006610350.tif 000028
【図２７】
0006610350.tif 000029
【図２８】
0006610350.tif 000030
【図２９】
0006610350.tif 000031
【図３０】
0006610350.tif 000032
【図３１】
0006610350.tif 000033
【図３２】
0006610350.tif 000034
【図３３】
0006610350.tif 000035
【図３４】
0006610350.tif 000036
【図３５】
0006610350.tif 000037