(57)【要約】【課題】多重格子法を用いて連立一次方程式を求解する場合に、格子の補間処理が不要な範囲でより多くの格子分割数を設定できるようにする。【解決手段】設定された格子の分割数を素因数分解し、その分解に係る素数に基づいて、連立一次方程式に係る線型方程式の係数、変数に対して粗視化作用素、精細化作用素子を作用させる。【選択図】図２
2007286738.tif 000002
【特許請求の範囲】
【請求項１】
多重格子法を用いて連立一次方程式を求解する数値計算装置であって、
格子の分割数を素因数分解する分解手段と、
前記分解に係る素数に対応する粗視化作用素を用いて連立一次方程式に係る線型方程式を粗視化する粗視化手段と、
前記分解に係る素数に対応する精細化作用素を用いて連立一次方程式に係る線型方程式を精細化する精細化手段と、
を有することを特徴とする数値計算装置。
【請求項２】
更に、前記分解に係る素数を小さい順に並べた素数ソーティングリストを作成する作成手段を有し、前記粗視化手段及び前記精細化手段は、それぞれ、該素数ソーティングリストから所定の順序で素数を１つずつ抽出して該素数に対応する粗視化作用素、精細化作用素を用いて前記粗視化、精細化を順次行なうことを特徴とする請求項１に記載の数値計算装置。
【請求項３】
前記粗視化手段は、前記素数ソーティングリストの最初の方から順に素数を１つずつ抽出して該素数に対応する粗視化作用素を用いて前記粗視化を順次行ない、前記精細化手段は、前記素数ソーティングリストの最後の方から順に素数を１つずつ抽出して該素数に対応する精細化作用素を用いて前記精細化を順次行なうことを特徴とする請求項２に記載の数値計算装置。
【請求項４】
前記作成手段は、所定の数値より大きな素数は、前記素数ソーティングリストにリストアップすることなく、該所定の数値より大きな素数を含む因数をリストアップすることを特徴とする請求項２又は３に記載の数値計算装置。
【請求項５】
前記分解手段は、多次元空間上の各次元に係る格子の分割数をそれぞれ素因数分解することを特徴とする請求項１〜４の何れかに記載の数値計算装置。
【請求項６】
前記作成手段は、多次元空間上の各次元に係る格子の分割数に対応する各素数ソーティングリストの要素の数を、前記分解に係る素数の数が最多数の素数ソーティングリストにおける素数の数に統一することを特徴とする請求項２〜５の何れかに記載の数値計算装置。
【請求項７】
前記作成手段は、前記分解に係る素数の数が最多数に満たない場合は、当該素数ソーティングリストに１を補充して前記統一を図ることを特徴とする請求項６に記載の数値計算装置。
【請求項８】
更に、前記粗視化された線型方程式及び前記精細化された線型方程式を反復法を用いて求解する求解手段を有することを特徴とする請求項１〜７の何れかに記載の数値計算装置。
【請求項９】
多重格子法を用いて連立一次方程式を求解する数値計算方法であって、
格子の分割数を素因数分解する分解工程と、
前記分解に係る素数に対応する粗視化作用素を用いて連立一次方程式に係る線型方程式を粗視化する粗視化工程と、
前記分解に係る素数に対応する精細化作用素を用いて連立一次方程式に係る線型方程式を精細化する精細化工程と、
を有することを特徴とする数値計算方法。
【請求項１０】
請求項９に記載の数値計算方法をコンピュータに実行させるコンピュータ読み取り可能なプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値計算装置、数値計算方法、及びプログラムに関し、特に、係数行列が大規模な連立一次方程式の解法技術に関する。
【背景技術】
【０００２】
複雑な偏微分方程式で表された各種の工学的な応用問題も、最終的には多元連立一次方程式で表されるため、連立一次方程式を高精度に求解することは、非常に重要である。従来、大規模な連立一次方程式の解法として反復法が有効であることが知られている。特に、零の要素が多いスパース（疎）な係数行列には反復法が有効である。
【０００３】
しかし、反復法として通常用いられるＧａｕｓｓ−Ｓｅｉｄｅｌ法やＳＯＲ法では、定義域を離散化して得られる格子のナイキスト周波数の近傍の周波数成分の誤差については、誤差の減少速度が速くなるが、それ以外の周波数成分の誤差では、誤差の減少速度が遅くなってしまう。このため、Ｇａｕｓｓ−Ｓｅｉｄｅｌ法やＳＯＲ法では、十分満足できる解へ収束するまでに長時間を要していた。
【０００４】
そこで、解への収束を加速させる方法として、多重格子法（マルチグリッド法：ＭＧ法）が提案されている。この多重格子法は、格子のナイキスト周波数と収束性との関係を逆手に取ったものであり、定義域としての格子の幅を変化させることによって、求解工程で高速に減少する誤差の周波数成分を変化させるものである（非特許文献１、非特許文献２参照）。すなわち、多重格子法では、格子の幅を複数設定することによって、誤差が高速に減少する周波数成分を誤差の全ての周波数帯域に亘って設定することができる。
【０００５】
この多重格子法を用いた従来の典型的な計算処理例を、図５に基づいて説明する。
【０００６】
まず、ステップＳ５０１では、誤差の周波数帯域を分割する分割数（格子の数に相当）として、新たな分割数を探す。この場合、各軸方向の与えられた分割数以上の数値であって、所定の数（例えば「２」）のべき乗で表せる数値を新たな分割数とする。例示した「２」のべき乗で表すものとすると、１２８＝２７と表せるので、この「１２８」が各軸方向の与えられた分割数以上であれば、「１２８」を新たな分割数とする。また、２５６＝２８と表せるので、この「１２８」が各軸方向の与えられた分割数以上であれば、「１２８」を新たな分割数とする。このようなべき乗で表せる数値を新たな分割数とするために、係数行列Ａ、右辺ベクトルｂを補間する。
【０００７】
ステップＳ５０２では、係数行列Ａ、右辺ベクトルｂ、変数ベクトルｘに粗視化作用素（ｒｅｓｔｒｉｃｔｉｏｎ）を作用させることにより、誤差の低周波数成分を減衰させ易くする。ステップＳ５０３では、上記の粗視化作用素で粗視化された方程式を解く。
【０００８】
ステップＳ５０４では、係数行列Ａ、右辺ベクトルｂ、変数ベクトルｘに精細化作用素（ｐｒｏｌｏｎｇａｔｉｏｎ）を作用させることにより、誤差の高周波数成分を減衰させ易くする。ステップＳ５０５では、上記の精細化作用素で精細化された線型方程式を解く。なお、ステップＳ５０３，Ｓ５０５で用いる求解法としては、Ｇａｕｓｓ−Ｓｅｉｄｅｌ法やＳＯＲ法等の反復法を適当に選ぶことができる。ステップＳ５０６では、補間により、分割数をステップＳ５０１における元の分割数に戻す。
【非特許文献１】建部修見：“ＭＧＣＧ法：ロバストで高効率な並列解法，”東京大学博士論文，１９９６年．
【非特許文献２】Ｗ．Ｌ．  Ｂｒｉｇｇｓ，  Ｖ．Ｅ．  Ｈｅｎｓｏｎ  ：  “Ａ  ｍｕｌｔｉｇｒｉｄ  Ｔｕｔｏｒｉａｌ，”  ＳＩＡＭ，  ２０００．
【発明の開示】
【発明が解決しようとする課題】
【０００９】
しかし、上記従来例では、元の分割数は、所定の数値（上記例では「２」）のべき乗に限定されてしまう。この場合、「２」のべき乗で表せない分割数、例えば「１００」については、「２」のべき乗で表せる分割数「１２８」となるように格子を補間してＭＧ法で求解し、その後、分割数「１２８」が「１００」となるように格子を補間する。
【００１０】
ところが、このような補間処理では、誤差が混入する可能性がある。何故なら，分割数「１２８」の格子上では、分割数「１００」の格子上よりも高周波数成分を持った関数を表現できるからである。このような高周波数成分は、最終的な解を得るべく最後にダウンサンプリングする際に、折り返しノイズとして最終的な解に悪影響を与えることが考えられる。
【００１１】
本発明は、このような背景の下になされたもので、その目的は、多重格子法を用いて連立一次方程式を求解する場合に、格子の補間処理が不要な範囲でより多くの格子分割数を設定できるようにすることにある。
【課題を解決するための手段】
【００１２】
上記目的を達成するため、本発明は、多重格子法を用いて連立一次方程式を求解する数値計算装置であって、格子の分割数を素因数分解する分解手段と、前記分解に係る素数に対応する粗視化作用素を用いて連立一次方程式に係る線型方程式を粗視化する粗視化手段と、前記分解に係る素数に対応する精細化作用素を用いて連立一次方程式に係る線型方程式を精細化する精細化手段とを有することを特徴とする。
【発明の効果】
【００１３】
本発明によれば、多重格子法を用いて連立一次方程式を求解する場合に、格子の補間処理が不要な範囲でより多くの格子分割数を設定することができ、解の精度を向上させることが可能となる。
【発明を実施するための最良の形態】
【００１４】
以下、本発明を実施するための最良の形態を説明する。
【００１５】
図１は、本発明の実施の形態に係る数値計算装置の概略構成を示すブロック図である。この数値計算装置は、コンピュータにより構成され、ＣＰＵ１０１、表示装置１０２、入力装置１０３、一次記憶装置１０４、二次記憶装置１０５、通信装置１０６を有している。また、これらデバイスは、相互に通信可能にバスライン１０７を介して接続されている。
【００１６】
本実施の形態に係る後述する数値計算手法（解法アルゴリズム）は、アプリケーションプログラムとして記述され、このアプリケーションプログラムは、予め二次記憶装置１０５に格納されている。そして、二次記憶装置１０５に格納された上記アプリケーションプログラムは、入力装置１０３、或いは通信装置１０６等からのコマンド入力によって、一次記憶装置１０４にロードされた上で、ＣＰＵ１０１により実行される。このＣＰＵ１０１によるアプリケーションプログラムの実行状況、すなわち、本実施の形態に係る数値計算の計算過程は、表示装置１０２に表示される。
【００１７】
次に、本実施の形態に係る数値計算手法を、図２のフローチャートに基づいて説明する。
【００１８】
なお、ここでは、連立一次方程式を解くことにより、３次元空間での物理現象のシミュレーションを行なうことを想定している。
【００１９】
ＣＰＵ１０１は、まず、３次元空間上のＸ軸，Ｙ軸，Ｘ軸方向の格子に対して予め設定された複数組の中の１組の分割数ＮＸ，ＮＹ，ＮＺを、それぞれ設定値ＰＮＭＡＸ以下の素数で分解し、小さい素数の順に並べた素数ソーティングリストを作成する（ステップＳ２０１）。
【００２０】
例えば、ＰＮＭＡＸ＝５であり、Ｘ軸方向の格子の分割数が「１００」、Ｙ軸方向の格子の分割数が「１０８」、Ｘ軸方向の格子の分割数が「１２８」であると仮定する。この場合、１００＝２２５２、１０８＝２２３３、１２８＝２７といったように、素因数分解できる。従って、格子の分割数が「１００」であるＸ軸方向の素数ソーティングリストは、｛２，２，５，５｝となり、格子の分割数が「１０８」であるＹ軸方向の素数ソーティングリストは、｛２，２，３，３，３｝となる。また、格子の分割数が「１２８」であるＺ軸方向の素数ソーティングリストは、｛２，２，２，２，２，２，２｝となる。
【００２１】
なお、格子の分割数が「２９４」の場合は、この「２９４」を素因数分解をすると２９４＝２１３１７２となるが、素数「７」はＰＮＭＡＸである「５」を超えた値なので、「７」は素数ソーティングリストには含めず、｛２，３，４９｝といった素数ソーティングリストを作成する。すなわち、ＰＮＭＡＸ素数「７」は素数ソーティングリストにリストアップせずに、代わりに素数「７」を含む因数である「４９」をリストアップする。
【００２２】
また、分割数それ自体が素数であり、且つＰＮＭＡＸを超える場合、例えば分割数が「５３」の場合は、素数ソーティングリストは、空の｛｝とする。
【００２３】
ＣＰＵ１０１は、このような素数ソーティングリストを、Ｘ方向，Ｙ方向，Ｚ方向の各軸方向に対して、それぞれ１つずつ作成する。また、ＣＰＵ１０１は、各軸方向に対する３つの素数ソーティングリストの間でリストアップされた素数の個数が異なる場合は、素数の個数が最多数の素数ソーティングリストにおける素数の個数に統一する。この際、素数の個数が不足する素数ソーティングリストには、「１」を補充する。
【００２４】
例えば、Ｘ軸方向の素数ソーティングリストが｛２，２，３，３，３，５，５，５｝、Ｙ軸方向の素数ソーティングリストが｛３，５，７｝、Ｚ軸方向の素数ソーティングリストが｛３，３，３，５，５｝であったとすると、下記のように修正する。すなわち、分解に係る素数として最多数の８個の素数を持つＸ軸方向の素数ソーティングリストは修正せずに｛２，２，３，３，３，５，５，５｝のままとし、分解に係る素数の個数が不足するＹ軸方向の素数ソーティングリストは｛３，５，７，１，１，１，１，１｝と修正する。同様に分解に係る素数の個数が不足するＺ軸方向の素数ソーティングリストは｛３，３，３，５，５，１，１，１｝と修正する。
【００２５】
次に、ＣＰＵ１０１は、各素数ソーティングリストに含まれる素数に対応する粗視化作用素と精細化作用素をそれぞれ作成し、メモリに格納する（ステップＳ２０２）。なお、使用する可能性のある各種の素数に係る粗視化作用素と精細化作用素とを予め作成し、データベース化して二次記憶装置１０５等に保存しておくのが望ましい。
【００２６】
このように、使用する可能性のある各種の素数に係る粗視化作用素と精細化作用素とを予めデータベース化しておけば、数値計算の過程で粗視化作用素、精細化作用素を作成する頻度を激減させることができ、計算速度を向上させることが可能となる。
【００２７】
図３は、粗視化作用素の粗視化作用を示す概念図であり、図３（a）は素数「２」、図３（b）は素数「３」、図３（c）は素数「５」の粗視化作用を示している。
【００２８】
素数「２」,「３」,「５」の粗視化作用素を行列形式で記述すると、下記のようになる。
【００２９】
すなわち、素数「２」に対する粗視化作用素Ｒ２は、
【００３０】
【数１】
2007286738.tif 000003
【００３１】
素数「３」に対する粗視化作用素はＲ３、
【００３２】
【数２】
2007286738.tif 000004
【００３３】
素数「５」に対する粗視化作用素Ｒ５は、
【００３４】
【数３】
2007286738.tif 000005
【００３５】
となる。
【００３６】
図４は、精細化作用素の精細化作用を示す概念図であり、図４（a）は素数「２」、図４（b）は素数「３」、図４（c）は素数「５」の精細化作用を示している。
【００３７】
素数「２」,「３」,「５」の精細化作用素を行列形式で記述すると、下記のようになる。
【００３８】
すなわち、素数「２」に対する精細化作用素Ｐ２は、
【００３９】
【数４】
2007286738.tif 000006
【００４０】
素数「３」に対する精細化作用素Ｐ３は、
【００４１】
【数５】
2007286738.tif 000007
【００４２】
素数「５」に対する精細化作用素Ｐ５は、
【００４３】
【数６】
2007286738.tif 000008
【００４４】
となる。
【００４５】
次に、ＣＰＵ１０１は、Ｘ軸方向，Ｙ軸方向，及びＺ軸方向に係る各素数ソーティングリストの最初の方から順に、それぞれ素数を１個ずつ抽出し、その素数に対応する粗視化作用素をデータベースから取得する（ステップＳ２０３）。
【００４６】
そして、ＣＰＵ１０１は、取得した粗視化作用素を、連立一次方程式に係る線型方程式Ａｘ＝ｂにおける係数Ａ、変数ｂ,ｘに作用させて、当該線型方程式を粗視化することにより、誤差の低周波数成分を減衰させ易くする（ステップＳ２０４）。ここで、変数ｂ,ｘは３次元空間上に定義されたスカラー関数であり、係数Ａはその関数に対する作用素である。
【００４７】
ＣＰＵ１０１は、上記ステップＳ２０３，Ｓ２０４の処理を、各素数ソーティングリストの全ての素数（補充した１を含む）に対して順次実行する。なお、素数ソーティングリストが空の場合は、すなわち、ステップＳ２０１で予め与えられた分割数ＮＸ，ＮＹ，ＮＺが素数であり、且つＰＮＭＡＸを超える場合は、ＣＰＵ１０１は、粗視化処理、及び後述の精細化処理を実行しないようにする。
【００４８】
次に、ＣＰＵ１０１は、線型方程式Ａｘ＝ｂを解く（ステップＳ２０５）。この解法としては、Ｇａｕｓｓ−Ｓｅｉｄｅｌ法、ＳＯＲ、ＣＧ法等の通常の反復法を使用することができる。
【００４９】
次に、ＣＰＵ１０１は、Ｘ軸方向，Ｙ軸方向，及びＺ軸方向に係る各素数ソーティングリストの最後の方から順に、それぞれ素数を１個ずつ抽出し、その素数に対応する精細化作用素をデータベースから取得する（ステップＳ２０６）。
【００５０】
次に、ＣＰＵ１０１は、取得した精細化作用素を、連立一次方程式に係る線型方程式Ａｘ＝ｂにおける係数Ａ、変数ｂ,ｘに作用させて、当該線型方程式を精細化することにより、誤差の高周波数成分を減衰させ易くする（ステップＳ２０７）。そして、ＣＰＵ１０１は、精細化された線型方程式を解く（ステップＳ２０５）。
【００５１】
ＣＰＵ１０１は、上記のステップＳ２０６→Ｓ２０７→Ｓ２０５のシーケンスを、各素数ソーティングリストの全ての要素（素数）に対して順次実行する。また、ＣＰＵ１０１は、図２に係る一連の処理を、３次元空間上のＸ軸，Ｙ軸，Ｘ軸方向の格子に対して予め設定された複数組の分割数(すなわち、格子の幅)について、順次実行していく。
【００５２】
以上説明したように、本実施の形態では、設定された格子の分割数を素因数分解し、その分解に係る素数に基づいて、連立一次方程式に係る線型方程式の係数、変数に対して粗視化作用素、精細化作用素子を作用させている。
【００５３】
従って、変更していく格子の分割数について、従来のように分割数を所定数値のべき乗法で表現する場合よりも、本実施の形態では、選択の自由度が大幅に高くなる。換言すれば、格子の補間処理を行なう必要のない範囲でより多くの格子分割数を設定することができ、連立一次方程式をノイズの無い形で高精度に計算することが可能となる。
【００５４】
図２に示した処理は、マルチグリッド法の中でも最も単純な方法であり、Ｆｕｌｌ  ＭｕｌｔｉＧｒｉｄ  Ｖ−Ｃｙｃｌｅと呼ばれるものである（非特許文献２参照）。なお、Ｆｕｌｌ  ＭｕｌｔｉＧｒｉｄ  Ｖ−Ｃｙｃｌｅと同様に、Ｖ−Ｃｙｃｌｅ，Ｗ−ＣｙｃｌｅおよびＦＭＧ−Ｃｙｃｌｅを適用することも可能である。
【００５５】
また、これら全ての方法は、共役勾配法、共役残差法を含む非定常反復解法の前処理として利用することが出来る。また、本発明は、３次元空間以外の２次元空間等の多次元空間上の各次元に係る格子の分割数に対して適用することも可能である。
【００５６】
さらに、本発明の目的は、前述した実施の形態の機能を実現するソフトウェアのプログラムコードを記憶した記憶媒体を、システム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）が記憶媒体に格納されたプログラムコードを読み出し実行することによっても達成される。
【００５７】
この場合、記憶媒体から読み出されたプログラムコード自体が前述した実施の形態の機能を実現することになり、そのプログラムコード及び該プログラムコードを記憶した記憶媒体は本発明を構成することになる。
【００５８】
また、プログラムコードを供給するための記憶媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ−ＲＯＭ、ＤＶＤ−ＲＡＭ、ＤＶＤ−ＲＷ、ＤＶＤ＋ＲＷ等の光ディスク、磁気テープ、不揮発性のメモリカード、ＲＯＭ等を用いることができる。または、プログラムコードをネットワークを介してダウンロードしてもよい。
【００５９】
また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施の形態の機能が実現されるだけではなく、そのプログラムコードの指示に基づき、コンピュータ上で稼動しているＯＳ（オペレーティングシステム）等が実際の処理の一部または全部を行い、その処理によって前述した実施の形態の機能が実現される場合も含まれる。
【００６０】
さらに、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づき、その拡張機能を拡張ボードや拡張ユニットに備わるＣＰＵ等が実際の処理の一部または全部を行い、その処理によって前述した実施の形態の機能が実現される場合も含まれる。
【図面の簡単な説明】
【００６１】
【図１】本発明の実施の形態に係る数値計算装置の概略構成を示すブロック図である。
【図２】本実施の形態に係る数値計算手法を示すフローチャートである。
【図３】素数２，３，５に対する粗視化作用素の粗視化作用を示す概念図である。
【図４】素数２，３，５に対する精細化作用素の精細化作用を示す概念図である。
【図５】多重格子法を用いた従来の計算処理例を示すフローチャートである。
【符号の説明】
【００６２】
１０１…ＣＰＵ、１０３…入力装置、１０４…一次記憶装置、１０５…二次記憶装置
【図１】
2007286738.tif 000009
【図２】
2007286738.tif 000010
【図３】
2007286738.tif 000011
【図４】
2007286738.tif 000012
【図５】
2007286738.tif 000013