(57)【特許請求の範囲】
【請求項１】
コンピュータにより、係数行列Ａ、変数ベクトルｘ、右辺ベクトルｂとして、Ａｘ＝ｂで表される連立一次方程式の解を前処理付き反復解法によって計算する計算方法であって、
前記コンピュータが、記憶装置から係数行列Ａのデータを読み込む工程と、
前記コンピュータが、係数行列Ａに対し、新たに発生した非零要素の少なくとも一部を棄却する不完全な三角分解を行うことによって、下三角行列Ｌ及び上三角行列Ｕを得る分解工程と、
前記コンピュータが、前記分解工程において分解された前記Ｌ及び前記Ｕの対角要素にスケーリングパラメタγを乗じ、非対角要素にスケーリングパラメタφを乗じることによって、下三角行列Ｌ′及び上三角行列Ｕ′を得るスケーリング工程と、
前記コンピュータが、前記Ｌ及び前記Ｕの代わりに前記Ｌ′及び前記Ｕ′を用いて、前処理付き反復解法を行う求解工程と、
を含むことを特徴とする連立一次方程式の計算方法。
【請求項２】
前記コンピュータが、前記スケーリングパラメタγ及びφを、前記係数行列Ａと前記Ｌ′及び前記Ｕ′より得られる前処理行列Ｍとより定義される評価関数を最小とするパラメタを推定することによって決定するパラメタ決定工程をさらに含むことを特徴とする請求項１に記載の連立一次方程式の計算方法。
【請求項３】
前記評価関数は、前記係数行列Ａと前記前処理行列Ｍの差行列の列毎の和を成分とするベクトル
【数１】
0005570228.tif 000007
のユークリッドノルムであることを特徴とする請求項２に記載の連立一次方程式の計算方
法。
【請求項４】
係数行列Ａ、変数ベクトルｘ、右辺ベクトルｂとして、Ａｘ＝ｂで表される連立一次方程式の解を前処理付き反復解法によって計算する数値計算装置であって、
係数行列Ａに対し、新たに発生した非零要素の少なくとも一部を棄却する不完全な三角分解を行うことによって、下三角行列Ｌ及び上三角行列Ｕを得る分解手段と、
前記分解手段によって分解された前記Ｌ及び前記Ｕの対角要素にスケーリングパラメタγを乗じ、非対角要素にスケーリングパラメタφを乗じることによって、下三角行列Ｌ′及び上三角行列Ｕ′を得るスケーリング手段と、
前記Ｌ及び前記Ｕの代わりに前記Ｌ′及び前記Ｕ′を用いて、前処理付き反復解法を行う求解手段と、
を備えることを特徴とする数値計算装置。
【請求項５】
前記スケーリングパラメタγ及びφを、前記係数行列Ａと前記Ｌ′及び前記Ｕ′より得られる前処理行列Ｍとより定義される評価関数を最小とするパラメタを推定することによって決定するパラメタ決定手段をさらに備えることを特徴とする請求項４に記載の数値計算装置。
【請求項６】
前記評価関数は、前記係数行列Ａと前記前処理行列Ｍの差行列の列毎の和を成分とするベクトル
【数２】
0005570228.tif 000008
のユークリッドノルムであることを特徴とする請求項５に記載の数値計算装置。
【請求項７】
請求項１〜３のうちいずれか１項に記載の計算方法の各工程をコンピュータに実行させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、連立一次方程式の計算方法に関し、特に、係数行列が大規模な疎行列である連立一次方程式の前処理付き反復解法に関する。
【背景技術】
【０００２】
物理現象の解析など各種の工学的な問題は、模擬的に表現された偏微分方程式を離散化することにより、最終的に多元連立一次方程式に帰着する。多くの場合、離散化の近似精度は高々数次程度なので方程式の係数は疎行列となる。一方問題を精度よく解くためには離散化の要素数を増やす必要があるため、係数行列は大規模な行列となることが多い。そのため大規模な疎行列を係数行列にもつ連立一次方程式を求解することは、非常に重要である。
【０００３】
反復解法とは連立一次方程式の求解手法の一つであり、特に係数行列が大規模疎行列の場合に最も有効な手法として知られている。具体的なアルゴリズムは、最初に近似解として適当な初期値を与え、より正確な解となるよう近似解の更新処理を繰り返すことで、徐々に真の解に近づけていくというものである。例えば正定値対称行列に特化したCG法（共役勾配法）や、非対称行列にも対応したBiCGStab法（安定化双共役勾配法）などがよく利用されている。
【０００４】
反復解法においてはその収束性、すなわちいかに少ない反復回数で近似解を真の解に十分近づけることができるか、が重要なポイントとなる。この収束性は係数行列の条件数、すなわち固有値の最大値を最小値で割った値に大きく左右され、条件数が小さいほど収束性が向上する。この性質を踏まえ、反復解法を高速化する手段として各種の前処理法が提案されている。
【０００５】
前処理の基本的な考え方は、解くべき方程式の係数行列を近似した前処理行列を用意し、その逆行列を方程式の両辺に乗じることで、方程式自体をより条件数の小さい係数行列で再構成し、反復解法における収束性を高めるというものである。実際は、方程式の再構成は負荷が大きいため陽には行わず、反復解法の反復ルーチン内に、前処理行列によるベクトルの修正処理（以降前処理と呼ぶ）を組み込む事で、再構成後の方程式に対する反復解法と数学的に同値となるアルゴリズムを構成する。
【０００６】
前処理行列に求められる性質は２つある。１つは係数行列に対する近似性が高いことである。近似性が高いほど、前処理によって係数行列が単位行列に近づき条件数が小さくなるため、反復解法の収束性が向上する。もう１つは、前処理行列の作成及び前処理に要する時間が少ないことである。前処理により収束性が向上しても、これらの処理に時間がかかり過ぎると解法全体として計算時間が増加する恐れがあるからである。一般にこの両者はトレードオフの関係にあり、バランスよく両立させるために様々な前処理法が提案されている。
【０００７】
その中でも最もポピュラーな方法の一つとして、不完全三角分解前処理（非特許文献１参照）が挙げられる。この方法は係数行列の三角分解を不完全に行い、その結果を前処理行列として用いる前処理法である。以下、その詳細について説明する。
三角分解とは、行列を下三角行列Ｌと上三角行列Ｕの積に分解するＬＵ分解と、下三角行列Ｌ、対角行列Ｄ、及び上三角行列Ｕとの積に分解するＬＤＵ分解の総称である。三角分解は、本来直接解法と呼ばれる連立一次方程式の求解手段の一種であるが、係数行列が大規模な疎行列の場合は通常使用されない。なぜならば、疎行列の三角分解には多くのｆ
ｉｌｌ−ｉｎ（フィルイン）が発生し、行列が大規模な場合に膨大な計算時間とメモリを必要とするからである。ここでｆｉｌｌ−ｉｎとは、分解中に元の行列では零であった位置へ新たに非零要素の書き込みが発生することを指す。
不完全三角分解前処理とは、係数行列の三角分解において、上記ｆｉｌｌ−ｉｎを棄却することで不完全に分解された行列を前処理行列として用いる前処理法である。前処理行列の係数行列に対する近似性は、例えば係数行列の対角成分のみを抜き出した行列を前処理行列とする対角スケーリング前処理に比べて高く、より良い収束性を示す。一方前処理行列の作成には多少の時間を要するが、ｆｉｌｌ−ｉｎを棄却するため比較的短時間で可能である。また前処理についても同様で、前処理行列がＬ＊ＵもしくはＬ＊Ｄ＊Ｕの形に分解されているため前進代入・後退代入という簡便な操作で済むことと、ＬやＵが係数行列の疎性を継承するため演算量が少ないことから、こちらも短時間で実行可能である。これらの特徴により、効率よく反復解法を高速化することができるため、最もポピュラーな前処理法の一つとして利用されている。
【０００８】
不完全三角分解前処理には、より収束性を高めるためにいくつかの派生法が提案されている。いずれも前処理行列の作成時に何らかの工夫を施すことを特徴とする。
シフトパラメタ付き不完全三角分解前処理（非特許文献１参照）と、補正パラメタ付き不完全三角分解前処理（非特許文献２参照）は、どちらも前処理行列作成の際に、元の係数行列の対角項にある変更を施してから不完全三角分解を行うことを特徴とする。前者は対角項にシフトパラメタを加算し、後者は対角項に補正パラメタを乗算する。特に係数行列の対角優位性が低い場合に、分解により対角項がゼロに近づくことを防ぎ、収束性の悪化を避けることを目的としている。
修正不完全三角分解前処理（非特許文献１参照）は、不完全三角分解においてｆｉｌｌ−ｉｎを単に棄却するのではなく、ｆｉｌｌ−ｉｎが発生した行の対角項からこの量を引くことを特徴とする。これは前処理行列が元の行列と同じ行和を持つよう強制するための操作である。実際は、ｆｉｌｌ−ｉｎに緩和パラメタを乗じた上で対角項から引くことが多く、この緩和パラメタ次第で収束性が向上することが知られている。
ｆｉｌｌ−ｉｎレベル付き不完全三角分解前処理（非特許文献１参照）は、不完全三角分解時のｆｉｌｌ−ｉｎを完全に棄却するのではなく、一定レベル以下のｆｉｌｌ−ｉｎを許容する手法である。ここでｆｉｌｌ−ｉｎのレベルとは、元の係数行列では非零であった位置をレベル０として、以降分解によってｋレベルの要素を参照するｆｉｌｌ−ｉｎはレベルｋ＋１と定義されるというものである。許容するレベルを上げるほど、不完全分解の“不完全さ”が低減するため、係数行列に対する近似性が高まり収束性が向上する。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】J.J.ドンガラ 他著／長谷川里美 長谷川秀彦 藤野清次 訳，「反復法Ｔｅｍｐｌａｔｅｓ」，朝倉書店，１９９６年，ｐｐ．５３−５７
【非特許文献２】片岡勲  安田秀辛  高野直樹  芝原正彦  著，「数値解析入門」，コロナ社，２００２年，ｐｐ．５７
【発明の概要】
【発明が解決しようとする課題】
【００１０】
以上説明した背景技術には２つの課題がある。１つは前処理による高速化の効果が低い点で、もう１つは各派生法に関し、それぞれ固有のパラメタについて、適切な値を決定する工程を持たない点である。
【００１１】
第１の課題である、前処理による高速化の効果が低い点について、その要因を以下に詳しく説明する。まず基本となる不完全三角分解前処理は、係数行列の三角分解を不完全な状態で中断し、そのまま前処理行列として利用する。しかし分解直後の前処理行列が、所
望の疎性を満たす範囲で係数行列に対し最良の近似性を持つという保障はなく、前処理による収束性の向上は限定的となる。またそれ以外の各派生法についても、基本となる不完全三角分解前処理に対して、不完全分解前および分解中に何らかの補正を施すものであり、不完全分解直後の行列をそのまま前処理行列として利用することに変わりはなく、やはり収束性の向上は限定的となる。
【００１２】
また各派生法はそれぞれ、高速化の効果を阻害する以下の要因を抱えている。まずシフトパラメタ付き及び補正パラメタ付き不完全三角分解前処理については、分解前の係数行列の対角項に変更を加えることが、前処理行列の係数行列に対する近似性を向上させる直接的な手段とはならないため、収束性の向上が限定的となる。また修正不完全三角分解については、本来計算する必要のないｆｉｌｌ−ｉｎを計算し、対角項から引くという操作を行うため、前処理行列の作成に多くの時間が必要となる。ｆｉｌｌ−ｉｎレベル付き不完全三角分解前処理についても、ｆｉｌｌ−ｉｎ許容量を増やすほど収束性はより向上するが、その代償として前処理行列の作成および前処理に要する計算時間が増加してしまう。
【００１３】
このように、上記で説明した全ての従来手法に関して、解法全体としての高速化の効果が低いという課題があり、より高速化の効果の高い前処理付き反復解法が望まれている。
【００１４】
次に第２の課題である、各派生法における固有のパラメタを適切に決定する工程を持たない点について、以下に詳しく説明する。背景技術に示した全ての派生法は、それぞれ固有のパラメタを持っている。いずれも各パラメタが高速化の度合いを大きく左右するが、係数行列毎に適切な値が異なるため注意が必要である。まずシフトパラメタ付き及び補正パラメタ付き不完全三角分解前処理については、それぞれパラメタが０及び１から離れすぎると、分解する行列が本来の係数行列とかけ離れてしまい、前処理行列の近似性が低下して収束性の悪化を招いてしまう。また修正不完全三角分解前処理についても、緩和パラメタが大きすぎると、分解により対角項がゼロに近づき収束性が悪化する事が知られている。ｆｉｌｌ−ｉｎレベル付き不完全三角分解については、ｆｉｌｌ−ｉｎ許容レベルを上げるほど収束性が向上するが、逆に前処理行列の作成および前処理に時間がかかってしまう。以上よりこれらの派生法に関しては、係数行列に応じて適切なパラメタを設定することが重要となる。しかし前処理による高速化の効果を損なわない程度の時間で、適切な値を推定する手段は明らかではない。そのためこれらの派生法は適切なパラメタを決定する工程を持たず、パラメタの決定は利用者に委ねられている。従って前処理による高速化の効果を確実に得るためには、係数行列に対していくつか適当なパラメタを用意して試行を繰り返し、適切なパラメタを手動で探索するしかない。
【００１５】
有限差分法や有限要素法といった数値解析手法においては、短い時間ステップごとに係数行列の異なる連立一次方程式を次々と解いていく必要がある。これに背景技術の各派生法を適用する場合、時間ステップごとに上記パラメタの手動探索を行うことは計算時間の点で非現実的であり、パラメタを適当な値に固定して計算を行うしかない。しかしその場合、時間ステップによって前処理の効果が変動し、最悪の場合収束が得られず計算が停止してしまう可能性もあり危険である。そのため、係数行列に応じた適切なパラメタを決定する工程を有する前処理付き反復解法が望まれている。
【００１６】
本発明は、このような背景の下になされたもので、その第１の目的は、前処理による高速化の効果が高い連立一次方程式の数値計算方法及び装置を提供することにある。また第２の目的は、第１の目的に加え、係数行列に応じて適切なパラメタを決定することが可能な連立一次方程式の数値計算方法及び装置を提供することにある。
【課題を解決するための手段】
【００１７】
上記第１の目的を達成するため、本発明は、コンピュータにより、係数行列Ａ、変数ベクトルｘ、右辺ベクトルｂとして、Ａｘ＝ｂで表される連立一次方程式の解を前処理付き反復解法によって計算する計算方法であって、前記コンピュータが、記憶装置から係数行列Ａのデータを読み込む工程と、前記コンピュータが、係数行列Ａに対し、新たに発生した非零要素の少なくとも一部を棄却する不完全な三角分解を行うことによって、下三角行列Ｌ及び上三角行列Ｕを得る分解工程と、前記コンピュータが、前記分解工程において分解された前記Ｌ及び前記Ｕの対角要素にスケーリングパラメタγを乗じ、非対角要素にスケーリングパラメタφを乗じることによって、下三角行列Ｌ′及び上三角行列Ｕ′を得るスケーリング工程と、前記コンピュータが、前記Ｌ及び前記Ｕの代わりに前記Ｌ′及び前記Ｕ′を用いて、前処理付き反復解法を行う求解工程と、を含むことを特徴とする。
【００１８】
また、上記第２の目的を達成するために、本発明は、前記コンピュータが、前記スケーリングパラメタγ及びφを、前記係数行列Ａと前記Ｌ′及び前記Ｕ′より得られる前処理行列Ｍとより定義される評価関数を最小とするパラメタを推定することによって決定するパラメタ決定工程をさらに含むことを特徴とする。この評価関数として、前記係数行列Ａと前記前処理行列Ｍの差行列の列毎の和を成分とするベクトル
【数１】
0005570228.tif 000002
のユークリッドノルムを採用することが好適である。
【発明の効果】
【００１９】
本発明によれば、前記スケーリング工程において下三角行列Ｌ及び上三角行列Ｕの対角項と非対角項のそれぞれにスケーリングパラメタを乗じることで、不完全三角分解後の前処理行列に対し、係数行列への近似性を向上させるための直接的な修正が可能となる。そのため収束性の向上が期待できる。一方前記スケーリング工程における演算量は微小であり、また前記スケーリング工程を経ても前処理行列の疎（スパース）性は保持される。そのため収束性の向上に対するトレードオフの要素、すなわち前処理行列の作成及び前処理に要する時間はほとんど増加せず、前処理による高速化の効果が高くなる。
【００２０】
また本発明によれば、前記パラメタ決定工程において、前処理の効果を損なわない程度の時間で、係数行列に応じた適切なスケーリングパラメタを推定することができる。
【００２１】
さらに背景技術のいくつかの派生法に対して本発明を適用することで、派生法固有のパラメタに対する感度を低減させ、実用性を向上させることができる。
【図面の簡単な説明】
【００２２】
【図１】本発明の数値解法を示すフローチャート
【図２】本発明の数値計算装置の構成を示す図
【図３】前処理つき共役勾配法（ＣＧ法）のアルゴリズム
【図４】本発明を適用した際のスケーリングパラメタと反復回数・計算時間の関係
【図５】従来手法と本発明との比較
【図６】従来手法に本発明を適用した場合の効果
【図７】従来手法に本発明を適用した場合の効果
【発明を実施するための形態】
【００２３】
図２は、本実施形態に係る数値計算装置として機能するコンピュータの基本構成を示すブロック図である。該コンピュータは、大略、バス２０７を介して接続されたＣＰＵ２０１、表示装置２０２，入力装置２０３、ＲＡＭ２０４、ハードディスク装置２０５、ＮＩ
Ｃ２０６等から構成される。
【００２４】
ＣＰＵ（中央演算装置）２０１で、ＲＡＭ（主記憶装置）２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。表示装置２０２は、液晶画面などにより構成されており、ＣＰＵ２０１による処理結果を文字が画像などでもって表示することができる。入力装置２０３は、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。ＲＡＭ２０４は、ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。ハードディスク装置２０５には、ＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されている。これらのプログラムはＲＡＭ２０４にロードされて、ＣＰＵ２０１によって実行されることで、本コンピュータが連立一次方程式を数値的に解くための数値計算装置として機能する。ＮＩＣ（ネットワークインターフェース）２０６は、外部装置とのデータ通信を行う際のインターフェースとして機能するものであり、本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対象としても良い。
【００２５】
次に、上記構成を備えるコンピュータが行う、連立一次方程式の解を求める為の数値解法について説明する。図１は、本実施形態に係るコンピュータが行う、連立一次方程式の解を求めるための数値計算処理のフローチャートである。なお、ここでは係数行列Ａ、変数ベクトルｘ、右辺ベクトルｂとして、Ａｘ＝ｂで表される連立一次方程式の求解処理に、本発明の前処理付き反復解法を適用する。
【００２６】
まず、係数行列Ａを不完全三角分解して下三角行列Ｌ、及び上三角行列Ｕを得る（ステップＳ１０１：分解工程）。なお、ここでは不完全三角分解としてＩＬＤＵ分解を例に説明するが、これは本発明を限定するものではなく、上記不完全三角分解はＩＬＵ分解やＩＣ（不完全コレスキー）分解であってもよい。また、背景技術に記載したシフトパラメタ付き不完全三角分解前処理、補正パラメタ付き不完全三角分解前処理、修正不完全三角分解前処理、及びｆｉｌｌ−ｉｎレベル付き不完全三角分解前処理のいずれの派生法における不完全三角分解処理であってもよい。すなわち、ここでの分解処理は、係数行列を不完全に分解し、Ｌ＊ＵもしくはＬ＊Ｄ＊Ｕの行列積を導出するいかなる処理であってよい。
【００２７】
次に、下三角行列Ｌの対角要素Ｌｉｉ及び上三角行列Ｕの対角要素Ｕｉｉに乗じるスケーリングパラメタγと、Ｌの非対角要素Ｌｉｊ（ｉ＞ｊ）及びＵの非対角要素Ｕｉｊ（ｉ＜ｊ）に乗じるスケーリングパラメタφを決定する（ステップＳ１０２）。スケーリングパラメタγ及びφは、それぞれＬ及びＵの対角項と非対角項に属する要素全てに対して一律に掛かる単一の定数である。なお、ｉ及びｊは行列の行方向及び列方向の添え字を表す。ここでのパラメタ決定方法は、決定に要する時間が、前処理による高速化の効果を損なわない程度であることが望ましい。またパラメタを作用させることで前処理行列の係数行列に対する近似性が向上するような、適切な値を導出できる方法が望ましい。そこで、スケーリングパラメタφ及びγを、係数行列Ａと前処理行列Ｍより定義される評価関数を最小とする値として決定する。この評価関数は、係数行列Ａと前処理行列Ｍの近似性が高いほど、値が小さくなる関数であることが好ましい。以下では、係数行列Ａと前処理行列Ｍとの差行列の列毎の和を成分とする下記ベクトルＲのユークリッドノルムを評価関数として、この評価関数を最小にするパラメタを推定する方法について説明する。すなわち評価関数は以下の式で表されるベクトルＲのユークリッドノルムとする。なお、nは行列の次
数を表す。
【数２】
0005570228.tif 000003
【００２８】
評価関数の最小化は、評価関数の一次微分をゼロにするγおよびφをニュートン法等の手法で数値的に解けばよいが、これは本発明を限定するものではなく、多変数関数の最適化手法のいずれを用いても良い。
【００２９】
次に、パラメタ決定工程で得られたスケーリングパラメタγｏｐｔおよびφｏｐｔを用いて、ＬおよびＵをスケーリングして、下三角行列Ｌ′と上三角行列Ｕ′を得る（ステップＳ１０３）。すなわち以下の式で表されるように、スケーリングパラメタγｏｐｔをＬの対角要素Ｌｉｉ及びＵの対角要素Ｕｉｉに乗じ、スケーリングパラメタφｏｐｔをＬの非対角要素Ｌｉｊ（ｉ＞ｊ）及びＵの非対角要素Ｕｉｊ（ｉ＜ｊ）に乗じる。
【数３】
0005570228.tif 000004
【００３０】
以上説明したステップＳ１０１からＳ１０３の操作により前処理行列Ｍ＝Ｌ′＊Ｄ＊Ｕ′が決定される。この前処理行列Ｍを用いて、以降のステップＳ１０４からＳ１０６では、前処理付き反復解法を行う。ここではＲｉｃｈａｒｄｓｏｎ反復法等の定常反復解法や、ＣＧ法、ＢｉＣＧ法、ＣＧＳ法、ＢｉＣＧＳｔａｂ法等の非定常反復解法をはじめ、その他いかなる前処理付き反復解法であってもよい。以下では前処理つきＣＧ法を例にとって説明する。前処理つきＣＧ法は公知であるため詳しい説明は省略するが、その実装法の一例を図３に示す。ただし、数学的に同値となるその他の手法によって計算しても構わない。
【００３１】
まずＣＧ法における残差ベクトルｒに対して前処理を行う（ステップＳ１０４）。すなわち、以下の式で表されるベクトルｚを、前進代入・後退代入を用いて求める。
【数４】
0005570228.tif 000005
なお、前進代入・後退代入の詳細については公知の技術であるためここでは記載しない。
【００３２】
次に、ステップＳ１０４で得られたベクトルｚを用いて、前処理付きＣＧ法における近似解の更新を行う（ステップＳ１０５）。なお、前処理付きＣＧ法の近似解更新は公知の技術であるため詳細な説明は省略するが、図３に示す近似解更新処理を行って、近似解を更新する。
【００３３】
次に、近似解の収束判定を行い、収束していれば処理を終了し、そうでなければステップＳ１０４に戻る（ステップＳ１０６）。ここでは、例えば、残差ベクトルｒのユークリッドノルムを右辺ベクトルｂのユークリッドノルムで割った値が収束判定用の閾値以下であるという収束判定条件を満たす場合に収束、そうでなければ未収束と判定する。ただし
これは本発明を限定するものではなく、前処理付き反復解法の収束判定であればいかなる方法を用いても良い。
【００３４】
なお、上記の説明では、前処理行列ＭをＬ＊ＵもしくはＬ＊Ｄ＊Ｕで構成する場合について説明したが、ここに任意のスカラー値ａを乗じ、Ｍ＝ａ＊Ｌ＊Ｕもしくはａ＊Ｌ＊Ｄ＊Ｕとして構成しても構わない。例えば、ａにスケーリングパラメタγの逆数を用いても良い。
【００３５】
また、上記の説明では、差行列の列毎の和を成分とするベクトルＲのユークリッドノルムを評価関数とし、これを最小とするスケーリングパラメタγ，φを採用しているが、その他の評価関数を用いても良い。例えば、差行列の行毎の和を成分とするベクトルのユークリッドノルムを評価関数としても良い。その場合、評価関数は以下の式で表されるベクトルＲのユークリッドノルムとなる。なお、nは行列の次数を表す。
【数５】
0005570228.tif 000006
また、スケーリングパラメタγ，φをこのような評価関数を最小とするものとして決定することが好適であるが、利用者が適当なパラメタを用意して試行を繰り返して適切なパラメタを手動で探索しても良い。
【００３６】
また、スケーリングパラメタ決定処理を簡易かつ高速なものとするために、スケーリングパラメタφを１．０で固定して、この条件の下で上記評価関数が最小となるスケーリングパラメタγを決定しても良い。この場合、スケーリング工程では、スケーリングパラメタγとφを不完全三角分解後の三角行列Ｌ，Ｕの対角要素と非対角要素に乗じていると捉えることもできるし、スケーリングパラメタγのみを三角行列Ｌ，Ｕの対角要素に乗じていると捉えることもできる。
【００３７】
以上説明した本発明の一実施形態に係る前処理付き反復解法の効果について、以下に詳しく説明する。なお求解の対象とする連立一次方程式には、流体解析で現れる３次元ポアソン方程式を有限差分法で解いた際の、ある時間ステップにおけるデータを用いた。係数行列は正定値対称行列で、規模は５４万元である。また本発明における非対角要素のスケーリングパラメタφは１．０に固定してある。なお、ベースとなる不完全三角分解にはＩＬＤＵ分解を用い、反復解法にはＣＧ法を用いる。
【００３８】
まず図４において、対角要素のスケーリングパラメタγの値を手動で変化させた場合の結果と、評価関数（ベクトルＲのユークリッドノルム）を最小にするようにスケーリングパラメタγを決定した場合の結果とを比較する。図４（Ａ）は反復回数を、図４（Ｂ）は計算時間を示すグラフである。なお、図中、白丸及び白三角がγを手動で変化させた場合、黒丸及び黒三角が評価関数を最小とするようにγを定めた場合である。
【００３９】
γを手動で変化させた結果、γ＝０．６３において反復回数と計算時間の両者が共に最小値を示した。γ＝１．０すなわち基本となるＩＬＤＵ分解前処理付きＣＧ法と比べると、反復回数は１７１回から９５回に低減し、収束性が大幅に向上している。一方計算時間も１２．７秒から７．３秒に低減しており、前処理行列の作成や前処理に要する時間はほとんど増加しないことがわかる。その結果、約１．７倍という大幅な高速化が得られている。
一方、評価関数を最小とするようにγを定めた場合はγ＝０．６４となり、手動探索による最適値に近い値が推定された。その結果、反復回数は９６回、計算時間は７．４秒となっており、手動でのパラメタ最適化と遜色ない結果が得られた。この結果は、上記のパラメタ決定方法が、本手法の高速化の効果を損なわない程度の時間で、係数行列に応じた適切なスケーリングパラメタを推定していることを示すものである。
【００４０】
次に、図５において、背景技術で説明した各種の従来手法（不完全三角分解前処理、その派生法であるシフトパラメタ付き、補正パラメタ付き、フィルインレベル付きの不完全三角分解前処理、修正不完全分解前処理）と本発明の手法とを比較する。図５（Ａ）では反復回数を、図５（Ｂ）では計算時間を比較してある。なお、本発明に関してはスケーリングパラメタを、上記のパラメタ決定手法により決定している。一方、背景技術の４つの派生法に関しては、それぞれの固有パラメタについて、様々な値を試して最も良い結果を与える最適値をあらかじめ探索しておき、そのパラメタを用いた結果を示す。図より、反復回数及び計算時間の両方で、本発明による結果が最も優れている。すなわち本発明の解法は、たった一度の計算で、何度も試行を繰り返してパラメタを最適化した背景技術の各派生法よりも優れたパフォーマンスを示し、実用性に優れることがわかる。
【００４１】
なお、下三角行列Ｌ及び上三角行列Ｕにスケーリングパラメタを乗じる本発明の手法は、不完全三角分解前処理の各派生法に対して同時に適用することができる。図６に、各派生法に対して本発明の手法を適用した場合の結果を示す。すなわち、本発明における不完全三角分解に、基本となるＩＬＤＵ分解だけでなく、各派生法で行われる独自のＩＬＤＵ分解を利用している。なお、スケーリングパラメタは上記パラメタ決定処理により評価関数を最小にするものを採用している。
【００４２】
図６（Ａ）では反復回数を、図６（Ｂ）では計算時間を示してある。また図５と同様、従来技術の４つの派生法については、固有パラメタを事前に最適化した結果を示している。図より、基本となるＩＬＤＵ分解前処理だけでなく、４つの派生法全てについて、本発明を適用する事で高速化されることがわかる。すなわち本発明はいずれの従来技術とも競合するものではなく、それらと組み合わせることでより効果を高めるという性質を持つ。
【００４３】
図６では各派生法に関し、固有パラメタを最適化した結果のみ示していたが、図７では固有パラメタに対して様々な値を試行した結果を示す。図７（Ａ）（Ｂ）がシフトパラメタ付きＩＬＤＵ分解前処理付きＣＧ法と、それに本発明を適用した場合の結果を示す。以降同様に、図７（Ｃ）（Ｄ）が補正パラメタ付きＩＬＤＵ分解、図７（Ｅ）（Ｆ）が修正ＩＬＤＵ分解、図７（Ｇ）（Ｈ）がｆｉｌｌ−ｉｎレベル付きＩＬＤＵ分解と続く。図中、白丸及び白三角は、従来手法のみによる解法での反復回数及び計算時間を示し、黒丸及び黒三角は、本発明と従来手法を組み合わせた場合の反復回数及び計算時間を示す。
【００４４】
まず全ての派生法に関して、それぞれ固有のパラメタ（シフトパラメタ、補正パラメタ緩和パラメタ、フィルインレベル）に対する感度が非常に高く、パラメタの設定が高速化の度合いを大きく左右することがわかる。一方、本発明を各派生法に適用すると、固有パラメタの値に関わらず反復回数と計算時間の両方が低減しており、固有パラメタに対するロバスト性が表れている。
【００４５】
また特に、補正パラメタ付きＩＬＤＵ分解前処理（図７（Ｃ）（Ｄ））、及び修正ＩＬＤＵ分解前処理（図７（Ｅ）（Ｆ））に本発明を適用することで、それぞれ補正パラメタと緩和パラメタに対する感度が低減していることがわかる。その結果ある範囲のパラメタが一律に最適値となっており、パラメタを固定して利用する際の安全性が増している。すなわち本発明の解法は、上記２つの派生法に対して、高速化と実用性の向上という２つの効果を与えることがわかる。
【符号の説明】
【００４６】
２０１…ＣＰＵ，２０２…表示装置，  ２０３…入力装置，  ２０４…ＲＡＭ
【図１】
0005570228.tif 000009
【図２】
0005570228.tif 000010
【図３】
0005570228.tif 000011
【図４】
0005570228.tif 000012
【図５】
0005570228.tif 000013
【図６】
0005570228.tif 000014
【図７】
0005570228.tif 000015