(57)【特許請求の範囲】
【請求項１】
磁場の連立方程式の解を反復計算で計算する磁化解析装置において、
前記連立方程式を構成する係数行列を分解して前処理行列を取得する取得部と、
取得された前記前処理行列の対角成分の値を調整するパラメータを示すシフトパラメータの値を２分探索に基づいて算出する算出部と、を有し、
前記算出部は、
前記シフトパラメータの値の第１の値および第２の値を予め設定する第１の設定部と、
設定された前記第１の値および設定された前記第２の値の中間値に対応する、前記前処理行列の対角成分の最小値を算出する第１の算出部と、
設定された前記第１の値と設定された前記第２の値との差分を算出する第２の算出部と、
算出された前記前処理行列の対角成分の最小値が零より大きい場合、設定された前記第１の値および設定された前記第２の値の中間値を新たな第２の値に設定し、前記最小値が零以下である場合、設定された前記第１の値および前記設定された第２の値との中間値を新たな第１の値に設定する第２の設定部と、を含み、
算出された前記差分が所定の閾値より小さくなるまで、前記第１の算出部による処理、前記第２の算出部による処理および前記第２の設定部による処理を繰り返す
ことを特徴とする磁化解析装置。
【請求項２】
前記第１の算出部は、前記連立方程式が回路に係る連立方程式を含む場合、前記前処理行列の対角成分のうち前記回路に係る対角成分を除外した、前記前処理行列の対角成分の最小値を算出する
ことを特徴とする請求項１に記載の磁化解析装置。
【請求項３】
前記算出部は、
算出された前記差分が前記所定の閾値より小さい場合、設定された前記第１の値と設定された前記第２の値との間の値に、１に前記シフトパラメータを調整するパラメータを示す調整パラメータの値を加算して得られた加算値を乗じて得た乗算値を前記シフトパラメータの値として算出する第３の算出部をさらに含む
ことを特徴とする請求項１に記載の磁化解析装置。
【請求項４】
前記第３の算出部は、
前記所定の閾値を、前記調整パラメータとして予め定められた一変数の線形式を用いて算出する
ことを特徴とする請求項３に記載の磁化解析装置。
【請求項５】
前記算出部は、
取得された前記前処理行列の第１の対角成分の値が正になるように、前記シフトパラメータの値を２分探索に基づいて算出する第１の算出制御部と、
取得された前記前処理行列の第２の対角成分の値が正になるように、前記シフトパラメータの値を２分探索に基づいて算出する第２の算出制御部と、
前記第１の算出制御部が算出したシフトパラメータの値と、前記第２の算出制御部が算出したシフトパラメータの値の最大値を共有する
ことを特徴とする請求項３または請求項４に記載の磁化解析装置。
【請求項６】
コンピュータが、
磁場の連立方程式を構成する係数行列を分解して前処理行列を取得し、
取得された前記前処理行列の対角成分の値を調整するパラメータを示すシフトパラメータの値の第１の値および第２の値を予め設定し、
設定された前記第１の値および設定された前記第２の値の中間値に対応する、前記前処理行列の対角成分の最小値を算出するとともに、設定された前記第１の値と設定された前記第２の値との差分を算出し、
算出された前記前処理行列の対角成分の最小値が零より大きい場合、設定された前記第１の値および設定された前記第２の値の中間値を新たな第２の値に設定し、前記最小値が零以下である場合、設定された前記第１の値および前記設定された第２の値との中間値を新たな第１の値に設定し、
算出された前記差分が所定の閾値より小さくなるまで、前記算出する処理、前記中間値を新たに設定する処理を繰り返す
処理を実行することを特徴とする磁化解析方法。
【請求項７】
コンピュータに、
磁場の連立方程式を構成する係数行列を分解して前処理行列を取得し、
取得された前記前処理行列の対角成分の値を調整するパラメータを示すシフトパラメータの値の第１の値および第２の値を予め設定し、
設定された前記第１の値および設定された前記第２の値の中間値に対応する、前記前処理行列の対角成分の最小値を算出するとともに、設定された前記第１の値と設定された前記第２の値との差分を算出し、
算出された前記前処理行列の対角成分の最小値が零より大きい場合、設定された前記第１の値および設定された前記第２の値の中間値を新たな第２の値に設定し、前記最小値が零以下である場合、設定された前記第１の値および前記設定された第２の値との中間値を新たな第１の値に設定し、
算出された前記差分が所定の閾値より小さくなるまで、前記算出する処理、前記中間値を新たに設定する処理を繰り返す
処理を実行させることを特徴とする磁化解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、磁化解析装置、磁化解析方法および磁化解析プログラムに関する。
【背景技術】
【０００２】
コンピュータを用いた数値シミュレーションでは、有限の数の物理量で場を表現するために離散化と呼ばれる技術が用いられる。この離散化手法には、例えば、有限差分法、有限要素法、有限体積法が挙げられる。いずれの離散化手法も空間に有限の数の物理量を配置し、物理現象を支配する方程式を離散化して連立方程式を作成する。そして、コンピュータで連立方程式を解くことで物理量が決定される。
【０００３】
連立方程式を行列とベクトルで表現するとき、この行列は係数行列と呼ばれる。数値シミュレーションで扱う連立方程式は様々であり、係数行列の対角成分に零が存在することがある。また、解が不能であったり、不定であったりすることがある。
【０００４】
辺要素を用いる磁場解析の過程で生成される連立方程式は不定である。辺要素は、例えば、有限要素法で用いられる要素のことである。不定の連立方程式を直接解法の線形ソルバーを適用して解を求めることは困難である。しかし、ＣＧ（conjugate  gradient）法またはＢｉＣＧ（Bi-Conjugate  Gradient）法のようなクロリフ部分空間を用いる反復法ソルバーであれば、不定の連立方程式の解を求めることが可能である。
【０００５】
ＣＧ法やＢｉＣＧ法に前処理を施す手法は、それぞれＩＣＣＧ法、ＩＬＵ  ＢｉＣＧ法と呼ばれる。ＩＣＣＧ（Incomplete  Cholesky  conjugate  gradient）法、ＩＬＵ  ＢｉＣＧ（Incomplete  LU  factorization）法は、前処理を施すことで、収束までに要する反復回数を減らすことが可能となる。
【０００６】
ＣＧ法の処理は、対称行列を反復計算で解を求める。図１１は、ＣＧ法の処理の一例を示す図である。図１１に示すように、ＣＧ法の処理は、前処理をしないで、対称行列を反復計算で解を求めている。図１１では、Ａが対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。
【０００７】
ＩＣＣＧ法の処理は、前処理を施した後に、対称行列を反復計算で解を求める。図１２は、前処理付きのＣＧ（ＩＣＣＧ）法の処理の一例を示す図である。図１２に示すように、ＩＣＣＧ法の処理は、前処理として前処理行列を計算した後に、対称行列を反復計算で解を求めている。図１２では、Ｃが前処理行列であり、Ａが対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。
【０００８】
Ａが対称行列である場合の前処理行列Ｃの計算は、以下の式（１）で表わされる。対称行列Ａを不完全コレスキー分解した行列Ｃは、対称行列Ａの近似である。
【数１】
0006337682.tif 000002
式（１）によると、前処理行列Ｃは、非対角成分の行列Ｌと対角成分の行列Ｄと非対角成分の行列Ｌの転置行列ＬＴにより表わされる。
【０００９】
行列Ｄの対角成分は、以下の式（２）で表わされる。
【数２】
0006337682.tif 000003
【００１０】
行列Ｌの非対角成分は、以下の式（３）で表わされる。
【数３】
0006337682.tif 000004
【００１１】
また、非対称行列を反復計算で解を求める手法には、ＢｉＣＧ法、ＢｉＣＧ  Ｓｔａｂ法、ＣＧＳ法やＣＲ法等が挙げられる。
【００１２】
ＣＲ法の処理は、非対称行列を反復処理で解を求める。図１３は、ＣＲ法の処理の一例を示す図である。図１３に示すように、ＣＲ法の処理は、前処理をしないで、非対称行列を反復計算で解を求めている。図１３では、Ａが非対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。
【００１３】
ＩＬＵ  ＣＲ法の処理は、前処理を施した後に、非対称行列を反復計算で解を求める。図１４は、前処理付きのＣＲ（ＩＬＵ  ＣＲ）法の処理の一例を示す図である。図１４に示すように、ＩＬＵ  ＣＲ法の処理は、前処理として前処理行列を計算した後に、非対称行列を反復計算で解を求めている。図１４では、Ｃが前処理行列であり、Ａが非対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。
【００１４】
Ａが非対称行列である場合の前処理行列Ｃの計算は、以下の式（４）で表わされる。非対称行列Ａを不完全ＬＵ分解した行列Ｃは、非対称行列Ａの近似である。
【数４】
0006337682.tif 000005
式（４）によると、前処理行列Ｃは、非対角成分の行列Ｌと対角成分の行列Ｄと非対角成分の行列Ｕにより表わされる。
【００１５】
行列Ｄの対角成分は、以下の式（５）で表わされる。
【数５】
0006337682.tif 000006
【００１６】
行列Ｌの非対角成分は、以下の式（６）で表わされる。
【数６】
0006337682.tif 000007
【００１７】
行列Ｕの非対角成分は、以下の式（７）で表わされる。
【数７】
0006337682.tif 000008
【００１８】
ところで、辺要素を用いる磁場解析の過程で生成される連立方程式の係数行列Ａを用いて前処理行列Ｃの対角成分Ｄｉｉを計算すると、零または負になる成分が発生する。零または負になる成分が発生すると、前処理が破綻するため、その後の反復計算が発散してしまう。前処理行列Ｃは、係数行列Ａの近似であるため、何らかのパラメータによって調整することが可能である。
【００１９】
従来では、以下の式（８）や式（９）のように１より大きな係数γによって、全てのｉに関して対角成分Ｄｉｉが零以下にならないように調整する方法が知られている。式（８）は、ＣＧ法の場合の調整する方法である。式（９）は、ＣＲ法の場合の調整する方法である。なお、係数γは、シフトパラメータと呼ばれる。
【数８】
0006337682.tif 000009
【００２０】
【数９】
0006337682.tif 000010
【００２１】
シフトパラメータγの値は、連立方程式を解くための反復計算の反復回数に影響を与えるため、適切に調整することで計算時間を大幅に削減することができる。シフトパラメータγの値には、適当に大きな値（１．０２〜１．２）が経験的に用いられている。
【先行技術文献】
【特許文献】
【００２２】
【特許文献１】特開２０１２−２４７９７３号公報
【特許文献２】国際公開第２００８／０２６２６１号
【特許文献３】特開２０１２−７３６８１号公報
【発明の概要】
【発明が解決しようとする課題】
【００２３】
磁場解析の連立方程式に対する前処理を安定させるために用いられるシフトパラメータに関し、適切な値を高速に算出することができないという問題がある。対角成分Ｄｉｉは全てのｉに関して正である必要があるが、対角成分Ｄｉｉが大きすぎるとその後の反復計算の収束の反復回数が増大してしまう。このため、対角成分Ｄｉｉが零よりやや大きな値となるシフトパラメータを高速に決めることが重要となる。
【００２４】
本発明は、１つの側面では、磁場解析の連立方程式に対する前処理で用いられるシフトパラメータに関し、適切な値を高速に算出することを目的とする。
【課題を解決するための手段】
【００２５】
第１の案では、磁化解析装置は、磁場の連立方程式の解を反復計算で計算する磁化解析装置において、前記連立方程式を構成する係数行列を分解して前処理行列を取得する取得部と、取得された前記前処理行列の対角成分の値を調整するパラメータを用いて、前記前処理行列の対角成分の値が正になるように、前記パラメータの値を２分探索に基づいて算出する算出部と、を有する。
【発明の効果】
【００２６】
１実施態様によれば、磁場解析の連立方程式に対する前処理で用いられるシフトパラメータに関し、適切な値を高速に算出することができる。
【図面の簡単な説明】
【００２７】
【図１】図１は、実施例１に係る磁化解析装置の構成を示す機能ブロック図である。
【図２】図２は、磁場の連立方程式を説明する図である。
【図３】図３は、前処理行列の対角成分Ｄｉｉとシフトパラメータγとの関係を示す図である。
【図４】図４は、回路無しの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。
【図５】図５は、回路有りの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。
【図６】図６は、実施例１に係るシフトパラメータ計算処理を示すフローチャートである。
【図７Ａ】図７Ａは、調整パラメータηと閾値εを独立した場合の反復回数の一例を示す図である。
【図７Ｂ】図７Ｂは、調整パラメータηと閾値εを連動させた場合の反復回数の一例を示す図である。
【図８】図８は、実施例２に係る磁化解析装置の構成を示す機能ブロック図である。
【図９】図９は、並列計算する場合のＣＰＵの分担の一例を示す図である。
【図１０】図１０は、磁化解析プログラムを実行するコンピュータの一例を示す図である。
【図１１】図１１は、ＣＧ法の処理の一例を示す図である。
【図１２】図１２は、前処理付きのＣＧ（ＩＣＣＧ）法の処理の一例を示す図である。
【図１３】図１３は、ＣＲ法の処理の一例を示す図である。
【図１４】図１４は、前処理付きのＣＲ（ＩＬＵ  ＣＲ）法の処理の一例を示す図である。
【発明を実施するための形態】
【００２８】
以下に、本願の開示する磁化解析装置、磁化解析方法および磁化解析プログラムの実施例を図面に基づいて詳細に説明する。なお、磁化解析装置は、磁場解析の過程で生成される連立方程式の解を反復計算で計算する装置であり、前処理付きのＣＧ法（ＩＣＣＧ法）や前処理付きのＣＲ法（ＩＬＵ  ＣＲ法）に適用した場合を示す。また、本発明は、実施例により限定されるものではなく、磁化解析に広く適用可能である。
【実施例１】
【００２９】
図１は、実施例１に係る磁化解析装置の構成を示す機能ブロック図である。図１に示すように、磁化解析装置１は、入力部２と、出力部３と、記憶部４と、制御部５とを有する。
【００３０】
入力部２は、解析を行うユーザが各種の情報や指示を磁化解析装置１に入力するための入力装置である。例えば、入力部２は、キーボード、マウス、タッチパネルに対応する。出力部３は、各種の情報を表示する出力装置である。例えば、出力部３は、ディスプレイ、タッチパネルに対応する。
【００３１】
記憶部４は、例えば、ＲＡＭ、フラッシュメモリ（Flash  Memory）などの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。記憶部４は、計算条件データ４１と、シフトパラメータデータ４２と、結果データ４３とを記憶する。
【００３２】
計算条件データ４１は、磁化解析の計算条件に関するデータである。計算条件データ４１には、例えば、前処理行列を生成するために用いられる係数行列、シフトパラメータを計算するために用いられる各種の定数が含まれる。
【００３３】
ここで、係数行列とは、有限要素法による磁場解析の過程で生成される連立方程式に含まれる疎行列のことをいう。前処理行列は、係数行列の近似である。係数行列は、対称行列や非対称行列を含む。シフトパラメータとは、前処理行列の対角成分の値を調整するパラメータのことをいう。すなわち、ＩＣＣＧ法では、シフトパラメータは、式（２）の前処理行列の対角成分Ｄｉｉが零よりやや大きな値となるように調整するパラメータのことをいう。ＩＬＵ  ＣＲ法では、シフトパラメータは、式（５）の前処理行列の対角成分Ｄｉｉが零よりやや大きな値となるように調整するパラメータのことをいう。なお、係数行列は、適宜、「Ａ」の符号で表わす。前処理行列は、適宜、「Ｃ」の符号で表わす。シフトパラメータは、適宜、「γ」の符号で表わす。
【００３４】
シフトパラメータデータ４２は、シフトパラメータに関するデータである。
【００３５】
結果データ４３は、磁場解析の過程で生成される連立方程式の解を反復計算で計算した結果を示すデータである。結果データ４３には、磁場の連立方程式の解が含まれる。
【００３６】
ここで、磁場の連立方程式について、図２を参照して説明する。図２は、磁場の連立方程式を説明する図である。図２に示すように、連立方程式ｒ０は、有限要素法による磁場解析の過程で生成される磁場の連立方程式ｒ１と、回路の連立方程式ｒ２とから構成される。係数行列Ａの四角の黒塗りの領域は、値が零でないことを示す。係数行列Ａの四角の白塗りの領域および何もない領域は、値が零であることを示す。連立方程式ｒ０の係数行列Ａの上部分は、磁場の連立方程式ｒ１の未知数であるベクトルポテンシャルやスカラーポテンシャルの領域（ｉ∈Ｘ）となるように並べ替えた状態である。連立方程式ｒ０の係数行列Ａの下部分は、回路の連立方程式ｒ２の未知数である電流の領域（ｉ∈Ｙ）となるように並べ替えた状態である。未知数が連立方程式の解となる。ここでいうｉは、全自由度の中の特定の行番号である。ここでいうＸは、ベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号の集合である。ここでいうＹは、電流を未知数とする行番号の集合である。
【００３７】
シフトパラメータγは、係数行列Ａの対角成分Ａｉｉを調整する（式（８）、式（９）参照）。連立方程式ｒ０に回路の連立方程式ｒ２を含む場合には、連立方程式ｒ０に含まれる係数行列Ａの電流を未知数とする行の対角成分Ａｉｉ（ｂ０）は零である。このような場合には、どんなに大きなシフトパラメータγを係数行列Ａの対角成分に乗じても前処理行列Ｃの対角成分Ｄｉｉは正にならない。これは、ＩＣＣＧ法の場合は式（８）、ＩＬＵ  ＣＲ法の場合は式（９）によると、前処理行列Ｃの対角成分Ｄｉｉは、負または零になってしまうからである。そこで、連立方程式ｒ０に回路の連立方程式ｒ２を含む場合には、電流を未知数とする行を除外して、残りのＤｉｉが正となるようにシフトパラメータγを計算する。
【００３８】
図１に戻って、制御部５は、ＣＰＵ（Central  Processing  Unit）などの電子回路に対応する。そして、制御部５は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。例えば、制御部５は、磁化解析処理を実行する。磁化解析処理は、記憶部４から計算条件データ４１を読み込んで計算を開始する。そして、磁化解析処理は、計算の過程で、計算条件データ４１の係数行列を用いて、シフトパラメータを計算し、計算したシフトパラメータを、前処理行列を調整するデータ（シフトパラメータデータ４２）として記憶部４に保存する。そして、磁化解析処理は、シフトパラメータを用いて、磁場の連立方程式の解を反復計算で求める。
【００３９】
制御部５は、シフトパラメータ計算部６と、前処理用行列計算部７と、反復解法計算部８とを有する。シフトパラメータ計算部６は、計算条件データ４の係数行列に近似する前処理行列の対角成分が零よりやや大きい正の値となるように調整すべく、シフトパラメータを計算する。シフトパラメータ計算部６は、計算されたシフトパラメータをシフトパラメータデータ４２として記憶部４に格納する。
【００４０】
シフトパラメータ計算部６は、回路無し計算部６１および回路有り計算部６２を有する。
【００４１】
回路無し計算部６１は、連立方程式ｒ０に回路の連立方程式ｒ２を含まない場合のシフトパラメータを計算する。例えば、回路無し計算部６１は、磁場の連立方程式を構成する係数行列を分解して前処理行列を取得する。この前処理行列を取得する形態には、磁場の連立方程式を構成する係数行列を計算することにより取得することも含まれる。係数行列の分解は、一例として、ＩＣＣＧ法の場合には、不完全コレスキー分解を指し、ＩＬＵ  ＣＲ法の場合には、不完全ＬＵ分解を指す。回路無し計算部６１は、前処理行列の対角成分が全て正になるように、シフトパラメータを２分探索法に基づいて算出する。
【００４２】
ここで、前処理行列の対角成分とシフトパラメータとの関係を、図３を参照して説明する。前処理行列の対角成分Ｄｉｉは、例えば、式（８）または式（９）で表わされ、シフトパラメータγに対して単調増加の関数で求められる。この関数をシフトパラメータγの関数としてＤｉｉ（γ）と表記するものとする。
【００４３】
図３は、前処理行列の対角成分Ｄｉｉとシフトパラメータγとの関係を示す図である。なお、Ｘ軸は、シフトパラメータγの値を示し、Ｙ軸は、対角成分Ｄｉｉの値を示す。図３に示すように、シフトパラメータγが増加するにつれて、対角成分Ｄｉｉは単調に増加する。
【００４４】
図１に戻って、一例として、回路無し計算部６１は、対角成分Ｄｉｉが負となる可能性のあるシフトパラメータγの値を最小値γｍｉｎとして設定する。例えば、図３が示す関数の場合であれば、シフトパラメータγが１のとき、Ｄｉｉが負になるので、最小値γｍｉｎを１．０としても良い。回路無し計算部６１は、対角成分Ｄｉｉが必ず正となるシフトパラメータγの値を最大値γｍａｘとして設定する。すなわち、Ｄｉｉ（γｍｉｎ）が負であって且つＤｉｉ（γｍａｘ）が正であれば、γｍｉｎとγｍａｘの間にＤｉｉが０となるシフトパラメータγが存在する。そこで、回路無し計算部６１は、Ｄｉｉが０に近似するシフトパラメータγを２分探索に基づいて算出する。回路無し計算部６１は、２分探索を、γｍｉｎおよびγｍａｘの差分が閾値εより小さくなったときに終了する。
【００４５】
ここで、閾値εは、シフトパラメータを調整するパラメータである調整パラメータηに対する線形の式ε＝kη（０＜ｋ＜１）に基づいて定められる。すなわち、ｋは、ε＜ηとなるように定められる。一例として、η、εの標準値は、それぞれ０．０２、１．０×１０−３と定めることができる。回路無し計算部６１は、それぞれの標準値を調整することで、さらに、計算を高速化することができる。
【００４６】
回路無し計算部６１は、２分探索が終了すると、調整パラメータη、γｍｉｎおよびγｍａｘを用いて、式（１０）に基づいて、シフトパラメータγを算出する。すなわち、回路無し計算部６１は、γｍｉｎとγｍａｘとの間のシフトパラメータに対して、１．０に調整パラメータηを加算して得られた値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する。なお、ｔは、０から１．０までの間の数値であり、例えば０．５であるが、これに限定しない。
【数１０】
0006337682.tif 000011
【００４７】
ここで、シフトパラメータγの算出式を式（１０）としたのは、以下のとおりである。閾値の線形の式ε＝ｋηのｋは、０＜ｋ＜１内のいずれかの値であるとする。調整パラメータηが規定値より大きい場合、閾値の線形の式ε＝ｋηにより閾値εはより大きい値となる。閾値εがより大きい値であれば、２分探索の回数が少なくなるので、式（１０）の右辺の第２式（ｔγｍｉｎ＋（１．０−ｔ）γｍａｘ）の値が目的とする値より小さくなる。そして、右辺の第１式（１．０＋η）の値が調整パラメータηによってより大きい値となるので、第２式の値が目的とする値より小さくなっても最終的にはシフトパラメータγの値は、小さい値とならず、零よりやや大きい正の値になり得る。一方、調整パラメータηが規定値より小さい場合、閾値の線形の式ε＝ｋηにより閾値εはより小さい値となる。閾値εがより小さい値であれば、第２式の値が目的とする値に近づく値となるので、第１式の値がより小さい値となっても、最終的にはシフトパラメータγの値は、小さい値とならず、零よりやや大きい正の値になり得る。したがって、シフトパラメータγの算出式を式（１０）とすることで、適切なシフトパラメータを算出することができ、その後の磁場の連立方程式の反復計算の収束性を保つことができる。
【００４８】
ここで、連立方程式ｒ０に回路の連立方程式ｒ２を含まない場合に、シフトパラメータγを計算する際の２分探索のアルゴリズムについて説明する。図４は、回路無しの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。なお、回路無し計算部６１は、最初に、Ｄｉｉが負となる可能性のあるシフトパラメータγの値を最小値γｍｉｎとして決定し、Ｄｉｉが必ず正となるシフトパラメータγの値を最大値γｍａｘとして決定したものとする。
【００４９】
図４に示すように、ステップ１（初期化処理）として、回路無し計算部６１は、調整パラメータηに初期値η０を設定し、閾値εに初期値ε０を設定する（Ｓ１００）。
【００５０】
ステップ２として、回路無し計算部６１は、現に設定された最小値γｍｉｎおよび現に設定された最大値γｍａｘの中間値（γｍｉｎ＋γｍａｘ）／２に対応するＤｉｉの最小値を算出する。回路無し計算部６１は、算出した最小値Ｄｉｉが零より大きければ（Ｓ１０１）、以下の処理を行う。すなわち、回路無し計算部６１は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの差分ｄγを算出する。回路無し計算部６１は、現に設定された最小値γｍｉｎおよび現に設定された最大値γｍａｘの中間値（γｍｉｎ＋γｍａｘ）／２を新たな最大値γｍａｘに設定する（Ｓ１０２）。すなわち、回路無し計算部６１は、シフトパラメータγの探索範囲を半分に狭めるべく、最大値γｍａｘをＤｉｉの値が小さくなる方向に設定する。つまり、Ｄｉｉの値が正であって零に近い方向に設定する。
【００５１】
一方、回路無し計算部６１は、算出した最小値Ｄｉｉが零以下であれば（Ｓ１０３）、以下の処理を行う。すなわち、回路無し計算部６１は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの差分ｄγを算出する。回路無し計算部６１は、現に設定された最小値γｍｉｎおよび現に設定された最大値γｍａｘの中間値（γｍｉｎ＋γｍａｘ）／２を新たな最小値γｍｉｎに設定する（Ｓ１０４）。すなわち、回路無し計算部６１は、シフトパラメータγの探索範囲を半分に狭めるべく、最小値γｍｉｎをＤｉｉの値が大きくなる方向に設定する。
【００５２】
そして、回路無し計算部６１は、差分ｄγが閾値εより小さければ（Ｓ１０５）、式（１０）に基づいて、シフトパラメータγを算出する。ここでは、式（１０）のｔを０．５としている。すなわち、回路無し計算部６１は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの中間値に、１より大きな値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する（Ｓ１０６）。
【００５３】
一方、回路無し計算部６１は、差分ｄγが閾値ε以上であれば（Ｓ１０７）、現に設定された最大値γｍａｘおよび最小値γｍｉｎを用いて２分探索を繰り返すべく、Ｓ１０１に遷移する。
【００５４】
図１に戻って、回路有り計算部６２は、連立方程式ｒ０に回路の連立方程式ｒ２を含む場合のシフトパラメータを計算する。例えば、回路有り計算部６２は、磁場の連立方程式を構成する係数行列を分解して前処理行列を取得する。この前処理行列を取得する形態には、磁場の連立方程式を構成する係数行列を計算することにより取得することも含まれる。係数行列の分解は、回路無し計算部６１におけるものと同じであるので、その説明を省略する。回路有り計算部６２は、前処理行列の対角成分のうち回路の行の対角成分を除外し、残った対角成分が全て正になるように、シフトパラメータを２分探索法に基づいて算出する。これは、回路の行の対角成分は零であるので（図２のｂ０参照）、ＩＣＣＧ法の場合は式（８）、ＩＬＵ  ＣＲ法の場合は式（９）によると、前処理行列Ｃの対角成分Ｄｉｉは、負または零になってしまうからである。
【００５５】
ここで、連立方程式ｒ０に回路の連立方程式ｒ２を含む場合に、シフトパラメータγを計算する際の２分探索のアルゴリズムについて説明する。図５は、回路有りの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。なお、連立方程式ｒ０に回路の連立方程式ｒ２を含まない場合の２分探索のアルゴリズムと同じ処理は、同じ符号で表わすものとする。また、回路有り計算部６２は、最初に、Ｄｉｉが負となる可能性のあるシフトパラメータγの値を最小値γｍｉｎとして決定し、Ｄｉｉが必ず正となるシフトパラメータγの値を最大値γｍａｘとして決定したものとする。
【００５６】
図５に示すように、ステップ１（初期化処理）として、回路有り計算部６２は、調整パラメータとしてのηに初期値η０を設定し、調整パラメータとしての閾値εに初期値ε０を設定する（Ｓ１００）。
【００５７】
ステップ２として、回路有り計算部６２は、現に設定された最小値γｍｉｎおよび現に設定された最大値γｍａｘの中間値（γｍｉｎ＋γｍａｘ）／２に対応するＤｉｉの最小値（Ｍｉｎ）を算出する。ここで、「ｉ∈Ｘ」は、行番号ｉはベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号の集合に属するという意味である。したがって、ベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号ｉのみが用いられ、回路の行、すなわち電流を未知数とする行番号ｉは除外される。回路有り計算部６２は、算出した最小値Ｄｉｉが零より大きければ（Ｓ２００）、以下の処理を行う。
【００５８】
すなわち、回路有り計算部６２は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの差分ｄγを算出する。回路有り計算部６２は、現に設定された最小値γｍｉｎおよび現に設定された最大値γｍａｘの中間値（γｍｉｎ＋γｍａｘ）／２を新たな最大値γｍａｘに設定する（Ｓ１０２）。すなわち、回路有り計算部６２は、シフトパラメータγの探索範囲を半分に狭めるべく、最大値γｍａｘをＤｉｉの値が小さくなる方向に設定する。つまり、Ｄｉｉの値が正であって零に近い方向に設定する。
【００５９】
一方、回路有り計算部６２は、算出した最小値Ｄｉｉが零以下であれば（Ｓ２０１）、以下の処理を行う。すなわち、回路有り計算部６２は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの差分ｄγを算出する。回路有り計算部６２は、現に設定された最小値γｍｉｎおよび現に設定された最大値γｍａｘの中間値（γｍｉｎ＋γｍａｘ）／２を新たな最小値γｍｉｎに設定する（Ｓ１０４）。すなわち、回路有り計算部６２は、シフトパラメータγの探索範囲を半分に狭めるべく、最大値γｍｉｎをＤｉｉの値が大きくなる方向に設定する。
【００６０】
そして、回路有り計算部６２は、差分ｄγが閾値εより小さければ（Ｓ１０５）、式（１０）に基づいて、シフトパラメータγを算出する。ここでは、式（１０）のｔを０．５としている。すなわち、回路有り計算部６２は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの中間値に、１より大きな値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する（Ｓ１０６）。
【００６１】
一方、回路有り計算部６２は、差分ｄγが閾値ε以上であれば（Ｓ１０７）、現に設定された最大値γｍａｘおよび最小値γｍｉｎを用いて２分探索を繰り返すべく、Ｓ２００に遷移する。
【００６２】
図１に戻って、前処理用行列計算部７は、シフトパラメータ計算部６によって計算されたシフトパラメータを基に、前処理用行列を計算する。例えば、前処理用行列計算部７は、記憶部４からシフトパラメータデータ４２を読み込む。
【００６３】
一例として、ＩＣＣＧ法の場合には、前処理用行列計算部７は、行列Ｄの対角成分Ｄｉｉを、式（８）を用いて計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄｉｉを式（３）に代入して、行列Ｌの非対角成分Ｌｉｊを計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄｉｉおよび行列Ｌの非対角成分Ｌｉｊを式（１）に代入して、前処理行列Ｃを計算する。
【００６４】
一例として、ＩＬＵ  ＣＲ法の場合には、前処理用行列計算部７は、行列Ｄの対角成分Ｄｉｉを、式（９）を用いて計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄｉｉを式（６）に代入して、非対角成分Ｌijを計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄｉｉを式（７）に代入して、非対角成分Ｕijを計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄｉｉ、非対角成分Ｌｉｊおよび非対角成分Ｕijを式（４）に代入して、前処理行列Ｃを計算する。
【００６５】
反復解法計算部８は、ＩＣＣＧ法またはＩＬＵ  ＣＲ法を用いた反復計算で磁場の連立方程式の解を計算する。
【００６６】
［シフトパラメータ計算処理の手順］
図６は、実施例１に係るシフトパラメータ計算処理を示すフローチャートである。なお、図６では、磁場の連立方程式の解を反復計算で計算する回路有りの場合のシフトパラメータ計算処理について説明する。係数行列Ａを含む磁場の連立方程式の解を求める方法は、ＩＣＣＧ法であり、シフトパラメータγの関数Ｄｉｉ（γ）は、式（８）で表わされる。また、図６では、係数行列Ａの全行数はＮであるとする。
【００６７】
図６に示すように、回路有り計算部６２は、シフトパラメータ計算処理で用いられる変数を初期化する（ステップＳ１１）。例えば、回路有り計算部６２は、計算条件データ４１に基づいて、変数ｋに初期値としてｋ０を設定し、調整パラメータηに初期値としてη０を設定し、調整パラメータとしての閾値εに変数ｋと調整パラメータηとを乗じた値を設定する。回路有り計算部６２は、シフトパラメータγｍｉｎに初期値として１．０を設定する。シフトパラメータγｍｉｎには、対角成分Ｄｉｉが負となる可能性のあるシフトパラメータγの値が設定される。回路有り計算部６２は、シフトパラメータγの初期値としてγｍｉｎの値を設定する。回路有り計算部６２は、シフトパラメータγｍａｘに初期値として２．０を設定する。シフトパラメータγｍａｘには、対角成分Ｄｉｉが必ず正となるシフトパラメータγの値が設定される。回路有り計算部６２は、変数ｔに初期値として０．５を設定する。なお、変数ｔは、０から１．０までの間の数値である。
【００６８】
回路有り計算部６２は、行番号ｉに初期値として１を設定し、対角成分Ｄｉｉの最小値Ｄｍｉｎに初期値として十分に大きな値である１．０Ｅ２０を設定する（ステップＳ１２）。回路有り計算部６２は、行番号ｉがベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号の集合Ｘに属するか否かを判定する（ステップＳ１３）。
【００６９】
回路有り計算部６２は、行番号ｉが行番号の集合Ｘに属しないと判定した場合には（ステップＳ１３；Ｎｏ）、次の行番号ｉに進むべく、ステップＳ１７に移行する。一方、回路有り計算部６２は、行番号ｉが行番号の集合Ｘに属すると判定した場合には（ステップＳ１３；Ｙｅｓ）、式（８）の関数Ｄｉｉ（γ）に基づいて対角成分Ｄｉｉの計算を行う（ステップＳ１４）。
【００７０】
回路有り計算部６２は、計算結果である対角成分ＤｉｉがＤｍｉｎ以下であるか否かを判定する（ステップＳ１５）。回路有り計算部６２は、計算結果である対角成分ＤｉｉがＤｍｉｎ以下であると判定した場合には（ステップＳ１５；Ｙｅｓ）、Ｄｍｉｎの値を計算結果である対角成分Ｄｉｉの値に置き換える（ステップＳ１６）。そして、回路有り計算部６２は、次の行番号ｉに進むべく、ステップＳ１７に移行する。
【００７１】
一方、回路有り計算部６２は、計算結果である対角成分ＤｉｉがＤｍｉｎより大きいと判定した場合には（ステップＳ１５；Ｎｏ）、次の行番号ｉに進むべく、ステップＳ１７に移行する。
【００７２】
ステップＳ１７において、回路有り計算部６２は、行番号ｉを１加算する(ステップＳ１７）。そして、回路有り計算部６２は、行番号ｉが全行数Ｎより大きいか否かを判定する（ステップＳ１８）。回路有り計算部６２は、行番号ｉが全行数Ｎより大きくないと判定した場合には（ステップＳ１８；Ｎｏ）、行番号ｉの対角成分Ｄｉｉにおける処理を行うべく、ステップＳ１３に遷移する。
【００７３】
一方、回路有り計算部６２は、行番号ｉが全行数Ｎより大きいと判定した場合には（ステップＳ１８；Ｙｅｓ）、現在のシフトパラメータγに対応するＤｍｉｎ、すなわちＤｉｉの最小値であるＤｍｉｎが０．０より大きいか否かを判定する（ステップＳ１９）。
【００７４】
回路有り計算部６２は、Ｄｍｉｎが０．０より大きくないと判定した場合には（ステップＳ１９；Ｎｏ）、以下の処理を行う。すなわち、回路有り計算部６２は、現に設定されたγｍａｘからγｍｉｎを減算し、差分ｄγを算出する。回路有り計算部６２は、現に設定されたγｍａｘおよびγｍｉｎの中間値（γｍｉｎ＋γｍａｘ）／２を算出し、新たなγｍｉｎに設定する。回路有り計算部６２は、シフトパラメータγにγｍｉｎの値を設定する（ステップＳ２１）。すなわち、回路有り計算部６２は、シフトパラメータγの値の探索範囲を半分に狭めるべく、γｍｉｎをＤｉｉの値が大きくなる方向に設定する。そして、回路有り計算部６２は、ステップＳ２２に移行する。
【００７５】
一方、回路有り計算部６２は、Ｄｍｉｎが０．０より大きいと判定した場合には（ステップＳ１９；Ｙｅｓ）、以下の処理を行う。すなわち、回路有り計算部６２は、現に設定されたγｍａｘからγｍｉｎを減算し、差分ｄγを算出する。回路有り計算部６２は、現に設定されたγｍａｘおよびγｍｉｎの中間値（γｍｉｎ＋γｍａｘ）／２を算出し、新たなγｍａｘに設定する。回路有り計算部６２は、シフトパラメータγにγｍａｘの値を設定する（ステップＳ２０）。すなわち、回路有り計算部６２は、シフトパラメータγの値の探索範囲を半分に狭めるべく、γｍａｘをＤｉｉの値が小さくなる方向に設定する。つまり、Ｄｉｉの値が正であって零に近い方向に設定する。そして、回路有り計算部６２は、ステップＳ２２に移行する。
【００７６】
ステップＳ２２において、回路有り計算部６２は、γｍａｘおよびγｍｉｎの差分ｄγが閾値εより小さいか否かを判定する（ステップＳ２２）。回路有り計算部６２は、差分ｄγが閾値εより小さくないと判定した場合には（ステップＳ２２；Ｎｏ）、新たに設定されたシフトパラメータγに対する対角成分Ｄｉｉにおける処理を行うべく、ステップＳ１２に遷移する。
【００７７】
一方、回路有り計算部６２は、差分ｄγが閾値εより小さいと判定した場合には（ステップＳ２２；Ｙｅｓ）、式（１０）に基づいて、シフトパラメータγを算出する（ステップＳ２３）。すなわち、回路有り計算部６２は、現に設定された最小値γｍｉｎと現に設定された最大値γｍａｘとの間のいずれかの値に、１より大きな値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する。そして、回路有り計算部６２は、シフトパラメータ計算処理を終了する。
【００７８】
[実施例１の効果]
上記実施例１によれば、磁化解析装置１は、磁場解析の過程で生成される連立方程式を構成する係数行列を分解して前処理行列を取得する。この前処理行列を取得する形態には、磁場の連立方程式を構成する係数行列を計算することにより取得することも含まれる。磁化解析装置１は、取得された前処理行列の対角成分の値を調整するシフトパラメータを用いて前処理行列の対角成分の値が正になるように、シフトパラメータの値を２分探索に基づいて算出する。かかる構成によれば、磁化解析装置１は、前処理行列の対角成分の値が正になるように、２分探索に基づいてシフトパラメータの値を算出するので、シフトパラメータの適切な値を高速に計算できる。この結果、磁化解析装置１は、磁場解析の過程で生成される連立方程式を安定且つ高速に計算できる。
【００７９】
また、上記実施例１によれば、磁化解析装置１は、連立方程式が回路に係る連立方程式を含む場合、回路に係る連立方程式に関係する、係数行列の行の対角成分を除外してシフトパラメータの値を算出する。かかる構成によれば、磁化解析装置１は、回路に関係する行の対角成分は零であるので、零であるとわかっている行の対角成分を予め除外してパラメータの値を算出するので、適切な値を算出できる。
【００８０】
また、上記実施例１によれば、磁化解析装置１は、シフトパラメータの値の最小値および最大値を予め設定する。磁化解析装置１は、現に設定されたシフトパラメータの最小値および現に設定されたシフトパラメータの最大値の中間値に対応する、前処理行列の対角成分の最小値を算出する。磁化解析装置１は、現に設定されたシフトパラメータの最小値と現に設定されたシフトパラメータの最大値との差分を算出する。磁化解析装置１は、前処理行列の対角成分の最小値が零より大きい場合に、現に設定されたシフトパラメータの最小値および現に設定されたシフトパラメータの最大値の中間値を新たなシフトパラメータの最大値に設定し、零以下である場合に、現に設定されたシフトパラメータの最小値および現に設定されたシフトパラメータの最大値との中間値を新たなシフトパラメータの最小値に設定する。磁化解析装置１は、算出された差分が閾値より小さくなるまで、シフトパラメータの算出処理を繰り返す。かかる構成によれば、磁化解析装置１は、シフトパラメータの探索範囲を半分に狭めながら、前処理行列の対角成分の値が正になるようにシフトパラメータの値を算出するので、シフトパラメータの適切な値を高速に計算できる。
【００８１】
また、上記実施例１によれば、磁化解析装置１は、算出された差分が閾値より小さい場合に、現に設定されたシフトパラメータの最小値と最大値との間の値に、１に調整パラメータηの値を加算して得られた値を乗じて得た値をシフトパラメータの値として算出する。かかる構成によれば、磁化解析装置１は、調整パラメータηに基づいて、シフトパラメータの値を微調整することができる。
【００８２】
また、上記実施例１によれば、磁化解析装置１は、閾値εを、調整パラメータとして予め定められた一変数ηの線形式ε＝kη（０＜ｋ＜１）を用いて算出する。かかる構成によれば、磁化解析装置１は、シフトパラメータの値を微調整する調整パラメータηを用いて閾値εを算出するので、より適切な値のシフトパラメータを算出できる。この結果、磁化解析装置１は、少ない反復回数で安定して磁場の連立方程式の解を求めることができる。
【００８３】
ここで、磁化解析装置１が、少ない反復回数で安定して磁場の連立方程式の解を求めることができるという効果について、図７Ａおよび図７Ｂを参照して説明する。図７Ａは、調整パラメータηと閾値εを独立に設定した場合の反復回数を示す図である。図７Ｂは、調整パラメータηと閾値εを連動させて設定した場合の反復回数を示す図である。なお、図７Ａおよび図７Ｂにおいて、「収束せず」とは、該当するシフトパラメータでは、磁場の連立方程式の反復計算が収束しないということを意味する。
【００８４】
調整パラメータηと閾値εを独立に設定した場合の反復回数は、以下のようになる。図７Ａに示すように、閾値εの設定によっては、収束しない場合が多くなる。ここでは、閾値εが０．１である場合には、調整パラメータηが０．０１、０．０２、０．０５、０．１である場合に収束しない。また、閾値εが０．００１の場合には、収束するが、調整パラメータηの値を増やすと、反復回数の増加が大きくなってしまう。
【００８５】
これに対して、調整パラメータηと閾値εを連動させて設定した場合の反復回数は、以下のようになる。図７Ｂに示すように、閾値εを調整パラメータηとｋの線形式から定めると、ｋが１．０である場合を除き、必ず収束している。少ない反復回数を１０未満とすると、ｋが０．５〜０．８とすることで、少ない反復回数で安定に磁場の連立方程式の解を求めることができる。
【実施例２】
【００８６】
ところで、実施例１に係る磁化解析装置１は、並列計算を用いないで、係数行列からシフトパラメータを算出する場合を説明した。しかしながら、磁化解析装置１は、並列計算を用いて、係数行列からシフトパラメータを算出しても良い。
【００８７】
そこで、実施例２では、磁化解析装置１は、並列計算を用いて、係数行列からシフトパラメータを算出する場合について説明する。
【００８８】
図８は、実施例２に係る磁化解析装置の構成を示す機能ブロック図である。なお、図１に示す磁化解析装置１と同一の構成については同一符号を示すことで、その重複する構成および動作の説明については省略する。実施例１と実施例２とが異なるところは、制御部５のシフトパラメータ計算部６Ａを変更した点にある。実施例１と実施例２とが異なるところは、制御部５にシフトパラメータ共有部９を追加した点にある。
【００８９】
シフトパラメータ計算部６Ａは、計算条件データ４１の係数行列Ａの中の、並列計算するＣＰＵ（Central  Processing  Unit）に対応させて分担する分担領域について、対角成分が零よりやや大きい正の値となるように調整すべく、シフトパラメータを計算する。シフトパラメータは、並列計算する各ＣＰＵによって、各担当領域について算出される。
【００９０】
ここで、係数行列Ａの中のＣＰＵの分担の一例を、図９を参照して説明する。図９は、並列計算する場合のＣＰＵの分担の一例を示す図である。図９では、２つのＣＰＵが２並列で並列計算する場合を示す。図９に示すように、係数行列Ａの四角の黒塗りの領域は、値が零でないことを示す。係数行列Ａの四角の白塗りの領域および何もない領域は、値が零であることを示す。２つのＣＰＵが２並列で並列計算する場合には、一般に、係数行列Ａは、零でない値の数がほぼ等しくなるように分割される。ＣＰＵ１およびＣＰＵ２は、それぞれ四角の黒塗りのある行を処理する。ここでは、ＣＰＵ１は、領域ａ１、領域ａ３を分担する。ＣＰＵ２は、領域ａ２、領域ａ４を担当する。
【００９１】
ＣＰＵ１およびＣＰＵ２は、それぞれ担当する領域について、それぞれシフトパラメータγを計算する。シフトパラメータの計算方法は、実施例１と同様に、電流を未知数とする行を除外して、残りのＤｉｉが正となるようにシフトパラメータγを計算する。
【００９２】
シフトパラメータ共有部９は、シフトパラメータ計算部６ＡによってＣＰＵ毎に計算されたシフトパラメータを取得し、取得したシフトパラメータのうち最大のシフトパラメータを共有する。すなわち、それぞれのＣＰＵで計算されるシフトパラメータは、同値とは限らない。そこで、シフトパラメータ共有部９は、各ＣＰＵで計算されたシフトパラメータのうち最大のシフトパラメータを共有し、共有するシフトパラメータをシフトパラメータデータ４２として記憶部４に格納する。最大のシフトパラメータの値を共有することで、前処理用行列計算部７および反復解法計算部８は、安定的に磁場の連立方程式の解を求めることができる。
【００９３】
なお、実施例２では、シフトパラメータ計算部６Ａは、２つのＣＰＵが並列計算する場合を説明した。すなわち、シフトパラメータ計算部６Ａは、係数行列Ａの中の２つのＣＰＵがそれぞれ分担する分担領域について、対角成分が零よりやや大きい正の値となるように調整すべく、シフトパラメータを計算する。しかしながら、シフトパラメータ計算部６Ａは、これに限定されず、３つ以上のＣＰＵが並列計算する場合であっても良い。また、シフトパラメータ計算部６Ａは、１つのＣＰＵであっても、１つのプロセスで複数のスレッドが並列計算する場合であっても良い。
【００９４】
［実施例２の効果］
上記実施例２によれば、磁化解析装置１は、前処理行列Ｃの第１の対角成分の値が正になるように、シフトパラメータの値を２分探索に基づいて算出する。磁化解析装置１は、前処理行列Ｃの第２の対角成分の値が正になるように、シフトパラメータの値を２分探索に基づいて算出する。磁化解析装置１は、算出された２つのパラメータの値の最大値を共有する。かかる構成によれば、磁化解析装置１は、前処理行列Ｃの第１の対角成分と第２の対角成分に担当を分けて、それぞれシフトパラメータを算出するので、高速にシフトパラメータを計算できる。また、磁化解析装置１は、算出された２つのパラメータの値の最大値を共有することで、磁場の連立方程式を安定且つ高速に計算できる。
【００９５】
［その他］
なお、図示した磁化解析装置１の各構成要素は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、磁化解析装置１の分散・統合の具体的態様は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。例えば、回路無し計算部６１と回路有り計算部６２とを１つの部として統合しても良い。また、記憶部４を磁化解析装置１の外部装置としてネットワーク経由で接続するようにしても良い。
【００９６】
また、上記実施例で説明した各種の処理は、予め用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータで実行することによって実現することができる。そこで、以下では、図１に示した磁化解析装置１と同様の機能を実現する磁化解析プログラムを実行するコンピュータの一例を説明する。図１０は、磁化解析プログラムを実行するコンピュータの一例を示す図である。
【００９７】
図１０に示すように、コンピュータ２００は、各種演算処理を実行するＣＰＵ２０３と、ユーザからのデータの入力を受け付ける入力装置２１５と、表示装置２０９を制御する表示制御部２０７とを有する。また、コンピュータ２００は、記憶媒体からプログラムなどを読取るドライブ装置２１３と、ネットワークを介して他のコンピュータとの間でデータの授受を行う通信制御部２１７とを有する。また、コンピュータ２００は、各種情報を一時記憶するメモリ２０１と、ＨＤＤ２０５を有する。そして、メモリ２０１、ＣＰＵ２０３、ＨＤＤ２０５、表示制御部２０７、ドライブ装置２１３、入力装置２１５、通信制御部２１７は、バス２１９で接続されている。
【００９８】
ドライブ装置２１３は、例えばリムーバブルディスク２１１用の装置である。ＨＤＤ２０５は、磁化解析プログラム２０５ａおよび磁化解析関連情報２０５ｂを記憶する。
【００９９】
ＣＰＵ２０３は、磁化解析プログラム２０５ａを読み出して、メモリ２０１に展開し、プロセスとして実行する。かかるプロセスは、磁化解析装置１の各機能部に対応する。磁化解析関連情報２０５ｂは、計算条件データ４１、シフトパラメータデータ４２および結果データ４３に対応する。そして、例えばリムーバブルディスク２１１が、磁化解析プログラム２０５ａなどの各情報を記憶する。
【０１００】
なお、磁化解析プログラム２０５ａについては、必ずしも最初からＨＤＤ２０５に記憶させておかなくても良い。例えば、コンピュータ２００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」に当該プログラムを記憶させておく。そして、コンピュータ２００がこれらから磁化解析プログラム２０５ａを読み出して実行するようにしても良い。
【符号の説明】
【０１０１】
１  磁化解析装置
２  入力部
３  出力部
４  記憶部
４１  計算条件データ
４２  シフトパラメータデータ
４３  結果データ
５  制御部
６  シフトパラメータ計算部
６１  回路無し計算部
６２  回路有り計算部
７  前処理用行列計算部
８  反復解法計算部
【図１】
0006337682.tif 000012
【図２】
0006337682.tif 000013
【図３】
0006337682.tif 000014
【図４】
0006337682.tif 000015
【図５】
0006337682.tif 000016
【図６】
0006337682.tif 000017
【図７Ａ】
0006337682.tif 000018
【図７Ｂ】
0006337682.tif 000019
【図８】
0006337682.tif 000020
【図９】
0006337682.tif 000021
【図１０】
0006337682.tif 000022
【図１１】
0006337682.tif 000023
【図１２】
0006337682.tif 000024
【図１３】
0006337682.tif 000025
【図１４】
0006337682.tif 000026