(57)【要約】【課題】Relaxed Supernode 融合による演算量の増加を、システム又はユーザが制御可能とし、かつ、高速なシンボリック分解手法を持つ連立一次方程式求解方法を得る。【解決手段】疎行列を係数行列とする連立一次方程式の解を求める連立一次方程式求解方法において、Supernode融合を実施した後にRelaxed Supernode融合を実施するステップと、Relaxed Supernode融合時に許容する演算量の増加率に任意の値を指定するステップと、指定された演算量増加許容率に基づいて、新たに追加できるfill-in数を算出しRelaxed Supernode融合を行うステップと、Relaxed Supernode融合によって変更された非ゼロ構造を再構築するステップとを有する。【選択図】図１
2009025962.tif 000002
【特許請求の範囲】
【請求項１】
システムパラメータを入力する処理ステップと、ユーザパラメータを入力する処理ステップと、オーダリングを行う処理ステップと、シンボリック分解を行う処理ステップと、数値分解を行う処理ステップと、前進・後退代入を行う処理ステップと、解の反復改良を行う処理ステップとを有し、係数行列が疎行列である連立一次方程式の解を計算する連立一次方程式求解方法において、
前記シンボリック分解を行う処理ステップは、入力された係数行列要素からfill-in を計算し、同一構造の節点をSupernode 融合する処理ステップと、Relaxed Supernode 融合による数値分解の演算量の増加許容率がユーザから指定されているか否かを確認し、指定されていない場合、システム既定の値を演算量の増加許容率に設定する処理ステップと、前記Supernode 融合された構造に対し前記演算量の増加許容率から算出した最大追加許容数を元にRelaxed Supernode 融合を行う処理ステップと、前記Relaxed Supernode 融合で変更された非ゼロ構造を再構築する処理ステップとを有することを特徴とする連立一次方程式求解方法。
【請求項２】
前記システムパラメータを入力する処理ステップは、ハードウェアに関するパラメータを入力する処理ステップと、Relaxed Supernode 融合による数値分解の演算量の増加許容率のシステム既定値を入力する処理ステップとを有することを特徴とする請求項１記載の連立一次方程式求解方法。
【請求項３】
前記ユーザパラメータを入力する処理ステップは、Relaxed Supernode 融合による数値分解の演算量の増加の上限値をユーザが指定する演算量の増加許容率を入力する処理ステップと、前記演算量の増加許容率でシステム既定値を上書きするか否かを選択するフラグを入力する処理ステップと、期待する解の精度を入力する処理ステップと、求解した解が期待する解の精度とならなかった場合に実行する解の反復改良の最大実行回数を示す最大反復改良回数を入力する処理ステップとを有することを特徴とする請求項１または２記載の連立一次方程式求解方法。
【請求項４】
システムパラメータを入力する手段と、ユーザパラメータを入力する手段と、オーダリングを行う手段と、シンボリック分解を行う手段と、数値分解を行う手段と、前進・後退代入を行う手段と、解の反復改良を行う手段とを有し、係数行列が疎行列である連立一次方程式の解を計算する連立一次方程式求解装置において、
前記シンボリック分解を行う手段は、入力された係数行列要素からfill-in を計算し、同一構造の節点をSupernode 融合する手段と、Relaxed Supernode 融合による数値分解の演算量の増加許容率がユーザから指定されているか否かを確認し、指定されていない場合、システム既定の値を演算量の増加許容率に設定する手段と、前記Supernode融合された構造に対し前記演算量の増加許容率から算出した最大追加許容数を元にRelaxed Supernode 融合を行う手段と、前記Relaxed Supernode 融合で変更された非ゼロ構造を再構築する手段とを有することを特徴とする連立一次方程式求解装置。
【請求項５】
システムパラメータを入力する処理と、ユーザパラメータを入力する処理と、オーダリングを行う処理と、シンボリック分解を行う処理と、数値分解を行う処理と、前進・後退代入を行う処理と、解の反復改良を行う処理とをＣＰＵに実行させることにより、係数行列が疎行列である連立一次方程式の解を計算する連立一次方程式求解プログラムにおいて、
前記シンボリック分解を行う処理は、入力された係数行列要素からfill-in を計算し、同一構造の節点をSupernode 融合する処理と、Relaxed Supernode 融合による数値分解の演算量の増加許容率がユーザから指定されているか否かを確認し、指定されていない場合、システム既定の値を演算量の増加許容率に設定する処理と、前記Supernode融合された構造に対し前記演算量の増加許容率から算出した最大追加許容数を元にRelaxed Supernode 融合を行う処理と、前記Relaxed Supernode 融合で変更された非ゼロ構造を再構築する処理とを有することを特徴とする連立一次方程式求解プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、連立一次方程式の数値解を計算する連立一次方程式求解方法及び装置に係り、特に、係数行列が疎行列である連立一次方程式を直接解法を用いて求解する連立一次方程式求解方法及び装置に関する。
【背景技術】
【０００２】
疎行列を係数行列に持つ連立一次方程式の解法の１つとして、ダイレクトスパースソルバと呼ばれる技術が知られている。一般に、ダイレクトスパースソルバにより連立一次方程式の数値解を求解する場合、その処理は、係数行列を並び替えるオーダリングの処理と、数値分解で新たに発生する非零要素（fill-in)の位置を前もって計算するシンボリック分解の処理、係数行列を三角分解する数値分解の処理、解を求める前進消去・後退代入の処理、求まった解の精度を改良する解の反復改良の処理の５つの処理に分割することができる（非特許文献１、非特許文献２参照）。
【０００３】
そして、前述した５つの処理の中で、最も処理に時間がかかるのは数値分解の処理である。この数値分解の処理を高速に行うことを可能にした手法として、シンボリック分解の処理時に、同一構造の節点同士を１つのノードとして処理するSupernode 融合と呼ばれる技術がある。Supernode 融合を行うことにより、データの局所性を高め、数値分解を高速化することができる。しかし、完全一致のSupernode 融合だけでは、数値分解の高速化には不充分であることが多い。そこで、意図的にfill-in を係数行列に追加することにより、Supernode 融合を促進するRelaxed Supernode 融合の技術が、非特許文献３等により提案されている。
【非特許文献１】小国  力編著「行列計算ソフトウェア−ＷＳ，スーパーコン，並列計算機」（丸善）
【非特許文献２】T.R.マッカーラ著「計算機のための数値計算法概論」（サイエンス社）
【非特許文献３】Cleave Ashcraft and Roger Grimes : ”The Influence of Relaxed Supernode Partitions on the Multifrontal Method” ACM Transactions on Mathematical Software, Vol.15,No.4,December 1989,Pages 291-309
【発明の開示】
【発明が解決しようとする課題】
【０００４】
前述したRelaxed Supernode 融合の従来技術は、Relaxed Supernode 融合時に追加するfill-in が、追加するfill-in の数だけでなく、追加する場所によっても数値分解の演算量を大きく変化させてしまうものである。そして、この従来技術は、Relaxed Supernode 融合による数値分解の演算量の増加をコントロール可能な高速なシンボリック分解手法について考慮されておらず、このため、多様な形状を持つ入力行列に対してRelaxed Supernode 融合を実施した場合に、数値分解でのデータの局所性が高まらない、あるいは、演算量が増大しすぎることにより数値分解での高速化を阻害してしまうという問題点を生じさせてしまう。
【０００５】
本発明の目的は、前述したような従来技術の問題点を解決し、Relaxed Supernode 融合による数値分解の演算量の増加を、システムまたはユーザが制御可能とし、数値分解処理の高速化を図ることができる連立一次方程式求解方法及び装置を提供することにある。
【課題を解決するための手段】
【０００６】
本発明によれば前記目的は、システムパラメータを入力する処理ステップと、ユーザパラメータを入力する処理ステップと、オーダリングを行う処理ステップと、シンボリック分解を行う処理ステップと、数値分解を行う処理ステップと、前進・後退代入を行う処理ステップと、解の反復改良を行う処理ステップとを有し、係数行列が疎行列である連立一次方程式の解を計算する連立一次方程式求解方法において、前記シンボリック分解を行う処理ステップは、入力された係数行列要素からfill-in を計算し、同一構造の節点をSupernode 融合する処理ステップと、Relaxed Supernode 融合による数値分解の演算量の増加許容率がユーザから指定されているか否かを確認し、指定されていない場合、システム既定の値を演算量の増加許容率に設定する処理ステップと、前記Supernode 融合された構造に対し前記演算量の増加許容率から算出した最大追加許容数を元にRelaxed Supernode 融合を行う処理ステップと、前記Relaxed Supernode 融合で変更された非ゼロ構造を再構築する処理ステップとを有することにより達成される。
【発明の効果】
【０００７】
本発明によれば、疎行列を係数行列とする連立一次方程式求解装置のシンボリック分解におけるRelaxed Supernode 融合による数値分解の演算量の増加割合をシステムまたはユーザにより制御可能となり、これにより、多様な形状の入力行列に対応したRelaxed Supernode 化の促進が可能となり、数値分解処理でのデータの局所性を高めることができるため、数値分解処理の高速化を図ることができる。
【発明を実施するための最良の形態】
【０００８】
以下、本発明による連立一次方程式求解方法及び装置の実施形態を図面により詳細に説明する。
【０００９】
図１は本発明の一実施形態による疎行列に対応した連立一次方程式求解装置の構成を示すブロック図である。
【００１０】
本発明の実施形態による連立一次方程式求解装置１００は、図１に示すように、ＣＰＵ１０１と、入力装置１０２と、出力装置１０３と、補助記憶装置１０４と、主記憶装置１０５とを有する一般的な情報処理装置内に構築されている。
【００１１】
ＣＰＵ１０１は、疎行列に対応した連立一次方程式求解装置１００全体の動作を制御する装置であり、入力装置１０２は、連立一次方程式の解を計算するための各種のデータ入力を行う装置である。出力装置１０３は、連立一次方程式の解の計算に伴う各種の処理結果の出力を行う表示装置、印刷装置等である。
【００１２】
補助記憶装置１０４は、ＨＤＤ等であり、疎行列に対応した連立一次方程式求解装置１００の動作を制御する各種処理プログラム及び入力データを格納しておく装置である。各種処理プログラムは、ＣＰＵ１０１により読み出されて、主記憶装置１０５上にロードされ実行される。補助記憶装置１０４には、Relaxed Supernode 融合による数値分解の演算量の増加の上限を指定する演算量増加許容率１０６、期待する解の精度１０７、求解した解が期待する解の精度１０７とならなかった場合に実行する解の反復改良の最大実行回数を示す最大反復改良回数１０８、求解時に主記憶装置１０５上に保持しているシステム既定値の値をユーザが指定した演算量許容増加率１０６に上書きするか否かを指定するシステム既定値設定フラグ１０９、係数行列の非ゼロ要素位置テーブル１１０、係数行列データ１１１、ユーザからの演算量許容増加率１０６が指定されなかった場合に使用されるシステム既定値１１２、ハードウェアの実行ＣＰＵ数、キャッシュサイズ、主記憶装置容量等を指定するハードウェアに関するパラメータ１１３の各データが入力データとして保持されている。
【００１３】
主記憶装置１０５は、前述した各種の処理プログラムや入力データがロードされる装置である。主記憶装置１０５にロードされる各種の処理プログラムは、システムパラメータ入力処理部１１４と、ユーザパラメータ入力処理部１１５と、オーダリング処理部１１６と、シンボリック分解処理部１１７と、数値分解処理部１１８と、前進・後退代入処理部１１９と、解の反復改良処理部１２０を構築するプログラムであり、これらの処理部は、ＣＰＵ１０１によりプログラムが実行されることにより構築される。
【００１４】
システムパラメータ入力処理部１１４及びユーザパラメータ入力処理部１１５は、補助記憶装置１０４から入力装置１０２を介して各種の入力情報の読み込みを行う。システムパラメータ入力処理部１１４及びユーザパラメータ入力処理部１１５の処理の詳細については図２により後述する。
【００１５】
シンボリック分解処理部１１７は、非ゼロ構造計算及びSupernode 融合処理部１２１とRelaxed Supernode 融合処理部１２２と非ゼロ構造再構築処理部１２３とにより構成される。非ゼロ構造計算及びSupernode 融合処理部１２１は、オーダリング処理部１１６で並べ替えられた係数行列の構造に対して、数値分解処理で発生するfill-in 要素の位置を特定し、構造が一致する節点同士をSupernode 融合する処理を行う。
【００１６】
Relaxed Supernode 融合処理部１２２は、非ゼロ構造計算及びSupernode 融合処理部１２１での処理で得られた非ゼロ構造に対して適切な位置にfill-in を追加し、構造が一致した節点同士をRelaxed Supernode 融合する処理を行う。Relaxed Supernode 融合処理部１２２は、最大追加許容数算出部１２４と、消去木作成部１２５と、消去木トレース処理部１２６と、Relax 化調査部１２７とにより構成され、Relax 化調査部１２７は、差分構造解析部１２８、差分構造融合部１２９、マージテーブル更新部１３０により構成される。このRelaxed Supernode 融合処理部１２２内の処理の詳細については、図４、図５、図６、図７、図８、図９を参照して後述する。
【００１７】
非ゼロ構造再構築処理部１２３は、Relaxed Supernode 融合処理部１２２で変更となった非ゼロ構造の再構築の処理を行う。この再構築の処理では、非ゼロ構造計算及びSupernode 融合処理部１２１で得られている非ゼロ構造を再利用して再構築を行うため、再構築にかかる処理時間は、非ゼロ構造計算及びSupernode 融合処理部１２１の処理時間と比較して無視できるものである。
【００１８】
数値分解処理部１１８は、シンボリック分解処理部１１７での処理結果をもとに、係数行列要素に対して数値分解を実施する。前進・後退代入処理部１１９は、数値分解処理部１１８での処理で数値分解された係数行列要素と、連立一次方程式の右辺ベクトルから前進消去、後退代入の処理を行って連立一次方程式の解を求める。また、解の反復改良処理部１２０は、前進・後退代入処理部１１９での処理で求められた解の精度が期待する解の精度１０７を満足するか否かをチェックし、充分な精度でなかった場合に、最大反復改良回数１０８を上限として解の反復改良を行う。
【００１９】
図２はシステムパラメータ入力処理部１１４及びユーザパラメータ入力処理部１１５での処理動作を説明するフローチャートと、その処理で作成されるパラメータ保持テーブルとを示す図である。まず、図２（ａ）を参照して、システムパラメータ入力処理部１１４での処理動作を説明する。
【００２０】
（１）システムパラメータ入力処理部１１４での処理が開始されると、システムパラメータ入力処理部１１４は、補助記憶装置１０４からシステム既定値１１２と、ハードウェアに関するパラメータ１１３を読み出す（ステップ２０１）。
【００２１】
（２）次に、システムパラメータ入力処理部１１４は、図２（ｃ）にパラメータ保持テーブル２０３として示しているテーブルのITB(2)２０５の位置に、ステップ２０１の処理で読み出したシステム既定値１１２を格納して、処理を終了する（ステップ２０２）。
【００２２】
次に、図２（ｂ）を参照して、ユーザパラメータ入力処理部１１５での処理動作を説明する。
【００２３】
（１）ユーザパラメータ入力処理部１１５での処理が開始されると、ユーザパラメータ入力処理部１１５は、補助記憶装置１０４から演算量増加許容率１０６、期待する解の精度１０７、最大反復改良回数１０８、システム既定値設定フラグ１０９、非ゼロ要素位置テーブル１１０、係数行列データ１１１を読み出す（ステップ２０６）。
【００２４】
（２）次に、パラメタ保持テーブル２０３のITB(1)２０４の位置にステップ２０６の処理で読み出した演算量増加許容率１０６を格納し、また、読み出したシステム既定値設定フラグ１０９が１となっているか否かを判定する（ステップ２０７、２０８）。
【００２５】
（３）ステップ２０８の判定で、システム既定値設定フラグ１０９が１となっていなかった場合、ここでの処理を終了し、システム既定値設定フラグ１０９が１となっていた場合、パラメータ保持テーブル２０３のITB(2)２０５の位置に演算量増加許容率１０６を格納する。すなわち、システムパラメータ入力処理部１１４でのステップ２０２の処理で保存したシステム既定値１１２を演算量増加許容率１０６で上書きする（ステップ２０９）。
【００２６】
図３は本発明の実施形態で取り扱う係数行列の形状及びその入力書式の例を示す図である。
【００２７】
図３（ａ）に示す行列３０１は、連立一時方程式の係数行列を示しており、次元数８（Ｎ＝８）、非ゼロ要素数２３（ＮＺ＝２３）の非対称疎行列の非ゼロ要素の位置関係を、非ゼロ要素を黒丸●３０２で示したものである。なお、解を求めることのできる連立一時方程式の係数の配列は、係数行列の１つの対角線上には、全て非ゼロ要素が配置された配列となる。このような行列を入力する場合、ＬＣ（Ｎ＋１）３０６として示すような一次元配列と、ＬＬ（ＮＺ）３０７として示すような一次元配列とを作成し、図３（ｂ）に示すような非ゼロ要素位置テーブル１１０として入力される。
【００２８】
一次元配列ＬＣ（Ｎ＋１）３０６の配列要素は、行列３０１の各行または各列の先頭の非ゼロ要素の列または行番号のインデックスを格納して作成したものである。なお、この一次元配列ＬＣ（Ｎ＋１）３０６のＮ＋１の位置には非ゼロ要素の全要素数＋１の値が格納される。図示例の場合、行列３０１の各行の先頭の非ゼロ要素が行順に数えたときに何番目の非ゼロ要素であるかを示す値が順に格納され、また、Ｎ＋１の位置、説明している例では、次元数８であるので９番目の位置には、非ゼロ要素の数に１を加えた２４が格納されている。この一次元配列ＬＣ（Ｎ＋１）３０６において、隣り合う値の右側の値から左側の値を引くと、左側の値が格納されている位置に対応する行の中に存在する非ゼロ要素の数となる。
【００２９】
また、一次元配列ＬＬ（ＮＺ）３０７の配列要素は、行列３０１の非ゼロ要素を行または列方向順に列または行番号のインデックスを格納して作成したものである。図示例の場合、行列３０１の行順に、その行に存在する非ゼロ要素の列番号が格納されている。
【００３０】
行列３０１中の白丸○３０３は、本来ゼロである要素が数値分解の処理時に非ゼロとなるfill-in 要素を示している。
【００３１】
行列３０１の中の太い黒線で囲んだ領域３０４は、その上三角部分の７行目と８行目、下三角部分の７列目と８列目の係数行列の構造が一致し、７段目が８段目にSupernode 融合された状態を示している。また、行列３０１の中の太い黒線で囲んだ領域３０５は、fill-in 要素が確定することにより、前述した対角線を挟んだ上側部分の４行目と５行目の係数行列の構造と、下側部分の４列目と５列目の係数行列の構造とが一致し、４段目が５段目にSupernode 融合された状態を示している。
【００３２】
図３（ｃ）に示す一次元配列ＭＴＢ（Ｎ）３０８は、Supernode 融合の融合情報を保持するテーブルである。ＭＴＢ（Ｎ）３０８の各配列要素１〜Ｎには行列３０１の各段でのSupernode 融合先の節点番号が格納される。図示例の場合、行列３０１の段（行列の非ゼロ要素が連なる対角線上で左上から数えた値）の１段目から３段目までは、前述で説明したようなSupernode 融合がされていないので、段番号が順に１、２、３と格納され、４段目、５段目は、前述で説明したようにSupernode 融合がされているので、５段側に融合されるとして、５、５が格納される。また、６段目は、Supernode 融合がされていないので、６が格納され、７段目、８段目は、前述で説明したようにSupernode 融合がされているので、８段側に融合されるとして、８、８が格納される。
【００３３】
図４はRelaxed Supernode 融合処理部１２２での処理動作を説明するフローチャート、及び、この処理で作成される消去木テーブル４０６の構成を示す図であり、次に、これについて説明する。
【００３４】
（１）図４（ａ）に示すように、Relaxed Supernode 融合処理部１２２での処理が開始されると、まず、最大追加許容数算出部１２４が、Relaxed Supernode 融合時に追加するfill-in の最大追加許容数を算出する。この最大追加許容数算出部１２４での処理の詳細は図５により後述する（ステップ４０１）。
【００３５】
（２）次に、消去木作成部１２５が、前述した非ゼロ構造計算及びSupernode 融合処理部１２１で確定した分解後の係数行列の下三角行列及び上三角行列の非ゼロ構造を元に消去木を作成する。消去木とは、Ａｉｊをｉ行ｊ列の行列要素としたとき、第ｋ列の親を、Ｍｉｎ（ｊ：ｋ＜ｊ，Ａｋｊ≠０）、あるいは、Ｍｉｎ（ｉ：ｋ＜ｉ，Ａｉｋ≠０）として定義する木構造である。ここで、第ｋ列の親となる節点が存在しない場合、［係数行列の次元数Ｎ＋１］を親節点とし、これがroot節点となる。また、Supernode 融合されている節点に関しては、Supernode の起点となっている節点番号のみを使用して消去木を作成する（ステップ４０２）。
【００３６】
前述の処理で作成される消去木テーブル４０６の例を図４（ｂ）に示しており、この例は、前述した行列３０１の例を元に、その下三角行列を用いて作成される消去木のデータを保持するテーブルの具体例を示している。この消去木テーブル４０６において、配列 ETB４０７は、［係数行列の次元数Ｎ＋１］だけ確保されるため、行列３０１の例の場合の配列 ETB４０７のサイズは９となる。配列 ETB４０７の各要素には、構造体フォーマット４０８として示される構造体へのポインタが格納され、ポインタがNULLの場合、その節点は子供節点を持たないことを示している。構造体フォーマット４０８は、節点番号４０９と、次へのポインタ４１０との２つの構成要素からなり、次へのポインタ４１０がNULLのときには、後に続く子供節点がないことを示している。
【００３７】
図４（ｂ）に示す消去木テーブル４０６の例では、配列 ETB４０７の節点５には、後に続く子供節点３があり、さらに、その節点３には、後に続く子供節点１があることを示しており、また、配列 ETB４０７の節点６には、後に続く子供節点２があることを示している。また、配列 ETB４０７の節点８には、子供節点６があり、さらに、その節点６には、後に続く子供節点５があることを示しており、配列 ETB４０７のルートノードとしての節点９には、後に続く子供節点８があることを示している。
【００３８】
（３）次に、消去木トレース処理部１２６が、ステップ４０２の処理で消去木作成部１２５が作成した消去木から、root節点が持つ全ての子供節点に対して、消去木トレースの処理を繰り返し実行し、root節点の子供節点を全て探索したか否かを判定し、全て探索した場合、子供節点探索ループ処理を終了して、Relaxed Supernode 融合処理部１２２での処理を終了する（ステップ４０３〜４０５）。
【００３９】
ステップ４０３〜４０５の処理において、消去木テーブル４０６の場合、root節点は「９」であるため、配列 ETB（９）のポインタが指し示す構造体４１１の節点番号「８」に対してステップ４０４の消去木トレース処理部１２６を実行する。節点番号「８」に対する消去木トレース処理部１２６の実行が完了した後、構造体４１１の次へのポインタが指し示す構造体に対して、同様の処理を実行する。但し、この例の場合、NULLであるため、ステップ４０５で子供節点探索ループ処理を終了する。
【００４０】
図５は図４（ａ）に示して説明したステップ４０１での最大追加許容数算出部１２４の処理の詳細を説明するフローチャート、及び、MAX_ZEROの算出方法を説明する図であり、次に、これについて説明する。
【００４１】
（１）図５（ａ）に示すように、最大追加許容数算出部１２４の処理が開始されると、最大追加許容数算出部１２４は、まず、図２に示したパラメータ保持テーブル２０３のITB(1)２０４に保持されている演算量増加許容率１０６の値が０以上であるか否か、ここでは、演算量増加許容率１０６の値として何らかの値が指定されているか否かを判定する（ステップ５０１）。
【００４２】
（２）ステップ５０１の判定で、ITB(1)２０４の値≧０、すなわち、何らかの値が指定されていた場合、演算量増加許容率をα（％）としたとき、αにユーザ指定の演算量増加許容率ITB(1)２０４を設定し、ユーザからの演算量増加許容率の指定がなかった場合、システムで前もって指定した演算量増加許容率のシステム既定値ITB(2)２０５をαに設定する（ステップ５０２、５０３）。
【００４３】
（３）次に、ステップ５０２あるいはステップ５０３の処理で設定したαを元に、Relaxed Supernode 融合時に追加するfill-in の最大追加許容数を算出し、ここでの処理を終了する。ここで、各段での最大追加許容数をMAX_ZEROとしたとき、MAX_ZEROには、前述の非ゼロ構造計算及びSupernode 融合処理部１２１で確定した分解後の係数行列から算出した平均帯幅とα／100 を乗じた値を設定する。各段において、このMAX_ZEROの範囲内でfill-in の追加を許容することにより、全体の演算量の増加を最大でα％以内とすることができる（ステップ５０５）。
【００４４】
図５（ｂ）に示すMAX_ZEROの算出の例５０５は、前述で説明した行列３０１を例に、ユーザからの演算量増加許容率１０６に６０（％）が指定されていた場合の、MAX_ZEROの算出の具体例を示したものである。MAX_ZEROは、α及び分解後の係数行列の平均帯幅から算出される。α＝６０（％）、行列３０１の分解後の係数行列の平均帯幅は３．３７５であるため、最大追加許容数  MAX_ZERO＝０．６×３．３７５≒２となる。なお、行列３０１の分解後の係数行列の平均帯幅は、図３（ａ）に示す係数行列の次元数が８であり、その行列に含まれる非ゼロ要素の数とfill-in 要素の数とを加えた数が２７であるため、２７／８＝３．３７５として求められたものである。
【００４５】
図６は図４（ａ）に示して説明したステップ４０４での消去木トレース処理部１２６の処理の詳細を説明するフローチャート、及び、図４（ｂ）に示して説明した消去木テーブル４０６の構造を木構造に書き換えて示した図である。
【００４６】
（１）図６（ａ）に示すように、消去木トレース処理部１２６の処理が開始されると、消去木トレース処理部１２６は、まず、図３（ｃ）に示して説明した一次元配列ＭＴＢ（Ｎ）３０８を参照し、当該節点がSupernode 融合されているか否かを確認する。この確認は、ＭＴＢ（当該節点番号）＝＝当該節点番号であるか否かを判定する処理である（ステップ６０１）。
【００４７】
（２）ステップ６０１の判定で、ＭＴＢ（当該節点番号）＝＝当該節点番号であった場合、当該節点はSupernode の起点であるか、あるいは、Supernode 融合されていない節点であるため、Relaxed Supernode 融合の対象となる節点を検索するために、Relax 化調査部１２７の処理を実行する。この処理の詳細については、図７により後述する（ステップ６０２）。
【００４８】
（３）ステップ６０２でのRelax 化調査部１２７の処理の終了後、あるいは、ステップ６０１の判定で、ＭＴＢ（当該節点番号）≠当該節点番号であった場合（この場合、当該節点はSupernode に融合済みの節点であり、Relax 化の調査の対象とならない）、図４（ａ）に示したフローのステップ４０２の処理で消去木作成部１２５が作成した消去木から、調査対象の当該節点番号が持つ全ての子供節点に対して、消去木トレース処理部１２６による処理を再帰的に実行し、当該節点の子供節点を全て探索したことを判定し、子供節点探索ループの処理を終了し、また、ここでの処理を終了する（ステップ６０３〜６０５）。
【００４９】
前述した消去木トレース処理部１２６での処理が、図３（ａ）に示した行列３０１の例において当該節点番号（８）で実行された場合、ＭＴＢ（Ｎ）３０８から、ＭＴＢ（８）==８であるため、ステップ６０１の判定からステップ６０２の処理に分岐し、節点番号（８）に対してRelax 化調査部１２７の処理が実行される。ステップ６０３の子供節点探索ループの処理では、図４（ｂ）に示した消去木テーブル４０６を参照し、配列 ETB４０７の要素（８）のポインタが指し示す構造体４１２の節点番号（６）に対してステップ６０４の消去木トレース処理部１２６の処理を実行する。節点番号（６）に対する消去木トレース処理部１２６の実行が完了後、構造体４１２の次へのポインタが指し示す構造体４１３の節点番号（５）に対して、同様の処理を次へのポインタがNULLとなるまで実行する。
【００５０】
図６（ｂ）に示す消去木６０６は、図４（ｂ）に示して説明した消去木テーブル４０６の構造を木構造に書き換えたものである。Supernode 融合されている節点に関しては、Supernode の起点となっている節点番号のみを使用して消去木テーブル４０６を作成しているため、節点番号（４）及び、節点番号（７）は消去木上から省略されている。消去木トレース処理部１２６の処理は、前述したような処理を各節点に対して実行することにより、節点番号（８）６０７、節点番号（６）６０８、節点番号（２）６０９、節点番号（５）６１０、節点番号（３）６１１、節点番号（１）６１２の順に消去木をトレースしている。
【００５１】
図７は図６に示して説明したステップ６０２でのRelax 化調査部１２７の処理の詳細を説明するフローチャートであり、次に、これについて説明する。
【００５２】
（１）Relax 化調査部１２７での処理が開始されると、まず、ステップ７０１の子供節点探索ループにより、図４（ａ）に示すフローのステップ４０２の処理で作成した消去木を元に、調査対象の節点番号を親節点としたとき、親節点が持つ全ての子供節点に対して、以降のステップ７０２〜７１１の処理を繰り返し実行することを設定する（ステップ７０１、７１２）。
【００５３】
（２）その後、まず、追加するfill-in 数をADD_FILLとしたとき、ADD_FILLを０で初期化し、差分構造解析部１２８が、調査対象の節点番号が持つ係数行列の下三角構造と、子供節点が持つ係数行列の下三角構造との差分を解析し、子供節点が持つ係数行列の下三角構造に対して差分要素を何個追加することにより構造が一致するかを求める（ステップ７０２、７０３）。
【００５４】
（３）次に、ステップ７０３の処理で求めた差分要素数をADD_FILLに加算して、新たなADD_FILLとしての値を求め、求めた新たなADD_FILLと、図４（ａ）に示すフローのステップ４０１の処理で求めたMAX_ZEROとを比較して、MAX_ZERO≧ADD_FILLであるか否かを判定する（ステップ７０４、７０５）。
【００５５】
（４）ステップ７０５の判定で、ADD_FILLに対してMAX_ZEROの方が小さかった場合、指定された演算量増加許容率１０６あるいはシステム既定値１１２の範囲内ではRelaxed Supernode 融合が不可であるため、次の子供節点に調査対象を移して、ステップ７０２からの処理を繰り返す（ステップ７１２）。
【００５６】
（５）ステップ７０５の判定で、ADD_FILLがMAX_ZERO以下であった場合、差分構造解析部１２５は、調査対象の節点番号が持つ係数行列の上三角構造と、子供節点が持つ係数行列の上三角構造との差分を解析し、子供節点が持つ係数行列の上三角構造に対して差分要素を何個追加することにより構造が一致するかを求める（ステップ７０６）。
【００５７】
（６）次に、ステップ７０６の処理で求めた差分要素数をADD_FILLに加算して、新たなADD_FILLとしての値を求め、求めた新たなADD_FILLと、図４（ａ）に示すフローのステップ４０１の処理で求めたMAX_ZEROとを比較して、MAX_ZERO≧ADD_FILLであるか否かを判定する（ステップ７０７、７０８）。
【００５８】
（７）ステップ７０８の判定で、ADD_FILLに対してMAX_ZEROの方が小さかった場合、指定された演算量増加許容率１０６あるいはシステム既定値１１２の範囲内ではRelaxed Supernode 融合が不可であるため、次の子供節点に調査対象を移して、ステップ７０２からの処理を繰り返す（ステップ７１２）。
【００５９】
（８）ステップ７０８の判定で、ADD_FILLがMAX_ZERO以下であった場合、子供節点と親節点とをRelaxed Supernode 融合するため、差分構造融合部１２９が、子供節点と親節点との差分構造の融合を行う（ステップ７０９）。
【００６０】
（９）次に、マージテーブル更新部１３０が、融合された子供節点の節点番号をｋとしたとき、図３（ｃ）に示すＭＴＢ（ｋ）を親節点の節点番号により更新し、その後、Relax 化調査部１２７が、融合された子供節点が持つ子供節点と、親節点とをRelaxed Supernode 融合することができるか否かの調査を再帰的に実行する（ステップ７１０、７１１）。
【００６１】
図８は図３（ａ）に示した行列３０１の例において当該節点番号（８）でRelax 化調査部１２７の処理が実行された場合に、差分構造解析部１２８及び差分構造融合部１３０で用いられるテーブル構造の具体例を示した図である。
【００６２】
親節点となる節点番号８の下三角構造テーブル８０１は、図８（ａ）に示すように、係数行列の次元数Ｎ個の配列NTBLが確保されて構成され、行列３０１の例では、NTBL(８)となる。配列NTBLには、節点番号（８）が持つ下三角要素の番号が配列の先頭要素からNTBL(１)＝７、NTBL(２)＝８が格納され、残りの要素はゼロで初期化されている。また、子供節点である節点番号６の下三角構造テーブル８０２は、当該子供節点が持つ下三角構造の要素数個の配列KTBLが確保されて構成され、行列３０１の例では、KTBL(２)となり、各要素には節点番号（６）が持つ下三角要素の番号が配列の先頭要素からKTBL(１)＝６、KTBL(２)＝７が格納されている。
【００６３】
前述で説明した図７に示すフローのステップ７０３の差分構造解析部１２８での処理の実行（下三角）では、配列KTBL(１)以外に格納されている要素と、配列NTBLに格納されているゼロ以外の要素との比較を順次実行し、差分要素数を算出する。
【００６４】
また、親節点となる節点番号８の上三角構造テーブル８０３は、図８（ｂ）に示すように、下三角構造テーブル８０１の場合と同様に配列NTBUが確保されて構成され、子供節点である節点番号６の上三角構造テーブル８０４も、下三角構造テーブル８０２の場合と同様に、配列KTBUが確保されて構成される。
【００６５】
そして、図７に示すフローのステップ７０６の差分解析部１２８での処理の実行（上三角）においても、同様の手順で節点番号８の上三角構造テーブル８０３と、節点番号６の下三角構造テーブル８０４から差分要素数を算出する。
【００６６】
前述した図７に示すフローのステップ７０９の差分構造融合部１２９での処理の実行では、親の構造に対して調査した子供の節点番号を付加する処理を行う。具体的には、下三角構造に関しては、節点番号８の下三角構造テーブル８０１に格納されているNTBL(２)をNTBL(３)に、NTBL(１)をNTBL(２)に格納し、節点番号６の下三角構造テーブル８０２のKTBL(１)の要素をNTBL(１)に格納することにより、融合後の下三角構造テーブル８０５を得る処理を行う。また、上三角構造に関しても、節点番号８の上三角構造テーブル８０３に格納されているNTBU(２)をNTBU(３)に、NTBU(１)をNTBU(２)に格納し、節点番号６の上三角構造テーブル８０４のKTBU(１)の要素をNTBU(１)に格納することにより、融合後の上三角構造テーブル８０６を得る処理を行う。
【００６７】
融合後の前述したテーブル８０５、８０６は、親となって、さらに、この子供の調査が行われることになる。
【００６８】
図９はユーザから演算量増加許容率１０６として６０％が指定されていた場合のRelax 後の係数行列の具体例を説明する図である。
【００６９】
図９（ａ）に示す係数行列９０１は、前述した行列３０１の例に対してRelax 処理を行った後の係数行列であり、この行列９０１において、消去木６０６の節点番号８（６０７）が持つ構造９０２に着目した場合、節点番号６（６０８）が持つ上三角構造との差分が黒星で示す位置９０３に１要素、下三角構造との差分が黒星で示す位置９０４に１要素の合計２つのfill-in を追加することにより上三角と下三角との構造が一致し、Relaxed Supernode 融合の対象となる。前述したMAX_ZEROの算出５０５から、α＝６０％の場合、MAX_ZERO＝２であるため、この追加fill-in 数は、MAX_ZEROの範囲内であり、節点番号６と節点番号８との構造がRelaxed Supernode 融合される。
【００７０】
図９（ｂ）、図９（ｃ）に示す部分行列９０５及び部分行列９０６は、Relaxed Supernode 融合前の節点番号６、節点番号８がそれぞれ持つ部分行列である。また、図９（ｄ）に示す部分行列９０７は、節点番号６と節点番号８とをRelaxed Supernode 融合した場合に節点番号８が持つ部分行列である。
【００７１】
前述したように、黒星の要素９０３黒星の要素９０４とを追加し、節点番号６と節点番号８の構造をRelaxed Supernode 融合することにより、数値分解処理において、節点番号６が持つ部分行列９０５及び節点番号８が持つ部分行列９０６のそれぞれが必要としていた数値分解処理を、部分行列９０７での１度の数値分解処理で求めることができることになる。
【００７２】
一般に、数値分解の高速化手法の１つにブロック化アルゴリズムがあるが、入力する係数行列が疎行列の場合、分解後の非ゼロ構造が疎であるため、適切なブロックサイズで処理できないことが多い。
【００７３】
前述した本発明の実施形態は、増加演算量を制御しつつRelaxed Supernode 融合を行い、より適切なブロックサイズでの処理を行うことができるため、数値分解の高速化を図ることができる。
【００７４】
前述した本発明の実施形態での各処理は、プログラムにより構成し、本発明が備えるＣＰＵに実行させることができ、また、それらのプログラムは、ＦＤ、ＣＤＲＯＭ、ＤＶＤ等の記録媒体に格納して提供することができ、また、ネットワークを介してディジタル情報により提供することができる。
【００７５】
以上説明したように本発明の実施形態によれば、疎行列を係数行列とした連立一次方程式求解装置のシンボリック分解におけるRelaxed Supernode 融合による数値分解の演算量の増加割合をシステムまたはユーザが制御することができるので、多様な形状の入力行列においてもRelaxed Supernode 化を促進することができ、データの局所性が高まり、数値分解処理にかかる時間を削減し、連立一次方程式求解装置の性能の向上を図ることができる。
【図面の簡単な説明】
【００７６】
【図１】本発明の一実施形態による疎行列に対応した連立一次方程式求解装置の構成を示すブロック図である。
【図２】システムパラメータ入力処理部及びユーザパラメータ入力処理部での処理動作を説明するフローチャートと、その処理で作成されるパラメータ保持テーブルとを示す図である。
【図３】本発明の実施形態で取り扱う係数行列の形状及びその入力書式の例を示す図である。
【図４】Relaxed Supernode 融合処理部１２２での処理動作を説明するフローチャート、及び、この処理で作成される消去木テーブル４０６の構成を示す図である。
【図５】図４（ａ）に示して説明したステップ４０１での最大追加許容数算出部１２４の処理の詳細を説明するフローチャート、及び、MAX_ZEROの算出方法を説明する図である。
【図６】図４（ａ）に示して説明したステップ４０４での消去木トレース処理部１２６の処理の詳細を説明するフローチャート、及び、図４（ｂ）に示して説明した消去木テーブル４０６の構造を木構造に書き換えて示した図である。
【図７】図６に示して説明したステップ６０２でのRelax 化調査部１２７の処理の詳細を説明するフローチャートである。
【図８】図３（ａ）に示した係数行列の例において当該節点番号（８）でRelax 化調査部１２７の処理が実行された場合に、差分構造解析部１２８及び差分構造融合部１３０で用いられるテーブル構造の具体例を示した図である。
【図９】ユーザから演算量増加許容率として６０％が指定されていた場合のRelax 後の係数行列の具体例を説明する図である。
【符号の説明】
【００７７】
１００  連立一次方程式求解装置
１０１  ＣＰＵ
１０２  入力装置
１０３  出力装置
１０４  補助記憶装置
１０５  主記憶装置
１０６  演算量増加許容率
１０７  期待する解の精度
１０８  最大反復改良回数
１０９  システム既定値設定フラグ
１１０  非ゼロ要素位置テーブル
１１１  係数行列データ
１１２  システム既定値
１１３  ハードウェアに関するパラメータ
１１４  システムパラメータ入力処理部
１１５  ユーザパラメータ入力処理部
１１６  オーダリング処理部
１１７  シンボリック分解処理部
１１８  数値分解処理部
１１９  前進・後退代入処理部
１２０  解の反復改良処理部
１２１  非ゼロ構造計算及びSupernode 融合処理部
１２２  Relaxed Supernode 融合処理部
１２３  非ゼロ構造再構築処理部
１２４  最大追加許容数算出部
１２５  消去木作成部
１２６  消去木トレース処理部
１２７  Relax 化調査部
１２８  差分構造解析部
１２９  差分構造融合部
１３０  マージテーブル更新部
【図１】
2009025962.tif 000003
【図２】
2009025962.tif 000004
【図３】
2009025962.tif 000005
【図４】
2009025962.tif 000006
【図５】
2009025962.tif 000007
【図６】
2009025962.tif 000008
【図７】
2009025962.tif 000009
【図８】
2009025962.tif 000010
【図９】
2009025962.tif 000011