(57)【要約】【課題】反復解法により連立一次方程式を解くメモリ分散型並列計算機において、多様なデータ格納方法に対応して効率的な並列処理行うことが課題である。【解決手段】反復解法の１つであるＭＧＣＲ法の探索ベクトルｓｉを格納する配列Ｗ１は４台のＰＥに分散配置され、係数行列Ａとｓｉを用いて計算された三角行列を格納する配列Ｗ２は各ＰＥに配置される。Ａ×ｓｉのような行列ベクトル積は４台のＰＥにより並列に計算され、Ａの格納方法に依存しない計算は各ＰＥ内で冗長に処理される。行列ベクトル積を反復解法のアルゴリズムから独立させることにより、各種の行列格納方法に対応することが可能になり、処理の汎用性が向上する。また、ＭＧＣＲ法を実装することで高速な反復解法が実現される。【選択図】図１５
2006164307.tif 000002
【特許請求の範囲】
【請求項１】
連立方程式を反復解法により解く並列計算機において、
係数行列と初期ベクトルから生成される探索ベクトルを格納する中間配列を、複数部分に分割して記憶する中間配列記憶手段と、
前記係数行列および探索ベクトルから生成された要素を含み、解ベクトルを求めるために用いられるＭＧＣＲ法の三角行列を記憶する行列記憶手段と、
前記係数行列を用いて行列ベクトル積を並列に計算し、該行列ベクトル積を用いて前記探索ベクトルを更新する並列演算手段と
を備えることを特徴とする並列処理装置。
【請求項２】
前記中間配列記憶手段は、さらに行列ベクトル積の結果を分割して記憶することを特徴とする請求項１記載の並列処理装置。
【請求項３】
前記行列記憶手段は、さらに前記三角行列の次元の定数ベクトルを記憶し、前記並列演算手段は、該三角行列と定数ベクトルとから形成される連立方程式を解き、得られた結果を用いて前記解ベクトルを求めることを特徴とする請求項１記載の並列処理装置。
【請求項４】
連立方程式を反復解法により解く計算機において、
係数行列と初期ベクトルから生成される探索ベクトルを、中間配列として記憶する中間配列記憶手段と、
前記係数行列および探索ベクトルから生成された要素を含み、解ベクトルを求めるために用いられるＭＧＣＲ法の三角行列を記憶する行列記憶手段と、
前記係数行列を用いて行列ベクトル積を計算し、該行列ベクトル積を用いて前記探索ベクトルを更新する演算手段と
を備えることを特徴とする処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、並列計算機を利用した連立一次方程式の反復解法に係り、係数行列を複数のプロセッシング・エレメントに分散配置して、並列処理により解を求める並列処理装置およびその方法に関する。
【背景技術】
【０００２】
物理現象の解析において現れる偏微分方程式の境界値問題や行列の固有値問題を解く時、一般に、大規模なスパース行列（疎行列）を係数行列とする次のような連立一次方程式を解く必要が生じる。
Ａｘ＝ｂ                                                      （１）
ここで、Ａは一般にｎ×ｎの非対称行列、ｘはｎ次元の変数ベクトル、ｂはｎ次元の定数ベクトルである。ｎの値が１００００以上になることも珍しくない。
【０００３】
大規模な連立一次方程式は、気象予測、原子炉設計、半導体の回路解析、航空工学における流体解析、構造物の構造解析等の多くの科学技術計算に用いられる。また、大規模な固有値問題は、構造物の構造解析、回路解析、地球科学における地震予知、原子炉の安全性解析、分子科学における多電子系のエネルギー計算、原子核の構造解析等の分野において、物理系の固有振動を記述するときに現れる。
【０００４】
したがって、（１）式のような大規模な連立一次方程式を効率よく高速に解くことは、科学技術計算の重要な問題の１つである。今日では、計算を高速化するために、複数のプロセッシング・エレメント（ＰＥ）を備えたメモリ分散型の並列計算機が多く用いられている。
【０００５】
計算機を用いて（１）式を解く１つの方法として、ＡをＬＵ分解するガウスの消去法に基づいた直接法がある。しかし、Ａが大きなスパース行列の場合、非零要素が各行に数個しかないこともあり、計算コストや記憶領域の面で無駄が多い。そこで、単純な行列ベクトル積を繰り返して近似解を求める反復解法が多く用いられている。
【０００６】
反復解法の多くはクリロフ（Krylov）部分空間法に帰着される。今、任意のベクトルｒ0 にＡを次々と乗じていくと、ｒ0，Ａｒ0 ，．．．，Ａk-1 ｒ0 のようなベクトル列が生成される。これらの一次独立なベクトルにより張られる空間はクリロフ部分空間と呼ばれ、（１）式の近似解ｘkをこれらのベクトルの一次結合で記述する一群の反復解法はクリロフ部分空間法と呼ばれる。
【０００７】
このクリロフ部分空間法としては、ＣＧ（Conjugate Gradient）法、ＢＣＧ（Bi-Conjugate Gradient ）法、ＣＲ（Conjugate Residuals ）法、ＧＣＲ（Generalized Conjugate Residuals ）法、ＭＧＣＲ（Modified Generalized Conjugate Residuals）法、ＧＭＲＥＳ（Generalized Mainimal RESidual ）法等がある。
【０００８】
ところで、（１）式のスパース行列Ａの形は与えられた問題によって様々であり、その非零要素を集めて配列に格納する方法を問題によって適当に選ぶ必要がある。しかし、データ格納方法が変わればその配置形態も変わるため、スパース行列とベクトルの積を並列に計算するアルゴリズムも変更する必要が生じる。このため、従来の反復解法では、多様なデータ格納方法の中から問題に適したもの選び、それに応じて演算アルゴリズムを個別に作成している。
【０００９】
したがって、同じ反復解法を用いていても、問題が変わる度に新しくプログラムを作成しなければならず、汎用性に欠けるという問題がある。また、この方法ではプログラマの負担も大きくなる。そこで、連立一次方程式を並列計算機で解く際に、各種のデータ格納方法をサポートしつつ、反復解法を実現することが望まれる。
【００１０】
また、反復解法で繰り返し現れる行列ベクトル積の並列演算を行う際、ＰＥ間のデータ転送をできるだけ少なくして、演算を高速化することが重要である。
【発明の開示】
【発明が解決しようとする課題】
【００１１】
本発明は、クリロフ部分空間を利用した反復解法により連立一次方程式を解くメモリ分散型並列計算機において、多様なデータ格納方法に対応して効率的な並列処理を行う並列処理装置およびその方法を提供することを目的とする。
【課題を解決するための手段】
【００１２】
図１は、本発明の並列処理装置の原理図である。図１の並列処理装置は、スパース行列を係数行列とする連立方程式を反復解法により解く並列計算機に設けられ、反復手段１、並列演算手段２、および中間配列記憶手段３を備える。
【００１３】
並列演算手段２は、上記係数行列を複数部分に分割して格納し、係数行列の格納形式に依存する計算を並列に行う。
中間配列記憶手段３は、並列演算手段２による計算結果を中間配列として、分割して記憶する。
【００１４】
反復手段１は、上記中間配列のデータを用いて、上記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する。
反復手段１は、例えば、上記並列計算機が有する各ＰＥに備えられ、与えられた初期ベクトルを用いて、行列ベクトル積を利用した反復処理を開始する。このとき、初期ベクトルに基づいて入力ベクトルを生成し、並列演算手段２に与える。
【００１５】
並列演算手段２は、例えば、複数のＰＥに対応し、係数行列の格納形式に応じて係数行列と入力ベクトルの積を並列に計算して、中間ベクトルを生成する。このような行列ベクトル積の演算は、反復解法の各繰り返し段階において、少なくとも１回以上必要になる。
【００１６】
生成された中間ベクトルは、中間配列記憶手段３内の中間配列に格納される。中間配列記憶手段３は、例えば、複数のＰＥの主記憶に対応し、中間ベクトルはこれらの主記憶に分割されて格納される。
【００１７】
反復手段１は、中間ベクトルを用いて次の入力ベクトルを生成し、それを並列演算手段２に与える。そして、このような処理の繰り返しにより解が収束すると、反復処理を終了して得られた解ベクトルを出力する。
【００１８】
中間配列を用いることにより、反復手段１は、係数行列の格納形式とは無関係に次の入力ベクトルを生成し、収束判定を行うことができる。したがって、反復手段１の処理を係数行列のデータ格納方法から完全に独立させることができる。これにより、特定のデータ格納方法を前提としない反復処理の実装が可能になり、反復解法の汎用性が高まる。
【００１９】
データ格納方法に依存する行列ベクトル積等の計算は、必要に応じて並列演算手段２に依頼し、反復手段１は与えられた反復解法の中核アルゴリズムのみを実行すればよい。一方、並列演算手段２は、採用された反復解法の種類とは無関係に、与えられた行列ベクトル積等の単純な計算を繰り返すだけである。こうして、効率のよい反復処理が実現される。
【００２０】
データ格納方法に依存する処理かどうかは、例えば、各処理の処理コードに付加されたパラメータの値等により識別することができる。
このような並列処理方法では、データ格納方法が変わっても、反復解法の中核アルゴリズムをプログラミングし直す必要がないので、プログラマの作業が大幅に削減される。
【００２１】
例えば、図１の反復手段は、実施形態の図２における個々のＰＥ３３または図１７における反復部５１に対応し、並列演算手段２は、図２における複数のＰＥ３３または図１７における行列ベクトル積演算部５２に対応する。
【発明の効果】
【００２２】
本発明によれば、スパース行列の連立一次方程式を解くメモリ分散型並列計算機において、同一の反復解法のアルゴリズムを、多様なスパース行列の格納法に適用することが可能になる。したがって、反復解法の汎用性が向上し、プログラマの負担が軽減される。
【００２３】
また、スパース行列をバンド行列とみなして処理することで、行列ベクトル積の演算に必要なデータ転送量が削減され、並列処理の台数効果が向上する。特に、反復解法として他の解法より高速なＭＧＣＲ法を用いた場合、さらに効率のよい並列処理が実現される。
【発明を実施するための最良の形態】
【００２４】
以下、図面を参照しながら、本発明の実施の形態を詳細に説明する。
本発明においては、反復解法において必要となる行列ベクトル積等の処理を行うために、ベクトルを格納する１次元配列を各ＰＥ（プロセッサ）で均等に分割し、この配列を中間インタフェースとして利用して、中核の計算アルゴリズムをスパース行列のデータ格納方法から独立させる。つまり、反復解法の各段階に現れる中間ベクトルを各ＰＥで均等に分割した中間領域に格納する。
【００２５】
これにより、データ格納方法に依存する行列ベクトル積等の処理が反復解法の中核アルゴリズムから分離され、多様なデータ格納方法に対処できるようになる。
また、代表的な行列のデータ格納法である対角形式格納法またはエルパック（Ellpack ）形式格納法を用いて、行列ベクトル積を並列に計算することができる。このとき、スパース行列を疑似的なバンド行列とみなしてそのバンド幅を求め、これを利用して行列ベクトル積を計算すれば、各ＰＥは近くの限られたＰＥとのみデータ転送を行えばよくなる。したがって、転送回数が削減され、色々なバンド幅のスパース行列について、行列ベクトル積の演算が効率よく実行される。
【００２６】
さらに、一般スパース行列の連立１次方程式の反復解法の１つであるＭＧＣＲ法を例に取り、本発明の反復処理をメモリ分散型のスーパーコンピュータＶＰＰシリーズ向けに構成した。この形態によれば、各種のデータ格納方法に対してＭＧＣＲ法を適用し、高速な計算処理を実現することが可能になる。
【００２７】
まず、図２および図３を参照しながら、実施形態で用いる並列計算機の構成を説明する。図２は、並列計算機の概略構成図である。図２の並列計算機は、入出力装置１１と、ＩＯバス１２と、クロスバ方式により互いに接続された複数のＰＥ１３とを備える。
【００２８】
入出力装置１１は、各ＰＥ１３との間でデータの入出力を行うための装置で、例えば、キーボードなどの入力機器とディスプレイやプリンタなどの出力機器を備えた計算機端末である。入出力装置１１は、ＩＯバス１２により各ＰＥ１３と接続されている。
【００２９】
図３は、ＰＥ１３の構成図である。図５のＰＥ１３は、主記憶２１、ＩＯポート２２、クロスバスイッチ２３、データ転送制御ユニット２４、メモリ制御ユニット２５、スカラユニット３８、およびベクトルユニット３９を備える。
【００３０】
主記憶２１は、与えられた問題の行列要素や計算の中間結果などを格納する。ＩＯポート２２はＩＯバス１２に接続され、入出力装置１１との間で入出力データの転送を行う。クロスバスイッチ２３は、主記憶２１内のデータを他のＰＥ１３に転送する際のデータの切り換え等を行う。データ転送制御ユニット２４は、クロスバスイッチ２３を介して他のＰＥ１３とのデータ転送および同期制御を行う。
【００３１】
メモリ制御ユニット２５は、データ転送制御ユニット２４、スカラユニット３８、およびベクトルユニット３９から発生するアクセス要求を受け取り、主記憶２１に対するアクセスを制御する。
【００３２】
スカラユニット３８は、キャッシュメモリ２６、ＧＰＲ（汎用レジスタ）／ＦＰＲ（浮動小数点レジスタ）２７、スカラ演算器２８を備え、スカラ演算を実行する。
ベクトルユニット３９は、マスクレジスタ２９とベクトルレジスタ３０の他に、ロードパイプライン３１、ストアパイプライン３２、マスクパイプライン３３、３４、乗算パイプライン３５、加算／論理演算パイプライン３６、および除算パイプライン３７を備える。
【００３３】
これらの各パイプラインは、それぞれ複数要素を同時に処理することができ、パイプライン３５、３６、３７のうち２本は同時に動作できる。また、他のパイプラインについては、すべて同時に動作可能である。マスクパイプライン３３は総和／検索処理用に使用され、マスクパイプライン３４は論理演算処理用に使用される。
【００３４】
ベクトルユニット３９は、これらのパイプラインを用いて複数のベクトル命令を並列に実行することができ、高速な行列演算が可能である。
次に、図４および図５を参照しながら、本発明の反復処理について説明する。図４は、多様なデータ格納方法に対応可能な反復処理のフローチャートである。図４において処理が開始されると、並列計算機の各ＰＥは、まずあらかじめ決められた反復法による処理コードを取り出す（ステップＳ１）。
【００３５】
そして、そのコードに記述された計算がデータ格納法に依存するかどうかを判定し（ステップＳ２）、データ格納法に依存しなければ各ＰＥで同じ処理を実行する（ステップＳ３）。データ格納法に依存するかどうかは、例えば、あらかじめ処理コード内に書き込まれたパラメータ値により判断する。
【００３６】
次に、解が十分に収束するなどの反復終了条件が満たされたかどうかを判定し（ステップＳ５）、それが満たされなければステップＳ１以降の処理を繰り返す。
また、ステップＳ２において計算がデータ格納法に依存する場合は、他のＰＥと並列処理を行って、計算結果を中間配列に格納する（ステップＳ４）。この中間配列は、各種のデータ格納法に使用できるように標準的に用意された領域で、各ＰＥにより分割されている。例えば、ＰＥ１、ＰＥ２、ＰＥ３、ＰＥ４の４台のＰＥを備える場合は、中間配列は図５に示すようになる。
【００３７】
その後、各ＰＥはステップＳ５以降の処理を繰り返し、反復終了条件が満たされれば処理を終了する。
このような反復処理によれば、例えば行列ベクトル積に代表されるデータ格納法依存性の高い並列演算の結果を中間配列に格納して、反復法における他の計算処理をデータ格納法と独立に行うことができる。これにより、各種スパース行列の格納法を用いて、スパース行列の連立１次方程式を容易に解くことが可能になる。
【００３８】
次に、図６から図１２までを参照しながら、具体的なデータ格納法に基づく効率のよい行列ベクトル積の演算方法を説明する。
図６は、スパース行列の対角形式格納法を示している。ここでは、ｎ次のスパース行列Ａをバンド行列の１種とみなし、対角線方向に並んだ非零要素から成る対角ベクトルを２次元の格納配列４１に格納する。
【００３９】
Ａの対角ベクトル部分を拡大すると図７に示すようになり、対角ベクトル間には対角方向の要素がすべて０の部分が存在する。このような部分は格納配列４１に格納しなくてもよいので、対角方向に集中して非零要素が存在するスパース行列に特に適した格納法であると言える。
【００４０】
図６および図７に示すように、Ａをバンド行列とみなしたときの上バンド幅ｂａｎｄｗｕと下バンド幅ｂａｎｄｗｌは、それぞれ対角要素から対角ベクトルまでの距離の最大値として求められる。このとき、格納配列４１の幅ｗは最大ｂａｎｄｗｕ＋ｂａｎｄｗｌ＋１となる。
【００４１】
また、図６において対角ベクトルの長さを揃えるために、斜線部分に０の要素を付け加えて格納配列４１に格納しておく。これにより、Ａの対角ベクトルを用いた演算を均一に行うことができる。
【００４２】
配列４１の１次元目（行の次元）を各ＰＥに均等に分割して格納し、さらにｎ次元のベクトルＸ、Ｙを格納する配列を各ＰＥに分割配置することにより、次式の行列ベクトル積を並列に計算することができる。
Ｙ＝ＡＸ                                                    （１１）
例えば、４台のＰＥで並列演算を行う場合は、Ａの配列４１、Ｘの配列４２、Ｙの配列４３は、それぞれ図８に示すように分割される。このうち、演算結果を格納する配列４３は図５の中間配列に相当する。
【００４３】
ところで、図６のＡは一定のバンド領域にのみ非零要素を持つスパース行列であるため、（１１）式の演算において、各ＰＥは必ずしもＸのすべての要素を持っている必要はない。例えば図８のＰＥ２が受け持つ行列ベクトル積の有効部分は、図９に示すようになる。
【００４４】
図９においてＰＥ２に割り当てられたＡの行数をｎ２とすると、その非零要素が存在する領域はＡの斜線部分である。この斜線部分に該当する列の範囲は、ｎ２本の行と同じ番号を持つｎ２本の列と、その左側のｂａｎｄｗｌの幅の各列と、右側のｂａｎｄｗｕの幅の各列となる。したがって、Ｘの要素のうち、これらの列に対応する（ｎ２＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）個の要素のみが乗算の結果に寄与する。他の要素には０が乗算されるので、その結果もまた０となる。
【００４５】
Ｘを格納する配列４２は、図８に示すように各ＰＥに分散されているため、ＰＥ２は乗算を行うために、少なくともｂａｎｄｗｌ個の要素をＰＥ１からコピーし、少なくともｂａｎｄｗｕ個の要素をＰＥ３からコピーする必要がある。他のＰＥについても同様のコピー処理が必要になる。
【００４６】
そこで、あらかじめＰＥ毎にワーク用の一次元配列Ｗを用意しておき、この配列に必要なＸの要素をコピーすることにする。例えばＰＥ２の場合、Ｗの長さとしては、最低（ｎ２＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）だけあれば十分である。しかし、ここではＸの各要素の論理的な添え字を、そのままＷの対応する要素の添え字として使用できるように、各ＰＥのＷの長さを（ｎ＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）に統一することにする。
【００４７】
こうして、図１０に示すようなＷが各ＰＥに用意される。ここで、ＰＥの台数をＰ、ｉ番目のＰＥに割り当てられたＡの行数をｎｉ（ｉ＝１，．．．，Ｐ）とすると、ｉ番目のＰＥでは、Ｘの必要な要素が図１０に示すようにＷの斜線部分にコピーされる。ここで、ｎｉはおおよそｎ／Ｐとなる。Ｗの斜線部分以外の領域は実際の処理には用いられないので、どんな値を格納していてもかまわない。
【００４８】
Ｘの斜線部分の要素のうち、元々自分が持っているｎｉ個の要素はそのまま主記憶２１上でコピーすることができ、その上下のｂａｎｄｗｌ個とｂａｎｄｗｕ個の要素は、クロスバスイッチ２３を介して隣接ＰＥまたは他の近くのＰＥからコピーする。配列４１に格納された対角ベクトルの数がＡの次元ｎに比べて十分に小さい場合は、高々隣接ＰＥとの通信だけで必要な要素をコピーすることができる。
【００４９】
このように、図６のような対角方向に非零要素が集中しているスパース行列の場合は、特定のＰＥ間の通信だけで行列ベクトル積を実行でき、通信に伴うオーバヘッドを大幅に削減することができる。
【００５０】
また、１番目のＰＥは、Ｗの上端部の長さｂａｎｄｗｌの部分にダミーデータを格納して乗算を行い、Ｐ番目のＰＥは、下端部の長さｂａｎｄｗｕの部分にダミーデータを格納して乗算を行う。これらのダミーデータには、図６の斜線部分の要素である０が乗算されるので、ダミーデータの値が行列ベクトル積の結果に影響を与えることはない。しかし、ダミーデータを用いることで、１番目とＰ番目のＰＥは他のＰＥと同様の乗算を行うことができる。
【００５１】
このような乗算方法によれば、各ＰＥ内の演算はベクトル演算器を用いて効率よく処理することができる。また、各ＰＥは演算の途中でベクトル長を変更する必要もない。
図１１は、バンド幅を利用したスパース行列の行列ベクトル積演算のフローチャートである。図１１において処理が開始されると、各ＰＥは、まずスパース行列をバンド行列とみなして、その上バンド幅と下バンド幅を求める（ステップＳ１１）。次に、分割されたベクトルＸの要素のうち、自分が保持する部分（図１０のｎｉ個の要素）をワーク用配列Ｗに並列にコピーする（ステップＳ１２）。
【００５２】
次に、各ＰＥは、必要な下バンド幅の長さの部分（図１０のｂａｎｄｗｌ個の要素）をＷに並列にコピーする（ステップＳ１３）。また、必要な上バンド幅の長さの部分（図１０のｂａｎｄｗｕ個の要素）をＷに並列にコピーする（ステップＳ１４）。
【００５３】
そして、行列ベクトル積を並列に計算して、計算結果をベクトルＹの各ＰＥに割り当てられた部分に格納し（ステップＳ１５）、処理を終了する。
この方法によれば、ステップＳ１３およびＳ１４において限られたＰＥとのみ通信を行えばよいので、データ転送量が少なくて済み、行列ベクトル積演算が高速化される。
【００５４】
このようなバンド幅を利用した行列ベクトル積の計算方法は、対角形式格納法に限らず、他の任意のデータ格納法にも適用可能である。他のデータ格納法の１つとして、エルパック形式格納法がある。
【００５５】
エルパック形式格納法では、スパース行列の各行ベクトル毎に非零要素を圧縮した形で、行列要素を格納する。このとき、各行内の非零要素に対応する元のスパース行列の列番号は、格納配列ｉｃｏｆに格納されている。この列番号の値を利用すれば、バンド幅を求めることができる。
【００５６】
図１２は、配列ｉｃｏｆの例を示している。図１２において、左のスパース行列は１００００×１００００の行列で、対角要素と他の対角方向の線分上に非零要素を持っている。また、ｉｃｏｆは、スパース行列の第ｉ行に含まれる非零要素の列番号ｉｃｏｆ（ｉ，＊）を、各行番号ｉ毎に格納している。ここで、＊は第ｉ行内の１つ以上の非零要素に付けられた識別番号である。
【００５７】
スパース行列のすべての非零要素のうち、対角要素から最も離れたものとその行の対角要素との距離を求めれば、バンド幅が得られる。
例えば、ｉ＝９０００の行内には、５０００番目、６０００番目、および９０００番目の各列に非零要素があり、これらの列番号がｉｃｏｆ（９０００，＊）としてｉｃｏｆに格納されている。このうち、９０００番目の列に対角要素があり、それから左に最も離れた要素は５０００番目の列にある。したがって、これらの要素間の距離は−（５０００−９０００）＝４０００となる。
【００５８】
また、この５０００番目の列の要素はバンドの最も外側に位置しているため、この場合の下バンド幅ｂａｎｄｗｌは４０００となる。上バンド幅ｂａｎｄｗｕもｉｃｏｆから同様にして求められる。
【００５９】
ｉｃｏｆ（ｉ，＊）を用いたバンド幅の計算式を一般化すると、次式のようになる。
【００６０】
【数１】
2006164307.tif 000003
【００６１】
（１２）、（１３）式において、＊に関するＭａｘは、第ｉ行内の＊番目の要素と対角要素との距離の最大値を表し、ｉに関するＭａｘは、各行の＊に関する最大値のうちの最大値を表す。
【００６２】
（１２）式で、（ｉｃｏｆ（ｉ，＊）−ｉ）の前に負の符号が付いているのは、左に位置する要素ほど対角要素までの距離を大きく評価するためである。この場合、対角要素より右に位置する要素までの距離は負の値を持つことになる。逆に、（１３）式では、右に位置する要素ほど対角要素までの距離が大きくなり、対角要素より左に位置する要素までの距離は負の値を持つ。
【００６３】
こうして、エルパック形式格納法においても上下バンド幅を計算することが可能になり、図１１の演算処理を適用することができる。
次に、図１３から図１７までを参照しながら、最も高速な反復解法の１つであるＭＧＣＲ法の実施形態について説明する。ＭＧＣＲ法は、各繰り返し段階における演算回数がＧＭＲＥＳ法より少ないため、この方法より高速であることが知られている（Z. Leyk,  “Modified generalized conjugate residuals for nonsymmetric systems of linear equations ” in “Proceeding of 6th Biennial Conference on Computational Techniques and Applications : CTAC93 ”, D. Stewart, H. Cardner and D. Singleton, eds., World Scientific, 1994, pp. 338-344）。
【００６４】
ＧＭＲＥＳ法は反復解法の中でもかなり速いアルゴリズムである。それより速いＭＧＣＲ法を実装することで、従来にない高速な反復処理装置を実現することができる。ここでは、一例としてスーパーコンピュータＶＰＰ５００を対象とした実装方法を説明する。
【００６５】
（１）式でｘ＝ｕ，ｂ＝ｆとおき直し、与えられた連立一次方程式を、
Ａｕ＝ｆ                                                    （１４）
と書くことにすると、ベクトルｕの近似解を求めるＭＧＣＲ法のアルゴリズムは次のようになる。
【００６６】
【数２】
2006164307.tif 000004
【００６７】
ここで、（１５）式のｕ0 は任意に与えられた初期ベクトルであり、｜ｒ1 ｜はベクトルｒ1の大きさを表す。また、ｓ1 ＝ｒ1 は１番目の探索ベクトルとなる。
（１６）〜（２３）式のループ処理では、ｉ＝１，．．．，ｋについてｋ個の探索ベクトルｓi+1 を生成し、（２３）式でｄi+1 が、与えられた収束判定値ε以内に収束すれば処理を終了する。そうでなければ、（２４）〜（２５）式によりｕ0を更新して、（１５）式からの始まる外側のループ処理を繰り返す。ｋの値は任意に指定可能であるが、通常は１０〜１００の範囲に設定すればよい。
【００６８】
（１６）、（２１）式において、（ｖ１，ｖ２）のような表記はベクトルｖ１とｖ２の内積を表し、（１９）、（２３）、（２４）式において、ｓｏｌｖｅ  Ｈｃ＝αのような表記は、連立方程式Ｈｃ＝αを解いてｃを求めることを意味する。
【００６９】
Ｈk は、次式で与えられるｋ次の上三角行列である。
【００７０】
【数３】
2006164307.tif 000005
【００７１】
Ｈk の要素βj,i-1 （ｉ＝１，．．．，ｋ，  ｊ＝１，．．．，ｉ）は、（１６）式または（１７）、（１８）式により、Ａおよびｓiから計算される。Ｈi-1 ，Ｈi についても同様である。
【００７２】
αk は、αk ＝（α1 ′，．．．，αk′）T のようなｋ次元の定数ベクトルである。このベクトルの要素αi ′は、（２１）式によりｓ1とｓi+1 から計算される。αi-1 ，αi についても同様である。
【００７３】
また、（２５）式におけるＮk は次式で与えられる。
Ｎk ＝（ｓ1 ，ｓ2 ，．．．，ｓk）                          （２８）
（１９）式のＮi-1 、（２３）式のＮi も同様である。
【００７４】
（１７）式の
2006164307.tif 000006
は実際には計算されず、βk,k-1 とαk ′は次式のようにして求められる。
【００７５】
【数４】
2006164307.tif 000007
【００７６】
このＭＧＣＲ法のアルゴリズムを並列計算機上に実装するとき、（１５）式のＡ×ｕ0 や（１６）、（１７）式のＡ×ｓi のようなスパース行列とベクトルの積を格納する領域が必要になる。そこで、各ＰＥにより分割された２次元の中間配列Ｗ１を用意し、Ｙ＝ＡＸの形の行列ベクトル積のＸとＹを格納する。また、（１７）式により生成されるベクトルｓi（ｉ＝１，．．．，ｋ）もＷ１に格納しておく。
【００７７】
４台のＰＥに分散配置されたＷ１は図１３に示すようになる。図１３において、Ｗ１はｎ×（ｋ＋２）の大きさを持ち、第１列にはＸが格納され、第２列にＹが格納される。また、残りのｋ本の列には各ｓi が格納される。第１列または第２列は、（１５）式のｒ1を格納する時にも使用される。
【００７８】
図１４は、各ＰＥ毎に設けられる格納配列を示している。図１４の配列Ｗ２はｋ×（ｋ＋１）の大きさを持ち、最初のｋ本の列にＨk が格納され、最後の列にαk が格納される。実際には、Ｈk用の領域にはＨi （ｉ＝１，．．．，ｋ）が順次格納されて処理に用いられ、αk 用の領域にはαi（ｉ＝１，．．．，ｋ）が順次格納されて処理に用いられる。
【００７９】
この他に、スパース行列Ａを格納する配列が各ＰＥに分散配置され、図１１のような行列ベクトル積の演算を行うために、図１０のようなワーク用配列が各ＰＥに１つ設けられる。スパース行列Ａは、例えば対角形式格納法やエルパック形式格納法により格納されるが、もちろん、他の任意のデータ格納法を採用することもできる。
【００８０】
図１５は、これらの格納配列を用いたＭＧＣＲ法による処理の概要を示している。図１５において、Ｗ１は４台のＰＥに分割されている。
これらのＰＥ１、ＰＥ２、ＰＥ３、ＰＥ４はＷ１を利用して、必要に応じてＡ×ｕ0 またはＡ×ｓi を並列に計算する。その後、各ＰＥはＷ２を利用して他の計算を同時に行い、ｕ0を更新していく。そして、（１９）式または（２３）式の条件が成り立てば、処理を終了する。
【００８１】
図１６は、このようなＭＧＣＲ法による処理のフローチャートである。図１６において処理が開始されると、並列計算機は、まず与えられた初期ベクトルｕ0 を用いた並列処理によりＡ×ｕ0 を計算し、（１５）式のｒ1とｄ1 を求める（ステップＳ２１）。Ａ×ｕ0 の行列ベクトル積は図１１の方法で実行される。
【００８２】
次に、ｉ＝１とおいて、ＭＧＣＲ法による処理を開始し（ステップＳ２２）、図１１の方法によりＡ×ｓi を並列に計算する（ステップＳ２３）。次に、各ＰＥでＡ×ｓiの結果を用いて、同じ（１６）、（１７）、（１８）式の計算を冗長に行い、βi,i-1 を求める（ステップＳ２４）。そして、βi,i-1が０かどうかを調べる（ステップＳ２５）。
【００８３】
βi,i-1 が０でなければ、各ＰＥでβi,i-1 の値を用いて、同じ（２０）、（２１）、（２２）式の計算を冗長に行い、ｓi+1およびｄi+1 を求める（ステップＳ２６）。そして、ｄi+1 がε以下になったかどうかを調べ、収束判定を行う（ステップＳ２７）。
【００８４】
ｄi+1 がεより大きければ、次にｉがｋに達したかどうかを調べる（ステップＳ２８）。そして、ｉ＝ｋでなければｉ＝ｉ＋１とおいて（ステップＳ２９）、ステップＳ２３以降の処理を繰り返す。また、ステップＳ２８においてｉ＝ｋとなった場合は、各ＰＥで（２４）、（２５）、（２６）式の再スタート処理を行ってｕ0を更新し（ステップＳ３０）、ステップＳ２１以降の処理を繰り返す。
【００８５】
そして、ステップＳ２５においてβi,i-1 ＝０となった場合は、各ＰＥは（１９）式の第１の終了処理を行って近似解ｕi-1を求め（ステップＳ３１）、処理を終了する。また、ステップＳ２７においてｄi+1 がε以下に収束した場合は、各ＰＥは（２３）式の第２の終了処理を行って近似解ｕiを求め（ステップＳ３２）、処理を終了する。
【００８６】
上述の第１の終了処理、第２の終了処理、および再スタート処理において、最大ｋ次の連立一次方程式を解く必要がある。しかし、ｋは高々１００の程度で、多くの場合ｎに比べてはるかに小さいので、各ＰＥ内で逐次的に解いてもそれほど時間はかからない。
【００８７】
以上のようなＭＧＣＲ法を実装した並列計算機を模式的に表すと、図１７のように書くことができる。図１７において、反復部５１は、スパース行列Ａのデータ格納法に依存しない計算処理や判定処理を行い、行列ベクトル積演算部５２は、Ａ×ｕ0 やＡ×ｓi の行列ベクトル積の計算処理を行う。
【００８８】
これらの行列ベクトル積の計算はＡの格納方法に依存するため、行列ベクトル積演算部５２に処理コードを与えて、入力ベクトルＸや出力ベクトルＹの格納場所等を指示する。このような制御用の処理コードは、コンパイラまたはプログラマにより作成される。
【００８９】
図１７のような構成を取ることにより、行列ベクトル積の演算をＭＧＣＲ法から分離して独立に実行することができ、各種のスパース行列格納法に対して、ＭＧＣＲ法を適用することが可能になる。
【００９０】
本発明で用いる反復解法はＭＧＣＲ法に限られるわけではなく、クリロフ部分空間を利用した任意の解法を用いることができる。また、図１７の構成において、反復部５１の処理をいずれかのＰＥに代表して行わせることも可能である。
【図面の簡単な説明】
【００９１】
【図１】本発明の原理図である。
【図２】並列計算機の構成図である。
【図３】プロセッシング・エレメントの構成図である。
【図４】本発明の反復処理のフローチャートである。
【図５】中間配列を示す図である。
【図６】スパース行列の対角形式格納法を示す図である。
【図７】スパース行列の対角ベクトルの拡大図である。
【図８】行列ベクトル積の並列処理を示す図である。
【図９】ＰＥ２による行列ベクトル積を示す図である。
【図１０】ワーク用配列を示す図である。
【図１１】バンド幅を利用した行列ベクトル積演算のフローチャートである。
【図１２】エルパック形式格納法における列番号配列を示す図である。
【図１３】ＭＧＣＲ法における中間配列を示す図である。
【図１４】各ＰＥが保持する格納配列を示す図である。
【図１５】ＭＧＣＲ法による処理の概要を示す図である。
【図１６】ＭＧＣＲ法による処理のフローチャートである。
【図１７】ＭＧＣＲ法の装置構成図である。
【符号の説明】
【００９２】
１  反復手段
２  第１の並列処理手段
３  中間配列記憶手段
１１  入出力装置
１２  ＩＯバス
１３  ＰＥ
２１  主記憶
２２  ＩＯポート
２３  クロスバスイッチ
２４  データ転送制御ユニット
２５  メモリ制御ユニット
２６  キャッシュメモリ
２７  ＧＰＲ／ＦＰＲ
２８  スカラ演算器
２９  マスクレジスタ
３０  ベクトルレジスタ
３１  ロードパイプライン
３２  ストアパイプライン
３３、３４  マスクパイプライン
３５  乗算パイプライン
３６  加算／論理演算パイプライン
３７  除算パイプライン
３８  スカラユニット
３９  ベクトルユニット
４１、４２、４３、Ｗ、Ｗ１、Ｗ２、ｉｃｏｆ  格納配列
５１  反復部
５２  行列ベクトル積演算部
【図１】
2006164307.tif 000008
【図２】
2006164307.tif 000009
【図３】
2006164307.tif 000010
【図４】
2006164307.tif 000011
【図５】
2006164307.tif 000012
【図６】
2006164307.tif 000013
【図７】
2006164307.tif 000014
【図８】
2006164307.tif 000015
【図９】
2006164307.tif 000016
【図１０】
2006164307.tif 000017
【図１１】
2006164307.tif 000018
【図１２】
2006164307.tif 000019
【図１３】
2006164307.tif 000020
【図１４】
2006164307.tif 000021
【図１５】
2006164307.tif 000022
【図１６】
2006164307.tif 000023
【図１７】
2006164307.tif 000024