JP3808925B2 - 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 - Google Patents
多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 Download PDFInfo
- Publication number
- JP3808925B2 JP3808925B2 JP01610596A JP1610596A JP3808925B2 JP 3808925 B2 JP3808925 B2 JP 3808925B2 JP 01610596 A JP01610596 A JP 01610596A JP 1610596 A JP1610596 A JP 1610596A JP 3808925 B2 JP3808925 B2 JP 3808925B2
- Authority
- JP
- Japan
- Prior art keywords
- coefficient matrix
- parallel
- processing
- vector
- storage format
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Complex Calculations (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、並列計算機を利用した連立一次方程式の反復解法に係り、係数行列を複数のプロセッシング・エレメントに分散配置して、並列処理により解を求める並列処理装置およびその方法に関する。
【０００２】
【従来の技術とその問題点】
物理現象の解析において現れる偏微分方程式の境界値問題や行列の固有値問題を解く時、一般に、大規模なスパース行列（疎行列）を係数行列とする次のような連立一次方程式を解く必要が生じる。
Ａｘ＝ｂ （１）
ここで、Ａは一般にｎ×ｎの非対称行列、ｘはｎ次元の変数ベクトル、ｂはｎ次元の定数ベクトルである。ｎの値が１００００以上になることも珍しくない。
【０００３】
大規模な連立一次方程式は、気象予測、原子炉設計、半導体の回路解析、航空工学における流体解析、構造物の構造解析等の多くの科学技術計算に用いられる。また、大規模な固有値問題は、構造物の構造解析、回路解析、地球科学における地震予知、原子炉の安全性解析、分子科学における多電子系のエネルギー計算、原子核の構造解析等の分野において、物理系の固有振動を記述するときに現れる。
【０００４】
したがって、（１）式のような大規模な連立一次方程式を効率よく高速に解くことは、科学技術計算の重要な問題の１つである。今日では、計算を高速化するために、複数のプロセッシング・エレメント（ＰＥ）を備えたメモリ分散型の並列計算機が多く用いられている。
【０００５】
計算機を用いて（１）式を解く１つの方法として、ＡをＬＵ分解するガウスの消去法に基づいた直接法がある。しかし、Ａが大きなスパース行列の場合、非零要素が各行に数個しかないこともあり、計算コストや記憶領域の面で無駄が多い。そこで、単純な行列ベクトル積を繰り返して近似解を求める反復解法が多く用いられている。
【０００６】
反復解法の多くはクリロフ（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 ）形式格納法を用いて、行列ベクトル積を並列に計算することができる。このとき、スパース行列を疑似的なバンド行列とみなしてそのバンド幅を求め、これを利用して行列ベクトル積を計算すれば、各ＰＥは近くの限られたＰＥとのみデータ転送を行えばよくなる。したがって、転送回数が削減され、色々なバンド幅のスパース行列について、行列ベクトル積の演算が効率よく実行される。
【００２４】
さらに、一般スパース行列の連立１次方程式の反復解法の１つであるＭＧＣＲ法を例に取り、本発明の反復処理をメモリ分散型のスーパーコンピュータＶＰＰシリーズ向けに構成した。この形態によれば、各種のデータ格納方法に対してＭＧＣＲ法を適用し、高速な計算処理を実現することが可能になる。
【００２５】
まず、図２および図３を参照しながら、実施形態で用いる並列計算機の構成を説明する。図２は、並列計算機の概略構成図である。図２の並列計算機は、入出力装置１１と、ＩＯバス１２と、クロスバ方式により互いに接続された複数のＰＥ１３とを備える。
【００２６】
入出力装置１１は、各ＰＥ１３との間でデータの入出力を行うための装置で、例えば、キーボードなどの入力機器とディスプレイやプリンタなどの出力機器を備えた計算機端末である。入出力装置１１は、ＩＯバス１２により各ＰＥ１３と接続されている。
【００２７】
図３は、ＰＥ１３の構成図である。図５のＰＥ１３は、主記憶２１、ＩＯポート２２、クロスバスイッチ２３、データ転送制御ユニット２４、メモリ制御ユニット２５、スカラユニット３８、およびベクトルユニット３９を備える。
【００２８】
主記憶２１は、与えられた問題の行列要素や計算の中間結果などを格納する。ＩＯポート２２はＩＯバス１２に接続され、入出力装置１１との間で入出力データの転送を行う。クロスバスイッチ２３は、主記憶２１内のデータを他のＰＥ１３に転送する際のデータの切り換え等を行う。データ転送制御ユニット２４は、クロスバスイッチ２３を介して他のＰＥ１３とのデータ転送および同期制御を行う。
【００２９】
メモリ制御ユニット２５は、データ転送制御ユニット２４、スカラユニット３８、およびベクトルユニット３９から発生するアクセス要求を受け取り、主記憶２１に対するアクセスを制御する。
【００３０】
スカラユニット３８は、キャッシュメモリ２６、ＧＰＲ（汎用レジスタ）／ＦＰＲ（浮動小数点レジスタ）２７、スカラ演算器２８を備え、スカラ演算を実行する。
【００３１】
ベクトルユニット３９は、マスクレジスタ２９とベクトルレジスタ３０の他に、ロードパイプライン３１、ストアパイプライン３２、マスクパイプライン３３、３４、乗算パイプライン３５、加算／論理演算パイプライン３６、および除算パイプライン３７を備える。
【００３２】
これらの各パイプラインは、それぞれ複数要素を同時に処理することができ、パイプライン３５、３６、３７のうち２本は同時に動作できる。また、他のパイプラインについては、すべて同時に動作可能である。マスクパイプライン３３は総和／検索処理用に使用され、マスクパイプライン３４は論理演算処理用に使用される。
【００３３】
ベクトルユニット３９は、これらのパイプラインを用いて複数のベクトル命令を並列に実行することができ、高速な行列演算が可能である。
次に、図４および図５を参照しながら、本発明の反復処理について説明する。図４は、多様なデータ格納方法に対応可能な反復処理のフローチャートである。図４において処理が開始されると、並列計算機の各ＰＥは、まずあらかじめ決められた反復法による処理コードを取り出す（ステップＳ１）。
【００３４】
そして、そのコードに記述された計算がデータ格納法に依存するかどうかを判定し（ステップＳ２）、データ格納法に依存しなければ各ＰＥで同じ処理を実行する（ステップＳ３）。データ格納法に依存するかどうかは、例えば、あらかじめ処理コード内に書き込まれたパラメータ値により判断する。
【００３５】
次に、解が十分に収束するなどの反復終了条件が満たされたかどうかを判定し（ステップＳ５）、それが満たされなければステップＳ１以降の処理を繰り返す。
【００３６】
また、ステップＳ２において計算がデータ格納法に依存する場合は、他のＰＥと並列処理を行って、計算結果を中間配列に格納する（ステップＳ４）。この中間配列は、各種のデータ格納法に使用できるように標準的に用意された領域で、各ＰＥにより分割されている。例えば、ＰＥ１、ＰＥ２、ＰＥ３、ＰＥ４の４台のＰＥを備える場合は、中間配列は図５に示すようになる。
【００３７】
その後、各ＰＥはステップＳ５以降の処理を繰り返し、反復終了条件が満たされれば処理を終了する。
このような反復処理によれば、例えば行列ベクトル積に代表されるデータ格納法依存性の高い並列演算の結果を中間配列に格納して、反復法における他の計算処理をデータ格納法と独立に行うことができる。これにより、各種スパース行列の格納法を用いて、スパース行列の連立１次方程式を容易に解くことが可能になる。
【００３８】
次に、図６から図１２までを参照しながら、具体的なデータ格納法に基づく効率のよい行列ベクトル積の演算方法を説明する。
図６は、スパース行列の対角形式格納法を示している。ここでは、ｎ次のスパース行列Ａをバンド行列の１種とみなし、対角線方向に並んだ非零要素から成る対角ベクトルを２次元の格納配列４１に格納する。
【００３９】
Ａの対角ベクトル部分を拡大すると図７に示すようになり、対角ベクトル間には対角方向の要素がすべて０の部分が存在する。このような部分は格納配列４１に格納しなくてもよいので、対角方向に集中して非零要素が存在するスパース行列に特に適した格納法であると言える。
【００４０】
図６および図７に示すように、Ａをバンド行列とみなしたときの上バンド幅ｂａｎｄｗｕと下バンド幅ｂａｎｄｗｌは、それぞれ対角要素から対角ベクトルまでの距離の最大値として求められる。このとき、格納配列４１の幅ｗは最大ｂａｎｄｗｕ＋ｂａｎｄｗｌ＋１となる。
【００４１】
また、図６において対角ベクトルの長さを揃えるために、斜線部分に０の要素を付け加えて格納配列４１に格納しておく。これにより、Ａの対角ベクトルを用いた演算を均一に行うことができる。
【００４２】
配列４１の１次元目（行の次元）を各ＰＥに均等に分割して格納し、さらにｎ次元のベクトルＸ、Ｙを格納する配列を各ＰＥに分割配置することにより、次式の行列ベクトル積を並列に計算することができる。
Ｙ＝ＡＸ （１１）
例えば、４台のＰＥで並列演算を行う場合は、Ａの配列４１、Ｘの配列４２、Ｙの配列４３は、それぞれ図８に示すように分割される。このうち、演算結果を格納する配列４３は図５の中間配列に相当する。
【００４３】
ところで、図６のＡは一定のバンド領域にのみ非零要素を持つスパース行列であるため、（１１）式の演算において、各ＰＥは必ずしもＸのすべての要素を持っている必要はない。例えば図８のＰＥ２が受け持つ行列ベクトル積の有効部分は、図９に示すようになる。
【００４４】
図９においてＰＥ２に割り当てられたＡの行数をｎ２とすると、その非零要素が存在する領域はＡの斜線部分である。この斜線部分に該当する列の範囲は、ｎ２本の行と同じ番号を持つｎ２本の列と、その左側のｂａｎｄｗｌの幅の各列と、右側のｂａｎｄｗｕの幅の各列となる。したがって、Ｘの要素のうち、これらの列に対応する（ｎ２＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）個の要素のみが乗算の結果に寄与する。他の要素には０が乗算されるので、その結果もまた０となる。
【００４５】
Ｘを格納する配列４２は、図８に示すように各ＰＥに分散されているため、ＰＥ２は乗算を行うために、少なくともｂａｎｄｗｌ個の要素をＰＥ１からコピーし、少なくともｂａｎｄｗｕ個の要素をＰＥ３からコピーする必要がある。他のＰＥについても同様のコピー処理が必要になる。
【００４６】
そこで、あらかじめＰＥ毎にワーク用の一次元配列Ｗを用意しておき、この配列に必要なＸの要素をコピーすることにする。例えばＰＥ２の場合、Ｗの長さとしては、最低（ｎ２＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）だけあれば十分である。しかし、ここではＸの各要素の論理的な添え字を、そのままＷの対応する要素の添え字として使用できるように、各ＰＥのＷの長さを（ｎ＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）に統一することにする。
【００４７】
こうして、図１０に示すようなＷが各ＰＥに用意される。ここで、ＰＥの台数をＰ、ｉ番目のＰＥに割り当てられたＡの行数をｎｉ（ｉ＝１，．．．，Ｐ）とすると、ｉ番目のＰＥでは、Ｘの必要な要素が図１０に示すようにＷの斜線部分にコピーされる。ここで、ｎｉはおおよそｎ／Ｐとなる。Ｗの斜線部分以外の領域は実際の処理には用いられないので、どんな値を格納していてもかまわない。
【００４８】
Ｘの斜線部分の要素のうち、元々自分が持っているｎｉ個の要素はそのまま主記憶２１上でコピーすることができ、その上下のｂａｎｄｗｌ個とｂａｎｄｗｕ個の要素は、クロスバスイッチ２３を介して隣接ＰＥまたは他の近くのＰＥからコピーする。配列４１に格納された対角ベクトルの数がＡの次元ｎに比べて十分に小さい場合は、高々隣接ＰＥとの通信だけで必要な要素をコピーすることができる。
【００４９】
このように、図６のような対角方向に非零要素が集中しているスパース行列の場合は、特定のＰＥ間の通信だけで行列ベクトル積を実行でき、通信に伴うオーバヘッドを大幅に削減することができる。
【００５０】
また、１番目のＰＥは、Ｗの上端部の長さｂａｎｄｗｌの部分にダミーデータを格納して乗算を行い、Ｐ番目のＰＥは、下端部の長さｂａｎｄｗｕの部分にダミーデータを格納して乗算を行う。これらのダミーデータには、図６の斜線部分の要素である０が乗算されるので、ダミーデータの値が行列ベクトル積の結果に影響を与えることはない。しかし、ダミーデータを用いることで、１番目とＰ番目のＰＥは他のＰＥと同様の乗算を行うことができる。
【００５１】
このような乗算方法によれば、各ＰＥ内の演算はベクトル演算器を用いて効率よく処理することができる。また、各ＰＥは演算の途中でベクトル長を変更する必要もない。
【００５２】
図１１は、バンド幅を利用したスパース行列の行列ベクトル積演算のフローチャートである。図１１において処理が開始されると、各ＰＥは、まずスパース行列をバンド行列とみなして、その上バンド幅と下バンド幅を求める（ステップＳ１１）。次に、分割されたベクトルＸの要素のうち、自分が保持する部分（図１０のｎｉ個の要素）をワーク用配列Ｗに並列にコピーする（ステップＳ１２）。
【００５３】
次に、各ＰＥは、必要な下バンド幅の長さの部分（図１０のｂａｎｄｗｌ個の要素）をＷに並列にコピーする（ステップＳ１３）。また、必要な上バンド幅の長さの部分（図１０のｂａｎｄｗｕ個の要素）をＷに並列にコピーする（ステップＳ１４）。
【００５４】
そして、行列ベクトル積を並列に計算して、計算結果をベクトルＹの各ＰＥに割り当てられた部分に格納し（ステップＳ１５）、処理を終了する。
この方法によれば、ステップＳ１３およびＳ１４において限られたＰＥとのみ通信を行えばよいので、データ転送量が少なくて済み、行列ベクトル積演算が高速化される。
【００５５】
このようなバンド幅を利用した行列ベクトル積の計算方法は、対角形式格納法に限らず、他の任意のデータ格納法にも適用可能である。他のデータ格納法の１つとして、エルパック形式格納法がある。
【００５６】
エルパック形式格納法では、スパース行列の各行ベクトル毎に非零要素を圧縮した形で、行列要素を格納する。このとき、各行内の非零要素に対応する元のスパース行列の列番号は、格納配列ｉｃｏｆに格納されている。この列番号の値を利用すれば、バンド幅を求めることができる。
【００５７】
図１２は、配列ｉｃｏｆの例を示している。図１２において、左のスパース行列は１００００×１００００の行列で、対角要素と他の対角方向の線分上に非零要素を持っている。また、ｉｃｏｆは、スパース行列の第ｉ行に含まれる非零要素の列番号ｉｃｏｆ（ｉ，＊）を、各行番号ｉ毎に格納している。ここで、＊は第ｉ行内の１つ以上の非零要素に付けられた識別番号である。
【００５８】
スパース行列のすべての非零要素のうち、対角要素から最も離れたものとその行の対角要素との距離を求めれば、バンド幅が得られる。
例えば、ｉ＝９０００の行内には、５０００番目、６０００番目、および９０００番目の各列に非零要素があり、これらの列番号がｉｃｏｆ（９０００，＊）としてｉｃｏｆに格納されている。このうち、９０００番目の列に対角要素があり、それから左に最も離れた要素は５０００番目の列にある。したがって、これらの要素間の距離は−（５０００−９０００）＝４０００となる。
【００５９】
また、この５０００番目の列の要素はバンドの最も外側に位置しているため、この場合の下バンド幅ｂａｎｄｗｌは４０００となる。上バンド幅ｂａｎｄｗｕもｉｃｏｆから同様にして求められる。
【００６０】
ｉｃｏｆ（ｉ，＊）を用いたバンド幅の計算式を一般化すると、次式のようになる。
【００６１】
【数１】
【００６２】
（１２）、（１３）式において、＊に関するＭａｘは、第ｉ行内の＊番目の要素と対角要素との距離の最大値を表し、ｉに関するＭａｘは、各行の＊に関する最大値のうちの最大値を表す。
【００６３】
（１２）式で、（ｉｃｏｆ（ｉ，＊）−ｉ）の前に負の符号が付いているのは、左に位置する要素ほど対角要素までの距離を大きく評価するためである。この場合、対角要素より右に位置する要素までの距離は負の値を持つことになる。逆に、（１３）式では、右に位置する要素ほど対角要素までの距離が大きくなり、対角要素より左に位置する要素までの距離は負の値を持つ。
【００６４】
こうして、エルパック形式格納法においても上下バンド幅を計算することが可能になり、図１１の演算処理を適用することができる。
次に、図１３から図１７までを参照しながら、最も高速な反復解法の１つであるＭＧＣＲ法の実施形態について説明する。ＭＧＣＲ法は、各繰り返し段階における演算回数がＧＭＲＥＳ法より少ないため、この方法より高速であることが知られている（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）。
【００６５】
ＧＭＲＥＳ法は反復解法の中でもかなり速いアルゴリズムである。それより速いＭＧＣＲ法を実装することで、従来にない高速な反復処理装置を実現することができる。ここでは、一例としてスーパーコンピュータＶＰＰ５００を対象とした実装方法を説明する。
【００６６】
（１）式でｘ＝ｕ，ｂ＝ｆとおき直し、与えられた連立一次方程式を、
Ａｕ＝ｆ （１４）
と書くことにすると、ベクトルｕの近似解を求めるＭＧＣＲ法のアルゴリズムは次のようになる。
【００６７】
【数２】
【００６８】
ここで、（１５）式のｕ0 は任意に与えられた初期ベクトルであり、｜ｒ1 ｜はベクトルｒ1 の大きさを表す。また、ｓ1 ＝ｒ1 は１番目の探索ベクトルとなる。
【００６９】
（１６）〜（２３）式のループ処理では、ｉ＝１，．．．，ｋについてｋ個の探索ベクトルｓi+1 を生成し、（２３）式でｄi+1 が、与えられた収束判定値ε以内に収束すれば処理を終了する。そうでなければ、（２４）〜（２５）式によりｕ0 を更新して、（１５）式からの始まる外側のループ処理を繰り返す。ｋの値は任意に指定可能であるが、通常は１０〜１００の範囲に設定すればよい。
【００７０】
（１６）、（２１）式において、（ｖ１，ｖ２）のような表記はベクトルｖ１とｖ２の内積を表し、（１９）、（２３）、（２４）式において、ｓｏｌｖｅ Ｈｃ＝αのような表記は、連立方程式Ｈｃ＝αを解いてｃを求めることを意味する。
【００７１】
Ｈk は、次式で与えられるｋ次の上三角行列である。
【００７２】
【数３】
【００７３】
Ｈ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 も同様である。
【００７６】
（１７）式の 外１ は実際には計算されず、βk,k-1 とαk ′は次式のよう
【００７７】
【外１】
【００７８】
にして求められる。
【００７９】
【数４】
【００８０】
このＭＧＣＲ法のアルゴリズムを並列計算機上に実装するとき、（１５）式のＡ×ｕ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 の行列ベクトル積の計算処理を行う。
【００９２】
これらの行列ベクトル積の計算はＡの格納方法に依存するため、行列ベクトル積演算部５２に処理コードを与えて、入力ベクトルＸや出力ベクトルＹの格納場所等を指示する。このような制御用の処理コードは、コンパイラまたはプログラマにより作成される。
【００９３】
図１７のような構成を取ることにより、行列ベクトル積の演算をＭＧＣＲ法から分離して独立に実行することができ、各種のスパース行列格納法に対して、ＭＧＣＲ法を適用することが可能になる。
【００９４】
本発明で用いる反復解法はＭＧＣＲ法に限られるわけではなく、クリロフ部分空間を利用した任意の解法を用いることができる。また、図１７の構成において、反復部５１の処理をいずれかのＰＥに代表して行わせることも可能である。
【００９５】
【発明の効果】
本発明によれば、スパース行列の連立一次方程式を解くメモリ分散型並列計算機において、同一の反復解法のアルゴリズムを、多様なスパース行列の格納法に適用することが可能になる。したがって、反復解法の汎用性が向上し、プログラマの負担が軽減される。
【００９６】
また、スパース行列をバンド行列とみなして処理することで、行列ベクトル積の演算に必要なデータ転送量が削減され、並列処理の台数効果が向上する。特に、反復解法として他の解法より高速なＭＧＣＲ法を用いた場合、さらに効率のよい並列処理が実現される。
【図面の簡単な説明】
【図１】本発明の原理図である。
【図２】並列計算機の構成図である。
【図３】プロセッシング・エレメントの構成図である。
【図４】本発明の反復処理のフローチャートである。
【図５】中間配列を示す図である。
【図６】スパース行列の対角形式格納法を示す図である。
【図７】スパース行列の対角ベクトルの拡大図である。
【図８】行列ベクトル積の並列処理を示す図である。
【図９】ＰＥ２による行列ベクトル積を示す図である。
【図１０】ワーク用配列を示す図である。
【図１１】バンド幅を利用した行列ベクトル積演算のフローチャートである。
【図１２】エルパック形式格納法における列番号配列を示す図である。
【図１３】ＭＧＣＲ法における中間配列を示す図である。
【図１４】各ＰＥが保持する格納配列を示す図である。
【図１５】ＭＧＣＲ法による処理の概要を示す図である。
【図１６】ＭＧＣＲ法による処理のフローチャートである。
【図１７】ＭＧＣＲ法の装置構成図である。
【符号の説明】
１ 反復手段
２ 第１の並列処理手段
３ 中間配列記憶手段
１１ 入出力装置
１２ ＩＯバス
１３ ＰＥ
２１ 主記憶
２２ ＩＯポート
２３ クロスバスイッチ
２４ データ転送制御ユニット
２５ メモリ制御ユニット
２６ キャッシュメモリ
２７ ＧＰＲ／ＦＰＲ
２８ スカラ演算器
２９ マスクレジスタ
３０ ベクトルレジスタ
３１ ロードパイプライン
３２ ストアパイプライン
３３、３４ マスクパイプライン
３５ 乗算パイプライン
３６ 加算／論理演算パイプライン
３７ 除算パイプライン
３８ スカラユニット
３９ ベクトルユニット
４１、４２、４３、Ｗ、Ｗ１、Ｗ２、ｉｃｏｆ 格納配列
５１ 反復部
５２ 行列ベクトル積演算部
【発明の属する技術分野】
本発明は、並列計算機を利用した連立一次方程式の反復解法に係り、係数行列を複数のプロセッシング・エレメントに分散配置して、並列処理により解を求める並列処理装置およびその方法に関する。
【０００２】
【従来の技術とその問題点】
物理現象の解析において現れる偏微分方程式の境界値問題や行列の固有値問題を解く時、一般に、大規模なスパース行列（疎行列）を係数行列とする次のような連立一次方程式を解く必要が生じる。
Ａｘ＝ｂ （１）
ここで、Ａは一般にｎ×ｎの非対称行列、ｘはｎ次元の変数ベクトル、ｂはｎ次元の定数ベクトルである。ｎの値が１００００以上になることも珍しくない。
【０００３】
大規模な連立一次方程式は、気象予測、原子炉設計、半導体の回路解析、航空工学における流体解析、構造物の構造解析等の多くの科学技術計算に用いられる。また、大規模な固有値問題は、構造物の構造解析、回路解析、地球科学における地震予知、原子炉の安全性解析、分子科学における多電子系のエネルギー計算、原子核の構造解析等の分野において、物理系の固有振動を記述するときに現れる。
【０００４】
したがって、（１）式のような大規模な連立一次方程式を効率よく高速に解くことは、科学技術計算の重要な問題の１つである。今日では、計算を高速化するために、複数のプロセッシング・エレメント（ＰＥ）を備えたメモリ分散型の並列計算機が多く用いられている。
【０００５】
計算機を用いて（１）式を解く１つの方法として、ＡをＬＵ分解するガウスの消去法に基づいた直接法がある。しかし、Ａが大きなスパース行列の場合、非零要素が各行に数個しかないこともあり、計算コストや記憶領域の面で無駄が多い。そこで、単純な行列ベクトル積を繰り返して近似解を求める反復解法が多く用いられている。
【０００６】
反復解法の多くはクリロフ（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 ）形式格納法を用いて、行列ベクトル積を並列に計算することができる。このとき、スパース行列を疑似的なバンド行列とみなしてそのバンド幅を求め、これを利用して行列ベクトル積を計算すれば、各ＰＥは近くの限られたＰＥとのみデータ転送を行えばよくなる。したがって、転送回数が削減され、色々なバンド幅のスパース行列について、行列ベクトル積の演算が効率よく実行される。
【００２４】
さらに、一般スパース行列の連立１次方程式の反復解法の１つであるＭＧＣＲ法を例に取り、本発明の反復処理をメモリ分散型のスーパーコンピュータＶＰＰシリーズ向けに構成した。この形態によれば、各種のデータ格納方法に対してＭＧＣＲ法を適用し、高速な計算処理を実現することが可能になる。
【００２５】
まず、図２および図３を参照しながら、実施形態で用いる並列計算機の構成を説明する。図２は、並列計算機の概略構成図である。図２の並列計算機は、入出力装置１１と、ＩＯバス１２と、クロスバ方式により互いに接続された複数のＰＥ１３とを備える。
【００２６】
入出力装置１１は、各ＰＥ１３との間でデータの入出力を行うための装置で、例えば、キーボードなどの入力機器とディスプレイやプリンタなどの出力機器を備えた計算機端末である。入出力装置１１は、ＩＯバス１２により各ＰＥ１３と接続されている。
【００２７】
図３は、ＰＥ１３の構成図である。図５のＰＥ１３は、主記憶２１、ＩＯポート２２、クロスバスイッチ２３、データ転送制御ユニット２４、メモリ制御ユニット２５、スカラユニット３８、およびベクトルユニット３９を備える。
【００２８】
主記憶２１は、与えられた問題の行列要素や計算の中間結果などを格納する。ＩＯポート２２はＩＯバス１２に接続され、入出力装置１１との間で入出力データの転送を行う。クロスバスイッチ２３は、主記憶２１内のデータを他のＰＥ１３に転送する際のデータの切り換え等を行う。データ転送制御ユニット２４は、クロスバスイッチ２３を介して他のＰＥ１３とのデータ転送および同期制御を行う。
【００２９】
メモリ制御ユニット２５は、データ転送制御ユニット２４、スカラユニット３８、およびベクトルユニット３９から発生するアクセス要求を受け取り、主記憶２１に対するアクセスを制御する。
【００３０】
スカラユニット３８は、キャッシュメモリ２６、ＧＰＲ（汎用レジスタ）／ＦＰＲ（浮動小数点レジスタ）２７、スカラ演算器２８を備え、スカラ演算を実行する。
【００３１】
ベクトルユニット３９は、マスクレジスタ２９とベクトルレジスタ３０の他に、ロードパイプライン３１、ストアパイプライン３２、マスクパイプライン３３、３４、乗算パイプライン３５、加算／論理演算パイプライン３６、および除算パイプライン３７を備える。
【００３２】
これらの各パイプラインは、それぞれ複数要素を同時に処理することができ、パイプライン３５、３６、３７のうち２本は同時に動作できる。また、他のパイプラインについては、すべて同時に動作可能である。マスクパイプライン３３は総和／検索処理用に使用され、マスクパイプライン３４は論理演算処理用に使用される。
【００３３】
ベクトルユニット３９は、これらのパイプラインを用いて複数のベクトル命令を並列に実行することができ、高速な行列演算が可能である。
次に、図４および図５を参照しながら、本発明の反復処理について説明する。図４は、多様なデータ格納方法に対応可能な反復処理のフローチャートである。図４において処理が開始されると、並列計算機の各ＰＥは、まずあらかじめ決められた反復法による処理コードを取り出す（ステップＳ１）。
【００３４】
そして、そのコードに記述された計算がデータ格納法に依存するかどうかを判定し（ステップＳ２）、データ格納法に依存しなければ各ＰＥで同じ処理を実行する（ステップＳ３）。データ格納法に依存するかどうかは、例えば、あらかじめ処理コード内に書き込まれたパラメータ値により判断する。
【００３５】
次に、解が十分に収束するなどの反復終了条件が満たされたかどうかを判定し（ステップＳ５）、それが満たされなければステップＳ１以降の処理を繰り返す。
【００３６】
また、ステップＳ２において計算がデータ格納法に依存する場合は、他のＰＥと並列処理を行って、計算結果を中間配列に格納する（ステップＳ４）。この中間配列は、各種のデータ格納法に使用できるように標準的に用意された領域で、各ＰＥにより分割されている。例えば、ＰＥ１、ＰＥ２、ＰＥ３、ＰＥ４の４台のＰＥを備える場合は、中間配列は図５に示すようになる。
【００３７】
その後、各ＰＥはステップＳ５以降の処理を繰り返し、反復終了条件が満たされれば処理を終了する。
このような反復処理によれば、例えば行列ベクトル積に代表されるデータ格納法依存性の高い並列演算の結果を中間配列に格納して、反復法における他の計算処理をデータ格納法と独立に行うことができる。これにより、各種スパース行列の格納法を用いて、スパース行列の連立１次方程式を容易に解くことが可能になる。
【００３８】
次に、図６から図１２までを参照しながら、具体的なデータ格納法に基づく効率のよい行列ベクトル積の演算方法を説明する。
図６は、スパース行列の対角形式格納法を示している。ここでは、ｎ次のスパース行列Ａをバンド行列の１種とみなし、対角線方向に並んだ非零要素から成る対角ベクトルを２次元の格納配列４１に格納する。
【００３９】
Ａの対角ベクトル部分を拡大すると図７に示すようになり、対角ベクトル間には対角方向の要素がすべて０の部分が存在する。このような部分は格納配列４１に格納しなくてもよいので、対角方向に集中して非零要素が存在するスパース行列に特に適した格納法であると言える。
【００４０】
図６および図７に示すように、Ａをバンド行列とみなしたときの上バンド幅ｂａｎｄｗｕと下バンド幅ｂａｎｄｗｌは、それぞれ対角要素から対角ベクトルまでの距離の最大値として求められる。このとき、格納配列４１の幅ｗは最大ｂａｎｄｗｕ＋ｂａｎｄｗｌ＋１となる。
【００４１】
また、図６において対角ベクトルの長さを揃えるために、斜線部分に０の要素を付け加えて格納配列４１に格納しておく。これにより、Ａの対角ベクトルを用いた演算を均一に行うことができる。
【００４２】
配列４１の１次元目（行の次元）を各ＰＥに均等に分割して格納し、さらにｎ次元のベクトルＸ、Ｙを格納する配列を各ＰＥに分割配置することにより、次式の行列ベクトル積を並列に計算することができる。
Ｙ＝ＡＸ （１１）
例えば、４台のＰＥで並列演算を行う場合は、Ａの配列４１、Ｘの配列４２、Ｙの配列４３は、それぞれ図８に示すように分割される。このうち、演算結果を格納する配列４３は図５の中間配列に相当する。
【００４３】
ところで、図６のＡは一定のバンド領域にのみ非零要素を持つスパース行列であるため、（１１）式の演算において、各ＰＥは必ずしもＸのすべての要素を持っている必要はない。例えば図８のＰＥ２が受け持つ行列ベクトル積の有効部分は、図９に示すようになる。
【００４４】
図９においてＰＥ２に割り当てられたＡの行数をｎ２とすると、その非零要素が存在する領域はＡの斜線部分である。この斜線部分に該当する列の範囲は、ｎ２本の行と同じ番号を持つｎ２本の列と、その左側のｂａｎｄｗｌの幅の各列と、右側のｂａｎｄｗｕの幅の各列となる。したがって、Ｘの要素のうち、これらの列に対応する（ｎ２＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）個の要素のみが乗算の結果に寄与する。他の要素には０が乗算されるので、その結果もまた０となる。
【００４５】
Ｘを格納する配列４２は、図８に示すように各ＰＥに分散されているため、ＰＥ２は乗算を行うために、少なくともｂａｎｄｗｌ個の要素をＰＥ１からコピーし、少なくともｂａｎｄｗｕ個の要素をＰＥ３からコピーする必要がある。他のＰＥについても同様のコピー処理が必要になる。
【００４６】
そこで、あらかじめＰＥ毎にワーク用の一次元配列Ｗを用意しておき、この配列に必要なＸの要素をコピーすることにする。例えばＰＥ２の場合、Ｗの長さとしては、最低（ｎ２＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）だけあれば十分である。しかし、ここではＸの各要素の論理的な添え字を、そのままＷの対応する要素の添え字として使用できるように、各ＰＥのＷの長さを（ｎ＋ｂａｎｄｗｌ＋ｂａｎｄｗｕ）に統一することにする。
【００４７】
こうして、図１０に示すようなＷが各ＰＥに用意される。ここで、ＰＥの台数をＰ、ｉ番目のＰＥに割り当てられたＡの行数をｎｉ（ｉ＝１，．．．，Ｐ）とすると、ｉ番目のＰＥでは、Ｘの必要な要素が図１０に示すようにＷの斜線部分にコピーされる。ここで、ｎｉはおおよそｎ／Ｐとなる。Ｗの斜線部分以外の領域は実際の処理には用いられないので、どんな値を格納していてもかまわない。
【００４８】
Ｘの斜線部分の要素のうち、元々自分が持っているｎｉ個の要素はそのまま主記憶２１上でコピーすることができ、その上下のｂａｎｄｗｌ個とｂａｎｄｗｕ個の要素は、クロスバスイッチ２３を介して隣接ＰＥまたは他の近くのＰＥからコピーする。配列４１に格納された対角ベクトルの数がＡの次元ｎに比べて十分に小さい場合は、高々隣接ＰＥとの通信だけで必要な要素をコピーすることができる。
【００４９】
このように、図６のような対角方向に非零要素が集中しているスパース行列の場合は、特定のＰＥ間の通信だけで行列ベクトル積を実行でき、通信に伴うオーバヘッドを大幅に削減することができる。
【００５０】
また、１番目のＰＥは、Ｗの上端部の長さｂａｎｄｗｌの部分にダミーデータを格納して乗算を行い、Ｐ番目のＰＥは、下端部の長さｂａｎｄｗｕの部分にダミーデータを格納して乗算を行う。これらのダミーデータには、図６の斜線部分の要素である０が乗算されるので、ダミーデータの値が行列ベクトル積の結果に影響を与えることはない。しかし、ダミーデータを用いることで、１番目とＰ番目のＰＥは他のＰＥと同様の乗算を行うことができる。
【００５１】
このような乗算方法によれば、各ＰＥ内の演算はベクトル演算器を用いて効率よく処理することができる。また、各ＰＥは演算の途中でベクトル長を変更する必要もない。
【００５２】
図１１は、バンド幅を利用したスパース行列の行列ベクトル積演算のフローチャートである。図１１において処理が開始されると、各ＰＥは、まずスパース行列をバンド行列とみなして、その上バンド幅と下バンド幅を求める（ステップＳ１１）。次に、分割されたベクトルＸの要素のうち、自分が保持する部分（図１０のｎｉ個の要素）をワーク用配列Ｗに並列にコピーする（ステップＳ１２）。
【００５３】
次に、各ＰＥは、必要な下バンド幅の長さの部分（図１０のｂａｎｄｗｌ個の要素）をＷに並列にコピーする（ステップＳ１３）。また、必要な上バンド幅の長さの部分（図１０のｂａｎｄｗｕ個の要素）をＷに並列にコピーする（ステップＳ１４）。
【００５４】
そして、行列ベクトル積を並列に計算して、計算結果をベクトルＹの各ＰＥに割り当てられた部分に格納し（ステップＳ１５）、処理を終了する。
この方法によれば、ステップＳ１３およびＳ１４において限られたＰＥとのみ通信を行えばよいので、データ転送量が少なくて済み、行列ベクトル積演算が高速化される。
【００５５】
このようなバンド幅を利用した行列ベクトル積の計算方法は、対角形式格納法に限らず、他の任意のデータ格納法にも適用可能である。他のデータ格納法の１つとして、エルパック形式格納法がある。
【００５６】
エルパック形式格納法では、スパース行列の各行ベクトル毎に非零要素を圧縮した形で、行列要素を格納する。このとき、各行内の非零要素に対応する元のスパース行列の列番号は、格納配列ｉｃｏｆに格納されている。この列番号の値を利用すれば、バンド幅を求めることができる。
【００５７】
図１２は、配列ｉｃｏｆの例を示している。図１２において、左のスパース行列は１００００×１００００の行列で、対角要素と他の対角方向の線分上に非零要素を持っている。また、ｉｃｏｆは、スパース行列の第ｉ行に含まれる非零要素の列番号ｉｃｏｆ（ｉ，＊）を、各行番号ｉ毎に格納している。ここで、＊は第ｉ行内の１つ以上の非零要素に付けられた識別番号である。
【００５８】
スパース行列のすべての非零要素のうち、対角要素から最も離れたものとその行の対角要素との距離を求めれば、バンド幅が得られる。
例えば、ｉ＝９０００の行内には、５０００番目、６０００番目、および９０００番目の各列に非零要素があり、これらの列番号がｉｃｏｆ（９０００，＊）としてｉｃｏｆに格納されている。このうち、９０００番目の列に対角要素があり、それから左に最も離れた要素は５０００番目の列にある。したがって、これらの要素間の距離は−（５０００−９０００）＝４０００となる。
【００５９】
また、この５０００番目の列の要素はバンドの最も外側に位置しているため、この場合の下バンド幅ｂａｎｄｗｌは４０００となる。上バンド幅ｂａｎｄｗｕもｉｃｏｆから同様にして求められる。
【００６０】
ｉｃｏｆ（ｉ，＊）を用いたバンド幅の計算式を一般化すると、次式のようになる。
【００６１】
【数１】
（１２）、（１３）式において、＊に関するＭａｘは、第ｉ行内の＊番目の要素と対角要素との距離の最大値を表し、ｉに関するＭａｘは、各行の＊に関する最大値のうちの最大値を表す。
【００６３】
（１２）式で、（ｉｃｏｆ（ｉ，＊）−ｉ）の前に負の符号が付いているのは、左に位置する要素ほど対角要素までの距離を大きく評価するためである。この場合、対角要素より右に位置する要素までの距離は負の値を持つことになる。逆に、（１３）式では、右に位置する要素ほど対角要素までの距離が大きくなり、対角要素より左に位置する要素までの距離は負の値を持つ。
【００６４】
こうして、エルパック形式格納法においても上下バンド幅を計算することが可能になり、図１１の演算処理を適用することができる。
次に、図１３から図１７までを参照しながら、最も高速な反復解法の１つであるＭＧＣＲ法の実施形態について説明する。ＭＧＣＲ法は、各繰り返し段階における演算回数がＧＭＲＥＳ法より少ないため、この方法より高速であることが知られている（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）。
【００６５】
ＧＭＲＥＳ法は反復解法の中でもかなり速いアルゴリズムである。それより速いＭＧＣＲ法を実装することで、従来にない高速な反復処理装置を実現することができる。ここでは、一例としてスーパーコンピュータＶＰＰ５００を対象とした実装方法を説明する。
【００６６】
（１）式でｘ＝ｕ，ｂ＝ｆとおき直し、与えられた連立一次方程式を、
Ａｕ＝ｆ （１４）
と書くことにすると、ベクトルｕの近似解を求めるＭＧＣＲ法のアルゴリズムは次のようになる。
【００６７】
【数２】
ここで、（１５）式のｕ0 は任意に与えられた初期ベクトルであり、｜ｒ1 ｜はベクトルｒ1 の大きさを表す。また、ｓ1 ＝ｒ1 は１番目の探索ベクトルとなる。
【００６９】
（１６）〜（２３）式のループ処理では、ｉ＝１，．．．，ｋについてｋ個の探索ベクトルｓi+1 を生成し、（２３）式でｄi+1 が、与えられた収束判定値ε以内に収束すれば処理を終了する。そうでなければ、（２４）〜（２５）式によりｕ0 を更新して、（１５）式からの始まる外側のループ処理を繰り返す。ｋの値は任意に指定可能であるが、通常は１０〜１００の範囲に設定すればよい。
【００７０】
（１６）、（２１）式において、（ｖ１，ｖ２）のような表記はベクトルｖ１とｖ２の内積を表し、（１９）、（２３）、（２４）式において、ｓｏｌｖｅ Ｈｃ＝αのような表記は、連立方程式Ｈｃ＝αを解いてｃを求めることを意味する。
【００７１】
Ｈk は、次式で与えられるｋ次の上三角行列である。
【００７２】
【数３】
Ｈ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 も同様である。
【００７６】
（１７）式の 外１ は実際には計算されず、βk,k-1 とαk ′は次式のよう
【００７７】
【外１】
にして求められる。
【００７９】
【数４】
このＭＧＣＲ法のアルゴリズムを並列計算機上に実装するとき、（１５）式のＡ×ｕ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 の行列ベクトル積の計算処理を行う。
【００９２】
これらの行列ベクトル積の計算はＡの格納方法に依存するため、行列ベクトル積演算部５２に処理コードを与えて、入力ベクトルＸや出力ベクトルＹの格納場所等を指示する。このような制御用の処理コードは、コンパイラまたはプログラマにより作成される。
【００９３】
図１７のような構成を取ることにより、行列ベクトル積の演算をＭＧＣＲ法から分離して独立に実行することができ、各種のスパース行列格納法に対して、ＭＧＣＲ法を適用することが可能になる。
【００９４】
本発明で用いる反復解法はＭＧＣＲ法に限られるわけではなく、クリロフ部分空間を利用した任意の解法を用いることができる。また、図１７の構成において、反復部５１の処理をいずれかのＰＥに代表して行わせることも可能である。
【００９５】
【発明の効果】
本発明によれば、スパース行列の連立一次方程式を解くメモリ分散型並列計算機において、同一の反復解法のアルゴリズムを、多様なスパース行列の格納法に適用することが可能になる。したがって、反復解法の汎用性が向上し、プログラマの負担が軽減される。
【００９６】
また、スパース行列をバンド行列とみなして処理することで、行列ベクトル積の演算に必要なデータ転送量が削減され、並列処理の台数効果が向上する。特に、反復解法として他の解法より高速なＭＧＣＲ法を用いた場合、さらに効率のよい並列処理が実現される。
【図面の簡単な説明】
【図１】本発明の原理図である。
【図２】並列計算機の構成図である。
【図３】プロセッシング・エレメントの構成図である。
【図４】本発明の反復処理のフローチャートである。
【図５】中間配列を示す図である。
【図６】スパース行列の対角形式格納法を示す図である。
【図７】スパース行列の対角ベクトルの拡大図である。
【図８】行列ベクトル積の並列処理を示す図である。
【図９】ＰＥ２による行列ベクトル積を示す図である。
【図１０】ワーク用配列を示す図である。
【図１１】バンド幅を利用した行列ベクトル積演算のフローチャートである。
【図１２】エルパック形式格納法における列番号配列を示す図である。
【図１３】ＭＧＣＲ法における中間配列を示す図である。
【図１４】各ＰＥが保持する格納配列を示す図である。
【図１５】ＭＧＣＲ法による処理の概要を示す図である。
【図１６】ＭＧＣＲ法による処理のフローチャートである。
【図１７】ＭＧＣＲ法の装置構成図である。
【符号の説明】
１ 反復手段
２ 第１の並列処理手段
３ 中間配列記憶手段
１１ 入出力装置
１２ ＩＯバス
１３ ＰＥ
２１ 主記憶
２２ ＩＯポート
２３ クロスバスイッチ
２４ データ転送制御ユニット
２５ メモリ制御ユニット
２６ キャッシュメモリ
２７ ＧＰＲ／ＦＰＲ
２８ スカラ演算器
２９ マスクレジスタ
３０ ベクトルレジスタ
３１ ロードパイプライン
３２ ストアパイプライン
３３、３４ マスクパイプライン
３５ 乗算パイプライン
３６ 加算／論理演算パイプライン
３７ 除算パイプライン
３８ スカラユニット
３９ ベクトルユニット
４１、４２、４３、Ｗ、Ｗ１、Ｗ２、ｉｃｏｆ 格納配列
５１ 反復部
５２ 行列ベクトル積演算部
Claims (12)
- スパース行列を係数行列とする連立方程式を、反復解法により解く並列計算機において、
前記係数行列を複数部分に分割して格納し、該係数行列の格納形式に依存する計算を並列に行う並列演算手段と、
前記並列演算手段による計算結果を中間配列として、該中間配列の第１次元を均等に分割して記憶する中間配列記憶手段と、
前記中間配列のデータを用いて、前記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する反復手段と、を備え
前記反復手段は、与えられた初期ベクトルを用いて行列ベクトル積を利用した前記反復処理を開始し、前記並列演算手段は、処理コードを介して、処理コード内のパラメータ値に基づき、前記係数行列の格納形式に依存する計算を識別し、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算して中間ベクトルを生成し、前記中間配列記憶手段は、該中間ベクトルを前記中間配列内に記憶し、前記反復手段は、前記反復処理が終了した後、得られた解ベクトルを出力し、前記係数行列の格納形式に依存する計算は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算することであることを特徴とするの並列処理装置。 - スパース行列を係数行列とする連立方程式を、反復解法により解く並列計算機において、
前記係数行列を複数部分に分割して格納し、該係数行列の格納形式に依存する計算を並列に行う並列演算手段と、
前記並列演算手段による計算結果を中間配列として、該中間配列の第１次元を均等に分割して記憶する中間配列記憶手段と、
前記中間配列のデータを用いて、前記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する反復手段と、を備え
前記反復手段は、与えられた初期ベクトルを用いて行列ベクトル積を利用した前記反復処理を開始し、前記並列演算手段は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算して中間ベクトルを生成し、前記中間配列記憶手段は、該中間ベクトルを前記中間配列内に記憶し、前記反復手段は、前記反復処理が終了した後、得られた解ベクトルを出力し、前記並列演算手段は、前記係数行列を分割して保持する複数のプロセッシング・エレメント手段を備え、各プロセッシング・エレメント手段は、前記係数行列の割り当てられた部分を用いて前記係数行列の格納形式に依存する計算を行い、前記係数行列の格納形式に依存する計算は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算することであることを特徴とする並列処理装置。 - スパース行列を係数行列とする連立方程式を、反復解法により解く並列計算機において、
前記係数行列を複数部分に分割して格納し、該係数行列の格納形式に依存する計算を並列に行う並列演算手段と、
前記並列演算手段による計算結果を中間配列として、該中間配列の第１次元を均等に分割して記憶する中間配列記憶手段と、
前記中間配列のデータを用いて、前記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する反復手段と、を備え
前記反復手段は、与えられた初期ベクトルを用いて行列ベクトル積を利用した前記反復処理を開始し、前記並列演算手段は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算して中間ベクトルを生成し、前記中間配列記憶手段は、該中間ベクトルを前記中間配列内に記憶し、前記反復手段は、前記反復処理が終了した後、得られた解ベクトルを出力し、前記並列演算手段は、前記係数行列の第１の部分に乗算する要素を格納する第１のワーク用配列を有し、前記係数行列のバンド幅に基づいて前記中間配列の必要な要素を決定し、該必要な要素を第１のワーク用配列に集めて、該第１の部分と第１のワーク用配列を用いた乗算を行い、前記係数行列の格納形式に依存する計算は、前記係 数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算することであることを特徴とする並列処理装置。 - スパース行列を係数行列とする連立方程式を、反復解法により解く並列計算機において、
前記係数行列を複数部分に分割して格納し、該係数行列の格納形式に依存する計算を並列に行う並列演算手段と、
前記並列演算手段による計算結果を中間配列として、該中間配列の第１次元を均等に分割して記憶する中間配列記憶手段と、
前記中間配列のデータを用いて、前記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する反復手段と、を備え
前記反復手段は、複数のプロセッシング・エレメント手段を備え、各プロセッシング・エレメント手段は、前記係数行列の格納形式に依存しない処理を冗長に実行し、前記並列演算手段は、処理コードを介して、処理コード内のパラメータ値に基づき、前記係数行列の格納形式に依存する計算を識別し、前記係数行列の格納形式に依存する計算は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算することであり、前記係数行列の格納形式に依存しない処理は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算すること以外の処理であることを特徴とする並列処理装置。 - スパース行列を係数行列とする連立方程式を、反復解法により解く並列計算機において、
前記係数行列を複数部分に分割して格納し、該係数行列の格納形式に依存する計算を並列に行う並列演算手段と、
前記並列演算手段による計算結果を中間配列として、該中間配列の第１次元を均等に分割して記憶する中間配列記憶手段と、
前記中間配列のデータを用いて、前記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する反復手段と、を備え
前記反復手段は、複数のプロセッシング・エレメント手段を備え、各プロセッシング・エレメント手段は、前記係数行列の格納形式に依存しない処理を冗長に実行し、前記並列演算手段は、前記係数行列を分割して保持する複数のプロセッシング・エレメント手段を備え、各プロセッシング・エレメント手段は、前記係数行列の割り当てられた部分を用いて前記係数行列の格納形式に依存する計算を行い、前記係数行列の格納形式に依存する計算は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算することであり、前記係数行列の格納形式に依存しない処理は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算すること以外の処理であることを特徴とする並列処理装置。 - スパース行列を係数行列とする連立方程式を、反復解法により解く並列計算機において、
前記係数行列を複数部分に分割して格納し、該係数行列の格納形式に依存する計算を並列に行う並列演算手段と、
前記並列演算手段による計算結果を中間配列として、該中間配列の第１次元を均等に分割して記憶する中間配列記憶手段と、
前記中間配列のデータを用いて、前記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する反復手段と、を備え
前記反復手段は、複数のプロセッシング・エレメント手段を備え、各プロセッシング・エレメント手段は、前記係数行列の格納形式に依存しない処理を冗長に実行し、前記並列演算手段は、前記係数行列の第１の部分に乗算する要素を格納する第１のワーク用配列を有し、前記係数行列のバンド幅に基づいて前記中間配列の必要な要素を決定し、該必要な要素を第１のワーク用配列に集めて、該第１の部分と第１のワーク用配列を用いた乗算を行い、前記係数行列の格納形式に依存する計算は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算することであり、前記係数行列の格納形式に依存しない処理は、前記係数行列の格納形式に応じて該係数行列とベクトルの積を並列に計算する こと以外の処理であることを特徴とする並列処理装置。 - 前記反復手段は、１つ前の繰り返し段階で生成された中間ベクトルを用いて、前記行列ベクトル積以外の計算処理または判定処理を行うことを特徴とする請求項１乃至３記載の並列処理装置。
- 前記並列演算手段は、前記係数行列を対角形式格納法またはエルパック形式格納法で分割して格納する複数のプロセッシング・エレメント手段を備え、プロセッシング・エレメント手段同士の通信により、前記必要な要素を第１のワーク用配列に集めることを特徴とする請求項３または６記載の並列処理装置。
- 前記中間配列記憶手段は、前記中間配列を前記複数のプロセッシング・エレメント手段に分割して記憶し、第１のプロセッシング・エレメント手段は、前記第１の部分および第１のワーク用配列を格納し、第２のプロセッシング・エレメント手段から転送された要素を該第１のワーク用配列に格納して、前記乗算を行うことを特徴とする請求項８記載の並列処理装置。
- 前記並列演算手段は、前記係数行列の上バンド幅または下バンド幅に基づいて前記必要な要素を決定することを特徴とする請求項３または６記載の並列処理装置。
- クリロフ部分空間を利用した反復解法を用いて前記連立方程式を解くことを特徴とする請求項１乃至６記載の並列処理装置。
- ＭＧＣＲ法を用いて前記連立方程式を解くことを特徴とする請求項１乃至６記載の並列処理装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP01610596A JP3808925B2 (ja) | 1996-01-31 | 1996-01-31 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP01610596A JP3808925B2 (ja) | 1996-01-31 | 1996-01-31 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Related Child Applications (2)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2006031724A Division JP2006164306A (ja) | 2006-02-08 | 2006-02-08 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
| JP2006031725A Division JP2006164307A (ja) | 2006-02-08 | 2006-02-08 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH09212483A JPH09212483A (ja) | 1997-08-15 |
| JP3808925B2 true JP3808925B2 (ja) | 2006-08-16 |
Family
ID=11907246
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP01610596A Expired - Fee Related JP3808925B2 (ja) | 1996-01-31 | 1996-01-31 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3808925B2 (ja) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE60037235T2 (de) * | 1999-06-14 | 2008-10-16 | Panalytical B.V. | Verfahren zur ermittlung eines wahren spektrums aus einem gemessenen spektrum |
| JP5262177B2 (ja) * | 2008-02-22 | 2013-08-14 | 富士通株式会社 | ベクトル積の並列処理方法 |
| JP5767576B2 (ja) | 2011-12-19 | 2015-08-19 | インターナショナル・ビジネス・マシーンズ・コーポレーションＩｎｔｅｒｎａｔｉｏｎａｌ Ｂｕｓｉｎｅｓｓ Ｍａｃｈｉｎｅｓ Ｃｏｒｐｏｒａｔｉｏｎ | システム同定のための行列格納方法、プログラム及びシステム |
-
1996
- 1996-01-31 JP JP01610596A patent/JP3808925B2/ja not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JPH09212483A (ja) | 1997-08-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US5887186A (en) | Method of solving simultaneous linear equations in a memory-distributed parallel computer | |
| US7337205B2 (en) | Matrix multiplication in a vector processing system | |
| US20030088600A1 (en) | Matrix transposition in a computer system | |
| EP0523544A2 (en) | Calculating equipment for solving systems of linear equation | |
| CN113326916A (zh) | 将卷积映射到分区通道卷积引擎 | |
| CN113344172A (zh) | 将卷积映射到通道卷积引擎 | |
| US20200226201A1 (en) | Methods and Apparatus for Constructing Digital Circuits for Performing Matrix Operations | |
| JP3808925B2 (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
| Geist et al. | Finding eigenvalues and eigenvectors of unsymmetric matrices using a distributed-memory multiprocessor | |
| CN113435569A (zh) | 使用每通道卷积运算的流水线逐点卷积 | |
| JP2006164307A (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
| Khan et al. | Parallel finite element analysis using Jacobi-conditioned conjugate gradient algorithm | |
| Aykanat et al. | Vectorization and parallelization of the conjugate gradient algorithm on hypercube-connected vector processors | |
| JP3659307B2 (ja) | ベクトル・コンピュータにおける演算方法、及び記録媒体 | |
| JP2006164306A (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
| CN111522776B (zh) | 一种计算架构 | |
| JP2018206078A (ja) | 並列処理装置、並列演算方法、及び並列演算プログラム | |
| Shewchuk et al. | A compiler for parallel finite element methods with domain-decomposed unstructured meshes | |
| BOSTIC et al. | A Lanczos eigenvalue method on a parallel computer | |
| Bischof et al. | A Cholesky up-and downdating algorithm for systolic and SIMD architectures | |
| JP2806262B2 (ja) | マルチプロセッサシステムのプロセス割当方法 | |
| Zubair et al. | Optimization of a solver for computational materials and structures problems on NVIDIA Volta and AMD Instinct GPUs | |
| Fatoohi et al. | Implementation of an ADI method on parallel computers | |
| Bischof et al. | A parallel algorithm for the multi-input Sylvester-observer equation | |
| JP3542184B2 (ja) | 線形計算方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20051101 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20051108 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060208 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20060516 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20060519 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090526 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100526 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100526 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110526 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120526 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130526 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140526 Year of fee payment: 8 |
|
| LAPS | Cancellation because of no payment of annual fees |