本発明の実施形態では、境界要素法で扱� ��面の幾何形状より、未知数間の係数が大き� ��もの同士が行列内で近い位置にくるような� ��べ替えを行うことで、反復による収束を速� ��するための前処理の効率化を図る。

 図2は、本発明の主要部の流れを示すフロー チャートである。
 ここでは、有限要素法と境界要素法とを用� ��て、解析対象を格子状にモデル化し、連立� ��程式を生成した際の連立方程式を構成する� ��数行列の高速演算処理方法における並べ替� ��処理方法を示す。

 まず、Step1で、境界要素法でモデル化され� �境界の格子で規定される各辺について、同� �平面状になく、かつ、互いの間の距離が小� �い方から順に2つの辺を取得する。
 Step2では、前記解析対象のモデルが、1つの� ��界によって2つの領域A及び領域Bに分割され� ��場合、前記係数行列の中央部分に、該1つの 境界に境界要素法を用いて得られた連立方程 式の係数を配置し、その前後に、該2つの領� �A及び領域Bに有限要素法を適用して得られた 連立方程式の係数を配置する。

 Step3では、前記係数行列において、境界� �素法による係数の前部分に配置される有限� �素法による前記領域Aについての係数は、領� ��Aを記述する未知数であって、かつ、前記1� �の境界を記述する未知数の内、最もインデ� �クス値が小さい未知数に関する連立方程式� �係数を開始位置として、カットヒル・マキ� �(Cuthill-Mckee)並べ替えを行う。また、前記係� �行列において、境界要素法による係数の後� �部分に配置される有限要素法による前記領� �Bについての係数は、領域Bを記述する未知数 であって、かつ、前記1つの境界を記述する� �知数の内、最もインデックス値が大きい未� �数に関する連立方程式の係数を開始位置と� �て、カットヒル・マキー並べ替えを行う。

 Step4では、境界要素法の係数について、� �互いの距離が小さい方から順に取得された2� ��の辺に対して得られる、連立方程式の係数� ��列の行列要素を、該係数行列の対角要素近� ��に並べ替える。

 図3は、本発明の実施形態の概略を説明する 図である。
 境界要素法で解く領域は密領域であるため� ��反復解法を用いて前処理に不完全コレスキ� ��分解を行うと、結果的に完全なコレスキー� ��解となってしまい、直接法による解法と同� ��になってしまう。そこで、密領域の全ての� ��数を使用せずに、前処理による反復法の収� ��性を高める必要がある。

 有限要素法と境界要素法の混合解法では� ��図3左図のように行列の対角要素付近にノン ゼロの値が集中する領域1と正方形に近い形� �領域にノンゼロの値が集中する領域2が存在� ��る。有限要素法解析では一般に行列全体に� ��して不完全コレスキー分解を行うが、この� ��うな行列に対して不完全コレスキー分解を� ��うと、領域2において計算時間の増加を招い てしまい非効率である。そこで領域2に対し� �図3右図のように対角要素の近傍の値だけを� ��象に不完全コレスキー分解を行うことで、� ��処理時間の短縮化を行う。

 すなわち、本発明の実施形態においては� ��有限要素法と境界要素法に基づく連立方程� ��の合成より成り立つ行列において、全ての� ��知数中の境界要素法で使用する未知数の並� ��替えを初めに行い、続いて残りの未知数の� ��べ替えを行う。

 図4及び図5は、数値解析すべきモデルとそ� �によってできる行列について説明する図で� �る。
 図4(a)のモデルは、HDDヘッドとHDD媒体との間 の動作状態を数値解析するためのモデルであ る。このHDDヘッド11と媒体10を含む空間をそ� �ぞれモデルBとモデルAとして、格子に区切り 、各格子点について物理量を考える。モデル AとモデルBの相互作用は、格子点間の相互作� ��として記述し、相互作用そのものは物理理� ��(たとえば、電磁気学)の法則にしたがって� �述する。媒体10を含むモデルA及び、HDDヘッ� �11を含むモデルBはそれぞれ、媒体10、HDDヘッ ド11のみではなく、その周りの空気領域も含� ��てモデル化したものである。このモデルAと Bの内部は、有限要素法を用いて計算される� �、モデルAとモデルBの両方の最外周領域は、 2つのモデルの最外周を1つの外周を表すモデ� ��としてモデル化し、外周を格子に区切って� ��物理量を計算する。図4(a)の場合、媒体10とH DDヘッド11が主に相互作用する境界要素面14が 、モデルAとモデルBの相対する外周面として� ��述される。このような外周境界の計算は、� ��界要素法を用いて計算される。

 以上のようにモデル化された解析対象に� ��いて、有限要素法と境界要素法を用いて連� ��方程式を立て、その係数を行列として整理� ��ると、図4(b)のような行列となる。すなわち 、モデルAの内部を記述する係数は、図4(b)の� ��列の領域1のように対角要素の周辺にノンゼ ロの係数が集中する。同様に、モデルBの内� �を記述する係数は、領域3のように対角要素� ��周辺にノンゼロの係数が集中する。一方、� ��界をモデル化し、境界要素法で連立方程式� ��つくり、係数を行列表示すると、領域2のよ うに、その部分だけ、密行列のようになって いる。

 境界要素法で使用する未知数の並べ替えは� ��代数的にモデル端部に位置する未知数(図4(a )のa点のような点における、求めるべき物理� ��をあらわすxなどの未知数。モデルが格子化 されているので、各格子点での物理量をx _i のようにインデックスiをつけて記述する)を� ��始No=1とし、それ以外の未知数との行列成分 (連立方程式は、a _ij x _j +a _ij+1 x _j+1 +・・・=b _i のような式からなっており、x _k をNo=1とした場合、a _ik )が最も大きい未知数をNo=2として登録する。� ��後に登録された未知数No=nに対して、まだ登 録されていない未知数との行列成分が最も大 きい未知数を最後のNoにn+1番目で登録する。� ��上の手順でナンバーリングされた番号を境� ��要素法領域における連立方程式の未知数の� ��番と定義する。

 図4(a)のように、境界要素法を使用する場 合、モデルの領域が複数の領域に分割されて いる。媒体移動を伴うHDDヘッドのライト・リ ードシミュレーションのように、ヘッドと媒 体というように2つの大きな領域に分割され� �いる場合に、以下の並べ替えを行うことで� �率的な領域分割を実現する。

 連立方程式より作成される行列の領域1は モデル左側の媒体+空気領域の境界要素面以� �の未知数から構成し、領域2は境界要素面の� ��知数から構成し、領域3はモデル右側のリー ドヘッド+空気領域の境界要素面以外の未知� �から構成する。

 領域2の並べ替え後に、媒体側未知数でか つ領域2の未知数のインデックス番号が最も� �さい番号を開始番号として領域1の並べ替え� ��行う。同様に、HDDヘッド側未知数でかつ領� ��2の未知数のインデックス番号が最も大きい 番号を開始番号として領域3の並べ替えを行� �。領域1~3の並べ替え後の順番に従って、最� �に領域1→領域2→領域3の通し番号でインデ� �クス番号付けを行って並べ替えを終了する� �

 境界要素法の未知数間の係数(未知数x _i と未知数x _j 間の係数はa _ij )は、3次元解析の場合は未知数が定義される� ��の単位法線ベクトルと、問題とする未知数� ��定義される位置から、相互作用を考慮すべ� ��未知数が定義される位置への位置ベクトル� ��の内積に比例するため、同一平面に存在す� ��面内の未知数間の係数は0となる。そこで、 連立方程式作成過程における計算では、境界 要素面を構成する全ての未知数の組み合わせ に対して成分計算は行わず、予め同一平面上 に存在するという情報を面要素グループに設 定し、異なる面要素グループとの組み合わせ に含まれる面に属する未知数間の係数だけを 成分計算することで計算時間の短縮化を行う 。

 なお、図4(b)に示されるように、対角要素 近傍のノンゼロ領域の幅をバンド幅と呼ぶ。 領域2のバンド幅は広いため、前処理におい� �、そのままのバンド幅を使用すると計算時� �が増加してしまう。前処理に使用するバン� �幅は、必ずしも各領域のバンド幅と同じで� �る必要はない。そこで、図5(a)の領域2に対し てはユーザーが任意にバンド幅を指定ユーザ 処理を行う領域を図5(b)の領域2のように制限� ��ることを可能とする。

 図6及び図7は、マルチCPU計算機での計算の� �方を説明する図である。
 複数のCPUを使用する並列計算では、ユーザ� ��使用する計算機ノード数を入力し、行列を� ��6のようにブロック分割する。また、各計算 機内に複数のCPUが存在する場合、ユーザが並 列スレッド数を入力すると、計算機ノード数 ×スレッド数で行列をブロック分割する。

 要素数を単純に分割するだけでは、各CPU� ��加わる計算負荷が均等にならない。計算負� ��は、モデル形状やメッシュ(格子)分割方法� �FEM(Finite Element Method:有限要素法)領域である かBEM(Boundary Element Method:境界要素法)領域で� �るかに依存する。そこで、各未知数に対す� �一回の反復計算あたりの計算時間をテスト� �算により予め算出し、計算付加が均等にな� �ように各CPUに未知数を割り当てる。

 図6では、CPU1~5に行列の係数を上から順に 割り当てている。図6のCPU3に割り当てられた� ��数は、ユーザが制限したバンド幅がCPU3に割 り当てられた係数のインデックス幅より大き い場合の前処理サイズである。一方、図7(a)� �CPU3に割り当てられた係数は、ユーザが制限� ��たバンド幅がCPU3に割り当てられた係数のイ ンデックス幅より小さい場合の前処理サイズ である。

 ユーザが設定したBEM領域のバンド幅制限� ��イズとCPUに割り当てられた未知数サイズの� ��小関係を自動的に判断し、図7(b)のようなBEM 領域の前処理ブロックを、図7(c)のように、� �ンド幅制限サイズに自動的に制限して計算� �行う。

 有限要素法で生成される行列のノンゼロ� ��分は、同じ材料物性グループ内では、大き� ��値と小さい値との差は小さい。しかしなが� ��、境界要素法によって生成される行列成分� ��、未知数と未知数の間の係数の値は、未知� ��が配置されている位置間の距離の2乗に反比 例する(3次元解析)ため、係数の大きさは幅広 く分布している。そこで未知数間の係数が大 きい値が行列対角部分に集中するような並べ 替えを行い、対角近傍成分だけを用いて反復 解法における不完全コレスキー分解を行う。 どの範囲までの対角近傍成分幅(制限バンド� �)を使用するかは、ユーザが予め数値入力し� ��指定する。また、最大係数に対してどの値� ��でを前処理に含めるかという閾値をユーザ� ��指定して、制限バンド幅を計算することを� ��能にする。

 HDDのヘッド-媒体解析やモータ解析のよう に、境界要素領域によって、解析すべきモデ ル全体が2つの領域に分割されるケースが多� �。行列としては、中央に境界要素法で得ら� �た係数がきて、前後に有限要素法で得られ� �係数を配置することで、FEM-BEM全体行列の前� ��理を効率的に行うことができ、かつ並列計� ��における計算機間の通信料を少なくするこ� ��ができる。

 有限要素法等を使用する数値解析では、� ��規模な連立方程式を作成し、解く必要があ� ��。連立方程式は、Aを係数行列、xを未知数� �列ベクトル、bを非斉次項の列ベクトルとし� ��、Ax=bの形で表現でき、Aの性質によって解� �が異なる。Aを行列として表現すると、対称� ��列の場合と非対称行列の場合の大きく分け� ��2つの場合がある。また、行列Aの成分が0の� ��域が多い疎行列、0以外の領域が多い密行列 、そして、疎と密がそれぞれ存在する疎密行 列の3種類が存在する。一般に有限要素法解� �では疎行列、境界要素法解析では密行列、� �限要素法と境界要素法の混合解析では疎密� �列を解く必要がある。

 連立方程式を行列で表現すると以下のよう� ��なる。
 図8は、FEMで得られた連立方程式とBEMで得ら れた連立方程式の連立方法を説明するフロー チャートである。

 以下のような手順により、連立方程式を得� ��ことができる。
ステップS10:連立方程式係数の初期化Aij=0、Bi= 0
ステップS11:FEMに基づく連立方程式の作成
解析領域全ての格子点に要素(インデックス)� ��号を付加する。
各要素の辺i-辺j番目の係数を計算してAijに追 加する。
磁性体要素の場合Biに値を追加する。
電流要素の場合Biに値を追加する。
ステップS12:BEMに基づく連立方程式の作成
BEM境界全ての辺iに対して、それ以外の辺jの� ��数を計算してAijに追加する。
ステップS13:連立方程式作成完了
 FEM連立方程式係数の作成過程における係数� ��算は以下のようにして行う。
透磁場率:μ、電流ベクトル:
、磁石の磁化ベクトル:
 各要素に対して
を計算しながら全ての要素に対しする大規模 なAijを計算する。ここで、Nは、モデルの形� �を現す補間関数であり、補間関数をモデル� �格子に沿って離散化して、列ベクトルとし� �ものである。dVは体積要素である。

 電流要素に対しては
を計算しながら大規模なBiを作成する。

 磁性体要素に対しては
を計算しながら大規模なBiを作成する。
 並べ替え処理は、以下の手順で行う。
(1)BEMでモデル化される領域で分離されるFEMで モデル化される領域AとBを作成する。
(2)BEMでモデル化された境界を構成する面領域 Cの辺をID=-1とする。
(3)面領域Cの重心を計算する。
(4)領域Aに属し、かつ、重心から最も離れた� �置の辺の番号Nを取得する。
(5)領域Bに属し、かつ、N番目の辺から最も離� ��た位置の辺の番号Mを取得する。
(6)全ての、辺番号が与えられている辺のIDを0 に初期化する。
(7)領域Aの辺番号Nを開始番号とし、この辺のI Dを1とする。
(8)ID=1の辺に隣接する辺番号が0の場合、この� ��をID=2とする。
(9)ID=nの辺に属する辺番号が0の場合、辺ID=n+1� ��する。
(10)領域Aの辺ID=0が無くなるまで繰り返し、最 後のID=Xを記憶する。
(11)領域Bの辺番号Mを開始番号とし、この辺の IDを1とする。
(12)ID=1の辺に隣接する辺番号が0の場合、この 辺をID=2とする。
(13)ID=nの辺に属する辺番号が0の場合、辺ID=n+1 とする。
(14)領域Bの辺ID=0が無くなるまで繰り返す。
(15)ID=Xを初期値として、逆の並べ替えをする( 領域A)。
(16)ID=-1の辺領域を順番に並べ替える(面領域C)
(17)領域Bを面領域Cの後に順番に並べる(領域B)
 並べ替えは、BEM領域とそれ以外のFEM領域に� ��して処理を行う。FEM領域に対しては、Cuthill -Mckeeの並べ替え手法を用いる。

 図9~図12は、Cuthill-Mckeeの並べ替え手法を説� �する図である。
 FEM+BEMの並べ替えでは、BEMで必要とする未知 数以外のデータに対しては、Cuthill-Mckee並べ� �えによって並べ替えを行う。図9に示される� ��うに、データ構造は、節点データは接続す� ��辺番号および辺で結ばれる隣接節点番号を� ��ポロジーデータとして持つ。磁場計算では� ��辺に未知数を配置するため、並べ替えに使� ��する並べ替え番号を辺に配置する。
<節点データ構造>
struct node_tpg{
  int node_num[4];   //点⇒隣接点
  int edge_num[4];   //点⇒辺
};
辺データは、辺両端の節点番号をトポロジー データとして持つ。
<辺データ構造>
struct edge_tpg
{
  int node_num[2];  //辺⇒2点
  int n_order;       //並べ替え番号
};
 Cuthill-Mckee並べ替えは、開始となる辺番号を 指定し、辺が持つ節点が接続する辺番号の順 で並べ替えを行う。近接する辺データを同心 円状に広げていく過程で、参照するリストと 蓄えるためのリストが必要となるため、リス トデータを2つ準備し、それぞれのリストを� �ックグラウンドリスト、フォアグラウンド� �ストとして切り替え可能とする。

 初期化では、全ての辺番号の並べ替え番� ��を-1とし、-1以外番号が設定されているデー タは、既に並べ替えされたデータとして判定 する。初めに開始となる一つの辺番号を入力 し、その辺番号を開始位置として、辺→節点 →辺→節点・・・と順番に辺を呼び出し、呼 び出される順に辺がナンバーリングされる。

 具体的には、既に追加されているフォア� ��ラウンドのリストの辺に対して、両端の節� ��番号N1,N2を取得し、それぞれの持つ辺番号(� ��だ並べ替えされていない)をバックグラウン ドリストに追加する。バックグラウンドリス トに追加する際に、n_order=Nedgeを代入し、Nedge =Nedge+1とする。

 図10~図12のフローチャートに従って説明す� �。
 まず、ステップS20において、すべての辺番� ��の並べ替え番号を-1とする。2つのリストを� ��ックグラウンドリストとする。ステップS21� ��おいて、エッジ番号を一つのバックグラウ� ��ドリストに追加し、別のリストをフォアグ� ��ウンドリストとする。また、作業変数Nallを 1に初期化する。ステップS22において、フォ� �グラウンドリストとバックグラウンドリス� �をバックグラウンドリストとし、リストデ� �タを空にする。バックグラウンドリストを� �ォアグラウンドリストとする。ステップS23� �おいて、フォアグラウンドイリスとからエ� �ジ番号を取得し、Nedgeに設定する(Nedge0とす� �)。

 ステップS24において、取得された辺の両� ��ノードN1、N2を取得する。ステップS25におい て、ノードN1の辺を取得し、I=0とする。ステ� ��プS26において、NedgeがIより大きいか否かを� ��断する。ステップS26の判断がNoの場合には� �図12に進む。ステップS26の判断がYesの場合に は、ステップS27において、ノードN1の辺(Nedge1 とする)を取得する。ステップS28において、� �(Nedge1)の並べ替え番号が0より小さいか否か� �すなわち、並べ替え済みか否かを判断する� �ステップS28の判断がNoの場合には、Iを1増加� ��てステップS26に戻る。ステップS28の判断がY esの場合には、ステップS29において、この辺� ��号をバックグラウンドリストに追加し、こ� ��辺(Nedge1)の並べ替え番号をNallにし、Nallを1� �加し、Iを1増加して、ステップS26に戻る。

 ステップS30において、ノードN1の辺を取� �し、I=0とする。ステップS31において、Nedgeが Iより大きいか否かを判断する。ステップS31� �判断がNoの場合には、ステップS35に進む。ス テップS31の判断がYesの場合には、ステップS32 において、ノードN1の辺(Nedge1とする)を取得� �る。ステップS33において、辺(Nedge1)の並べ替 え番号が0より小さいか否か、すなわち、並� �替え済みか否かを判断する。ステップS33の� �断がNoの場合には、Iを1増加してステップS31� ��戻る。ステップS33の判断がYesの場合には、� ��テップS34において、この辺番号をバックグ� ��ウンドリストに追加し、この辺(Nedge1)の並� �替え番号をNallにし、Nallを1増加し、Iを1増加 して、ステップS31に戻る。

 ステップS35では、バックグラウンドリス� ��が空か否かを判断する。ステップS35の判断� ��Noの場合には、図10のAに戻る。ステップS35� �判断がYesの場合には、処理を終了する。

 図13~図17は、BEM領域に対する並べ替え処理� �説明する図である。
 この処理は、境界要素法の処理を行うプロ� ��スが連立方程式を生成する前に、境界要素� ��を適用する境界部分の格子の配列番号を並� ��替え、この並べ替えられた配列番号に基づ� ��てプロセスが連立方程式を生成すれば、係� ��値の大きな連立方程式の係数が行列の対角� ��素付近に集まるようになるというものであ� ��。

 BEMで作成される方程式は、辺同士が互い� ��近い場合に係数が大きくなる。連立方程式� ��対角部分に大きな数値を配置するためには� ��幾何学的に近い位置にある辺番号が近くな� ��ように工夫を行う必要がある。そこで、BEM� ��界面に対して背景にセル状の格子をつくり� ��各セルに対応する辺を登録可能なリストを� ��成し、このリストのセルの中に辺の中心位� ��が存在する辺をリストの当該セルに代入す� ��ことで並べ替えを行う。

 図13では、BEM境界は、モデルAとモデルBのY-Z 面に存在している。そのBEM境界に背景セルを 配置した様子が図14である。
 図14においては、グリッド座標の原点は(x0,y 0)とし、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の整数部
 図15は、BEM境界を構成する辺の並べ替えの� �体フローチャートである。

 ステップS40で、上記のように配列リスト(背 景のセル状リスト)へ辺番号を代入する。ス� �ップS41で、配列リストを用いた並べ替えを� �う。
 図16は、配列リストへの辺番号代入の処理� �詳細に示したフローチャートである。

 ステップS45において、コンピュータのキ� ��ッシュにImax×Jmaxのセルからなる配列リスト のデータエリアを確保し、その値を空にし、 変数Nを格納する領域をレジスタに確保し、� �のレジスタ値を0にする。ステップS46におい� ��、BEM境界を構成する辺の数を、境界要素法� ��実行するプロセスから取得し、他のレジス� ��にNedge_allとして設定する。ステップS47にお� ��て、境界要素法を実行するプロセスが生成� ��た、格子化されたBEM境界を構成する辺番号� ��取得する。通常は、この辺番号はばらばら� ��配列となっている。この辺番号をNedgeとす� �。ステップS48において、Nedgeの辺番号を有す る辺の座標(z、y)を、境界要素法を実行する� �ロセスが生成したモデルから取得し、配列� �ストを有する背景セルの座標(I、J)を図14で� �明した計算式を演算器で計算することによ� �て取得する。ステップS49において、背景セ� �座標(I、J)の配列リストに辺番号Nedgeを追加� �る。ステップS50において、NがNedge_allより小� ��いか否かを判断する。ステップS50の判断がY esの場合には、Nを1増加し、ステップS47に戻� �。ステップS50の判断がNoの場合には、処理を 終了する。

 図17は、配列リストを用いた並べ替えの処� �を詳細に示したフローチャートである。
 ステップS55において、3つのレジスタに変数 N、I、Jを設定し、N=0、I=0、J=0と初期化する。 ステップS56において、キャッシュの配列(I、J )の配列リストが空か否かを判断する。ステ� �プS56の判断がYesの場合には、ステップS58に� �む。ステップS56の判断がNoの場合には、ステ ップS57において、キャッシュの配列リストか ら辺番号Nedgeを取得し、NedgeをNに置き換える� ��そして、Nを1増加してステップS56に戻る。Ne dgeは、一般に境界要素法を行うプロセスがラ ンダムに境界の格子点に与えた番号である。 これを順次1ずつ増加するNに置き換えること� ��より、格子点に与えられる番号が順序良く� ��べられるようになる。

 ステップS58において、レジスタのIを1増� �し、ステップS59において、IがImaxより小さい か否かを判断する。ステップS59の判断がYesの 場合、ステップS56に戻る。ステップS59の判断 がNoの場合、ステップS60において、レジスタ� ��Jを1増加し、ステップS61において、JがJmaxよ り小さいか否かを判断する。ステップS61の判 断がYesの場合には、ステップS56に戻る。ステ ップS61の判断がNoの場合には、処理を終了す� ��。

 図18は、前処理を説明するための図である� �
 非対称行列を解くための反復手法にBiCG法が 公知の技術として存在する。反復回数を減ら すために前処理を行う手法が一般に用いられ ている。Ax=bという方程式に対して、正方行� �AをLDUという行列の積で分解すると以下のよ� ��になる。

ここで、L,D,Uは以下の性質を持つ。

ここで対角成分は、
非対角成分(i<j)は、
また、非対角成分(i>j)は、
 このような演算により完全にL,D,U行列を計� �してしまうと、直接解法のように演算量が� �くなってしまう。そこで、不完全コレスキ� �分解では以下のような処理を行う。
Aij=0 ⇒ Lij=0、Uij=0
 行列Aが疎行列であれば、L、Uの成分が0の領 域がおくなり、L、Uを求める過程の演算量が� ��なくなる。しかしながら、BEMの連立方程式� ��行列Aは完全に密(≠0)になるため、前処理が 完全なコレスキー分解となってしまうために 、演算量が多くなる。そこで、予め設定した 前処理用のバンド幅2S領域を用いて以下のよ� ��な不完全なLDU分解を前処理として採用する� ��

 Aij=0またはj-i>Sまたはi-j>Sのとき、Uij=0,L ij=0
すなわち、図18で、演算をバンド幅2Sに限定� �ることを意味する。

 なお、以上の前処理では、バンド幅を規定� ��て、対角要素を挟んだバンド幅内の係数に� ��いてのみ前処理を行っているが、別の方法� ��して、係数の値に閾値を設け、係数の値が� ��値より大きい係数についてのみ前処理を行� ��ことも可能である。

 図19及び図20は、マルチCPU型の並列計算機に おいて本発明の実施形態の計算を行う場合の 各CPUへのデータの振り分け方法を説明するフ ローチャートである。
 計算負荷は連立方程式Ax=bの行列Aのゼロ以� �成分の大きさで決まる。各CPUには、Aijのi成� ��の範囲を割り当てる。

 図19、及び、図20のフローチャートでは、Aij の≠0成分をNでカウントしながら、AN[N]=Iに縦 の番号を代入する。
 ステップS70において、I=0、J=0、N=0と変数を� ��期化する。ステップS71において、CPUの数(N_C PU)と、行列のサイズ(Imax、Jmax)を計算機に入� �する。ステップS72において、行列のある要� �A _ij が0か否かを判断する。ステップS72の判断がYe sの場合には、ステップS74に進む。ステップS7 2の判断がNoの場合には、変数AN[N]にIを代入し 、Nを1増加して、ステップS74に進む。ステッ� ��S74においては、IがImaxより小さいか否かを� �断する。ステップS74の判断がYesの場合には� �ステップS72に進む。ステップS74の判断がNoの 場合には、ステップS75において、JがJmaxより� ��さいか否かを判断する。ステップS75の判断� ��Yesの場合には、ステップS72に戻る。ステッ� ��S75の判断がNoの場合には、ステップS76に進� �。ステップS76では、D=N/N_CPUを計算する。こ� �は、各CPUの演算量を簡易的に求めたもので� �る。

 ステップS77において、I=0、J=0、M=0、CPU[0]= 0と、変数を初期化する。ステップS78におい� �、IがDより小さいか否かを判断する。ステッ プS78の判断がYesの場合には、ステップS80に進 む。ステップS78の判断がNoの場合には、ステ� ��プS79において、I=0、J=J+1、CPU[J]=AN[N]を計算� �る。ステップS80において、IとNを1ずつ増加� �、ステップS81において、NがMより小さいか否 かを判断する。ステップS81の判断がYesの場合 には、ステップS78に戻る。ステップS81の判断 がNoの場合には、処理を終了する。

 ここで、J番目のCPUにはCPU[J]~CPU[J+1]の成分 が割り当てられる。