JP6384331B2 - 情報処理装置、情報処理方法、および情報処理プログラム - Google Patents

Description

本発明は、情報処理装置、情報処理方法、および情報処理プログラムに関する。

従来、連立１次方程式の直接解法の一つに、連立１次方程式の係数行列になる複素非対称スパース行列をＬＵ分解することによって連立１次方程式の解を求める方法がある。関連する技術としては、例えば、連立１次方程式の解法として、前処理による高速化の効果が高い前処理付き反復解法を提供するためのものがある。また、例えば、主方程式と拘束方程式から形成される全体系方程式の収束解を適度な演算量で得るための技術がある。また、例えば、圧縮列格納法を用いたスパース行列とベクトルとの積を効率よく並列に処理するための技術がある。また、例えば、連立１次方程式の解法として、不完全ＬＵ分解前処理付きクリロフ部分空間法においてベクトル演算機能を有する並列計算機上で前処理を高速に実行するための技術がある。また、例えば、スパース行列の成分のうち、非零の成分のみを１列ごとに格納するための技術がある。

特開２０１１−１４５９９９号公報 特開２００５−１１５４９７号公報 特開２００９−１９９４３０号公報 特開２００７−１３３７１０号公報 特開２０１０−１２２８５０号公報

しかしながら、上述した従来技術では、複素非対称スパース行列をＬＵ分解するにあたって、複素非対称スパース行列をＬＵ分解した結果を格納する領域として、確保しなくてもよい領域までも確保してしまう場合がある。この場合、例えば、複素非対称スパース行列をＬＵ分解する際のメモリ使用量が増大するとともに、行わなくてもよい演算を行うことになり処理効率が低下してしまう。

１つの側面では、本発明は、複素非対称スパース行列のＬＵ分解を効率的に行う情報処理装置、情報処理方法、および情報処理プログラムを提供することを目的とする。

本発明の一側面によれば、複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する情報処理装置、情報処理方法、および情報処理プログラムが提案される。

本発明の一態様によれば、複素非対称スパース行列のＬＵ分解を効率的に行うことができるという効果を奏する。

図１は、ＬＵ分解するスパース行列Ａの非零要素のパターンの一例を示す説明図である。 図２は、下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの非零要素のパターンの一例を示す説明図である。 図３は、エリミネーションツリー３０１を示す説明図である。 図４は、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との近似された非零要素のパターンの一例を示す説明図である。 図５は、ＬＵ分解するスパース行列Ａの非零要素のパターンの別の例を示す説明図である。 図６は、下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの非零要素のパターンの別の例を示す説明図である。 図７は、エリミネーションツリー７０１を示す説明図である。 図８は、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との近似された非零要素のパターンの別の例を示す説明図である。 図９は、実施の形態にかかる情報処理装置１００のハードウェアの一例を示すブロック図である。 図１０は、情報処理装置１００の機能的構成例を示すブロック図である。 図１１は、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）との非零要素のパターンを示す説明図である。 図１２は、下三角行列Ｃ（Ｐ）の非零要素のパターンを示す説明図である。 図１３は、下三角行列Ｌ（Ｐ）の非零要素のパターンを示す説明図である。 図１４は、エリミネーションツリー１４０１を示す説明図である。 図１５は、７行目の非零要素を表現する部分木１５０１の一例を示す説明図である。 図１６は、１１行目の非零要素を表現する部分木１６０１の一例を示す説明図である。 図１７は、圧縮列格納法の一例を示す説明図である。 図１８は、実施例１にかかる算出処理手順の一例を示すフローチャートである。 図１９は、実施例２にかかる算出処理手順の一例を示すフローチャートである。 図２０は、ｃｏｌｕｍｎ ｃｏｕｎｔの計数処理手順の一例を示すフローチャート（その１）である。 図２１は、ｃｏｌｕｍｎ ｃｏｕｎｔの計数処理手順の一例を示すフローチャート（その２）である。 図２２は、ｃｏｌｕｍｎ ｃｏｕｎｔの計数処理手順の一例を示すフローチャート（その３）である。

以下に、図面を参照して、本発明にかかる情報処理装置、情報処理方法、および情報処理プログラムの実施の形態を詳細に説明する。

（情報処理方法の一実施例）
まず、図１〜図４を用いて、本実施の形態にかかる情報処理方法の一実施例について説明する。図１〜図４において、情報処理装置１００は、複素非対称スパース行列ＡをＬＵ分解した結果を格納するメモリ領域量を決定するコンピュータである。

スパース行列とは、行列の要素として零である要素を多く含む行列である。スパース行列は、疎行列とも呼ばれる。複素非対称スパース行列は、当該複素スパース行列の転置行列と一致しない行列である。ＬＵ分解とは、ある行列を、下三角行列Ｌと上三角行列Ｕとの積で表現することである。下三角行列とは、対角要素より上にある要素がすべて零である行列である。対角要素とは、行番号と列番号が一致する位置にある要素である。上三角行列とは、対角要素より下にある要素がすべて零である行列である。以下の説明では、零である要素を「零要素」と表記する場合がある。また、以下の説明では、零ではない要素を「非零要素」と表記する場合がある。

複素非対称スパース行列ＡをＬＵ分解した結果を格納する領域として、複素非対称スパース行列Ａと複素非対称スパース行列Ａの転置行列Ａ＾Ｔとの和になる対称行列ＰをＬＬ＾Ｔ分解した結果を格納するための領域を確保することが考えられる。この場合、例えば、演算装置は、対称行列ＰをＬＬ＾Ｔ分解する場合の各列の依存関係を示すエリミネーションツリー（ｅｌｉｍｉｎａｔｉｏｎ ｔｒｅｅ）を生成して、対称行列ＰをＬＬ＾Ｔ分解した結果を格納するための領域を算出することになる。

ＬＬ＾Ｔ分解とは、ＬＵ分解の一つである。ＬＬ＾Ｔ分解とは、ある対称行列を、下三角行列Ｌと、下三角行列Ｌの転置行列Ｌ＾Ｔとの積で表現することである。ＬＬ＾Ｔ分解は、コレスキー分解とも呼ばれる。エリミネーションツリーとは、対称行列Ｐの各列の依存関係を示すツリーである。エリミネーションツリーは、各列に関するノード（ｎｏｄｅ）を含むツリーである。エリミネーションツリーは、消去木とも呼ばれる。ノードは、節点とも呼ばれる。以下の説明では、ｊ列目に関するノードを「ノード［ｊ］」と表記する場合がある。また、以下の説明では、ノード［ｊ］における「ｊ」を「インデックス」と表記する場合がある。

しかしながら、この場合、複素非対称スパース行列ＡをＬＵ分解した結果を格納する領域は、格納しなくてもよい零要素を格納するための領域までも含んでしまう。すなわち、複素非対称スパース行列Ａが大きくなるにつれて、確保する領域が足りなくなって、複素非対称スパース行列ＡをＬＵ分解することができなくなることがある。また、ＬＵ分解において行わなくてもよい零要素に関する演算を演算装置が行ってしまい、ＬＵ分解にかかる処理時間の増大を招いてしまう。

特に、非対称スパース行列Ａの非対称度合いが大きくなると、複素非対称スパース行列Ａの互いに対称位置にある要素の組み合わせの多くが非零要素と零要素との組み合わせになってしまう場合がある。この場合、複素非対称スパース行列Ａから生成した対称行列Ｐにおいて、複素非対称スパース行列Ａでは零要素があった位置に対応する位置にも、非零要素が出現してしまう。このため、対称行列Ｐの非零要素が対称行列ＰをＬＬ＾Ｔ分解した結果に影響し、対称行列ＰをＬＬ＾Ｔ分解した結果に含まれる非零要素の数を、複素非対称スパース行列ＡをＬＵ分解した結果に含まれる非零要素よりも増大させてしまう。結果として、対称行列ＰをＬＬ＾Ｔ分解した結果に基づいて確保された、複素非対称スパース行列ＡをＬＵ分解した結果を格納する領域は、格納しなくてもよい零要素を格納するための領域までも含んでしまう。

そこで、本実施の形態では、複素非対称スパース行列のＬＵ分解を効率的に行うことができる情報処理方法について説明する。この情報処理方法によれば、格納しなくてもよい零要素を格納する領域を確保してしまうことを抑制し、複素非対称スパース行列のＬＵ分解を効率的に行うことができる。以下の説明では、複素非対称スパース行列Ａを「スパース行列Ａ」と表記する場合がある。

まず、情報処理装置１００は、ＬＵ分解するスパース行列Ａを取得する。ＬＵ分解するスパース行列Ａは、例えば、１１行１１列の行列である。以下の説明では、スパース行列Ａのｉ行ｊ列にある要素を「要素ａ［ｉ，ｊ］」と表記する場合がある。ここで、図１を用いて、ＬＵ分解するスパース行列Ａの非零要素のパターンの一例について説明する。

図１は、ＬＵ分解するスパース行列Ａの非零要素のパターンの一例を示す説明図である。図１の方眼１０１のｉ行ｊ列目の升目は、スパース行列Ａのｉ行ｊ列目の要素に対応し、スパース行列Ａのｉ行ｊ列目の要素が対角要素、零要素、および非零要素のいずれであるかを示す。例えば、対角要素は、「対角要素があるスパース行列Ａの行番号ｉ（＝列番号ｊ）」で示される。また、零要素は、「空白」で示される。また、非零要素は、「●」で示される。図１に示すように、スパース行列Ａの７，１０行目には、対角要素を除いて非零要素はない。

次に、情報処理装置１００は、スパース行列Ａとスパース行列Ａ＾Ｔとの和となる、スパース行列Ａを対称化した対称行列ＰをＬＬ＾Ｔ分解して、下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの積で表現する。

以下の説明では、下三角行列Ｌ（Ｐ）のｉ行ｊ列にある要素を「要素ｌ［ｉ，ｊ］」と表記する場合がある。また、下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔのｉ行ｊ列にある要素を「要素ｌｔ［ｉ，ｊ］」と表記する場合がある。ここで、図２を用いて、対称行列ＰをＬＬ＾Ｔ分解して得られる、下三角行列Ｌ（Ｐ）と、下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの非零要素のパターンの一例について説明する。

図２は、下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの非零要素のパターンの一例を示す説明図である。図２の方眼２０１は、対角要素で分割された下三角部分によって下三角行列Ｌ（Ｐ）の非零要素のパターンを示し、対角要素で分割された上三角部分によって下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔの非零要素のパターンを示す。

図２の方眼２０１のｉ行ｊ列目の升目は、ｉ＞ｊの場合には、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素に対応し、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素が対角要素、零要素、非零要素、およびフィルインのいずれであるかを示す。フィルインとは、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素ｌ［ｉ，ｊ］であって、対称行列Ｐにおいて同じ位置にあるｉ行ｊ列目の要素ｐ［ｉ，ｊ］が零要素であるのに、非零要素になってしまう要素である。例えば、フィルインは、「○」で示される。

一方で、図２の方眼２０１のｉ行ｊ列目の升目は、ｉ＜ｊの場合には、下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔのｉ行ｊ列目の要素に対応し、転置行列Ｌ（Ｐ）＾Ｔのｉ行ｊ列目の要素が対角要素、零要素、非零要素、およびフィルインのいずれであるかを示す。

ここで、対称行列Ｐの非零要素のパターンは、スパース行列Ａの非零要素のパターンを包含する。このため、スパース行列Ａの非零要素が影響して下三角行列Ｌ（Ａ）において非零要素が出現する位置は、スパース行列Ａの非零要素と同じ位置にある対称行列Ｐの非零要素が影響して下三角行列Ｌ（Ｐ）において非零要素が出現する位置と一致する。同様に、上三角行列Ｕ（Ａ）において非零要素が出現する位置は、転置行列Ｌ（Ｐ）＾Ｔにおいて非零要素が出現する位置と一致する。

一方で、スパース行列Ａにおいて零要素がある位置にも対称行列Ｐにおいては非零要素があるため、当該非零要素が影響して、下三角行列Ｌ（Ａ）において非零要素が出現しない位置にも下三角行列Ｌ（Ｐ）においては非零要素が出現する。同様に、上三角行列Ｕ（Ａ）において非零要素が出現しない位置にも転置行列Ｌ（Ｐ）＾Ｔにおいては非零要素が出現する。これらのことから、対称行列ＰをＬＬ＾Ｔ分解した場合の下三角行列Ｌ（Ｐ）や転置行列Ｌ（Ｐ）＾Ｔの非零要素のパターンは、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）や上三角行列Ｕ（Ａ）の非零要素のパターンを包含する。

次に、情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解した結果に基づいて、対称行列Ｐのエリミネーションツリーを生成する。ここで、図３を用いて、図２に示した下三角行列Ｌ（Ｐ）や下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔの非零要素のパターンに基づくエリミネーションツリーについて説明する。

図３は、エリミネーションツリー３０１を示す説明図である。情報処理装置１００は、ｉ＞ｊであるノード［ｉ］とノード［ｊ］とがある場合に、ｍｉｎ｛ｉ｜ｉ＞ｊかつｌ［ｉ，ｊ］≠０｝を満たせば、ノード［ｉ］をノード［ｊ］の親ノード（ｐａｒｅｎｔ）として、エリミネーションツリー３０１を生成する。エリミネーションツリー３０１において、ノード［ｉ］がノード［ｊ］の祖先でなければ、ＬＬ＾Ｔ分解した下三角行列Ｌ（Ｐ）のｉ列目の要素を算出する際に、ＬＬ＾Ｔ分解した下三角行列Ｌ（Ｐ）のｊ列目の要素が影響することはない。

ここでは、情報処理装置１００が、エリミネーションツリー３０１を、対称行列ＰをＬＬ＾Ｔ分解した結果に基づいて生成する場合について説明したが、これに限らない。例えば、情報処理装置１００は、エリミネーションツリー３０１を、対称行列ＰをＬＬ＾Ｔ分解しなくても、スパース行列Ａに基づいて生成することができるし、ＬＬ＾Ｔ分解する前の対称行列Ｐの下三角行列Ｃ（Ｐ）に基づいて生成することもできる。

次に、情報処理装置１００は、エリミネーションツリー３０１に基づいて、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との非零要素のパターンを特定する。ここで、下三角行列Ｌ（Ａ）のｉ行目の非零要素のパターンは、下三角行列Ｌ（Ａ）のｉ行目において対角要素がある列に関するノードを根ノード（ｒｏｏｔ）として含む、エリミネーションツリー３０１のロウサブツリー（ｒｏｗ ｓｕｂｔｒｅｅ）で、近似して表現される。同様に、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔのｉ行目の非零要素のパターンは、転置行列Ｕ（Ａ）＾Ｔのｉ行目において対角要素がある列に関するノードを根ノードとして含む、エリミネーションツリー３０１のロウサブツリーで、近似して表現される。以下の説明では、ロウサブツリーを「部分木」と表記する場合がある。

このため、情報処理装置１００は、下三角行列Ｃ（Ａ）のｉ行目において対角要素がある列に関するノードを根ノードとして特定し、対角要素以外の非零要素がある列に関するノードを葉ノード（ｌｅａｆ）として特定する。次に、情報処理装置１００は、下三角行列Ｌ（Ａ）のｉ行目の非零要素のパターンを、特定した根ノードと葉ノードとを含むエリミネーションツリー３０１の部分木によって近似する。そして、情報処理装置１００は、下三角行列Ｌ（Ａ）の各行の非零要素のパターンに基づいて、下三角行列Ｌ（Ａ）の各列にある非零要素の数を算出する。

同様に、情報処理装置１００は、上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔのｉ行目において対角要素がある列に関するノードを根ノードとして特定し、対角要素以外の非零要素がある列に関するノードを葉ノードとして特定する。次に、情報処理装置１００は、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔのｉ行目の非零要素のパターンを、特定した根ノードと葉ノードとを含むエリミネーションツリー３０１の部分木によって近似する。そして、情報処理装置１００は、転置行列Ｕ（Ａ）＾Ｔの各行の非零要素のパターンに基づいて、転置行列Ｕ（Ａ）＾Ｔの各列にある非零要素の数を算出する。

情報処理装置１００は、例えば、図１のスパース行列Ａの下三角行列Ｃ（Ａ）の７行目の非零要素のパターンから、ＬＵ分解した下三角行列Ｌ（Ａ）の７行目に葉ノードはないと特定する。そして、情報処理装置１００は、ノード［７］を含む部分木によって下三角行列Ｌ（Ａ）の７行目の非零要素のパターンを近似する。

また、情報処理装置１００は、図１のスパース行列Ａの下三角行列Ｃ（Ａ）の１０行目の非零要素のパターンから、下三角行列Ｌ（Ａ）の１０行目に葉ノードはないと特定する。そして、情報処理装置１００は、ノード［１０］を含む部分木によって下三角行列Ｌ（Ａ）の１０行目の非零要素のパターンを近似する。ここで、図４を用いて、スパース行列ＡをＬＵ分解して得られた下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との近似された非零要素のパターンの一例について説明する。

図４は、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との近似された非零要素のパターンの一例を示す説明図である。図４の方眼４０１では、対角要素で分割された下三角部分によって下三角行列Ｌ（Ａ）の近似された非零要素のパターンを示し、対角要素で分割された上三角部分によって上三角行列Ｕ（Ａ）の近似された非零要素のパターンを示す。

図４の方眼４０１のｉ行ｊ列目の升目は、ｉ＞ｊの場合には、下三角行列Ｌ（Ａ）のｉ行ｊ列目の要素に対応し、下三角行列Ｌ（Ａ）のｉ行ｊ列目の要素が対角要素、零要素、非零要素、フィルイン、および偽フィルインのいずれであるかを示す。偽フィルインとは、実際にはフィルインにならないが、非零要素のパターンを近似したためにフィルインになった要素である。例えば、偽フィルインは、「◎」で示される。

一方で、図４の方眼４０１のｉ行ｊ列目の升目は、ｉ＜ｊの場合には、上三角行列Ｕ（Ａ）のｉ行ｊ列目の要素に対応し、上三角行列Ｕ（Ａ）のｉ行ｊ列目の要素が対角要素、零要素、非零要素、フィルイン、および偽フィルインのいずれであるかを示す。

ここで、対称行列Ｐにおいて非零要素がある位置であってもスパース行列Ａにおいては非零要素がない場合がある。このため、下三角行列Ｌ（Ａ）の非零要素のパターンを、対称行列Ｐとエリミネーションツリー３０１の組み合わせから近似するよりも、スパース行列Ａとエリミネーションツリー３０１の組み合わせから近似する方が、非零要素の影響が少なくなる。結果として、実際にＬＵ分解した下三角行列Ｌ（Ａ）では非零要素にならない要素を、非零要素になると判定することが抑制される。同様に、実際にＬＵ分解した上三角行列Ｕ（Ａ）では非零要素にならない要素を、非零要素になると判定することが抑制される。

換言すれば、図４の非零要素のパターンは、対称行列ＰをＬＬ＾Ｔ分解した場合の図２に示した非零要素のパターンに包含され、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンよりもフィルインの数が少なくなることがある。図４の例では、図４の非零要素のパターンにおける７，１０行目にあるフィルインの数は、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンにおける７，１０行目にあるフィルインの数よりも少なくなる。

また、エリミネーションツリー３０１は、スパース行列ＡをＬＵ分解する場合の各列の依存関係を示したツリーではないため、実際にＬＵ分解する場合には依存関係のない列同士が依存関係のある列同士とされることがある。しかしながら、少なくとも、実際にＬＵ分解する場合に依存関係のある列同士は、依存関係のある列同士として示される。このため、スパース行列Ａの非零要素が影響して、非零要素が出現する下三角行列Ｌ（Ａ）の位置については、少なくとも、非零要素が出現する位置として判定されることになる。同様に、スパース行列Ａの非零要素が影響して、非零要素が出現する上三角行列Ｕ（Ａ）の位置については、少なくとも、非零要素が出現する位置として判定されることになる。

換言すれば、図４の非零要素のパターンは、スパース行列ＡをＬＵ分解した場合の非零要素のパターンよりもフィルインの数が多くなることがあるが、少なくともスパース行列ＡをＬＵ分解した場合の非零要素のパターンを包含することになる。

これにより、情報処理装置１００は、下三角行列Ｌ（Ａ）を格納する領域の大きさを精度よく算出することができるようになる。情報処理装置１００は、例えば、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンに基づく領域よりも小さく、かつ、実際にスパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）や上三角行列Ｕ（Ａ）を格納可能な領域の大きさを算出することができる。情報処理装置１００は、具体的には、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）を格納する領域として、下三角行列Ｌ（Ａ）の各列と上三角行列Ｕ（Ａ）の各行との非零要素があるインデックスを格納する領域と非零要素を格納する領域とを用意する。

このように、情報処理装置１００は、ＬＵ分解した結果を格納する領域の大きさを低減することができ、スパース行列Ａが大きくなってもＬＵ分解した結果を格納することができるようになる。例えば、スパース行列Ａの非対称度合いが大きく、下三角行列Ｃ（Ａ）のみに非零要素があるような場合がある。この場合であれば、情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンに基づいて領域の大きさを算出する場合に比べて、領域の大きさを約半分まで低減することができる可能性がある。

また、情報処理装置１００は、ＬＵ分解において行わなくてもよい演算を省略することができ、効率よくＬＵ分解することができる可能性がある。例えば、スパース行列Ａの非対称度合いが大きく、下三角行列Ｃ（Ａ）のみに非零要素があるような場合がある。この場合であれば、情報処理装置１００は、演算量も約半分に低減することができる可能性がある。

ここで、情報処理装置１００が、実際にＬＵ分解する場合を例に挙げて、演算量を低減することについて説明する。情報処理装置１００は、実際にＬＵ分解する際には、エリミネーションツリー３０１を深さ優先探索（ｄｅｐｔｈ ｆｉｒｓｔ ｓｅａｒｃｈ）して、ポストオーダー（ｐｏｓｔｏｒｄｅｒ）を付与する。次に、情報処理装置１００は、付与したポストオーダーの順にｌｅｆｔ ｌｏｏｋｉｎｇおよびｕｐｗａｒｄ ｌｏｏｋｉｎｇによって、下三角行列Ｌ（Ａ）の各列および上三角行列Ｕ（Ａ）の各行を更新する。

ｌｅｆｔ ｌｏｏｋｉｎｇとは、スパース行列ＡのＬＵ分解において、下三角行列Ｌ（Ａ）のｊ列にある要素を、下三角行列Ｌ（Ａ）のｊ列よりも左側にある列の要素を参照して更新することである。また、ｕｐｗａｒｄ ｌｏｏｋｉｎｇとは、スパース行列ＡのＬＵ分解において、上三角行列Ｕ（Ａ）のｉ行にある要素を、上三角行列Ｕ（Ａ）のｉ行よりも上側の行の要素を参照して更新することである。以下の説明では、下三角行列Ｌ（Ａ）のｉ行ｊ列目の要素を「ｌａ［ｉ，ｊ］」と表記する場合がある。また、以下の説明では、上三角行列Ｕ（Ａ）のｉ行ｊ列目の要素を「ｕａ［ｉ，ｊ］」と表記する場合がある。

例えば、情報処理装置１００は、上三角行列Ｕ（Ａ）の７行目の要素ｕａ［７，ｊ］を更新する場合がある。この場合には、情報処理装置１００は、下三角行列Ｌ（Ａ）の７行目の非零要素ｌａ［７，ｉ］（ｉ＜７）と、当該非零要素の列番号と同じ値を行番号として有する上三角行列Ｕ（Ａ）の要素ｕａ［ｉ，ｊ］を乗算して、ａ［７，ｊ］から減算する。しかしながら、情報処理装置１００は、下三角行列Ｌ（Ａ）の要素ｌａ［７，ｉ］（ｉ＜７）に非零要素がなければ、上三角行列Ｕ（Ａ）の７行目を更新する場合の演算を行わなくてもよいことになる。

このため、情報処理装置１００は、下三角行列Ｌ（Ａ）の非零要素のパターンが図４の例になる場合であれば、上三角行列Ｕ（Ａ）の７行目を更新する演算を省略することにより、ＬＵ分解の際に演算量を低減することができる。具体的には、情報処理装置１００は、下三角行列Ｌ（Ａ）の７行目に関して領域が確保されているか否かを判定し、確保されていなければ下三角行列Ｌ（Ａ）の７行目は非零要素であるため、上三角行列Ｕ（Ａ）の７行目についての演算を省略する。

ここで、情報処理装置１００が、上述した演算の中で、非零要素がある位置をどのように特定するのかについて説明する。情報処理装置１００は、非零要素がある位置を特定する際には、各ノードを、付与されたポストオーダーの順に辿ればよい。情報処理装置１００は、例えば、各ノードを、付与されたポストオーダーの順に辿り、当該ノードの子ノード（ｃｈｉｌｄ ｎｏｄｅ）のインデックスの集合と、当該ノードに応じてＬＵ分解する列の非零要素があるインデックスの集合の和を算出する。これにより、情報処理装置１００は、下三角行列Ｌ（Ａ）の各列の非零要素のインデックスを特定することができる。

同様に、情報処理装置１００は、各ノードを、付与されたポストオーダーの順に辿り、当該ノードの子ノードのインデックスの集合と、当該ノードに応じてＬＵ分解する行の非零要素があるインデックスの集合の和を算出する。これにより、情報処理装置１００は、上三角行列Ｕ（Ａ）の各行の非零要素のインデックスを特定することができる。

これらのことから、情報処理装置１００によれば、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを削減することができる。そして、情報処理装置１００によれば、確保しなくてもよい零要素を格納する領域を確保しないため、ＬＵ分解にかかる処理量を抑えて、処理時間の増大を防ぎ、ＬＵ分解を効率的に行うことができる。

これにより、例えば、電磁場、音響、量子力学、および回路などの解析における複素非対称スパース行列を用いた連立１次方程式を解く際に、ＬＵ分解した結果を格納する領域の大きさを削減して、大規模問題に対応することが可能となる。

（情報処理方法の他の実施例）
まず、情報処理装置１００は、図１に示す非零要素のパターンとは異なる非零要素のパターンを有する、ＬＵ分解するスパース行列Ａを取得する。ＬＵ分解するスパース行列Ａは、例えば、１１行１１列の行列である。ここで、図５を用いて、ＬＵ分解するスパース行列Ａの非零要素のパターンの別の例を示す。

図５は、ＬＵ分解するスパース行列Ａの非零要素のパターンの別の例を示す説明図である。図５の方眼５０１のｉ行ｊ列目の升目は、スパース行列Ａのｉ行ｊ列目の要素に対応し、スパース行列Ａのｉ行ｊ列目の要素が対角要素、零要素、および非零要素のいずれであるかを示す。図５に示すように、スパース行列Ａは、対称度合いが大きい行列であって、スパース行列Ａの下三角行列Ｃ（Ａ）にある非零要素が上三角行列Ｒ（Ａ）にある非零要素よりも多い。

次に、情報処理装置１００は、スパース行列Ａを対称化した対称行列Ｐを生成し、生成した対称行列ＰをＬＬ＾Ｔ分解して下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの積で表現する。ここで、図６を用いて、対称行列ＰをＬＬ＾Ｔ分解して得られた下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの非零要素のパターンの別の例について説明する。

図６は、下三角行列Ｌ（Ｐ）と下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔとの非零要素のパターンの別の例を示す説明図である。図６の方眼６０１では、対角要素で分割された下三角部分によって下三角行列Ｌ（Ｐ）の非零要素のパターンを示し、対角要素で分割された上三角部分によって下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔの非零要素のパターンを示す。

図６の方眼６０１のｉ行ｊ列目の升目は、ｉ＞ｊの場合には、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素に対応し、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素が対角要素、零要素、非零要素、およびフィルインのいずれであるかを示す。一方で、図６の方眼６０１のｉ行ｊ列目の升目は、ｉ＜ｊの場合には、下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔのｉ行ｊ列目の要素に対応し、転置行列Ｌ（Ｐ）＾Ｔのｉ行ｊ列目の要素が対角要素、零要素、非零要素、およびフィルインのいずれであるかを示す。

ここで、対称行列ＰをＬＬ＾Ｔ分解した場合の下三角行列Ｌ（Ｐ）や転置行列Ｌ（Ｐ）＾Ｔの非零要素のパターンは、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）や上三角行列Ｕ（Ａ）の非零要素のパターンを包含する。

次に、情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解した結果に基づいて、対称行列Ｐのエリミネーションツリー７０１を生成する。ここで、図７を用いて、図６に示した下三角行列Ｌ（Ｐ）や下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔの非零要素のパターンに基づくエリミネーションツリー７０１について説明する。

図７は、エリミネーションツリー７０１を示す説明図である。情報処理装置１００は、ｉ＞ｊであるノード［ｉ］とノード［ｊ］とがある場合に、ｍｉｎ｛ｉ｜ｉ＞ｊかつｌ［ｉ，ｊ］≠０｝を満たせば、ノード［ｉ］をノード［ｊ］の親ノード（親ノード）として、エリミネーションツリー７０１を生成する。

次に、情報処理装置１００は、エリミネーションツリー７０１に基づいて、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との非零要素のパターンを特定する。

例えば、情報処理装置１００は、下三角行列Ｃ（Ａ）のｉ行目において対角要素がある列に関するノードを根ノードとして特定し、対角要素以外の非零要素がある列に関するノードを葉ノード（ｌｅａｆ）として特定する。次に、情報処理装置１００は、下三角行列Ｌ（Ａ）のｉ行目の非零要素のパターンを、特定した根ノードと葉ノードとを含むエリミネーションツリー７０１の部分木によって近似する。そして、情報処理装置１００は、下三角行列Ｌ（Ａ）の各行の非零要素のパターンに基づいて、下三角行列Ｌ（Ａ）の各列にある非零要素の数を算出する。

同様に、情報処理装置１００は、上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔのｉ行目において対角要素がある列に関するノードを根ノードとして特定し、対角要素以外の非零要素がある列に関するノードを葉ノードとして特定する。次に、情報処理装置１００は、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔのｉ行目の非零要素のパターンを、特定した根ノードと葉ノードとを含むエリミネーションツリー７０１の部分木によって近似する。そして、情報処理装置１００は、転置行列Ｕ（Ａ）＾Ｔの各行の非零要素のパターンに基づいて、転置行列Ｕ（Ａ）＾Ｔの各列にある非零要素の数を算出する。

情報処理装置１００は、例えば、図５のスパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの１，２，４，６〜９，１１行目の非零要素のパターンを特定する。次に、情報処理装置１００は、特定した非零要素のパターンから、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの１，２，４，６〜９，１１行目に葉ノードはないと特定する。そして、情報処理装置１００は、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの１，２，４，６〜９，１１行目に対応する部分木によって上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの１，２，４，６〜９，１１行目の非零要素のパターンを近似する。

また、情報処理装置１００は、図５のスパース行列Ａの下三角行列Ｃ（Ａ）の転置行列Ｃ（Ａ）＾Ｔの３，５，１０行目の非零要素のパターンから、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの３，５，１０行目にある葉ノードを特定する。そして、情報処理装置１００は、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの３，５，１０行目に対応する部分木によって上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの３，５，１０行目の非零要素のパターンを近似する。ここで、図８を用いて、スパース行列ＡをＬＵ分解して得られた下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との近似された非零要素のパターンの別の例について説明する。

図８は、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）との近似された非零要素のパターンの別の例を示す説明図である。図８の方眼８０１では、対角要素で分割された下三角部分によって下三角行列Ｌ（Ａ）の近似された非零要素のパターンを示し、対角要素で分割された上三角部分によって上三角行列Ｕ（Ａ）の近似された非零要素のパターンを示す。

図８の方眼８０１のｉ行ｊ列目の升目は、ｉ＞ｊの場合には、下三角行列Ｌ（Ａ）のｉ行ｊ列目の要素に対応し、下三角行列Ｌ（Ａ）のｉ行ｊ列目の要素が対角要素、零要素、非零要素、フィルイン、および偽フィルインのいずれであるかを示す。一方で、図８の方眼８０１のｉ行ｊ列目の升目は、ｉ＜ｊの場合には、上三角行列Ｕ（Ａ）のｉ行ｊ列目の要素に対応し、上三角行列Ｕ（Ａ）のｉ行ｊ列目の要素が対角要素、零要素、非零要素、フィルイン、および偽フィルインのいずれであるかを示す。

ここで、図８の非零要素のパターンは、スパース行列ＡをＬＵ分解した場合の非零要素のパターンよりもフィルインの数が多くなることがあるが、スパース行列ＡをＬＵ分解した場合の非零要素のパターンを包含することになる。一方で、図８の非零要素のパターンは、対称行列ＰをＬＬ＾Ｔ分解した場合の図６に示した非零要素のパターンに包含され、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンよりもフィルインの数が少なくなることがある。図８の例では、図８の非零要素のパターンにおける各行にあるフィルインの数は、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンにおける各行にあるフィルインの数よりも少なくなる。

これにより、情報処理装置１００は、下三角行列Ｌ（Ａ）を格納する領域の大きさを精度よく算出することができるようになる。情報処理装置１００は、例えば、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンに基づく領域よりも小さく、かつ、実際にスパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）や上三角行列Ｕ（Ａ）を格納可能な領域の大きさを算出することができる。情報処理装置１００は、具体的には、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）を格納する領域として、下三角行列Ｌ（Ａ）の各列と上三角行列Ｕ（Ａ）の各行との非零要素があるインデックスを格納する領域と非零要素を格納する領域とを用意する。

このように、情報処理装置１００は、ＬＵ分解した結果を格納する領域の大きさを低減することができ、スパース行列Ａが大きくなってもＬＵ分解した結果を格納することができるようになる。例えば、スパース行列Ａの非対称度合いが大きく、下三角行列Ｃ（Ａ）のみに非零要素があるような場合がある。この場合であれば、情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解した場合の非零要素のパターンに基づいて領域の大きさを算出する場合に比べて、領域の大きさを約半分に低減することができる可能性がある。

また、情報処理装置１００は、ＬＵ分解において行わなくてもよい演算を省略することができ、効率よくＬＵ分解することができる。例えば、スパース行列Ａの非対称度合いが大きく、下三角行列Ｃ（Ａ）のみに非零要素があるような場合がある。この場合であれば、情報処理装置１００は、演算量も約半分に低減することができる可能性がある。

（情報処理装置１００のハードウェア）
次に、図９を用いて、実施の形態にかかる情報処理装置１００のハードウェアの一例について説明する。

図９は、実施の形態にかかる情報処理装置１００のハードウェアの一例を示すブロック図である。図９において、情報処理装置１００は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）９０１と、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）９０２と、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）９０３と、を有する。また、情報処理装置１００は、さらに、ディスクドライブ９０４と、ディスク９０５と、インターフェース（Ｉ／Ｆ：Ｉｎｔｅｒｆａｃｅ）９０６と、を有する。

また、ＣＰＵ９０１と、ＲＯＭ９０２と、ＲＡＭ９０３と、ディスクドライブ９０４と、Ｉ／Ｆ９０６とは、バス９００によってそれぞれ接続されている。情報処理装置１００は、例えば、サーバ、ＰＣ（Ｐｅｒｓｏｎａｌ Ｃｏｍｐｕｔｅｒ）、ノートＰＣ、タブレット型ＰＣなどである。

ＣＰＵ９０１は、情報処理装置１００の全体の制御を司る。ＲＯＭ９０２は、ブートプログラム、情報処理プログラムなどの各種プログラムを記憶する。ＲＡＭ９０３は、ＣＰＵ９０１のワークエリアとして使用される。また、ＲＡＭ９０３は、各種プログラムの実行により得られたデータなどの各種データを記憶する。ディスクドライブ９０４は、ＣＰＵ９０１の制御により、ディスク９０５に対するデータのリード／ライトを制御する。ディスク９０５は、ディスクドライブ９０４の制御により書き込まれたデータを記憶する。

Ｉ／Ｆ９０６は、通信回線を通じてネットワーク９１０に接続され、このネットワーク９１０を介して他の装置に接続される。ネットワーク９１０は、例えば、ＬＡＮ（Ｌｏｃａｌ Ａｒｅａ Ｎｅｔｗｏｒｋ）、ＷＡＮ（Ｗｉｄｅ Ａｒｅａ Ｎｅｔｗｏｒｋ）、インターネットなどである。Ｉ／Ｆ９０６は、ネットワーク９１０と内部のインターフェースを司り、外部装置からのデータの入出力を制御する。Ｉ／Ｆ９０６は、例えば、モデムやＬＡＮアダプタなどである。

情報処理装置１００は、ディスクドライブ９０４とディスク９０５との代わりに、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）と半導体メモリとを有していてもよい。また、情報処理装置１００は、光ディスク、ディスプレイ、キーボード、マウス、スキャナ、およびプリンタの少なくともいずれか一つを有してもよい。

（情報処理装置１００の機能的構成例）
次に、図１０を用いて、情報処理装置１００の機能的構成例について説明する。

図１０は、情報処理装置１００の機能的構成例を示すブロック図である。情報処理装置１００は、制御部となる機能として、取得部１００１と、算出部１００２と、分解部１００３とを含む。

取得部１００１は、ＬＵ分解する行列を取得する。ＬＵ分解する行列は、例えば、図１や図５に示した非零要素のパターンを有するスパース行列Ａである。これにより、取得部１００１は、算出部１００２にスパース行列ＡをＬＵ分解した結果を格納する領域の大きさを算出させるために、算出部１００２にスパース行列Ａを入力することができる。

取得されたデータは、例えば、ＲＡＭ９０３、ディスク９０５などの記憶領域に記憶される。取得部１００１は、例えば、図９に示したＲＯＭ９０２、ＲＡＭ９０３、ディスク９０５などの記憶装置に記憶されたプログラムをＣＰＵ９０１に実行させることにより、または、Ｉ／Ｆ９０６により、その機能を実現する。

算出部１００２は、スパース行列Ａから、スパース行列Ａの対称行列Ｐのエリミネーションツリーを生成する。算出部１００２は、例えば、スパース行列Ａの対称行列Ｐの非零要素のパターンを特定して、対称行列Ｐのエリミネーションツリーを生成する。

ここで、算出部１００２は、対称行列Ｐの非零要素のパターンを特定することができれば、対称行列Ｐを生成しなくてもよい。また、算出部１００２は、例えば、スパース行列Ａの互いに対称位置にある要素ａ［ｉ，ｊ］とａ［ｊ，ｉ］とから、対称行列Ｐを生成せずに、対称行列Ｐのエリミネーションツリーを生成してもよい。対称位置とは、対角要素に対して対称な位置であり、ｉ行ｊ列目の位置とｊ行ｉ列目の位置とである。エリミネーションツリーを生成する詳細は、後述する実施例１の第１の工程〜第４の工程において説明する。

算出部１００２は、生成したエリミネーションツリーに基づいて、スパース行列Ａの下三角行列Ｃ（Ａ）の各行の部分木を抽出する。各行の部分木とは、各行に対応するロウサブツリーである。各行の部分木は、当該行の行番号と同じ値をインデックスとして有するノードを根ノードとし、当該行において非零要素がある列の列番号と同じ値をインデックスとして有するノードを葉ノードとする、エリミネーションツリーの部分木である。そして、算出部１００２は、抽出した下三角行列Ｃ（Ａ）の各行の部分木のうち、エリミネーションツリーの各ノードを含む部分木の数を算出する。

算出部１００２は、例えば、スパース行列Ａの下三角行列Ｃ（Ａ）の非零要素のパターンを特定して、下三角行列Ｃ（Ａ）の各行の部分木を抽出する。そして、算出部１００２は、エリミネーションツリーのノードごとに、当該ノードを含む部分木がいくつあるかを計数し、当該ノードを含む部分木の数を算出する。

ここで、算出部１００２は、下三角行列Ｃ（Ａ）の非零要素のパターンを特定することができれば、下三角行列Ｃ（Ａ）を生成しなくてもよい。部分木の数を算出する詳細は、後述する実施例１の第５の工程において説明する。

算出部１００２は、生成したエリミネーションツリーに基づいて、スパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの各行の部分木を抽出する。そして、算出部１００２は、抽出した転置行列Ｒ（Ａ）＾Ｔの各行の部分木のうち、エリミネーションツリーの各ノードを含む部分木の数を算出する。

算出部１００２は、例えば、スパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの非零要素のパターンを特定して、転置行列Ｒ（Ａ）＾Ｔの各行の部分木を抽出する。そして、算出部１００２は、エリミネーションツリーのノードごとに、当該ノードを含む部分木がいくつあるかを計数し、当該ノードを含む部分木の数を算出する。

ここで、算出部１００２は、転置行列Ｒ（Ａ）＾Ｔの非零要素のパターンを特定することができれば、上三角行列Ｒ（Ａ）や転置行列Ｒ（Ａ）＾Ｔを生成しなくてもよい。部分木の数を算出する詳細は、後述する実施例１の第６の工程において説明する。

算出部１００２は、生成したエリミネーションツリーの各ノードを含む部分木の数に基づいて、スパース行列ＡのＬＵ分解結果を格納するメモリ領域量を算出する。メモリ領域量とは、ＬＵ分解した結果を格納する領域の大きさである。算出部１００２は、例えば、生成したエリミネーションツリーの各ノードを含む、下三角行列Ｃ（Ａ）から得られた部分木の数を、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）の各列の非零要素の数とする。

また、算出部１００２は、生成したエリミネーションツリーの各ノードを含む、転置行列Ｒ（Ａ）＾Ｔから得られた部分木の数を、スパース行列ＡをＬＵ分解した場合の上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数とする。転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数は、上三角行列Ｕ（Ａ）の各行の非零要素の数に対応する。そして、算出部１００２は、下三角行列Ｌ（Ａ）の各列の非零要素の数と上三角行列Ｕ（Ａ）の各行の非零要素の数とに基づいて、ＬＵ分解した結果を格納する領域の大きさを算出する。メモリ領域量を算出する詳細は、後述する実施例１の第７の工程において説明する。

これにより、算出部１００２は、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを低減することができる。算出結果は、例えば、ＲＡＭ９０３、ディスク９０５などの記憶領域に記憶される。算出部１００２は、例えば、図９に示したＲＯＭ９０２、ＲＡＭ９０３、ディスク９０５などの記憶装置に記憶されたプログラムをＣＰＵ９０１に実行させることにより、その機能を実現する。

分解部１００３は、算出部１００２が算出した領域の大きさを確保して、スパース行列ＡをＬＵ分解する。分解部１００３は、例えば、ＲＡＭ９０３、ディスク９０５などの記憶領域に、算出部１００２が算出した大きさ分の領域を確保して、スパース行列ＡをＬＵ分解し、スパース行列Ａを分解した結果を確保した領域に格納する。分解部１００３は、例えば、図９に示したＲＯＭ９０２、ＲＡＭ９０３、ディスク９０５などの記憶装置に記憶されたプログラムをＣＰＵ９０１に実行させることにより、その機能を実現する。

（実施例１）
次に、図１１〜１７を用いて、実施例１について説明する。実施例１において、情報処理装置１００は、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを、スパース行列ＡをＬＵ分解するのに先立って近似的に算出し、ＬＵ分解した結果を格納する領域を確保してからＬＵ分解を行う。

ここでは、図５に示した非零要素のパターンを有するスパース行列Ａを例に挙げて、情報処理装置１００が、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）とを格納する領域の大きさを算出する各種工程について説明する。

＜第１の工程＞
まず、第１の工程について説明する。第１の工程は、情報処理装置１００が、スパース行列Ａから、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）とを生成し、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）との非零要素のパターンを特定する工程である。ここで、情報処理装置１００は、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）とのそれぞれを、対角要素が非零要素である行列として生成する。

情報処理装置１００は、例えば、スパース行列Ａと同じ大きさであって、スパース行列Ａの対角要素よりも右側および上側にある要素ａ［ｉ，ｊ］（ｉ＜ｊ）を零要素に変更した行列を、下三角行列Ｃ（Ａ）として生成する。以下の説明では、下三角行列Ｃ（Ａ）のｉ行ｊ列にある要素を「要素ｃ［ｉ，ｊ］」と表記する場合がある。そして、情報処理装置１００は、生成した下三角行列Ｃ（Ａ）を、圧縮列格納法を用いて格納しておく。圧縮列格納法については、図１７を用いて後述する。

また、情報処理装置１００は、スパース行列Ａと同じ大きさであって、スパース行列Ａの対角要素よりも左側および下側にある要素ａ［ｉ，ｊ］（ｉ＞ｊ）を零要素にした行列を、上三角行列Ｒ（Ａ）として生成する。以下の説明では、上三角行列Ｒ（Ａ）のｉ行ｊ列にある要素を「要素ｒ［ｉ，ｊ］」と表記する場合がある。そして、情報処理装置１００は、生成した上三角行列Ｒ（Ａ）を、圧縮行格納法を用いて格納しておく。圧縮行格納法については、図１７を用いて後述する。

換言すれば、情報処理装置１００は、実質的に、生成した上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔを、圧縮列格納法で格納しておくことになる。ここで、図１１を用いて、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）との非零要素のパターンについて説明する。

図１１は、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）との非零要素のパターンを示す説明図である。図１１の方眼１１０１は、下三角行列Ｃ（Ａ）の非零要素のパターンを示す。図１１の方眼１１０１のｉ行ｊ列目の升目は、下三角行列Ｃ（Ａ）のｉ行ｊ列目の要素ｃ［ｉ，ｊ］に対応し、下三角行列Ｃ（Ａ）のｉ行ｊ列目の要素ｃ［ｉ，ｊ］が対角要素、零要素、および非零要素のいずれであるかを示す。

図１１の例では、図１１の方眼１１０１の１行１列目の升目は、対応する下三角行列Ｃ（Ａ）の１行１列目の要素ｃ［１，１］が対角要素であるため、行番号「１」で示される。また、例えば、図１１の方眼１１０１の２行３列目の升目は、対応する下三角行列Ｃ（Ａ）の２行３列目の要素ｃ［２，３］が、スパース行列Ａの２行３列目の要素ａ［２，３］に関わらず零要素に変更されたため、「空白」で示される。また、例えば、図１１の方眼１１０１の３行２列目の升目は、対応する下三角行列Ｃ（Ａ）の３行２列目の要素ｃ［３，２］が、スパース行列Ａの３行２列目の要素ａ［３，２］そのものであり、非零要素であるため、「●」で示される。

一方で、図１１の方眼１１０２は、上三角行列Ｒ（Ａ）の非零要素のパターンを示す。図１１の方眼１１０２のｉ行ｊ列目の升目は、上三角行列Ｒ（Ａ）のｉ行ｊ列目の要素ｒ［ｉ，ｊ］に対応し、上三角行列Ｒ（Ａ）のｉ行ｊ列目の要素ｒ［ｉ，ｊ］が対角要素、零要素、および非零要素のいずれであるかを示す。

図１１の例では、図１１の方眼１１０２の１行１列目の升目は、対応する上三角行列Ｒ（Ａ）の１行１列目の要素ｒ［１，１］が対角要素であるため、行番号「１」で示される。また、例えば、図１１の方眼１１０２の２行３列目の升目は、対応する上三角行列Ｒ（Ａ）の２行３列目の要素ｒ［２，３］が、スパース行列Ａの２行３列目の要素ａ［２，３］そのものであり、非零要素であるため、「●」で示される。また、例えば、図１１の方眼１１０２の３行２列目の升目は、対応する上三角行列Ｒ（Ａ）の３行２列目の要素ｒ［３，２］が、スパース行列Ａの３行２列目の要素ａ［３，２］に関わらず零要素に変更されたため、「空白」で示される。

ここでは、情報処理装置１００が、下三角行列Ｃ（Ａ）と、上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔとを生成する場合について説明したが、これに限らない。例えば、情報処理装置１００は、下三角行列Ｃ（Ａ）と転置行列Ｒ（Ａ）＾Ｔとの非零要素のパターンを特定することができれば、下三角行列Ｃ（Ａ）と転置行列Ｒ（Ａ）＾Ｔを生成しなくてもよい。次に、情報処理装置１００は、第２の工程に移行する。

＜第２の工程＞
次に、第２の工程について説明する。第２の工程は、情報処理装置１００が、下三角行列Ｃ（Ａ）と、上三角行列Ｒ（Ａ）とから、スパース行列Ａを対称化した対称行列Ｐ＝Ａ＋Ａ＾Ｔの下三角行列Ｃ（Ｐ）を生成する工程である。

情報処理装置１００は、例えば、下三角行列Ｃ（Ａ）の非零要素のパターンと、上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの非零要素のパターンとを結合して、スパース行列Ａを対称化した対称行列Ｐの下三角行列Ｃ（Ｐ）を生成する。以下の説明では、対称行列Ｐのｉ行ｊ列にある要素を「要素ｐ［ｉ，ｊ］」と表記する場合がある。また、以下の説明では、下三角行列Ｃ（Ｐ）のｉ行ｊ列にある要素を「要素ｃｐ［ｉ，ｊ］」と表記する場合がある。ｉ＞ｊであれば、要素ｐ［ｉ，ｊ］＝要素ｃｐ［ｉ，ｊ］である。そして、情報処理装置１００は、生成した下三角行列Ｃ（Ｐ）の非零要素のパターンを特定する。ここで、図１２を用いて、下三角行列Ｃ（Ｐ）の非零要素のパターンについて説明する。

図１２は、下三角行列Ｃ（Ｐ）の非零要素のパターンを示す説明図である。図１２の方眼１２０１のｉ行ｊ列目の升目は、下三角行列Ｃ（Ｐ）のｉ行ｊ列目の要素ｃｐ［ｉ，ｊ］に対応し、下三角行列Ｃ（Ｐ）のｉ行ｊ列目の要素ｃｐ［ｉ，ｊ］が対角要素、零要素、および非零要素のいずれであるかを示す。

図１２の例では、図１２の方眼１２０１の２行３列目の升目は、対応する下三角行列Ｃ（Ｐ）の２行３列目の要素ｃｐ［２，３］が、零要素であるため、「空白」で示される。また、例えば、図１２の方眼１２０１の３行２列目の升目は、対応する下三角行列Ｃ（Ｐ）の３行２列目の要素ｃｐ［３，２］が、非零要素であるため、「●」で示される。下三角行列Ｃ（Ｐ）の３行２列目の要素ｃｐ［３，２］は、スパース行列Ａの下三角行列Ｃ（Ａ）の３行２列目の要素ｃ［３，２］と上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの３行２列目の要素ｒ［３，２］との和である。

ここでは、情報処理装置１００が、下三角行列Ｃ（Ｐ）を生成する場合について説明したが、これに限らない。例えば、情報処理装置１００は、下三角行列Ｃ（Ｐ）の非零要素のパターンを特定することができれば、下三角行列Ｃ（Ｐ）を生成しなくてもよい。次に、情報処理装置１００は、第３の工程に移行する。

＜第３の工程＞
次に、第３の工程について説明する。第３の工程は、情報処理装置１００が、対称行列ＰをＬＬ＾Ｔ分解した場合の下三角行列Ｌ（Ｐ）に基づいて、対称行列Ｐのエリミネーションツリーを生成する工程である。まず、情報処理装置１００は、対称行列Ｐの下三角行列Ｃ（Ｐ）の非零要素のパターンに基づいて、対称行列ＰをＬＬ＾Ｔ分解した場合の下三角行列Ｌ（Ｐ）を生成する。以下の説明では、下三角行列Ｌ（Ｐ）のｉ行ｊ列にある要素を「要素ｌ［ｉ，ｊ］」と表記する場合がある。

ここで、ＬＬ＾Ｔ分解の定義より、「対称行列Ｐ＝下三角行列Ｌ（Ｐ）・下三角行列Ｌ（Ｐ）の転置行列Ｌ（Ｐ）＾Ｔ」であるため、対称行列Ｐの各要素ｐ［ｉ，ｊ］について、下記式（１）が成立することになる。

さらに、上記式（１）を変形することにより、下三角行列Ｌ（Ｐ）の各要素ｌ［ｉ，ｊ］（ｉ＞ｊ）について、下記式（２）および下記式（３）が成立することになる。

上記式（２）および上記式（３）が成立するため、情報処理装置１００は、下三角行列Ｌ（Ｐ）を生成する際には、下三角行列Ｌ（Ｐ）の対角要素ｌ［ｊ，ｊ］をｊ＝１から順番に決定することになる。さらに、情報処理装置１００は、下三角行列Ｌ（Ｐ）の対角要素ｌ［ｊ，ｊ］を決定すると、下三角行列Ｌ（Ｐ）のｊ列目の各要素ｌ［１，ｊ］〜ｌ［ｊ−１，ｊ］をｊ＝１から順番に決定することになる。ここで、図１３を用いて、下三角行列Ｌ（Ｐ）の非零要素のパターンについて説明する。

図１３は、下三角行列Ｌ（Ｐ）の非零要素のパターンを示す説明図である。図１３の方眼１３０１のｉ行ｊ列目の升目は、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素に対応し、下三角行列Ｌ（Ｐ）のｉ行ｊ列目の要素が対角要素、零要素、および非零要素のいずれであるかを示す。図１３の例では、図１３の方眼１３０１の２行３列目の升目は、対応する下三角行列Ｌ（Ｐ）の２行３列目の要素ｌ［２，３］が、零要素であるため、「空白」で示される。また、例えば、図１３の方眼１３０１の３行２列目の升目は、対応する下三角行列Ｌ（Ｐ）の３行２列目の要素ｌ［３，２］が、非零要素であるため、「●」で示される。また、例えば、図１３の方眼１３０１の７行６列目の升目は、対応する下三角行列Ｌ（Ｐ）の７行６列目の要素ｌ［７，６］が、フィルインであるため、「○」で示される。

次に、情報処理装置１００は、下三角行列Ｌ（Ｐ）の非零要素のパターンに基づいて、対称行列ＰをＬＬ＾Ｔ分解してＬ（Ｐ）・Ｌ（Ｐ）＾Ｔで表現する場合の対称行列Ｐのエリミネーションツリーを生成する。

ここで、上記式（２）および上記式（３）によれば、下三角行列Ｌ（Ｐ）のｉ行ｋ列目が非零要素であれば、下三角行列Ｌ（Ｐ）のｉ列目の要素を算出する際に、下三角行列Ｌ（Ｐ）のｋ列目の要素が影響することになる。一方で、下三角行列Ｌ（Ｐ）のｉ行ｋ列目が零要素であれば、下三角行列Ｌ（Ｐ）のｉ列目の要素を算出する際に、下三角行列Ｌ（Ｐ）のｋ列目の要素が影響することはない。

エリミネーションツリーは、下三角行列Ｌ（Ｐ）の各列に関するノードを含み、ｉ列目の要素についてｊ列目の要素が影響する場合に、ｉ列目に関するノードとｊ列目に関するノードとを連結して親子関係とする。エリミネーションツリーは、例えば、ｉ列目の要素についてｊ列目の要素が影響すれば、ｉ列目に関するノードを親ノードとし、ｊ列目に関するノードを子ノードとする。

このため、ｊ列目に関するノード［ｊ］の親ノードは、ｍｉｎ｛ｉ｜ｉ＞ｊかつＬ［ｉ，ｊ］≠０｝を満たすｉ列目に関するノード［ｉ］になる。換言すれば、ｊ列目に関するノード［ｊ］の親ノードは、下三角行列Ｌ（Ｐ）のｊ列目において、対角要素以外の非零要素であって、最も対角要素の近くにある要素がある行の行番号と一致する列番号の列に関するノード［ｉ］になる。

生成したエリミネーションツリーは、ＬＬ＾Ｔ分解によってフィルインが発生する箇所を表現するツリーになる。また、ノード［ｊ］とノード［ｉ］との親子関係を、例えば、配列ｎｐａｒｅｎｔ［ｊ］を用いて、ｎｐａｒｅｎｔ［ｊ］＝ｉとして表現する場合がある。ここで、図１４を用いて、エリミネーションツリーについて説明する。

図１４は、エリミネーションツリー１４０１を示す説明図である。エリミネーションツリー１４０１は、１列目〜１１列目に関するノード［１］〜ノード［１１］を含む。図１４の例では、エリミネーションツリー１４０１において、ノード［１］の親ノードは、ノード［６］である。ノード［１］とノード［６］との親子関係は、ＬＬ＾Ｔ分解の過程で、下三角行列Ｌ（Ｐ）の１列目の要素が下三角行列Ｌ（Ｐ）の６列目の要素に影響し、下三角行列Ｌ（Ｐ）の６列目にフィルインが発生する可能性があることを示している。

また、エリミネーションツリー１４０１において、ノード［６］の親ノードは、ノード［７］である。ノード［６］とノード［７］との親子関係は、ＬＬ＾Ｔ分解の過程で、下三角行列Ｌ（Ｐ）の６列目の要素が下三角行列Ｌ（Ｐ）の７列目の要素に影響し、下三角行列Ｌ（Ｐ）の７列目にフィルインが発生する可能性があることを示している。このとき、下三角行列Ｌ（Ｐ）の１列目の要素が間接的に下三角行列Ｌ（Ｐ）の７列目の要素に影響し、下三角行列Ｌ（Ｐ）の７列目にフィルインが発生する可能性がある。次に、情報処理装置１００は、第４の工程に移行する。

＜第４の工程＞
次に、第４の工程について説明する。第４の工程は、情報処理装置１００が、エリミネーションツリー１４０１に関するパラメータを設定する工程である。情報処理装置１００は、例えば、親ノードを辿り、親ノードが辿れなくなったノードをエリミネーションツリー１４０１の根ノードとして設定する。また、情報処理装置１００は、あるノード［ｑ］の親ノードをノード［ｐ］としたとき、ノード［ｑ］をノード［ｐ］の子ノードとして設定する。

また、情報処理装置１００は、エリミネーションツリー１４０１の根ノードから深さ優先探索した場合に各ノードが探索された順番を、各ノードのポストオーダーとして付与する。図１４の例では、情報処理装置１００は、エリミネーションツリー１４０１のノード［２，３，１，６，７，４，５，８，９，１０，１１］の順に、ポストオーダー「１〜１１」のそれぞれを割り振る。また、情報処理装置１００は、エリミネーションツリー１４０１の根ノードから深さ優先探索した場合に、あるノードよりも深く探索することができなければ、当該ノードを葉ノードとして設定する。

また、情報処理装置１００は、各葉ノードから親ノードを辿り、エリミネーションツリー１４０１の根ノードまで遡った経路にある各ノードに、当該葉ノードを対応付けておく。また、情報処理装置１００は、各ノードに対応付ける葉ノードが複数ある場合には、複数の葉ノードの中で付与されたポストオーダーの最も小さいものを対応付ける。また、情報処理装置１００は、あるノードに対応付けた葉ノードを、あるノードの第１子孫（ｆｉｒｓｔ ｄｅｓｃｅｎｄａｎｔ）として設定する。次に、情報処理装置１００は、第５の工程に移行する。

＜第５の工程＞
次に、第５の工程について説明する。第５の工程は、情報処理装置１００が、スパース行列ＡをＬＵ分解したときの下三角行列Ｌ（Ａ）の各列の非零要素の数を算出する工程である。

情報処理装置１００は、例えば、エリミネーションツリー１４０１と、スパース行列Ａの下三角行列Ｃ（Ａ）の非零要素のパターンとに基づいて、スパース行列ＡをＬＵ分解したときの下三角行列Ｌ（Ａ）の各列の非零要素の数を算出する。

ここで、対称行列Ｐの下三角行列Ｃ（Ｐ）のｊ列目のベクトルをｂｊとしたとき、ＬＵ分解した下三角行列Ｌ（Ａ）の非零要素のパターンは、ｂｊと、ノード［ｊ］の子ノード［ｋ］に関する下三角行列Ｌ（Ａ）のｋ列目のベクトルｂｋとの和集合になる。このため、下三角行列Ｌ（Ａ）のｉ行目の非零要素は、エリミネーションツリー１４０１の部分木として表現することができる。例えば、下三角行列Ｌ（Ａ）のｉ行目の非零要素は、ノード［ｉ］を根ノードとする部分木として表現することができる。ここで、図１５および図１６を用いて、部分木の一例について説明する。

図１５は、７行目の非零要素を表現する部分木１５０１の一例を示す説明図である。ここで、情報処理装置１００は、対角要素以外の非零要素がある列番号を特定して、部分木１５０１を抽出する。図１５の例では、情報処理装置１００は、スパース行列Ａの７行目の非零要素が、対角要素を除いて、１，３列目にあると特定する。

７行目の非零要素を表現する部分木１５０１であるため、７列目に関するノード［７］が根ノードになる。また、非零要素がある１列目に関するノード［１］が、葉ノードになり、ノード［６］を経由して、根ノードになるノード［７］まで連結されている。また、非零要素がある３列目に関するノード［３］が、葉ノードになり、ノード［７］まで連結されている。

このため、ノード［７］を根ノードとする部分木１５０１は、エリミネーションツリー１４０１のうちのノード［１，６，３，７］を含むツリーになる。ノード［１，３］は、葉ノードである。ここで、図１５の部分木１５０１は、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）の７行目において、１，３，６列目にフィルインが発生する可能性があることを表現している。

図１６は、１１行目の非零要素を表現する部分木１６０１の一例を示す説明図である。ここで、情報処理装置１００は、対角要素以外の非零要素がある列番号を特定して、部分木１６０１を抽出する。図１６の例では、情報処理装置１００は、スパース行列Ａの１１行目の非零要素が、対角要素を除いて、１，５，８列目にあると特定する。

１１行目の非零要素を表現する部分木１６０１であるため、１１列目に関するノード［１１］が根ノードになる。また、非零要素がある１列目に関するノード［１］が、葉ノードになり、ノード［６，７，８，９，１０］を経由して、根ノードになるノード［１１］まで連結されている。また、非零要素がある５列目に関するノード［５］が、葉ノードになり、ノード［８，９，１０］を経由して、ノード［１１］まで連結されている。また、非零要素がある８列目に関するノード［８］が、ノード［９，１０］を経由して、ノード［１１］まで連結されている。

このため、ノード［１１］を根ノードとする部分木１６０１は、エリミネーションツリー１４０１のうちのノード［１，５〜１１］を含むツリーになる。ノード［１，５］は、葉ノードである。ここで、図１６の部分木１６０１は、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）の１１行目において、１，５〜１１列目にフィルインが発生する可能性があることを表現している。

情報処理装置１００は、具体的には、下三角行列Ｃ（Ａ）の１１行目にある非零要素がある列に関するノードを、付与されたポストオーダーの順に取り出す。ここで、１，５，８列目の要素に関するノード［１，５，８］のそれぞれの第１子孫はノード［１，４，２］が設定されている。情報処理装置１００は、ノード［１］を最初に取り出したため、ノード［１］を葉ノードにする。

また、情報処理装置１００は、ノード［５］を取り出すと、ノード［５］の第１子孫となるノード［４］に付与されたポストオーダーが、ノード［１］に付与されたポストオーダーより大きいことを検出する。このため、情報処理装置１００は、分枝ノードで枝分かれしているとして、ノード［５］も葉ノードにする。また、情報処理装置１００は、ノード［８］の第１子孫となるノード［２］に付与されたポストオーダーが、ノード［５］に付与されたポストオーダーより小さいため、ノード［５，８］の間で枝分かれはないとして、ノード［８］を葉ノードにしない。そして、情報処理装置１００は、エリミネーションツリー１４０１から、各葉ノードから根ノードまでを含む部分木１６０１を抽出する。

このように、情報処理装置１００は、スパース行列Ａの行の非零要素がある列に関するノードを、付与されたポストオーダーの順に取り出す。次に、情報処理装置１００は、一つ前に取り出したノードに付与されたポストオーダーと、現在取り出したノードの第１子孫に付与されたポストオーダーを比較する。ここで、一つ前に取り出したノードと、現在取り出したノードとの２つのノードは、深さ優先探索でポストオーダーを付与したため、現在取り出したノードの第１子孫に付与されたポストオーダーの方が大きければ、共通の祖先で分枝していることになる。そして、情報処理装置１００は、比較した結果、現在取り出したノードの第１子孫のほうが大きければ、現在取り出したノードを葉ノードにして、エリミネーションツリー１４０１から部分木１６０１を抽出する。

また、情報処理装置１００は、行の非零要素がある列に関するノードを、付与されたポストオーダーの順に取り出して、一つ前のノードを記憶しておく代わりに、一つ前の葉ノードを記憶しておき新たな葉ノードが見つかったときに更新するようにしてもよい。

ここで、部分木１５０１，１６０１などは、各行の非零要素を表現している。このため、下三角行列Ｌ（Ａ）のｊ番目の列の非零要素の数は、部分木１５０１，１６０１などといったエリミネーションツリー１４０１の部分木のうちの、ノード［ｊ］を含む部分木の数になる。これにより、情報処理装置１００は、Ｏ（｜Ｌ（Ａ）｜）の演算量で、非零要素の数を算出することができる。ここで、｜Ｌ（Ａ）｜は行列Ｌ（Ａ）の非零要素の数を表す。

ここで、下三角行列Ｃ（Ａ）および上三角行列Ｒ（Ａ）の非零要素のパターンは、対称行列Ｐの非零要素のパターンの部分集合になる。このため、下三角行列Ｃ（Ａ）および上三角行列Ｒ（Ａ）の非零要素のパターンと、対称行列Ｐの非零要素のパターンとには包含関係があり、下三角行列Ｃ（Ａ）⊆対称行列Ｐ、かつ、上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔ⊆対称行列Ｐが成立する。

また、下三角行列Ｃ（Ａ）または上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの非零要素を調べて抽出した部分木は、対称行列Ｐから抽出した部分木よりもノードの数が少なくなる。また、対称行列Ｐのエリミネーションツリー１４０１を使っているので、上述した第５の工程において算出した非零要素の数は、実際の非零要素の数より多くなる可能性がある。すなわち、下三角行列Ｃ（Ａ）および上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの各列の非零要素の数（ｃｏｌｕｍｎ ｃｏｕｎｔ）は、対応する対称行列Ｐの各列の非零要素の数以下となる。

＜第５の工程の他の例＞
次に、第５の工程の他の例について説明する。第５の工程の他の例は、情報処理装置１００が、特性関数を用いて、スパース行列ＡをＬＵ分解したときの下三角行列Ｌ（Ａ）の各列の非零要素の数を算出する一例である。

ここで、下記式（４）および下記式（５）に示す特性関数を用意する。ｒｏｗ ｓｕｂｔｒｅｅ ｉは、ｉ行目のロウサブツリーである。

特性関数は、部分木１６０１の葉ノードに１を設定しておき、葉ノード以外であれば０を設定しておく。そして、特性関数は、情報処理装置１００によって、各ノードに付与されたポストオーダーの順にエリミネーションツリー１４０１が辿られ、子ノードの特性関数の値が伝播され加算されることにより、更新される。

例えば、特性関数は、部分木１６０１の葉ノードのときは、１を設定しておく。また、特性関数は、部分木１６０１の構成ノードであって、部分木１６０１の葉ノードに行き当たる子ノードがｄ個あるとき、１−ｄを加算する。また、特性関数は、部分木１６０１の根ノードの親ノードなら−１を加算する。

情報処理装置１００は、実際には、部分木１６０１の根ノードに対応する行をポストオーダーの順にスキャンして、一つ前の葉ノードとペアにして、ペアのｃｏｍｍｏｎ ａｎｃｅｓｔｏｒのノードに−１を加えることで算出することができる。ａｎｃｅｓｔｏｒは、祖先のノードである。ｃｏｍｍｏｎ ａｎｃｅｓｔｏｒは、ペアのノードに共通する祖先のノードであって、ペアのノードから最も近いノードである。

情報処理装置１００は、各ノードに対して変数ｃｏｌｕｍｎ ｃｏｕｎｔを用意し、ポストオーダーの順にノードを辿り、親ノードのｃｏｌｕｍｎ ｃｏｕｎｔに子ノードのｃｏｌｕｍｎ ｃｏｕｎｔを加算していく。これにより、各部分木１６０１に含まれれば、葉ノードに設定した「１」がエリミネーションツリー１４０１の枝を伝播していき、特性関数を実現することができる。子ノードが多いノードでは値が調整され、部分木１６０１の根ノードの親ノードで値の伝播がキャンセルされる。

＜第６の工程＞
次に、第６の工程について説明する。第６の工程は、情報処理装置１００が、スパース行列ＡをＬＵ分解したときの上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数を算出する工程である。情報処理装置１００は、例えば、エリミネーションツリー１４０１と、スパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの非零要素のパターンとから、転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数を算出する。

ここで、情報処理装置１００は、第５の工程におけるスパース行列Ａの下三角行列Ｃ（Ａ）を、スパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔに置き換えれば、同様に、転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数を算出することができる。このため、転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数を算出する説明については省略する。

＜第７の工程＞
次に、第７の工程について説明する。第７の工程は、情報処理装置１００が、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを算出する工程である。

情報処理装置１００は、例えば、スパース行列ＡをＬＵ分解したときの下三角行列Ｌ（Ａ）の非零要素の総数に基づいて、スパース行列ＡをＬＵ分解したときの下三角行列Ｌ（Ａ）を格納する領域の大きさを算出する。また、情報処理装置１００は、スパース行列ＡをＬＵ分解したときの上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの非零要素の総数に基づいて、スパース行列ＡをＬＵ分解したときの上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔを格納する領域の大きさを算出する。

情報処理装置１００は、具体的には、下三角行列Ｌ（Ａ）の各列のｃｏｌｕｍｎ ｃｏｕｎｔをすべて加算することにより下三角行列Ｌ（Ａ）の列を圧縮列格納法で格納するときに必要な領域の大きさを算出する。上三角行列Ｕ（Ａ）は、単位上三角行列であるため、上三角行列Ｕ（Ａ）の対角要素は１である。すなわち、上三角行列Ｕ（Ａ）の対角要素は、格納しなくてもよい。ここで、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの各列のｃｏｌｕｍｎ ｃｏｕｎｔは、上三角行列Ｕ（Ａ）の各行の非零要素の数になる。このため、情報処理装置１００は、上三角行列Ｕ（Ａ）の各行の非零要素の数から１を減算した値をすべて加算することにより、対角要素を除いた非零要素を圧縮行格納法で格納するときに必要な領域の大きさを算出する。

・圧縮列格納法の一例
ここで、図１７を用いて、圧縮列格納法の一例について説明する。図１７は、圧縮列格納法の一例を示す説明図である。

情報処理装置１００は、図１７の行列ｍａｔの各列の非零要素を圧縮して、配列ａに順次格納する。次に、情報処理装置１００は、配列ａに格納された要素が、何行目に位置する要素であるかを示す情報を、配列ｎｒｏｗに同じ順序で格納する。そして、情報処理装置１００は、各列の最初の非零要素が、何番目の配列ａに格納されるかを示す情報を、配列ｎｆｃｎｚに格納する。

ここで、行列ｍａｔの次数をｎ、非零要素の総数をｎｚとしたとき、１次元配列ｎｆｃｎｚの大きさはｎ＋１になり、配列ｎｆｃｎｚの６つ目の要素にはｎｚ＋１となる仮想位置が格納される。配列ｎｆｃｎｚは、配列ｎｒｏｗに対しても、配列ａと同様に位置を示す。配列ａおよび配列ｎｒｏｗは大きさｎｚの１次元配列である。

配列ａは、例えば、倍精度複素数型である。配列ｎｆｃｎｚや配列ｎｒｏｗは、例えば、整数型である。ここで、圧縮行格納法は、格納する行列を転置して圧縮列格納法で格納する場合と同様であるため、説明を省略する。

（算出処理手順の一例）
次に、図１８を用いて、実施例１にかかる算出処理手順の一例について説明する。

図１８は、実施例１にかかる算出処理手順の一例を示すフローチャートである。図１８において、まず、情報処理装置１００は、スパース行列Ａから、スパース行列Ａの下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）とを生成する（ステップＳ１８０１）。

次に、情報処理装置１００は、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔとの非零要素のパターンをマージして、スパース行列Ａの対称行列Ｐの下三角行列Ｃ（Ｐ）の非零要素のパターンを生成する（ステップＳ１８０２）。そして、情報処理装置１００は、対称行列Ｐの下三角行列Ｃ（Ｐ）の非零要素のパターンから、対称行列Ｐのエリミネーションツリー１４０１を生成する（ステップＳ１８０３）。

次に、情報処理装置１００は、エリミネーションツリー１４０１に関するパラメータを特定する（ステップＳ１８０４）。そして、情報処理装置１００は、エリミネーションツリー１４０１と、スパース行列Ａの下三角行列Ｃ（Ａ）の非零要素のパターンとから、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）の各列の非零要素の数を近似する（ステップＳ１８０５）。

次に、情報処理装置１００は、エリミネーションツリー１４０１と、スパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの非零要素のパターンとから、スパース行列ＡをＬＵ分解した場合の上三角行列Ｕ（Ａ）の各行の非零要素の数を近似する（ステップＳ１８０６）。そして、情報処理装置１００は、下三角行列Ｌ（Ａ）の各列の非零要素の数の総和を算出し、下三角行列Ｌ（Ａ）の非零要素を格納する領域の大きさを算出する（ステップＳ１８０７）。

次に、情報処理装置１００は、上三角行列Ｕ（Ａ）の各行の非零要素の数の総和を算出し、上三角行列Ｕ（Ａ）の非零要素を格納する領域の大きさを算出する（ステップＳ１８０８）。そして、情報処理装置１００は、算出処理を終了する。これにより、情報処理装置１００は、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを算出することができる。

（実施例２）
次に、実施例２について説明する。実施例２は、情報処理装置１００が、スーパーノードを用いて下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）とを格納する領域の大きさを算出する一例である。実施例２において、情報処理装置１００は、実施例１と同様に、第１の工程〜第６の工程によって、スパース行列ＡをＬＵ分解したときの下三角行列Ｌ（Ａ）の各列の非零要素の数と、上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔの各列の非零要素の数とを算出する。

実施例２では、情報処理装置１００は、対称行列Ｐの非零要素のパターンから、対称行列Ｐの各行のｃｏｌｕｍｎ ｃｏｕｎｔを算出して、対称行列ＰをＬＬ＾Ｔ分解するときの複数のノードを纏めたスーパーノードを特定する。スーパーノードは、エリミネーションツリー１４０１の連続する複数のノードを纏めたものである。スーパーノードは、インデックスが大きい方のノードに対応する列にある非零要素のパターンが、インデックスが小さい方のノードに対応する列にある非零要素のパターンと一致する場合に、複数のノードを纏めたものである。

また、スーパーノードは、インデックスが大きい方のノードに対応する列にある非零要素のパターンが、インデックスが小さい方のノードに対応する列にある非零要素のパターンと類似する場合に、複数のノードを纏めたものであってもよい。スーパーノードとして纏められた複数のノードに対応する複数の列をパネル（ｐａｎｅｌ）とする。

ここで、スーパーノードとして纏められる複数のノードが満たす条件は、複数のノードのうちの親ノードが、スーパーノードとして纏められる複数のノードに対応する複数の列を纏めたパネルの右終端の列になることである。そして、スーパーノードとして纏められる複数のノードが満たす条件は、複数のノードのうちの他ノードが当該親ノードの子孫になることである。さらに、スーパーノードとして纏められる複数のノードが満たす条件は、当該親ノードと当該子孫との間にある他ノードが含まれることである。子ノードをマージしてスーパーノードを生成するとき、親ノードがパネルの右終端の列に対応するようにすれば、対角部分以外での非零要素の数は、親ノードの「ｃｏｌｕｍｎ ｃｏｕｎｔ−１」となる。

情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解するときのスーパーノードとして纏められた複数のノードを、スパース行列ＡをＬＵ分解するときの下三角行列Ｌ（Ａ）についてのスーパーノードとして纏める複数のノードとして採用する。情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解するときのスーパーノードとして纏められた複数のノードを、スパース行列ＡをＬＵ分解するときの上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔについてのスーパーノードとして纏める複数のノードとして採用する。これにより、情報処理装置１００は、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）の転置行列Ｕ（Ａ）＾Ｔとについて、スーパーノードを格納する領域の大きさを計算できる。

スーパーノードの大きさとしてスーパーノードに含まれるノードの数をｂとし、右端のノードをｅとすれば、スーパーノードに対応する下三角行列Ｌ（Ａ）の複数の列は、纏めて、（ｂ＋ｌｃｃ（ｅ）−１）×ｂの大きさのパネルに格納される。ここで、ｌｃｃ（ｅ）は、下三角行列Ｌ（Ａ）のノードｅのｃｏｌｕｍｎ ｃｏｕｎｔである。また、スーパーノードに対応する上三角行列Ｕ（Ａ）の複数の列は、対角要素は下三角行列Ｌ（Ａ）のｐａｎｅｌに格納されるため、残りの要素を（ｕｔｃｃ（ｅ）−１）×ｂのパネルに格納される。ここで、ｕｔｃｃ（ｅ）は、上三角行列Ｕ（Ａ）のノードｅのｃｏｌｕｍｎ ｃｏｕｎｔである。

情報処理装置１００は、エリミネーションツリー１４０１において連続するノードの集合であり、ノードに対応する列にある非零要素のパターンが一致するものを纏めたスーパーノードについて、非零要素がある列または行を圧縮した形式で格納することができる。このため、情報処理装置１００は、非零要素がある列または行のデータを格納し、非零要素がない列または行のデータを格納しなくてもよくなり、格納する領域の大きさを低減することができる。

（算出処理手順の一例）
次に、図１９を用いて、実施例２にかかる算出処理手順の一例について説明する。

図１９は、実施例２にかかる算出処理手順の一例を示すフローチャートである。図１９において、まず、情報処理装置１００は、スパース行列Ａから、スパース行列Ａの下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）とを生成する（ステップＳ１９０１）。

次に、情報処理装置１００は、下三角行列Ｃ（Ａ）と上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔとの非零要素のパターンをマージして、スパース行列Ａの対称行列Ｐの下三角行列Ｃ（Ｐ）の非零要素のパターンを生成する（ステップＳ１９０２）。そして、情報処理装置１００は、対称行列Ｐの下三角行列Ｃ（Ｐ）の非零要素のパターンから、対称行列Ｐのエリミネーションツリー１４０１を生成する（ステップＳ１９０３）。

次に、情報処理装置１００は、エリミネーションツリー１４０１に関するパラメータを特定する（ステップＳ１９０４）。そして、情報処理装置１００は、対称行列Ｐの下三角行列Ｃ（Ｐ）から、対称行列ＰをＬＬ＾Ｔ分解した場合の下三角行列Ｌ（Ｐ）の各列の非零要素の数を算出する（ステップＳ１９０５）。

次に、情報処理装置１００は、エリミネーションツリー１４０１と、スパース行列Ａの下三角行列Ｃ（Ａ）の非零要素のパターンとから、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）の各列の非零要素の数を近似する（ステップＳ１９０６）。そして、情報処理装置１００は、エリミネーションツリー１４０１と、スパース行列Ａの上三角行列Ｒ（Ａ）の転置行列Ｒ（Ａ）＾Ｔの非零要素のパターンとから、スパース行列ＡをＬＵ分解した場合の上三角行列Ｕ（Ａ）の各行の非零要素の数を近似する（ステップＳ１９０７）。ステップＳ１９０６とステップＳ１９０７との処理は、より具体的には、後述するｃｏｌｕｍｎ ｃｏｕｎｔを算出する処理を行うことにより実現される。

次に、情報処理装置１００は、対称行列ＰをＬＬ＾Ｔ分解した場合のスーパーノードを特定する（ステップＳ１９０８）。そして、情報処理装置１００は、下三角行列Ｌ（Ａ）と上三角行列Ｕ（Ａ）とのスーパーノードに対応する部分を格納するパネル（ｐａｎｅｌ）の大きさを算出する（ステップＳ１９０９）。

次に、情報処理装置１００は、下三角行列Ｌ（Ａ）のスーパーノードに対応するパネルの大きさを加えて、下三角行列Ｌ（Ａ）を格納する領域の大きさを算出する（ステップＳ１９１０）。そして、情報処理装置１００は、上三角行列Ｕ（Ａ）のスーパーノードに対応するパネルの大きさを加えて、上三角行列Ｕ（Ａ）を格納する領域の大きさを算出する（ステップＳ１９１１）。その後、情報処理装置１００は、算出処理を終了する。これにより、情報処理装置１００は、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを算出することができる。

（ｃｏｌｕｍｎ ｃｏｕｎｔを算出する詳細）
次に、情報処理装置１００が、ｃｏｌｕｍｎ ｃｏｕｎｔを算出する詳細について説明する。ｃｏｌｕｍｎ ｃｏｕｎｔを算出する処理は、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）の各列の非零要素の数と、上三角行列Ｕ（Ａ）の各行の非零要素の数とを近似する、ステップＳ１９０６とステップＳ１９０７との処理に対応する。

ここで、ノードの数をｎとする。ノード［ｊ］がノード［ｉ］の親ノードであることをｊ＝ｎｐａｒｅｎｔ（ｉ）とする。情報処理装置１００は、１次元配列ｎｒｏｗ（ｎｚ）を用意する。ｎｚは下三角行列の非零要素の総数である。１次元配列ｎｒｏｗ（ｎｚ）は、各列の非零要素の行番号が格納される配列である。また、情報処理装置１００は、１次元配列ｎｐａｒｅｎｔ（ｎ）を用意する。１次元配列ｎｐａｒｅｎｔ（ｎ）は、エリミネーションツリー１４０１を表現する配列である。また、情報処理装置１００は、１次元配列ｎｐｏｓｔｏ（ｎ）を用意する。１次元配列ｎｐｏｓｔｏ（ｎ）は、ポストオーダーが格納される配列である。

また、情報処理装置１００は、１次元配列ｎｐｏｓｔｏｉｎｖ（ｎ）を用意する。１次元配列ｎｐｏｓｔｏｉｎｖ（ｎ）は、ポストオーダー順にノードが格納される配列である。また、情報処理装置１００は、１次元配列ｎｆｉｒｓｔｄｅｓｃｅｎｄａｎｔ（ｎ）を用意する。１次元配列ｎｆｉｒｓｔｄｅｓｃｅｎｄａｎｔ（ｎ）は、各ノードのｆｉｒｓｔ ｄｅｓｃｅｎｄａｎｔが格納される配列である。

また、情報処理装置１００は、１次元配列ｎｐｒｅｖｐ（ｎ）を用意する。１次元配列ｎｐｒｅｖｐ（ｎ）は、ロウサブツリーの一つ前に検出された葉ノードが格納される配列である。１次元配列ｎｐｒｅｖｐ（ｎ）の初期値は０である。また、情報処理装置１００は、１次元配列ｎｒｏｗｓｕｂｆｌａｇ（ｎ）を用意する。１次元配列ｎｒｏｗｓｕｂｆｌａｇ（ｎ）は、ロウサブツリーが葉ノードを持つか否かを示す情報が格納される配列である。１次元配列ｎｒｏｗｓｕｂｆｌａｇ（ｎ）の初期値は０である。

また、情報処理装置１００は、１次元配列ｎｓｋｅｌｅｔｏｎｍａｔ（ｎ）を用意する。１次元配列ｎｓｋｅｌｅｔｏｎｍａｔ（ｉ）は、ノード［ｉ］がロウサブツリーの葉ノードであれば１が加算される。１次元配列ｎｓｋｅｌｅｔｏｎｍａｔ（ｉ）の初期値は０である。また、情報処理装置１００は、１次元配列ｎｄｅｌｔａ（ｎ）を用意する。１次元配列ｎｄｅｌｔａ（ｎ）は、葉ノードのペアにとってのｃｏｍｍｏｎ ａｎｃｅｓｔｏｒ ｉｃａに対応する要素が格納され、−１を加算される。１次元配列ｎｄｅｌｔａ（ｎ）の初期値は０である。

また、情報処理装置１００は、１次元配列ｎｃｃｏｕｎｔ（ｎ）を用意する。１次元配列ｎｃｃｏｕｎｔ（ｎ）は、各ノード対する列の非零要素数ｃｏｌｕｍｎ ｃｏｕｎｔが格納される配列である。また、情報処理装置１００は、１次元配列ｎａｎｃｅｓｔｏｒ（ｎ）を用意する。１次元配列ｎａｎｃｅｓｔｏｒ（ｎ）は、ポストオーダー順に処理した親ノードを格納する。１次元配列ｎａｎｃｅｓｔｏｒ（ｎ）の初期値はｉである。

＜ｃｏｌｕｍｎ ｃｏｕｎｔの計数処理手順＞
次に、図２０〜図２２を用いて、ｃｏｌｕｍｎ ｃｏｕｎｔの計数処理手順について説明する。以下の説明では、「ａ＝＝ｂ」はａとｂが一致することを示す。「ａ．ｎｅ．ｂ」はａとｂが一致しないことを示す。「ａ＝ｂ」はａにｂを代入することを示す。「ａ：ｂ」はａ〜ｂを示す。

図２０〜図２２は、ｃｏｌｕｍｎ ｃｏｕｎｔの計数処理手順の一例を示すフローチャートである。図２０において、情報処理装置１００は、配列を初期化する（ステップＳ２００１）。情報処理装置１００は、例えば、１次元配列ｎｐｒｅｖｐ（１：ｎ）＝０、１次元配列ｎｓｋｅｌｅｔｏｎｍａｔ（１：ｎ）＝０、１次元配列ｎｄｅｌｔａ（１：ｎ）＝０、１次元配列ｎｒｏｗｓｕｂｆｌａｇ（１：ｎ）＝０を設定する。また、情報処理装置１００は、例えば、１次元配列ｎａｎｃｅｓｔｏｒ（ｋ）＝ｋ、ｋ＝１，・・・，ｎ、ｉ＝１を設定する。

次に、情報処理装置１００は、ｎｏｄｅｐ＝ｎｐｏｓｔｏ（ｉ）として、ポストオーダーが「ｉ」であるノードのインデックスを取得する（ステップＳ２００２）。そして、情報処理装置１００は、ｊ＝ｎｆｃｎｚ（ｎｏｄｅｐ）とする（ステップＳ２００３）。

次に、情報処理装置１００は、ｎｏｄｅｕ＝ｎｒｏｗ（ｊ）とする（ステップＳ２００４）。そして、情報処理装置１００は、ｎｏｄｅｕ＞ｎｏｄｅｐであるか否かを判定する（ステップＳ２００５）。ここで、ｎｏｄｅｕ＞ｎｏｄｅｐではない場合（ステップＳ２００５：Ｎｏ）、情報処理装置１００は、ステップＳ２１０６の処理に移行する。

一方で、ｎｏｄｅｕ＞ｎｏｄｅｐである場合（ステップＳ２００５：Ｙｅｓ）、情報処理装置１００は、ｎｒｏｗｓｕｂｆｌａｇ（ｎｏｄｅｕ）＝＝０であるか否かを判定する（ステップＳ２００６）。ここで、ｎｒｏｗｓｕｂｆｌａｇ（ｎｏｄｅｕ）＝＝０ではない場合（ステップＳ２００６：Ｎｏ）、情報処理装置１００は、ステップＳ２００８の処理に移行する。

一方で、ｎｒｏｗｓｕｂｆｌａｇ（ｎｏｄｅｕ）＝＝０である場合（ステップＳ２００６：Ｙｅｓ）、情報処理装置１００は、ｎｒｏｗｓｕｂｆｌａｇ（ｎｏｄｅｕ）＝１とする（ステップＳ２００７）。次に、情報処理装置１００は、ｎｐｒｅｖｎｂｒｎｏｄｅｕ＝ｎｐｒｅｖｐ（ｎｏｄｅｕ）とする（ステップＳ２００８）。

そして、情報処理装置１００は、ｎｐｒｅｖｎｂｒｎｏｄｅｕ≠０であるか否かを判定する（ステップＳ２００９）。ここで、ｎｐｒｅｖｎｂｒｎｏｄｅｕ≠０ではない場合（ステップＳ２００９：Ｎｏ）、情報処理装置１００は、ステップＳ２１０１の処理に移行する。

一方で、ｎｐｒｅｖｎｂｒｎｏｄｅｕ≠０である場合（ステップＳ２００９：Ｙｅｓ）、情報処理装置１００は、ｎｐｒｅｖｎｂｒｎｏｄｅｕ＝ｎｐｏｓｔｏｉｎｖ（ｎｐｒｅｖｎｂｒｎｏｄｅｕ）として、ｎｐｒｅｖｎｂｒｎｏｄｅｕのポストオーダーを取得する（ステップＳ２０１０）。次に、情報処理装置１００は、図２１のステップＳ２１０１の処理に移行する。

図２１において、情報処理装置１００は、葉ノードであるか否かをチェックするために、ｎｐｏｓｔｏｉｎｖ（ｎｆｉｒｓｔｄｅｓｃｅｎｄａｎｔ（ｎｏｄｅｐ））＞ｎｐｒｅｖｎｂｒｎｏｄｅｕであるか否かを判定する（ステップＳ２１０１）。ここで、ｎｐｏｓｔｏｉｎｖ（ｎｆｉｒｓｔｄｅｓｃｅｎｄａｎｔ（ｎｏｄｅｐ））＞ｎｐｒｅｖｎｂｒｎｏｄｅｕではない場合（ステップＳ２１０１：Ｎｏ）、情報処理装置１００は、ステップＳ２１０６の処理に移行する。

一方で、ｎｐｏｓｔｏｉｎｖ（ｎｆｉｒｓｔｄｅｓｃｅｎｄａｎｔ（ｎｏｄｅｐ））＞ｎｐｒｅｖｎｂｒｎｏｄｅｕである場合（ステップＳ２１０１：Ｙｅｓ）、情報処理装置１００は、ステップＳ２１０２の処理に移行する。ステップＳ２１０２において、情報処理装置１００は、ｎｓｋｅｌｅｔｏｎｍａｔ（ｎｏｄｅｐ）＝ｎｓｋｅｌｅｔｏｎｍａｔ（ｎｏｄｅｐ）＋１とし、ｎｏｄｅｐｐ＝ｎｐｒｅｖｐ（ｎｏｄｅｕ）とする（ステップＳ２１０２）。

次に、情報処理装置１００は、ｎｏｄｅｐｐ≠０であるか否かを判定する（ステップＳ２１０３）。ここで、ｎｏｄｅｐｐ≠０ではない場合（ステップＳ２１０３：Ｎｏ）、情報処理装置１００は、ステップＳ２１０５の処理に移行する。

一方で、ｎｏｄｅｐｐ≠０である場合（ステップＳ２１０３：Ｙｅｓ）、情報処理装置１００は、ｃｏｍｍｏｎ ａｎｃｅｓｔｏｒを探索し、ｎｏｄｅｑ＝ｎｏｄｅｐｐとし、ｎｏｄｅｑ＝ｎａｎｃｅｓｔｏｒ（ｎｏｄｅｑ）とする。そして、情報処理装置１００は、条件（ｎｏｄｅｑ．ｎｅ．ｎａｎｃｅｓｔｏｒ（ｎｏｄｅｑ））を満たすまでｎｏｄｅｑ＝ｎａｎｃｅｓｔｏｒ（ｎｏｄｅｑ）を繰り返し、ｎｄｅｌｔａ（ｎｏｄｅｑ）＝ｎｄｅｌｔａ（ｎｏｄｅｑ）−１とする（ステップＳ２１０４）。

次に、情報処理装置１００は、ｎｐｒｅｖｐ（ｎｏｄｅｕ）＝ｎｏｄｅｐとする（ステップＳ２１０５）。そして、情報処理装置１００は、ｊ＝ｊ＋１とする（ステップＳ２１０６）。次に、情報処理装置１００は、ｊ＞ｎｆｃｎｚ（ｎｏｄｅ＋１）−１であるか否かを判定する（ステップＳ２１０７）。ここで、ｊ＞ｎｆｃｎｚ（ｎｏｄｅ＋１）−１ではない場合（ステップＳ２１０７：Ｎｏ）、情報処理装置１００は、ステップＳ２００４の処理に移行する。

一方で、ｊ＞ｎｆｃｎｚ（ｎｏｄｅ＋１）−１である場合（ステップＳ２１０７：Ｙｅｓ）、情報処理装置１００は、ｎｐａｒｅｎｔ（ｎｏｄｅｐ）≠０であるか否かを判定する（ステップＳ２１０８）。ここで、ｎｐａｒｅｎｔ（ｎｏｄｅｐ）≠０ではない場合（ステップＳ２１０８：Ｎｏ）、情報処理装置１００は、ステップＳ２１１０の処理に移行する。

一方で、ｎｐａｒｅｎｔ（ｎｏｄｅｐ）≠０である場合（ステップＳ２１０８：Ｙｅｓ）、情報処理装置１００は、ｎａｎｃｅｓｔｏｒ（ｎｏｄｅｐ）＝ｎｐａｒｅｎｔ（ｎｏｄｅｐ）とする（ステップＳ２１０９）。次に、情報処理装置１００は、ｉ＝ｉ＋１とする（ステップＳ２１１０）。そして、情報処理装置１００は、ｉ＞Ｎであるか否かを判定する（ステップＳ２１１１）。ここで、ｉ＞Ｎではない場合（ステップＳ２１１１：Ｎｏ）、情報処理装置１００は、ステップＳ２００２の処理に移行する。一方で、ｉ＞Ｎである場合（ステップＳ２１１１：Ｙｅｓ）、情報処理装置１００は、図２２のステップＳ２２０１の処理に移行する。

図２２において、情報処理装置１００は、ｎｃｃｏｕｎｔ（ｉ）＝ｎｓｋｅｌｅｔｏｎｍａｔ（ｉ）＋ｎｄｅｌｔａ（ｉ）とし、ｎｒｏｗｓｕｂｆｌａｇ（ｉ）＝＝０であればｎｃｃｏｕｎｔ（ｉ）＝ｎｃｃｏｕｎｔ（ｉ）＋１とする処理をｉ＝１〜ｎまで繰り返す（ステップＳ２２０１）。

次に、情報処理装置１００は、ｊ＝ｎｐｏｓｔｏ（ｉ）とし、ｎｐａｒｅｎｔ（ｊ）．ｎｅ．０であればｎｃｃｏｕｎｔ（ｎｐａｒｅｎｔ（ｉ））＝ｎｃｃｏｕｎｔ（ｎｐａｒｅｎｔ（ｊ））＋（ｎｃｃｏｕｎｔ（ｊ）−１）とする処理をｉ＝１〜ｎまで繰り返す（ステップＳ２２０２）。

そして、情報処理装置１００は、ｎｃｃｏｕｎｔ（ｉ）を返値としてｒｅｔｕｒｎし（ステップＳ２２０３）、計数処理を終了する。これにより、情報処理装置１００は、分枝するノードがロウサブツリーにある場合には、複数の子ノードからの伝播をキャンセルして、１つの子ノードからの特性関数の値を伝播させることができる。そして、情報処理装置１００は、ｃｏｌｕｍｎ ｃｏｕｎｔを計数することができる。

ここで、情報処理装置１００が、２つのノードのｃｏｍｍｏｎ ａｎｃｅｓｔｏｒを探索することについて説明する。情報処理装置１００は、例えば、各ノードのａｎｃｅｓｔｏｒを示す情報をｎａｎｃｅｓｔｏｒに設定する。情報処理装置１００は、各ノード自体をａｎｃｅｓｔｏｒとして初期化する。情報処理装置１００は、ポストオーダー順にノードｉを取り出して、ノードｉのａｎｃｅｓｔｏｒ（ｉ）＝ｎｐａｒｅｎｔ（ｉ）を設定する。情報処理装置１００は、ノードｉに関する列にある非零要素の行番号ｒに関して、ノードｉがｒ番目の行に対応するロウサブツリーの葉ノードか否かを判定する。

情報処理装置１００は、ノードｒのひとつ前の葉ノードｕをｎｐｒｅｖｐ（ｒ）に記憶する。情報処理装置１００は、ノードｉが葉ノードであれば、ノードｕのａｎｃｅｓｔｏｒを辿ったときにａｎｃｅｓｔｏｒがノードｉでない最後のノードを、ｃｏｍｍｏｎ ａｎｃｅｓｔｏｒのノードｉｃａとする。

以上説明したように、情報処理装置１００によれば、スパース行列Ａを対称化した対称行列Ｐのエリミネーションツリー１４０１から、スパース行列ＡをＬＵ分解した場合の下三角行列Ｌ（Ａ）や上三角行列Ｕ（Ａ）の非零要素の数を算出することができる。これにより、情報処理装置１００は、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを低減することができ、スパース行列Ａが大きくなってもＬＵ分解した結果を格納する領域を確保しやすくすることができる。また、情報処理装置１００は、ＬＵ分解において行わなくてもよい演算を省略することができ、効率よくスパース行列ＡをＬＵ分解することができる。

また、情報処理装置１００によれば、スパース行列Ａから対称行列Ｐを生成しなくても、スパース行列Ａの互いに対称位置にある要素から、エリミネーションツリー１４０１を生成することができる。これにより、情報処理装置１００は、エリミネーションツリー１４０１の生成にかかる演算を簡略化することができ、効率よくエリミネーションツリー１４０１を生成することができる。

また、情報処理装置１００によれば、スーパーノードを用いてスパース行列ＡをＬＵ分解した結果を格納する領域を算出することができる。これにより、情報処理装置１００は、スパース行列ＡをＬＵ分解した結果を格納する領域の大きさを低減することができる。

なお、本実施の形態で説明した情報処理方法は、予め用意されたプログラムをパーソナル・コンピュータやワークステーション等のコンピュータで実行することにより実現することができる。本情報処理プログラムは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行される。また本情報処理プログラムは、インターネット等のネットワークを介して配布してもよい。

上述した実施の形態に関し、さらに以下の付記を開示する。

（付記１）複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、
生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、
抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、
制御部を有することを特徴とする情報処理装置。

（付記２）前記制御部は、前記複素非対称スパース行列の対称位置にある要素の組み合わせから、前記対称行列のエリミネーションツリーを生成する、ことを特徴とする付記１に記載の情報処理装置。

（付記３）前記制御部は、前記対称行列のＬＵ分解結果において非零要素のパターンが共通する複数の列または行をまとめたスーパーノードに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、ことを特徴とする付記１または２に記載の情報処理装置。

（付記４）コンピュータが、
複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、
生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、
抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、
処理を実行することを特徴とする情報処理方法。

（付記５）コンピュータに、
複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、
生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、
抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、
処理を実行させることを特徴とする情報処理プログラム。

１００ 情報処理装置
１００１ 取得部
１００２ 算出部
１００３ 分解部

Claims (4)

1. 複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、
   生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、
   抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、
   制御部を有することを特徴とする情報処理装置。
2. 前記制御部は、前記対称行列のＬＵ分解結果において非零要素のパターンが共通する複数の列または行をまとめたスーパーノードに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、ことを特徴とする請求項１に記載の情報処理装置。
3. コンピュータが、
   複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、
   生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、
   抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、
   処理を実行することを特徴とする情報処理方法。
4. コンピュータに、
   複素非対称スパース行列から、前記複素非対称スパース行列の対称行列のエリミネーションツリーを生成し、
   生成した前記エリミネーションツリーに基づいて、前記複素非対称スパース行列の下三角行列の各行のロウサブツリーと、前記複素非対称スパース行列の上三角行列の転置行列の各行のロウサブツリーとを抽出し、
   抽出した前記下三角行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数と、抽出した前記上三角行列の転置行列の各行のロウサブツリーのうち、前記エリミネーションツリーの各ノードを含むロウサブツリーの数とに基づいて、前記複素非対称スパース行列のＬＵ分解結果を格納するメモリ領域量を決定する、
   処理を実行させることを特徴とする情報処理プログラム。

Images