JP2006164307A - 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 - Google Patents
多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 Download PDFInfo
- Publication number
- JP2006164307A JP2006164307A JP2006031725A JP2006031725A JP2006164307A JP 2006164307 A JP2006164307 A JP 2006164307A JP 2006031725 A JP2006031725 A JP 2006031725A JP 2006031725 A JP2006031725 A JP 2006031725A JP 2006164307 A JP2006164307 A JP 2006164307A
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- vector
- parallel
- iterative
- processing
- 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.)
- Pending
Links
Images
Landscapes
- Complex Calculations (AREA)
Abstract
【課題】反復解法により連立一次方程式を解くメモリ分散型並列計算機において、多様なデータ格納方法に対応して効率的な並列処理行うことが課題である。
【解決手段】反復解法の1つであるMGCR法の探索ベクトルsiを格納する配列W1は4台のPEに分散配置され、係数行列Aとsiを用いて計算された三角行列を格納する配列W2は各PEに配置される。A×siのような行列ベクトル積は4台のPEにより並列に計算され、Aの格納方法に依存しない計算は各PE内で冗長に処理される。行列ベクトル積を反復解法のアルゴリズムから独立させることにより、各種の行列格納方法に対応することが可能になり、処理の汎用性が向上する。また、MGCR法を実装することで高速な反復解法が実現される。
【選択図】図15
【解決手段】反復解法の1つであるMGCR法の探索ベクトルsiを格納する配列W1は4台のPEに分散配置され、係数行列Aとsiを用いて計算された三角行列を格納する配列W2は各PEに配置される。A×siのような行列ベクトル積は4台のPEにより並列に計算され、Aの格納方法に依存しない計算は各PE内で冗長に処理される。行列ベクトル積を反復解法のアルゴリズムから独立させることにより、各種の行列格納方法に対応することが可能になり、処理の汎用性が向上する。また、MGCR法を実装することで高速な反復解法が実現される。
【選択図】図15
Description
本発明は、並列計算機を利用した連立一次方程式の反復解法に係り、係数行列を複数のプロセッシング・エレメントに分散配置して、並列処理により解を求める並列処理装置およびその方法に関する。
物理現象の解析において現れる偏微分方程式の境界値問題や行列の固有値問題を解く時、一般に、大規模なスパース行列(疎行列)を係数行列とする次のような連立一次方程式を解く必要が生じる。
Ax=b (1)
ここで、Aは一般にn×nの非対称行列、xはn次元の変数ベクトル、bはn次元の定数ベクトルである。nの値が10000以上になることも珍しくない。
Ax=b (1)
ここで、Aは一般にn×nの非対称行列、xはn次元の変数ベクトル、bはn次元の定数ベクトルである。nの値が10000以上になることも珍しくない。
大規模な連立一次方程式は、気象予測、原子炉設計、半導体の回路解析、航空工学における流体解析、構造物の構造解析等の多くの科学技術計算に用いられる。また、大規模な固有値問題は、構造物の構造解析、回路解析、地球科学における地震予知、原子炉の安全性解析、分子科学における多電子系のエネルギー計算、原子核の構造解析等の分野において、物理系の固有振動を記述するときに現れる。
したがって、(1)式のような大規模な連立一次方程式を効率よく高速に解くことは、科学技術計算の重要な問題の1つである。今日では、計算を高速化するために、複数のプロセッシング・エレメント(PE)を備えたメモリ分散型の並列計算機が多く用いられている。
計算機を用いて(1)式を解く1つの方法として、AをLU分解するガウスの消去法に基づいた直接法がある。しかし、Aが大きなスパース行列の場合、非零要素が各行に数個しかないこともあり、計算コストや記憶領域の面で無駄が多い。そこで、単純な行列ベクトル積を繰り返して近似解を求める反復解法が多く用いられている。
反復解法の多くはクリロフ(Krylov)部分空間法に帰着される。今、任意のベクトルr0 にAを次々と乗じていくと、r0,Ar0 ,...,Ak-1 r0 のようなベクトル列が生成される。これらの一次独立なベクトルにより張られる空間はクリロフ部分空間と呼ばれ、(1)式の近似解xkをこれらのベクトルの一次結合で記述する一群の反復解法はクリロフ部分空間法と呼ばれる。
このクリロフ部分空間法としては、CG(Conjugate Gradient)法、BCG(Bi-Conjugate Gradient )法、CR(Conjugate Residuals )法、GCR(Generalized Conjugate Residuals )法、MGCR(Modified Generalized Conjugate Residuals)法、GMRES(Generalized Mainimal RESidual )法等がある。
ところで、(1)式のスパース行列Aの形は与えられた問題によって様々であり、その非零要素を集めて配列に格納する方法を問題によって適当に選ぶ必要がある。しかし、データ格納方法が変わればその配置形態も変わるため、スパース行列とベクトルの積を並列に計算するアルゴリズムも変更する必要が生じる。このため、従来の反復解法では、多様なデータ格納方法の中から問題に適したもの選び、それに応じて演算アルゴリズムを個別に作成している。
したがって、同じ反復解法を用いていても、問題が変わる度に新しくプログラムを作成しなければならず、汎用性に欠けるという問題がある。また、この方法ではプログラマの負担も大きくなる。そこで、連立一次方程式を並列計算機で解く際に、各種のデータ格納方法をサポートしつつ、反復解法を実現することが望まれる。
また、反復解法で繰り返し現れる行列ベクトル積の並列演算を行う際、PE間のデータ転送をできるだけ少なくして、演算を高速化することが重要である。
本発明は、クリロフ部分空間を利用した反復解法により連立一次方程式を解くメモリ分散型並列計算機において、多様なデータ格納方法に対応して効率的な並列処理を行う並列処理装置およびその方法を提供することを目的とする。
図1は、本発明の並列処理装置の原理図である。図1の並列処理装置は、スパース行列を係数行列とする連立方程式を反復解法により解く並列計算機に設けられ、反復手段1、並列演算手段2、および中間配列記憶手段3を備える。
並列演算手段2は、上記係数行列を複数部分に分割して格納し、係数行列の格納形式に依存する計算を並列に行う。
中間配列記憶手段3は、並列演算手段2による計算結果を中間配列として、分割して記憶する。
中間配列記憶手段3は、並列演算手段2による計算結果を中間配列として、分割して記憶する。
反復手段1は、上記中間配列のデータを用いて、上記係数行列の格納形式とは独立に処理を実行し、反復処理を制御する。
反復手段1は、例えば、上記並列計算機が有する各PEに備えられ、与えられた初期ベクトルを用いて、行列ベクトル積を利用した反復処理を開始する。このとき、初期ベクトルに基づいて入力ベクトルを生成し、並列演算手段2に与える。
反復手段1は、例えば、上記並列計算機が有する各PEに備えられ、与えられた初期ベクトルを用いて、行列ベクトル積を利用した反復処理を開始する。このとき、初期ベクトルに基づいて入力ベクトルを生成し、並列演算手段2に与える。
並列演算手段2は、例えば、複数のPEに対応し、係数行列の格納形式に応じて係数行列と入力ベクトルの積を並列に計算して、中間ベクトルを生成する。このような行列ベクトル積の演算は、反復解法の各繰り返し段階において、少なくとも1回以上必要になる。
生成された中間ベクトルは、中間配列記憶手段3内の中間配列に格納される。中間配列記憶手段3は、例えば、複数のPEの主記憶に対応し、中間ベクトルはこれらの主記憶に分割されて格納される。
反復手段1は、中間ベクトルを用いて次の入力ベクトルを生成し、それを並列演算手段2に与える。そして、このような処理の繰り返しにより解が収束すると、反復処理を終了して得られた解ベクトルを出力する。
中間配列を用いることにより、反復手段1は、係数行列の格納形式とは無関係に次の入力ベクトルを生成し、収束判定を行うことができる。したがって、反復手段1の処理を係数行列のデータ格納方法から完全に独立させることができる。これにより、特定のデータ格納方法を前提としない反復処理の実装が可能になり、反復解法の汎用性が高まる。
データ格納方法に依存する行列ベクトル積等の計算は、必要に応じて並列演算手段2に依頼し、反復手段1は与えられた反復解法の中核アルゴリズムのみを実行すればよい。一方、並列演算手段2は、採用された反復解法の種類とは無関係に、与えられた行列ベクトル積等の単純な計算を繰り返すだけである。こうして、効率のよい反復処理が実現される。
データ格納方法に依存する処理かどうかは、例えば、各処理の処理コードに付加されたパラメータの値等により識別することができる。
このような並列処理方法では、データ格納方法が変わっても、反復解法の中核アルゴリズムをプログラミングし直す必要がないので、プログラマの作業が大幅に削減される。
このような並列処理方法では、データ格納方法が変わっても、反復解法の中核アルゴリズムをプログラミングし直す必要がないので、プログラマの作業が大幅に削減される。
例えば、図1の反復手段は、実施形態の図2における個々のPE33または図17における反復部51に対応し、並列演算手段2は、図2における複数のPE33または図17における行列ベクトル積演算部52に対応する。
本発明によれば、スパース行列の連立一次方程式を解くメモリ分散型並列計算機において、同一の反復解法のアルゴリズムを、多様なスパース行列の格納法に適用することが可能になる。したがって、反復解法の汎用性が向上し、プログラマの負担が軽減される。
また、スパース行列をバンド行列とみなして処理することで、行列ベクトル積の演算に必要なデータ転送量が削減され、並列処理の台数効果が向上する。特に、反復解法として他の解法より高速なMGCR法を用いた場合、さらに効率のよい並列処理が実現される。
以下、図面を参照しながら、本発明の実施の形態を詳細に説明する。
本発明においては、反復解法において必要となる行列ベクトル積等の処理を行うために、ベクトルを格納する1次元配列を各PE(プロセッサ)で均等に分割し、この配列を中間インタフェースとして利用して、中核の計算アルゴリズムをスパース行列のデータ格納方法から独立させる。つまり、反復解法の各段階に現れる中間ベクトルを各PEで均等に分割した中間領域に格納する。
本発明においては、反復解法において必要となる行列ベクトル積等の処理を行うために、ベクトルを格納する1次元配列を各PE(プロセッサ)で均等に分割し、この配列を中間インタフェースとして利用して、中核の計算アルゴリズムをスパース行列のデータ格納方法から独立させる。つまり、反復解法の各段階に現れる中間ベクトルを各PEで均等に分割した中間領域に格納する。
これにより、データ格納方法に依存する行列ベクトル積等の処理が反復解法の中核アルゴリズムから分離され、多様なデータ格納方法に対処できるようになる。
また、代表的な行列のデータ格納法である対角形式格納法またはエルパック(Ellpack )形式格納法を用いて、行列ベクトル積を並列に計算することができる。このとき、スパース行列を疑似的なバンド行列とみなしてそのバンド幅を求め、これを利用して行列ベクトル積を計算すれば、各PEは近くの限られたPEとのみデータ転送を行えばよくなる。したがって、転送回数が削減され、色々なバンド幅のスパース行列について、行列ベクトル積の演算が効率よく実行される。
また、代表的な行列のデータ格納法である対角形式格納法またはエルパック(Ellpack )形式格納法を用いて、行列ベクトル積を並列に計算することができる。このとき、スパース行列を疑似的なバンド行列とみなしてそのバンド幅を求め、これを利用して行列ベクトル積を計算すれば、各PEは近くの限られたPEとのみデータ転送を行えばよくなる。したがって、転送回数が削減され、色々なバンド幅のスパース行列について、行列ベクトル積の演算が効率よく実行される。
さらに、一般スパース行列の連立1次方程式の反復解法の1つであるMGCR法を例に取り、本発明の反復処理をメモリ分散型のスーパーコンピュータVPPシリーズ向けに構成した。この形態によれば、各種のデータ格納方法に対してMGCR法を適用し、高速な計算処理を実現することが可能になる。
まず、図2および図3を参照しながら、実施形態で用いる並列計算機の構成を説明する。図2は、並列計算機の概略構成図である。図2の並列計算機は、入出力装置11と、IOバス12と、クロスバ方式により互いに接続された複数のPE13とを備える。
入出力装置11は、各PE13との間でデータの入出力を行うための装置で、例えば、キーボードなどの入力機器とディスプレイやプリンタなどの出力機器を備えた計算機端末である。入出力装置11は、IOバス12により各PE13と接続されている。
図3は、PE13の構成図である。図5のPE13は、主記憶21、IOポート22、クロスバスイッチ23、データ転送制御ユニット24、メモリ制御ユニット25、スカラユニット38、およびベクトルユニット39を備える。
主記憶21は、与えられた問題の行列要素や計算の中間結果などを格納する。IOポート22はIOバス12に接続され、入出力装置11との間で入出力データの転送を行う。クロスバスイッチ23は、主記憶21内のデータを他のPE13に転送する際のデータの切り換え等を行う。データ転送制御ユニット24は、クロスバスイッチ23を介して他のPE13とのデータ転送および同期制御を行う。
メモリ制御ユニット25は、データ転送制御ユニット24、スカラユニット38、およびベクトルユニット39から発生するアクセス要求を受け取り、主記憶21に対するアクセスを制御する。
スカラユニット38は、キャッシュメモリ26、GPR(汎用レジスタ)/FPR(浮動小数点レジスタ)27、スカラ演算器28を備え、スカラ演算を実行する。
ベクトルユニット39は、マスクレジスタ29とベクトルレジスタ30の他に、ロードパイプライン31、ストアパイプライン32、マスクパイプライン33、34、乗算パイプライン35、加算/論理演算パイプライン36、および除算パイプライン37を備える。
ベクトルユニット39は、マスクレジスタ29とベクトルレジスタ30の他に、ロードパイプライン31、ストアパイプライン32、マスクパイプライン33、34、乗算パイプライン35、加算/論理演算パイプライン36、および除算パイプライン37を備える。
これらの各パイプラインは、それぞれ複数要素を同時に処理することができ、パイプライン35、36、37のうち2本は同時に動作できる。また、他のパイプラインについては、すべて同時に動作可能である。マスクパイプライン33は総和/検索処理用に使用され、マスクパイプライン34は論理演算処理用に使用される。
ベクトルユニット39は、これらのパイプラインを用いて複数のベクトル命令を並列に実行することができ、高速な行列演算が可能である。
次に、図4および図5を参照しながら、本発明の反復処理について説明する。図4は、多様なデータ格納方法に対応可能な反復処理のフローチャートである。図4において処理が開始されると、並列計算機の各PEは、まずあらかじめ決められた反復法による処理コードを取り出す(ステップS1)。
次に、図4および図5を参照しながら、本発明の反復処理について説明する。図4は、多様なデータ格納方法に対応可能な反復処理のフローチャートである。図4において処理が開始されると、並列計算機の各PEは、まずあらかじめ決められた反復法による処理コードを取り出す(ステップS1)。
そして、そのコードに記述された計算がデータ格納法に依存するかどうかを判定し(ステップS2)、データ格納法に依存しなければ各PEで同じ処理を実行する(ステップS3)。データ格納法に依存するかどうかは、例えば、あらかじめ処理コード内に書き込まれたパラメータ値により判断する。
次に、解が十分に収束するなどの反復終了条件が満たされたかどうかを判定し(ステップS5)、それが満たされなければステップS1以降の処理を繰り返す。
また、ステップS2において計算がデータ格納法に依存する場合は、他のPEと並列処理を行って、計算結果を中間配列に格納する(ステップS4)。この中間配列は、各種のデータ格納法に使用できるように標準的に用意された領域で、各PEにより分割されている。例えば、PE1、PE2、PE3、PE4の4台のPEを備える場合は、中間配列は図5に示すようになる。
また、ステップS2において計算がデータ格納法に依存する場合は、他のPEと並列処理を行って、計算結果を中間配列に格納する(ステップS4)。この中間配列は、各種のデータ格納法に使用できるように標準的に用意された領域で、各PEにより分割されている。例えば、PE1、PE2、PE3、PE4の4台のPEを備える場合は、中間配列は図5に示すようになる。
その後、各PEはステップS5以降の処理を繰り返し、反復終了条件が満たされれば処理を終了する。
このような反復処理によれば、例えば行列ベクトル積に代表されるデータ格納法依存性の高い並列演算の結果を中間配列に格納して、反復法における他の計算処理をデータ格納法と独立に行うことができる。これにより、各種スパース行列の格納法を用いて、スパース行列の連立1次方程式を容易に解くことが可能になる。
このような反復処理によれば、例えば行列ベクトル積に代表されるデータ格納法依存性の高い並列演算の結果を中間配列に格納して、反復法における他の計算処理をデータ格納法と独立に行うことができる。これにより、各種スパース行列の格納法を用いて、スパース行列の連立1次方程式を容易に解くことが可能になる。
次に、図6から図12までを参照しながら、具体的なデータ格納法に基づく効率のよい行列ベクトル積の演算方法を説明する。
図6は、スパース行列の対角形式格納法を示している。ここでは、n次のスパース行列Aをバンド行列の1種とみなし、対角線方向に並んだ非零要素から成る対角ベクトルを2次元の格納配列41に格納する。
図6は、スパース行列の対角形式格納法を示している。ここでは、n次のスパース行列Aをバンド行列の1種とみなし、対角線方向に並んだ非零要素から成る対角ベクトルを2次元の格納配列41に格納する。
Aの対角ベクトル部分を拡大すると図7に示すようになり、対角ベクトル間には対角方向の要素がすべて0の部分が存在する。このような部分は格納配列41に格納しなくてもよいので、対角方向に集中して非零要素が存在するスパース行列に特に適した格納法であると言える。
図6および図7に示すように、Aをバンド行列とみなしたときの上バンド幅bandwuと下バンド幅bandwlは、それぞれ対角要素から対角ベクトルまでの距離の最大値として求められる。このとき、格納配列41の幅wは最大bandwu+bandwl+1となる。
また、図6において対角ベクトルの長さを揃えるために、斜線部分に0の要素を付け加えて格納配列41に格納しておく。これにより、Aの対角ベクトルを用いた演算を均一に行うことができる。
配列41の1次元目(行の次元)を各PEに均等に分割して格納し、さらにn次元のベクトルX、Yを格納する配列を各PEに分割配置することにより、次式の行列ベクトル積を並列に計算することができる。
Y=AX (11)
例えば、4台のPEで並列演算を行う場合は、Aの配列41、Xの配列42、Yの配列43は、それぞれ図8に示すように分割される。このうち、演算結果を格納する配列43は図5の中間配列に相当する。
Y=AX (11)
例えば、4台のPEで並列演算を行う場合は、Aの配列41、Xの配列42、Yの配列43は、それぞれ図8に示すように分割される。このうち、演算結果を格納する配列43は図5の中間配列に相当する。
ところで、図6のAは一定のバンド領域にのみ非零要素を持つスパース行列であるため、(11)式の演算において、各PEは必ずしもXのすべての要素を持っている必要はない。例えば図8のPE2が受け持つ行列ベクトル積の有効部分は、図9に示すようになる。
図9においてPE2に割り当てられたAの行数をn2とすると、その非零要素が存在する領域はAの斜線部分である。この斜線部分に該当する列の範囲は、n2本の行と同じ番号を持つn2本の列と、その左側のbandwlの幅の各列と、右側のbandwuの幅の各列となる。したがって、Xの要素のうち、これらの列に対応する(n2+bandwl+bandwu)個の要素のみが乗算の結果に寄与する。他の要素には0が乗算されるので、その結果もまた0となる。
Xを格納する配列42は、図8に示すように各PEに分散されているため、PE2は乗算を行うために、少なくともbandwl個の要素をPE1からコピーし、少なくともbandwu個の要素をPE3からコピーする必要がある。他のPEについても同様のコピー処理が必要になる。
そこで、あらかじめPE毎にワーク用の一次元配列Wを用意しておき、この配列に必要なXの要素をコピーすることにする。例えばPE2の場合、Wの長さとしては、最低(n2+bandwl+bandwu)だけあれば十分である。しかし、ここではXの各要素の論理的な添え字を、そのままWの対応する要素の添え字として使用できるように、各PEのWの長さを(n+bandwl+bandwu)に統一することにする。
こうして、図10に示すようなWが各PEに用意される。ここで、PEの台数をP、i番目のPEに割り当てられたAの行数をni(i=1,...,P)とすると、i番目のPEでは、Xの必要な要素が図10に示すようにWの斜線部分にコピーされる。ここで、niはおおよそn/Pとなる。Wの斜線部分以外の領域は実際の処理には用いられないので、どんな値を格納していてもかまわない。
Xの斜線部分の要素のうち、元々自分が持っているni個の要素はそのまま主記憶21上でコピーすることができ、その上下のbandwl個とbandwu個の要素は、クロスバスイッチ23を介して隣接PEまたは他の近くのPEからコピーする。配列41に格納された対角ベクトルの数がAの次元nに比べて十分に小さい場合は、高々隣接PEとの通信だけで必要な要素をコピーすることができる。
このように、図6のような対角方向に非零要素が集中しているスパース行列の場合は、特定のPE間の通信だけで行列ベクトル積を実行でき、通信に伴うオーバヘッドを大幅に削減することができる。
また、1番目のPEは、Wの上端部の長さbandwlの部分にダミーデータを格納して乗算を行い、P番目のPEは、下端部の長さbandwuの部分にダミーデータを格納して乗算を行う。これらのダミーデータには、図6の斜線部分の要素である0が乗算されるので、ダミーデータの値が行列ベクトル積の結果に影響を与えることはない。しかし、ダミーデータを用いることで、1番目とP番目のPEは他のPEと同様の乗算を行うことができる。
このような乗算方法によれば、各PE内の演算はベクトル演算器を用いて効率よく処理することができる。また、各PEは演算の途中でベクトル長を変更する必要もない。
図11は、バンド幅を利用したスパース行列の行列ベクトル積演算のフローチャートである。図11において処理が開始されると、各PEは、まずスパース行列をバンド行列とみなして、その上バンド幅と下バンド幅を求める(ステップS11)。次に、分割されたベクトルXの要素のうち、自分が保持する部分(図10のni個の要素)をワーク用配列Wに並列にコピーする(ステップS12)。
図11は、バンド幅を利用したスパース行列の行列ベクトル積演算のフローチャートである。図11において処理が開始されると、各PEは、まずスパース行列をバンド行列とみなして、その上バンド幅と下バンド幅を求める(ステップS11)。次に、分割されたベクトルXの要素のうち、自分が保持する部分(図10のni個の要素)をワーク用配列Wに並列にコピーする(ステップS12)。
次に、各PEは、必要な下バンド幅の長さの部分(図10のbandwl個の要素)をWに並列にコピーする(ステップS13)。また、必要な上バンド幅の長さの部分(図10のbandwu個の要素)をWに並列にコピーする(ステップS14)。
そして、行列ベクトル積を並列に計算して、計算結果をベクトルYの各PEに割り当てられた部分に格納し(ステップS15)、処理を終了する。
この方法によれば、ステップS13およびS14において限られたPEとのみ通信を行えばよいので、データ転送量が少なくて済み、行列ベクトル積演算が高速化される。
この方法によれば、ステップS13およびS14において限られたPEとのみ通信を行えばよいので、データ転送量が少なくて済み、行列ベクトル積演算が高速化される。
このようなバンド幅を利用した行列ベクトル積の計算方法は、対角形式格納法に限らず、他の任意のデータ格納法にも適用可能である。他のデータ格納法の1つとして、エルパック形式格納法がある。
エルパック形式格納法では、スパース行列の各行ベクトル毎に非零要素を圧縮した形で、行列要素を格納する。このとき、各行内の非零要素に対応する元のスパース行列の列番号は、格納配列icofに格納されている。この列番号の値を利用すれば、バンド幅を求めることができる。
図12は、配列icofの例を示している。図12において、左のスパース行列は10000×10000の行列で、対角要素と他の対角方向の線分上に非零要素を持っている。また、icofは、スパース行列の第i行に含まれる非零要素の列番号icof(i,*)を、各行番号i毎に格納している。ここで、*は第i行内の1つ以上の非零要素に付けられた識別番号である。
スパース行列のすべての非零要素のうち、対角要素から最も離れたものとその行の対角要素との距離を求めれば、バンド幅が得られる。
例えば、i=9000の行内には、5000番目、6000番目、および9000番目の各列に非零要素があり、これらの列番号がicof(9000,*)としてicofに格納されている。このうち、9000番目の列に対角要素があり、それから左に最も離れた要素は5000番目の列にある。したがって、これらの要素間の距離は−(5000−9000)=4000となる。
例えば、i=9000の行内には、5000番目、6000番目、および9000番目の各列に非零要素があり、これらの列番号がicof(9000,*)としてicofに格納されている。このうち、9000番目の列に対角要素があり、それから左に最も離れた要素は5000番目の列にある。したがって、これらの要素間の距離は−(5000−9000)=4000となる。
また、この5000番目の列の要素はバンドの最も外側に位置しているため、この場合の下バンド幅bandwlは4000となる。上バンド幅bandwuもicofから同様にして求められる。
icof(i,*)を用いたバンド幅の計算式を一般化すると、次式のようになる。
(12)、(13)式において、*に関するMaxは、第i行内の*番目の要素と対角要素との距離の最大値を表し、iに関するMaxは、各行の*に関する最大値のうちの最大値を表す。
(12)式で、(icof(i,*)−i)の前に負の符号が付いているのは、左に位置する要素ほど対角要素までの距離を大きく評価するためである。この場合、対角要素より右に位置する要素までの距離は負の値を持つことになる。逆に、(13)式では、右に位置する要素ほど対角要素までの距離が大きくなり、対角要素より左に位置する要素までの距離は負の値を持つ。
こうして、エルパック形式格納法においても上下バンド幅を計算することが可能になり、図11の演算処理を適用することができる。
次に、図13から図17までを参照しながら、最も高速な反復解法の1つであるMGCR法の実施形態について説明する。MGCR法は、各繰り返し段階における演算回数がGMRES法より少ないため、この方法より高速であることが知られている(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)。
次に、図13から図17までを参照しながら、最も高速な反復解法の1つであるMGCR法の実施形態について説明する。MGCR法は、各繰り返し段階における演算回数がGMRES法より少ないため、この方法より高速であることが知られている(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)。
GMRES法は反復解法の中でもかなり速いアルゴリズムである。それより速いMGCR法を実装することで、従来にない高速な反復処理装置を実現することができる。ここでは、一例としてスーパーコンピュータVPP500を対象とした実装方法を説明する。
(1)式でx=u,b=fとおき直し、与えられた連立一次方程式を、
Au=f (14)
と書くことにすると、ベクトルuの近似解を求めるMGCR法のアルゴリズムは次のようになる。
Au=f (14)
と書くことにすると、ベクトルuの近似解を求めるMGCR法のアルゴリズムは次のようになる。
ここで、(15)式のu0 は任意に与えられた初期ベクトルであり、|r1 |はベクトルr1の大きさを表す。また、s1 =r1 は1番目の探索ベクトルとなる。
(16)〜(23)式のループ処理では、i=1,...,kについてk個の探索ベクトルsi+1 を生成し、(23)式でdi+1 が、与えられた収束判定値ε以内に収束すれば処理を終了する。そうでなければ、(24)〜(25)式によりu0を更新して、(15)式からの始まる外側のループ処理を繰り返す。kの値は任意に指定可能であるが、通常は10〜100の範囲に設定すればよい。
(16)〜(23)式のループ処理では、i=1,...,kについてk個の探索ベクトルsi+1 を生成し、(23)式でdi+1 が、与えられた収束判定値ε以内に収束すれば処理を終了する。そうでなければ、(24)〜(25)式によりu0を更新して、(15)式からの始まる外側のループ処理を繰り返す。kの値は任意に指定可能であるが、通常は10〜100の範囲に設定すればよい。
(16)、(21)式において、(v1,v2)のような表記はベクトルv1とv2の内積を表し、(19)、(23)、(24)式において、solve Hc=αのような表記は、連立方程式Hc=αを解いてcを求めることを意味する。
Hk は、次式で与えられるk次の上三角行列である。
Hk の要素βj,i-1 (i=1,...,k, j=1,...,i)は、(16)式または(17)、(18)式により、Aおよびsiから計算される。Hi-1 ,Hi についても同様である。
αk は、αk =(α1 ′,...,αk′)T のようなk次元の定数ベクトルである。このベクトルの要素αi ′は、(21)式によりs1とsi+1 から計算される。αi-1 ,αi についても同様である。
また、(25)式におけるNk は次式で与えられる。
Nk =(s1 ,s2 ,...,sk) (28)
(19)式のNi-1 、(23)式のNi も同様である。
Nk =(s1 ,s2 ,...,sk) (28)
(19)式のNi-1 、(23)式のNi も同様である。
このMGCR法のアルゴリズムを並列計算機上に実装するとき、(15)式のA×u0 や(16)、(17)式のA×si のようなスパース行列とベクトルの積を格納する領域が必要になる。そこで、各PEにより分割された2次元の中間配列W1を用意し、Y=AXの形の行列ベクトル積のXとYを格納する。また、(17)式により生成されるベクトルsi(i=1,...,k)もW1に格納しておく。
4台のPEに分散配置されたW1は図13に示すようになる。図13において、W1はn×(k+2)の大きさを持ち、第1列にはXが格納され、第2列にYが格納される。また、残りのk本の列には各si が格納される。第1列または第2列は、(15)式のr1を格納する時にも使用される。
図14は、各PE毎に設けられる格納配列を示している。図14の配列W2はk×(k+1)の大きさを持ち、最初のk本の列にHk が格納され、最後の列にαk が格納される。実際には、Hk用の領域にはHi (i=1,...,k)が順次格納されて処理に用いられ、αk 用の領域にはαi(i=1,...,k)が順次格納されて処理に用いられる。
この他に、スパース行列Aを格納する配列が各PEに分散配置され、図11のような行列ベクトル積の演算を行うために、図10のようなワーク用配列が各PEに1つ設けられる。スパース行列Aは、例えば対角形式格納法やエルパック形式格納法により格納されるが、もちろん、他の任意のデータ格納法を採用することもできる。
図15は、これらの格納配列を用いたMGCR法による処理の概要を示している。図15において、W1は4台のPEに分割されている。
これらのPE1、PE2、PE3、PE4はW1を利用して、必要に応じてA×u0 またはA×si を並列に計算する。その後、各PEはW2を利用して他の計算を同時に行い、u0を更新していく。そして、(19)式または(23)式の条件が成り立てば、処理を終了する。
これらのPE1、PE2、PE3、PE4はW1を利用して、必要に応じてA×u0 またはA×si を並列に計算する。その後、各PEはW2を利用して他の計算を同時に行い、u0を更新していく。そして、(19)式または(23)式の条件が成り立てば、処理を終了する。
図16は、このようなMGCR法による処理のフローチャートである。図16において処理が開始されると、並列計算機は、まず与えられた初期ベクトルu0 を用いた並列処理によりA×u0 を計算し、(15)式のr1とd1 を求める(ステップS21)。A×u0 の行列ベクトル積は図11の方法で実行される。
次に、i=1とおいて、MGCR法による処理を開始し(ステップS22)、図11の方法によりA×si を並列に計算する(ステップS23)。次に、各PEでA×siの結果を用いて、同じ(16)、(17)、(18)式の計算を冗長に行い、βi,i-1 を求める(ステップS24)。そして、βi,i-1が0かどうかを調べる(ステップS25)。
βi,i-1 が0でなければ、各PEでβi,i-1 の値を用いて、同じ(20)、(21)、(22)式の計算を冗長に行い、si+1およびdi+1 を求める(ステップS26)。そして、di+1 がε以下になったかどうかを調べ、収束判定を行う(ステップS27)。
di+1 がεより大きければ、次にiがkに達したかどうかを調べる(ステップS28)。そして、i=kでなければi=i+1とおいて(ステップS29)、ステップS23以降の処理を繰り返す。また、ステップS28においてi=kとなった場合は、各PEで(24)、(25)、(26)式の再スタート処理を行ってu0を更新し(ステップS30)、ステップS21以降の処理を繰り返す。
そして、ステップS25においてβi,i-1 =0となった場合は、各PEは(19)式の第1の終了処理を行って近似解ui-1を求め(ステップS31)、処理を終了する。また、ステップS27においてdi+1 がε以下に収束した場合は、各PEは(23)式の第2の終了処理を行って近似解uiを求め(ステップS32)、処理を終了する。
上述の第1の終了処理、第2の終了処理、および再スタート処理において、最大k次の連立一次方程式を解く必要がある。しかし、kは高々100の程度で、多くの場合nに比べてはるかに小さいので、各PE内で逐次的に解いてもそれほど時間はかからない。
以上のようなMGCR法を実装した並列計算機を模式的に表すと、図17のように書くことができる。図17において、反復部51は、スパース行列Aのデータ格納法に依存しない計算処理や判定処理を行い、行列ベクトル積演算部52は、A×u0 やA×si の行列ベクトル積の計算処理を行う。
これらの行列ベクトル積の計算はAの格納方法に依存するため、行列ベクトル積演算部52に処理コードを与えて、入力ベクトルXや出力ベクトルYの格納場所等を指示する。このような制御用の処理コードは、コンパイラまたはプログラマにより作成される。
図17のような構成を取ることにより、行列ベクトル積の演算をMGCR法から分離して独立に実行することができ、各種のスパース行列格納法に対して、MGCR法を適用することが可能になる。
本発明で用いる反復解法はMGCR法に限られるわけではなく、クリロフ部分空間を利用した任意の解法を用いることができる。また、図17の構成において、反復部51の処理をいずれかのPEに代表して行わせることも可能である。
1 反復手段
2 第1の並列処理手段
3 中間配列記憶手段
11 入出力装置
12 IOバス
13 PE
21 主記憶
22 IOポート
23 クロスバスイッチ
24 データ転送制御ユニット
25 メモリ制御ユニット
26 キャッシュメモリ
27 GPR/FPR
28 スカラ演算器
29 マスクレジスタ
30 ベクトルレジスタ
31 ロードパイプライン
32 ストアパイプライン
33、34 マスクパイプライン
35 乗算パイプライン
36 加算/論理演算パイプライン
37 除算パイプライン
38 スカラユニット
39 ベクトルユニット
41、42、43、W、W1、W2、icof 格納配列
51 反復部
52 行列ベクトル積演算部
2 第1の並列処理手段
3 中間配列記憶手段
11 入出力装置
12 IOバス
13 PE
21 主記憶
22 IOポート
23 クロスバスイッチ
24 データ転送制御ユニット
25 メモリ制御ユニット
26 キャッシュメモリ
27 GPR/FPR
28 スカラ演算器
29 マスクレジスタ
30 ベクトルレジスタ
31 ロードパイプライン
32 ストアパイプライン
33、34 マスクパイプライン
35 乗算パイプライン
36 加算/論理演算パイプライン
37 除算パイプライン
38 スカラユニット
39 ベクトルユニット
41、42、43、W、W1、W2、icof 格納配列
51 反復部
52 行列ベクトル積演算部
Claims (4)
- 連立方程式を反復解法により解く並列計算機において、
係数行列と初期ベクトルから生成される探索ベクトルを格納する中間配列を、複数部分に分割して記憶する中間配列記憶手段と、
前記係数行列および探索ベクトルから生成された要素を含み、解ベクトルを求めるために用いられるMGCR法の三角行列を記憶する行列記憶手段と、
前記係数行列を用いて行列ベクトル積を並列に計算し、該行列ベクトル積を用いて前記探索ベクトルを更新する並列演算手段と
を備えることを特徴とする並列処理装置。 - 前記中間配列記憶手段は、さらに行列ベクトル積の結果を分割して記憶することを特徴とする請求項1記載の並列処理装置。
- 前記行列記憶手段は、さらに前記三角行列の次元の定数ベクトルを記憶し、前記並列演算手段は、該三角行列と定数ベクトルとから形成される連立方程式を解き、得られた結果を用いて前記解ベクトルを求めることを特徴とする請求項1記載の並列処理装置。
- 連立方程式を反復解法により解く計算機において、
係数行列と初期ベクトルから生成される探索ベクトルを、中間配列として記憶する中間配列記憶手段と、
前記係数行列および探索ベクトルから生成された要素を含み、解ベクトルを求めるために用いられるMGCR法の三角行列を記憶する行列記憶手段と、
前記係数行列を用いて行列ベクトル積を計算し、該行列ベクトル積を用いて前記探索ベクトルを更新する演算手段と
を備えることを特徴とする処理装置。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2006031725A JP2006164307A (ja) | 2006-02-08 | 2006-02-08 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2006031725A JP2006164307A (ja) | 2006-02-08 | 2006-02-08 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Related Parent Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP01610596A Division JP3808925B2 (ja) | 1996-01-31 | 1996-01-31 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2006164307A true JP2006164307A (ja) | 2006-06-22 |
Family
ID=36666161
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2006031725A Pending JP2006164307A (ja) | 2006-02-08 | 2006-02-08 | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP2006164307A (ja) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106897163A (zh) * | 2017-03-08 | 2017-06-27 | 郑州云海信息技术有限公司 | 一种基于knl平台的代数系统求解方法及系统 |
CN112130998A (zh) * | 2020-09-23 | 2020-12-25 | 中国核动力研究设计院 | 适用于多环路压水堆的核反应堆系统分析程序的优化方法 |
-
2006
- 2006-02-08 JP JP2006031725A patent/JP2006164307A/ja active Pending
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106897163A (zh) * | 2017-03-08 | 2017-06-27 | 郑州云海信息技术有限公司 | 一种基于knl平台的代数系统求解方法及系统 |
CN112130998A (zh) * | 2020-09-23 | 2020-12-25 | 中国核动力研究设计院 | 适用于多环路压水堆的核反应堆系统分析程序的优化方法 |
CN112130998B (zh) * | 2020-09-23 | 2022-02-01 | 中国核动力研究设计院 | 适用于多环路压水堆的核反应堆系统分析程序的优化方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US5887186A (en) | Method of solving simultaneous linear equations in a memory-distributed parallel computer | |
CN114897132A (zh) | 在硬件中执行核心跨越 | |
CN115271049A (zh) | 硬件中转置神经网络矩阵 | |
CN111310904A (zh) | 一种用于执行卷积神经网络训练的装置和方法 | |
CN113344172A (zh) | 将卷积映射到通道卷积引擎 | |
US10867008B2 (en) | Hierarchical Jacobi methods and systems implementing a dense symmetric eigenvalue solver | |
US20200226201A1 (en) | Methods and Apparatus for Constructing Digital Circuits for Performing Matrix Operations | |
KR20220038579A (ko) | 데이터 처리 | |
JP2023070746A (ja) | 情報処理プログラム、情報処理装置、及び情報処理方法 | |
JP2006164307A (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
JP3808925B2 (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
JP2006164306A (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
US20180349321A1 (en) | Parallel processing apparatus, parallel operation method, and parallel operation program | |
JP3659307B2 (ja) | ベクトル・コンピュータにおける演算方法、及び記録媒体 | |
CN111522776B (zh) | 一种计算架构 | |
JP6666548B2 (ja) | 並列計算機、fft演算プログラムおよびfft演算方法 | |
JP2008269329A (ja) | 連立一次方程式の解を反復的に決定する方法 | |
Kamra et al. | A stable parallel algorithm for block tridiagonal toeplitz–block–toeplitz linear systems | |
Tokura et al. | Gpu-accelerated bulk computation of the eigenvalue problem for many small real non-symmetric matrices | |
Kuznetsov | An approach of the QR factorization for tall-and-skinny matrices on multicore platforms | |
JP2806262B2 (ja) | マルチプロセッサシステムのプロセス割当方法 | |
Brown | A comparison of techniques for solving the Poisson equation in CFD | |
JP6511937B2 (ja) | 並列計算機システム、演算方法、演算プログラム、及び情報処理装置 | |
Gurung et al. | Simultaneous Solving of Linear Programming Problems in GPU | |
Yamada et al. | Optimization of reordering procedures in hotrg for distributed parallel computing |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20071004 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20071016 |
|
A02 | Decision of refusal |
Effective date: 20080311 Free format text: JAPANESE INTERMEDIATE CODE: A02 |