(57)【要約】【課題】連立方程式の係数行列の分解を効率化する。【解決手段】記憶部１１は、連立方程式に含まれる変数を分類したグループを示すグループ情報１３と、変数にそれぞれ割り当てられた変数番号に応じて連立方程式に含まれる係数が配置された行列１４とを記憶する。処理部１２は、グループ情報１３が示すグループ毎に、当該グループに属する変数の間で変数番号を入れ替えることで、行列１４から係数の配置を変更した行列１５を生成するオーダリング処理を実行する。処理部１２は、行列１５を行列１６と行列１７の積に分解する行列分解処理を実行する。処理部１２は、行列１６と行列１７とを用いて、変数の解１８を算出する。【選択図】図１
2019219811.tif 000002
【特許請求の範囲】
【請求項１】
連立方程式に含まれる複数の変数を分類したグループを示すグループ情報と、前記連立方程式に含まれる複数の係数を示す行列であって、前記複数の変数にそれぞれ割り当てられた変数番号に応じて前記複数の係数が配置された第１の行列とを記憶する記憶部と、
前記グループ情報が示すグループ毎に、前記複数の変数のうち当該グループに属する変数の間で変数番号を入れ替えることで、前記第１の行列から前記複数の係数の配置を変更した第２の行列を生成するオーダリング処理を実行し、前記第２の行列を第３の行列と第４の行列の積に分解する行列分解処理を実行し、前記第３の行列と前記第４の行列とを用いて、前記複数の変数の解を算出する処理部と、
を有する連立方程式処理装置。
【請求項２】
前記複数の変数は、物体の形状を示すモデルに含まれる辺に対応付けた第１のグループの変数と、前記モデルに含まれる節点に対応付けた第２のグループの変数とを含み、
前記第１の行列では、前記第１のグループの変数に対応する行が前記第２のグループの変数に対応する行よりも上側に配置されるように、前記第１のグループの変数および前記第２のグループの変数に対して変数番号が割り当てられており、
前記オーダリング処理では、前記処理部は、前記第１のグループの変数の変数番号と前記第２のグループの変数の変数番号との間の大小関係を維持する、
請求項１記載の連立方程式処理装置。
【請求項３】
前記第１のグループの変数は前記モデルに含まれる辺のベクトルポテンシャルを示し、前記第２のグループの変数は前記モデルに含まれる節点のスカラポテンシャルを示す、
請求項２記載の連立方程式処理装置。
【請求項４】
前記行列分解処理では、前記処理部は、前記第２の行列に含まれる行を入れ替えるピボッティング処理を抑制する、
請求項１記載の連立方程式処理装置。
【請求項５】
前記オーダリング処理では、前記処理部は、ある変数に割り当てられた第１の変数番号を第２の変数番号に変更する場合、前記第１の行列に含まれる前記第１の変数番号に対応する第１の行を前記第２の変数番号に対応する第２の行に移動させ、前記第１の変数番号に対応する第１の列を前記第２の変数番号に対応する第２の列に移動させる、
請求項１記載の連立方程式処理装置。
【請求項６】
前記オーダリング処理では、前記処理部は、前記グループ情報が示すグループ毎に、前記第１の行列に含まれる複数の行のうち当該グループに属する変数に対応する行の間で、非ゼロ要素が少ない行が上側に配置されるように変数番号を入れ替える、
請求項１記載の連立方程式処理装置。
【請求項７】
コンピュータが、
連立方程式に含まれる複数の変数を分類したグループを示すグループ情報と、前記連立方程式に含まれる複数の係数を示す行列であって、前記複数の変数にそれぞれ割り当てられた変数番号に応じて前記複数の係数が配置された第１の行列とを取得し、
前記グループ情報が示すグループ毎に、前記複数の変数のうち当該グループに属する変数の間で変数番号を入れ替えることで、前記第１の行列から前記複数の係数の配置を変更した第２の行列を生成するオーダリング処理を実行し、
前記第２の行列を第３の行列と第４の行列の積に分解する行列分解処理を実行し、
前記第３の行列と前記第４の行列とを用いて、前記複数の変数の解を算出する、
連立方程式処理方法。
【請求項８】
コンピュータに、
連立方程式に含まれる複数の変数を分類したグループを示すグループ情報と、前記連立方程式に含まれる複数の係数を示す行列であって、前記複数の変数にそれぞれ割り当てられた変数番号に応じて前記複数の係数が配置された第１の行列とを取得し、
前記グループ情報が示すグループ毎に、前記複数の変数のうち当該グループに属する変数の間で変数番号を入れ替えることで、前記第１の行列から前記複数の係数の配置を変更した第２の行列を生成するオーダリング処理を実行し、
前記第２の行列を第３の行列と第４の行列の積に分解する行列分解処理を実行し、
前記第３の行列と前記第４の行列とを用いて、前記複数の変数の解を算出する、
処理を実行させる連立方程式処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は連立方程式処理装置、連立方程式処理方法および連立方程式処理プログラムに関する。
【背景技術】
【０００２】
シミュレーションや科学技術計算などの計算分野では、多数の変数を含む大規模連立方程式の解をコンピュータに計算させることがある。その場合、連立方程式の係数を示す係数行列を生成し、係数行列に対して所定の行列演算を行うことで解を算出する。連立方程式の解を求める方法の１つとして、ＬＵ分解を用いた直接法が挙げられる。
【０００３】
ＬＵ分解を用いた直接法では、係数行列を、対角線より上側の要素がゼロである下三角行列（行列Ｌ）と、対角線より下側の要素がゼロである上三角行列（行列Ｕ）との積に分解する。そして、行列Ｌを左辺から右辺に移す行列演算に相当する前進消去と、行列Ｕを左辺から右辺に移す行列演算に相当する後退代入とを通じて、変数の解を求める。
【０００４】
元の係数行列が、値がゼロでない要素（非ゼロ要素）が少ない疎行列である場合、ＬＵ分解によって生成される行列Ｌや行列Ｕも非ゼロ要素が少ないことが好ましい。しかし、ＬＵ分解では、ある行番号および列番号によって特定される要素が、元の係数行列では値がゼロの要素（ゼロ要素）であっても、行列Ｌまたは行列Ｕでは非ゼロ要素になってしまい、非ゼロ要素が増加するというフィルインが発生することがある。フィルインの発生は、元の係数行列における非ゼロ要素の配置パターンに依存する。例えば、係数行列の上側に非ゼロ要素が多い場合、フィルインが発生しやすいことがある。
【０００５】
そこで、ＬＵ分解の前処理として、フィルインを低減するために変数の順序を入れ替えるオーダリングを行うことがある。例えば、オーダリングによって、非ゼロ要素が少ない行が上側に配置され非ゼロ要素が多い行が下側に配置されるように、変数の順序を入れ替えることで、ＬＵ分解後の行列Ｌや行列Ｕの非ゼロ要素を削減できることがある。オーダリング方法としては、ＭＤ（Minimum Degree）オーダリングやＡＭＤ（Approximate Minimum Degree）オーダリングなどが提案されている。
【０００６】
なお、複数の連立方程式を逐次的に処理する連立方程式求解装置が提案されている。提案の連立方程式求解装置は、前回の係数行列と今回の係数行列との間の差分を算出し、差分に基づいて、今回の係数行列の中でオーダリングを行うべき部分オーダリング領域を決定する。連立方程式求解装置は、決定した部分オーダリング領域以外の領域に対しては前回のオーダリング結果を流用して、連立方程式の解を求める。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００５−１８３６６号公報
【非特許文献】
【０００８】
【非特許文献１】William F. Tinney and John W. Walker, "Direct Solutions of Sparse Network Equations by Optimally Ordered Triangular Factorization", Proceedings of the IEEE, Volume 55 Issue 11, pp.1801-1809, November 1967
【非特許文献２】Patrick R. Amestoy, Timothy A. Davis and Iain S. Duff, "An Approximate Minimum Degree Ordering Algorithm", SIAM (Society for Industrial and Applied Mathematics) Journal on Matrix Analysis and Applications, Volume 17 Issue 4, pp.886-905, December 1996
【発明の概要】
【発明が解決しようとする課題】
【０００９】
ＬＵ分解は、例えば、元の係数行列の要素に基づいて、行列Ｕの上側の行から下側の行に向かって要素を確定していき、行列Ｕと並行して行列Ｌの左側の列から右側の列に向かって要素を確定していく。このとき、行列Ｌの要素の計算において、確定済みの行列Ｕの対角要素を分母として使用する除算が行われることがある。分母として使用する行列Ｕの対角要素がゼロである場合、そのままでは除算を行うことができず行列Ｌの要素が計算不能になってしまう。そこで、ゼロによる除算を避けるため、ＬＵ分解の途中で、次に処理すべき行とそれより下側にある未処理の行とを入れ替えるピボッティングが行われることがある。ピボッティングは、ピボット選択や軸選択と言うこともある。
【００１０】
ＬＵ分解の計算順序の関係上、元の係数行列の上側の対角要素に値がゼロ以下の要素が存在する場合、行列Ｕの対角要素がゼロになるリスクが大きくなる。そこで、ＬＵ分解の途中でピボッティングが行われることを抑制する方法として、ゼロ以下の対角要素ができる限り係数行列の下側に配置されるように変数の順序を決定し、決定した変数の順序に従って係数行列を作成しておくことが考えられる。
【００１１】
しかし、ＬＵ分解の前処理として前述のオーダリングを行うと、元の係数行列とは変数の順序が変わってしまい、オーダリング後の係数行列の上側にゼロ以下の対角要素が出現しやすくなってしまう。よって、係数の配置パターンを調整する上で、フィルインの抑制とピボッティングの抑制をどの様に両立させるかが問題となる。
【００１２】
１つの側面では、本発明は、係数行列の分解を効率化する連立方程式処理装置、連立方程式処理方法および連立方程式処理プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１３】
１つの態様では、記憶部と処理部とを有する連立方程式処理装置が提供される。記憶部は、連立方程式に含まれる複数の変数を分類したグループを示すグループ情報と、連立方程式に含まれる複数の係数を示す行列であって、複数の変数にそれぞれ割り当てられた変数番号に応じて複数の係数が配置された第１の行列とを記憶する。処理部は、グループ情報が示すグループ毎に、複数の変数のうち当該グループに属する変数の間で変数番号を入れ替えることで、第１の行列から複数の係数の配置を変更した第２の行列を生成するオーダリング処理を実行し、第２の行列を第３の行列と第４の行列の積に分解する行列分解処理を実行し、第３の行列と第４の行列とを用いて、複数の変数の解を算出する。
【００１４】
また、１つの態様では、連立方程式処理方法が提供される。また、１つの態様では、コンピュータに実行させる連立方程式処理プログラムが提供される。
【発明の効果】
【００１５】
１つの側面では、係数行列の分解を効率化できる。
【図面の簡単な説明】
【００１６】
【図１】第１の実施の形態の連立方程式処理装置の例を説明する図である。
【図２】シミュレーション装置のハードウェア例を示すブロック図である。
【図３】シミュレーション装置の機能例を示すブロック図である。
【図４】ＬＵ分解前の係数行列のオーダリング例を示す図である。
【図５】係数行列データの例を示す図である。
【図６】シミュレーション対象の電気回路の例を示す図である。
【図７】モデルの領域分割の例を示す図である。
【図８】モデルデータと物理データの例を示す図である。
【図９】領域内でのオーダリングの例を示す図である。
【図１０】パラメータデータと種類ＩＤベクトルの例を示す図である。
【図１１】結果データの例を示す図である。
【図１２】シミュレーションの手順例を示すフローチャートである。
【図１３】シミュレーションの手順例を示すフローチャート（続き）である。
【発明を実施するための形態】
【００１７】
以下、本実施の形態を図面を参照して説明する。
［第１の実施の形態］
第１の実施の形態を説明する。
【００１８】
図１は、第１の実施の形態の連立方程式処理装置の例を説明する図である。
第１の実施の形態の連立方程式処理装置１０は、ＬＵ分解を用いた直接法により連立方程式の解を算出する。連立方程式処理装置１０は、電磁界シミュレーションなどのシミュレーションや科学技術計算に使用することができる。連立方程式処理装置１０は、コンピュータや情報処理装置と言うこともでき、以下で説明する情報処理を並列に実行する並列処理装置であってもよい。また、連立方程式処理装置１０は、ユーザが操作するクライアント装置でもよいし、クライアント装置からアクセスされるサーバ装置でもよい。
【００１９】
連立方程式処理装置１０は、記憶部１１および処理部１２を有する。記憶部１１は、ＲＡＭ（Random Access Memory）などの揮発性の半導体メモリでもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性ストレージでもよい。処理部１２は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、ＤＳＰ（Digital Signal Processor）などのプロセッサである。ただし、処理部１２は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などの特定用途の電子回路を含んでもよい。複数のプロセッサの集合を「マルチプロセッサ」または単に「プロセッサ」と言うことがある。
【００２０】
プロセッサは、ＲＡＭなどのメモリ（記憶部１１でもよい）に記憶されたプログラムを実行する。プログラムには、後述する処理を規定した連立方程式処理プログラムが含まれる。連立方程式処理プログラムは、ユーザプログラムやアプリケーションプログラムから呼び出される数値計算ライブラリとして実装してもよい。
【００２１】
記憶部１１は、グループ情報１３および行列１４を記憶する。
グループ情報１３は、連立方程式に含まれる複数の変数を分類したグループを示す。変数は、連立方程式を解くことで値を算出すべき項目であり、未知数や未知変数などと言うこともできる。例えば、行列１４において同一のグループに属する変数が連続しており、グループ情報１３はグループ毎に当該グループに属する変数の範囲を示すものである。例えば、変数ｘ１，ｘ２がグループＧ１に属し、変数ｘ３，ｘ４がグループＧ２に属する。
【００２２】
変数のグループは、変数の性質に基づいて予め指定される。グループ情報１３は、ユーザが手作業で作成してもよいし、行列１４を生成するユーザプログラムやアプリケーションプログラムが自動的に作成してもよい。２以上の節点および２以上の辺を含むモデルを用いた有限要素法のシミュレーションの場合、例えば、節点に割り当てた２以上の変数から１つのグループが形成され、辺に割り当てた２以上の変数から別のグループが形成される。電磁界シミュレーションの場合、例えば、節点に割り当てた変数がスカラポテンシャルを示し、辺に割り当てた変数がベクトルポテンシャルを示す。
【００２３】
行列１４は、連立方程式に含まれる複数の係数を示す行列（第１の行列）である。行列１４では、複数の変数にそれぞれ割り当てられた変数番号に応じて複数の係数が配置されている。変数番号は、複数の変数の順番を示す一意な識別番号である。例えば、Ｎ個（Ｎは２以上の整数）の変数が定義されている場合、１からＮまでの整数がそれら変数に割り振られ、Ｎ行Ｎ列の行列１４が作成される。変数番号は行番号および列番号に対応する。
【００２４】
行列１４は、ユーザが手作業で作成してもよいし、ユーザプログラムやアプリケーションプログラムが自動的に作成してもよい。行列１４は、グループ情報１３が示すグループを考慮して作成される。例えば、同一のグループに属する変数に対応する行や列が連続するように行列１４が作成される。グループ間の変数番号の大小関係は、グループの性質に基づいて予め決定される。好ましくは、対角要素がゼロ以下になる行が、行列１４の上側では相対的に少なくなるようにし、行列１４の下側では相対的に多くなるようにする。
【００２５】
電磁界シミュレーションの場合、例えば、ベクトルポテンシャルを示す変数に小さい変数番号を割り当て、当該変数に対応する行が行列１４の上側に配置されるようにする。一方、スカラポテンシャルを示す変数に大きい変数番号を割り当て、当該変数に対応する行が行列１４の下側に配置されるようにする。これにより、行列１４において、上側の対角要素にゼロ以下の要素が少なくなり下側の対角要素にゼロ以下の要素が多くなる。
【００２６】
処理部１２は、グループ情報１３および行列１４を取得する。処理部１２は、行列１４に対してオーダリング処理を行って行列１５（第２の行列）を生成する。行列１５は、行列１４から複数の係数の配置を変更したものである。オーダリング処理は、後の行列分解処理においてフィルインが抑制されるようにすることを目的としている。
【００２７】
例えば、処理部１２は、オーダリング処理において、上側の行の非ゼロ要素が少なくなり下側の行の非ゼロ要素が多くなるように、変数の順序を入れ替える。変数の順序の入れ替えは、変数に割り当てた変数番号を入れ替えることで行う。ある変数に割り当てた変数番号＃１を変数番号＃２に変更する場合、例えば、処理部１２は、行列１４の行＃１を行＃２に移動させ、行列１４の列＃１を列＃２に移動させる。
【００２８】
このとき、処理部１２は、グループ情報１３が示すグループ毎に、複数の変数のうち当該グループに属する変数の間で変数番号を入れ替えることを許容し、異なるグループに属する変数の間で変数番号を入れ替えることを制限する。これにより、行列１４の作成時に決定されたグループ間の変数番号の大小関係が、オーダリング処理後の行列１５においても維持される。よって、行列１４において、あるグループに属する変数に対応する行が他のグループに属する変数に対応する行よりも上側に配置されている場合、行列１５においてもその上下関係が維持されて逆転してしまうことを抑制できる。対角要素がゼロ以下になる行が行列１４の上側で少ない場合、行列１５の上側でも少なくなる。
【００２９】
例えば、処理部１２は、変数ｘ１に対応する行＃１と変数ｘ２に対応する行＃２との入れ替えを許容し、変数ｘ３に対応する行＃３と変数ｘ４に対応する行＃４との入れ替えを許容する。一方、処理部１２は、変数ｘ１に対応する行＃１と変数ｘ３，ｘ４に対応する行＃３，＃４との入れ替えを禁止し、変数ｘ２に対応する行＃２と変数ｘ３，ｘ４に対応する行＃３，＃４との入れ替えを禁止する。
【００３０】
処理部１２は、行列１５を行列１６（第３の行列）と行列１７（第４の行列）の積に分解する行列分解処理を実行する。例えば、行列１６は対角線より上の要素がゼロである下三角行列であり、行列１７は対角線より下の要素がゼロである上三角行列である。この行列分解処理をＬＵ分解と言うことがある。行列分解処理には、外積形式ガウス法、内積形式ガウス法、クラウト法などの各種の分解方法を用いることができる。
【００３１】
行列分解処理の過程で、処理部１２は、行列１７の上側の行から下側の行に向かって要素を確定していき、行列１７と並行して行列１６の左側の列から右側の列に向かって要素を確定していくことが考えられる。行列１６の要素の計算にあたって、処理部１２は、確定した行列１７の対角要素を分母に使用した除算を行うことがある。そこで、処理部１２は行列分解処理の中で、ゼロによる除算を避けるため、次に処理する行とそれより下側にある未処理の行とを入れ替えるピボッティング処理を行うことがある。
【００３２】
ただし、上側の行にゼロ以下の対角要素が少なくなるように行列１５が作成されている場合、行列１７の対角要素がゼロになる可能性を低減できる。前述のように、上側の行にゼロ以下の対角要素が少なくなるように行列１４が作成されている場合、オーダリング処理を行ってもその性質が行列１５に引き継がれる。よって、処理部１２は、行列分解処理におけるピボッティング処理を抑制することができる。
【００３３】
また、行列分解処理によって生成される行列１６，１７の非ゼロ要素が行列１５よりも増加するというフィルインが発生することがある。ただし、上側の行に非ゼロ要素が少なく下側の行に非ゼロ要素が多くなるように行列１５が作成されている場合、フィルインを抑制できる。よって、上記のオーダリング処理によってフィルインが抑制される。
【００３４】
処理部１２は、行列１６，１７を用いて変数の解１８を算出する。処理部１２は、前進消去と後退代入の二段階の行列演算によって解１８を算出することができる。行列１６をＬ、行列１７をＵ、変数ベクトルをｘ、連立方程式の右辺ベクトルをｂとおくと、前進消去はＬＵｘ＝ｂをＵｘ＝Ｌ−１ｂに変形する行列演算に相当し、後退代入はＵｘ＝Ｌ−１ｂをｘ＝Ｕ−１Ｌ−１ｂに変形する行列演算に相当する。前進消去によって行列Ｌが右辺から左辺に移動し、後退代入によって行列Ｕが右辺から左辺に移動する。
【００３５】
処理部１２は、解１８を出力する。例えば、処理部１２は、グループ情報１３および行列１４を入力したユーザプログラムやアプリケーションプログラムに解１８を出力する。また、例えば、処理部１２は、連立方程式処理装置１０が有する記憶装置に解１８を保存する。また、例えば、処理部１２は、連立方程式処理装置１０が有するディスプレイに解１８を表示させるなど、出力デバイスに解１８を出力する。
【００３６】
第１の実施の形態の連立方程式処理装置１０によれば、直接法による連立方程式の求解において、行列分解処理の前にオーダリング処理が行われる。よって、フィルインが抑制され、行列分解処理によって生成される行列１６，１７の非ゼロ要素が少なくなる。このため、行列１６，１７を保持する記憶領域を削減でき、解１８を算出する処理時間を短縮することができる。また、オーダリング処理は同一のグループに属する変数の間で行われ、異なるグループに属する変数の間では制限される。よって、変数のグループを特定の順序に並べて行列１４を作成しておけば、行列分解処理においてゼロによる除算が発生する可能性を低減でき、ピボッティング処理を抑制することが可能となる。このように、オーダリングによるフィルインの抑制と変数のグループの順序付けによるピボッティングの抑制とを両立させることができ、行列分解処理を効率化することができる。
【００３７】
［第２の実施の形態］
次に、第２の実施の形態を説明する。
第２の実施の形態のシミュレーション装置１００は、有限要素法による電磁界シミュレーションを行う。電磁界シミュレーションにあたり、シミュレーション装置１００は、大規模連立一次方程式を生成し、その連立方程式の解を直接法によって算出する。シミュレーション装置１００は、コンピュータや情報処理装置と言うこともできる。また、後述するように、シミュレーション装置１００は複数のプロセッサを有しており、並列処理装置と言うこともできる。シミュレーション装置１００は、ユーザが操作するクライアント装置でもよいし、クライアント装置からアクセスされるサーバ装置でもよい。
【００３８】
図２は、シミュレーション装置のハードウェア例を示すブロック図である。
シミュレーション装置１００は、ＣＰＵ１０１，１０１−１，１０１−２，１０１−３を含む複数のＣＰＵ、ＲＡＭ１０２、ＨＤＤ１０３、画像信号処理部１０４、入力信号処理部１０５、媒体リーダ１０６および通信インタフェース１０７を有する。上記ユニットはバスに接続されている。シミュレーション装置１００は、第１の実施の形態の連立方程式処理装置１０に対応する。ＣＰＵ１０１は、第１の実施の形態の処理部１２に対応する。ＲＡＭ１０２またはＨＤＤ１０３は、第１の実施の形態の記憶部１１に対応する。
【００３９】
ＣＰＵ１０１，１０１−１，１０１−２，１０１−３は、プログラムの命令を実行するプロセッサである。ＣＰＵ１０１，１０１−１，１０１−２，１０１−３は並列にスレッドを実行することができる。ＣＰＵ１０１，１０１−１，１０１−２，１０１−３は、ＨＤＤ１０３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０２にロードし、プログラムを実行する。なお、シミュレーション装置１００は、複数のプロセッサコアを有してもよく、多数のプロセッサコアをもつＧＰＵを有してもよい。
【００４０】
ＲＡＭ１０２は、プロセッサが実行するプログラムやプロセッサが演算に使用するデータを一時的に記憶する揮発性の半導体メモリである。なお、シミュレーション装置１００は、ＲＡＭ以外の種類のメモリを備えてもよく、複数のメモリを備えてもよい。
【００４１】
ＨＤＤ１０３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、および、データを記憶する不揮発性ストレージである。なお、シミュレーション装置１００は、フラッシュメモリやＳＳＤ（Solid State Drive）など他の種類の記憶装置を備えてもよく、複数の記憶装置を備えてもよい。
【００４２】
画像信号処理部１０４は、ＣＰＵ１０１からの命令に従って、シミュレーション装置１００に接続されたディスプレイ１０４ａに画像を出力する。ディスプレイ１０４ａとしては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなど、任意の種類のディスプレイを使用することができる。
【００４３】
入力信号処理部１０５は、シミュレーション装置１００に接続された入力デバイス１０５ａから入力信号を受信する。入力デバイス１０５ａとして、マウス、タッチパネル、タッチパッド、キーボードなど、任意の種類の入力デバイスを使用できる。また、シミュレーション装置１００に複数の種類の入力デバイスが接続されてもよい。
【００４４】
媒体リーダ１０６は、記録媒体１０６ａに記録されたプログラムやデータを読み取る読み取り装置である。記録媒体１０６ａとして、例えば、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤなどの磁気ディスク、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）などの光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。媒体リーダ１０６は、例えば、記録媒体１０６ａから読み取ったプログラムやデータをＲＡＭ１０２またはＨＤＤ１０３に格納する。
【００４５】
通信インタフェース１０７は、ネットワーク１０７ａに接続され、ネットワーク１０７ａを介して他の情報処理装置と通信を行うインタフェースである。通信インタフェース１０７は、スイッチやルータなどの有線通信装置に接続される有線通信インタフェースでもよいし、基地局やアクセスポイントに接続される無線通信インタフェースでもよい。
【００４６】
図３は、シミュレーション装置の機能例を示すブロック図である。
シミュレーション装置１００は、モデル記憶部１１１、連立方程式生成部１１２、連立方程式記憶部１１３、線形ソルバ１１４、解記憶部１１５および結果表示部１１６を有する。モデル記憶部１１１、連立方程式記憶部１１３および解記憶部１１５は、例えば、ＲＡＭ１０２またはＨＤＤ１０３の記憶領域を用いて実現される。連立方程式生成部１１２および結果表示部１１６は、例えば、電磁界シミュレーションを行うアプリケーションプログラムを用いて実現される。線形ソルバ１１４は、例えば、アプリケーションプログラムから呼び出される数値計算用ライブラリプログラムを用いて実現される。
【００４７】
モデル記憶部１１１は、シミュレーション対象の物体の形状を示すモデルを記憶する。モデルは二次元または三次元の複数の要素に細分化されている。モデルは複数の節点（ノード）を含む。節点が要素の頂点になり、節点間を結合する辺（エッジ）が要素の外辺になる。また、モデル記憶部１１１は、シミュレーション対象の物体の材質に応じた物理特性を示す物理データを記憶する。例えば、モデルは電磁コイルの形状を示し、物理データは電磁コイルの電気伝導率や磁気抵抗率などの電磁的特性を示す。モデルおよび物理データは、シミュレーション装置１００で作成されてもよいし、他の情報処理装置で作成されてシミュレーション装置１００に入力されてもよい。
【００４８】
連立方程式生成部１１２は、モデルに対して複数の変数（未知数）を設定し、モデル記憶部１１１に記憶されたモデルおよび物理データと電磁気学の支配方程式とから、変数の解を求めるための連立方程式を生成する。例えば、連立方程式生成部１１２は、モデルの節点にスカラポテンシャルを示す変数を割り当て、モデルの辺にベクトルポテンシャルを示す変数を割り当てる。連立方程式を解くことで、各節点のスカラポテンシャルおよび各辺のベクトルポテンシャルを算出することができる。連立方程式は係数行列によって表現される。係数行列は多くの場合、次元数が多く非ゼロ要素が少ない大規模疎行列になる。
【００４９】
また、連立方程式生成部１１２は、連立方程式に含まれる変数の種類を示すパラメータデータを生成する。変数の種類は、モデルから係数行列を生成する際に連立方程式生成部１１２が判断することができる。ただし、シミュレーション装置１００を使用するユーザが変数の種類を指定するようにしてもよい。例えば、節点のスカラポテンシャルを示す変数と辺のベクトルポテンシャルを示す変数は、異なる種類の変数である。
【００５０】
係数行列は、同じ種類の変数に対応する行や列が連続するように生成される。係数行列を用いた連立方程式の求解の過程でも、異なる種類の変数に対応する行や列が混ざらないようにする。例えば、パラメータデータは、係数行列の１行目から１２行目および１列目から１２列目が第１の種類の変数に対応し、係数行列の１３行目から２１行目および１３列目から２１列目が第２の種類の変数に対応することを示す。例えば、上記の第１の種類の変数がスカラポテンシャルを示す変数であり、上記の第２の種類の変数がベクトルポテンシャルを示す変数である。連立方程式生成部１１２は、生成した係数行列およびパラメータデータを連立方程式記憶部１１３に格納する。
【００５１】
連立方程式記憶部１１３は、連立方程式生成部１１２によって生成された連立方程式の係数行列およびパラメータデータを記憶する。
線形ソルバ１１４は、連立方程式生成部１１２からの呼び出しに応じて、ＬＵ分解を用いた直接法により連立方程式の解を算出する。線形ソルバ１１４は、連立方程式記憶部１１３から係数行列およびパラメータデータを読み出す。線形ソルバ１１４は、係数行列に対してオーダリングを行い、係数行列に含まれる係数の配置を変更する。オーダリングには、ＭＤオーダリングやＡＭＤオーダリングなどの各種のオーダリング方法を用いることができる。ただし、線形ソルバ１１４は、同じ種類の変数の間では順序の入れ替えを許容する一方、異なる種類の変数の間では順序の入れ替えを禁止する。
【００５２】
線形ソルバ１１４は、オーダリング後の係数行列を、下三角行列と上三角行列の積に分解するＬＵ分解を行う。ＬＵ分解には、外積形式ガウス法、内積形式ガウス法、クラウト法などの各種のＬＵ分解方法を用いることができる。線形ソルバ１１４は、変数ベクトルを含む左辺からの下三角行列の消去に相当する前進消去と、左辺からの上三角行列の消去に相当する後退代入とを行って、変数の解を算出する。例えば、線形ソルバ１１４は、モデルに含まれる各節点のスカラポテンシャルと各辺のベクトルポテンシャルとを算出する。線形ソルバ１１４は、解を示す結果データを解記憶部１１５に格納する。
【００５３】
解記憶部１１５は、線形ソルバ１１４が出力した結果データを記憶する。
結果表示部１１６は、線形ソルバ１１４による求解が終了すると、解記憶部１１５に記憶された結果データをディスプレイ１０４ａに表示させる。ただし、結果表示部１１６は、シミュレーション装置１００が備える他の出力デバイスに結果データを出力してもよいし、他の情報処理装置に結果データを送信してもよい。このとき、結果表示部１１６は、変数の値を文字情報として出力してもよいし、図形要素に変換して可視化してもよい。例えば、結果表示部１１６は、変数の値を矢印などの図形要素に変換してモデルにマッピングし、当該モデルをディスプレイ１０４ａに表示させることが考えられる。電磁界シミュレーションの場合、これにより物体の電磁的性質を可視化することができる。
【００５４】
次に、連立方程式の解法について説明する。
連立方程式は、既知数である係数行列Ａおよび右辺ベクトルｂと未知数である変数ベクトルｘとを用いて、数式（１）のように定義される。係数行列Ａは、行数および列数がｎ（ｎは２以上の整数）の正方行列である。変数ベクトルｘおよび右辺ベクトルｂは、行数がｎの列ベクトルである。係数行列ＡはＬＵ分解によって、下三角行列である行列Ｌと上三角行列である行列Ｕの積に分解される。行列Ｌおよび行列Ｕは、行数および列数がｎの正方行列である。ＬＵ分解後の連立方程式からは、数式（２）に示すように、行列Ｌを左辺から右辺に移動させることに相当する前進消去と行列Ｕを左辺から右辺に移動させることに相当する後退代入とを通じて、変数ベクトルｘの解を算出できる。
【００５５】
【数１】
2019219811.tif 000003
【００５６】
【数２】
2019219811.tif 000004
【００５７】
数式（３）に示すように、行列Ｌは、対角要素が「１」であり、対角要素より上側の要素が全てゼロであり、対角要素より下側に非ゼロ要素が存在する正方行列である。行列Ｕは、対角要素より下側の要素が全てゼロであり、対角要素および対角要素より上側に非ゼロ要素が存在する正方行列である。行列Ｕの対角要素がゼロになる可能性もある。ここで、行列Ｌの対角要素が１であり、行列Ｌの対角要素より上側の要素がゼロであることから、係数行列Ａと、行列Ｌと行列Ｕの積との間の関係は数式（４）のようになる。そこで、行列Ｌの対角要素より下側の要素は数式（５）のように算出することができ、行列Ｕの対角要素より上側の要素は数式（６）のように算出することができる。
【００５８】
【数３】
2019219811.tif 000005
【００５９】
【数４】
2019219811.tif 000006
【００６０】
【数５】
2019219811.tif 000007
【００６１】
【数６】
2019219811.tif 000008
【００６２】
ＬＵ分解では、係数行列Ａを上側の行から下側の行に向かって走査する間に、行列Ｕの要素を上側の行から下側の行に向かって順に確定していき、行列Ｌの要素を左側の列から右側の列に向かって順に確定していく。このとき、数式（５）に示すように、行列Ｕの対角要素を分母として使用する除算を行うことがある。もし、分母として使用する行列Ｕの対角要素がゼロ要素である場合、そのままでは行列Ｌの要素を算出できない。そこで、ゼロによる除算が発生し得るときは、線形ソルバ１１４は、ＬＵ分解の過程でピボッティングを行う。ピボッティングでは、処理対象の行とそれより下側にある未処理の行とを入れ替えることで、行列Ｕの処理対象の行の対角要素がゼロになることを回避する。
【００６３】
ただし、分母として使用する行列Ｕの対角要素がゼロになることは、ＬＵ分解の前に係数行列Ａの作成方法を工夫することで抑制できる。係数行列Ａおよび行列Ｕは上側の行から順に処理していくという性質上、係数行列Ａの上側にゼロまたは負の対角要素が多く存在する場合、分母として使用する行列Ｕの対角要素がゼロになるリスクが大きくなる。一方、係数行列Ａの下側にゼロまたは負の対角要素が多く存在しても、行列Ｕを更新していく過程で行列Ｕの対角要素の値が大きくなっていくことが多いため、分母として使用する行列Ｕの対角要素がゼロになるリスクは小さい。そこで、係数行列Ａにおいてゼロ以下の対角要素はできる限り下側に集まるように予め変数の順序を決めておくことで、ＬＵ分解中のピボッティングを削減できＬＵ分解の負荷を低減できる。
【００６４】
ところで、係数行列Ａはゼロ要素が多い疎行列であり、行列Ｌおよび行列Ｕも疎行列であることが期待される。ただし、係数行列Ａにおける非ゼロ要素の配置パターンによっては、行列Ｌおよび行列Ｕそれぞれの非ゼロ要素が係数行列Ａよりも多くなることがある。係数行列Ａではゼロ要素であった位置に、行列Ｌまたは行列Ｕでは非ゼロ要素が配置されてしまうことを、フィルインと言うことがある。
【００６５】
フィルインは、非ゼロ要素が多い行が係数行列Ａの上側に集まり、非ゼロ要素が少ない行が係数行列Ａの下側に集まっている場合に発生しやすい。例えば、数式（７）の係数行列Ａは、対角線および第１列目に非ゼロ要素が集まっている。この係数行列Ａからは数式（８）のような行列Ｌが生成されることがある。数式（８）の行列Ｌは、第２列目や第３列目などでフィルインが発生して非ゼロ要素が増加している。一方、数式（９）の係数行列Ａは、対角線および第ｎ行目に非ゼロ要素が集まっている。この係数行列Ａからは数式（１０）のような行列Ｌが生成されることがある。数式（１０）の行列Ｌは、フィルインが発生しておらず非ゼロ要素が増加していない。
【００６６】
【数７】
2019219811.tif 000009
【００６７】
【数８】
2019219811.tif 000010
【００６８】
【数９】
2019219811.tif 000011
【００６９】
【数１０】
2019219811.tif 000012
【００７０】
フィルインが多く発生すると、行列Ｌや行列Ｕを保持する記憶領域が増大し、また、行列Ｌおよび行列Ｕを用いた前進消去および後退代入の負荷が増大する。そこで、線形ソルバ１１４は、ＬＵ分解の前処理としてオーダリングを行う。オーダリングでは、非ゼロ要素が少ない行ができる限り係数行列Ａの上側に存在し、非ゼロ要素が多い行ができる限り係数行列Ａの下側に存在するように、変数の順序を入れ替える。変数の順序の入れ替えによって、行の入れ替えおよび列の入れ替えが発生し、係数の配置が変わる。
【００７１】
図４は、ＬＵ分解前の係数行列のオーダリング例を示す図である。
係数行列Ａの一例である行列１３１は、６行６列の正方行列であり１３個の非ゼロ要素を含む。変数ｘ１は１行目および１列目に対応し、変数ｘ２は２行目および２列目に対応し、変数ｘ３は３行目および３列目に対応し、変数ｘ４は４行目および４列目に対応し、変数ｘ５は５行目および５列目に対応し、変数ｘ６は６行目および６列目に対応する。１行目の非ゼロ要素は２個、２行目の非ゼロ要素は２個、３行目の非ゼロ要素は３個、４行目の非ゼロ要素は３個、５行目の非ゼロ要素は１個、６行目の非ゼロ要素は３個である。
【００７２】
ベクトルＧは、変数ｘ１〜ｘ６のうち未処理の変数の変数番号を列挙したベクトルである。変数番号は行番号および列番号に対応する。初期状態では変数ｘ１〜ｘ６が順に並んでいるため、ベクトルＧの初期値は｛１，２，３，４，５，６｝である。ベクトルＨは、並び替え後の変数の変数番号を列挙したベクトルであり、その初期値は｛｝（空）である。
【００７３】
線形ソルバ１１４は、行列１３１から最も非ゼロ要素が少ない行を検索する。すると、行列１３１の５行目が検索される。線形ソルバ１１４は、ベクトルＧから変数番号＃５を削除し、ベクトルＨの末尾に変数番号＃５を挿入する。また、線形ソルバ１１４は、行列１３１から５行目と５列目を削除して、５行５列の行列１３２を生成する。行列１３２の１行目の非ゼロ要素は２個、２行目の非ゼロ要素は２個、３行目の非ゼロ要素は３個、４行目の非ゼロ要素は３個、５行目（元の６行目）の非ゼロ要素は３個である。
【００７４】
次に、線形ソルバ１１４は、行列１３２から最も非ゼロ要素が少ない行を検索する。最も非ゼロ要素が少ない行が複数ある場合、線形ソルバ１１４は、行番号の小さい行を優先的に選択することとする。すると、行列１３２の１行目が検索される。線形ソルバ１１４は、ベクトルＧから変数番号＃１を削除し、ベクトルＨの末尾に変数番号＃１を挿入する。また、線形ソルバ１１４は、行列１３２から１行目と１列目を削除して、４行４列の行列１３３を生成する。行列１３３の１行目（元の２行目）の非ゼロ要素は２個、２行目（元の３行目）の非ゼロ要素は２個、３行目（元の４行目）の非ゼロ要素は３個、４行目（元の６行目）の非ゼロ要素は３個である。
【００７５】
次に、線形ソルバ１１４は、行列１３３から最も非ゼロ要素が少ない１行目（元の２行目）を選択する。線形ソルバ１１４は、ベクトルＧから変数番号＃２を削除し、ベクトルＨの末尾に変数番号＃２を挿入する。また、線形ソルバ１１４は、行列１３３から１行目（元の２行目）と１列目（元の２列目）を削除して、３行３列の行列１３４を生成する。行列１３４の１行目（元の３行目）の非ゼロ要素は２個、２行目（元の４行目）の非ゼロ要素は２個、３行目（元の６行目）の非ゼロ要素は３個である。
【００７６】
次に、線形ソルバ１１４は、行列１３４から最も非ゼロ要素が少ない１行目（元の３行目）を選択する。線形ソルバ１１４は、ベクトルＧから変数番号＃３を削除し、ベクトルＨの末尾に変数番号＃３を挿入する。また、線形ソルバ１１４は、行列１３４から１行目（元の３行目）と１列目（元の３列目）を削除して、２行２列の行列１３５を生成する。
【００７７】
線形ソルバ１１４は、行列１３５から最も非ゼロ要素が少ない１行目（元の４行目）を選択する。線形ソルバ１１４は、ベクトルＧから変数番号＃４を削除し、ベクトルＨの末尾に変数番号＃４を挿入する。また、線形ソルバ１１４は、行列１３５から１行目（元の４行目）と１列目（元の４列目）を削除して、１行１列の行列１３６を生成する。最後に、ベクトルＧから変数番号＃６を削除し、ベクトルＨの末尾に変数番号＃６を挿入する。
【００７８】
その結果、ベクトルＨは｛５，１，２，３，４，６｝となる。これは、変数ｘ１〜ｘ６の順序をｘ５，ｘ１，ｘ２，ｘ３，ｘ４，ｘ６に変更することを示す。すなわち、変数ｘ５の変数番号を１に変換し、変数ｘ１の変数番号を２に変換し、変数ｘ２の変数番号を３に変換し、変数ｘ３の変数番号を４に変換し、変数ｘ４の変数番号を５に変換することを示している。よって、行列１３１の５行目が１行目に移動し、１行目が２行目に移動し、２行目が３行目に移動し、３行目が４行目に移動し、４行目が５行目に移動する。また、行列１３１の５列目が１列目に移動し、１列目が２列目に移動し、２列目が３列目に移動し、３列目が４列目に移動し、４列目が５列目に移動する。
【００７９】
これによって、係数行列Ａである行列１３１が行列１３７に変換される。行列１３７の１行目の非ゼロ要素は１個、２行目の非ゼロ要素は２個、３行目の非ゼロ要素は２個、４行目の非ゼロ要素は３個、５行目の非ゼロ要素は３個、６行目の非ゼロ要素は３個である。このようにして、係数行列Ａのオーダリングが行われる。ただし、上記のオーダリング方法は一例であり、他の方法でオーダリングを行うことも可能である。
【００８０】
以上、ピボッティングを抑制するための係数行列Ａの作成方法およびフィルインを抑制するための前処理であるオーダリングについて説明した。しかし、これら２つの係数配置方法は両立しないおそれがあるという問題がある。ゼロ以下の対角要素が上側に出現しないように予め係数行列Ａを作成しておいても、係数行列Ａの全体に対してオーダリングを行うと、ゼロ以下の対角要素が係数行列Ａの上側に出現してしまうことがある。その結果、ＬＵ分解の過程でピボッティングが発生してしまうことがある。
【００８１】
そこで、第２の実施の形態では、ピボッティングの抑制とフィルインの抑制のバランスを図るようにオーダリングを行うこととする。すなわち、ピボッティングが抑制されるように決定された変数の順序を尊重し、ピボッティングが十分に抑制される範囲で、フィルインができる限り少なくなるようにオーダリングを行う。
【００８２】
線形ソルバ１１４は、前述のパラメータデータに基づいて変数の種類を把握し、同じ種類の変数の間では変数番号の入れ替えによる順序変更を許容し、異なる種類の変数の間では変数番号の入れ替えによる順序変更を禁止する。これにより、係数行列Ａの下側に存在する対角要素がゼロ以下である行が、係数行列Ａの上側に移動してピボッティングが発生することを抑制できる。また、同じ種類の変数の間では、非ゼロ要素が少ない行が上方に移動するため、オーダリングを行わない場合よりもフィルインが減少する。
【００８３】
例えば、数式（１１）に示すような６行６列の係数行列Ａを考える。変数ｘ１と変数ｘ２が同一の種類であり、変数ｘ３と変数ｘ４が同一の種類であり、変数ｘ５と変数ｘ６が同一の種類である。この場合、線形ソルバ１１４は、変数ｘ１と変数ｘ２の入れ替えを許容し、変数ｘ３と変数ｘ４の入れ替えを許容し、変数ｘ５と変数ｘ６の入れ替えを許容する。一方、線形ソルバ１１４は、ピボッティングを抑制するために変数ｘ１，ｘ２のグループと変数ｘ３，ｘ４のグループと変数ｘ５，ｘ６のグループがその順に並べられていることを重視し、異なる種類の変数を入れ替えることを禁止する。
【００８４】
【数１１】
2019219811.tif 000013
【００８５】
なお、係数行列Ａは疎行列であることが多いため、ゼロ要素を省略する圧縮行格納法（ＣＲＳ：Compressed Row Storage）や圧縮列格納法（ＣＣＳ：Compressed Column Storage）などの圧縮形式で保存しておく。例えば、数式（１２）に示すような４行４列の係数行列Ａを考える。数式（１２）の係数行列Ａは、９個の非ゼロ要素を含む。４個の変数は３種類に分類されている。１番目の変数が１番目のグループに属し、２番目と３番目の変数が２番目のグループに属し、４番目の変数が３番目のグループに属する。
【００８６】
【数１２】
2019219811.tif 000014
【００８７】
図５は、係数行列データの例を示す図である。
連立方程式記憶部１１３は、圧縮行格納法で表現した係数行列データを記憶する。図５の係数行列データの例は、上記の数式（１２）の係数行列Ａに対応する。係数行列データは、値ベクトル１２３、列番号ベクトル１２４および行位置ベクトル１２５を含む。
【００８８】
値ベクトル１２３は、係数行列Ａに含まれる非ゼロ要素を列挙したベクトルである。値ベクトル１２３の長さは、係数行列Ａに含まれる非ゼロ要素の数と同じである。値ベクトル１２３は、行番号が小さい行の非ゼロ要素を先に列挙し、同じ行に複数の非ゼロ要素が存在する場合は列番号が小さい列の非ゼロ要素を先に列挙する。数式（１２）の例の場合、値ベクトル１２３は、１行１列、１行３列、１行４列、２行２列、３行１列、３行３列、３行４列、４行２列、４行４列の順に非ゼロ要素を列挙する。
【００８９】
列番号ベクトル１２４は、値ベクトル１２３に含まれる非ゼロ要素の列番号を列挙したベクトルである。列番号ベクトル１２４の長さは、係数行列Ａに含まれる非ゼロ要素の数と同じである。数式（１２）の例の場合、列番号ベクトル１２４は、列＃１、列＃３、列＃４、列＃２、列＃１、列＃３、列＃４、列＃２、列＃４の順に列番号を列挙する。
【００９０】
行位置ベクトル１２５は、値ベクトル１２３に含まれる非ゼロ要素のうち、各行の最初に出現する非ゼロ要素を指すポインタを列挙したベクトルである。ポインタは、値ベクトル１２３の要素番号によって表現される。行位置ベクトル１２５の長さは、係数行列Ａの行数より１だけ大きい。数式（１２）の例の場合、係数行列Ａの１行目は値ベクトル１２３の要素＃０から始まり、２行目は値ベクトル１２３の要素＃３から始まり、３行目は値ベクトル１２３の要素＃４から始まり、４行目は値ベクトル１２３の要素＃７から始まる。よって、行位置ベクトル１２５は、０，３，４，７の順に要素番号を列挙する。ただし、係数行列Ａの最後の非ゼロ要素を特定できるように、値ベクトル１２３の最大の要素番号より１だけおおきい番号を行位置ベクトル１２５の末尾に記載しておく。
【００９１】
次に、電磁界シミュレーションに使用する連立方程式の生成例を説明する。
図６は、シミュレーション対象の電気回路の例を示す図である。
シミュレーション対象の電気回路の例は、交流電源１４１、抵抗器１４２、コア１４３およびコイル１４４を含む。交流電源１４１および交流電源１４１に接続された抵抗器１４２は、コア１４３に巻かれたコイル１４４に電気を供給する入力回路である。交流電源１４１の入力電圧がＶ、交流電源１４１の近くの入力電流がＩ、抵抗器１４２の入力抵抗がＲである。入力電圧Ｖと入力抵抗Ｒは所与の既知数であり、入力電流Ｉは電磁界シミュレーションによって算出される未知数（変数）である。
【００９２】
電磁界シミュレーションは、コア１４３およびコア１４３の近くの空間を要素に細分化したモデルを使用した有限要素法によって行う。モデル上の各節点におけるスカラポテンシャルφと、モデル上の各辺におけるベクトルポテンシャルＰとが、電磁界シミュレーションによって算出される未知数（変数）である。よって、ベクトルポテンシャルＰとスカラポテンシャルφと入力電流Ｉが、連立方程式の変数になる。
【００９３】
ここで、モデルの節点および辺には変数番号が割り当てられる。また、複数のプロセッサを用いて並列処理できるように、モデルが複数の領域に分割される。
図７は、モデルの領域分割の例を示す図である。
【００９４】
モデル１５０は、節点１５１〜１５９を含む。
節点１５２，１５１の間の辺に変数ｘ１を割り当て、節点１５３，１５２の間の辺に変数ｘ２を割り当て、節点１５３，１５６の間の辺に変数ｘ３を割り当て、節点１５２，１５５の間の辺に変数ｘ４を割り当て、節点１５１，１５４の間の辺に変数ｘ５を割り当てている。また、節点１５５，１５４の間の辺に変数ｘ６を割り当て、節点１５６，１５５の間の辺に変数ｘ７を割り当て、節点１５４，１５７の間の辺に変数ｘ８を割り当て、節点１５５，１５８の間の辺に変数ｘ９を割り当て、節点１５６，１５９の間の辺に変数ｘ１０を割り当てている。また、節点１５８，１５７の間の辺に変数ｘ１１を割り当て、節点１５９，１５８の間の辺に変数ｘ１２を割り当てている。また、節点１５３に変数ｘ１３を割り当て、節点１５２に変数ｘ１４を割り当て、節点１５１に変数ｘ１５を割り当て、節点１５４に変数ｘ１６を割り当て、節点１５５に変数ｘ１７を割り当てている。また、節点１５６に変数ｘ１８を割り当て、節点１５７に変数ｘ１９を割り当て、節点１５８に変数ｘ２０を割り当て、節点１５９に変数ｘ２１を割り当てている。
【００９５】
変数ｘ１〜ｘ１２はベクトルポテンシャルを示す変数であり、変数ｘ１３〜ｘ２１はスカラポテンシャルを示す変数である。ベクトルポテンシャルを示す変数が１つのグループを形成し、スカラポテンシャルを示す変数が１つのグループを形成する。ベクトルポテンシャルを示す変数には、スカラポテンシャルを示す変数よりも小さい変数番号が割り当てられている。よって、変数ベクトルｘでは、ベクトルポテンシャルを示す変数が、スカラポテンシャルを示す変数よりも先に列挙されることになる。
【００９６】
このモデル１５０は、領域１５０−１と領域１５０−２に分割することが可能である。領域１５０−１と領域１５０−２は辺に沿って分割される。境界上の節点および辺は、領域１５０−１と領域１５０−２の両方に含まれる。領域１５０−１は節点１５１〜１５６の６個の節点を含み、領域１５０−２は節点１５４〜１５９の６個の節点を含む。
【００９７】
線形ソルバ１１４は、モデル１５０に対応する大きな係数行列から、領域１５０−１，１５０−２に対応する複数の小さな係数行列を生成する。線形ソルバ１１４は、小さな係数行列毎に独立に変数の解を直接法によって算出する。異なる小さな係数行列は、異なるプロセッサを用いて並列に処理することができる。線形ソルバ１１４は、複数の小さな係数行列に対応する変数の解を合成して、モデル１５０に対応する解ベクトルを生成する。１つの小さな係数行列にのみ属する変数については、解ベクトルの中の適切な位置にその解を配置すればよい。２以上の小さな係数行列に属する変数については、例えば、異なる小さな係数行列から算出された解を加算して解ベクトルの中の適切な位置に配置する。
【００９８】
図８は、モデルデータと物理データの例を示す図である。
モデルデータ１２１は、モデル記憶部１１１に記憶されている。モデルデータ１２１は、シミュレーション対象の物体の形状を表した二次元または三次元のモデルを示す。モデルデータ１２１は、節点を識別する節点番号と節点の座標とを対応付けた節点データを含む。また、モデルデータ１２１は、辺を識別する辺番号と辺の両端にある節点を列挙した節点リストとを対応付けた辺データを含む。また、モデルデータ１２１は、要素を識別する要素番号と要素を囲む辺を列挙した辺リストとを対応付けた要素データを含む。
【００９９】
物理データ１２２は、モデル記憶部１１１に記憶されている。物理データ１２２は、電磁界シミュレーションに使用する物理パラメータを示す。物理データ１２２が示す物理パラメータは、シミュレーション対象の物体の材質やその周辺の空間の材質に基づいて予め与えられる既知数である。例えば、物理データ１２２は、電気伝導率σ、磁気抵抗率ν、コイル電流Ｊ、入力電圧Ｖおよび入力抵抗Ｒを含む。
【０１００】
電磁界シミュレーションの係数行列Ａ、変数ベクトルｘおよび右辺ベクトルｂは、数式（１３）〜（１５）の支配方程式に基づいて決定することができる。数式（１３）は磁場の方程式であり、数式（１４）は渦電流の保存に関する電流保存の拘束条件であり、数式（１５）は磁束変化に伴う誘導起電力と連携する入力回路の方程式である。Ｐはベクトルポテンシャル、φはスカラポテンシャル、σは電気伝導率、νは磁気抵抗率、Ｊはコイルの電流、Φは磁束、Ｒは入力抵抗、Ｉは入力電流、Ｖは入力電圧、ｔは時間である。
【０１０１】
【数１３】
2019219811.tif 000015
【０１０２】
【数１４】
2019219811.tif 000016
【０１０３】
【数１５】
2019219811.tif 000017
【０１０４】
電磁界シミュレーションの連立方程式は、数式（１６）のように生成される。数式（１６）の１行目は数式（１３）に対応し、数式（１６）の２行目は数式（１４）に対応し、数式（１６）の３行目は数式（１５）に対応する。｛Ｐ｝はベクトルポテンシャルの変数を列挙した列ベクトルである。｛φ｝はスカラポテンシャルの変数を列挙した列ベクトルである。［Ｋ］は電気伝導率σ、磁気抵抗率νおよび辺の間の位置関係に応じた補間関数から算出される係数を列挙した小行列である。［Ｂ］は電気伝導率σおよび節点の間の位置関係に応じた補間関数から算出される係数を列挙した小行列である。｛Ｃ｝は入力回路とベクトルポテンシャルの関係を表す係数を列挙した行ベクトルである。｛Ｆ｝は電気伝導率σおよびコイル電流Ｊから算出される値を列挙した列ベクトルである。
【０１０５】
【数１６】
2019219811.tif 000018
【０１０６】
ベクトルポテンシャルを示す変数の個数（ベクトルポテンシャルの自由度）をｎｐ、スカラポテンシャルを示す変数の個数（スカラポテンシャルの自由度）をｎφとおく。すると、数式（１６）の係数行列Ａのサイズは（ｎｐ＋ｎφ＋１）×（ｎｐ＋ｎφ＋１）であり、変数ベクトルｘおよび右辺ベクトルｂのサイズは（ｎｐ＋ｎφ＋１）×１である。図７のモデル１５０の場合、ｎｐ＝１２、ｎφ＝９であるため、係数行列Ａのサイズは２１行２１列、変数ベクトルｘおよび右辺ベクトルｂのサイズは２１行１列となる。
【０１０７】
変数ベクトルｘは、ベクトルポテンシャルＰとスカラポテンシャルφと入力電流Ｉをその順に列挙している。これは、対角要素がゼロ以下である行を係数行列Ａの下側に集めるためである。もし、スカラポテンシャルφや入力電流ＩをベクトルポテンシャルＰよりも上側に配置してしまうと、対角要素がゼロ以下である行が係数行列Ａの上側に集まってしまう。よって、数式（１６）のように３種類の変数の順序を決定することで、ＬＵ分解におけるピボッティングを抑制することができる。
【０１０８】
線形ソルバ１１４は、オーダリングにおいて、ベクトルポテンシャルＰを示す変数の中で順序を入れ替えることは許容する。また、線形ソルバ１１４は、スカラポテンシャルφを示す変数の中で順序を入れ替えることは許容する。しかし、線形ソルバ１１４は、ベクトルポテンシャルＰを示す変数とスカラポテンシャルφを示す変数の間で順序を入れ替えることは禁止する。また、線形ソルバ１１４は、ベクトルポテンシャルＰを示す変数と入力電流Ｉを示す変数の間で順序を入れ替えることは禁止する。
【０１０９】
複数のプロセッサを用いた並列処理では、線形ソルバ１１４は、全体のモデルに対応する数式（１６）の連立方程式から、各領域に対応する数式（１７）の連立方程式を生成する。｛Ｐ（１）｝は特定の領域に含まれるベクトルポテンシャルの変数を列挙した列ベクトルである。｛φ（１）｝は特定の領域に含まれるスカラポテンシャルの変数を列挙した列ベクトルである。［Ｋ（１）］は［Ｋ］の中から特定の領域に関連する係数を抽出した小行列である。［Ｂ（１）］は［Ｂ］の中から特定の領域に関連する係数を抽出した小行列である。｛Ｆ（１）｝は｛Ｆ｝の中から特定の領域に関連する値を抽出した列ベクトルである。
【０１１０】
【数１７】
2019219811.tif 000019
【０１１１】
特定の領域に含まれるベクトルポテンシャルを示す変数の個数をｎｐ（１）、特定の領域に含まれるスカラポテンシャルを示す変数の個数をｎφ（１）とおく。数式（１７）の係数行列Ａのサイズは（ｎｐ（１）＋ｎφ（１）＋１）×（ｎｐ（１）＋ｎφ（１）＋１）であり、変数ベクトルｘおよび右辺ベクトルｂのサイズは（ｎｐ（１）＋ｎφ（１）＋１）×１である。図７の領域１５０−１の場合、ｎｐ（１）＝７、ｎφ（１）＝６であるため、係数行列Ａのサイズは１４行１４列、変数ベクトルｘおよび右辺ベクトルｂのサイズは１４行１列となる。
【０１１２】
領域分割を行った場合、オーダリングは各領域で行われる。
図９は、領域内でのオーダリングの例を示す図である。
ここでは、図７の領域１５０−１のオーダリングを考える。節点１５２，１５１の間の辺の変数ｘ１を変数ｘ２に変更し、節点１５３，１５２の間の辺の変数ｘ２を変数ｘ４に変更し、節点１５３，１５６の間の辺の変数ｘ３を変数ｘ５に変更している。また、節点１５２，１５５の間の辺の変数ｘ４を変数ｘ３に変更し、節点１５１，１５４の間の辺の変数ｘ５を変数ｘ１に変更している。また、節点１５３の変数ｘ１３を変数ｘ１５に変更し、節点１５１の変数ｘ１５を変数ｘ１３に変更している。
【０１１３】
領域１５０−２との境界に位置する辺の変数および節点の変数については、領域１５０−２の計算結果との整合性を維持するため、オーダリングの対象から除外している。また、辺に割り当てた変数同士は入れ替え可能とし、節点に割り当てた変数同士は入れ替え可能とする一方、辺に割り当てた変数と節点に割り当てた変数は入れ替え不可としている。よって、オーダリングを行っても、辺に割り当てた変数の変数番号が節点に割り当てた変数の変数番号よりも小さいという状況が維持される。
【０１１４】
線形ソルバ１１４は、変数の種類を認識するためにパラメータデータを取得する。
図１０は、パラメータデータと種類ＩＤベクトルの例を示す図である。
パラメータデータ１２６は、連立方程式記憶部１１３に係数行列データと併せて格納される。パラメータデータ１２６は、変数の種類の個数である種類数Ｎｐ、変数の総数Ｎａｌｌおよび各種類の変数の数Ｎ１，Ｎ２，…を含む。数式（１６）の連立方程式の場合、Ｎｐ＝３である。図７のモデル１５０および数式（１６）の連立方程式の場合、Ｎａｌｌ＝２２、Ｎ１＝１２、Ｎ２＝９、Ｎ３＝１である。
【０１１５】
種類ＩＤベクトル１２７は、領域分割およびオーダリングの前に、元の係数行列データおよびパラメータデータ１２６から線形ソルバ１１４が生成するベクトルである。種類ＩＤベクトル１２７は、係数行列データが示す各変数に割り当てた種類ＩＤを列挙する。種類ＩＤは、変数の種類を識別する識別番号である。１番目の種類の変数には種類ＩＤ＝１が割り当てられ、２番目の種類の変数には種類ＩＤ＝２が割り当てられ、３番目の種類の変数には種類ＩＤ＝３が割り当てられる。例えば、ベクトルポテンシャルＰを示す変数に種類ＩＤ＝１が割り当てられ、スカラポテンシャルφを示す変数に種類ＩＤ＝２が割り当てられ、入力電流Ｉを示す変数に種類ＩＤ＝３が割り当てられる。
【０１１６】
線形ソルバ１１４は、変数番号の小さい方から順に、パラメータデータ１２６が示す各種類の変数の個数だけ同じ種類ＩＤを割り当てていくことで、種類ＩＤベクトル１２７を生成することができる。例えば、線形ソルバ１１４は、変数番号の小さい方からＮ１個の変数に種類ＩＤ＝１を割り当て、次のＮ２個の変数に種類ＩＤ＝２を割り当て、次のＮ３個の変数に種類ＩＤ＝３を割り当てる。
【０１１７】
このようにして、ＬＵ分解が行われて解ベクトルが算出される。
図１１は、結果データの例を示す図である。
結果データ１２８は、解記憶部１１５に格納される。
【０１１８】
結果データ１２８は、モデルに含まれる辺を識別する辺番号と辺に対して算出されたベクトルポテンシャルＰとを対応付けたベクトルポテンシャルデータを含む。また、結果データ１２８は、モデルに含まれる節点を識別する節点番号と節点に対して算出されたスカラポテンシャルφとを対応付けたスカラポテンシャルデータを含む。また、結果データ１２８は、算出された入力電流Ｉを示す入力電流データを含む。
【０１１９】
なお、電磁界の時間変化のシミュレーションを行った場合、結果データ１２８は、異なる複数の時刻のベクトルポテンシャルデータ、異なる複数の時刻のスカラポテンシャルデータおよび異なる複数の時刻の入力電流データを含んでもよい。また、結果データ１２８に含まれる数値がディスプレイ１０４ａに表示されてもよいし、結果データ１２８に含まれる数値が矢印などの図形要素に変換されて可視化されてもよい。
【０１２０】
次に、シミュレーション装置１００の処理手順について説明する。
図１２は、シミュレーションの手順例を示すフローチャートである。
（Ｓ１０）連立方程式生成部１１２は、予め作成されたモデルデータ１２１および物理データ１２２をモデル記憶部１１１から読み込む。
【０１２１】
（Ｓ１１）連立方程式生成部１１２は、モデルデータ１２１が示す複数の辺にそれぞれベクトルポテンシャルＰを示す変数を割り当て、モデルデータ１２１が示す複数の節点にそれぞれスカラポテンシャルφを示す変数を割り当てる。また、連立方程式生成部１１２は、入力電流Ｉを示す変数を割り当てる。
【０１２２】
（Ｓ１２）連立方程式生成部１１２は、ステップＳ１１で割り当てた変数から変数ベクトルｘを生成し、変数ベクトルｘと物理データ１２２と所定の支配方程式から、連立方程式の係数行列Ａと右辺ベクトルｂを生成する。このとき、連立方程式生成部１１２は、ベクトルポテンシャルＰを示す変数、スカラポテンシャルφを示す変数、入力電流Ｉを示す変数の順に並ぶように変数ベクトルｘを生成する。また、連立方程式生成部１１２は、係数行列Ａを圧縮行格納法の係数行列データとして表現する。
【０１２３】
（Ｓ１３）連立方程式生成部１１２は、変数の個数に関するパラメータデータ１２６を生成する。種類数Ｎｐは、ベクトルポテンシャルＰの変数、スカラポテンシャルφの変数および入力電流Ｉの変数の３種類である。変数の数Ｎ１はモデルに含まれる辺の数、変数の数Ｎ２はモデルに含まれる節点の数、変数の数Ｎ３は１である。変数の総数Ｎａｌｌは、変数の数Ｎ１と変数の数Ｎ２と変数の数Ｎ３の合計である。
【０１２４】
（Ｓ１４）線形ソルバ１１４は、ｉ＝１に初期化する。
（Ｓ１５）線形ソルバ１１４は、パラメータデータ１２６から変数の数Ｎｉを読み込む。線形ソルバ１１４は、種類ＩＤをまだ割り当てていない変数のうち変数番号が小さい方からＮｉ個の変数に種類ＩＤ＝ｉを割り当てる。
【０１２５】
（Ｓ１６）線形ソルバ１１４は、ｉが種類数Ｎｐより小さいか判断する。ｉ＜Ｎｐの場合はステップＳ１７に進み、ｉ≧Ｎｐの場合はステップＳ１８に進む。
（Ｓ１７）線形ソルバ１１４は、ｉに１を加算する。そして、ステップＳ１５に進む。
【０１２６】
（Ｓ１８）線形ソルバ１１４は、モデル全体に対応する係数行列Ａに基づいて、モデルから分割された各領域に対応する係数行列Ａを生成する。モデルの分割方法は連立方程式生成部１１２が決定する。例えば、連立方程式生成部１１２から線形ソルバ１１４に、各領域に属する変数を通知するなど、係数行列Ａの分割方法が通知される。ただし、連立方程式生成部１１２から線形ソルバ１１４に、節点の隣接関係を示す隣接行列を通知し、線形ソルバ１１４が変数の集合を適切なサイズに分割してもよい。線形ソルバ１１４は、領域分割後の複数の係数行列Ａを複数のスレッドに割り振り、複数のプロセッサを用いて以下のステップＳ１９〜Ｓ２５を異なるスレッドに並列に実行させる。
【０１２７】
図１３は、シミュレーションの手順例を示すフローチャート（続き）である。
（Ｓ１９）線形ソルバ１１４は、ｉ＝１に初期化する。
（Ｓ２０）線形ソルバ１１４は、種類ＩＤ＝ｉの変数の間でオーダリングを実行する。例えば、線形ソルバ１１４は、種類ＩＤ＝ｉの変数に対応する行の中から最も非ゼロ要素が少ない行を選択し、選択した行を未処理の行よりも上側に移動させ、選択した行の行番号と同じ列番号をもつ列を未処理の列よりも左側に移動させる。線形ソルバ１１４は、上記を未処理の行がなくなるまで繰り返す。ただし、線形ソルバ１１４は、他のスレッドと共有する変数（領域の境界にある節点の変数や辺の変数）については移動対象から除外する。また、種類ＩＤ＝ｉの変数に対応する行の範囲を超えた行の移動や、種類ＩＤ＝ｉの変数に対応する列の範囲を超えた列の移動は行わない。
【０１２８】
（Ｓ２１）線形ソルバ１１４は、ｉが種類数Ｎｐより小さいか判断する。ｉ＜Ｎｐの場合はステップＳ２２に進み、ｉ≧Ｎｐの場合はステップＳ２３に進む。
（Ｓ２２）線形ソルバ１１４は、ｉに１を加算する。そして、ステップＳ２０に進む。
【０１２９】
（Ｓ２３）線形ソルバ１１４は、係数行列Ａを下三角行列である行列Ｌと上三角行列である行列Ｕの積に分解するＬＵ分解を実行する。このとき、線形ソルバ１１４は、係数行列Ａの上側の対角要素はゼロより大きいためピボッティングを抑制する。
【０１３０】
（Ｓ２４）線形ソルバ１１４は、行列Ｌと右辺ベクトルｂを用いて、連立方程式の左辺から行列Ｌを消去することに相当する前進消去を実行する。
（Ｓ２５）線形ソルバ１１４は、ステップＳ２４の前進消去の結果と行列Ｕを用いて、連立方程式の左辺から行列Ｕを消去することに相当する後退代入を実行する。これにより、変数ベクトルｘの解である解ベクトルが生成される。線形ソルバ１１４は、ステップＳ２０のオーダリングによって変更された変数の順序を元に戻す。すなわち、線形ソルバ１１４は、解ベクトルに含まれる値を、元の変数の順序に合わせて並べ替える。
【０１３１】
（Ｓ２６）線形ソルバ１１４は、複数のスレッドから解ベクトルを収集する。線形ソルバ１１４は、複数の解ベクトルを合成して全体の解ベクトルを生成する。このとき、線形ソルバ１１４は、１つのスレッドのみが値を計算した変数については、当該１つのスレッドが計算した値を解ベクトルの中の適切な位置（その変数に対応する行）に配置する。一方、線形ソルバ１１４は、２以上のスレッドが値を計算した変数（領域の境界にある辺の変数や節点の変数）については、それら２以上のスレッドが計算した値を合計して解ベクトルの中の適切な位置（その変数に対応する行）に配置する。
【０１３２】
（Ｓ２７）線形ソルバ１１４は、ベクトルポテンシャルＰの値とスカラポテンシャルφの値と入力電流Ｉの値を含む結果データ１２８を生成する。結果表示部１１６は、結果データ１２８をディスプレイ１０４ａに表示させる。例えば、結果表示部１１６は、結果データ１２８の数値をディスプレイ１０４ａに表示させる。また、例えば、結果表示部１１６は、結果データ１２８の数値を図形要素に変換してモデルにマッピングし、マッピング後のモデルを可視化データとしてディスプレイ１０４ａに表示させる。
【０１３３】
第２の実施の形態のシミュレーション装置１００によれば、対角要素がゼロ以下になる行が下側に集まるように係数行列が生成される。よって、係数行列のＬＵ分解の際のピボッティングが抑制され、ＬＵ分解の計算時間が短縮される。また、ＬＵ分解の前処理としてオーダリングが行われる。よって、ＬＵ分解によるフィルインが抑制され、ＬＵ分解後の行列を記憶する記憶領域を削減できると共に、前進消去および後退代入の計算時間が短縮される。また、オーダリングは、同じ種類の変数の間に限定して行われ、異なる種類の変数の間では禁止される。よって、オーダリングの副作用としてピボッティングが発生してしまうことを抑制でき、ピボッティングの抑制とフィルインの抑制のバランスを図ることができる。その結果、直接法による連立方程式の求解を効率化できる。
【符号の説明】
【０１３４】
１０  連立方程式処理装置
１１  記憶部
１２  処理部
１３  グループ情報
１４，１５，１６，１７  行列
１８  解
【図１】
2019219811.tif 000020
【図２】
2019219811.tif 000021
【図３】
2019219811.tif 000022
【図４】
2019219811.tif 000023
【図５】
2019219811.tif 000024
【図６】
2019219811.tif 000025
【図７】
2019219811.tif 000026
【図８】
2019219811.tif 000027
【図９】
2019219811.tif 000028
【図１０】
2019219811.tif 000029
【図１１】
2019219811.tif 000030
【図１２】
2019219811.tif 000031
【図１３】
2019219811.tif 000032