(57)【特許請求の範囲】
【請求項１】
対称な非零構造を持つＮ（Ｎ：２以上の整数）次元の三重対角行列Ａを係数行列とし、ベクトルｂを定数項とする連立一次方程式Ax=bを、ｐ（ｐ：２以上の整数）台のプロセッサとメモリと入力装置と出力装置とを備えた並列計算機システムにＬＵ分解を実行させることにより解く求解方法であって、
（１）前記係数行列Ａ、右辺ベクトルbを前記入力装置により前記ｐ台のプロセッサに入力するステップと、
（２）前記ｐ台のプロセッサの少なくとも１台を用いて、下記条件１〜３を満たすように１〜Ｎまでの数の集合を互いに重なりのない２p＋１個の集合に分割し、１．それぞれが３以上、Ｎ−２以下であり、かつ互いに４以上離れたp-1個の増加する整数の列IC１、IC２、・・・、ICｐ−１を選んで集合IC＝{IC１、IC２、・・・、ICｐ−１}とし、
２．IBをIB１={IC１-1}、IBｐ={ICｐ−１+1}、２≦k≦p-1のときIBｋ={ICｋ−１+1、ICｋ-1}、すなわち上記ICの各整数に値が隣接する整数の集合とし、
３．IAをIA１={1、...、IB１-1}、IAｐ ={IBｐ+1、...、N}、2≦k≦p-1のときIA={ICｋ−１+2、ICｋ-2}、すなわち上記集合ICの各整数に値が隣接しない整数の集合とし、
（３）前記ｐ台のプロセッサの少なくとも１台を用いて、前記三重対角行列Ａの行及び列を、並び替え前の行番号と列番号がそれぞれIA１、IA２、...、IAp、IB１、IB２、...、IBp、ICと並ぶ集合の整数の順になるように並び替えるステップと、
（４）並べ替えで得た行列の列集合IA１、IA２、...、IAp(並び替え前からの列番号の集合で表示)のそれぞれに対し、該行列を上三角行列と下三角行列との積の形に分解するＬＵ分解処理を前記ｐ台のプロセッサを用いて並列に実行し、該ＬＵ分解処理の際に、該ＬＵ分解処理の計算対象が第Ｋ列であるとき第Ｋ列の第Ｋ行目以降の要素の中で最大の絶対値を持つ要素を選び当該最大の絶対値を持つ要素の行と第Ｋ行とを入れ替えて以降の計算を行うステップと、
（５）前記列集合IA１、IA２、...、IApのそれぞれに対するＬＵ分解処理の結果に基づき、前記並べ替えで得た行列の列集合IB１、IB２、...、IBp、IC(並び替え前からの列番号の集合で表示)に対し、ＬＵ分解処理を前記ｐ台のプロセッサの少なくとも１台を用いて実行するステップと、
（６）前記ｐ台のプロセッサの少なくとも１台を用いて、ＬＵ分解処理が実行された行列および右辺ベクトルbに対して前進消去演算と後退代入演算とを実行して前記解ベクトルxを求めるステップと、
（７）前記解ベクトルxを上記出力装置より出力するステップとを有することを特徴とする並列計算機システムを用いた連立一次方程式の求解方法。
【請求項２】
請求項１に記載の連立一次方程式の求解方法において、前記集合ICを構成する整数列が等差数列となるように前記分割ステップを実行することを特徴とする連立一次方程式の求解方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は並列計算機システムを用いた連立一次方程式の求解方法とこれを用いた数値シミュレーション実行方法等に関する。
【０００２】
【従来の技術】
連立一次方程式の解を求める計算は、科学技術計算におけるもっとも基本的な計算の一つであり、様々な分野で利用される。たとえば半導体デバイスの特性を求めるための電子状態計算では、個々の電子の軌道を求めるために大規模行列の固有値・固有ベクトル計算が必要であるが、この固有ベクトル計算では、次数が数千以上の三重対角行列（図１の１）を係数とする連立一次方程式を多数回解くことが必要となる。また、建造物などの構造解析を有限要素法または差分法を用いて行う場合には、各節点での応力や変位を求めるために、剛性行列と呼ばれる行列を係数とする連立一次方程式を解くことが必要である。差分法における解析メッシュの例と、それに対応する剛性行列とをそれぞれ図２、図３に示す。剛性行列の次数はメッシュ内の節点4の個数と等しく、その各行および各列にそれぞれ節点が対応し、節点I（図２の47）と節点J（図２の48）とがメッシュ上で隣接する場合に限り、剛性行列の第（I、 J）要素（図３の49）および第（J、 I）要素（図３の50）は零でない値を持つ。節点の数、すなわち連立一次方程式の次数は数万から数十万に上ることが普通である。衝撃が加わった場合の変形など時間変化を求める場合、あるいは強い非線形性を持つ材料を扱う非線形解析の場合などでは、時間変化の１ステップごと、あるいは収束計算の１ステップごとにこの連立一次方程式を解く必要があり、極めて多数の連立一次方程式を解くことが必要となる。また、LSI設計などに用いられる回路シミュレーションでも、回路の動的特性を計算するため、回路網行列と呼ばれる行列を係数行列とする連立一次方程式を、時間変化の１ステップごとに解くことが必要である。回路シミュレーションにおける回路の例と、それに対応する回路網行列とをそれぞれ図４、図５に示す。回路網行列の次数は回路内の複数の素子（図４の素子７、51など）が接続する節点（図４の47、49など）の個数と等しく、その各行および各列にそれぞれ節点が対応し、節点i（図４の節点47）と節点j（図４の節点49）とが回路上で１個の素子を通して結合している場合に限り、回路網行列の第（i、 j）要素（図５の要素53）および第（j、 i）要素（図５の要素54）は零でない値を持つ。回路内の節点の数、すなわち連立一次方程式の次数は数万以上にのぼることが普通である。この場合も、強い非線形性を持つ回路素子を扱う場合には、各時間ステップにおいて収束計算を行う必要が生じ、連立一次方程式を解く回数は更に多くなる。
【０００３】
以上に挙げた計算に共通していることとして、解くべき連立一次方程式の次数が数千から数十万と極めて大きいこと、一つの問題の中で連立一次方程式を多数回繰り返して解く必要があること、解を求める際に誤差が入り込みやすい悪条件の問題が多いこと、係数行列が密行列ではなく疎行列であり、かつ三重対角行列、剛性行列、回路網行列などの特殊な非零構造を持つ行列であることなどが挙げられる。従来、このような連立一次方程式を解く方法としては、逐次型計算機上でのLU分解（LUはLower Upperの略で、行列を上三角行列と下三角行列の積の形に分解すること）による方法があった（たとえばI.S.Duff & J.K.Reid: Some Design Features of a Sparse Matrix Code、 ACM Transactions on Mathematical Software、 Vol. 5、 No. 1、 pp. 18-35 (1979)参照）。LU分解による方法は、悪条件の問題でも精度良く解けること、同じ非零構造を持つ連立一次方程式を多数回解く場合に計算量を削減できる効率的な方法があること、係数行列の特殊な非零構造を利用して計算量を削減できることなどから、これらの問題に対して広く使われてきた。
【０００４】
近年、電子状態計算、構造解析、回路シミュレーションなどの数値シミュレーションの利用が本格化し、シミュレーションにより実験を代替できるような、より精密なモデル化の必要性が高まってきた。そのためには、非線形性を取り入れたより複雑なモデル、より多くの節点を持つ詳細なメッシュ、より多くの回路素子を考慮した回路シミュレーションなどが必要であり、シミュレーション計算の中核である連立一次方程式を解くために必要な計算量も、飛躍的に大きくなりつつある。
【０００５】
このような大規模な連立一次方程式を解くための手段としては、並列計算機が有力である。並列計算機は数十個から数万個に上る多数の高速プロセッサをネットワークで結んだシステムであり、従来の逐次型計算機に比べ、プロセッサ台数を増やすことでピーク性能をいくらでも高めることができるという利点を持つ。しかし、並列計算機の性能を引き出すには、プロセッサ間に処理を均等に分割でき、プロセッサ間での通信がなるべく少なくて済む並列性の高い計算方法が必要である。
【０００６】
そのため、LU分解による連立一次方程式の解法を並列計算機向けに改良しようという研究が、従来から行われてきた。その中でよく使われている方法として、Nested Dissection法（以下、ND法と略す。たとえばK.A.Gallivan et.al.: Parallel Algorithms for Matrix Computations、 SIAM、 Philadelphia、 1990参照）が挙げられる。ND法は、三重対角行列、剛性行列、回路網行列などの特殊な非零構造を持つ行列を、その非零構造を利用して少ない計算量で解く方法であり、計算の並列性を高めて並列計算機のプロセッサを有効に利用できるようになっている。以下、三重対角行列の場合を例に取り、逐次型計算機向けのアルゴリズムと対比させながらND法のアルゴリズムを述べる。
【０００７】
まず、図１の三重対角行列に対応する格子を考える。いま、1から行列の次数Nまでの番号が付いたN個の節点を考え、行列の第（I、 J）要素および第（J、 I）要素が零でない値を持つときに限り、節点Iと節点Jとを結ぶ。すると、三重対角行列に対応する格子は、図６に示す１次元格子になる。
【０００８】
逐次型計算機では、三重対角行列に対し、第１列目から順にLU分解を行っていく。第K列に対するLU分解の演算は、第K列より後のすべての行Iおよび列Jに対して、AIJ=AIJ-AIK*AKJ/AKKという演算を行うことであるが、三重対角行列の場合には、三重対角以外の部分が零であることから、実際に計算するのはAK+1、K+1のみとなる。このアルゴリズムでは、A11、 A22、 ... 、 ANNを
順に計算して行かなければならないため、並列性はない。
【０００９】
これに対してND法では、並列計算機のプロセッサ台数をpとするとき、行列の列集合をp個のほぼ等しい大きさの集合A1、 A2、 ... 、 Apと、それらの境界に当たるp-1本の列の集合Cに分割する。そして、Cに属する(p-1)本の列が行列の最後に来るよう、行列の行と列をそれぞれ入れ替える。p=2の場合に、入れ替えを行った結果の行列を図７に示す。この行列では、入れ替えにより、もともとの三重対角行列にあった非零要素の一部（図７の要素14）がなくなった結果、A1に属する列のLU分解により、A2に属する列が影響を受けなくなっている。したがって、A1の部分のLU分解の終了を待つことなく、A2の部分のLU分解を独立に行うことができる。集合Cに属する列の本数はp-1本であるから、三重対角行列の次数Nがプロセッサ台数pに比べて十分大きければ、ND法は極めて高い並列性を達成できる。ND法のアルゴリズムは、対応する１次元格子（図８）上で見た場合、格子をプロセッサ台数に等しい数の部分領域15と、その間の境界点16とに分割し、部分領域15に属する節点に最初に番号を振り、境界点16に最後に番号を振るように番号の付け替えを行って、最初に部分領域内部でそれぞれ独立にLU分解を行い、その後に境界点の部分に対する処理を行って分解を完成させるアルゴリズムと見ることができる。
【００１０】
ND法は三重対角行列だけでなく、より広い範囲の疎行列に対して適用可能である。たとえば、有限要素法または差分法で生じる剛性行列に対しては、その行列に対応するメッシュを考え、それをプロセッサ台数と等しい数の部分領域と、それら部分領域の間の境界点の集合に分割する。ここで、部分領域への分割は、２つの異なる部分領域に属する点が、メッシュ上で隣接することがないように取るものとする。このとき、各部分領域に属する節点の番号の集合をA1、 A2、 ... 、 Apとし、境界点に属する節点の集合をCとして、A1、 A2、 ... 、 Apの節点に最初に番号を振り、Cの節点に最後に番号を振るようにして、これに対応した行列の行と列の並べ替えを行った後でLU分解を行えば、各部分領域に属する列に対するLU分解は、お互いに影響を及ぼさないので独立に行うことができ、これらすべての分解が終わった後にCの部分の分解のみを逐次的に行えばよいことになる。この場合も、プロセッサ台数pに対して節点の数Nが十分に大きければ、境界点の数は相対的に少なくなり、LU分解のほとんどの計算を並列に行うことができる。
【００１１】
また、回路シミュレーションで生じる回路網行列に対しては、その行列に対応する回路を考え、それをプロセッサ台数と等しい数の部分回路と、それら部分回路の間の境界部分の回路に分割する。ここで、部分回路への分割は、２つの異なる部分回路に属する節点が、回路上で１つの素子で結ばれることのないように行うものとする。このとき、各部分回路に属する節点の番号の集合をA1、 A2、 ... 、 Apとし、境界部分の回路に属する節点の集合をCとして、A1、 A2、 ... 、 Apの節点に最初に番号を振り、Cの節点に最後に番号を振るようにして、これに対応した行列の行と列の並べ替えを行った後でLU分解を行えば、各部分回路に属する列に対するLU分解は、お互いに影響を及ぼさないので独立に行うことができ、これらすべての分解が終わった後にCの部分の分解のみを逐次的に行えばよいことになる。この場合も、プロセッサ台数pに対して回路内の節点の数Nが十分に大きければ、境界部分の回路内の節点の数は相対的に少なくなり、LU分解のほとんどの計算を並列に行うことができ、並列計算機のプロセッサを有効に利用することができる。
【００１２】
【発明が解決しようとする課題】
上記のND法による連立一時方程式の解法では、三重対角行列の行と列とを入れ替えて図７のような形に変形することにより、並列性を確保していた。しかし、LU分解を用いて連立一次方程式を解く場合、解に誤差が入りやすい悪条件の問題でも正しく解くためには、一般にピボット選択という操作を取り入れなければならない。ピボット選択とは、第K列に対するLU分解の演算AIJ=AIJ-AIK*AKJ/AKKにおいて、右辺の分母に来る値AKKがなるべく大きくなるように、第K列の第K行目以降の要素の中で最大の絶対値を持つ要素A K′Kを選び、第K行と第K′行とを入れ替えて以降の計算を行う方法である。
【００１３】
ところが、ND法においてピボット選択を行うと、行の入れ替えを行うために、図７の形が崩れてしまう。その結果、部分領域間の並列性は失われてしまい、ND法のメリットがなくなる。プロセッサ台数ｐが２の場合を例に取り、この様子を図９に示す。ピボット選択を行った結果、非零要素の数は増大し、かつ、列集合IA1の部分のLU分解の結果として、列集合IA2の最初の列の要素（図９の要素17）が影響を受けてしまう。その結果、IA1のLU分解が終了しないとIA2のLU分解が開始できなくなり、独立性が失われる。
【００１４】
ここでは三重対角行列を例に取って示したが、有限要素法や差分法で生じる剛性行列の場合、回路シミュレーションで生じる回路網行列の場合にも、ND法にピボット選択を取り入れると、行列の形が崩れてしまい、部分領域間の並列性は失われる。
【００１５】
このような問題から、従来のND法は、悪条件でないことが予めわかっており、ピボット選択なしでも正しい解が得られることが理論的に保証されている一部の問題にしか使えず、その適用範囲が大きく限定されていた。特に、分子計算のための固有値・固有ベクトル計算で出てくる連立一次方程式、回路シミュレーションで生じる連立一次方程式は、典型的な悪条件の問題であり、従来これらの問題をND法を用いて解くことはできなかった。
【００１６】
本発明では、以上に述べたようなND法の問題点を解決し、悪条件の問題にも適用できて、かつ並列計算機の多数のプロセッサを有効に利用できるような、新しい連立一次方程式の解法を提供することを目的とする。
【００１７】
【課題を解決するための手段】
上記目的を達成するため、本発明では、ND法における行と列の入れ替え方式を改良し、ピボット選択を行っても、並列性が保たれるようにした。
【００１８】
ND法でピボット選択を取り入れたために部分領域間の演算の独立性が崩れる過程をよく見ると、部分領域に含まれる点のうち、境界点に接する点に対応する列のLU分解で問題が生じていることがわかる。そこで本発明では、ND法における各部分領域の点を、境界点に接する内点と、それ以外の純粋な内点とに分割し、各部分領域の純粋な内点のみが最初に来て、次に各部分領域内の境界点に接する内点、最後に境界点が来るように行列の行と列を並べ替える。これにより、ピボット交換で独立性を失わせる可能性のある列は純粋な内点の後に来るため、最初の各部分領域の純粋な内点の部分は、部分領域ごとに独立にLU分解することが可能となる。それ以降の点では並列性は失われるが、プロセッサ台数pに比べて行列の次数Nが十分に大きければ、境界点の数も、境界点に接する内点の数も、純粋な内点の数に比べれば十分小さくなり、本発明の方法ではLU分解の演算の大部分において高い並列性が得られる。
【００１９】
【発明の実施の形態】
（実施例１）
以下、本発明の原理および実施例を、図面により詳細に説明する。ここで実施例として挙げるのは、分散メモリ型の並列計算機上で、三重対角行列を係数とする連立一次方程式を解く場合である。
【００２０】
本方法により得られる並列プログラムを実行する並列計算機システムを図10に示す。本システムは係数行列Aおよび右辺ベクトルbを入力するための入力装置18、それぞれがメモリ26を備えたp台のプロセッサ19を持つ処理装置20、計算結果を出力するための出力装置21、プログラム22およびデータ23を格納するための外部記憶装置24から構成される。
【００２１】
本発明の方法により連立一次方程式の解を求める手順を図11に示す。ただし、並列計算機のプロセッサの数をp個とする。まず、係数行列Aおよび右辺ベクトルbを入力する（処理27）。次に、三重対角行列のプロセッサへの分割を行うために、それぞれが３以上N-2以下で、かつ互いに４以上離れたp-1個の増加する整数の列IC1、 IC2、 ... 、 ICp-1を選び、 IC={IC1、 IC2、 ... ICp-1}とする。（処理28）。ここで、 ICkは、プロセッサk-1の担当する部分領域とプロセッサkの担当する部分領域との境界となるので、IC1、 IC2、 ... 、 ICp-1がなるべく等差数列に近くなるように選ぶと、各プロセッサの負荷が均等になり、連立一次方程式を解くための処理時間が短くて済む。次に、IB1={IC1-1}、IBp={ICp-1+1}、２≦k≦p-1のときIBk={ICk-1+1、 ICk-1}と設定する（処理29）。 IBkは、プロセッサkの担当する部分領域のうちで、境界点の集合ICに接する点の集合である。次に、 IA1={1、 ... 、 IB1-1}、 IAp={ IBp+1、 ... 、 N}、 ２≦k≦p-1のときIAk={ICk-1+2、 ICk-2}と設定する（処理30）。 IAkは、プロセッサkの担当する部分領域のうちで、境界点の集合に接しない点の集合である。以上のように集合IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICをを選んだ後、行列Aの行と列とを、 IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICの順に並べ替える（処理31）。そして、 列集合IA1、 IA2、 ... 、 IApに対するピボット選択付きLU分解を、それぞれ第1、 2、 ...  p番目のプロセッサで独立に行う（処理32）。なお、これら各列集合のLU分解の結果、列集合IB1、 IB2、 ... 、 IBp、 ICに属する列の要素が影響を受けることがあるが、それもこの段階で担当のプロセッサが処理を行っておく。列集合IA1、 IA2、 ... 、 IApに対するLU分解が終了したら、次に、その結果に基づいて、列集合IB1、 IB2、 ... 、 IBp、 ICに対するピボット選択付きLU分解を順に行う（処理33）。このLU分解は、単一のプロセッサが行ってもよいし、また、上記のp個のプロセッサのうちの複数個が協調して行ってもよい。以上により、行列AのLU分解が終了したら、その結果を用いて、前進消去演算（処理34）および後退代入演算（処理35）を行って解ベクトルを求め、出力装置より解を出力する（処理36）。
【００２２】
以上の処理において、プロセッサの台数が２の場合を例に取り、本発明の方法ではピボット選択を行っても各IA1、 IA2、 ... 、 IApに対するLU分解が独立に行えることを示す。まず、行列Aの行と列とを、 IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICの順に並べ替えた結果を図12に示す。また、列集合IA1に対するLU分解が終わった段階での、行列の非零構造および影響を受けた要素を図13に示す。図12よりわかるように、境界点に接する部分領域の点を後ろに持っていった結果、ピボット選択付きのND法において並列性を損なう原因となっていた要素38が、IA2に属する列の後ろに行く。この結果、図13から明らかなように、列集合IA1に対するLU分解によって、列集合IA2に属する列の要素は影響を受けなくなる。したがって、列集合IA2に属する列のLU分解は、 IA2に属する列のLU分解と独立に並行して行うことができる。以上ではpが２のときを示したが、２以上のpについても、各IAkはそれ以前のIAkのLU分解の影響を受けず、p個の列集合IA1、 IA2、 ... 、 IApに対するLU分解はそれぞれ独立に行ってよいことがわかる。
【００２３】
以上では、本発明の解法を分散メモリ型計算機上で実現する方法を示したが、本発明の方法では列集合IA1、 IA2、 ... 、 IApに対するLU分解はそれぞれ独立に行えるため、本発明の方法を共有メモリ型並列計算機にも同様に適用できることは明らかである。
【００２４】
（実施例２）
別の非零構造を持つ疎行列の例として、有限要素法もしくは差分法で生じる剛性行列に対して本発明を適用する方法を説明する。本方法により得られる並列プログラムを実行する並列計算機システムは、実施例１のものと同じである。
【００２５】
本発明の方法により剛性行列を係数とする連立一次方程式の解を求める手順を図14に示す。ただし、pは並列計算機のプロセッサの数を示す。まず、係数行列Aおよび右辺ベクトルbを入力する（処理27）。次に、剛性行列のプロセッサへの分割を行うために、行列Aに対応するメッシュを、p個の部分領域およびそれらの境界の領域ICに分割する（処理42）。ここで、各部分領域に含まれる節点の数がなるべく均等になるように部分領域を取ると、LU分解において各プロセッサの負荷が均等になり、連立一次方程式を解くための処理時間が短くて済む。また、これと同時に、境界の領域ICに属する節点の数がなるべく少なくなるように部分領域を取ると、並列に処理できない境界点の部分の割合が減少するため、並列計算機の効率的な利用が可能となり、連立一次方程式を解くための処理時間がさらに短くて済む。次に、上記p個の部分領域のそれぞれにおいて、境界点に接する節点のみからなる集合をそれぞれIB1、 IB2、 ... IBpに設定し（処理43）、それ以外の節点の集合をそれぞれIA1、 IA2、 ... 、 IApに設定する（処理44）。以上のように集合IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICをを選んだ後、行列Aの行と列とを、 IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICの順に並べ替える（処理31）。列集合IA1、 IA2、 ... 、 IApに対するピボット選択付きLU分解を、それぞれ第1、 2、 ...  p番目のプロセッサで独立に行う（処理32）。なお、これら各列集合のLU分解の結果、列集合IB1、 IB2、 ... 、 IBp、 ICに属する列の要素が影響を受けることがあるが、それもこの段階で担当のプロセッサが処理を行っておく。処理32で、列集合IA1、 IA2、 ... 、 IApに対するLU分解が終了したら、次に、その結果に基づいて、列集合IB1、 IB2、 ... 、 IBp、 ICに対するピボット選択付きLU分解を順に行う（処理33）。このLU分解は、単一のプロセッサが行ってもよいし、p個のプロセッサのうちの複数個が協調して行ってもよい。以上により、行列AのLU分解が終了したら、その結果を用いて、前進消去演算（処理34）および後退代入演算（処理35）を行って解ベクトルを求め、出力装置21より解を出力する（処理36）。
【００２６】
以上の処理において、プロセッサの台数が２の場合を例に取り、本発明の方法ではピボット選択を行っても各IA1、 IA2、 ... 、 IApに対するLU分解が独立に行えることを示す。まず、行列Aの行と列とを、 IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICの順に並べ替えた結果を図15に示す。また、列集合IA1に対するLU分解が終わった段階での、行列の非零構造および影響を受けた要素を図16に示す。なお、図15では、これまでの図と異なり、非零要素の一個一個を四角で示さずに、非零要素の存在する可能性のある領域をまとめて網掛けで、それ以外の部分を白抜きで示してある。また、図16では、IA1に対するLU分解の影響を受けた非零要素を濃い網掛け41で、受けていない非零要素を薄い網掛け2で、零要素を白抜き3で示してある。図16から明らかなように、列集合IA1に対するLU分解によって、列集合IA2に属する列の要素は影響を受けない。したがって、列集合IA2に属する列のLU分解は、 IA2に属する列のLU分解と独立に並行して行うことができる。以上ではp=2の場合を説明したが、p>2の場合についても、各IAkはそれ以前のIAkのLU分解の影響を受けず、p個の列集合IA1、 IA2、 ... 、 IApに対するLU分解はそれぞれ独立に行ってよい。
【００２７】
本実施例においても、分散メモリ型計算機上での実現方法を示したが、本発明の方法では列集合IA1、 IA2、 ... 、 IApに対するLU分解はそれぞれ独立に行えるため、本発明の方法を共有メモリ型並列計算機にも同様に適用できることは明らかである。
【００２８】
（実施例３）
更に別の非零構造を持つ疎行列の例として、回路シミュレーションで生じる回路網行列に対して本発明を適用する方法を説明する。本方法により得られる並列プログラムを実行する並列計算機システムは、実施例１のものと同じである。
【００２９】
一般に回路シミュレーションでは、回路（図４）中において複数の素子（素子7、51など）が接続する点を節点（節点47、49など）とし、各節点iにおける電圧Viを未知数として方程式を立てる。この手順を図17に示す。まず、回路中の各素子のインピーダンスΩ、素子間の接続関係を表す情報、および外部からの入力電流を入力する（処理59）。節点iにおける方程式を立てるには、まずiに接続する各素子に対して、素子のもう一方の端点での電圧Vjおよび素子のインピーダンスΩより、その素子を通して節点iに流れ込む電流(Vj-Vi)/Ωを求める（処理63）。次に、これをiに接続するすべての素子に対して足し合わせ（処理64）、それを０と置く。これは節点iにおける電流保存則であり、このようにして各節点においてそれぞれ１つの方程式ができる。これらをすべての節点について連立させることにより、Viに関する連立一次方程式ができる。これを行列の形で書くことにより、回路網行列を作る（処理65）。行列の第（i、 j）要素が非零ならば、これは節点iと節点jとが素子を通してつながっているということだから、第（j、 i）要素も非零となる。したがって、回路網行列は対称な非零構造を持つ。また、右辺ベクトルの各要素を、外部からの電流入力がある節点に対してはその電流の値、ない節点に対しては０と設定し（処理66）、回路網行列および右辺ベクトルを出力する（処理67）。
【００３０】
本発明の方法により回路網行列を係数とする連立一次方程式の解を求める手順を図18に示す。ただし、pは並列計算機のプロセッサの数を示す。まず、回路網行列Aおよび右辺ベクトルbを入力する（処理55）。次に、回路網行列のプロセッサへの分割を行うために、行列Aに対応する回路を、p個の部分回路およびそれらの境界部の回路ICに分割する（処理56）。なお、境界部の回路は、２つの異なる部分回路に含まれる節点が１つの素子で結ばれることのないように取るものとする。ここで、各部分回路に含まれる節点の数がなるべく均等になるように部分回路を取ると、LU分解において各プロセッサの負荷が均等になり、連立一次方程式を解くための処理時間が短くて済む。また、これと同時に、境界部の回路ICに属する節点の数がなるべく少なくなるように部分回路を取ると、並列に処理できない境界点の部分の割合が減少するため、並列計算機の効率的な利用が可能となり、連立一次方程式を解くための処理時間がさらに短くて済む。次に、上記p個の部分回路のそれぞれにおいて、境界部の回路に属する節点に１つの素子を通して接続する節点のみからなる集合をそれぞれIB1、 IB2、 ... IBpに設定し（処理57）、それ以外の節点の集合をそれぞれIA1、 IA2、 ... 、 IApに設定する（処理58）。以上のように集合IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICをを選んだ後、行列Aの行と列とを、 IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICの順に並べ替える（処理31）。列集合IA1、 IA2、 ... 、 IApに対するピボットLU分解を、それぞれ第1、 2、 ...  p番目のプロセッサで独立に行う（処理32）。なお、これら各列集合のLU分解の結果、列集合IB1、 IB2、 ... 、 IBp、 ICに属する列の要素が影響を受けることがあるが、それもこの段階で担当のプロセッサが処理を行っておく。処理32で、列集合IA1、 IA2、 ... 、 IApに対するLU分解が終了したら、その結果に基づいて、列集合IB1、 IB2、 ... 、 IBp、 ICに対するピボット選択付きLU分解を順に行う（処理33）。このLU分解は、単一のプロセッサが行ってもよいし、p個のプロセッサのうちの複数個が協調して行ってもよい。以上により、行列AのLU分解が終了したら、その結果を用いて、前進消去演算（処理34）および後退代入演算（処理35）を行って解ベクトルを求め、出力装置21より解を出力する（処理36）。この解ベクトルの各要素が、各節点での電圧を表す。
【００３１】
以上の処理において、プロセッサの台数が２の場合を例に取り、本発明の方法ではピボット選択を行っても各IA1、 IA2、 ... 、 IApに対するLU分解が独立に行えることは、実施例２の場合と全く同様にして示せる。まず、行列Aの行と列とを、 IA1、 IA2、 ... 、 IAp、 IB1、 IB2、 ... 、 IBp、 ICの順に並べ替えた結果を図15に示す。また、列集合IA1に対するLU分解が終わった段階での、行列の非零構造および影響を受けた要素を図16に示す。なお、図15では、これまでの図と異なり、非零要素の一個一個を四角で示さずに、非零要素の存在する可能性のある領域をまとめて網掛けで、それ以外の部分を白抜きで示してある。また、図16では、IA1に対するLU分解の影響を受けた非零要素を濃い網掛け41で、受けていない非零要素を薄い網掛け2で、零要素を白抜き3で示してある。図16から明らかなように、列集合IA1に対するLU分解によって、列集合IA2に属する列の要素は影響を受けない。したがって、列集合IA2に属する列のLU分解は、 IA2に属する列のLU分解と独立に並行して行うことができる。以上、p=2の場合を説明したが、p>2の場合についても、各IAkはそれ以前のIAkのLU分解の影響を受けず、p個の列集合IA1、 IA2、 ... 、 IApに対するLU分解はそれぞれ独立に行ってよい。
【００３２】
本実施例においても、分散メモリ型計算機上での実現方法を示したが、本発明の方法では列集合IA1、 IA2、 ... 、 IApに対するLU分解はそれぞれ独立に行えるため、本発明の方法を共有メモリ型並列計算機にも同様に適用できることは明らかである。
【００３３】
（実施例４）
本発明の他の利用例として、並列計算機を用いて電子状態計算を行うシステムに対して本発明を適用する方法を説明する。本発明により得られる並列プログラムを実行する並列計算機システムは、実施例１のものと同じである。
【００３４】
本発明の方法により電子状態計算を行う手順を図19に示す。まず、電子状態計算を行うべき系において、系を構成する原子の種類および座標、電子の個数Nを入力する（処理68）。次に、通常の量子力学の手順に基づき、各電子の波動関数ψjを複素指数関数exp(ikr)によりψj=Σcjk exp(ikr)と展開し、系の全エネルギーHを、各電子の運動エネルギーT({cjk})、電子間のクーロン相互作用エネルギーEc({cjk})、電子間の交換相関相互作用エネルギーExc({cjk})、電子と原子間の相互作用エネルギーU({cjk})の和として、H({cjk})=T({cjk})+Ec({cjk})+Exc({cjk})+ U({cjk})と表す（処理69）。量子力学の定理より、求める波動関数は、関数H({cjk})を最小化し、かつ異なる電子の波動関数が直交するような{cjk}により与えられ、この２つの条件を満たす{cjk}は、H({cjk})を{cjk}の２次形式と見たときの２次形式を表す対称行列Hの、下からN個の固有値に属する固有ベクトルとして与えられる。そこでHより、２次形式を表す対称行列Hを求める（処理70）。次に、固有値と固有ベクトルを求めるため、Hをハウスホルダー変換と呼ばれる通常の手続きにより、三重対角行列Tに相似変換する（処理71）。次に三重対角行列の下からN個の固有値λ1、...、λNを求め（処理72）、j=1、...、Nについて、三重対角行列Tを係数とする連立一次方程式Tvj=λjvjを解いて固有ベクトルvjを求める（処理73）。この求解処理では、実施例１の方法を適用する（処理74）。こうして固有ベクトルが求められると、{v1、...、vN}が求める{cjk}であり、これをψj=Σcjk exp(ikr)に代入して波動関数ψjを求め（処理75）、H({cjk})に代入して系の全エネルギーHを求める（処理76）。こうして求めたψjとHを出力装置から出力する（処理77）。
【００３５】
本方法において、入力から三重対角行列Tおよび固有値λ1、...、λNを求める部分、および{v1、...、vN}からψjとHを求める部分は、単一のプロセッサが行ってもよいし、複数のプロセッサが行ってもよい。
【００３６】
以上では、本発明の解法を図10の分散メモリ型計算機上で実現する方法を示したが、実施例１の場合と同様、本発明の方法を共有メモリ型並列計算機にも適用できることは明らかである。
【００３７】
【発明の効果】
以上説明したように、本発明によれば、並列計算機を用いて三重対角行列、剛性行列、回路網行列などを係数とする連立一次方程式を解く場合において、並列性を保ったままピボット選択を行うことが可能になり、従来並列計算機の適用が難しかった悪条件の問題に対しても、並列計算機を有効に利用した解法が実現できる。また、これにより、連立一次方程式の解を求める計算が大きな部分を占める、電子状態計算、有限要素法または差分法の計算、回路シミュレーションの計算などにおいて、計算に要する時間を大幅に削減できる。
【図面の簡単な説明】
【図１】三重対角行列を示す図。
【図２】差分法における解析メッシュを示す図。
【図３】差分法における剛性行列を示す図。
【図４】回路の例を示す図。
【図５】回路網行列を示す図。
【図６】三重対角行列に対応する１次元格子を示す図。
【図７】  ND法により行と列の入れ替えを行った三重対角行列を示す図。
【図８】  ND法により行と列の入れ替えを行った三重対角行列に対応する１次元格子を示す図。
【図９】  ND法においてピボット選択を行った結果を示す図。
【図１０】本発明の実施例１および実施例２を適用すべき並列計算機の図。
【図１１】本発明の方法により三重対角行列を係数とする連立一次方程式の解を求める手順を示す図。
【図１２】本発明の方法により行と列の入れ替えを行った三重対角行列を示す図。
【図１３】本発明の方法においてピボット選択を行った結果を示す図。
【図１４】本発明の方法により剛性行列を係数とする連立一次方程式の解を求める手順を示す図。
【図１５】本発明の方法により行と列の入れ替えを行った剛性行列を示す図。
【図１６】本発明の方法においてピボット選択を行った結果を示す図。
【図１７】回路網行列を作成する手順を示す図。
【図１８】本発明の方法により回路網行列を係数とする連立一次方程式の解を求める手順を示す図。
【図１９】本発明の方法により電子状態計算を行う手順を示す図。
【符号の説明】
１：三重対角行列、２：非零要素、３：零要素、４：節点、５：剛性行列、６：回路、７：コンデンサー、８：抵抗、９：トランジスタ、10：回路網行列、11：ND法により行と列の入れ替えを行った三重対角行列、12：A1のLU分解により影響を受ける要素、13：A1のLU分解により影響を受けない要素、14：もともとの三重対角行列にあった要素、15：各プロセッサが担当する部分領域、16：境界点、17：列集合IA1のLU分解の影響を受ける列集合IA2の要素、18：入力装置、19：プロセッサ、20：処理装置、21：出力装置、22：プログラム、23：データ、24：外部記憶装置、25：メモリ、26：スタート、27：係数行列Aおよび右辺ベクトルbの入力、28：集合ICの設定、29：集合IB1、 IB2、 ... 、 IBpの設定、30：集合IA1、 IA2、 ... 、 IApの設定、31：行列の並べ替え、32：各プロセッサで列集合IA1、 IA2、 ... 、 IAp、 ICに対するLU分解を独立に行う、33：列集合IB1、 IB2、 ... 、 IBpに対するLU分解、34：前進消去演算、35：後退代入演算、36：解ベクトルxの出力、37：終了、38：ND法でピボット選択を行ったときに並列性を損なう原因となっていた要素、39：本発明の方法により行と列の入れ替えを行った三重対角行列、40：本発明の方法においてピボット選択を行った結果、41：列集合IA1のLU分解の影響を受ける要素、42：列集合IA1のLU分解の影響を受ける要素、43：境界点に接する内点の集合をIB1、 IB2、 ... 、 IBpに設定、44：それ以外の内点の集合をIA1、 IA2、 ... 、 IApに設定、45：本発明の方法により行と列の入れ替えを行った剛性行列を示す図、46：本発明の方法においてピボット選択を行った結果を示す図、47：節点、48：節点番号、49：節点、50：節点、51：素子、52：素子、53：回路網行列の第(I、J)要素、54：回路網行列の第(J、I)要素、55：回路網行列および右辺ベクトルbの入力、56：回路の部分回路と境界部の回路への分割、57：IB1、 IB2、 ... 、 IBpの設定、58： IA1、 IA2、 ...、 IApの設定、59： 各素子のインピーダンスおよび素子間の接続関係を表す情報の入力、60： 各素子のインピーダンス、素子間の接続関係を表す情報、および入力電流の入力、61： すべての節点iについて繰り返す、62： iに接続するすべての素子について繰り返す、63： 素子を通して節点jからiに流れ込む電流(Vj-Vi)/Ωを計算、64： (Vj-Vi)/Ωをjについて足し合わせ、和を０とおいて方程式を立てる、65：すべての節点iに関する方程式を連立させ、回路網行列を作る、66：入力電流から右辺ベクトルを作る、67：回路網行列および右辺ベクトルの出力、68：原子の種類と座標および電子の個数Nの入力、69：系の全エネルギーHを{cjk}の関数として表す、70：２次形式を表す行列Hを求める、71：H'を三重対角行列Tに相似変換する、72：Tの固有値λ1、...、λNを求める、73：j=1、 ...、 Nについて繰り返す、74：実施例１の方法によりTvj=λjvjを解く、75：波動関数ψjを求める、76：全エネルギーHを求める、77：ψjとHの出力。
【図１】
0004292604.tif 000002
【図２】
0004292604.tif 000003
【図３】
0004292604.tif 000004
【図４】
0004292604.tif 000005
【図５】
0004292604.tif 000006
【図６】
0004292604.tif 000007
【図７】
0004292604.tif 000008
【図８】
0004292604.tif 000009
【図９】
0004292604.tif 000010
【図１０】
0004292604.tif 000011
【図１１】
0004292604.tif 000012
【図１２】
0004292604.tif 000013
【図１３】
0004292604.tif 000014
【図１４】
0004292604.tif 000015
【図１５】
0004292604.tif 000016
【図１６】
0004292604.tif 000017
【図１７】
0004292604.tif 000018
【図１８】
0004292604.tif 000019
【図１９】
0004292604.tif 000020