(57)【要約】【課題】ＳＯＲ法もしくはＳＳＯＲ法を用いた連立一次方程式の解法において、問題に応じて緩和パラメタの適切な値を高速かつ自動的に決定する技術を提供する。【解決手段】本発明は、数値計算装置における連立一次方程式の計算方法であって、前記数値計算装置が、解くべき連立一次方程式の係数行列に基づき緩和パラメタを決定する緩和パラメタ決定工程と、前記数値計算装置が、決定した緩和パラメタを用いてＳＯＲ法またはＳＳＯＲ法により前記連立一次方程式の解を求める求解工程と、を有し、前記緩和パラメタ決定工程では、前記係数行列と、前記係数行列に対する前処理行列との線形結合により定義される評価関数を用いて、前記緩和パラメタを決定する方法を提供する。【選択図】図１
2015049724.tif 000002
【特許請求の範囲】
【請求項１】
数値計算装置における連立一次方程式の計算方法であって、
前記数値計算装置が、解くべき連立一次方程式の係数行列に基づき緩和パラメタを決定する緩和パラメタ決定工程と、
前記数値計算装置が、決定した緩和パラメタを用いてＳＯＲ法またはＳＳＯＲ法により前記連立一次方程式の解を求める求解工程と、を有し、
前記緩和パラメタ決定工程では、前記係数行列と、前記係数行列に対する前処理行列との線形結合により定義される評価関数を用いて、前記緩和パラメタを決定する
ことを特徴とする連立一次方程式の計算方法。
【請求項２】
前記評価関数は、前記前処理行列の前記係数行列に対する近似の度合いを評価するための関数であり、
前記緩和パラメタ決定工程では、前記前処理行列の前記係数行列に対する近似の度合いが最も高くなるように、前記緩和パラメタを決定する
ことを特徴とする請求項１に記載の連立一次方程式の計算方法。
【請求項３】
前記係数行列の下三角要素行列をＬ、対角要素行列をＤ、上三角要素行列をＵ、前記緩和パラメタをωとしたときに、
前記前処理行列Ｍは、
【数１】
2015049724.tif 000009
である
ことを特徴とする請求項１または２に記載の連立一次方程式の計算方法。
【請求項４】
前記評価関数は、前記係数行列と前記前処理行列の差行列の列毎の和を成分とするベクトル、もしくは、前記差行列の行毎の和を成分とするベクトルのユークリッドノルムである
ことを特徴とする請求項１〜３のうちいずれか１項に記載の連立一次方程式の計算方法。
【請求項５】
前記数値計算装置が、マルチグリッド法もしくはマルチグリッド前処理の少なくとも一部の粗視化レベルのスムージング処理における近似解の更新に、前記緩和パラメタ決定工程および前記求解工程を適用する
ことを特徴とする請求項１〜４のうちいずれか１項に記載の連立一次方程式の計算方法。
【請求項６】
連立一次方程式を計算するための数値計算装置であって、
解くべき連立一次方程式の係数行列に基づき緩和パラメタを決定する緩和パラメタ決定手段と、
決定した緩和パラメタを用いてＳＯＲ法またはＳＳＯＲ法により前記連立一次方程式の解を求める求解手段と、を有し、
前記緩和パラメタ決定手段は、前記係数行列と、前記係数行列に対する前処理行列との
線形結合により定義される評価関数を用いて、前記緩和パラメタを決定する
ことを特徴とする数値計算装置。
【請求項７】
請求項１〜５のうちいずれか１項に記載の連立一次方程式の計算方法の各工程をコンピュータに実行させる
ことを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、連立一次方程式の計算方法及び数値計算装置に関し、特に、係数行列が大規模な疎行列である連立一次方程式の解法技術に関する。
【背景技術】
【０００２】
物理現象の解析など各種の工学的な問題は、模擬的に表現された偏微分方程式を離散化することにより、最終的に、大規模な疎行列を係数にもつ多元連立一次方程式に帰着する。従来より、この種の大規模係数行列をもつ連立一次方程式を計算機を用いて高速に且つ精度よく解くためのアルゴリズムが数多く提案されている（例えば特許文献１、非特許文献１参照）。
【０００３】
連立一次方程式の解法の一つである定常反復法は、マルチグリッド法の適用や、クリロフ部分空間法の前処理（内部ソルバとも呼ばれる）への応用が可能であり、有用な解法技術である。
定常反復法の一つであるＧＳ法（ガウス−ザイデル法）は、解くべき連立一次方程式をＡｘ＝ｂとした場合、次のいずれかの漸化式に従って近似解を更新する。
【数１】
2015049724.tif 000003
Ｌ、Ｄ、及びＵはそれぞれ、係数行列Ａの下三角要素行列、対角要素行列、及び上三角要素行列である。ｋは漸化式のステップ数を表す。
ＳＧＳ法（対称ガウス−ザイデル法）は式（１）と式（２）を交互に繰り返す手法である。共役勾配法等、係数行列に対称性を求めるタイプのクリロフ部分空間法の前処理に応用する場合に好適である。
【０００４】
ＳＯＲ法（逐次的過剰緩和法）およびＳＳＯＲ法（対称逐次的過剰緩和法）は、それぞれ、ＧＳ法およびＳＧＳ法に対して緩和パラメタωを導入し、収束性の向上を狙ったものである。式（１）、式（２）それぞれの処理の直後に以下の処理を実施する。
【数２】
2015049724.tif 000004
緩和パラメタωが問題に応じた適切な値に設定されていれば、元のＧＳ法およびＳＧＳ法より収束性が向上する。
【０００５】
マルチグリッド法は、問題を粗い解像度に粗視化していき、粗視化レベル毎に定常反復法（スムーサと呼ばれる）により誤差を推算し除去する解法である。各解像度に特徴的な誤差の周波数成分をまとめて除去することができるため収束性が高い。このスムーサにＳＯＲ法およびＳＳＯＲ法を用いる場合も、粗視化レベル毎の緩和パラメタωが適切な値であれば、ＧＳ法およびＳＧＳ法をスムーサに用いる場合より収束性が向上する。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１１−１４５９９９号公報
【非特許文献】
【０００７】
【非特許文献１】Yuichiro Miki and Teruyoshi Washizawa, "A2ILU: AUTO-ACCELERATED ILU PRECONDITIONER FOR SPARSE LINEAR SYSTEMS", SIAM J. SCI. COMPUT., Vol. 35, No. 2, pp. A1212-A1232, 2013
【発明の概要】
【発明が解決しようとする課題】
【０００８】
ＳＯＲ法およびＳＳＯＲ法では、利用者自身が緩和パラメタの値を決定しなければならない。緩和パラメタの値が適切でなければ、それぞれＧＳ法およびＳＧＳ法と比較して収束性が悪化し、解が得られなくなる可能性もある。また緩和パラメタの値が適切であっても、それを推算するのに多大な時間を要しては本末転倒である。
本発明はこのような背景の下になされたものであり、その目的は、ＳＯＲ法もしくはＳＳＯＲ法を用いた連立一次方程式の解法において、問題に応じて緩和パラメタの適切な値を高速かつ自動的に決定する技術を提供することにある。
【課題を解決するための手段】
【０００９】
本発明の第一態様は、数値計算装置における連立一次方程式の計算方法であって、前記数値計算装置が、解くべき連立一次方程式の係数行列に基づき緩和パラメタを決定する緩和パラメタ決定工程と、前記数値計算装置が、決定した緩和パラメタを用いてＳＯＲ法またはＳＳＯＲ法により前記連立一次方程式の解を求める求解工程と、を有し、前記緩和パラメタ決定工程では、前記係数行列と、前記係数行列に対する前処理行列との線形結合により定義される評価関数を用いて、前記緩和パラメタを決定することを特徴とする連立一次方程式の計算方法である。
【００１０】
本発明の第二態様は、連立一次方程式を計算するための数値計算装置であって、解くべき連立一次方程式の係数行列に基づき緩和パラメタを決定する緩和パラメタ決定手段と、決定した緩和パラメタを用いてＳＯＲ法またはＳＳＯＲ法により前記連立一次方程式の解を求める求解手段と、を有し、前記緩和パラメタ決定手段は、前記係数行列と、前記係数行列に対する前処理行列との線形結合により定義される評価関数を用いて、前記緩和パラメタを決定することを特徴とする数値計算装置である。
【００１１】
本発明の第三態様は、本発明の第一態様に係る連立一次方程式の計算方法の各工程をコンピュータに実行させることを特徴とするプログラムである。
【発明の効果】
【００１２】
本発明によれば、ＳＯＲ法もしくはＳＳＯＲ法を用いた連立一次方程式の解法において、問題に応じて緩和パラメタの適切な値を高速かつ自動的に決定することができる。
【図面の簡単な説明】
【００１３】
【図１】本発明の実施の形態に係る数値計算方法を示すフローチャート
【図２】本発明の実施の形態に係る数値計算装置の構成を示すブロック図
【図３】マルチグリッド法の計算処理例を示すフローチャート
【発明を実施するための形態】
【００１４】
以下に述べる本発明の実施形態は、数値計算装置における連立一次方程式の計算方法に関し、特に、ＳＯＲ法またはＳＳＯＲ法の改良に関するものである。具体的には、ＳＯＲ
法またはＳＳＯＲ法で用いる緩和パラメタ（加速パラメタとも呼ばれる）の値を、連立一次方程式の係数行列と該係数行列に対する前処理行列との線形結合により定義される評価関数を用いて決定する、という緩和パラメタ決定処理を導入する。この処理により、緩和パラメタの値が自動で決められるので、利用者によるパラメタの設定および入力の手間を省くことができる。また、解くべき連立一次方程式の係数行列に基づく評価関数を用いて緩和パラメタが決定されるので、解くべき問題に応じた適切な（最適または好適な）緩和パラメタの値を設定でき、従来の定常反復法よりも高速に求解することができる。しかも、評価関数が係数行列とその前処理行列との線形結合により定義される形式のものであるため、緩和パラメタの計算を高速にでき、推算に要する時間を実用的な時間に収めることができる。評価関数としては、前処理行列の係数行列に対する近似の度合いを評価するための関数を好ましく採用でき、その場合は近似の度合いが最も高くなるように緩和パラメタの決定を行えばよい。
本実施形態の方法は、大規模な疎行列を係数行列にもつ多元連立一次方程式の求解に好ましく適用でき、例えば、流体解析、電磁場解析、構造解析、挙動解析などの各種シミュレーションシステムに応用可能である。
以下、図面を参照して、本実施形態に係る数値計算装置の具体的な構成を例示的に説明する。
【００１５】
まず、本実施形態に係る数値計算装置として機能するコンピュータの基本構成を図２のブロック図に示す。該コンピュータは、大略、バス２０７を介して接続されたＣＰＵ２０１、表示装置２０２、入力装置２０３、ＲＡＭ２０４、ハードディスク装置２０５、ＮＩＣ２０６等から構成される。
【００１６】
ＣＰＵ（中央演算装置）２０１で、ＲＡＭ（主記憶装置）２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。表示装置２０２は、液晶画面などにより構成されており、ＣＰＵ２０１による処理結果を文字や画像などでもって表示することができる。入力装置２０３は、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。ＲＡＭ２０４は、ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。ハードディスク装置２０５には、ＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されている。これらのプログラムはＲＡＭ２０４にロードされて、ＣＰＵ２０１によって実行されることで、本コンピュータが連立一次方程式を数値的に解くための数値計算装置として機能する。ＮＩＣ（ネットワークインターフェース）２０６は、外部装置とのデータ通信を行う際のインターフェースとして機能するものであり、本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対象としても良い。また、クライアント−サーバシステム、グリッドコンピューティング、クラウドコンピューティングなどを利用してもよい。
【００１７】
次に、本実施形態に係る数値計算方法を、図１のフローチャートに基づいて説明する。図１のフローチャートに示す各工程は、数値計算装置（ＣＰＵ２０１）がプログラムに従って実行する処理を示す。なお、ここでは係数行列Ａ、解ベクトルｘ、右辺ベクトルｂとして、Ａｘ＝ｂで表される連立一次方程式の求解処理に、本実施形態の計算方法を適用する。
【００１８】
まず、数値計算装置が緩和パラメタωを決定する（緩和パラメタ決定ステップＳ１０１）。ここで、定常反復法のアルゴリズムは、係数行列Ａに対する前処理行列Ｍを定義する
ことにより得られる。ＳＳＯＲ法における前処理行列Ｍは係数行列Ａの下三角要素行列Ｌ、対角要素行列Ｄ、上三角要素行列Ｕ、及び緩和パラメタωを用いて以下の式で表わされる。
【数３】
2015049724.tif 000005
【００１９】
本実施形態では、この前処理行列Ｍが係数行列Ａを良く近似するほど定常反復法の収束性が向上することに着目する（以下、前処理行列Ｍを近似行列Ｍと呼ぶ。）。すなわち、近似行列Ｍの係数行列Ａに対する近似の度合いを表すなんらかの評価関数を定義し、近似の度合いが最も高くなるωを推算する。ただし、ωの推算に要する時間は、緩和パラメタの導入による収束性向上の効果を阻害しない程度に抑えなければ、実用的でない。そこで本実施の形態においては、評価関数に、係数行列Ａと近似行列Ｍとの差行列の列毎の和を成分とするベクトルのユークリッドノルムを用いる。すなわち評価関数ｆ（ω）は以下の式で表される。なお、ａ，ｌ，ｄ，ｕはそれぞれＡ，Ｌ，Ｄ，Ｕの成分を表し、ｎは行列の次数を表す。
【数４】
2015049724.tif 000006
この評価関数の場合、関数の値が小さいほど、行列ＡとＭの近似の度合いが高いといえる。したがって、評価関数の値を最小化する（最もゼロに近づける）ωが緩和パラメタの最適値である。
ただし上記式は本発明を限定するものではなく、評価関数は係数行列Ａと近似行列Ｍの線形結合により定義される何らかの関数であればよい。例えば係数行列Ａと近似行列Ｍとの差行列の行毎の和を成分とするベクトルのユークリッドノルムであってもよい。
【００２０】
ここで、上記の近似行列Ｍは係数行列Ａと比較し疎性が劣化している。そのためＡとＭの線形結合をメインメモリに陽に保存するには多大な時間を要する。そこで、評価関数を構築する際は、上記線形結合を要素毎にレジスタやキャッシュメモリ等の高速な記憶領域に展開し、逐一評価関数に反映させていく。構築された評価関数を最小化するωは、例えば評価関数の一次微分をゼロにするωをニュートン法等の手法により数値的に導くことにより推算できるが、これは本発明を限定するものではなく、多変数関数の最適化手法のいずれを用いても良い。なお、必ずしもωの最適値を求める必要はなく、評価関数を十分に小さくできる（例えば所定の基準値よりも小さくできる）ωを緩和パラメタに採用してもよい。
【００２１】
次に、数値計算装置が上記ステップＳ１０１により決定された緩和パラメタωを用いて、ＳＳＯＲ法により連立一次方程式を解く（求解ステップＳ１０２）。ＳＳＯＲ法は公知であるため詳しい説明は省略するが、式（１）と式（３）による近似解の更新と式（２）と式（３）による近似解の更新を、近似解の要求精度が満たされるまで交互に繰り返すものである。
【００２２】
本実施形態においては、本発明をＳＳＯＲ法に適用する場合について説明したが、本発明をＳＯＲ法に適用する場合についても同様である。
ただしその場合の近似行列Ｍは、昇順のＳＯＲ法の場合は
【数５】
2015049724.tif 000007
となり、降順のＳＯＲ法の場合は
【数６】
2015049724.tif 000008
となる。ＳＯＲ法も公知であるため詳しい説明は省略するが、昇順のＳＯＲ法は式（１）と式（３）による近似解の更新を、降順のＳＯＲ法は式（２）と式（３）による近似解の更新を、近似解の要求精度が満たされるまで交互に繰り返すものである。
【００２３】
これらの方法は、クリロフ部分空間法の前処理として利用することも可能である。特に共役勾配法等、係数行列に対称性を求めるタイプのクリロフ部分空間法の前処理に応用する場合は、本発明を適用したＳＳＯＲ法が好適である。
【００２４】
またこれらの方法は、マルチグリッド法もしくはマルチグリッド前処理のスムージング処理における近似解の更新に応用することも可能である。
マルチグリッド法の実施の一例について図３のフローチャートに基づいて説明する。図３のフローチャートに示す各工程は、数値計算装置（ＣＰＵ２０１）がプログラムに従って実行する処理を示す。
まず数値計算装置がプレスムーサにより近似解ｘを更新する（ステップＳ３０１）。次に、数値計算装置が係数行列Ａおよび残差（ｒ＝ｂ−Ａｘ）を粗い格子に粗視化する（ステップＳ３０２）。次に、数値計算装置が粗い格子で残差方程式（Ａｅ＝ｒ）を解き、当該格子に特徴的な周波数成分を含む誤差ｅを推算する（ステップＳ３０３）。次に、数値計算装置が得られた誤差ｅを精細化して元の格子に戻し、近似解ｘから除去する（ステップＳ３０４）。最後に数値計算装置がポストスムーサにより近似解ｘの更新を行う（ステップＳ３０５）。なお、ステップＳ３０３の線形方程式の求解処理に、ステップＳ３０１〜Ｓ３０５までを再帰的に用いれば、粗視化レベルを深くすることができる。
【００２５】
本発明の方法をマルチグリッド法もしくはマルチグリッド前処理に応用するためには、一部もしくは全ての粗視化レベルにおけるスムーサに、本発明を適用したＳＯＲ法もしくはＳＳＯＲ法を用いればよい。なお、共役勾配法等、係数行列に対称性を求めるタイプのクリロフ部分空間法のマルチグリッド前処理に応用する場合は、特に本発明を適用したＳＳＯＲ法が好適である。その場合、プレスムーサとポストスムーサの両方にＳＳＯＲ法を用いてもよい。あるいは、プレスムーサにＳＯＲ法の昇順処理（式（１）と式（３））を用い、ポストスムーサにＳＯＲ法の降順処理（式（２）と式（３））を用いてもよい。またその逆でもよい。これはオーダリングを変えた場合も同様である。すなわち、例えばＲｅｄ−Ｂｌａｃｋオーダリングを施した場合は、プレスムーサにＲｅｄ、Ｂｌａｃｋの順で処理を割り当て、ポストスムーサにはその逆の順序で処理を割り当てればよい。またその逆でもよい。
【００２６】
以上述べた本実施形態によれば、ＳＯＲ法またはＳＳＯＲ法で用いる緩和パラメタの値が自動で決められるので、利用者によるパラメタの設定および入力の手間を省くことができる。また、解くべき問題に応じた適切な（最適または好適な）緩和パラメタの値を設定でき、従来の定常反復法よりも高速に求解することができる。しかも、緩和パラメタの計
算を高速にでき、推算に要する時間を実用的な時間に収めることができる。したがって、本実施形態の計算方法（アルゴリズム）は、各種シミュレーションシステムにおける大規模な連立一次方程式の求解エンジンとして好ましく実装することができる。
【００２７】
記憶装置に記録されたプログラムを読み込み実行することで前述した実施形態の機能を実現するシステムや装置のコンピュータ（又はＣＰＵ、ＭＰＵ等のデバイス）によっても、本発明を実施することができる。また、例えば、記憶装置に記録されたプログラムを読み込み実行することで前述した実施形態の機能を実現するシステムや装置のコンピュータによって実行されるステップからなる方法によっても、本発明を実施することができる。この目的のために、上記プログラムは、例えば、ネットワークを通じて、又は、上記記憶装置となり得る様々なタイプの記録媒体（つまり、非一時的にデータを保持するコンピュータ読取可能な記録媒体）から、上記コンピュータに提供される。したがって、上記コンピュータ（ＣＰＵ、ＭＰＵ等のデバイスを含む）、上記方法、上記プログラム（プログラムコード、プログラムプロダクトを含む）、上記プログラムを非一時的に保持するコンピュータ読取可能な記録媒体は、いずれも本発明の範疇に含まれる。
【符号の説明】
【００２８】
２０１…ＣＰＵ，２０２…表示装置，２０３…入力装置，２０４…ＲＡＭ，２０５…ハードディスク装置，２０６…ＮＩＣ，２０７…バス
【図１】
2015049724.tif 000010
【図２】
2015049724.tif 000011
【図３】
2015049724.tif 000012