(57)【特許請求の範囲】
【請求項１】
親プロセッサと、各々にメモリを有する複数の子プロセッサと、これらプロセッサに共有される入力用外部メモリ及び出力用外部メモリとを有する分散メモリ型並列計算機システムにおいて大規模行列の固有ベクトルの数値計算を行う並列機向け固有ベクトル求解方法であって、
（１）行列データと固有値データを上記入力用外部メモリから各子プロセッサに付随のメモリに入力する入力処理、
（２）個々の固有値について他の全ての固有値との距離がある基準値以下かどうかを上記親プロセッサによって判定する固有値距離解析処理、
（３）上記固有値距離解析処理において求められた結果を用いて固有ベクトルを解く順序を上記親プロセッサによって決定する固有ベクトル計算優先度決定処理、
（４）上記固有値距離解析処理の結果と個々のプロセッサの性能に基づき、上記親プロセッサによって各固有ベクトル計算を各子プロセッサに分配する固有ベクトル計算分配処理、
（５）上記固有ベクトル計算優先度決定処理並びに固有ベクトル計算分配処理の結果に基づいて、必要があれば固有ベクトルのプロセッサ間通信をしながら、上記固有値距離解析処理の結果の範囲にある求解済みの固有ベクトルと直交化をしながら、逆反復法によって固有ベクトルを上記各子プロセッサによって求める固有ベクトル計算処理、
（６）上記求められた固有ベクトルを上記親プロセッサによって統合し、出力用外部メモリへ出力する出力処理、
との処理手順に従うことを特徴とする並列機向け固有ベクトル求解方法。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
本発明は分散メモリ型並列計算機（以下、単に並列計算機）における大規模行列の固有ベクトル計算方法に係り、特に、弾性体解析、流体解析などの数値シミュレーションにあらわれる大型行列の固有ベクトル計算の高速実行方法に関する。
【０００２】
【従来の技術】
固有値eがわかっている場合に、適当な初期ベクトルから逆反復法（数１）によって固有ベクトルvを求めることができる。
【０００３】
【数１】
0003648286.tif 000002
【０００４】
ここで、この逆反復法をより高速にかつ高精度にするため、１回の反復の度に、すでに求められた他の固有ベクトルｗと数２のように直交化させることにより固有ベクトルを計算することが考えられていた。
【０００５】
【数２】
0003648286.tif 000003
【０００６】
しかし、この直交化の過程は、非常に計算量が多いため、特に近くの固有値に属する固有ベクトルとだけ直交化させることによって、数値計算上現れる大部分の誤差を除去することが一般的に行われている。この方法による固有ベクトル計算を、以下、直交化付き逆反復法と呼ぶことにする。
【０００７】
この直交化付き逆反復法による固有ベクトル計算を、高速処理のために、並列計算機上で実行することが、従来行われきた（Wilkinson, Reinsch, Linear Algebra, Sprnger-Verlag pp418-pp439）が、この処理方法について図１３と図１４を用いて説明する。
【０００８】
まず、図１３において固有値のグループ化が行われる。数直線131上に今、小さいものから順に、固有値e1、e2、などがあるとする。図１３では説明のためその位置を矢印で示している。ここで、ある基準epsを用い、隣り合う固有値がeps以下ならその二つの固有値は同じグループとするという規則に従い、全ての固有値をグループに分け、グループ132、グループ133、などとする。なお、この処理を以下、固有値のグループ化と呼ぶ。
【０００９】
次に図１４を用いてグループ化された固有値に属する固有ベクトル計算のプロセッサへの分配方法を説明する。ここでは、グループ数を使用するプロセッサ数で割った数のグループ数を一個のプロセッサが担当し、あまったグループは一個ずつどれかのプロセッサが担当することにする。図１４では５個に別れたグループを４個のプロセッサが担当する場合を示している。この場合４個のプロセッサのうち、３個が１グループを担当し、一個のプロセッサが２グループを担当する。そこで、プロセッサ１がグループ141と142を担当し、プロセッサ２がグループ143を担当し、プロセッサ３がグループ144を担当し、プロセッサ４がグループ145を担当する。
【００１０】
各プロセッサは、担当するグループ内で直交化付き逆反復法を実行する。ここでは、一個のグループ内の最も小さい固有値に属する固有ベクトルを直交化なしで計算を行った後には、順次、一回の逆反復計算の度に、グループ内の自分の値より小さい固有値に属する固有ベクトルとはすべて直交化を行う、という方法によって直交化付き逆反復法を実現する。
【００１１】
【発明が解決しようとする課題】
並列計算機において、固有値が与えられた場合に、大型行列の固有ベクトルを、直交化付き逆反復法によって高速かつ高精度に求める時、従来の技術は以下に示すような問題がある。
【００１２】
第一の問題は、固有値のグループ化によって、実際の距離は基準値以上でも、同じグループに属したがために必要以上に直交化を行っていることである。即ち、大量の不必要な計算を行っていることである。
【００１３】
第二の問題は、固有値のグループ化によって利用すべきプロセッサの台数に制限があることである。従って、固有値の分布によってはグループの個数が少ないために並列計算実行がほとんど有効に機能しない場合もある。このことは、特にプロセッサ台数が多い超並列計算機のハードウェアの有効利用を考えた場合、非常に大きな問題である。
【００１４】
第三の問題は、固有値のグループ化によって、それぞれのグループを担当するプロセッサエレメントの計算量が決まってしまうため、全体の計算量と使用する個々のプロセッサの性能を鑑みたうえで各プロセッサエレメントへの計算量を割り付けることが不可能になっていることである。即ち、プロセッサエレメントのロードバランスの調整が出来ず、特にプロセッサ台数が多い超並列計算機のハードウェア並びに、個々のプロセッサの性能が異なるような並列計算機のハードウェアの有効利用を考えた場合、非常に大きな問題である。
【００１５】
【課題を解決するための手段】
上記の課題を解決するために、本発明では、固有値のグループ化を行わず、個々の固有値について、自分の固有値からの距離がある基準以下の範囲にある他の固有値に属する固有ベクトルと直交させる、という方針で直交化付き逆反復法を実現する。
【００１６】
この方針に基づく本発明の並列機向け固有ベクトル解法は次の六つの処理から構成される。
【００１７】
（１）行列データと固有値データの入力処理。
【００１８】
（２）個々の固有値についてその他の固有値との距離がある基準値以下のものがどの固有値かを判定し、直交化させる範囲を決定する、固有値距離解析処理。
【００１９】
（３）固有値距離解析処理において求められた直交化させる範囲を鑑み、優先度を全ての固有値に与えていく。ここでは、自分より小さくかつ直交化の範囲にある固有値の優先度と重ならないよう、できるだけ高い優先度を与える。その後に、優先度の高い固有値に属する固有ベクトルの計算から順に、優先度が同じ場合は値の小さい固有値に属する固有ベクトルから順に計算を行う、という順序を決定する。このように固有ベクトルの計算の順序を決定する、固有ベクトル計算優先度決定処理。
【００２０】
（４）各固有ベクトル計算において、固有値距離解析処理の結果から導かれる直交化の計算量を鑑み、個々のプロセッサにその性能に応じて各固有ベクトル計算を分配する、固有ベクトル分配処理。
【００２１】
（５）固有ベクトル計算優先度決定処理並びに固有ベクトル計算分配処理の結果に基づいて、必要があれば固有ベクトルのプロセッサ間通信をしながら、固有値距離解析処理の結果の範囲の、既に求められた固有ベクトルと直交化をしながら、逆反復法によって固有ベクトルを求める、固有ベクトル計算処理。
【００２２】
（６）各プロセッサにおいて求められた固有ベクトルを統合し出力する、出力処理。
【００２３】
【作用】
上記（２）の固有値距離解析処理によって、固有値同士の距離がある基準以下の固有値に属する固有ベクトルとのみ直交させるため、従来の方法で行われていた必要以上の直交化の計算がなくなり、第一の問題が解決される。
【００２４】
また、従来の方法でのグループ化を回避し、上記（２）の固有値距離解析処理の結果に基づいて個々の固有ベクトルの計算量を決定した後に、上記（４）の固有ベクトル計算分配処理において各プロセッサに固有ベクトル計算を分配しているため、解くべき固有ベクトルの個数までの並列化が可能であり、従って、従来の方法におけるような、固有値の分布に左右されて殆ど並列化が望めないような場合が回避され、第二の問題が解決される。
【００２５】
さらに、従来の方法でのグループ化を回避し、上記（４）の固有ベクトル計算分配処理において、各固有ベクトルに要する計算量と各プロセッサの性能に応じて、一つ一つの固有ベクトルの計算を個々のプロセッサに割り付けることができるため、使用するプロセッサ全体のロードバランスの最適な調整が達成され、第三の問題が解決される。
【００２６】
【実施例】
以下、本発明の第一の実施例として並列機向け固有ベクトル解法の実施例を図１から図１２に基づき説明する。
【００２７】
まず、図１を用いて実施例で利用される全体処理装置の一例について説明する。この全体処理装置はプロセッサ101、入力用外部メモリ102、出力用外部メモリ103、プロセッサ104、プロセッサ105、プロセッサ106、プロセッサ107、メモリ108、メモリ109、メモリ110、メモリ111から構成され、メモリ108はプロセッサ104に、メモリ109はプロセッサ105に、メモリ110はプロセッサ106に、メモリ111はプロセッサ107に付随しているとする。このように、以下の説明では、固有ベクトル計算に使用するプロセッサの数を４として説明する。また以下の処理内容は、全体処理装置においてライブラリ登録されるものとする。
【００２８】
図２は本発明実施例の処理全体の手順を示すフローチャートであり、入力処理21、固有値距離解析処理22、固有ベクトル計算優先度決定処理23、固有ベクトル計算分配処理24、固有ベクトル計算処理25、出力処理26から構成される。
【００２９】
次に、図３を用いて入力処理の具体的な処理内容を説明する。ここでは、入力用外部メモリにある、行列要素、行列の固有値、プロセッサの性能比、固有値距離解析に用いる基準値epsのデータを、各プロセッサｋを通して、そのプロセッサｋに付随するメモリｋに入力する（ブロック31）、という操作をｋが１から４まで繰り返す（ブロック32）。
【００３０】
次に、図４を用いて固有値距離解析処理について説明する。ここでは、各々の固有値について距離が基準値eps未満の他の固有値がどこからどこまでかを解析する。
【００３１】
今、固有値e1、e2、e3、e4、e5、e6、e7、e8、e9が小さい順に数直線41の上にあるとする。各固有値の位置について理解のために数直線の上に矢印を引いてある。固有値e1についてはeps未満の距離にあるのは、e2、e3である。図４では、理解のためにe1からe2、e3に折れ線42、43を引いた。他の固有値についても同様である。
【００３２】
次に、図５を用いて固有値解析処理の具体的な処理内容について説明する。ここでは、j番目の固有値に対して、自分も含めて距離が基準値eps未満である固有値で最も小さい固有値の番号をfir(j)、最も大きい固有値をend(j)とする。
【００３３】
まず、「ｊ＝１であるかどうかを判定し（ブロック501）、そうであればfir(j)＝１とし（ブロック502）、そうでなければ以下に進む。まずk=j-1とする（ブロック503）。次に｜e(k)-e(j)|<epsであるかどうかを判定し（ブロック504）、そうでなければfir(j)=k-1とし（ブロック505）、そうであれば以下に進む。まず、k=k-1とし（ブロック506）、ブロック501に戻る」という操作をｊ＝１から９まで実行する（ブロック507）。
【００３４】
次に、「j=9であるかどうかを判定し（ブロック508）、そうであればend(j)＝9とし
（ブロック509）、そうでなければ以下に進む。まずk=j+1とする（ブロック510）。次に｜e(k)-e(j)|<epsであるかどうかを判定し（ブロック511）、そうでなければend(j)=k+1とし（ブロック512）、そうであれば以下に進む。まず、k=k-1とし（ブロック513）、ブロック501に戻る」という操作をｊ＝１から９まで実行する（ブロック514）。
【００３５】
次に、図６を用いて固有ベクトル計算優先度決定処理について説明する。
【００３６】
ここではまず、各固有ベクトル計算の優先度を定め、これに基づいて全体の固有ベクトル計算の順序を決定する。ここでは次の規則に従って各固有ベクトルに優先度を定める。（規則１）e(1)の優先度は１とする。（規則２）各固有ベクトルについて、出来る限り高い優先度を付けるようにする。（規則３）距離がeps未満であって、自分よりも小さい固有値の優先度と同じ優先度になってはいけない。以上の規則によってe(1)からe(9)の固有値の優先度を決定する。61は、e(1)の優先度を（規則１）に従って定め、次にe(2)の優先度を（規則２）と（規則３）によって定めた結果を表したものである。62、63、64は（規則２）と（規則３）により残りの固有値について優先度を定めたものである。
【００３７】
このようにして、優先度が決まった後に、「各プロセッサで担当する固有ベクトル計算について、まず優先度の高い固有ベクトル計算から、同じ優先度の固有ベクトル計算がある場合は、固有値の値が小さいものから順に、固有ベクトル計算を行う」という順序を定める。
【００３８】
次に、図７のフローチャートを用いて固有ベクトル計算優先度決定処理の具体的な処理内容について説明する。ここでは、各固有値e(j)に属する固有ベクトル計算の優先度をprio(j)と表示する。
【００３９】
「j=1であるかを判定し（ブロック71）、そうであればprio(j)=1とし（ブロック72）、そうでなければ以下に進む。まず、k=1、m=fir(j)とする（ブロック73）。次にk=prio(m)かどうかを判定し（ブロック74）、そうでなければk=k+1とし（ブロック75）、ブロック74に戻る。そうであれば、以下に進む。まず、m=m+1とし（ブロック76）、次にm=jであるかどうかを判定し（ブロック77）、そうでなければブロック74に戻り、そうであればprio(j)=kとする（ブロック78）」という操作についてｊを１から９まで繰り返す（ブロック79）。
【００４０】
次に、図８を用いて固有ベクトル計算分配処理について説明する。  ここでは固有ベクトル計算優先度決定処理において求められた優先度に基づいて、直交化に関する計算量を図り、各プロセッサの性能比に応じて、固有ベクトル計算を各プロセッサに割り当てる。
【００４１】
まず、各固有ベクトル計算での直交化の計算量81を、直交化の時に用いる他の固有ベクトルの個数、即ち、距離がeps未満の範囲にあり、優先度の低い他の固有値の個数から定める。例えば、優先度が１である固有ベクトル計算は直交化がないので直交化計算量は０である。e(2)は、距離がeps未満の範囲にあり、優先度の低い他の固有値はe(1)のみなので、１である。他の固有値についても同様である。
【００４２】
次に、直交化計算量とプロセッサ性能比がバランスよくなるよう、各プロセッサに固有ベクトル計算を線82に示すように分配する。即ち、（各プロセッサが担当する直交化計算量/ 各プロセッサの性能比）ができるだけ等しくなるように分配する。
【００４３】
次に、図９のフローチャートを用いて固有ベクトル計算分配処理の前半の具体的な処理内容について説明する。ここでは各固有ベクトル計算における直交化計算量を定めるが、これを、固有値e(j)に属する固有ベクトル計算についての直交化計算量をmas(j)と表示する。
【００４４】
「まずmas(j)=0かつm=fir(j)とする（ブロック91)。次にprio(m)<prio(j)かどうかを判定し（ブロック92）、そうでなければm=m+1とし（ブロック93）、ブロック92に戻り、そうであれば以下に進む。まずmas(j)=mas(j)+1とし（ブロック94）、次にm=end(j)であるかどうかを判定し（ブロック95）、そうでなければブロック93に進んでブロック92へ戻り、そうであればｊのループを終了する（ブロック96）」という操作をｊが１から９まで繰り返す（ブロック97）。
【００４５】
次に、図１０を用いて固有ベクトル計算分配処理の後半の具体的な処理内容について説明する。ここでは直交化計算量mas(j)を用いて固有ベクトル計算の各プロセッサへの分配を定めるが、これを、固有値e(j)に属する固有ベクトル計算を担当するプロセッサ番号（本実施例では１、２、３、４）をmem(j)と表示する。また、各プロセッサｋの性能比をpro(k)とする。
【００４６】
まず全てのmas(j)を加えmassumとする（ブロック1001）。次に全てのpro(k)を加えprosumとする（ブロック1002）。そしてj=1とする（ブロック1003）。次に「まずgpro(k)=0とする（ブロック1004）。次にgpro(k)>massum/prosumであるかどうかを判定し（ブロック1005）、そうであればj=j+1とし（ブロック1006）、そうでなければ次に進む。まず、gpro(k)=gpro(k)+(mas(j)/pro(k))とし（ブロック1007）、次にmem(j)=kとし（ブロック1008）、j=j+1とし（ブロック1009）、ブロック1005に戻る」という操作をｋが１から４まで繰り返す（ブロック1010）。
【００４７】
次に、図１１のフローチャートを用いて固有ベクトル計算処理の具体的な処理内容について説明する。ここでは求めるべき固有ベクトルをv(i,j);iは各固有ベクトルの要素のための番号でありi=1,...,n、またjは各固有ベクトルの番号であり、j=1,...,9とする。また、逆反復計算の反復回数は１回とは限らないが、多くの場合１回で済むことから本実施例では反復計算１回である処理の内容について説明する。
【００４８】
prio(j)の最大値をmaxprioとする（ブロック1101）。次に「「プロセッサｋにおいてmem(j)=kであるかどうかを判定し（ブロック1102）、そうでなければ次のjに進み（ブロック1112）、そうであれば以下に進む。「prio(j)=k2であるかどうかを判定し（ブロック1103）、そうでなければ次のk2へ進み（ブロック1104）、そうであれば以下に進む。まずm=fir(j)とし（ブロック1105）、prio(m)<prio(j)であるかどうかを判定し（ブロック1106）、そうでなければブロック1110に進み、そうであれば以下に進む。次にmem(m)=mem(j)であるかどうかを判定し（ブロック1107）、そうでなければv(i,m);i=1,....,nをmem(m)からmem(j)に転送し（ブロック1108）してブロック1109に進み、そうであれば、近似の固有ベクトルv(i,j);i=1,...,nからv(i,m);i=1,...,nの成分を直交化の手法により取り除く（ブロック1109）。次にm=m+1とし（ブロック1110）、m=end(j)+1とし（ブロック1111）、次のjへいく（ブロック1112）」という操作をk2が１からmaxprioまで繰り返す（ブロック1113）」という操作をプロセッサｋにおいてｊが１から９まで繰り返す（ブロック1114）」という操作をｋが１から４まで繰り返す（ブロック1115）。
【００４９】
最後に、図１２のフローチャートを用いて出力処理の具体的な処理内容について説明する。ここでは「プロセッサｋに付随するメモリｋに保有されている固有ベクトルv(i,j);i=1,...,n,（jはmem(j)=kなるjについて）のデータを出力用外部メモリに出力する（ブロック121）」という操作をkが１から４まで繰り返す（ブロック122）。
【００５６】
【発明の効果】
本発明によれば、個々の固有ベクトル計算において必要最小限の直交化を行うことで全体の演算量の削減を図り、各固有ベクトル計算に要する計算量と各プロセッサの性能とに応じて個々の固有ベクトルの計算を個々のプロセッサに割り付けることでプロセッサ全体のロードバランスの調整を可能にし、固有ベクトル計算を高精度かつ高速に実行できる。これによって、さらに、構造解析シミュレーションが高速に処理できる。
【図面の簡単な説明】
【図１】  全体処理装置のブロック図。
【図２】  全体処理のブロック図。
【図３】  入力処理のフローチャート。
【図４】  固有値距離解析処理の説明図。
【図５】  固有値距離解析処理のフローチャート。
【図６】  固有ベクトル計算優先度決定処理の説明図。
【図７】  固有ベクトル計算優先度決定処理のフローチャート。
【図８】  固有ベクトル計算分配処理の説明図。
【図９】  固有ベクトル分配処理前半のフローチャート。
【図１０】  固有ベクトル計算分配処理後半のフローチャート。
【図１１】  固有ベクトル計算処理のフローチャート。
【図１２】  出力処理のフローチャート。
【図１３】  従来の固有値のグループ化の説明図。
【図１４】  従来のプロセッサへの配分の説明図。
【符号の説明】
21…出力用外部メモリ101を用いる入力処理、22…プロセッサ101による固有値距離解析処理、23…プロセッサ101による固有ベクトル計算優先度決定処理、24…プロセッサ101による固有ベクトル計算分配処理、25…プロセッサ104、105、106、107、メモリ108、109、110、111を用いる固有ベクトル計算処理、26…出力用外部メモリ103への出力処理。
【図１】
0003648286.tif 000004
【図２】
0003648286.tif 000005
【図３】
0003648286.tif 000006
【図４】
0003648286.tif 000007
【図５】
0003648286.tif 000008
【図６】
0003648286.tif 000009
【図７】
0003648286.tif 000010
【図８】
0003648286.tif 000011
【図９】
0003648286.tif 000012
【図１０】
0003648286.tif 000013
【図１１】
0003648286.tif 000014
【図１２】
0003648286.tif 000015
【図１３】
0003648286.tif 000016
【図１４】
0003648286.tif 000017