JP6090602B2 - 連立１次方程式の解の算出装置及び算出方法、並びにそれに適用されるプログラム - Google Patents

Description

本発明は、連立１次方程式の解の算出装置及び算出方法、並びにそれに適用されるプログラムに関し、特に複数の演算器で並列処理をする連立１次方程式の解の算出装置及び算出方法、並びにそれに適用されるプログラムに関するものである。

近年、例えば有限要素法を用いたシミュレーション装置等のいわゆる工業数学の手法を応用した装置においては、コンピュータを用いて大規模な連立１次方程式を解く必要が生じることが多くなっている。

このような大規模連立１次方程式の求解を行う方法は、まとめれば直接法と反復法とに分類される。

直接法は、連立１次方程式の未知数を一つずつ求めていく方法であり、代表的な直接法としては、いわゆるＬＵ分解を利用したものが挙げられる。一方で、反復法は、漸化式を用いて、暫定解ベクトルを真の解へと漸近させていく方法である。直接法では、決まった計算ステップによって確実に解に到達できるのに対して、反復法では、解に到達するまでの計算量が予測できないという側面がある。

しかしながら、本発明者の検討によれば、例えば有限要素法を用いたシミュレーション方法で現れる連立１次方程式を含め、工業数学の手法を応用した分野で求解対象となる大規模の連立１次方程式Ａｘ＝ｂにおける係数行列Ａは疎行列であることが多いが、その場合、直接法が係数行列の次元数の二乗に比例するメモリを必要とするのに対して、反復法は係数行列の次元数の一乗に比例するメモリしか必要としない。従って、非常に大規模な連立方程式を扱う場合において、反復法の使用が有利であり、現実的には反復法しか使えない場合がほとんどであると考えられる。

このような、反復法として用いられるものには、いわゆるクリノフ部分空間法に属する反復法が挙げられる。クリノフ部分空間法に属する反復法としては、共役勾配法、共役残差法、双共役勾配法、二乗共役勾配法、一般化共役残差法、一般化最小残差法、及びＭＲＴＲ法等が挙げられる。かかるクリノフ部分空間法は、他の反復法と比較してはるかに少ない反復回数で解に漸近するという有意性を有する。

更に、Ｍ^−１≒Ａ^−１の成立する前処理行列Ｍ^−１を漸化式の中に取り入れることで、クリノフ部分空間法に属する反復法の必要反復回数が減る傾向があるため、前処理行列Ｍ^−１が漸化式の中に取り入れられることがある。前処理行列Ｍ^−１を漸化式の中に取り入れる手法としては、いわゆる不完全ＬＵ分解によって、Ａ≒ＬＵが成立する下三角行列Ｌと上三角行列Ｕを計算して、（ＬＵ）^−１、即ちＵ^−１Ｌ^−１を前処理行列として用いる手法が挙げられる。

ここで、不完全ＬＵ分解というのは、Ａの非零要素が存在する場所にしか、ＬとＵとが非零要素を持たないという制約の下で、Ａ≒ＬＵが成立するようＬ及びＵを求める計算手法のことである。

また、前処理行列としてＵ^−１Ｌ^−１を用いるというのは、Ｕ^−１やＬ^−１を実際に求めてそれらの積を用いることを意味しない。というのは、Ｕ及びＬは疎行列であるが、Ｕ^−１及びＬ^−１は密行列となるため、Ｕ^−１及びＬ^−１を実際に求めようとすると計算量や必要メモリ量が膨大となり、反復法を用いるメリットがなくなってしまうからである。そこで、実際には、Ｕ及びＬを計算しておき、それらを以下の手法で用いることになる。具体的には、漸化式の中で行われる操作は、前処理行列をベクトルに作用させる操作であり、これは結局、Ｌ^−１をベクトルに作用させた後、この結果にＵ^−１を作用させる操作となる。つまり、Ｌ^−１をベクトルに作用させる操作は、Ｌを係数行列とする１次方程式を前進代入法によって求める操作で実現し、Ｕ^−１をベクトルに作用させる操作は、Ｕを係数行列とする１次方程式を後退代入法によって求める操作で実現することが可能である。

ここで、本発明者の更なる検討によれば、不完全ＬＵ分解で得られたＬ及びＵを用いて、前処理行列としてＵ^−１Ｌ^−１を使用するという手法は、反復回数の削減に非常に効果的であるものの、並列計算に向かないという傾向を有する。

というのは、漸化式の中にＬを係数行列とする１次方程式を前進代入法によって求める操作が入っているが、前進代入法はベクトルの成分を上から順番に確定させていくプロセスであるため、ベクトル内の複数の場所を並列に計算することが非常に難しい傾向があるからである。併せて、漸化式の中にＵを係数行列とする１次方程式を後退代入法によって求める操作も入っているが、後退代入法はベクトルの成分を下から順番に確定させていくプロセスであるため、ベクトル内の複数の場所を並列に計算することが非常に難しい傾向があるからである。つまり、前処理行列としてＵ^−１Ｌ^−１を使用する手法は、前進代入及び後退代入というベクトル内の計算順序の制約が厳しい演算プロセスを含むため、いわゆるベクトルコンピューティングの手法による並列計算の適用が非常に困難な傾向が強い。

そのため、前処理行列としてＵ^−１Ｌ^−１を使用するという手法は、演算処理の並列化による高速化が困難であるという課題を有するものであり、特に近年、マルチコアＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）やＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）による並列計算が待望されている現状を鑑みると、これは非常に大きな解決すべき課題であるといえる。

本発明は、以上の検討を経てなされたもので、簡便な構成で、並列計算による高速な連立１次方程式の解法を実現できる連立１次方程式の解の算出装置及び算出方法、並びにそれに適用されるプログラムを提供することを目的とする。

以上の目的を達成すべく、本発明の第１の局面においては、クリノフ部分空間法に属する反復法を用いた連立１次方程式の解法を実現する連立１次方程式の解の算出装置であって、複数の演算器を備え、前記連立１次方程式の解に漸近していくベクトル列ｘ_ｋ（ｋは０を含む自然数）は、その次数に応じた複数の成分からなり、前記ベクトル列ｘ_ｋを前記解に漸近させていく反復計算において、前記ベクトル列ｘ_ｋを前記複数の成分に対応した互いに異なる複数の領域に分割すると共に、前記複数の演算器の各々に対して、前記互いに異なる複数の領域に対応した演算処理を並行して実行させるに際して、前記反復計算において用いられる前処理行列が、対角行列として与えられ、前記反復計算において、前記対角行列を前記反復計算で用いられる対応したベクトル列に作用させる際に、前記対角行列の第ｎ（ｎは自然数）成分と前記対応したベクトル列の第ｎ成分との積の演算処理を並行して実行する算出装置である。

かかる構成においては、反復計算において用いられる前処理行列が対角行列であるため、前処理行列をベクトルに作用させる計算操作において前進代入や後退代入を使用する必要性はなくなる。また、対角行列である前処理行列の対角成分を並べて作成したベクトル（以下、前処理用ベクトルという）を用意しておけば、前処理行列を被作用ベクトルに作用させるという演算プロセスは、前処理用ベクトルと被作用ベクトルとの対応する成分毎の積を計算して、それを並べたものを計算結果のベクトルとすることで実現できることになる。ここで、前処理用ベクトルと被作用ベクトルとの成分毎の積の計算は、成分毎に独立して実行することができるため、簡便なベクトルコンピューティングの手法で並列計算ができる。ここでいう「簡便なベクトルコンピューティングの手法」とは、計算対象となるベクトルを複数の領域に分割して、その相異なる領域の計算を、異なる演算器が並行して行う計算手法のことである。また、演算器とは、四則演算が可能で電力で動作する論理回路である。

前記複数の演算器が、同一のマルチコアＣＰＵの中に設けられた複数のコア同士の組み合わせ、同一のＧＰＵの中に設けられた複数のコア同士の組み合わせ、複数のＣＰＵ同士の組み合わせ、複数のＧＰＵ同士の組み合わせ、ＣＰＵ及びＧＰＵの組み合わせ、又は複数のコンピュータ同士の組み合わせであることを第２の局面とする。

かかる構成においては、共にメモリの共有が可能である同一ＣＰＵの中の複数コアを使用する構成、及び同一ＧＰＵの中の複数コアを使用する構成は、データ転送の時間を短縮し易いため、高速演算に有利であるといえる。一方で、複数のＣＰＵ及び複数のＧＰＵ、更には複数のコンピュータを用いて並列演算を行う構成も考えられる。かかる場合は、より大規模な連立方程式を扱えるようになるという利点や、より並列度を上げることができるようになるという利点がある。もちろん、これらは、必要とされる性能や、実際に確保できる演算器の種類や台数に応じて、自由度高く組み合わせることができるものである。

また、本発明は、かかる第１又は第２の局面に加え、前記反復法が、対角行列Ｍ^−１を前記前処理行列とする前処理付共役勾配法であって、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、前記ベクトル列ｘ_ｋ、ベクトル列ｒ_ｋ及びｐ_ｋ、並びにスカラー列α_ｋ及びβ_ｋを準備して、数式（数１）で示される漸化式を用いることで、前記ベクトル列ｘ_ｋを前記解ベクトルｘに漸近させていくことを第３の局面とする。

かかる構成においては、共役勾配法を用いるものであるため、一回の反復計算の中の計算量が減少するものである。

また、本発明は、かかる第１又は第２の局面に加え、前記反復法が、対角行列Ｍ^−１を前記前処理行列とする前処理付ＭＲＴＲ法であって、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、前記ベクトル列ｘ_ｋ、ベクトル列ｙ_ｋ、ｒ_ｋ、ｕ_ｋ及びｐ_ｋ、並びにスカラー列ζ_ｋ、η_ｋ及びγ_ｋを準備して、数式（数２）で示される漸化式を用いることで、前記ベクトル列ｘ_ｋを前記解ベクトルｘに漸近させていくことを第４の局面とする。

かかる構成においては、ＭＲＴＲ法を用いるものであるため、共役勾配法と比較して一回の反復計算の中の計算量が多くなるが、共役勾配法よりも少ない反復回数で解に到達できることが多くなる。また、Ａが非対称行列である場合、共役勾配法では解への漸近が保障されないのに対して、ＭＲＴＲ法であればＡが非対称行列であっても正定値行列でさえあれば、解への漸近が保障されている。更に、いわゆる丸め誤差が蓄積しにくく、収束性に優れたものといえる。

また、本発明は、かかる第１から４のいずれかの局面に加え、前記前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように前記対角行列Ｍ^−１が定められていることを第５の局面とする。

かかる構成においては、前処理行列として、行列積Ｍ^−１Ａの対角成分の絶対値が全て１の近傍、理想的には０．９以上１．１以下となるような対角行列Ｍ^−１を用いているため、前処理に必要な計算量が削減されるものである。例えば、行列積Ｍ^−１Ａの対角成分の絶対値が全て１になるような対角行列Ｍ^−１は、Ａの対角成分の逆数を対角成分として並べていけば作成されることになる。

また、本発明は、かかる第１から第４のいずれかの局面に加え、前記前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように前記対角行列Ｍ^−１が定められていることを第６の局面とする。

かかる構成においては、前処理行列として、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるような対角行列Ｍ^−１を用いているため、前処理に必要な計算量が削減されるものである。例えば、Ａを構成する行ベクトルのノルムの逆数を対角成分として並べていけば、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１になるような対角行列Ｍ^−１を作成されることになる。

また、本発明は、かかる第１から第４のいずれかの局面に加え、前記前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ_１ ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように前記対角行列Ｍ_１ ^−１を定め、行列積Ｍ_２ ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように対角行列Ｍ_２ ^−１を定め、かつ、任意の複素数ｍを用いてｍＭ_１ ^−１＋（１−ｍ）Ｍ_２ ^−１を計算し、これを前記前処理行列として用いることを第７の局面とする。

かかる構成においては、前処理行列Ｍ^−１として、行列積Ｍ_１ ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように定めた対角行列Ｍ_１ ^−１、行列積Ｍ_２ ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように定めた対角行列Ｍ_２ ^−１、及び任意の複素数ｍを用いて規定されるｍＭ_１ ^−１＋（１−ｍ）Ｍ_２ ^−１を用いるため、前処理に必要な計算量が削減されるものである。

また、本発明は、別の局面において、第１から第７のいずれかの局面における前記算出装置を使用して、前記対角行列である前記前処理行列を用いて前記ベクトル列ｘ_ｋを前記連立１次方程式の前記解に漸近させていく前記反復計算を実行させる際に、前記複数の演算器の各々に対して前記ベクトル列ｘ_ｋにおける前記次数に対応して前記並列の演算処理を実行させるプログラムである。

かかる構成のプログラムをコンピュータ上で実行させることにより、簡便な構成で、並列計算しながら高速に連立１次方程式の解が求まる。

また、本発明は、更に別の局面において、第１から第７のいずれかの局面における前記算出装置を使用して、前記対角行列である前記前処理行列を用いて前記ベクトル列ｘ_ｋを前記連立１次方程式の前記解に漸近させていく前記反復計算を実行させる際に、前記複数の演算器の各々に対して前記ベクトル列ｘ_ｋにおける前記次数に対応して前記並列の演算処理を実行させる算出方法である。

かかる構成の算出方法により、簡便な構成で、並列計算しながら高速に連立１次方程式の解が求まる。

本発明の第１の局面における構成によれば、クリノフ部分空間法に属する反復法を用いた連立１次方程式の解法を実現する連立１次方程式の解の算出装置であって、複数の演算器を備え、連立１次方程式の解に漸近していくベクトル列ｘ_ｋ（ｋは０を含む自然数）が、その次数に応じた複数の成分からなり、ベクトル列ｘ_ｋを解に漸近させていく反復計算において、ベクトル列ｘ_ｋを複数の成分に対応した互いに異なる複数の領域に分割すると共に、複数の演算器の各々に対して、互いに異なる複数の領域に対応した演算処理を並行して実行させるに際して、反復計算において用いられる前処理行列が、対角行列として与えられ、反復計算において、対角行列を反復計算で用いられる対応したベクトル列に作用させる際に、対角行列の第ｎ（ｎは自然数）成分と対応したベクトル列の第ｎ成分との積の演算処理を並行して実行するため、簡便な構成で、並列計算により連立１次方程式の解を高速に算出することができる。

また、本発明の第２の局面における構成によれば、複数の演算器が、同一のマルチコアＣＰＵの中に設けられた複数のコア同士の組み合わせ、同一のＧＰＵの中に設けられた複数のコア同士の組み合わせ、複数のＣＰＵ同士の組み合わせ、複数のＧＰＵ同士の組み合わせ、ＣＰＵ及びＧＰＵの組み合わせ、又は複数のコンピュータ同士の組み合わせであることにより、共にメモリの共有が可能である同一ＣＰＵの中の複数コアを使用する場合、及び同一ＧＰＵの中の複数コアを使用する場合には、データ転送の時間を短縮できて、高速演算を実現できる。また、複数のＣＰＵ及び複数のＧＰＵ、更には複数のコンピュータを用いて並列演算を行う場合には、より大規模な連立方程式を扱うことができると共に、より並列度を上げることができる。

また、本発明の第３局面における構成によれば、反復法が、対角行列Ｍ^−１を前処理行列とする前処理付共役勾配法であって、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、ベクトル列ｘ_ｋ、ベクトル列ｒ_ｋ及びｐ_ｋ、並びにスカラー列α_ｋ及びβ_ｋを準備して、前述の数式（数１）で示される漸化式を用いることで、ベクトル列ｘ_ｋを解ベクトルｘに漸近させていくことにより、一回の反復計算の中の計算量を減少することができる。

また、本発明の第４の局面における構成によれば、反復法が、対角行列Ｍ^−１を前処理行列とする前処理付ＭＲＴＲ法であって、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、ベクトル列ｘ_ｋ、ベクトル列ｙ_ｋ、ｒ_ｋ、ｕ_ｋ及びｐ_ｋ、並びにスカラー列ζ_ｋ、η_ｋ及びγ_ｋを準備して、前述した数式（数２）で示される漸化式を用いることで、ベクトル列ｘ_ｋを解ベクトルｘに漸近させていくことにより、共役勾配法と比較して一回の反復計算の中の計算量が多くなるが、共役勾配法よりも少ない反復回数で解に到達できる。また、Ａが非対称行列である場合、共役勾配法では解への漸近が保障されないのに対して、ＭＲＴＲ法であればＡが非対称行列であっても正定値行列でさえあれば、解への漸近を保障することができる。更に、いわゆる丸め誤差が蓄積しにくく、解への収束性を向上することができる。

また、本発明の第５の局面における構成によれば、前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように対角行列Ｍ^−１が定められていることにより、前処理に必要な計算量を削減することができる。

また、本発明の第６の局面における構成によれば、前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように対角行列Ｍ^−１が定められていることにより、前処理に必要な計算量を削減することができる。

また、本発明の別な局面における構成によれば、第１から第７のいずれかの局面における算出装置を使用して、対角行列である前処理行列を用いてベクトル列ｘ_ｋを連立１次方程式の解に漸近させていく反復計算を実行させる際に、複数の演算器の各々に対してベクトル列ｘ_ｋにおける次数に対応して並列の演算処理を実行させるプログラムであるため、かかる構成のプログラムをコンピュータ上で実行させることにより、簡便な構成で、並列計算しながら高速に連立１次方程式の解を算出することができる。

また、本発明の更に別の局面における構成によれば、第１から第７のいずれかの局面における算出装置を使用して、対角行列である前処理行列を用いてベクトル列ｘ_ｋを連立１次方程式の解に漸近させていく反復計算を実行させる際に、複数の演算器の各々に対してベクトル列ｘ_ｋにおける次数に対応して並列の演算処理を実行させる算出方法であることにより、簡便な構成で、並列計算しながら高速に連立１次方程式の解を算出することができる。

本発明における実施形態の連立１次方程式の解の算出装置の構成を示すブロック図である。 本実施形態の変形例における連立１次方程式の解の算出装置の構成を示すブロック図である。 本実施形態の別の変形例における連立１次方程式の解の算出装置の構成を示すブロック図である。 本実施形態の更に別の変形例における連立１次方程式の解の算出装置の構成を示すブロック図である。 本実施形態における連立１次方程式の解の算出装置を用いて行った有限要素法によるシミュレーションの求解演算結果を、比較例と共に示す図であり、横軸が反復回数で、縦軸が残差である。 本実施形態における連立１次方程式の解の算出装置を用いて行った有限要素法によるシミュレーションの別の求解演算結果を、比較例と共に示す図であり、横軸が反復回数で、縦軸が残差である。

以下、図面を適宜参照して、本発明の実施形態における連立１次方程式の解の算出装置及び算出方法、並びにそれに適用されるプログラムにつき、詳細に説明する。

［算出装置の構成］
まず、本実施形態における連立１次方程式の解の算出装置につき、その各種変形例を含め、図１から図４を参照して、詳細に説明する。

図１は、本実施形態の連立１次方程式の解の算出装置の構成を示すブロック図である。また、図２は、本実施形態の変形例における連立１次方程式の解の算出装置の構成を示すブロック図である。図３は、本実施形態の別の変形例における連立１次方程式の解の算出装置の構成を示すブロック図である。また、図４は、本実施形態の更に別の変形例における連立１次方程式の解の算出装置の構成を示すブロック図である。

図１に示すように、本実施形態における連立１次方程式の解の算出装置は、コンピュータ１０を有する。

コンピュータ１０は、コア２１〜２４を内蔵するＣＰＵ２０と、ＣＰＵ２０に電気的に接続する通信バス３０と、通信バス３０に電気的に接続するメモリ４０と、を備える。

コア２１〜２４は、各々、連立１次方程式の解の算出を実行するにあたり、並列の演算処理を実行自在な演算器である。ここでは、複数のコアとして、商用品の一例として計４個のコア２１〜２４を示すが、その総数は限定的なものではなく、原理的には複数個設けられていればよい。

メモリ４０は、コア２１〜２４に連立１次方程式の解の算出を実行させるにあたり、並列の演算処理を実行させるアルゴリズムを有するプログラムを格納しており、通信バス３０を介して、ＣＰＵ２０内のコア２１〜２４に電気的に接続している。ここでは、メモリ４０を１個のみ示すが、その総数は限定的なものではなく、複数個設けられていてもよい。また、メモリ４０は、典型的には不揮発性のメモリである。

このように、本実施形態におけるコア２１〜２４を内蔵するＣＰＵ２０を備えたコンピュータ１０を用いた算出装置では、同一ＣＰＵ２０の中の複数のコア２１〜２４を複数の演算器として用いて連立１次方程式の解の算出を並行して実行させることが可能になると共に、メモリ４０の共有が可能となり、通信バス３０を介したデータ転送の時間を短縮することができるため、高速演算に有利であるといえる。
［各種変形例における算出装置の構成］

ここで、本実施形態における連立１次方程式の解の算出装置においては、連立１次方程式の解の算出を実行する際に、その複数の演算器に並列の演算処理を行わせるものであるから、種々の変形例が考えられる。

かかる連立１次方程式の解の算出装置の変形例としては、図２から図４に示すものが挙げられる。

まず、図２に示すように、本実施形態の変形例における連立１次方程式の解の算出装置は、コンピュータ５０を有する。

コンピュータ５０は、コア６１〜７６を内蔵するＧＰＵ６０と、ＧＰＵ６０に電気的に接続する通信バス８０と、通信バス８０に電気的に接続するメモリ９０と、を備える。

ＧＰＵ６０においても、コア６１〜７６は、各々、連立１次方程式の解の算出を実行するにあたり、並列の演算処理を実行自在な演算器である。ここでは、複数のコアとして、商用品の一例として計１６個のコア６１〜７６を示すが、その総数は限定的なものではなく、原理的には複数個設けられていればよい。

メモリ９０は、コア６１〜７６に連立１次方程式の解の算出を実行させるにあたり、並列の演算処理を実行させるアルゴリズムを有するプログラムを格納しており、通信バス８０を介して、ＧＰＵ６０内のコア６１〜７６に電気的に接続している。ここでは、メモリ９０を１個のみ示すが、その総数は限定的なものではなく、複数個設けられていてもよい。また、メモリ９０は、典型的には不揮発性のメモリである。

このように、本変形例におけるコア６１〜７６を内蔵するＧＰＵ６０を備えたコンピュータ５０を用いた算出装置においても、同一ＧＰＵ６０の中の複数のコア６１〜７６を複数の演算器として用いて連立１次方程式の解の算出を並行して実行させることが可能になると共に、メモリ９０の共有が可能となり、通信バス８０を介したデータ転送の時間を短縮することができるため、高速演算に有利であるといえる。更に、一般に、ＧＰＵが内蔵するコアの数は、ＣＰＵが内蔵するコアの数よりも多いため、本変形例におけるＧＰＵ６０を備えた算出装置は、図１に示す本実施形態におけるＣＰＵ２０を備えた算出装置よりも高価にはなるが、並列して行い得る演算処理の数をより増やすことが可能となり、高速演算にはより有利であるといえる。

次に、図３に示すように、本実施形態の別の変形例における連立１次方程式の解の算出装置は、コンピュータ１１０、１２０、１３０及び１４０を連携させたコンピュータ群１００を有する。

コンピュータ１１０は、図示を省略する複数のコアを内蔵するＣＰＵ１１１と、ＣＰＵ１１１に、図示を省略する通信バスを介して、電気的に接続するメモリ１１２と、を備える。ＣＰＵ１１１の構成は、図１に示す本実施形態におけるＣＰＵ２０のものと同様である。但し、本変形例では、複数のコンピュータを演算器として用いるものであるため、ＣＰＵ１１１が内蔵するコアの数は１個であってもよい。

コンピュータ１２０は、図示を省略する複数のコアを内蔵するＣＰＵ１２１と、ＣＰＵ１２１に、図示を省略する通信バスを介して、電気的に接続するメモリ１２２と、を備える。かかるコンピュータ１２０の構成は、コンピュータ１１０のものと同様である。

コンピュータ１３０は、図示を省略する複数のコアを内蔵するＧＰＵ１３１と、ＧＰＵ１３１に、図示を省略する通信バスを介して、電気的に接続するメモリ１３２と、を備える。ＧＰＵ１３１の構成は、図２に示す変形例におけるＧＰＵ６０のものと同様である。かかる構成は、ＧＰＵ１４１と、ＧＰＵ１４１に、図示を省略する通信バスを介して、電気的に接続するメモリ１４２と、を備えたコンピュータ１４０についても同様である。

これらのコンピュータ１１０、１２０、１３０及び１４０は、各々、連立１次方程式の解の算出を実行するにあたり、並列の演算処理を実行自在な演算器である。ここでは、複数のコンピュータとして、一例として計４個のコンピュータ１１０、１２０、１３０及び１４０を示すが、その総数は限定的なものではなく、原理的には複数個設けられていればよい。また、かかる複数のコンピュータは、ＣＰＵのみを内蔵するコンピュータのみから構成されていてもよいし、ＧＰＵのみを内蔵するコンピュータのみから構成されていてもよい。

メモリ１１２及び１２２の構成は、図１に示す本実施形態におけるメモリ４０のものと同様であり、メモリ１３２及び１４２の構成は、図２に示す変形例におけるメモリ９０のものと同様である。ここで、コンピュータ１１０、１２０、１３０及び１４０に連立１次方程式の解の算出を実行させるにあたり、並列の演算処理を実行させるアルゴリズムを有するプログラムを格納しておくメモリは、典型的にはメモリ１１２、１２２、１３２及び１４２のいずれかに設定しておけばよい。かかる場合には、そのプログラムを格納したコンピュータが、マスター側となって残余のコンピュータをスレーブ側として制御し、コンピュータ１１０、１２０、１３０及び１４０が協働して、それら全体で並列の演算処理を実行すればよい。

コンピュータ１１０、１２０、１３０及び１４０は、通信網１５０に電気的に接続され、それを介して、互いに必要な信号を送受信する。通信網１５０は、典型的には、ＬＡＮ（Ｌｏｃａｌ Ａｒｅａ Ｎｅｔｗｏｒｋ）、インターネット網、又はそれらの組み合わせから構成される。通信網１５０に図示を省略するメモリが電気的に接続される場合には、そのメモリ内に、連立１次方程式の解を算出するために並列の演算処理を実行させるアルゴリズムを有するプログラムを格納しておいてもよい。

このように、本変形例におけるコンピュータ１１０、１２０、１３０及び１４０を連携させたコンピュータ群１００を用いた算出装置においては、通信網１５０に電気的に接続した複数のコンピュータ１１０、１２０、１３０及び１４０を複数の演算器として用いて連立１次方程式の解の算出を並行して実行させることが可能になる。更に、その各演算器の演算処理能力が格段に増強されると共に、より拡張的に演算器の並列度を上げることが可能になるため、より大規模な連立１次方程式の解の算出を扱うことが可能になる。

次に、図４に示すように、本実施形態の更に別の変形例における連立１次方程式の解の算出装置は、コンピュータ２００を有する。

コンピュータ２００は、ＣＰＵ２１１及び２２１と、ＧＰＵ２３１及び２４１と、ＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１に電気的に接続する通信バス２５０と、通信バス２５０に電気的に接続するメモリ２６０と、を備える。

ＣＰＵ２１１及び２２１の構成は、図１に示す本実施形態におけるＣＰＵ２０のものと同様である。但し、本変形例では、複数のＣＰＵ等を演算器として用いるものであるため、ＣＰＵ２１１及び２２１が内蔵するコアの数は１個であってもよい。

ＧＰＵ２３１及び２４１の構成は、図２に示す変形例におけるＧＰＵ６０のものと同様である。

コンピュータ２００においては、ＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１は、各々、連立１次方程式の解の算出を実行するにあたり、並列の演算処理を実行自在な演算器である。ここでは、複数の演算器として、一例として各々２個のＣＰＵ及びＧＰＵを示すが、その個数は限定的なものではなく、原理的にはトータルで複数個設けられていればよい。また、かかるコンピュータ２００は、ＣＰＵのみを内蔵していてもよいし、ＧＰＵのみを内蔵していてもよい。

メモリ２６０は、ＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１に連立１次方程式の解の算出を実行させるにあたり、並列の演算処理を実行させるアルゴリズムを有するプログラムを格納しており、通信バス２５０を介して、ＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１に電気的に接続している。かかるメモリ２６０の構成は、図１に示す本実施形態におけるメモリ４０、又は図２に示す変形例におけるメモリ９０と同様である。

このように、本変形例における同一のコンピュータ２００の中のＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１を内蔵したコンピュータ２００を用いた算出装置においては、複数のＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１を複数の演算器として用いて連立１次方程式の解の算出を並行して実行させることが可能になる。更に、その各演算器の演算処理能力が、外部の通信網を利用することがない態様で格段に増強されることが可能になるため、より大規模な連立１次方程式の解の算出を扱うことが可能になる。但し、本変形例では、ＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１を動作させるＯＳ（Ｏｐｅｒａｔｉｏｎ Ｓｙｓｔｅｍ）の構成が煩雑なるため、商用品としての選択自由度には限界がある。

［算出装置に適用されるプログラムの内容］
次に、各種変形例を含む本実施形態における連立１次方程式の解の算出装置に適用されるプログラムにつき、詳細に説明する。かかるプログラムは、各種変形例を含む本実施形態における連立１次方程式の解の算出装置のメモリに格納された後、必要に応じて読み出されて実行されるものである。また、本実施形態における連立１次方程式の解の算出方法は、かかるプログラムを実行することにより、各種変形例を含む本実施形態における連立１次方程式の解の算出装置を用いてなされるものである。

本実施形態における連立１次方程式の解の算出装置が、以上説明した同一ＣＰＵ２０の中の複数のコア２１〜２４、同一ＧＰＵ６０の中の複数のコア６１〜７６、通信網１５０に電気的に接続した複数のコンピュータ１１０、１２０、１３０及び１４０、及び同一のコンピュータ２００の中のＣＰＵ２１１及び２２１並びにＧＰＵ２３１及び２４１を複数の演算器として用いて、連立１次方程式の解の算出を並行して実行させるための算出処理で用いる解法は、クリノフ部分空間法に属する反復法である。

具体的には、クリノフ部分空間法に属する反復法においては、ベクトル列ｘ_ｋ（ｋは０を含む自然数）を求解方程式である連立１次方程式の解に漸近させていく反復計算において、前処理行列が用いられる。

ここで、各種変形例を含む本実施形態では、複数の演算器に連立１次方程式の解の算出を並行して実行させるために、反復計算において用いられる前処理行列が、対角行列に設定される。

このように、前処理行列が対角行列に設定されることにより、ベクトル列ｘ_ｋを、その次数に応じた複数の成分に対応して、互いに異なる複数の領域に分割すると共に、複数の演算器の各々に対して、互いに異なる複数の領域に対応した演算処理を並行して実行させることが可能となる。

以下、求解方程式である連立１次方程式をＡｘ＝ｂと表現して、前処理行列及びそれを用いたアルゴリズムにつき説明を進めることとする。ここで、ｘが求解対象となる未知ベクトルで、Ａは既知（値の確定した）の係数行列、及びｂは既知（値の確定した）のベクトルである。

前処理行列として使用できるＭ^−１としては、以下のようなものを挙げることができる。

まず、行列積Ｍ^−１Ａの対角成分の絶対値が全て１の近傍、理想的には０．９以上１．１以下となるような対角行列Ｍ^−１である。例えば、Ａの対角成分の逆数を対角成分として並べていけば、行列積Ｍ^−１Ａの対角成分の絶対値が全て１になるような対角行列Ｍ^−１を作成することができる。

次に、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１の近傍、理想的には０．９以上１．１以下となるような対角行列Ｍ^−１である。例えば、Ａを構成する行ベクトルのノルムの逆数を対角成分として並べていけば、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１になるような対角行列Ｍ^−１を作成することができる。

次に前処理行列として使用できる別のＭ^−１として、先に挙げた２つの例から合成される、以下のようなものを挙げることができる。即ち、行列積Ｍ_１ ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように対角行列Ｍ_１ ^−１を定め、行列積Ｍ_２ ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように対角行列Ｍ_２ ^−１を定め、任意の複素数ｍを用いてＭ^−１＝ｍＭ_１ ^−１＋（１−ｍ）Ｍ_２ ^−１として作成したＭ^−１である。

以上、前処理行列として使用できる対角行列Ｍ^−１の例をいくつか挙げたが、いずれが優れているかに関しては、与えられた１次方程式に依存するため一概にはいえない。しかし、いずれを用いた場合にも、前処理を全く行わないのと比較すれば、かなり大きな反復回数削減効果を得ることができるものである。

但し、前処理行列として不完全ＬＵ分解の結果を用いる従来の手法と比較すると、各種変形例を含む本実施形態における対角行列を前処理行列とする手法は、解に到達するまでの反復回数が多く必要となることが多い。即ち、反復回数の削減という観点では、対角行列を前処理行列とする手法は、不完全ＬＵ分解を用いる従来の手法に劣ることが多い。

しかしながら、不完全ＬＵ分解を用いる従来の手法と比較して、対角行列を前処理行列とする手法では、前処理に必要な計算量が削減されるため、反復回数が多くなっても解に到達するまでの計算量の合計はそれほど多くならないか、又は削減される。更に、対角行列を前処理行列とする手法によって反復計算の並列化が効率的にできるようになるため、解に到達するまでの計算時間そのものは大幅に短縮されることになる。

対角行列を前処理行列として使用する反復法のアルゴリズムとしては、原理的には、クリノフ部分空間法に属する手法を任意に選択して、ベクトル列ｘ_ｋを求解方程式である連立１次方程式Ａｘ＝ｂの解に漸近させていく反復計算を実行する際に、ベクトル列ｘ_ｋをその次数に応じた複数の成分に対応した互いに異なる複数の領域に分割すると共に、複数の演算器の各々に対して、かかる互いに異なる複数の領域に対応した演算処理を並行して実行させればよい。

ここで、クリノフ部分空間法に属する手法としては、共役勾配法、共役残差法、双共役勾配法、二乗共役勾配法、一般化共役残差法、一般化最小残差法、及びＭＲＴＲ法等をその例として挙げることができる。

その中で、特に、一回の反復計算の中の計算量が少ない共役勾配法は好適な手法として挙げることができる。共役勾配法は、ベクトル列ｘ_ｋ、ｒ_ｋ、ｐ_ｋ及びスカラー列α_ｋ、β_ｋを準備して、数式（数３）で示される漸化式を用いてベクトル列ｘ_ｋを解ベクトルに漸近させていく手法である。ここで、数式（数３）において、前処理行列としての対角行列Ｍ^−１をベクトル列に作用させる際には、ｘ_ｋ＋１の式の右辺のα_ｋｐ_ｋの項がα_ｋ等の式から計算されるため、例えば、η_ｋを計算する際のＭ^−１ｒ_ｋに着目するときにベクトル列ｒ_ｋの次数がｎ（自然数）であるとすれば、Ｍ^−１ｒ_ｋの第１成分は、Ｍ^−１の対角第１成分とｒ_ｋの第１成分との積、Ｍ^−１ｒ_ｋの第２成分は、Ｍ^−１の対角第２成分とｒ_ｋの第２成分との積、Ｍ^−１ｒ_ｋの第３成分は、Ｍ^−１の対角第３成分とｒ_ｋの第３成分との積…Ｍ^−１ｒ_ｋの第ｎ成分は、Ｍ^−１の対角第ｎ成分とｒ_ｋの第ｎ成分との積として、求められる。また、例えば、複数の演算器の数がｎであれば、かかるｎ個の積の計算は、ｎ個の演算器の各々が１つの積の計算を実行することが最も並列化の高い演算処理となる。また、この際、例えば、複数の演算器の数がｎ個よりも少なければ、同一の演算ルーチンのループの中で、複数の演算器の内で複数回の積の演算をするものが出てくることになる。なお、数式（数３）において、ｘ_０の値は、ベクトル列の初期値として与えられる。

また、使用する反復法のアルゴリズムの別の例としてＭＲＴＲ法を挙げることもできる。ＭＲＴＲ法は、ベクトル列ｘ_ｋ、ｙ_ｋ、ｒ_ｋ、ｕ_ｋ、ｚ_ｋ、ｐ_ｋ及びスカラー列ζ_ｋ、η_ｋ、γ_ｋを準備して、数式（数４）で示される漸化式を用いて、ベクトル列ｘ_ｋを解ベクトルに漸近させていく手法である。ここで、数式（数４）において、前処理行列としての対角行列Ｍ^−１をベクトル列に作用させる際には、ｘ_ｋ＋１の式の右辺のζ_ｋｐ_ｋの項がＭ^−１Ａｕ_ｋ等の式から計算されるため、例えば、Ｍ^−１Ａｕ_ｋに着目するときに行列Ａがｎ次の正方行列でベクトル列ｕ_ｋの次数がｎであるとすれば（ｎは自然数）、Ｍ^−１Ａｕ_ｋの第１成分は、Ｍ^−１の対角第１成分とＡｕ_ｋ第１成分との積、Ｍ^−１Ａｕ_ｋの第２成分は、Ｍ^−１の対角第２成分とＡｕ_ｋの第２成分との積、Ｍ^−１Ａｕ_ｋの第３成分は、Ｍ^−１の対角第３成分とＡｕ_ｋの第３成分との積…Ｍ^−１Ａｕ_ｋの第ｎ成分は、Ｍ^−１の対角第ｎ成分とＡｕ_ｋの第ｎ成分との積として、求められる。また、例えば、複数の演算器の数がｎであれば、かかるｎ個の積の計算は、ｎ個の演算器の各々が１つの積の計算を実行することが最も並列化の高い演算処理となる。また、この際、例えば、複数の演算器の数がｎ個よりも少なければ、同一の演算ルーチンのループの中で、複数の演算器の内で複数回の積の演算をするものが出てくることになる。なお、数式（数４）において、ｘ_０及びｒ_０の値は、ベクトル列の初期値として与えられる。

かかるＭＲＴＲ法は、共役勾配法と比較して一回の反復計算の中の計算量が多いが、共役勾配法よりも少ない反復回数で解に到達できることが多い。また、Ａが非対称行列である場合、共役勾配法では解への漸近が保障されないのに対して、ＭＲＴＲ法であればＡが非対称行列であっても正定値行列でさえあれば、解への漸近が保障されている。共役勾配法とＭＲＴＲ法とのいずれが有利であるかは、演算処理をしようとする対象に依存する。

数式（数４）で示される漸化式を用いた場合にも、一回の反復計算の中の計算量を減少させることはできる。しかしながら、各種変形例を含む本実施形態においてＭＲＴＲ法を使用するにあたっては、数式（数５）で示される漸化式を用いるのが更に好ましい。ここで、数式（数５）において、前処理行列としての対角行列Ｍ^−１をベクトル列に作用させる際には、ｘ_ｋ＋１の式の右辺のζ_ｋｐ_ｋの項がＭ^−１Ａｕ_ｋ等の式から計算されるため、例えば、Ｍ^−１Ａｕ_ｋに着目するときに行列Ａがｎ次の正方行列でベクトル列ｕ_ｋの次数がｎであるとすれば（ｎは自然数）、Ｍ^−１Ａｕ_ｋの第１成分は、Ｍ^−１の対角第１成分とＡｕ_ｋ第１成分との積、Ｍ^−１Ａｕ_ｋの第２成分は、Ｍ^−１の対角第２成分とＡｕ_ｋの第２成分との積、Ｍ^−１Ａｕ_ｋの第３成分は、Ｍ^−１の対角第３成分とＡｕ_ｋの第３成分との積…Ｍ^−１Ａｕ_ｋの第ｎ成分は、Ｍ^−１の対角第ｎ成分とＡｕ_ｋの第ｎ成分との積として、求められる。また、例えば、複数の演算器の数がｎであれば、かかるｎ個の積の計算は、ｎ個の演算器の各々が１つの積の計算を実行することが最も並列化の高い演算処理となる。また、この際、例えば、複数の演算器の数がｎ個よりも少なければ、同一の演算ルーチンのループの中で、複数の演算器の内で複数回の積の演算をするものが出てくることになる。なお、数式（数５）において、ｘ_０及びｒ_０の値は、ベクトル列の初期値として与えられる。

ここで、数式（数４）で示される漸化式と比較すると、数式（数５）で示される漸化式にはベクトル列ｚ_ｋが含まれないと共に、ｕ_ｋ＋１を求める式が異なる。ベクトル列ｚ_ｋ以外の部分に関しては、これらの漸化式は数学的に同値であって、一方から他方を互いに演繹できるが、コンピュータを用いて実際の数値計算に適用した場合には、いわゆる丸め誤差の蓄積に対するふるまいが異なる。

結論としては、数式（数５）で示される漸化式を用いた方が、丸め誤差が蓄積しにくく、収束性が優れる。両者を比較すると、数式（数５）で示される漸化式は、ベクトル列ｕ_ｋを毎回本来の定義式ｕ_ｋ＝Ｍ^−１ｒ_ｋで再定義するのに対して、数式（数４）で示される漸化式はそれを避けて単純なｚ_ｋ＋１との１次結合によってｕ_ｋ＋１を演繹している。結局、数式（数４）で示される漸化式はＭ^−１の作用を一回避けるかわりに、反復計算による丸め誤差の蓄積を甘受しているといえる。Ｍ^−１の作用、つまり前処理として、不完全ＬＵ分解によって得たＬ、Ｕを使用する場合には、Ｍ^−１の作用の計算量が多いため、最初の漸化式を用いてこれを１回避けることにメリットがあるが、各種変形例を含む本実施形態のように、Ｍ^−１の作用が前処理ベクトルとの間の単純な成分毎積計算である場合には、わざわざこれを１回避けることと引き換えに、丸め誤差の蓄積を甘受するメリットはない。従って、各種変形例を含む本実施形態の中でＭＲＴＲ法を使用する場合には、数式（数５）で示される漸化式を用いることが好ましいことになる。
（算出例１）
以下、本実施形態における連立１次方程式の解の算出装置を用いて行った有限要素法によるシミュレーションの求解演算結果について説明する。

図５は、本実施形態における連立１次方程式の解の算出装置を用いて行った有限要素法によるシミュレーションの求解演算結果を、比較例と共に示す図であり、横軸が反復回数で、縦軸が残差である。

本算出例において、残差は、｜ｂ−Ａｘ_ｋ｜／｜ｂ｜で定義される量であり、ｘ_ｋがどの程度解から乖離しているかを示している。一般的には、残差の大きさが１×１０^−８程度になったところで計算を打ち切る。反復計算のアルゴリズムとしては、共役勾配法を用いた。

図５においては、前処理を行わなかった場合のデータをＡ、前処理行列として不完全ＬＵ分解から得られたＵ^−１Ｌ^−１を使用した場合のデータをＢ、数式（数３）の漸化式を用いて、前処理行列として行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１になるように定められた対角行列Ｍ^−１を使用した場合のデータをＣとして、各々示している。

図５に示すように、データＢと比較するとデータＣは、解に到達するまでに２．５倍程度の反復回数を必要としている。

しかしながら、前処理行列をベクトルに作用させるときの計算量が大きく異なるため、データＢと比較すると、データＣは、１回の反復計算の中の計算量が１／２．５程度になっている。従って、解に到達するまでに必要なデータＣの計算量は、データＢのものと実質的に同じである。

ここで、前処理行列をベクトルに作用させる計算ステップにおいて、データＢを得る場合に対応する計算アルゴリズムでは並列化が困難であるのに対して、データＣを得る場合に対応する計算アルゴリズムでは高効率の並列化が容易である。よって、並列計算を行った状態では、データＢを得る場合に比較して、データＣを得る場合には、はるかに計算時間が短くなる。

一方で、データＡで示す前処理なしの場合には、解に到達するまでの反復回数が多くなり過ぎる。もちろん、前処理なしにすれば並列計算が容易であるし、一回の反復計算の中での計算量も少なくなるが、データＣの場合に比較すると、一回の反復計算の中での計算量は９割程度（１割減る程度）であり、並列計算における並列化効率もデータＣの場合を得る場合と変わらない。

よって、総合的に考えれば、以上の３者の中では、反復計算のアルゴリズムが共役勾配法の場合に、数式（数３）の漸化式を用いて、前処理行列として行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１になるように定められた対角行列Ｍ^−１を使用し、かつ、前処理行列をベクトルに作用させる計算ステップで並列計算を実行するデータＣに対応した場合が、３者の中では最も高速に計算ができる手法といえる。
（算出例２）
次に、本実施形態における連立１次方程式の解の算出装置を用いて行った有限要素法によるシミュレーションの別の求解演算結果について説明する。

図６は、本実施形態における連立１次方程式の解の算出装置を用いて行った有限要素法によるシミュレーションの別の求解演算結果を、比較例と共に示す図であり、横軸が反復回数で、縦軸が残差である。

本算出例においては、反復計算のアルゴリズムとしてはＭＲＴＲ法を用いた。なお、残差は、算出例１と同じく、｜ｂ−Ａｘ_ｋ｜／｜ｂ｜で定義される量である。

図６において、比較例として示したデータは、前処理行列として不完全ＬＵ分解から得られたＵ^−１Ｌ^−１を使用した場合の結果である。また、算出例及び別の算出例として示した各々のデータは、いずれも前処理行列として行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１になるように定められた対角行列Ｍ^−１を使用した場合であるが、使用している漸化式が異なり、前者では、数式（数４）で前述した漸化式を用い、後者では、数式（数５）で前述した漸化式を用いた。

図６に示すように、算出例における反復回数は、比較例における反復回数よりも減少されているが、算出例では、図中で解に到達していない。一方で、別の算出例におけるものは、反復回数が一番少なくなり、図中で解に実質的に到達していると評価できる。

よって、反復計算のアルゴリズムとしてＭＲＴＲ法を用いた場合に、数式（数５）の漸化式を用いて、前処理行列として行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て１になるように定められた対角行列Ｍ^−１を使用し、かつ、前処理行列をベクトルに作用させる計算ステップで並列計算を実行する場合が、３者の中では最も高速に計算ができる手法といえる。

以上の各種変形例を含む本実施形態における連立１次方程式の解の算出装置、それに用いられる算出方法、及びそれに適用されるプログラムの構成によれば、簡便な構成で、並列計算により連立１次方程式の解を高速に算出することができる。

なお、本発明は、構成要素の形状、配置、個数等は前述の実施形態に限定されるものではなく、かかる構成要素を同等の作用効果を奏するものに適宜置換する等、発明の要旨を逸脱しない範囲で適宜変更可能であることはもちろんである。

以上のように、本発明においては、簡便な構成で、並列計算による高速な連立１次方程式の解法を実現できる連立１次方程式の解の算出装置及び算出方法、並びにそれに適用するプログラムを提供することができるものであるため、その汎用普遍的な性格から広範にシミュレーション装置等の工業数学を応用した分野に適用され得るものと期待される。

Claims (9)

1. クリノフ部分空間法に属する反復法を用いた連立１次方程式の解法を実現する連立１次方程式の解の算出装置であって、
   複数の演算器を備え、
   前記連立１次方程式の解に漸近していくベクトル列ｘ_ｋ（ｋは０を含む自然数）は、その次数に応じた複数の成分からなり、
   前記ベクトル列ｘ_ｋを前記解に漸近させていく反復計算において、前記ベクトル列ｘ_ｋを前記複数の成分に対応した互いに異なる複数の領域に分割すると共に、前記複数の演算器の各々に対して、前記互いに異なる複数の領域に対応した演算処理を並行して実行させるに際して、
   前記反復計算において用いられる前処理行列が、対角行列として与えられ、
   前記反復計算において、前記対角行列を前記反復計算で用いられる対応したベクトル列に作用させる際に、前記対角行列の第ｎ（ｎは自然数）成分と前記対応したベクトル列の第ｎ成分との積の演算処理を並行して実行する算出装置。
2. 前記複数の演算器が、同一のマルチコアＣＰＵの中に設けられた複数のコア同士の組み合わせ、同一のＧＰＵの中に設けられた複数のコア同士の組み合わせ、複数のＣＰＵ同士の組み合わせ、複数のＧＰＵ同士の組み合わせ、ＣＰＵ及びＧＰＵの組み合わせ、又は複数のコンピュータ同士の組み合わせである請求項１に記載の算出装置。
3. 前記反復法が、対角行列Ｍ ^−１ を前記前処理行列とする前処理付共役勾配法であって、
   既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、前記ベクトル列ｘｋ、ベクトル列ｒｋ及びｐｋ、並びにスカラー列αｋ及びβｋを準備して、数式（数１）で示される漸化式を用いることで、前記ベクトル列ｘｋを前記解ベクトルｘに漸近させていく請求項１又は２に記載の算出装置。
   

   
4. 前記反復法が、対角行列Ｍ^−１を前記前処理行列とする前処理付ＭＲＴＲ法であって、
   既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、前記ベクトル列ｘ_ｋ、ベクトル列ｙ_ｋ、ｒ_ｋ、ｕ_ｋ及びｐ_ｋ、並びにスカラー列ζ_ｋ、η_ｋ及びγ_ｋを準備して、数式（数２）で示される漸化式を用いることで、前記ベクトル列ｘ_ｋを前記解ベクトルｘに漸近させていく請求項１又は２に記載の算出装置。
   

   
5. 前記前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように前記対角行列Ｍ^−１が定められている請求項１から４のいずれかに記載の算出装置。
6. 前記前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように前記対角行列Ｍ^−１が定められている請求項１から４のいずれかに記載の算出装置。
7. 前記前処理行列を対角行列Ｍ^−１として、既知の行列Ａ、既知のベクトルｂからなる連立１次方程式Ａｘ＝ｂの解ベクトルｘを求めるときに、行列積Ｍ_１ ^−１Ａの対角成分の絶対値が全て０．９以上１．１以下となるように前記対角行列Ｍ_１ ^−１を定め、行列積Ｍ_２ ^−１Ａを構成する行ベクトルのノルムが全て０．９以上１．１以下となるように対角行列Ｍ_２ ^−１を定め、かつ、任意の複素数ｍを用いてｍＭ_１ ^−１＋（１−ｍ）Ｍ_２ ^−１を計算し、これを前記前処理行列として用いる請求項１から４のいずれかに記載の算出装置。
8. 前記請求項１から７のいずれかに記載の前記算出装置を使用して、前記対角行列である前記前処理行列を用いて前記ベクトル列ｘ_ｋを前記連立１次方程式の前記解に漸近させていく前記反復計算を実行させる際に、前記複数の演算器の各々に対して、前記ベクトル列ｘ_ｋにおける前記互いに異なる複数の領域に対応した演算処理を、並行して実行させるプログラム。
9. 前記請求項１から７のいずれかに記載の前記算出装置を使用して、前記対角行列である前記前処理行列を用いて前記ベクトル列ｘ_ｋを前記連立１次方程式の前記解に漸近させていく前記反復計算を実行させる際に、前記複数の演算器の各々に対して、前記ベクトル列ｘ_ｋにおける前記互いに異なる複数の領域に対応した演算処理を、並行して実行させる算出方法。

Images