トップ :: G 物理学 :: G06 計算;計数

【発明の名称】 連立一次方程式に関する計算装置
【発明者】 【氏名】江澤 良孝

【氏名】マーチン・フィールド

【氏名】後 保範

【要約】 【課題】短時間で高精度の計算を導く前処理行程を有する連立一次方程式に関する計算装置を提供する。

【解決手段】連立一次方程式の係数行列Aに関する情報が入力されるデータ入力部、前記係数行列Aに対する前処理行列Cを求める前処理部、前記係数行列Aと前記求めた前処理行列Cを用いて演算する演算部、前記演算結果の出力部を備え、前記前処理部は、前記係数行列Aを部分行列Dに分解する第一の処理部、前記分解された部分行列Dのノルムを計算する第二の処理部、前記計算された部分行列Dのノルムの値のうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する第三の処理部、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する第四の処理部、前記追加された行列を演算して前記前処理行列Cを求める第五の処理部を有する。
【特許請求の範囲】
【請求項1】連立一次方程式の係数行列Aに関する情報が入力されるデータ入力部と、前記係数行列Aに対する前処理行列Cを求める前処理部と、前記入力された係数行列Aと前記求めた前処理行列Cの値を用いて連立一次方程式を解く演算部と、前記演算結果を出力する出力部とを備え、前記前処理部は、前記係数行列Aを部分行列Dに分解する第一の処理部と、前記分解された部分行列Dのノルムを計算する第二の処理部と、前記計算された部分行列Dのノルムの値のうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する第三の処理部と、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する第四の処理部と、前記追加された行列を演算して前記前処理行列Cを求める第五の処理部と、を有することを特徴とする連立一次方程式に関する計算装置。
【請求項2】連立一次方程式の係数行列Aに関する情報が入力されるデータ入力部と、前記係数行列Aの逆行列の近似行列C~1を求める前処理部と、前記入力された係数行列Aと前記求めた前近似行列C~1を用いて連立一次方程式を解く演算部と、前記演算結果を出力する出力部とを備え、前記前処理部は、前記係数行列Aのマトリックスを構成する部分行列Dの演算値が予め記憶された設定値より小さい成分を削除して第二のマトリックスに相当するスパースパターンを作成する第一の演算部と、前記第二のマトリックスを計算した場合に零から非零に変わる成分を前記スパースパターンに追加して第二のスパースパターンを作成する第二の演算部と、前記第二のスパースパターンを用いて前記近似行列C~1を求める第三の演算部と、を有することを特徴とする連立一次方程式に関する計算装置。
【請求項3】連立方程式の係数行列Aに対する前処理行列Cを計算する装置であって、係数行列A、計算結果の前処理行列Cおよび途中計算結果を記憶する手段と、前記記憶手段を参照しながら、係数行列Aを部分行列Dの集合に分解し、前記部分行列Dの成分からノルムを計算し、前記ノルムが予め設定された設定値以上のもののみの部分行列の位置を記憶する手段と、前記位置を示す行列Sを対角項から上側の成分がゼロある下三角行列Lと対角項から下側の成分がゼロである上三角行列Uの積に分解する三角分解を行った場合、元々非零の成分の影響でゼロから非零に変化する前記部分行列Dの成分の位置を計算し、前記部分行列Dの成分の位置情報を前記位置を示す行列Sに追加した結果を記憶する手段と、前記追加された結果を反映した位置を示す行列S’が示す非零部分行列を用いて不完全三角分解を行う処理手段と、を有することを特徴とする連立方程式に関する計算装置。
【請求項4】請求項1乃至3において、前記部分行列Dは、サイズが1×1、2×2或いは3×3から選択されることを特徴とする連立方程式に関する計算装置。
【請求項5】請求項1乃至3において、前記部分行列Dに関する設定値は、0.05以上0.25以下であることを特徴とする連立方程式に関する計算装置。
【請求項6】請求項1乃至3において、前記連立方程式は、有限要素法による数値シミュレーションにおいて生成される連立一次方程式、或いは差分法による数値シミュレーションにおいて生成される連立一次方程式であることを特徴とする連立一次方程式に関する計算装置。
【請求項7】入力される連立一次方程式の係数行列Aに関するデータに基づき、前記係数行列Aに対する前処理行列Cを求める前処理行程と、前記入力された係数行列Aと前記求めた前処理行列Cの値を用いて連立一次方程式を解く演算行程と、前記演算結果を出力する出力行程と、を実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体であって、前記前処理行程は、前記係数行列Aを部分行列Dに分解する行程と、前記分解された部分行列Dのノルムを計算する行程と、前記計算された部分行列Dのノルムのうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する行程と、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する行程と、前記追加された行列を演算して前記前処理行列Cを求める行程と、を実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項8】入力される連立一次方程式の係数行列Aに関するデータに基づき、前記係数行列Aに対する前処理行列Cを求める前処理行程と、前記入力された係数行列Aと前記求めた前処理行列Cの値を用いて連立一次方程式を解く演算行程と、前記演算結果を出力する出力行程と、を有するデータ処理方法であって、前記前処理行程は、前記係数行列Aを部分行列Dに分解する行程と、前記分解された部分行列Dのノルムを計算する行程と、前記計算された部分行列Dのノルムのうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する行程と、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する行程と、前記追加された行列を演算して前記前処理行列Cを求める行程と、を有することを特徴とするデータ処理方法。
【発明の詳細な説明】【0001】
【発明の属する技術分野】本発明は、連立一次方程式の数値解を計算する装置に関し、連立一次方程式の係数行列に対する前処理行列を計算する行程を有する装置に係る。
【0002】
【従来の技術】従来の技術としては、村田健郎他著「スーパーコンピュータ、科学技術計算への応用」(丸善)に論じられている前処理つき共役勾配法(Preconditioning CGmethod、PCG法)がある。
【0003】共役勾配法とは、Aを既知のマトリックス、xを未知ベクトル、bを既知ベクトルとしたとき、連立方程式Ax=bの解が関数f(x)=1/2(x,Ax)−(x,b)
を最小にするという事実に基づいている。実際、f(x+h)=f(x)+(h,Ax−b)+1/2(h,Ah)
であり、Aが正定値だから(h,Ah)>0となる。したがって、Ax=bの解x^に対してf(x^+h)>f(x^)
すなわち、どんなh(≠0)を選んでも、f(x^+h)はf(x)^となるので、f(x^)は最小となる。
【0004】共役勾配法は任意の初期値x0から始めて、f(x)の最小点を逐次探索する方法で、いま、第k番目の近似値xkと探索方向pkが決まったとするとxk+1=xk+αpkとしてf(xk+1)が最小になるようにαkを定めればよい。それには【0005】
【数1】

【0006】数式(1)をαkの関数とみて、それが最小になるようにαkを定める。そうすると【0007】
【数2】

【0008】となる。ただし、rkは第k近似解に対する残差であり、次のように定義される。
【0009】rk+1=b−Axkつぎに次のステップの探索方向をpk+1=rk+1+βkkによって定める。 ただし、p0=r0である。ここで、βkを定める方法として、pk+1とApkが直交するようにβkを定める。
【0010】
【数3】

【0011】から【0012】
【数4】

【0013】となる。以上が、共役勾配法の大筋である。
【0014】このような近似解から正解を反復計算で求める手法を反復解法と呼ぶ。これに対し、連立方程式の未知数を、順に式の和や差によって消去して求める方法が直接解法と呼ばれるもので、その代表的なものにガウス消去法(村田健郎、名取亮、唐木幸比古著「大型数値シミュレーション」岩波書店)などがある。直接法と反復解法を比べると、反復解法の方が必要とするメモリ量が圧倒的に小さいという利点をもっている。
【0015】
【発明が解決しようとする課題】しかしながら、反復解法には解の収束性の問題がつきまとっている。
【0016】すなわち、共役勾配法では、{ri}の直交性から、数式の上ではあるm≦n(n:マトリックスAのサイズ)で残差はゼロになる。すなわち、たかだかn回の反復で正解が得られる。しかしながら、実際には丸め誤差があるのでn回以下の反復で正解が得られないことがある。具体的には、構造物の変形や応力を有限要素法で求める際には、材料が降伏したり、変形が大きくなることによりこの収束性が悪化することがあり、シミュレーションに適用する際の問題となっていた。
【0017】共役勾配法の収束の速さは、Aの固有値分布と深い関係がある。すなわち、Aの固有値に重複があるときには収束が速い。完全に重複していなくても、固有値が密集していると収束が速い。その理由は以下の通りである。いま、Aの固有値、固有ベクトルをλi、νiとし、初期残差r0を固有ベクトルで展開して、【0018】
【数5】

【0019】としよう。ただし、(νij)=δijとする。ここに【0020】
【数6】

【0021】
【数7】

【0022】
【数8】

【0023】となるので、rk=Rk(A)r0と書けることが分かる。ここに、R0(0)=1である。
【0024】f(x)を残差を用いて書き直すとf(x)=1/2(x,Ax)−(x,b)
=1/2(r,A~1r)−1/2(b,A~1b)
となるから、f(x)最小にすることと(r,A~1r)を最小にすることは同じである。したがって、多項式Rk(A)は各kにおいて(r,A~1r)を最小にするように決められていると考えてよい。すなわち、【0025】
【数9】

【0026】となる。ここで、П1kはk次多項式でRk(0)=1を満たすものの集合である。
【0027】ここで、【0028】
【数10】

【0029】となるので、もし、λp、λp+1、…、λp+qが重複していればRk(λp)=Rk(λp+1)=…=R(λp+q
であるから、rkの動きうる空間の次元がqだけ落ちる。したがって、たかだかn-q回の反復で収束する。
【0030】このような原理を用いて、この共役勾配法の収束を速くするために、事前に方程式を変形しておいてから共役勾配法を適用するものがある。これを前処理付き共役勾配法(PCG法)と呼ぶ。
【0031】具体的には、Aに近い適当な行列Cを選んで、Ax=bの代わりにC~1x=C~1bに対して共役勾配法(CG法)を適用するものである。CがAに近ければ、C~1Aは単位行列に近くなり、固有値が密集するので解の収束は速くなる。
【0032】実際のPCG法のアルゴリズムは、C~1Aを直接計算しないで、下記のように残差ベクトルrにC~1を掛けるようにする。
【0033】
【数11】

【0034】これはC~1Aの計算を直接やるとサイズn×nの行列の積になり、nの三乗の回数の積を計算しなければならず、計算量が膨大になるからである。
【0035】したがって、PCG法のアルゴリズムは以下のようになる。
【0036】(1) xの近似解x0を用意(2) r0=b−Ax0 ; p0=C~10 ; v=0(3)【0037】
【数12】

【0038】
【数13】

【0039】
【数14】

【0040】
【数15】

【0041】
【数16】

【0042】
【数17】

【0043】この前処理つき共役勾配法の代表的なものが、C~1の計算に、三角分解の一種であるコレスキー分解を使用するICCG法である。この方法は、不完全コレスキー分解(IC)をまず解くべき方程式の行列に施し、行列を単位行列に近づけてから共役勾配法(CG法)を適用するものである。三角分解とは、行列Aを対角項の上側の成分がすべてゼロである下三角行列Lと、対角項の下側の成分がすべてゼロである上三角行列Uの積LU(=A)に分解するものである。また、不完全三角分解とはA=LU+Rのように非零の残差Rが生じるものをいう。不完全コレスキー分解ではA=LL~1+Rとなる。不完全コレスキー分解の計算では、通常のコレスキー分解手続きを用いるが、すべてのコレスキー分解手続きをすべて行わないで、適当に手抜きをして行う。すなわち、Aの非零要素の場所の集合(スパースパターン)をSA、すなわちSA={(i,j):aij≠0}とし、S⊇SAなる集合を予めきめておいて、Aをコレスキー分解するときに、L、LTの要素中Sに属する場所のものだけを計算し、他のものはゼロにする。これを不完全コレスキー分解という。
【0044】ICCG法における不完全コレスキー分解の具体的手順は以下の通りである。
【0045】
【数18】

【0046】この方法の問題点は、不完全コレスキー分解(IC)の処理が、複数の計算処理装置(プロセッサ)をもつ並列計算機への適用が困難なことである。すなわち、上記手続きの中の【0047】
【数19】

【0048】は、処理の完全な並列化が困難であり、あるプロセッサでの計算に他のプロセッサの計算結果が必要になる。したがって、使用するプロセッサが増えるに従い、プロセッサ間のデータ転送の増加が避けられない。そのため、サイズの大きい係数行列を計算するためにプロセッサ数を多くすると、データ転送による計算時間が多くなり、並列化効率が上がらないという問題が生じる。この問題を解決するため、CG法に代わって、非対称行列用の反復解法であるBi-CGSTAB法を使う方法等が提案されてきた。Bi-CGSTAB法は積型の反復解法であり、前処理は規則シャドウと呼ばれる、隣接するプロセッサ間のみでデータ転送を行う手法で並列化できる。しかしながら、本来非対称行列用の手法であるため反復ごとの計算量が多いのと必要とする記憶容量が大きいという問題を抱えており、その計算量の増加を打ち消すだけの前処理の効果が必要になる。
【0049】有限要素法で使われる係数行列は大部分対称行列であり、その対称行列用の反復解法であるCG法を用いながら並列化処理のできる前処理の方法が望まれる。このような目的で開発された手法にFSAI(Factorized Sparse Approximate Inverse)法がある。この手法はKolotilinaとYeremin(ロシア、モスクワ大学)によって提案されたもので、その詳細は文献「L.Yu. Kolotilina and A.Yu. Yeremin(ロシア、モスクワ大): Factorised sparse approximate inversepreconditioning I:Theory; SIAM Journal of Matrix Analysis andApplications, 14, pp.45-58(1993)」に記載されている。この方法は以下のようなものである。
【0050】SLを下記で定義されるAの下三角スパースパターン(行列の成分の非零の位置を示す行列)とする。
【0051】
【数20】

【0052】スパースパターンSLは具体的には行列Aの非零の成分の位置に対応する行列Sの成分には、非零であることを示す数値(1でも0でもよい)なり記号なりを割り当て、行列Aの零の成分の位置に対応する行列Sの成分には零であることを示す数値なり記号なりを割り当てる。
【0053】さらに、MSLを下記のようにこのスパースパターンSLに含まれない成分をゼロとし、それ以外の成分は未知数xであるN×Nの未知マトリックスとする。
【0054】
【数21】

【0055】ここで関数【0056】
【数22】

【0057】を考える。この関数を下記のように最小化することにより、コレスキー分解における下三角マトリックスLの逆マトリックスのスパース近似(Sparse ApproximateInverse)【0058】
【数23】

【0059】を求めることができる。
【0060】
【数24】

【0061】ここに、【0062】
【数25】

【0063】であることから、これをXで偏微分し、【0064】
【数26】

【0065】とすることにより【0066】
【数27】

【0067】なる連立方程式が得られるが、上記の式ではコレスキー分解後のマトリックスLの対角項の値を必要とする。Lの計算をするということは、コレスキー分解をするということであり、通常の直接解法を使ったのと変わりなくなくなる。そこで、上記のG〜の代わりに、以下の式で定義されるマトリックスG∈MSLを用いることにする。
【0068】
【数28】

【0069】ここで以下のような対角マトリックスDを導入する。
【0070】D=Diag(G)
このDを用いると【0071】
【数29】

【0072】から【0073】
【数30】

【0074】ここで【0075】
【数31】

【0076】だから【0077】
【数32】

【0078】となる。これを前処理C~1AにおけるC~1として用いる。すなわち【0079】
【数33】

【0080】となる。このC~1を既に述べた前処理付きCG法に適用すればよい。
【0081】Gを求めるために解くべき式は【0082】
【数34】

【0083】であるが、マトリックスGのi行目の非ゼロ成分だけからなるベクトルをgiとし、Ji={j:(i,j)∈SL} (ni:Jiのサイズ)とすると、GA=I (I :単位マトリックス)は【0084】
【数35】

【0085】
【数36】

【0086】
【数37】

【0087】となることから、ベクトルgiはni×niの小さな連立方程式AJii=eni (1)
の解となることが分かる。ここに、AJiはni×niのAのサブマトリックス、eniはni×niの単位マトリックスのni列目のベクトルである。このni×niの連立方程式生成の様子を示したのが、図1である。図1の太枠で囲まれた部分の対でひとつの方程式が生成される。式(1)はn個の独立な式であり、n個のプロセッサで相互のデータ転送なしに完全な並列処理が可能である。
【0088】前記のFSAI法は完全な並列処理が可能であるが、問題はこの手法でどれだけ前処理の効果が得られるかという点にある。図2は187,280節点、561,840自由度の有限要素法構造解析の問題に適用した結果であるが、行列Aの対角項のオーダをそろえるスケーリングと呼ばれる処理のみをした場合と比べると計算時間の短縮は約1割から2割であり、前処理の効果は十分とは言い難い。
【0089】そこで、本発明は、短時間で高精度の計算を導く前処理行程を有する連立一次方程式に関する計算装置を提供することにある。
【0090】
【課題を解決するための手段】前記の前処理効果の問題を本発明により解決することができる。前処理の効果は、前処理行列C~1を計算するときフィルイン(零の項が非零に変わること)をどれだけ考慮するかで決まってくる。すべてのフィルインを考慮すると、完全な三角分解をすることに相当し、直接解法をしたことと同じになる。逆にフィルインの考慮が少なすぎると、前処理に要する計算時間は短くなるが、前処理の効果は小さくなる。
【0091】よって、本発明の連立一次方程式に関する計算装置は、連立一次方程式の係数行列Aに関する情報が入力されるデータ入力部と、前記係数行列Aに対する前処理行列Cを求める前処理部と、前記入力された係数行列Aと前記求めた前処理行列Cの値を用いて連立一次方程式を解く演算部と、前記演算結果を出力する出力部とを備え、前記前処理部は、前記係数行列Aを部分行列Dに分解する第一の処理部と、前記分解された部分行列Dのノルムを計算する第二の処理部と、前記計算された部分行列Dのノルムの値のうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する第三の処理部と、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する第四の処理部と、前記追加された行列を演算して前記前処理行列Cを求める第五の処理部と、を有することを特徴とする。
【0092】または、連立一次方程式の係数行列Aに関する情報が入力されるデータ入力部と、前記係数行列Aの逆行列の近似行列C~1を求める前処理部と、前記入力された係数行列Aと前記求めた前近似行列C~1を用いて連立一次方程式を解く演算部と、前記演算結果を出力する出力部とを備え、前記前処理部は、前記係数行列Aのマトリックスを構成する部分行列Dの演算値が予め記憶された設定値より小さい成分を削除して第二のマトリックスに相当するスパースパターンを作成する第一の演算部と、前記第二のマトリックスを計算した場合に零から非零に変わる成分を前記スパースパターンに追加して第二のスパースパターンを作成する第二の演算部と、前記第二のスパースパターンを用いて前記近似行列C~1を求める第三の演算部と、を有するようにしてもよい。
【0093】または、本発明の連立一次方程式に関する計算装置は、連立方程式の係数行列Aに対する前処理行列Cを計算する装置であって、係数行列A、計算結果の前処理行列Cおよび途中計算結果を記憶する手段と、前記記憶手段を参照しながら、係数行列Aを部分行列Dの集合に分解し、前記部分行列Dの成分からノルムを計算し、前記ノルムが予め設定された設定値以上のもののみの部分行列の位置を記憶する手段と、前記位置を示す行列Sを対角項から上側の成分がゼロある下三角行列Lと対角項から下側の成分がゼロである上三角行列Uの積に分解する三角分解を行った場合、元々非零の成分の影響でゼロから非零に変化する前記部分行列Dの成分の位置を計算し、前記部分行列Dの成分の位置情報を前記位置を示す行列Sに追加した結果を記憶する手段と、前記追加された結果を反映した位置を示す行列S’が示す非零部分行列を用いて不完全三角分解を行う処理手段と、備えるようにしてもよい。
【0094】さらに、例えば、前記前処理行列Cを計算する手順において、上記手段により作成し、記憶装置に記憶された位置情報Sをもとに、連立一次方程式の係数行列Aの成分から位置情報Sに対応する成分のみを残して、他の成分はゼロとした行列A’を作成して記憶装置に記憶し、位置情報Sに対応した成分以外はゼロとした未知行列Gのi行目に対応する行ベクトルgiを行列A’の右側から掛けたものを左辺とし、i行目が1で他をゼロとしたベクトルを右辺とした式を各i行ごとに別々のプロセッサで解き、その結果をひとつのプロセッサに集めて未知行列Gを定め、そのGを係数行列Aの三角分解における下三角行列の逆マトリックスの近似値とみなすことにより前処理行列Cを計算するようにすることができる。
【0095】前記の連立一次方程式に関する計算装置においては、前記部分行列Dは、サイズが1×1、2×2或いは3×3から選択されることが好ましい。計算対象となるデータが3次元である場合或いは前記係数行列Aに関する情報が3次元のものに関する場合には前記3×3を選択することが好ましい。一方、計算対象となるデータが2次元である場合或いは前記係数行列Aに関する情報が2次元のものに関する場合には前記2×2を選択することが好ましい。
【0096】前記の連立一次方程式に関する計算装置においては、前記部分行列Dに関する設定値は、0.05以上0.25以下であってもよい。
【0097】また、例えば、具体的には、解かれる方程式が三次元の構造物の変位と力の関係を表す場合には、部分行列Dの成分から計算するノルムが一定値以上のものを抽出するノルムの基準値を、0.1以上0.2以下とすることにより高速化が図れて好ましい。
【0098】また、前記発明において、前記連立方程式は、有限要素法による数値シミュレーションにおいて生成される連立一次方程式、或いは差分法による数値シミュレーションにおいて生成される連立一次方程式であることを特徴とする。
【0099】また、他の本発明は、入力される連立一次方程式の係数行列Aに関するデータに基づき、前記係数行列Aに対する前処理行列Cを求める前処理行程と、前記入力された係数行列Aと前記求めた前処理行列Cの値を用いて連立一次方程式を解く演算行程と、前記演算結果を出力する出力行程と、を実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体であって、前記前処理行程は、前記係数行列Aを部分行列Dに分解する行程と、前記分解された部分行列Dのノルムを計算する行程と、前記計算された部分行列Dのノルムのうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する行程と、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する行程と、前記追加された行列を演算して前記前処理行列Cを求める行程と、を実行させるものである。
【0100】また、他の本発明は、入力される連立一次方程式の係数行列Aに関するデータに基づき、前記係数行列Aに対する前処理行列Cを求める前処理行程と、前記入力された係数行列Aと前記求めた前処理行列Cの値を用いて連立一次方程式を解く演算行程と、前記演算結果を出力する出力行程と、を有するデータ処理方法であって、前記前処理行程は、前記係数行列Aを部分行列Dに分解する行程と、前記分解された部分行列Dのノルムを計算する行程と、前記計算された部分行列Dのノルムのうち予め設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する行程と、前記記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する行程と、前記追加された行列を演算して前記前処理行列Cを求める行程と、を有するものである。
【0101】なお、前記FSAI法を用いる際、本発明において改良したスパースパターンSLの決め方の内容を以下に各ステップと共に詳述する。
【0102】(1)第1ステップ:値の小さい係数を削除したマトリックスを作成し、そのスパースパターンを作成。
【0103】(2)第2ステップ:上記のマトリックスを三角分解したときのフィルインを、上記のスパースパターンに追加する。このときフィルインは、計算時間を短縮するため一次のフィルイン、すなわち、それ自身初期に零で、もともと非零だった成分からの演算で生成される非零成分のみに限定する。
【0104】(3)第3ステップ:前期で生成されたスパースパターンを用いてFSAI法を適用する。
【0105】従来のFSAI法では下記のように剛性マトリックスKの非零のパターンをそのまま用いていた。
【0106】SL={(i,j):kij≠0,i≧j}収束までの反復数を低減するには、非零の項を増やせばよい。しかし、単純に増やすと前処理の計算時間が増えてしまうので、そのバランスを考えなければならない。上式の問題点は、LDLT分解時に発生するフィルインを無視していることである。前処理の効果が小さかったのは、そのためと考えられる。そこで、フィルインで非零に変わる項もスパースパターンに加えることを考える。
【0107】ところで、フィルインには、もともと非零の項の影響で零から非零になるもの(これを1次のフィルインと呼ぶ)、さらにその1次のフィルインの影響で発生するフィルイン(これを2次のフィルインと呼ぶ)などがある。すべてのフィルインを含めてしまうと先に述べたように直接法で解いているのと変わりなくなってしまう。そこで、本発明では、一次のフィルインのみを考える。
【0108】ただし、できるだけ無駄なフィルインの考慮を避けるため、値の小さな項は無視する。ここに、三次元構造物の解析では、1節点の自由度は3(シェル要素ではその倍の6)になることから、剛性マトリックスの小さな成分の削除は、3×3の部分マトリックスの単位で行うことにした。すなわち、まず、全体剛性マトリックスAをスケーリングする。例えば【0109】
【数38】

【0110】にスケーリングを行うと次のようになる。
【0111】
【数39】

【0112】次に、縦横3行ずつの単位で分割し、3×3の部分マトリックスA〜の集合に変える。
【0113】次にこの個別の部分マトリックスA〜の成分のノルムを計算する。すなわち、【0114】
【数40】

【0115】ただし、上式では、アインシュタインの総和規約を用いている。したがって、もしA〜が単位行列であれば、【0116】
【数41】

【0117】となる。
【0118】ここで、許容値εを定め、この値以下のK〜をスパースパターンSLから除去する。除去後のスパースパターンをStolLとする。
【0119】次に、1次のフィルインの追加を行う。
【0120】これは以下のようにすればよい。
【0121】
【数42】

【0122】以上の手順を図示すると図3のようになる。
【0123】このように定めたスパースパターンStolLをスパースパターンSLの代わりに用いてFSAI法を適用する。
【0124】すなわち、このスパースパターンStolLをもとに、連立方程式の係数行列Aの成分から位置情報Sに対応する成分のみを残して、他の成分はゼロとした行列A’を作成して記憶装置に記憶し、位置情報Sに対応した成分以外はゼロとした未知行列Gのi行目に対応する行ベクトルgiを行列A’の右側から掛けたものを左辺とし、i行目が1で他をゼロとしたベクトルを右辺とした式A’gi=Iiを各i行ごとに別々のプロセッサで解き、その結果をひとつのプロセッサに集めて【0125】
【数43】

【0126】であることから未知行列Gを定め、そのGを係数行列Aの三角分解における下三角行列の逆マトリックスの近似値とみなすことにより【0127】
【数44】

【0128】から、前処理行列Cを計算する。このCを親プロセッサの記憶装置に記憶する。つぎに親プロセッサは、このCを次の前処理つき共役勾配法のアルゴリズムの計算を開始する。
【0129】(1) xの近似解x0を用意(2) r0=b−Ax0 ; p0=C~10 ; v=0(3)【0130】
【数45】

【0131】
【数46】

【0132】
【数47】

【0133】
【数48】

【0134】
【数49】

【0135】
【数50】

【0136】この計算のなかでApvは、Aが各子プロセッサのメモリ上にある部分領域の係数行列の和であることを利用し、各部分領域の係数行列にかかる計算は各プロセッサで計算し、その計算結果を親プロセッサに転送する。親プロセッサは子プロセッサから送られてきたApvの結果を用いて、上記の計算を行う。
【0137】以上を図示すると、図4のようになる。
【0138】ここで問題となるのは許容値εをどのように定めるかということである。そこで、いくつかの例題について、εを変えて計算時間の比較を行った。
【0139】図5は有限要素法による15,492自由度の比較的小さな熱応力解析問題に適用した場合である。図6は、69,579自由度のより大きい熱応力解析問題の場合である。両方の場合とも許容値εは0.1から0.2のときに計算時間が最小となっている。許容値εを小さくすると、前処理に要する計算量は減るが、前処理の効果は小さくなるので、必要とする反復回数は増加する。逆に、許容値εを小さくすると、前処理の効果は大きくなるが前処理に要する計算量は増加する。そのため、このような極小値が出てくることになる。そこで、この数値実験から、0.1以上0.2以下を標準値として採用する。ただし、解きたい係数行列の行列式の値がゼロに近い場合は反復法では収束がし難いので、許容値εをもっと大きくして収束性を良くすることも考えられる。
【0140】
【発明の実施の形態】本発明の一実施例を以下に示す。
【0141】図7は、本発明の計算装置の概要を示す。(a)が概要を示し、(b)が(a)の内前処理部の内容の内容を示す。 本願発明の計算装置は例えば有限要素方による三次元や二次元の応力、変形状態、温度、他の状態変化をシミュレーション等する計算機に使用することができる。もっともそれ以外であっても、連立一次方程式を用いて解く事象について適応することができる。
【0142】本発明の一実施例の計算装置は、図7(a)に示すように、連立一次方程式の係数行列Aに関する情報が入力されるデータ入力部1と、その係数行列Aに対する前処理行列Cを求める前処理部2と、入力された係数行列Aと前処理部で求めた前処理行列Cの値を用いて連立一次方程式を解く演算部3と、この演算部での演算結果を出力する出力部4とを備える。
【0143】図示したように各行程での演算結果は一旦記憶部5に記憶されるようにして、次の行程で記憶した情報を使用する際には、記憶部から呼び出すようにすることが好ましい。
【0144】また、図7(b)に示すように、前記前処理部では、まず、データ入力部で入力されたデータを記憶部から取り出し、第一の処理部6において係数行列Aを部分行列Dに分解する。次に、第二の処理部7において第一の処理部で分解された部分行列Dのノルムが計算される。次に、第三の処理部8において第二の処理部で計算された部分行列Dのノルムの値と記憶部に記憶されているノルムの設定値とを比較する。さらに、第二の処理部で計算された部分行列Dのノルムのうち前記設定されたノルムの値より大きい前記部分行列Dを構成する成分の位置を記憶する。次に第四の処理部9において、第三の処理部での処理により記憶された成分の位置に基づいて形成される行列Sを演算した際にゼロから非ゼロになる行列成分の位置情報を前記行列Sに追加する。そして、第五の処理部10において第四の処理部で追加された行列を演算して前記前処理行列Cを求める。そして、それを記憶部に記憶する。
【0145】前記(a)にも述べたように、(b)においても 各処理部での処理結果は一旦記憶部に記憶されるようにして、次の行程で記憶した情報を使用する際には、記憶部から呼び出すようにすることが好ましい。
【0146】このように構成して、計算が短時間で済む程度であって精度の低下を抑制するように計算対象を絞ることにより、短時間で高精度の計算ができる計算装置を提供することができる。
【0147】次に本願発明の一実施例における第二の例を以下に示す。
【0148】図8は、有限要素法により三次元の構造物の応力、変形状態を並列計算機でシミュレーションするシステムのフローチャートを示したものである。プログラムは、図9に示すように販売時点では、通常、CD−ROMか、磁気テープ、フロッピディスク、MOディスク等の記録媒体に記録された状態になっており、使用する前にプログラムをハードディスク上にインストールする。起動時には、ハードディスク上にあるプログラムは、各プロセッサのメモリ上にロードされる。プログラムはすべてのプロセッサのメモリに同じプログラムが分配される。ただし、実際に動作するときは、各プロセッサに分配されたプログラムは、自分が何番のプロセッサ上にいるのかを検知して、そのプロセッサ番号に応じた別々の動作をする。
【0149】複数のプロセッサでどのようにロードを分散させるかを、図10を用いて説明する。
【0150】まず、図10の左に示す全体領域を、各プロセッサが担当する複数の部分領域に分ける。この場合、分けかたとして節点数ができるだけ同じになるように分ける方法と、要素数が同じになるように分ける方法が考えられる。有限要素法による応力解析においては、要素生成に要する時間より、連立方程式を解く時間の方が大部分を占めるのが普通である。そこでプロセッサの負荷をできるだけ均等にするため、節点数を基準にして領域分割を行う。各プロセッサはその部分領域の要素剛性マトリックスを生成する。この場合問題となるのは、各部分領域にまたがる要素であるが、各プロセッサはその境界部の要素の構成節点の情報も含めてもつようにし、その境界部の要素の剛性マトリックスも生成するようにする。図10の右側が分割例で、同じ濃さの部分がひとつのプロセッサで担当する部分である。
【0151】実際の製品の解析モデルの領域分割例を図11に示す。濃さで分けられた各部分が各プロセッサで担当する部分である。
【0152】図12に示すように、まず親プロセッサのメモリにロードされたプログラムはハードディスク上にある入力データファイルを読み込む。入力データには、シミュレーションしたい構造物の形状データ、物性データ、荷重データ、構造を複数の細かい領域(有限要素)に分けるデータなどが含まれている。つぎに、分割された部分領域のデータをそれぞれを担当する子プロセッサに転送する。
【0153】各プロセッサ上のプログラムは入力データを用いて、担当する部分領域における有限要素の変位と荷重の関係を示す下記のような関係式(要素剛性マトリックス)を作成する。
【0154】[Ae]{u}={f}
ここに{u}は要素を構成する節点における変位で{f}はその節点における荷重を表している。つぎに、この各要素における関係式をつなげて下記のような構造全体の関係式(全体剛性マトリックス)を組み立てる。
【0155】[A]{u}={f}
本発明はこの全体剛性方程式を解く段階で適用する。全体剛性マトリックスを3行3列ごとに分割し、各3×3の部分マトリックスごとにノルムを計算し、そのノルムが0.15以下の部分マトリックスは削除する。残った非零成分の位置から、スパースパターンを作成する。つぎに、この削除した後のマトリックスから一次のフィルインを計算し、その位置をスパースパターンに追加する。このようにして定めたStolLを用いて、【0156】
【数51】

【0157】を解く。具体的には、マトリックスGのi行目の非ゼロ成分だけからなるベクトルをgiとし、ベクトルgiをni×niの小さな連立方程式AJii=eni (1)
の解として求める。ここに、AJiはni×niのAのサブマトリックス、eniはni×niの単位マトリックスのni列目のベクトルである。
【0158】N個の独立な式(1)は完全に並列に解くことができる。しかも、この方程式は図1から分かるように一般に非常に小さなサイズになっており、少ない計算量で解くことが出来る。例えばガウスの消去法などの直接法で解けばよい。また、ICCG法においては前処理における丸め誤差の蓄積が問題となるが、本手法ではAJiのサイズが小さく、しかも独立に計算できるということから、前処理において丸め誤差の蓄積が起きにくいという特長がある。条件数の悪いマトリックス(例えば、拘束の弱い構造物)は丸め誤差の影響を受け易いので、その点でもICCG法より優れている。
【0159】子プロセッサで剛性マトリックスが生成されたのを親プロセッサが確認すると、親プロセッサは全体剛性方程式の求解を開始する。ただし、この際、実際に全体剛性マトリックス[A]を組みたてることはしないで、共役勾配法の計算手順のなかで全体剛性マトリックス[A]の関係する計算部分[A]{pv}は[A]が部分領域の剛性マトリックス[Ai]の総和であることから [A]{pv}≡([A1]+[A2]+…+[Ak]){pv} (k:部分領域の数)
=[A1]{pv1}+[A2]{pv2}+…+[Ak]{pvk}として、各プロセッサで[Ai]{pvi}を計算して、その結果のベクトルを親プロセッサに送るようにする。このようにすることにより、計算とメモリを子プロセッサに分担させることができる。前処理用のマトリックスCを計算するための式AJii=eniも同様に子プロセッサに分担させて計算する。親プロセッサは子プロセッサが計算した結果を集めて、共役勾配法の計算を行う。親プロセッサがこのようにして求めた節点変位、節点荷重は、それぞれの点が所属する領域を担当する子プロセッサに送り返す。子プロセッサは送られてきた節点変位データ、節点荷重データを用いて、担当領域の有限要素の応力、ひずみを計算する。その計算結果は、また親プロセッサに送り返す。親プロセッサは収集した応力、ひずみデータと自分が所有している節点変位、節点荷重データを合わせて、出力データファイルを作成し、ハードディスク上に出力する。この様子は以前に示した図4と同じである。
【0160】本発明の一実施例における第三の例を以下に示す。
【0161】図13は、並列計算機上で動作する差分法による構造解析プログラムのフローチャートを示したものである。これは温度分布を差分法で解く例である。まず、温度分布を求めたい領域を正方格子に分割する。次に、支配方程式を格子間の差分方程式に変換し、それを領域全体にわたって結合する。最終的には下記のような差分方程式が得られる。
【0162】[K]{θ}={q}ここに{θ}は格子点における温度であり、{q}は格子点を通過する熱流速である。
【0163】マトリックス[K]は一般に図14に示すように非零の項が対角線に平行に現れる。温度分布をシミュレーションする場合、ひとつの格子点の自由度は1である。構造の応力、変位を求めるシミュレーションでは節点の自由度は3であったが、この場合は1なので、スパースパターンの最初の小さな成分の除去は、3×3の部分マトリックスの単位ではなく、成分そのものの大きさによって行う。したがって、プロセスの順番としては、(1)マトリックスを各対角項の大きさで各行、各列を割ることによって、スケーリングする。(2)マトリックスの成分のうち、既定値以下のものを削除したスパースパターンを作成する。この既定値は0から1の間にとる。(3)残った成分から一次のフィルインを計算し、先のスパースパターンに追加する。(4)FSAI法により、コレスキー分解における下三角マトリックスLの逆マトリックスのスパース近似(Sparse Approximate Inverse)G〜を並列計算機で並列に計算する。(5)このG〜を使って前処理付き共役勾配法を適用する。
【0164】これにより各格子点における温度が求められるので、親プロセッサは、その結果を各部分領域を担当する子プロセッサに転送する。子プロセッサは、その値から差分格子における熱流速を計算する。以上から求まった熱流速を親プロセッサは再度集めてひとつのファイルを作成し、ハードディスク上に出力する。
【0165】本発明の一実施例における第四の例を以下に示す。
【0166】平面応力状態(紙面に垂直な方向の応力が零)や平面ひずみ状態(紙面に垂直な方向のひずみが零)にある二次元形状の構造物の応力、変位をシミュレーションする場合は、フローチャートは実施例1の場合と同じ図8になる。それ以下の手順も基本的に同じであるが、三次元構造物の実施例1では節点の自由度が3であったのに対し、本実施例では2となるので、図15のようになる。すなわち、全体剛性マトリックスは2行2列ごとに分割し、サイズ2×2の部分マトリックスを作成する。部分マトリックスの削除は、判定の基準値を0から2の間にとり、ノルムがその基準値以下のものをスパースパターンから削除する。つぎに一次のフィルインの位置を計算し、スパースパターンに追加する。このスパースパターンを使って、コレスキー分解における下三角マトリックスLの逆マトリックスのスパース近似G〜を並列計算機で並列に計算する。このG〜を使って前処理付き共役勾配法を適用する。
【0167】本発明の一実施例(例として第二の例を示す)における効果を、図2と同じ187,280節点、561,840自由度の有限要素法構造解析の問題に適用した結果を図16に示す。ここでは、許容値εは0.15とした。この図から、新しいスパースパターンを使うと従来のスケーリングのみを行った場合やオリジナルのFSAI法に比べて約1/2の計算時間になっていることが分かる。
【0168】本発明を使えば、同じ並列計算機を使っても、より速いシミュレーションを実現できる。また、プロセッサ数の多い超並列計算機でシミュレーションを行うときには、通常、プロセッサ間のデータ転送回数の増加が問題となるが、本発明では前処理用行列の計算が完全に並列に計算できるので、前処理用行列の計算でプロセッサ間のデータ転送は不要であり、超並列計算機でも高速化の効果が期待できる。
【0169】また、前処理の効果が大きいところから、シェル構造のシミュレーションのように、従来の反復解法では解の収束が悪く、反復解法が適用できなかったようなものでも、本発明を適用することにより反復解法が適用できるようになる。すなわち、並列計算機でない通常のワークステーションやパソコンの計算でも、本発明を用いることにより、省メモリという利点をもっている反復解法で解の収束性を向上させることができる。
【0170】
【発明の効果】本発明により、短時間で高精度の計算を導く前処理行程を有する連立一次方程式に関する計算装置を提供することができる。
【出願人】 【識別番号】000005108
【氏名又は名称】株式会社日立製作所
【出願日】 平成12年11月10日(2000.11.10)
【代理人】 【識別番号】100075096
【弁理士】
【氏名又は名称】作田 康夫
【公開番号】 特開2002−149630(P2002−149630A)
【公開日】 平成14年5月24日(2002.5.24)
【出願番号】 特願2000−349306(P2000−349306)