(57)【要約】  有限要素法と境界要素法の両方を用いて生成された連立方程式を行列形式であらわしたとき、その係数行列は、対角要素付近にノンゼロの値が集中する領域１と正方形に近い形の領域にノンゼロの値が集中する領域２が存在する。領域２に対して、対角要素の近傍の値だけを対象に不完全コレスキー分解を行うことで、前処理時間の短縮化を行う。有限要素法と境界要素法に基づく連立方程式の合成より成り立つ行列において、全ての未知数中の境界要素法で使用する未知数の並べ替えを初めに行い、続いて残りの未知数の並べ替えを行う。
2008026261.tif 000002
【特許請求の範囲】
【請求項１】
有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理方法であって、
境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得し、
該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替えることを特徴とする高速演算処理方法。
【請求項２】
前記解析対象のモデルが、１つの境界によって２つの領域Ａ及び領域Ｂに分割される場合、前記係数行列の中央部分に、該１つの境界に境界要素法を用いて得られた連立方程式の係数を配置し、その前後に、該２つの領域Ａ及び領域Ｂに有限要素法を適用して得られた連立方程式の係数を配置することを特徴とする請求項１に記載の高速演算処理方法。
【請求項３】
前記係数行列において、境界要素法による係数の前部分に配置される有限要素法による前記領域Ａについての係数は、領域Ａを記述する未知数であって、かつ、前記１つの境界を記述する未知数の内、最もインデックス値が小さい未知数に関する連立方程式の係数を開始位置として、カットヒル・マキー並べ替えを行うことを特徴とする請求項２に記載の高速演算処理方法。
【請求項４】
前記係数行列において、境界要素法による係数の後ろ部分に配置される有限要素法による前記領域Ｂについての係数は、領域Ｂを記述する未知数であって、かつ、前記１つの境界を記述する未知数の内、最もインデックス値が大きい未知数に関する連立方程式の係数を開始位置として、カットヒル・マキー並べ替えを行うことを特徴とする請求項２に記載の高速演算処理方法。
【請求項５】
前記係数行列において、境界要素法による係数は、対角要素を中心とするユーザが入力した値の幅に含まれるものについてのみ前処理が施されることを特徴とする請求項１に記載の高速演算処理方法。
【請求項６】
前記係数行列において、境界要素法による係数は、ユーザが入力した値より大きな係数値を有する係数についてのみ前処理が施されることを特徴とする請求項１に記載の高速演算処理方法。
【請求項７】
前記前処理は、不完全コレスキー分解であることを特徴とする請求項５または６に記載の高速演算処理方法。
【請求項８】
複数のＣＰＵを有する並列計算機を用いて、各ＣＰＵに前記係数行列の係数を均等に配分して計算されることを特徴とする請求項１に記載の高速演算処理方法。
【請求項９】
有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理方法をコンピュータに実現させるプログラムであって、
境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得し、
該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替えることを特徴とする高速演算処理方法をコンピュータに実現させるプログラム。
【請求項１０】
有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理装置であって、
境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得する辺取得手段と、
該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替える並べ替え手段と、
を備えることを特徴とする高速演算処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、有限要素法と境界要素法による結合方程式を高速に解いて、電子デバイス等の性能評価などを行うための高速演算処理方法に関する。
【背景技術】
【０００２】
有限要素法で得られた方程式と境界要素法で得られた方程式を結合した結合方程式を生成し、これを数値解析する場合（静電場、静磁場、動磁場、電磁場のシミュレーション）、大規模連立方程式を並列計算機環境で高速に解く必要がある。【０００３】
計算速度を向上するためには、計算機間の通信量を少なくするためのデータ並べ替えが必用である。また、境界要素法による連立方程式の行列は大きくなるため、前処理を行う領域を制限する必要があるが、領域の制限のしかたによって計算時間が大きく変わる。【０００４】
有限要素法と境界要素法を結合する数値解析分野では、従来、大規模連立方程式を直接的に解く手法（直説法）が使用されている。直接法は反復法に比べて計算負荷が大きく、膨大な計算時間を要する。【０００５】
有限要素法と境界要素法を結合する数値解析（静電場、静磁場、動磁場、電磁場）は、モデルの移動を考慮した解析が可能となるため、有効な解析手法である。通常、境界要素法で作成される連立方程式は０成分が無い密行列となり、解法には多くの時間を要することが分かっている。このような特徴を持つ連立方程式は反復法で解くことができるが、不完全コレスキー分解による前処理を行うことができないため収束効率が悪い。そのため、直接法が使用されている。また、有限要素法と結合しても、行列の密な領域が存在するため、同様に直接法による解法を用いるのが従来の解法である。【０００６】
図１は、行列の種類を説明する図である。
図１（ａ）は、疎行列である。疎行列とは、対角成分の近傍成分以外の係数が０となっている行列である。図１（ｂ）は、密行列である。密行列とは、０となる係数がない行列のことである。図１（ｃ）は、疎密行列である。疎密行列とは、行列のある部分は疎行列としての性質を持ち、他の部分は密行列としての性質を持つ、疎行列と密行列が混ざったような行列のことである。【０００７】
特許文献１においては、並列計算機の計算効率を上げるため、連立一次方程式の前処理行列を作成し、数値解の収束性を向上させる技術が開示されている。
【特許文献１】特開平５−８１３１０号公報
【発明の開示】
【０００８】
本発明の課題は、有限要素法と境界要素法とを結合する数値解析において演算すべき行列の演算を高速に行うことができる高速演算処理方法を提供することである。
本発明の高速演算処理方法は、有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理方法であって、境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得し、該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替えることを特徴とする。
【０００９】
本発明に基づく並べ替えを行うことで、1万自由度以上の未知数の一般的なモデルを用いた数値解析において、100倍以上高速に計算することが可能である。直接法では計算に必要なメモリ容量が多く、2GBで5万自由度が限界であるが、本発明を用いると100万自由度まで計算可能である。
【００１０】
従来の手法の手法では、直接法（計算時間増大、メモリ容量の増大化）を用いており、従来の不完全コレスキー分解の前処理付BiCG法反復法（並列化効率の低下）を使ったとしても、並べ替え処理無しでは並列化効率低下し、バンド幅制限をしなければ、計算時間増加、メモリ使用量が増加を招いてしまう。これに対し、本発明（並べ替え処理+バンド幅制限）では、並列計算の高速化を実現可能である。
【図面の簡単な説明】
【００１１】
【図１】行列の種類を説明する図である。
【図２】本発明の主要部の流れを示すフローチャートである。
【図３】本発明の実施形態の概略を説明する図である。
【図４】数値解析すべきモデルとそれによってできる行列について説明する図（その１）である。
【図５】数値解析すべきモデルとそれによってできる行列について説明する図（その２）である。
【図６】マルチCPU計算機での計算の仕方を説明する図（その１）である。
【図７】マルチCPU計算機での計算の仕方を説明する図（その２）である。
【図８】FEMで得られた連立方程式とBEMで得られた連立方程式の連立方法を説明するフローチャートである。
【図９】Cuthill-Mckeeの並べ替え手法を説明する図（その１）である。
【図１０】Cuthill-Mckeeの並べ替え手法を説明する図（その２）である。
【図１１】Cuthill-Mckeeの並べ替え手法を説明する図（その３）である。
【図１２】Cuthill-Mckeeの並べ替え手法を説明する図（その４）である。
【図１３】BEM領域に対する並べ替え処理を説明する図（その１）である。
【図１４】BEM領域に対する並べ替え処理を説明する図（その２）である。
【図１５】BEM領域に対する並べ替え処理を説明する図（その３）である。
【図１６】BEM領域に対する並べ替え処理を説明する図（その４）である。
【図１７】BEM領域に対する並べ替え処理を説明する図（その５）である。
【図１８】前処理を説明するための図である。
【図１９】マルチCPU型の並列計算機において本発明の実施形態の計算を行う場合の各CPUへのデータの振り分け方法を説明するフローチャート（その１）である。
【図２０】マルチCPU型の並列計算機において本発明の実施形態の計算を行う場合の各CPUへのデータの振り分け方法を説明するフローチャート（その２）である。
【発明を実施するための最良の形態】
【００１２】
本発明の実施形態では、境界要素法で扱う面の幾何形状より、未知数間の係数が大きいもの同士が行列内で近い位置にくるような並べ替えを行うことで、反復による収束を速くするための前処理の効率化を図る。【００１３】
図２は、本発明の主要部の流れを示すフローチャートである。
ここでは、有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理方法における並べ替え処理方法を示す。【００１４】
まず、Ｓｔｅｐ１で、境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得する。
Ｓｔｅｐ２では、前記解析対象のモデルが、１つの境界によって２つの領域Ａ及び領域Ｂに分割される場合、前記係数行列の中央部分に、該１つの境界に境界要素法を用いて得られた連立方程式の係数を配置し、その前後に、該２つの領域Ａ及び領域Ｂに有限要素法を適用して得られた連立方程式の係数を配置する。【００１５】
Ｓｔｅｐ３では、前記係数行列において、境界要素法による係数の前部分に配置される有限要素法による前記領域Ａについての係数は、領域Ａを記述する未知数であって、かつ、前記１つの境界を記述する未知数の内、最もインデックス値が小さい未知数に関する連立方程式の係数を開始位置として、カットヒル・マキー（Cuthill-Mckee）並べ替えを行う。また、前記係数行列において、境界要素法による係数の後ろ部分に配置される有限要素法による前記領域Ｂについての係数は、領域Ｂを記述する未知数であって、かつ、前記１つの境界を記述する未知数の内、最もインデックス値が大きい未知数に関する連立方程式の係数を開始位置として、カットヒル・マキー並べ替えを行う。【００１６】
Ｓｔｅｐ４では、境界要素法の係数について、該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替える。【００１７】
図３は、本発明の実施形態の概略を説明する図である。
境界要素法で解く領域は密領域であるため、反復解法を用いて前処理に不完全コレスキー分解を行うと、結果的に完全なコレスキー分解となってしまい、直接法による解法と同じになってしまう。そこで、密領域の全ての係数を使用せずに、前処理による反復法の収束性を高める必要がある。【００１８】
有限要素法と境界要素法の混合解法では、図３左図のように行列の対角要素付近にノンゼロの値が集中する領域１と正方形に近い形の領域にノンゼロの値が集中する領域２が存在する。有限要素法解析では一般に行列全体に対して不完全コレスキー分解を行うが、このような行列に対して不完全コレスキー分解を行うと、領域２において計算時間の増加を招いてしまい非効率である。そこで領域２に対して図３右図のように対角要素の近傍の値だけを対象に不完全コレスキー分解を行うことで、前処理時間の短縮化を行う。【００１９】
すなわち、本発明の実施形態においては、有限要素法と境界要素法に基づく連立方程式の合成より成り立つ行列において、全ての未知数中の境界要素法で使用する未知数の並べ替えを初めに行い、続いて残りの未知数の並べ替えを行う。【００２０】
図４及び図５は、数値解析すべきモデルとそれによってできる行列について説明する図である。
図４（ａ）のモデルは、ＨＤＤヘッドとＨＤＤ媒体との間の動作状態を数値解析するためのモデルである。このＨＤＤヘッド１１と媒体１０を含む空間をそれぞれモデルＢとモデルＡとして、格子に区切り、各格子点について物理量を考える。モデルＡとモデルＢの相互作用は、格子点間の相互作用として記述し、相互作用そのものは物理理論（たとえば、電磁気学）の法則にしたがって記述する。媒体１０を含むモデルＡ及び、ＨＤＤヘッド１１を含むモデルＢはそれぞれ、媒体１０、ＨＤＤヘッド１１のみではなく、その周りの空気領域も含めてモデル化したものである。このモデルＡとＢの内部は、有限要素法を用いて計算されるが、モデルＡとモデルＢの両方の最外周領域は、２つのモデルの最外周を１つの外周を表すモデルとしてモデル化し、外周を格子に区切って、物理量を計算する。図４（ａ）の場合、媒体１０とＨＤＤヘッド１１が主に相互作用する境界要素面１４が、モデルＡとモデルＢの相対する外周面として記述される。このような外周境界の計算は、境界要素法を用いて計算される。【００２１】
以上のようにモデル化された解析対象について、有限要素法と境界要素法を用いて連立方程式を立て、その係数を行列として整理すると、図４（ｂ）のような行列となる。すなわち、モデルＡの内部を記述する係数は、図４（ｂ）の行列の領域１のように対角要素の周辺にノンゼロの係数が集中する。同様に、モデルＢの内部を記述する係数は、領域３のように対角要素の周辺にノンゼロの係数が集中する。一方、境界をモデル化し、境界要素法で連立方程式をつくり、係数を行列表示すると、領域２のように、その部分だけ、密行列のようになっている。【００２２】
境界要素法で使用する未知数の並べ替えは、代数的にモデル端部に位置する未知数（図４（ａ）のａ点のような点における、求めるべき物理量をあらわすｘなどの未知数。モデルが格子化されているので、各格子点での物理量をｘｉのようにインデックスｉをつけて記述する）を開始No=１とし、それ以外の未知数との行列成分（連立方程式は、ａｉｊｘｊ＋ａｉｊ＋１ｘｊ＋１＋・・・＝ｂｉのような式からなっており、ｘｋをNo=1とした場合、ａｉｋ）が最も大きい未知数をNo=2として登録する。最後に登録された未知数No=nに対して、まだ登録されていない未知数との行列成分が最も大きい未知数を最後のNoにn+1番目で登録する。以上の手順でナンバーリングされた番号を境界要素法領域における連立方程式の未知数の順番と定義する。【００２３】
図４（ａ）のように、境界要素法を使用する場合、モデルの領域が複数の領域に分割されている。媒体移動を伴うHDDヘッドのライト・リードシミュレーションのように、ヘッドと媒体というように２つの大きな領域に分割されている場合に、以下の並べ替えを行うことで効率的な領域分割を実現する。【００２４】
連立方程式より作成される行列の領域１はモデル左側の媒体+空気領域の境界要素面以外の未知数から構成し、領域２は境界要素面の未知数から構成し、領域３はモデル右側のリードヘッド+空気領域の境界要素面以外の未知数から構成する。【００２５】
領域２の並べ替え後に、媒体側未知数でかつ領域２の未知数のインデックス番号が最も小さい番号を開始番号として領域１の並べ替えを行う。同様に、HDDヘッド側未知数でかつ領域２の未知数のインデックス番号が最も大きい番号を開始番号として領域３の並べ替えを行う。領域１〜３の並べ替え後の順番に従って、最後に領域１→領域２→領域３の通し番号でインデックス番号付けを行って並べ替えを終了する。【００２６】
境界要素法の未知数間の係数（未知数ｘｉと未知数ｘｊ間の係数はａｉｊ）は、３次元解析の場合は未知数が定義される面の単位法線ベクトルと、問題とする未知数が定義される位置から、相互作用を考慮すべき未知数が定義される位置への位置ベクトルとの内積に比例するため、同一平面に存在する面内の未知数間の係数は０となる。そこで、連立方程式作成過程における計算では、境界要素面を構成する全ての未知数の組み合わせに対して成分計算は行わず、予め同一平面上に存在するという情報を面要素グループに設定し、異なる面要素グループとの組み合わせに含まれる面に属する未知数間の係数だけを成分計算することで計算時間の短縮化を行う。【００２７】
なお、図４（ｂ）に示されるように、対角要素近傍のノンゼロ領域の幅をバンド幅と呼ぶ。領域２のバンド幅は広いため、前処理において、そのままのバンド幅を使用すると計算時間が増加してしまう。前処理に使用するバンド幅は、必ずしも各領域のバンド幅と同じである必要はない。そこで、図５（ａ）の領域２に対してはユーザーが任意にバンド幅を指定ユーザ処理を行う領域を図５（ｂ）の領域２のように制限することを可能とする。【００２８】
図６及び図７は、マルチCPU計算機での計算の仕方を説明する図である。
複数のCPUを使用する並列計算では、ユーザが使用する計算機ノード数を入力し、行列を図６のようにブロック分割する。また、各計算機内に複数のCPUが存在する場合、ユーザが並列スレッド数を入力すると、計算機ノード数×スレッド数で行列をブロック分割する。【００２９】
要素数を単純に分割するだけでは、各CPUに加わる計算負荷が均等にならない。計算負荷は、モデル形状やメッシュ（格子）分割方法、FEM（Finite Element Method：有限要素法）領域であるかBEM(Boundary Element Method：境界要素法）領域であるかに依存する。そこで、各未知数に対する一回の反復計算あたりの計算時間をテスト計算により予め算出し、計算付加が均等になるように各CPUに未知数を割り当てる。【００３０】
図６では、CPU1〜5に行列の係数を上から順に割り当てている。図６のCPU3に割り当てられた係数は、ユーザが制限したバンド幅がCPU3に割り当てられた係数のインデックス幅より大きい場合の前処理サイズである。一方、図７（ａ）のCPU3に割り当てられた係数は、ユーザが制限したバンド幅がCPU3に割り当てられた係数のインデックス幅より小さい場合の前処理サイズである。【００３１】
ユーザが設定したBEM領域のバンド幅制限サイズとCPUに割り当てられた未知数サイズの大小関係を自動的に判断し、図７（ｂ）のようなBEM領域の前処理ブロックを、図７（ｃ）のように、バンド幅制限サイズに自動的に制限して計算を行う。【００３２】
有限要素法で生成される行列のノンゼロ成分は、同じ材料物性グループ内では、大きい値と小さい値との差は小さい。しかしながら、境界要素法によって生成される行列成分は、未知数と未知数の間の係数の値は、未知数が配置されている位置間の距離の２乗に反比例する(３次元解析)ため、係数の大きさは幅広く分布している。そこで未知数間の係数が大きい値が行列対角部分に集中するような並べ替えを行い、対角近傍成分だけを用いて反復解法における不完全コレスキー分解を行う。どの範囲までの対角近傍成分幅(制限バンド幅)を使用するかは、ユーザが予め数値入力して指定する。また、最大係数に対してどの値までを前処理に含めるかという閾値をユーザが指定して、制限バンド幅を計算することを可能にする。【００３３】
HDDのヘッド-媒体解析やモータ解析のように、境界要素領域によって、解析すべきモデル全体が２つの領域に分割されるケースが多い。行列としては、中央に境界要素法で得られた係数がきて、前後に有限要素法で得られた係数を配置することで、FEM-BEM全体行列の前処理を効率的に行うことができ、かつ並列計算における計算機間の通信料を少なくすることができる。【００３４】
有限要素法等を使用する数値解析では、大規模な連立方程式を作成し、解く必要がある。連立方程式は、Aを係数行列、xを未知数の列ベクトル、bを非斉次項の列ベクトルとして、Ax=bの形で表現でき、Aの性質によって解法が異なる。Aを行列として表現すると、対称行列の場合と非対称行列の場合の大きく分けて２つの場合がある。また、行列Aの成分が０の領域が多い疎行列、０以外の領域が多い密行列、そして、疎と密がそれぞれ存在する疎密行列の３種類が存在する。一般に有限要素法解析では疎行列、境界要素法解析では密行列、有限要素法と境界要素法の混合解析では疎密行列を解く必要がある。【００３５】
連立方程式を行列で表現すると以下のようになる。
【数１】
2008026261.tif 000003
図８は、FEMで得られた連立方程式とBEMで得られた連立方程式の連立方法を説明するフローチャートである。【００３６】
以下のような手順により、連立方程式を得ることができる。
ステップＳ１０：連立方程式係数の初期化Aij＝０、Bi＝０
ステップＳ１１：FEMに基づく連立方程式の作成
解析領域全ての格子点に要素（インデックス）番号を付加する。
各要素の辺i-辺j番目の係数を計算してAijに追加する。
磁性体要素の場合Biに値を追加する。
電流要素の場合Biに値を追加する。
ステップＳ１２：BEMに基づく連立方程式の作成
BEM境界全ての辺iに対して、それ以外の辺ｊの係数を計算してAijに追加する。
ステップＳ１３：連立方程式作成完了
FEM連立方程式係数の作成過程における係数計算は以下のようにして行う。
透磁場率：μ、電流ベクトル：
2008026261.tif 000004
、磁石の磁化ベクトル：
2008026261.tif 000005
各要素に対して
【数２】
2008026261.tif 000006
を計算しながら全ての要素に対しする大規模なAijを計算する。ここで、Ｎは、モデルの形状を現す補間関数であり、補間関数をモデルの格子に沿って離散化して、列ベクトルとしたものである。ｄＶは体積要素である。【００３７】
電流要素に対しては
【数３】
2008026261.tif 000007
を計算しながら大規模なBiを作成する。【００３８】
磁性体要素に対しては
【数４】
2008026261.tif 000008
を計算しながら大規模なBiを作成する。
並べ替え処理は、以下の手順で行う。
（１）BEMでモデル化される領域で分離されるFEMでモデル化される領域AとBを作成する。
（２）BEMでモデル化された境界を構成する面領域Cの辺をID＝−１とする。
（３）面領域Cの重心を計算する。
（４）領域Aに属し、かつ、重心から最も離れた位置の辺の番号Nを取得する。
（５）領域Bに属し、かつ、N番目の辺から最も離れた位置の辺の番号Mを取得する。
（６）全ての、辺番号が与えられている辺のIDを０に初期化する。
（７）領域Aの辺番号Nを開始番号とし、この辺のIDを１とする。
（８）ID＝１の辺に隣接する辺番号が０の場合、この辺をID＝２とする。
（９）ID=nの辺に属する辺番号が０の場合、辺ID＝n+1とする。
（１０）領域Aの辺ID＝０が無くなるまで繰り返し、最後のID＝Xを記憶する。
（１１）領域Bの辺番号Mを開始番号とし、この辺のIDを１とする。
（１２）ID＝１の辺に隣接する辺番号が０の場合、この辺をID＝２とする。
（１３）ID=nの辺に属する辺番号が０の場合、辺ID＝n+1とする。
（１４）領域Bの辺ID＝０が無くなるまで繰り返す。
（１５）ID=Xを初期値として、逆の並べ替えをする（領域A）。
（１６）ID＝−１の辺領域を順番に並べ替える（面領域C）
（１７）領域Bを面領域Cの後に順番に並べる（領域B）
並べ替えは、BEM領域とそれ以外のFEM領域に対して処理を行う。FEM領域に対しては、Cuthill-Mckeeの並べ替え手法を用いる。【００３９】
図９〜図１２は、Cuthill-Mckeeの並べ替え手法を説明する図である。
FEM+BEMの並べ替えでは、BEMで必要とする未知数以外のデータに対しては、Cuthill-Mckee並べ替えによって並べ替えを行う。図９に示されるように、データ構造は、節点データは接続する辺番号および辺で結ばれる隣接節点番号をトポロジーデータとして持つ。磁場計算では、辺に未知数を配置するため、並べ替えに使用する並べ替え番号を辺に配置する。
＜節点データ構造＞
struct node_tpg{
int node_num[4];   //点⇒隣接点
int edge_num[4];   //点⇒辺
};
辺データは、辺両端の節点番号をトポロジーデータとして持つ。
＜辺データ構造＞
struct edge_tpg
{
int node_num[2];  //辺⇒2点
int n_order;       //並べ替え番号
};
Cuthill-Mckee並べ替えは、開始となる辺番号を指定し、辺が持つ節点が接続する辺番号の順で並べ替えを行う。近接する辺データを同心円状に広げていく過程で、参照するリストと蓄えるためのリストが必要となるため、リストデータを２つ準備し、それぞれのリストをバックグラウンドリスト、フォアグラウンドリストとして切り替え可能とする。【００４０】
初期化では、全ての辺番号の並べ替え番号を-１とし、-１以外番号が設定されているデータは、既に並べ替えされたデータとして判定する。初めに開始となる一つの辺番号を入力し、その辺番号を開始位置として、辺→節点→辺→節点・・・と順番に辺を呼び出し、呼び出される順に辺がナンバーリングされる。【００４１】
具体的には、既に追加されているフォアグラウンドのリストの辺に対して、両端の節点番号N1,N2を取得し、それぞれの持つ辺番号(未だ並べ替えされていない)をバックグラウンドリストに追加する。バックグラウンドリストに追加する際に、n_order=Nedgeを代入し、Nedge=Nedge+1とする。【００４２】
図１０〜図１２のフローチャートに従って説明する。
まず、ステップＳ２０において、すべての辺番号の並べ替え番号を−１とする。２つのリストをバックグラウンドリストとする。ステップＳ２１において、エッジ番号を一つのバックグラウンドリストに追加し、別のリストをフォアグラウンドリストとする。また、作業変数Ｎａｌｌを１に初期化する。ステップＳ２２において、フォアグラウンドリストとバックグラウンドリストをバックグラウンドリストとし、リストデータを空にする。バックグラウンドリストをフォアグラウンドリストとする。ステップＳ２３において、フォアグラウンドイリスとからエッジ番号を取得し、Ｎｅｄｇｅに設定する（Ｎｅｄｇｅ０とする）。【００４３】
ステップＳ２４において、取得された辺の両端ノードＮ１、Ｎ２を取得する。ステップＳ２５において、ノードＮ１の辺を取得し、Ｉ＝０とする。ステップＳ２６において、ＮｅｄｇｅがＩより大きいか否かを判断する。ステップＳ２６の判断がＮｏの場合には、図１２に進む。ステップＳ２６の判断がＹｅｓの場合には、ステップＳ２７において、ノードＮ１の辺（Ｎｅｄｇｅ１とする）を取得する。ステップＳ２８において、辺（Ｎｅｄｇｅ１）の並べ替え番号が０より小さいか否か、すなわち、並べ替え済みか否かを判断する。ステップＳ２８の判断がＮｏの場合には、Ｉを１増加してステップＳ２６に戻る。ステップＳ２８の判断がＹｅｓの場合には、ステップＳ２９において、この辺番号をバックグラウンドリストに追加し、この辺（Ｎｅｄｇｅ１）の並べ替え番号をＮａｌｌにし、Ｎａｌｌを１増加し、Ｉを１増加して、ステップＳ２６に戻る。【００４４】
ステップＳ３０において、ノードＮ１の辺を取得し、Ｉ＝０とする。ステップＳ３１において、ＮｅｄｇｅがＩより大きいか否かを判断する。ステップＳ３１の判断がＮｏの場合には、ステップＳ３５に進む。ステップＳ３１の判断がＹｅｓの場合には、ステップＳ３２において、ノードＮ１の辺（Ｎｅｄｇｅ１とする）を取得する。ステップＳ３３において、辺（Ｎｅｄｇｅ１）の並べ替え番号が０より小さいか否か、すなわち、並べ替え済みか否かを判断する。ステップＳ３３の判断がＮｏの場合には、Ｉを１増加してステップＳ３１に戻る。ステップＳ３３の判断がＹｅｓの場合には、ステップＳ３４において、この辺番号をバックグラウンドリストに追加し、この辺（Ｎｅｄｇｅ１）の並べ替え番号をＮａｌｌにし、Ｎａｌｌを１増加し、Ｉを１増加して、ステップＳ３１に戻る。【００４５】
ステップＳ３５では、バックグラウンドリストが空か否かを判断する。ステップＳ３５の判断がＮｏの場合には、図１０のＡに戻る。ステップＳ３５の判断がＹｅｓの場合には、処理を終了する。【００４６】
図１３〜図１７は、BEM領域に対する並べ替え処理を説明する図である。
この処理は、境界要素法の処理を行うプロセスが連立方程式を生成する前に、境界要素法を適用する境界部分の格子の配列番号を並べ替え、この並べ替えられた配列番号に基づいてプロセスが連立方程式を生成すれば、係数値の大きな連立方程式の係数が行列の対角要素付近に集まるようになるというものである。【００４７】
BEMで作成される方程式は、辺同士が互いに近い場合に係数が大きくなる。連立方程式の対角部分に大きな数値を配置するためには、幾何学的に近い位置にある辺番号が近くなるように工夫を行う必要がある。そこで、BEM境界面に対して背景にセル状の格子をつくり、各セルに対応する辺を登録可能なリストを作成し、このリストのセルの中に辺の中心位置が存在する辺をリストの当該セルに代入することで並べ替えを行う。【００４８】
図１３では、BEM境界は、モデルＡとモデルＢのＹ−Ｚ面に存在している。そのBEM境界に背景セルを配置した様子が図１４である。
図１４においては、グリッド座標の原点は(x0,y0)とし、Z軸方向にdz刻み、Y軸方向にdy刻みでグリッド座標値が増加する。位置(i,j)のグリッド座標値は、(z0+i*dz, y0+j*y0)となる。また、セル（I,J）は、Z座標がz0+I*dz〜z0+(I+1)*dz、Y座標がy0+j*y0〜y0+(j+1)*y0の範囲の領域である。
ゆえに、追加する辺の座標値(Z,Y)が分かれば、セルの(I,J)が以下の式で計算できる。
I=(Z-z0)/dzの整数部
J=(Y-z0)/dyの整数部
図１５は、BEM境界を構成する辺の並べ替えの全体フローチャートである。【００４９】
ステップＳ４０で、上記のように配列リスト（背景のセル状リスト）へ辺番号を代入する。ステップＳ４１で、配列リストを用いた並べ替えを行う。
図１６は、配列リストへの辺番号代入の処理を詳細に示したフローチャートである。【００５０】
ステップＳ４５において、コンピュータのキャッシュにＩｍａｘ×Ｊｍａｘのセルからなる配列リストのデータエリアを確保し、その値を空にし、変数Ｎを格納する領域をレジスタに確保し、このレジスタ値を０にする。ステップＳ４６において、BEM境界を構成する辺の数を、境界要素法を実行するプロセスから取得し、他のレジスタにＮｅｄｇｅ＿ａｌｌとして設定する。ステップＳ４７において、境界要素法を実行するプロセスが生成した、格子化されたBEM境界を構成する辺番号を取得する。通常は、この辺番号はばらばらな配列となっている。この辺番号をＮｅｄｇｅとする。ステップＳ４８において、Ｎｅｄｇｅの辺番号を有する辺の座標（ｚ、ｙ）を、境界要素法を実行するプロセスが生成したモデルから取得し、配列リストを有する背景セルの座標（Ｉ、Ｊ）を図１４で説明した計算式を演算器で計算することによって取得する。ステップＳ４９において、背景セル座標（Ｉ、Ｊ）の配列リストに辺番号Ｎｅｄｇｅを追加する。ステップＳ５０において、ＮがＮｅｄｇｅ＿ａｌｌより小さいか否かを判断する。ステップＳ５０の判断がＹｅｓの場合には、Ｎを１増加し、ステップＳ４７に戻る。ステップＳ５０の判断がＮｏの場合には、処理を終了する。【００５１】
図１７は、配列リストを用いた並べ替えの処理を詳細に示したフローチャートである。
ステップＳ５５において、３つのレジスタに変数Ｎ、Ｉ、Ｊを設定し、Ｎ＝０、Ｉ＝０、Ｊ＝０と初期化する。ステップＳ５６において、キャッシュの配列（Ｉ、Ｊ）の配列リストが空か否かを判断する。ステップＳ５６の判断がＹｅｓの場合には、ステップＳ５８に進む。ステップＳ５６の判断がＮｏの場合には、ステップＳ５７において、キャッシュの配列リストから辺番号Ｎｅｄｇｅを取得し、ＮｅｄｇｅをＮに置き換える。そして、Ｎを１増加してステップＳ５６に戻る。Ｎｅｄｇｅは、一般に境界要素法を行うプロセスがランダムに境界の格子点に与えた番号である。これを順次１ずつ増加するＮに置き換えることにより、格子点に与えられる番号が順序良く並べられるようになる。【００５２】
ステップＳ５８において、レジスタのＩを１増加し、ステップＳ５９において、ＩがＩｍａｘより小さいか否かを判断する。ステップＳ５９の判断がＹｅｓの場合、ステップＳ５６に戻る。ステップＳ５９の判断がＮｏの場合、ステップＳ６０において、レジスタのＪを１増加し、ステップＳ６１において、ＪがＪｍａｘより小さいか否かを判断する。ステップＳ６１の判断がＹｅｓの場合には、ステップＳ５６に戻る。ステップＳ６１の判断がＮｏの場合には、処理を終了する。【００５３】
図１８は、前処理を説明するための図である。
非対称行列を解くための反復手法にBiCG法が公知の技術として存在する。反復回数を減らすために前処理を行う手法が一般に用いられている。Ax=bという方程式に対して、正方行列AをLDUという行列の積で分解すると以下のようになる。【００５４】
【数５】
2008026261.tif 000009
ここで、L,D,Uは以下の性質を持つ。【００５５】
【数６】
2008026261.tif 000010
ここで対角成分は、
【数７】
2008026261.tif 000011
非対角成分(i<j)は、
【数８】
2008026261.tif 000012
また、非対角成分(i>j)は、
【数９】
2008026261.tif 000013
このような演算により完全にL,D,U行列を計算してしまうと、直接解法のように演算量が多くなってしまう。そこで、不完全コレスキー分解では以下のような処理を行う。
Aij=0  ⇒  Lij=0、Uij=0
行列Aが疎行列であれば、L、Uの成分が０の領域がおくなり、L、Uを求める過程の演算量が少なくなる。しかしながら、BEMの連立方程式の行列Aは完全に密(≠0)になるため、前処理が完全なコレスキー分解となってしまうために、演算量が多くなる。そこで、予め設定した前処理用のバンド幅2S領域を用いて以下のような不完全なLDU分解を前処理として採用する。【００５６】
Aij=0またはj-i>Sまたはi-j>Sのとき、Uij=0,Lij=0
すなわち、図１８で、演算をバンド幅２Ｓに限定することを意味する。【００５７】
【数１０】
2008026261.tif 000014
なお、以上の前処理では、バンド幅を規定して、対角要素を挟んだバンド幅内の係数についてのみ前処理を行っているが、別の方法として、係数の値に閾値を設け、係数の値が閾値より大きい係数についてのみ前処理を行うことも可能である。【００５８】
図１９及び図２０は、マルチCPU型の並列計算機において本発明の実施形態の計算を行う場合の各CPUへのデータの振り分け方法を説明するフローチャートである。
計算負荷は連立方程式Ax=bの行列Aのゼロ以外成分の大きさで決まる。各CPUには、Aijのｉ成分の範囲を割り当てる。【００５９】
図１９、及び、図２０のフローチャートでは、Aijの≠０成分をＮでカウントしながら、AN[N]=Iに縦の番号を代入する。
ステップＳ７０において、Ｉ＝０、Ｊ＝０、Ｎ＝０と変数を初期化する。ステップＳ７１において、CPUの数（Ｎ＿ＣＰＵ）と、行列のサイズ（Ｉｍａｘ、Ｊｍａｘ）を計算機に入力する。ステップＳ７２において、行列のある要素Aijが０か否かを判断する。ステップＳ７２の判断がＹｅｓの場合には、ステップＳ７４に進む。ステップＳ７２の判断がＮｏの場合には、変数AN[N]にＩを代入し、Ｎを１増加して、ステップＳ７４に進む。ステップＳ７４においては、ＩがＩｍａｘより小さいか否かを判断する。ステップＳ７４の判断がＹｅｓの場合には、ステップＳ７２に進む。ステップＳ７４の判断がＮｏの場合には、ステップＳ７５において、ＪがＪｍａｘより小さいか否かを判断する。ステップＳ７５の判断がＹｅｓの場合には、ステップＳ７２に戻る。ステップＳ７５の判断がＮｏの場合には、ステップＳ７６に進む。ステップＳ７６では、Ｄ＝Ｎ／Ｎ＿ＣＰＵを計算する。これは、各ＣＰＵの演算量を簡易的に求めたものである。【００６０】
ステップＳ７７において、Ｉ＝０、Ｊ＝０、Ｍ＝０、ＣＰＵ［０］＝０と、変数を初期化する。ステップＳ７８において、ＩがＤより小さいか否かを判断する。ステップＳ７８の判断がＹｅｓの場合には、ステップＳ８０に進む。ステップＳ７８の判断がＮｏの場合には、ステップＳ７９において、Ｉ＝０、Ｊ＝Ｊ＋１、ＣＰＵ［Ｊ］＝ＡＮ［Ｎ］を計算する。ステップＳ８０において、ＩとＮを１ずつ増加し、ステップＳ８１において、ＮがＭより小さいか否かを判断する。ステップＳ８１の判断がＹｅｓの場合には、ステップＳ７８に戻る。ステップＳ８１の判断がＮｏの場合には、処理を終了する。【００６１】
ここで、Ｊ番目のＣＰＵにはＣＰＵ［Ｊ］〜ＣＰＵ［Ｊ＋１］の成分が割り当てられる。
【図２】
2008026261.tif 000015
【図８】
2008026261.tif 000016
【図９】
2008026261.tif 000017
【図１０】
2008026261.tif 000018
【図１１】
2008026261.tif 000019
【図１２】
2008026261.tif 000020
【図１４】
2008026261.tif 000021
【図１５】
2008026261.tif 000022
【図１６】
2008026261.tif 000023
【図１７】
2008026261.tif 000024
【図１９】
2008026261.tif 000025
【図２０】
2008026261.tif 000026
【図１】
2008026261.tif 000027
【図３】
2008026261.tif 000028
【図４】
2008026261.tif 000029
【図５】
2008026261.tif 000030
【図６】
2008026261.tif 000031
【図７】
2008026261.tif 000032
【図１３】
2008026261.tif 000033
【図１８】
2008026261.tif 000034
【手続補正書】
【提出日】平成20年9月9日(2008.9.9)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理方法であって、
境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得し、
該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替えることをコンピュータに実行させることを特徴とする高速演算処理方法。
【請求項２】
前記解析対象のモデルが、１つの境界によって２つの領域Ａ及び領域Ｂに分割される場合、前記係数行列の中央部分に、該１つの境界に境界要素法を用いて得られた連立方程式の係数を配置し、その前後に、該２つの領域Ａ及び領域Ｂに有限要素法を適用して得られた連立方程式の係数を配置することを特徴とする請求項１に記載の高速演算処理方法。
【請求項３】
前記係数行列において、境界要素法による係数の前部分に配置される有限要素法による前記領域Ａについての係数は、領域Ａを記述する未知数であって、かつ、前記１つの境界を記述する未知数の内、最もインデックス値が小さい未知数に関する連立方程式の係数を開始位置として、カットヒル・マキー並べ替えを行うことを特徴とする請求項２に記載の高速演算処理方法。
【請求項４】
前記係数行列において、境界要素法による係数の後ろ部分に配置される有限要素法による前記領域Ｂについての係数は、領域Ｂを記述する未知数であって、かつ、前記１つの境界を記述する未知数の内、最もインデックス値が大きい未知数に関する連立方程式の係数を開始位置として、カットヒル・マキー並べ替えを行うことを特徴とする請求項２に記載の高速演算処理方法。
【請求項５】
前記係数行列において、境界要素法による係数は、対角要素を中心とするユーザが入力した値の幅に含まれるものについてのみ前処理が施されることを特徴とする請求項１に記載の高速演算処理方法。
【請求項６】
前記係数行列において、境界要素法による係数は、ユーザが入力した値より大きな係数値を有する係数についてのみ前処理が施されることを特徴とする請求項１に記載の高速演算処理方法。
【請求項７】
前記前処理は、不完全コレスキー分解であることを特徴とする請求項５または６に記載の高速演算処理方法。
【請求項８】
複数のＣＰＵを有する並列計算機を用いて、各ＣＰＵに前記係数行列の係数を均等に配分して計算されることを特徴とする請求項１に記載の高速演算処理方法。
【請求項９】
有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理方法をコンピュータに実現させるプログラムであって、
境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得し、
該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替えることを特徴とする高速演算処理方法をコンピュータに実現させるプログラム。
【請求項１０】
有限要素法と境界要素法とを用いて、解析対象を格子状にモデル化し、連立方程式を生成した際の連立方程式を構成する係数行列の高速演算処理装置であって、
境界要素法でモデル化された境界の格子で規定される各辺について、同一平面状になく、かつ、互いの間の距離が小さい方から順に２つの辺を取得する辺取得手段と、
該互いの距離が小さい方から順に取得された２つの辺に対して得られる、連立方程式の係数行列の行列要素を、該係数行列の対角要素近傍に並べ替える並べ替え手段と、
を備えることを特徴とする高速演算処理装置。
【国際調査報告】
2008026261.tif 000035
2008026261.tif 000036
2008026261.tif 000037
2008026261.tif 000038
2008026261.tif 000039
2008026261.tif 000040
2008026261.tif 000041