JP6337682B2 - 磁化解析装置、磁化解析方法および磁化解析プログラム - Google Patents

Description

本発明は、磁化解析装置、磁化解析方法および磁化解析プログラムに関する。

コンピュータを用いた数値シミュレーションでは、有限の数の物理量で場を表現するために離散化と呼ばれる技術が用いられる。この離散化手法には、例えば、有限差分法、有限要素法、有限体積法が挙げられる。いずれの離散化手法も空間に有限の数の物理量を配置し、物理現象を支配する方程式を離散化して連立方程式を作成する。そして、コンピュータで連立方程式を解くことで物理量が決定される。

連立方程式を行列とベクトルで表現するとき、この行列は係数行列と呼ばれる。数値シミュレーションで扱う連立方程式は様々であり、係数行列の対角成分に零が存在することがある。また、解が不能であったり、不定であったりすることがある。

辺要素を用いる磁場解析の過程で生成される連立方程式は不定である。辺要素は、例えば、有限要素法で用いられる要素のことである。不定の連立方程式を直接解法の線形ソルバーを適用して解を求めることは困難である。しかし、ＣＧ（conjugate gradient）法またはＢｉＣＧ（Bi-Conjugate Gradient）法のようなクロリフ部分空間を用いる反復法ソルバーであれば、不定の連立方程式の解を求めることが可能である。

ＣＧ法やＢｉＣＧ法に前処理を施す手法は、それぞれＩＣＣＧ法、ＩＬＵ ＢｉＣＧ法と呼ばれる。ＩＣＣＧ（Incomplete Cholesky conjugate gradient）法、ＩＬＵ ＢｉＣＧ（Incomplete LU factorization）法は、前処理を施すことで、収束までに要する反復回数を減らすことが可能となる。

ＣＧ法の処理は、対称行列を反復計算で解を求める。図１１は、ＣＧ法の処理の一例を示す図である。図１１に示すように、ＣＧ法の処理は、前処理をしないで、対称行列を反復計算で解を求めている。図１１では、Ａが対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。

ＩＣＣＧ法の処理は、前処理を施した後に、対称行列を反復計算で解を求める。図１２は、前処理付きのＣＧ（ＩＣＣＧ）法の処理の一例を示す図である。図１２に示すように、ＩＣＣＧ法の処理は、前処理として前処理行列を計算した後に、対称行列を反復計算で解を求めている。図１２では、Ｃが前処理行列であり、Ａが対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。

Ａが対称行列である場合の前処理行列Ｃの計算は、以下の式（１）で表わされる。対称行列Ａを不完全コレスキー分解した行列Ｃは、対称行列Ａの近似である。

式（１）によると、前処理行列Ｃは、非対角成分の行列Ｌと対角成分の行列Ｄと非対角成分の行列Ｌの転置行列Ｌ^Ｔにより表わされる。

行列Ｄの対角成分は、以下の式（２）で表わされる。

行列Ｌの非対角成分は、以下の式（３）で表わされる。

また、非対称行列を反復計算で解を求める手法には、ＢｉＣＧ法、ＢｉＣＧ Ｓｔａｂ法、ＣＧＳ法やＣＲ法等が挙げられる。

ＣＲ法の処理は、非対称行列を反復処理で解を求める。図１３は、ＣＲ法の処理の一例を示す図である。図１３に示すように、ＣＲ法の処理は、前処理をしないで、非対称行列を反復計算で解を求めている。図１３では、Ａが非対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。

ＩＬＵ ＣＲ法の処理は、前処理を施した後に、非対称行列を反復計算で解を求める。図１４は、前処理付きのＣＲ（ＩＬＵ ＣＲ）法の処理の一例を示す図である。図１４に示すように、ＩＬＵ ＣＲ法の処理は、前処理として前処理行列を計算した後に、非対称行列を反復計算で解を求めている。図１４では、Ｃが前処理行列であり、Ａが非対称行列であり、Ａ以外はベクトルである。そして、ｘが解である。

Ａが非対称行列である場合の前処理行列Ｃの計算は、以下の式（４）で表わされる。非対称行列Ａを不完全ＬＵ分解した行列Ｃは、非対称行列Ａの近似である。

式（４）によると、前処理行列Ｃは、非対角成分の行列Ｌと対角成分の行列Ｄと非対角成分の行列Ｕにより表わされる。

行列Ｄの対角成分は、以下の式（５）で表わされる。

行列Ｌの非対角成分は、以下の式（６）で表わされる。

行列Ｕの非対角成分は、以下の式（７）で表わされる。

ところで、辺要素を用いる磁場解析の過程で生成される連立方程式の係数行列Ａを用いて前処理行列Ｃの対角成分Ｄ_ｉｉを計算すると、零または負になる成分が発生する。零または負になる成分が発生すると、前処理が破綻するため、その後の反復計算が発散してしまう。前処理行列Ｃは、係数行列Ａの近似であるため、何らかのパラメータによって調整することが可能である。

従来では、以下の式（８）や式（９）のように１より大きな係数γによって、全てのｉに関して対角成分Ｄ_ｉｉが零以下にならないように調整する方法が知られている。式（８）は、ＣＧ法の場合の調整する方法である。式（９）は、ＣＲ法の場合の調整する方法である。なお、係数γは、シフトパラメータと呼ばれる。

シフトパラメータγの値は、連立方程式を解くための反復計算の反復回数に影響を与えるため、適切に調整することで計算時間を大幅に削減することができる。シフトパラメータγの値には、適当に大きな値（１．０２〜１．２）が経験的に用いられている。

特開２０１２−２４７９７３号公報 国際公開第２００８／０２６２６１号 特開２０１２−７３６８１号公報

磁場解析の連立方程式に対する前処理を安定させるために用いられるシフトパラメータに関し、適切な値を高速に算出することができないという問題がある。対角成分Ｄ_ｉｉは全てのｉに関して正である必要があるが、対角成分Ｄ_ｉｉが大きすぎるとその後の反復計算の収束の反復回数が増大してしまう。このため、対角成分Ｄ_ｉｉが零よりやや大きな値となるシフトパラメータを高速に決めることが重要となる。

本発明は、１つの側面では、磁場解析の連立方程式に対する前処理で用いられるシフトパラメータに関し、適切な値を高速に算出することを目的とする。

第１の案では、磁化解析装置は、磁場の連立方程式の解を反復計算で計算する磁化解析装置において、前記連立方程式を構成する係数行列を分解して前処理行列を取得する取得部と、取得された前記前処理行列の対角成分の値を調整するパラメータを用いて、前記前処理行列の対角成分の値が正になるように、前記パラメータの値を２分探索に基づいて算出する算出部と、を有する。

１実施態様によれば、磁場解析の連立方程式に対する前処理で用いられるシフトパラメータに関し、適切な値を高速に算出することができる。

図１は、実施例１に係る磁化解析装置の構成を示す機能ブロック図である。 図２は、磁場の連立方程式を説明する図である。 図３は、前処理行列の対角成分Ｄ_ｉｉとシフトパラメータγとの関係を示す図である。 図４は、回路無しの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。 図５は、回路有りの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。 図６は、実施例１に係るシフトパラメータ計算処理を示すフローチャートである。 図７Ａは、調整パラメータηと閾値εを独立した場合の反復回数の一例を示す図である。 図７Ｂは、調整パラメータηと閾値εを連動させた場合の反復回数の一例を示す図である。 図８は、実施例２に係る磁化解析装置の構成を示す機能ブロック図である。 図９は、並列計算する場合のＣＰＵの分担の一例を示す図である。 図１０は、磁化解析プログラムを実行するコンピュータの一例を示す図である。 図１１は、ＣＧ法の処理の一例を示す図である。 図１２は、前処理付きのＣＧ（ＩＣＣＧ）法の処理の一例を示す図である。 図１３は、ＣＲ法の処理の一例を示す図である。 図１４は、前処理付きのＣＲ（ＩＬＵ ＣＲ）法の処理の一例を示す図である。

以下に、本願の開示する磁化解析装置、磁化解析方法および磁化解析プログラムの実施例を図面に基づいて詳細に説明する。なお、磁化解析装置は、磁場解析の過程で生成される連立方程式の解を反復計算で計算する装置であり、前処理付きのＣＧ法（ＩＣＣＧ法）や前処理付きのＣＲ法（ＩＬＵ ＣＲ法）に適用した場合を示す。また、本発明は、実施例により限定されるものではなく、磁化解析に広く適用可能である。

図１は、実施例１に係る磁化解析装置の構成を示す機能ブロック図である。図１に示すように、磁化解析装置１は、入力部２と、出力部３と、記憶部４と、制御部５とを有する。

入力部２は、解析を行うユーザが各種の情報や指示を磁化解析装置１に入力するための入力装置である。例えば、入力部２は、キーボード、マウス、タッチパネルに対応する。出力部３は、各種の情報を表示する出力装置である。例えば、出力部３は、ディスプレイ、タッチパネルに対応する。

記憶部４は、例えば、ＲＡＭ、フラッシュメモリ（Flash Memory）などの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。記憶部４は、計算条件データ４１と、シフトパラメータデータ４２と、結果データ４３とを記憶する。

計算条件データ４１は、磁化解析の計算条件に関するデータである。計算条件データ４１には、例えば、前処理行列を生成するために用いられる係数行列、シフトパラメータを計算するために用いられる各種の定数が含まれる。

ここで、係数行列とは、有限要素法による磁場解析の過程で生成される連立方程式に含まれる疎行列のことをいう。前処理行列は、係数行列の近似である。係数行列は、対称行列や非対称行列を含む。シフトパラメータとは、前処理行列の対角成分の値を調整するパラメータのことをいう。すなわち、ＩＣＣＧ法では、シフトパラメータは、式（２）の前処理行列の対角成分Ｄ_ｉｉが零よりやや大きな値となるように調整するパラメータのことをいう。ＩＬＵ ＣＲ法では、シフトパラメータは、式（５）の前処理行列の対角成分Ｄ_ｉｉが零よりやや大きな値となるように調整するパラメータのことをいう。なお、係数行列は、適宜、「Ａ」の符号で表わす。前処理行列は、適宜、「Ｃ」の符号で表わす。シフトパラメータは、適宜、「γ」の符号で表わす。

シフトパラメータデータ４２は、シフトパラメータに関するデータである。

結果データ４３は、磁場解析の過程で生成される連立方程式の解を反復計算で計算した結果を示すデータである。結果データ４３には、磁場の連立方程式の解が含まれる。

ここで、磁場の連立方程式について、図２を参照して説明する。図２は、磁場の連立方程式を説明する図である。図２に示すように、連立方程式ｒ０は、有限要素法による磁場解析の過程で生成される磁場の連立方程式ｒ１と、回路の連立方程式ｒ２とから構成される。係数行列Ａの四角の黒塗りの領域は、値が零でないことを示す。係数行列Ａの四角の白塗りの領域および何もない領域は、値が零であることを示す。連立方程式ｒ０の係数行列Ａの上部分は、磁場の連立方程式ｒ１の未知数であるベクトルポテンシャルやスカラーポテンシャルの領域（ｉ∈Ｘ）となるように並べ替えた状態である。連立方程式ｒ０の係数行列Ａの下部分は、回路の連立方程式ｒ２の未知数である電流の領域（ｉ∈Ｙ）となるように並べ替えた状態である。未知数が連立方程式の解となる。ここでいうｉは、全自由度の中の特定の行番号である。ここでいうＸは、ベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号の集合である。ここでいうＹは、電流を未知数とする行番号の集合である。

シフトパラメータγは、係数行列Ａの対角成分Ａ_ｉｉを調整する（式（８）、式（９）参照）。連立方程式ｒ０に回路の連立方程式ｒ２を含む場合には、連立方程式ｒ０に含まれる係数行列Ａの電流を未知数とする行の対角成分Ａ_ｉｉ（ｂ０）は零である。このような場合には、どんなに大きなシフトパラメータγを係数行列Ａの対角成分に乗じても前処理行列Ｃの対角成分Ｄ_ｉｉは正にならない。これは、ＩＣＣＧ法の場合は式（８）、ＩＬＵ ＣＲ法の場合は式（９）によると、前処理行列Ｃの対角成分Ｄ_ｉｉは、負または零になってしまうからである。そこで、連立方程式ｒ０に回路の連立方程式ｒ２を含む場合には、電流を未知数とする行を除外して、残りのＤ_ｉｉが正となるようにシフトパラメータγを計算する。

図１に戻って、制御部５は、ＣＰＵ（Central Processing Unit）などの電子回路に対応する。そして、制御部５は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。例えば、制御部５は、磁化解析処理を実行する。磁化解析処理は、記憶部４から計算条件データ４１を読み込んで計算を開始する。そして、磁化解析処理は、計算の過程で、計算条件データ４１の係数行列を用いて、シフトパラメータを計算し、計算したシフトパラメータを、前処理行列を調整するデータ（シフトパラメータデータ４２）として記憶部４に保存する。そして、磁化解析処理は、シフトパラメータを用いて、磁場の連立方程式の解を反復計算で求める。

制御部５は、シフトパラメータ計算部６と、前処理用行列計算部７と、反復解法計算部８とを有する。シフトパラメータ計算部６は、計算条件データ４の係数行列に近似する前処理行列の対角成分が零よりやや大きい正の値となるように調整すべく、シフトパラメータを計算する。シフトパラメータ計算部６は、計算されたシフトパラメータをシフトパラメータデータ４２として記憶部４に格納する。

シフトパラメータ計算部６は、回路無し計算部６１および回路有り計算部６２を有する。

回路無し計算部６１は、連立方程式ｒ０に回路の連立方程式ｒ２を含まない場合のシフトパラメータを計算する。例えば、回路無し計算部６１は、磁場の連立方程式を構成する係数行列を分解して前処理行列を取得する。この前処理行列を取得する形態には、磁場の連立方程式を構成する係数行列を計算することにより取得することも含まれる。係数行列の分解は、一例として、ＩＣＣＧ法の場合には、不完全コレスキー分解を指し、ＩＬＵ ＣＲ法の場合には、不完全ＬＵ分解を指す。回路無し計算部６１は、前処理行列の対角成分が全て正になるように、シフトパラメータを２分探索法に基づいて算出する。

ここで、前処理行列の対角成分とシフトパラメータとの関係を、図３を参照して説明する。前処理行列の対角成分Ｄ_ｉｉは、例えば、式（８）または式（９）で表わされ、シフトパラメータγに対して単調増加の関数で求められる。この関数をシフトパラメータγの関数としてＤ_ｉｉ（γ）と表記するものとする。

図３は、前処理行列の対角成分Ｄ_ｉｉとシフトパラメータγとの関係を示す図である。なお、Ｘ軸は、シフトパラメータγの値を示し、Ｙ軸は、対角成分Ｄ_ｉｉの値を示す。図３に示すように、シフトパラメータγが増加するにつれて、対角成分Ｄ_ｉｉは単調に増加する。

図１に戻って、一例として、回路無し計算部６１は、対角成分Ｄ_ｉｉが負となる可能性のあるシフトパラメータγの値を最小値γ_ｍｉｎとして設定する。例えば、図３が示す関数の場合であれば、シフトパラメータγが１のとき、Ｄ_ｉｉが負になるので、最小値γ_ｍｉｎを１．０としても良い。回路無し計算部６１は、対角成分Ｄ_ｉｉが必ず正となるシフトパラメータγの値を最大値γ_ｍａｘとして設定する。すなわち、Ｄ_ｉｉ（γ_ｍｉｎ）が負であって且つＤ_ｉｉ（γ_ｍａｘ）が正であれば、γ_ｍｉｎとγ_ｍａｘの間にＤ_ｉｉが０となるシフトパラメータγが存在する。そこで、回路無し計算部６１は、Ｄ_ｉｉが０に近似するシフトパラメータγを２分探索に基づいて算出する。回路無し計算部６１は、２分探索を、γ_ｍｉｎおよびγ_ｍａｘの差分が閾値εより小さくなったときに終了する。

ここで、閾値εは、シフトパラメータを調整するパラメータである調整パラメータηに対する線形の式ε＝kη（０＜ｋ＜１）に基づいて定められる。すなわち、ｋは、ε＜ηとなるように定められる。一例として、η、εの標準値は、それぞれ０．０２、１．０×１０^−３と定めることができる。回路無し計算部６１は、それぞれの標準値を調整することで、さらに、計算を高速化することができる。

回路無し計算部６１は、２分探索が終了すると、調整パラメータη、γ_ｍｉｎおよびγ_ｍａｘを用いて、式（１０）に基づいて、シフトパラメータγを算出する。すなわち、回路無し計算部６１は、γ_ｍｉｎとγ_ｍａｘとの間のシフトパラメータに対して、１．０に調整パラメータηを加算して得られた値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する。なお、ｔは、０から１．０までの間の数値であり、例えば０．５であるが、これに限定しない。

ここで、シフトパラメータγの算出式を式（１０）としたのは、以下のとおりである。閾値の線形の式ε＝ｋηのｋは、０＜ｋ＜１内のいずれかの値であるとする。調整パラメータηが規定値より大きい場合、閾値の線形の式ε＝ｋηにより閾値εはより大きい値となる。閾値εがより大きい値であれば、２分探索の回数が少なくなるので、式（１０）の右辺の第２式（ｔγ_ｍｉｎ＋（１．０−ｔ）γ_ｍａｘ）の値が目的とする値より小さくなる。そして、右辺の第１式（１．０＋η）の値が調整パラメータηによってより大きい値となるので、第２式の値が目的とする値より小さくなっても最終的にはシフトパラメータγの値は、小さい値とならず、零よりやや大きい正の値になり得る。一方、調整パラメータηが規定値より小さい場合、閾値の線形の式ε＝ｋηにより閾値εはより小さい値となる。閾値εがより小さい値であれば、第２式の値が目的とする値に近づく値となるので、第１式の値がより小さい値となっても、最終的にはシフトパラメータγの値は、小さい値とならず、零よりやや大きい正の値になり得る。したがって、シフトパラメータγの算出式を式（１０）とすることで、適切なシフトパラメータを算出することができ、その後の磁場の連立方程式の反復計算の収束性を保つことができる。

ここで、連立方程式ｒ０に回路の連立方程式ｒ２を含まない場合に、シフトパラメータγを計算する際の２分探索のアルゴリズムについて説明する。図４は、回路無しの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。なお、回路無し計算部６１は、最初に、Ｄ_ｉｉが負となる可能性のあるシフトパラメータγの値を最小値γ_ｍｉｎとして決定し、Ｄ_ｉｉが必ず正となるシフトパラメータγの値を最大値γ_ｍａｘとして決定したものとする。

図４に示すように、ステップ１（初期化処理）として、回路無し計算部６１は、調整パラメータηに初期値η_０を設定し、閾値εに初期値ε_０を設定する（Ｓ１００）。

ステップ２として、回路無し計算部６１は、現に設定された最小値γ_ｍｉｎおよび現に設定された最大値γ_ｍａｘの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２に対応するＤ_ｉｉの最小値を算出する。回路無し計算部６１は、算出した最小値Ｄ_ｉｉが零より大きければ（Ｓ１０１）、以下の処理を行う。すなわち、回路無し計算部６１は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの差分ｄγを算出する。回路無し計算部６１は、現に設定された最小値γ_ｍｉｎおよび現に設定された最大値γ_ｍａｘの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２を新たな最大値γ_ｍａｘに設定する（Ｓ１０２）。すなわち、回路無し計算部６１は、シフトパラメータγの探索範囲を半分に狭めるべく、最大値γ_ｍａｘをＤ_ｉｉの値が小さくなる方向に設定する。つまり、Ｄ_ｉｉの値が正であって零に近い方向に設定する。

一方、回路無し計算部６１は、算出した最小値Ｄ_ｉｉが零以下であれば（Ｓ１０３）、以下の処理を行う。すなわち、回路無し計算部６１は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの差分ｄγを算出する。回路無し計算部６１は、現に設定された最小値γ_ｍｉｎおよび現に設定された最大値γ_ｍａｘの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２を新たな最小値γ_ｍｉｎに設定する（Ｓ１０４）。すなわち、回路無し計算部６１は、シフトパラメータγの探索範囲を半分に狭めるべく、最小値γ_ｍｉｎをＤ_ｉｉの値が大きくなる方向に設定する。

そして、回路無し計算部６１は、差分ｄγが閾値εより小さければ（Ｓ１０５）、式（１０）に基づいて、シフトパラメータγを算出する。ここでは、式（１０）のｔを０．５としている。すなわち、回路無し計算部６１は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの中間値に、１より大きな値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する（Ｓ１０６）。

一方、回路無し計算部６１は、差分ｄγが閾値ε以上であれば（Ｓ１０７）、現に設定された最大値γ_ｍａｘおよび最小値γ_ｍｉｎを用いて２分探索を繰り返すべく、Ｓ１０１に遷移する。

図１に戻って、回路有り計算部６２は、連立方程式ｒ０に回路の連立方程式ｒ２を含む場合のシフトパラメータを計算する。例えば、回路有り計算部６２は、磁場の連立方程式を構成する係数行列を分解して前処理行列を取得する。この前処理行列を取得する形態には、磁場の連立方程式を構成する係数行列を計算することにより取得することも含まれる。係数行列の分解は、回路無し計算部６１におけるものと同じであるので、その説明を省略する。回路有り計算部６２は、前処理行列の対角成分のうち回路の行の対角成分を除外し、残った対角成分が全て正になるように、シフトパラメータを２分探索法に基づいて算出する。これは、回路の行の対角成分は零であるので（図２のｂ０参照）、ＩＣＣＧ法の場合は式（８）、ＩＬＵ ＣＲ法の場合は式（９）によると、前処理行列Ｃの対角成分Ｄ_ｉｉは、負または零になってしまうからである。

ここで、連立方程式ｒ０に回路の連立方程式ｒ２を含む場合に、シフトパラメータγを計算する際の２分探索のアルゴリズムについて説明する。図５は、回路有りの場合にシフトパラメータを計算する際の２分探索のアルゴリズムを示す図である。なお、連立方程式ｒ０に回路の連立方程式ｒ２を含まない場合の２分探索のアルゴリズムと同じ処理は、同じ符号で表わすものとする。また、回路有り計算部６２は、最初に、Ｄ_ｉｉが負となる可能性のあるシフトパラメータγの値を最小値γ_ｍｉｎとして決定し、Ｄ_ｉｉが必ず正となるシフトパラメータγの値を最大値γ_ｍａｘとして決定したものとする。

図５に示すように、ステップ１（初期化処理）として、回路有り計算部６２は、調整パラメータとしてのηに初期値η_０を設定し、調整パラメータとしての閾値εに初期値ε_０を設定する（Ｓ１００）。

ステップ２として、回路有り計算部６２は、現に設定された最小値γ_ｍｉｎおよび現に設定された最大値γ_ｍａｘの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２に対応するＤ_ｉｉの最小値（Ｍｉｎ）を算出する。ここで、「ｉ∈Ｘ」は、行番号ｉはベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号の集合に属するという意味である。したがって、ベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号ｉのみが用いられ、回路の行、すなわち電流を未知数とする行番号ｉは除外される。回路有り計算部６２は、算出した最小値Ｄ_ｉｉが零より大きければ（Ｓ２００）、以下の処理を行う。

すなわち、回路有り計算部６２は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの差分ｄγを算出する。回路有り計算部６２は、現に設定された最小値γ_ｍｉｎおよび現に設定された最大値γ_ｍａｘの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２を新たな最大値γ_ｍａｘに設定する（Ｓ１０２）。すなわち、回路有り計算部６２は、シフトパラメータγの探索範囲を半分に狭めるべく、最大値γ_ｍａｘをＤ_ｉｉの値が小さくなる方向に設定する。つまり、Ｄ_ｉｉの値が正であって零に近い方向に設定する。

一方、回路有り計算部６２は、算出した最小値Ｄ_ｉｉが零以下であれば（Ｓ２０１）、以下の処理を行う。すなわち、回路有り計算部６２は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの差分ｄγを算出する。回路有り計算部６２は、現に設定された最小値γ_ｍｉｎおよび現に設定された最大値γ_ｍａｘの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２を新たな最小値γ_ｍｉｎに設定する（Ｓ１０４）。すなわち、回路有り計算部６２は、シフトパラメータγの探索範囲を半分に狭めるべく、最大値γ_ｍｉｎをＤ_ｉｉの値が大きくなる方向に設定する。

そして、回路有り計算部６２は、差分ｄγが閾値εより小さければ（Ｓ１０５）、式（１０）に基づいて、シフトパラメータγを算出する。ここでは、式（１０）のｔを０．５としている。すなわち、回路有り計算部６２は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの中間値に、１より大きな値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する（Ｓ１０６）。

一方、回路有り計算部６２は、差分ｄγが閾値ε以上であれば（Ｓ１０７）、現に設定された最大値γ_ｍａｘおよび最小値γ_ｍｉｎを用いて２分探索を繰り返すべく、Ｓ２００に遷移する。

図１に戻って、前処理用行列計算部７は、シフトパラメータ計算部６によって計算されたシフトパラメータを基に、前処理用行列を計算する。例えば、前処理用行列計算部７は、記憶部４からシフトパラメータデータ４２を読み込む。

一例として、ＩＣＣＧ法の場合には、前処理用行列計算部７は、行列Ｄの対角成分Ｄ_ｉｉを、式（８）を用いて計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄ_ｉｉを式（３）に代入して、行列Ｌの非対角成分Ｌ_ｉｊを計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄ_ｉｉおよび行列Ｌの非対角成分Ｌ_ｉｊを式（１）に代入して、前処理行列Ｃを計算する。

一例として、ＩＬＵ ＣＲ法の場合には、前処理用行列計算部７は、行列Ｄの対角成分Ｄ_ｉｉを、式（９）を用いて計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄ_ｉｉを式（６）に代入して、非対角成分Ｌ_ijを計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄ_ｉｉを式（７）に代入して、非対角成分Ｕ_ijを計算する。前処理用行列計算部７は、計算された行列Ｄの対角成分Ｄ_ｉｉ、非対角成分Ｌ_ｉｊおよび非対角成分Ｕ_ijを式（４）に代入して、前処理行列Ｃを計算する。

反復解法計算部８は、ＩＣＣＧ法またはＩＬＵ ＣＲ法を用いた反復計算で磁場の連立方程式の解を計算する。

［シフトパラメータ計算処理の手順］
図６は、実施例１に係るシフトパラメータ計算処理を示すフローチャートである。なお、図６では、磁場の連立方程式の解を反復計算で計算する回路有りの場合のシフトパラメータ計算処理について説明する。係数行列Ａを含む磁場の連立方程式の解を求める方法は、ＩＣＣＧ法であり、シフトパラメータγの関数Ｄ_ｉｉ（γ）は、式（８）で表わされる。また、図６では、係数行列Ａの全行数はＮであるとする。

図６に示すように、回路有り計算部６２は、シフトパラメータ計算処理で用いられる変数を初期化する（ステップＳ１１）。例えば、回路有り計算部６２は、計算条件データ４１に基づいて、変数ｋに初期値としてｋ_０を設定し、調整パラメータηに初期値としてη_０を設定し、調整パラメータとしての閾値εに変数ｋと調整パラメータηとを乗じた値を設定する。回路有り計算部６２は、シフトパラメータγ_ｍｉｎに初期値として１．０を設定する。シフトパラメータγ_ｍｉｎには、対角成分Ｄ_ｉｉが負となる可能性のあるシフトパラメータγの値が設定される。回路有り計算部６２は、シフトパラメータγの初期値としてγ_ｍｉｎの値を設定する。回路有り計算部６２は、シフトパラメータγ_ｍａｘに初期値として２．０を設定する。シフトパラメータγ_ｍａｘには、対角成分Ｄ_ｉｉが必ず正となるシフトパラメータγの値が設定される。回路有り計算部６２は、変数ｔに初期値として０．５を設定する。なお、変数ｔは、０から１．０までの間の数値である。

回路有り計算部６２は、行番号ｉに初期値として１を設定し、対角成分Ｄ_ｉｉの最小値Ｄ_ｍｉｎに初期値として十分に大きな値である１．０Ｅ２０を設定する（ステップＳ１２）。回路有り計算部６２は、行番号ｉがベクトルポテンシャルやスカラーポテンシャルを未知数とする行番号の集合Ｘに属するか否かを判定する（ステップＳ１３）。

回路有り計算部６２は、行番号ｉが行番号の集合Ｘに属しないと判定した場合には（ステップＳ１３；Ｎｏ）、次の行番号ｉに進むべく、ステップＳ１７に移行する。一方、回路有り計算部６２は、行番号ｉが行番号の集合Ｘに属すると判定した場合には（ステップＳ１３；Ｙｅｓ）、式（８）の関数Ｄ_ｉｉ（γ）に基づいて対角成分Ｄ_ｉｉの計算を行う（ステップＳ１４）。

回路有り計算部６２は、計算結果である対角成分Ｄ_ｉｉがＤ_ｍｉｎ以下であるか否かを判定する（ステップＳ１５）。回路有り計算部６２は、計算結果である対角成分Ｄ_ｉｉがＤ_ｍｉｎ以下であると判定した場合には（ステップＳ１５；Ｙｅｓ）、Ｄ_ｍｉｎの値を計算結果である対角成分Ｄ_ｉｉの値に置き換える（ステップＳ１６）。そして、回路有り計算部６２は、次の行番号ｉに進むべく、ステップＳ１７に移行する。

一方、回路有り計算部６２は、計算結果である対角成分Ｄ_ｉｉがＤ_ｍｉｎより大きいと判定した場合には（ステップＳ１５；Ｎｏ）、次の行番号ｉに進むべく、ステップＳ１７に移行する。

ステップＳ１７において、回路有り計算部６２は、行番号ｉを１加算する(ステップＳ１７）。そして、回路有り計算部６２は、行番号ｉが全行数Ｎより大きいか否かを判定する（ステップＳ１８）。回路有り計算部６２は、行番号ｉが全行数Ｎより大きくないと判定した場合には（ステップＳ１８；Ｎｏ）、行番号ｉの対角成分Ｄ_ｉｉにおける処理を行うべく、ステップＳ１３に遷移する。

一方、回路有り計算部６２は、行番号ｉが全行数Ｎより大きいと判定した場合には（ステップＳ１８；Ｙｅｓ）、現在のシフトパラメータγに対応するＤ_ｍｉｎ、すなわちＤ_ｉｉの最小値であるＤ_ｍｉｎが０．０より大きいか否かを判定する（ステップＳ１９）。

回路有り計算部６２は、Ｄ_ｍｉｎが０．０より大きくないと判定した場合には（ステップＳ１９；Ｎｏ）、以下の処理を行う。すなわち、回路有り計算部６２は、現に設定されたγ_ｍａｘからγ_ｍｉｎを減算し、差分ｄγを算出する。回路有り計算部６２は、現に設定されたγ_ｍａｘおよびγ_ｍｉｎの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２を算出し、新たなγ_ｍｉｎに設定する。回路有り計算部６２は、シフトパラメータγにγ_ｍｉｎの値を設定する（ステップＳ２１）。すなわち、回路有り計算部６２は、シフトパラメータγの値の探索範囲を半分に狭めるべく、γ_ｍｉｎをＤ_ｉｉの値が大きくなる方向に設定する。そして、回路有り計算部６２は、ステップＳ２２に移行する。

一方、回路有り計算部６２は、Ｄ_ｍｉｎが０．０より大きいと判定した場合には（ステップＳ１９；Ｙｅｓ）、以下の処理を行う。すなわち、回路有り計算部６２は、現に設定されたγ_ｍａｘからγ_ｍｉｎを減算し、差分ｄγを算出する。回路有り計算部６２は、現に設定されたγ_ｍａｘおよびγ_ｍｉｎの中間値（γ_ｍｉｎ＋γ_ｍａｘ）／２を算出し、新たなγ_ｍａｘに設定する。回路有り計算部６２は、シフトパラメータγにγ_ｍａｘの値を設定する（ステップＳ２０）。すなわち、回路有り計算部６２は、シフトパラメータγの値の探索範囲を半分に狭めるべく、γ_ｍａｘをＤ_ｉｉの値が小さくなる方向に設定する。つまり、Ｄ_ｉｉの値が正であって零に近い方向に設定する。そして、回路有り計算部６２は、ステップＳ２２に移行する。

ステップＳ２２において、回路有り計算部６２は、γ_ｍａｘおよびγ_ｍｉｎの差分ｄγが閾値εより小さいか否かを判定する（ステップＳ２２）。回路有り計算部６２は、差分ｄγが閾値εより小さくないと判定した場合には（ステップＳ２２；Ｎｏ）、新たに設定されたシフトパラメータγに対する対角成分Ｄ_ｉｉにおける処理を行うべく、ステップＳ１２に遷移する。

一方、回路有り計算部６２は、差分ｄγが閾値εより小さいと判定した場合には（ステップＳ２２；Ｙｅｓ）、式（１０）に基づいて、シフトパラメータγを算出する（ステップＳ２３）。すなわち、回路有り計算部６２は、現に設定された最小値γ_ｍｉｎと現に設定された最大値γ_ｍａｘとの間のいずれかの値に、１より大きな値（１．０＋η）を乗じて得た値をシフトパラメータγの値として算出する。そして、回路有り計算部６２は、シフトパラメータ計算処理を終了する。

[実施例１の効果]
上記実施例１によれば、磁化解析装置１は、磁場解析の過程で生成される連立方程式を構成する係数行列を分解して前処理行列を取得する。この前処理行列を取得する形態には、磁場の連立方程式を構成する係数行列を計算することにより取得することも含まれる。磁化解析装置１は、取得された前処理行列の対角成分の値を調整するシフトパラメータを用いて前処理行列の対角成分の値が正になるように、シフトパラメータの値を２分探索に基づいて算出する。かかる構成によれば、磁化解析装置１は、前処理行列の対角成分の値が正になるように、２分探索に基づいてシフトパラメータの値を算出するので、シフトパラメータの適切な値を高速に計算できる。この結果、磁化解析装置１は、磁場解析の過程で生成される連立方程式を安定且つ高速に計算できる。

また、上記実施例１によれば、磁化解析装置１は、連立方程式が回路に係る連立方程式を含む場合、回路に係る連立方程式に関係する、係数行列の行の対角成分を除外してシフトパラメータの値を算出する。かかる構成によれば、磁化解析装置１は、回路に関係する行の対角成分は零であるので、零であるとわかっている行の対角成分を予め除外してパラメータの値を算出するので、適切な値を算出できる。

また、上記実施例１によれば、磁化解析装置１は、シフトパラメータの値の最小値および最大値を予め設定する。磁化解析装置１は、現に設定されたシフトパラメータの最小値および現に設定されたシフトパラメータの最大値の中間値に対応する、前処理行列の対角成分の最小値を算出する。磁化解析装置１は、現に設定されたシフトパラメータの最小値と現に設定されたシフトパラメータの最大値との差分を算出する。磁化解析装置１は、前処理行列の対角成分の最小値が零より大きい場合に、現に設定されたシフトパラメータの最小値および現に設定されたシフトパラメータの最大値の中間値を新たなシフトパラメータの最大値に設定し、零以下である場合に、現に設定されたシフトパラメータの最小値および現に設定されたシフトパラメータの最大値との中間値を新たなシフトパラメータの最小値に設定する。磁化解析装置１は、算出された差分が閾値より小さくなるまで、シフトパラメータの算出処理を繰り返す。かかる構成によれば、磁化解析装置１は、シフトパラメータの探索範囲を半分に狭めながら、前処理行列の対角成分の値が正になるようにシフトパラメータの値を算出するので、シフトパラメータの適切な値を高速に計算できる。

また、上記実施例１によれば、磁化解析装置１は、算出された差分が閾値より小さい場合に、現に設定されたシフトパラメータの最小値と最大値との間の値に、１に調整パラメータηの値を加算して得られた値を乗じて得た値をシフトパラメータの値として算出する。かかる構成によれば、磁化解析装置１は、調整パラメータηに基づいて、シフトパラメータの値を微調整することができる。

また、上記実施例１によれば、磁化解析装置１は、閾値εを、調整パラメータとして予め定められた一変数ηの線形式ε＝kη（０＜ｋ＜１）を用いて算出する。かかる構成によれば、磁化解析装置１は、シフトパラメータの値を微調整する調整パラメータηを用いて閾値εを算出するので、より適切な値のシフトパラメータを算出できる。この結果、磁化解析装置１は、少ない反復回数で安定して磁場の連立方程式の解を求めることができる。

ここで、磁化解析装置１が、少ない反復回数で安定して磁場の連立方程式の解を求めることができるという効果について、図７Ａおよび図７Ｂを参照して説明する。図７Ａは、調整パラメータηと閾値εを独立に設定した場合の反復回数を示す図である。図７Ｂは、調整パラメータηと閾値εを連動させて設定した場合の反復回数を示す図である。なお、図７Ａおよび図７Ｂにおいて、「収束せず」とは、該当するシフトパラメータでは、磁場の連立方程式の反復計算が収束しないということを意味する。

調整パラメータηと閾値εを独立に設定した場合の反復回数は、以下のようになる。図７Ａに示すように、閾値εの設定によっては、収束しない場合が多くなる。ここでは、閾値εが０．１である場合には、調整パラメータηが０．０１、０．０２、０．０５、０．１である場合に収束しない。また、閾値εが０．００１の場合には、収束するが、調整パラメータηの値を増やすと、反復回数の増加が大きくなってしまう。

これに対して、調整パラメータηと閾値εを連動させて設定した場合の反復回数は、以下のようになる。図７Ｂに示すように、閾値εを調整パラメータηとｋの線形式から定めると、ｋが１．０である場合を除き、必ず収束している。少ない反復回数を１０未満とすると、ｋが０．５〜０．８とすることで、少ない反復回数で安定に磁場の連立方程式の解を求めることができる。

ところで、実施例１に係る磁化解析装置１は、並列計算を用いないで、係数行列からシフトパラメータを算出する場合を説明した。しかしながら、磁化解析装置１は、並列計算を用いて、係数行列からシフトパラメータを算出しても良い。

そこで、実施例２では、磁化解析装置１は、並列計算を用いて、係数行列からシフトパラメータを算出する場合について説明する。

図８は、実施例２に係る磁化解析装置の構成を示す機能ブロック図である。なお、図１に示す磁化解析装置１と同一の構成については同一符号を示すことで、その重複する構成および動作の説明については省略する。実施例１と実施例２とが異なるところは、制御部５のシフトパラメータ計算部６Ａを変更した点にある。実施例１と実施例２とが異なるところは、制御部５にシフトパラメータ共有部９を追加した点にある。

シフトパラメータ計算部６Ａは、計算条件データ４１の係数行列Ａの中の、並列計算するＣＰＵ（Central Processing Unit）に対応させて分担する分担領域について、対角成分が零よりやや大きい正の値となるように調整すべく、シフトパラメータを計算する。シフトパラメータは、並列計算する各ＣＰＵによって、各担当領域について算出される。

ここで、係数行列Ａの中のＣＰＵの分担の一例を、図９を参照して説明する。図９は、並列計算する場合のＣＰＵの分担の一例を示す図である。図９では、２つのＣＰＵが２並列で並列計算する場合を示す。図９に示すように、係数行列Ａの四角の黒塗りの領域は、値が零でないことを示す。係数行列Ａの四角の白塗りの領域および何もない領域は、値が零であることを示す。２つのＣＰＵが２並列で並列計算する場合には、一般に、係数行列Ａは、零でない値の数がほぼ等しくなるように分割される。ＣＰＵ１およびＣＰＵ２は、それぞれ四角の黒塗りのある行を処理する。ここでは、ＣＰＵ１は、領域ａ１、領域ａ３を分担する。ＣＰＵ２は、領域ａ２、領域ａ４を担当する。

ＣＰＵ１およびＣＰＵ２は、それぞれ担当する領域について、それぞれシフトパラメータγを計算する。シフトパラメータの計算方法は、実施例１と同様に、電流を未知数とする行を除外して、残りのＤ_ｉｉが正となるようにシフトパラメータγを計算する。

シフトパラメータ共有部９は、シフトパラメータ計算部６ＡによってＣＰＵ毎に計算されたシフトパラメータを取得し、取得したシフトパラメータのうち最大のシフトパラメータを共有する。すなわち、それぞれのＣＰＵで計算されるシフトパラメータは、同値とは限らない。そこで、シフトパラメータ共有部９は、各ＣＰＵで計算されたシフトパラメータのうち最大のシフトパラメータを共有し、共有するシフトパラメータをシフトパラメータデータ４２として記憶部４に格納する。最大のシフトパラメータの値を共有することで、前処理用行列計算部７および反復解法計算部８は、安定的に磁場の連立方程式の解を求めることができる。

なお、実施例２では、シフトパラメータ計算部６Ａは、２つのＣＰＵが並列計算する場合を説明した。すなわち、シフトパラメータ計算部６Ａは、係数行列Ａの中の２つのＣＰＵがそれぞれ分担する分担領域について、対角成分が零よりやや大きい正の値となるように調整すべく、シフトパラメータを計算する。しかしながら、シフトパラメータ計算部６Ａは、これに限定されず、３つ以上のＣＰＵが並列計算する場合であっても良い。また、シフトパラメータ計算部６Ａは、１つのＣＰＵであっても、１つのプロセスで複数のスレッドが並列計算する場合であっても良い。

［実施例２の効果］
上記実施例２によれば、磁化解析装置１は、前処理行列Ｃの第１の対角成分の値が正になるように、シフトパラメータの値を２分探索に基づいて算出する。磁化解析装置１は、前処理行列Ｃの第２の対角成分の値が正になるように、シフトパラメータの値を２分探索に基づいて算出する。磁化解析装置１は、算出された２つのパラメータの値の最大値を共有する。かかる構成によれば、磁化解析装置１は、前処理行列Ｃの第１の対角成分と第２の対角成分に担当を分けて、それぞれシフトパラメータを算出するので、高速にシフトパラメータを計算できる。また、磁化解析装置１は、算出された２つのパラメータの値の最大値を共有することで、磁場の連立方程式を安定且つ高速に計算できる。

［その他］
なお、図示した磁化解析装置１の各構成要素は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、磁化解析装置１の分散・統合の具体的態様は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。例えば、回路無し計算部６１と回路有り計算部６２とを１つの部として統合しても良い。また、記憶部４を磁化解析装置１の外部装置としてネットワーク経由で接続するようにしても良い。

また、上記実施例で説明した各種の処理は、予め用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータで実行することによって実現することができる。そこで、以下では、図１に示した磁化解析装置１と同様の機能を実現する磁化解析プログラムを実行するコンピュータの一例を説明する。図１０は、磁化解析プログラムを実行するコンピュータの一例を示す図である。

図１０に示すように、コンピュータ２００は、各種演算処理を実行するＣＰＵ２０３と、ユーザからのデータの入力を受け付ける入力装置２１５と、表示装置２０９を制御する表示制御部２０７とを有する。また、コンピュータ２００は、記憶媒体からプログラムなどを読取るドライブ装置２１３と、ネットワークを介して他のコンピュータとの間でデータの授受を行う通信制御部２１７とを有する。また、コンピュータ２００は、各種情報を一時記憶するメモリ２０１と、ＨＤＤ２０５を有する。そして、メモリ２０１、ＣＰＵ２０３、ＨＤＤ２０５、表示制御部２０７、ドライブ装置２１３、入力装置２１５、通信制御部２１７は、バス２１９で接続されている。

ドライブ装置２１３は、例えばリムーバブルディスク２１１用の装置である。ＨＤＤ２０５は、磁化解析プログラム２０５ａおよび磁化解析関連情報２０５ｂを記憶する。

ＣＰＵ２０３は、磁化解析プログラム２０５ａを読み出して、メモリ２０１に展開し、プロセスとして実行する。かかるプロセスは、磁化解析装置１の各機能部に対応する。磁化解析関連情報２０５ｂは、計算条件データ４１、シフトパラメータデータ４２および結果データ４３に対応する。そして、例えばリムーバブルディスク２１１が、磁化解析プログラム２０５ａなどの各情報を記憶する。

なお、磁化解析プログラム２０５ａについては、必ずしも最初からＨＤＤ２０５に記憶させておかなくても良い。例えば、コンピュータ２００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」に当該プログラムを記憶させておく。そして、コンピュータ２００がこれらから磁化解析プログラム２０５ａを読み出して実行するようにしても良い。

１ 磁化解析装置
２ 入力部
３ 出力部
４ 記憶部
４１ 計算条件データ
４２ シフトパラメータデータ
４３ 結果データ
５ 制御部
６ シフトパラメータ計算部
６１ 回路無し計算部
６２ 回路有り計算部
７ 前処理用行列計算部
８ 反復解法計算部

Claims (7)

1. 磁場の連立方程式の解を反復計算で計算する磁化解析装置において、
   前記連立方程式を構成する係数行列を分解して前処理行列を取得する取得部と、
   取得された前記前処理行列の対角成分の値を調整するパラメータを示すシフトパラメータの値を２分探索に基づいて算出する算出部と、を有し、
   前記算出部は、
   前記シフトパラメータの値の第１の値および第２の値を予め設定する第１の設定部と、
   設定された前記第１の値および設定された前記第２の値の中間値に対応する、前記前処理行列の対角成分の最小値を算出する第１の算出部と、
   設定された前記第１の値と設定された前記第２の値との差分を算出する第２の算出部と、
   算出された前記前処理行列の対角成分の最小値が零より大きい場合、設定された前記第１の値および設定された前記第２の値の中間値を新たな第２の値に設定し、前記最小値が零以下である場合、設定された前記第１の値および前記設定された第２の値との中間値を新たな第１の値に設定する第２の設定部と、を含み、
   算出された前記差分が所定の閾値より小さくなるまで、前記第１の算出部による処理、前記第２の算出部による処理および前記第２の設定部による処理を繰り返す
   ことを特徴とする磁化解析装置。
2. 前記第１の算出部は、前記連立方程式が回路に係る連立方程式を含む場合、前記前処理行列の対角成分のうち前記回路に係る対角成分を除外した、前記前処理行列の対角成分の最小値を算出する
   ことを特徴とする請求項１に記載の磁化解析装置。
3. 前記算出部は、
   算出された前記差分が前記所定の閾値より小さい場合、設定された前記第１の値と設定された前記第２の値との間の値に、１に前記シフトパラメータを調整するパラメータを示す調整パラメータの値を加算して得られた加算値を乗じて得た乗算値を前記シフトパラメータの値として算出する第３の算出部をさらに含む
   ことを特徴とする請求項１に記載の磁化解析装置。
4. 前記第３の算出部は、
   前記所定の閾値を、前記調整パラメータとして予め定められた一変数の線形式を用いて算出する
   ことを特徴とする請求項３に記載の磁化解析装置。
5. 前記算出部は、
   取得された前記前処理行列の第１の対角成分の値が正になるように、前記シフトパラメータの値を２分探索に基づいて算出する第１の算出制御部と、
   取得された前記前処理行列の第２の対角成分の値が正になるように、前記シフトパラメータの値を２分探索に基づいて算出する第２の算出制御部と、
   前記第１の算出制御部が算出したシフトパラメータの値と、前記第２の算出制御部が算出したシフトパラメータの値の最大値を共有する
   ことを特徴とする請求項３または請求項４に記載の磁化解析装置。
6. コンピュータが、
   磁場の連立方程式を構成する係数行列を分解して前処理行列を取得し、
   取得された前記前処理行列の対角成分の値を調整するパラメータを示すシフトパラメータの値の第１の値および第２の値を予め設定し、
   設定された前記第１の値および設定された前記第２の値の中間値に対応する、前記前処理行列の対角成分の最小値を算出するとともに、設定された前記第１の値と設定された前記第２の値との差分を算出し、
   算出された前記前処理行列の対角成分の最小値が零より大きい場合、設定された前記第１の値および設定された前記第２の値の中間値を新たな第２の値に設定し、前記最小値が零以下である場合、設定された前記第１の値および前記設定された第２の値との中間値を新たな第１の値に設定し、
   算出された前記差分が所定の閾値より小さくなるまで、前記算出する処理、前記中間値を新たに設定する処理を繰り返す
   処理を実行することを特徴とする磁化解析方法。
7. コンピュータに、
   磁場の連立方程式を構成する係数行列を分解して前処理行列を取得し、
   取得された前記前処理行列の対角成分の値を調整するパラメータを示すシフトパラメータの値の第１の値および第２の値を予め設定し、
   設定された前記第１の値および設定された前記第２の値の中間値に対応する、前記前処理行列の対角成分の最小値を算出するとともに、設定された前記第１の値と設定された前記第２の値との差分を算出し、
   算出された前記前処理行列の対角成分の最小値が零より大きい場合、設定された前記第１の値および設定された前記第２の値の中間値を新たな第２の値に設定し、前記最小値が零以下である場合、設定された前記第１の値および前記設定された第２の値との中間値を新たな第１の値に設定し、
   算出された前記差分が所定の閾値より小さくなるまで、前記算出する処理、前記中間値を新たに設定する処理を繰り返す
   処理を実行させることを特徴とする磁化解析プログラム。

Images