JP3639206B2 - 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 - Google Patents
共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 Download PDFInfo
- Publication number
- JP3639206B2 JP3639206B2 JP2000358232A JP2000358232A JP3639206B2 JP 3639206 B2 JP3639206 B2 JP 3639206B2 JP 2000358232 A JP2000358232 A JP 2000358232A JP 2000358232 A JP2000358232 A JP 2000358232A JP 3639206 B2 JP3639206 B2 JP 3639206B2
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- block
- decomposition
- blocks
- recording medium
- 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
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Complex Calculations (AREA)
Description
【発明の属する技術分野】
本発明は、共有メモリ型スカラ並列計算機における並列行列処理方法に関する。
【0002】
【従来の技術】
連立1次方程式を計算機によって解く場合には、連立1次方程式を行列表示し、この行列について処理を施すことによって、解の求めやすい形に変形し、このような形にしてから方程式の解を求める方法が採られている。
【0003】
すなわち、連立1次方程式は、係数を表す行列と、変数を表す列ベクトルとの積が所定の列ベクトルに等しくなるというように記述することが出来る。ここで、LU分解(LU Factorization)という方法によれば、係数を表す行列を上三角行列(upper-triangular matrix )と下三角行列(lower-triangular matrix )に分解することによって、連立1次方程式の解を求める。したがって、この場合においては、係数行列をLU分解することが連立1次方程式の解を得るために重要な処理となる。また、LU分解の特別な場合として(変形)コレスキー分解(Cholesky Factorization)という行列の分解方法がある。
a)実行列の連立1次方程式の解法
実行列の連立1次方程式の解法に関して、ベクトル並列計算機での連立1次方程式は、ブロック化した外積型のLU分解をベースに並列化を行っている。列ベクトルを何本か束ねたブロックをLU分解(1)した後、対応した行ベクトルを束ねたブロックを更新して(2)から、正方小行列を更新(3)する処理を繰り返す。
【0004】
従来、(1)の処理は、一つのプロセッサで逐次的に行っていた。並列効率を高めるためにブロック幅は12(列あるいは行:係数行列の行幅あるいは列幅を示す)程度の比較的小さい値としていた。この結果(2)及び(3)の部分の更新も幅12程度の行列演算となった。
【0005】
最もコストの大きい(3)の計算で、幅が12程度と小さくても、効率の良い方法があった。共有メモリ型スカラ並列計算機(SMP)では、幅が小さいと性能を引き出せない。これは以下の理由による。
【0006】
すなわち、(3)の計算は、行列積である。幅が小さいと更新する行列の要素(メモリに格納されている)をロードして、更新結果をストアするというメモリにアクセスするためのコストが行列の更新を行う演算に比べて大きくなり性能を引き出せない。
【0007】
このため、ブロック幅を大きくする必要があるが、ブロック幅を大きくすると、このブロックのLU分解のコストが大きくなり並列化効率が落ちる。
b)正値対称行列の連立1次方程式の解法
正値対称行列の連立1次方程式の解法に関して、下三角行列部分のみを利用してコレスキー分解を行うときに、分散メモリ型並列計算機では小行列ブロックをcyclicに各プロセッサに分散して負担させ、各プロセッサの負荷を均等にして解いていた。実行列の連立1次方程式と同じようにブロック化するときブロック幅を比較的小さくでき、並列化効率を高めることが可能であった。SMPでは、上記(3)の更新で現れる行列積でブロック幅が大きい方が性能が良いため、ブロック幅を大きくする必要がある。
【0008】
【発明が解決しようとする課題】
すなわち、共有メモリ型スカラ並列計算機では、LU分解あるいはコレスキー分解における行列積による行列の更新に必要とされるコストよりも、共有メモリにアクセスするためのコストが大きくなってしまい、従来ベクトル並列計算機で行っていた方法をそのまま共有メモリ型スカラ計算機に適用しても十分な性能を引き出せない。
【0009】
本発明の課題は、共有メモリ型スカラ計算機に適した並列行列処理方法を提供することである。
【0010】
【課題を解決するための手段】
並列行列処理方法は、複数のプロセッサモジュールを持つ共有メモリ型スカラ並列計算機の行列演算において、行列を小行列ブロックに分けるブロック化ステップと、該小行列ブロックの内、対角ブロックと対角ブロックでない小行列ブロックとを該複数のプロセッサモジュールのローカルメモリ領域に格納する格納ステップと、該複数のプロセッサモジュールが並列に、それぞれ有するブロックを演算することにより、複数のプロセッサモジュールで、対角ブロックを冗長に演算する演算ステップと、該演算ステップで得られた小行列ブロックの演算結果を使って、該行列を更新する更新ステップとを備えることを特徴とする。
【0011】
本発明によれば、複数のプロセッサモジュールがそれぞれ有するローカルメモリ領域を有効に利用して、効率的な行列演算を行うことが出来る。
【0012】
【発明の実施の形態】
図1は、共有メモリ型スカラ並列計算機のハードウェア構成例を示す図である。
【0013】
共有メモリ型スカラ並列計算機は、複数のプロセッサ10−1、10−2、・・・10−nが2次キャッシュメモリ13−1、13−2、・・・13−nを介して相互結合網12に接続される。各プロセッサ10−1、10−2、・・・10−nは、その内部あるいは、2次キャッシュメモリ13−1、13−2、・・・13−nよりプロセッサ側に1次キャッシュメモリが設けられる。また、各プロセッサ10−1、10−2、・・・10−nに共有となっているメモリモジュール11−1、11−2、・・・11−nは、相互結合網12を介してプロセッサ10−1、10−2、・・・10−nがアクセス可能となってる。プロセッサ10−1、10−2、・・・10−nがデータ処理を行う場合には、まず、メモリモジュール11−1、11−2、・・・11−nから1つのプロセッサが担当するデータを2次キャッシュメモリ13−1、13−2、・・・13−nに格納し、更に、2次キャッシュメモリから処理単位となるデータを1次キャッシュメモリにコピーして処理を行う。
【0014】
処理が終わると、1次キャッシュメモリから2次キャッシュメモリに処理データが格納され、2次キャッシュメモリ内のデータが全て処理し終わると、メモリモジュール11−1、11−2、・・・11−nの内、最初にデータを持ってきたメモリモジュールに対してデータの更新を行う。また、次のデータ処理を行う場合には、上述したように、メモリモジュールから各プロセッサが担当する分のデータを2次キャッシュメモリに格納し、1次キャッシュメモリに処理単位のデータを持ってきて、プロセッサが処理を行う。このような処理を繰り返して、並列にデータ処理を完了する。このとき、各プロセッサが処理した後のデータをメモリモジュールに書き込み、次の処理のために、再びメモリモジュールからデータを読み込む際、各プロセッサが自分のタイミングでデータの読み込みを行っていたのでは、データ更新された後のデータを読み込むべきところを、データ更新される前のデータを読み込んでしまう可能性が有る。したがって、このときには、全てのプロセッサがメモリモジュールにデータ更新し終わるまで、他のプロセッサがメモリモジュールからデータを読み込まないようにする必要がある。このように、プロセッサのメモリモジュールからのデータの読み込みを制限して、全体のプロセッサの処理の同期をとることをバリア同期(Barrier Synchronization)を取るという。
【0015】
図2及び図3は、本発明の実施形態に従ったLU分解の並列処理の概念を説明する図である。
図2は、処理すべき行列の模式図であり、図3は、処理単位となるデータの構成を説明する図である。
a)本実施形態に従った実行列の連立1次方程式の解法
本実施形態においては、プロセッサが処理を担当する小行列ブロックの幅を従来より大きくしてLU分解する部分を並列化する。すなわち、図2において、従来では、L1〜L3の部分は、ブロック幅を小さくして、1つのプロセッサ(PE、あるいは、スレッド)で処理していたが、本実施形態においては、このブロック幅を大きくすると共に、L1〜L3をそれぞれ別のスレッドに割り当てて、並列に処理させる。なお、ここでは、スレッドの数を3つとしている。
【0016】
共有メモリ型スカラ並列計算機の各プロセッサは独立に1次及び2次キャッシュが備わっている。特に1次キャッシュ上のデータに載る範囲で、計算を行うことが高性能を引き出す上で重要である。
【0017】
データ量の大きな問題を多数PEで解く場合、各PEでデータを局所化して全体としてはブロックのLU分解を並列に計算し、できるだけ大きなブロック幅を確保することが必要となる。
【0018】
このために、図3に示されるように、各プロセッサでローカルに作業域を確保してL2キャッシュ(2次キャッシュ)に載る大きさで計算を行う。このとき、並列に更新するとき必要なブロック対角部分Dは、各PE(各スレッド)でコピーして各プロセッサ(各スレッド)で冗長に計算する。また、LU分解において、ピボットを取った後、行ベクトルの入れ替えを行う場合は、いずれかのプロセッサの2次キャッシュメモリ上に設けられた共用メモリ域を介して、各PEで通信して、必要な情報の共用を行うようにする。
【0019】
b:ブロック幅、k:各プロセッサが分担する列ブロックの1次元目の大きさ(ブロック幅が小行列ブロックの列方向であるので、小行列ブロックの行の数、すなわち、図2におけるL1〜L3のブロックの合計の行の数)としたとき、b×(b+k)×8〜8Mbyteを満たすものをブロック幅として採用する。
【0020】
そして、作業域(1次キャッシュメモリ)にコピーした部分に関して、キャッシュメモリ上のデータを利用したLU分解を行う。
また、処理すべき行列が占有するメモリ量が大きく(行列が大きく)、それに比べてプロセッサ数が少ないため、並列処理向けの列ブロックのブロック幅が小さくなるときは、ブロック幅を分割して必要なブロックを確保してから、行ブロックの更新と行列積の更新を並列に行う。
【0021】
更に、各プロセッサ(各スレッド)での作業域上でのLU分解は、内積法など更新部分の内積ベクトルの長さが大きな方法を、例えば、アルゴリズムの再帰的な呼び出しで、更新部分の性能を引き出しながらLU分解を行う方法を利用することで、キャッシュ上のデータを効率よく利用する。
【0022】
図4は、本実施形態のLU分解の処理の流れを示す概略フローチャートである。
まず、ステップS1において、スレッド数及び問題の大きさ(処理すべき行列の大きさ)からブロック幅を決定する。次に、ステップS2において、各スレッドが処理するブロックを決定し、各スレッドで処理するブロック(図2のD及びLi)を作業域にコピーする。そして、ステップS3において、各スレッドでピボットを決定し、その中での最大値を示すピボットを共用域を使って決定し、最大値を示すピボットを用いて、行ベクトルを入れ替え、各スレッドで、上記ブロックDとブロックLiとをLU分解する。
【0023】
ステップS4においては、処理が終わりか否かを判断し、終わりの場合には、処理を終了するが、終わりでない場合には、ステップS5において、各スレッドで並列にブロックLL(図2参照)を使って、ブロックUi(図2参照)をLU分解のアルゴリズムに従って更新する。そして、ステップS6において、各スレッドで並列に、ブロックCi(図2参照)をブロックLiとブロックU(Uiを組み合わせたもの:図2参照)の積で更新する。そして、ステップS2に戻り、次のブロック(Ciからなるブロック)を、同様の方法でLU分解する。そして、次第に小さくなる未処理ブロックが最後にブロックDに対応する部分のみになり、1つのスレッドでLU分解が完了すると、行列全体についてLU分解が終了する。
【0024】
図5〜図10は、本実施形態のLU分解の方法をより詳細に説明する図である。
ここでは、2048×2048の行列を4スレッドでLU分解する場合を例にとって説明する。
【0025】
まず、2048×2048の行列をブロック幅256でLU分解するものとする。各ブロックの区分けは、図5に示すとおりである。そして、4つのCPU(スレッド)で処理を実行する場合、各スレッドで連続な領域((256+448)×256:8MB(L2キャッシュの大きさ)より小さい領域)を確保し、図6に示すように、各スレッドの各領域にD1+L1、D1+L2、D1+L3、D1+L4をコピーする。
【0026】
なお、ブロック幅は例えば、以下のようにして決める。
問題の大きさ(行列の大きさ:行列の次数)をn、スレッドの数を#THREADとして、
【0027】
【数1】
【0028】
とおいて、
nb≧512なら、ブロック幅=512
nb≧256なら、ブロック幅=256
nb≧128なら、ブロック幅=128
それ以外、ブロック幅=64
というように、メニュー化しておき、この中から選ぶようにする。
【0029】
すなわち、LU分解のコストは、2n3 /3(nは行列の次数)で、n3 に比例する。したがって、全体のコストを#THREADで並列化し、その1%が最後に1スレッドで行うブロック幅程度となるように決める。
【0030】
ここで、理解を助けるため、並列化しない場合のLU分解のアルゴリズムを図7に示す。図7においては、LT=D1+L1+L2+L3+L4の部分をLU分解するアルゴリズムを示している。LTは、2048×256のブロックとなっている。
【0031】
まず、(1)の部分において、ピボットを決定する。iblks はLTの幅であり、今の場合、256である。また、lengは、LTの長さであり、今の場合、2048である。jjには、ピボットの存在する行番号が、TMPには、ピボットの絶対値が設定される。
【0032】
そして、(2)の部分において、現在処理しているLT内の列番号iより、ピボットの存在する行番号jjが大きい場合に、i番の行のデータをjj番の行のデータと入れ替える。次に、(3)の部分で、列iのピボットを使って、LU分解の演算を行う。
【0033】
この(1)〜(3)をiが1〜iblks にわたって繰り返し演算する。
ここで、LTの長さである、lengがもっと大きくなると、これらの処理はL2キャッシュメモリのデータを入れ替えてしまい、著しい性能低下を引き起こす。そこで、図5のように、データを分散し、L2キャッシュメモリにデータを保持したまま処理を行うようにする。各スレッドにおけるアルゴリズムは図8に示すとおりである。
【0034】
なお、図8においては、LTiは、ローカル域、pivot(4)、GPIVOT、ROW (iblks )は、共用域に格納されるデータである。
まず、(4)において、各スレッドでピボットを取る。そして、(5)で最大値を配列pivot(4)のスレッド番号の要素に格納する。(5)の後に、バリア同期を取って、(6)において、最大ピボットを持つスレッド番号をGPIVOTに格納する。そして、(6)の後で、再びバリア同期を取る。次に、(7)において、最大ピボットを持つスレッドが共用域ROWに最大ピボットの行ベクトルを格納し、バリア同期を取る。(8)においては、GPIVOTが0のときは、最大ピボットはD1の内にあるか、入れ替え不要であり、ローカル域に入れ替える。GPIVOTが#THREADに等しいとき、すなわち、GPIVOTが0より大きい場合には、最大ピボットを持たないスレッドは、ROWの内容とi行目の行ベクトルを入れ替える。そして、(9)、(10)において、LU分解のための演算を行う。そして、上記(4)〜(10)を処理するブロックの全ての列について行う。すなわち、1〜iblks までのiについて処理を繰り返す。
【0035】
ここで、LTiのLU分解の最後でバリア同期を取る。そのあと、各スレッドのD1の部分はLLとUUにLU分解されている。そして、各スレッドでU1←LL-1U1、U2←LL-1U2、U3←LL-1U3、U4←LL-1U4を各スレッドで並列に計算する。この計算の後で、D1、L1、L2、L3、L4をローカル域から行列Aにコピーバックし、バリア同期を取る。更に、C1←C1−L1×U、C2←C2−L2×U、C3←C3−L3×U、C4←C4−L4×Uを並列に書くスレッドで行い、最後にバリア同期を取る。
【0036】
図9は、上記処理によって、1段階の処理が終わった後の行列の様子を説明する図である。
図9に示されるように、上記処理をすることによって、行列の外側の行及び列が処理されたので、次に、残された左下の部分を同様の方法によって順次処理する。すなわち、ブロック幅iblksを縮小した部分を同じように分割し、図9に示すように、D1、L1、L2、L3、L4のブロックに分けて、各スレッドにコピーし、上記と同じ処理を行う。このように、処理を繰り返していくと、最後に256×256のブロックが残る。この部分は1つのスレッドでLU分解して処理を終了する。
【0037】
なお、上記処理では、LTiをLU分解するとき、キャッシュメモリ上のデータを効率よく利用するため再帰的なLU分解を利用している。また、Ci←Ci−Li×Uの演算は、キャッシュメモリ上のデータを有効に利用した方法が既知の技術として知られている。
【0038】
図10は、再帰的LU分解アルゴリズムを説明する図である。
再帰的LU分解のアルゴリズムはサブルーチンLUとして与えられる。LUの取る変数は、図9のアルゴリズムで出てきた、LTi(各スレッドでD1+Liを格納)、k(LTiの1次元目の大きさ)、iblks (ブロック幅)の他に、LU分解を始める位置を示すist、LU分解を行う幅であるnwidである。
【0039】
まず、サブルーチンの最初で、nwidが8、すなわち、LU分解を行う幅が8であるか否かを判断する。YESの場合には、LTi(ist :k、ist :ist +nwid−1)を並列にLU分解する。ここで、ist :kと言う表記は、変数がist からkまでのLTiを示す意味で、ist :ist +nwid−1というのは、変数がist からist +nwid−1までのLTiを示す意味である。以下においても同様である。
【0040】
LTiのLU分解においては、上記(4)〜(10)の処理を行う。ただし、行の入れ替え部分は、長さiblks で、LTi(i、1:iblks )を入れ替える。また、上記判断がNOの場合には、LU分解を行う幅nwidを2で割った値のLU分解をのLU分解のサブルーチンを再帰的に呼び出して行う。その後、TRSというルーチンを呼び出す。このルーチンは、LTi(ist :ist +nwid/2−1、ist +nwid/2:ist +nwid)を更新する。更に、LTi(ist :ist +nwid/2−1、ist :ist +nwid/2−1)の下三角行列LLを利用して、LL-1をCiに左からかけて更新する。次に、MMというルーチンを呼び出す。このルーチンでは、
LTi(ist +nwid/2:k、ist +nwid/2:ist +nwid)=LTi(ist +nwid/2:k、ist +nwid/2:ist +nwid)−LTi(ist +nwid/2:k、ist :ist +nwid/2−1)×LTi(ist :ist +nwid/2−1、ist +nwid/2:ist +nwid)
を演算する。
そして、その後、バリア同期を取り、LU分解のサブルーチンを再帰的に呼び出し、処理した後、処理が終わると、サブルーチンを抜ける。
b)正値対称行列の連立1次方程式の解法
図11〜図13は、正値対称行列の場合にコレスキー分解を行う処理の概念を説明する図である。
【0041】
実行列の連立1次方程式と同様に、行列をD、L1、L2、L3に分割して、各スレッドに対角行列Dと更新する列ブロック部分L1、L2、L3を作業域にコピーして(図12参照)、独立に並列にして列ブロック部分をコレスキー分解する。なお、この場合、ピボットを取る必要がない。この分解された列ブロックを利用して、小下三角行列(C1〜C6からなる)を更新する。この更新部分を、並列に行うとき負荷を均等にするために、更新するべき下三角行列をスレッド数を#Tとしたとき、2×#Tに同じブロック幅に分ける(すなわち、C1とC6、C2とC5、C3とC4を組み合わせて処理を行うようにする)。各スレッドはi番目及び2×#T+1−i番目のブロックを更新することで負荷を均等にする。
【0042】
現在考えている行列が正値対称行列であるので、図2のUに対応する部分は、L1+L2+L3からなる列ブロックの転置LT となっているので、この場合には、LT 部分は、演算を行う必要が無く、L1、L2、L3をそれぞれ転置してコピーすればよい。
【0043】
図13は、再帰的コレスキー分解の処理の進行の状況を説明する図である。
まず、(1)において、行列の一番左のブロックをコレスキー分解し、30の部分を31にコピーする。そして、31の下の点線で囲まれている部分を斜線部分を用いて更新する。次に、(2)において、処理する列の幅を2倍にして32の部分を33にコピーし、33の下の点線で囲まれた部分を斜線部分を用いて、更新する。そして、(3)に進んで、34の部分を35にコピーし、35の下の点線で示された部分を、斜線部分を用いて更新する。更に、(4)において、36で示される部分を37にコピーし、37の下の部分を斜線部分を用いて更新する。更に、(5)において、38を39にコピーし、39の下の点線で示されている部分を斜線部分を用いて更新し、(6)において、40を41にコピーし、斜線部分を用いて、41の下の部分を更新する。更に、(7)において、42の部分を43にコピーし、斜線部分を用いて43の下の部分を更新する。このように、行列の一部にコレスキー分解を再帰的に繰り返し適用しながら、最終的には行列全体をコレスキー分解する。
【0044】
図14〜図16は、変形コレスキー分解のアルゴリズムをより詳細に説明する図である。
ここでは、説明を簡略化するため4スレッドを使って処理を行う場合を説明する。
【0045】
まず、各スレッドに図14に示されるDx及びLiを連続的な領域LTiにコピーする。そして、LTiをLDLT 分解する。LDLT 分解は再帰的な方法で行う。そして、DLiへ以下の計算で値を並列にコピーする。
【0046】
DLi←D×LiT 、ここで、DはDxの対角要素であり、また、右辺は各スレッドのローカル域。
そして、Dx(スレッド1番から)他のLiは並列に、もとの領域にコピーバックする。そして、ここで、バリア同期を取り、図15に示されるように、C1とC8、C2とC7、C3とC6、C4とC5をペアとして、各スレッドで並列に更新する。すなわち、
・スレッド1では、
C1←C1−L11×DL11
C8←C8−L42×DL42
・スレッド2では、
C2←C2−L12×DL12
C7←C7−L41×DL41
・スレッド3では、
C3←C3−L21×DL21
C6←C6−L32×DL32
・スレッド4では、
C4←C4−L22×DL22
C5←C5−L31×DL31
という演算を行う。
【0047】
そして、ここまでの演算が終わった時点でバリア同期を取り、1周り小さくなった行列の領域に関して同じ処理を行う。以上を繰り返し、最後は1つのスレッドで、LDLT 分解して処理を終了する。
【0048】
上記、各スレッドで行うCiの更新処理をより一般的な言葉で述べると、スレッド数を#THREADとして、Lを2×#THREAD個に分割し、Cの下三角部分も同じく2×#THREAD個に分割する。そして、上と下から#THREAD個のペアを作り、このペアでCの分割した部分を更新するという処理になる。
【0049】
図16は、LDLT 分解の再帰的アルゴリズムを示す図である。
LDLT 分解のアルゴリズムは、サブルーチンLDLとして実現される。サブルーチンが取る変数は、前述のLU分解の場合と同様である。
【0050】
まず、nwidが8の場合には、(20)の部分で、直接LDLT 分解を行う。そして、LTi(ist +8:k、ist :ist +7)を更新する。このとき、LTi(ist :ist +7、ist :ist +7)の上三角部分にDLT が入っているので、(DLT )-1を右からかけることによって更新する。
【0051】
(20)の最初のIF文で、nwidが8でないと判断された場合には、nwid/2を新たなnwidとしてサブルーチンLDLを呼び出し、実行する。ここで、LTi(ist :ist +nwid/2−1、ist +nwid/2:ist +nwid−1)にDLT をコピーする。Dは、LTi(ist :ist +nwid/2−1、ist :ist +nwid/2−1)の対角要素であり、Lは、LTi(ist +nwid/2:ist +nwid−1、ist :ist +nwid/2−1)であり、このLを転置する。
【0052】
そして、LTi(ist +nwid/2:k、ist +nwid/2:ist +nwid−1)を更新する。すなわち、
LTi(ist +nwid/2:k、ist +nwid/2:ist +nwid−1)=LTi(ist +nwid/2:k、ist +nwid/2:ist +nwid−1)−LTi(ist +nwid/2:k、ist :ist +nwid−1)×LTi(ist :ist +nwid/2−1、ist +nwid/2:ist +nwid−1)
の演算を行う。
【0053】
次に、LDLT 分解のサブルーチンLDLを再帰的に呼び出す。そして、処理が終わったら、サブルーチンを抜ける。
なお、本発明の実施形態は、上記説明から分かるように、共有メモリ型スカラ並列計算機のアルゴリズムとして与えられるので、このアルゴリズムをプログラムとして実現することになる。あるいは、該並列計算機をLU分解専用機あるいはコレスキー分解専用機として使用する場合には、ROMなどにプログラムを書き込んでおくことも可能であるが、汎用の並列計算機として使用する場合には、本発明の実施形態のアルゴリズムは、CD−ROM等の可搬記録媒体や、ハードディスクなどの記録媒体にプログラムとして記録しておき、必要に応じて、プログラムをプロセッサにロードして使用する形態が考えられる。
【0054】
このような場合、本発明の実施形態のアルゴリズムを実現するプログラムは、可搬記録媒体などを使って、ユーザに配布が可能である。
(参考文献)
・LU分解の文献は、1)、2)、変形コレスキー分解の文献は2)
1)P.AMESTOY, M.DAYDE, and I.DUFF,“Use of computational kernels in the solution of full and sparse linear equations”, M.COSNARD, Y.ROBERT, Q.QUINTON, and M.RAYNAL, PARALLEL&DISTRIBUTED ALGORITHMS, North-Holland, 1989, pp,13-19
2)G.H.Golub, C.F.van Loan, “Matrix Computations”, second edition, The Johns Hopkins University Press, 1989
・日本語の文献では、以下のものにLU分解及びLDLT 分解の解説がある。
・「数値解析」森 正武著、共立出版会社
・「スーパーコンピュータとプログラミング」島崎眞昭著、共立出版社
(付記1)複数のプロセッサモジュールを持つ共有メモリ型スカラ並列計算機の行列演算において、
行列を小行列ブロックに分けるブロック化ステップと、
該小行列ブロックの内、対角ブロックと対角ブロックでない小行列ブロックとを該複数のプロセッサモジュールのローカルメモリ領域に格納する格納ステップと、
該複数のプロセッサモジュールが並列に、それぞれ有するブロックを演算することにより、複数のプロセッサモジュールで、対角ブロックを冗長に演算する演算ステップと、
該演算ステップで得られた小行列ブロックの演算結果を使って、該行列を更新する更新ステップと、
を備えることを特徴とする並列行列処理方法を情報装置に実現させるプログラムを格納した、情報装置読み取り可能な記録媒体。
【0055】
(付記2)該行列演算は、行列のLU分解であることを特徴とする付記1に記載の記録媒体。
(付記3)前記複数のプロセッサモジュールが有する小行列ブロックのデータからそれぞれがピボットの候補を抽出する抽出ステップと、
該ピボットの候補から該複数のプロセッサモジュールに共通のメモリ領域において、データ値が最大値を示すピボット候補を最終的なピボットと決定するピボット決定ステップと、
を更に備え、該決定されたピボットを用いてLU分解を行うことを特徴とする付記2に記載の記録媒体。
【0056】
(付記4)前記LU分解は、前記行列の外側から再帰的なアルゴリズムによって順次行列の更新を行い、前記行列の内、最後に更新し残った部分を、1つのプロセッサモジュールでLU分解することにより、該行列全体についてLU分解を完了することを特徴とする付記2に記載の記録媒体。
【0057】
(付記5)該行列演算は、行列のコレスキー分解あるいは変形コレスキー分解であることを特徴とする付記1に記載の記録媒体。
(付記6)前記コレスキー分解あるいは変形コレスキー分解は、前記行列の外側から再帰的なアルゴリズムによって順次行列の更新を行い、前記行列の内、最後に更新し残った部分を、1つのプロセッサモジュールでLU分解することにより、該行列全体についてLU分解を完了することを特徴とする付記2に記載の記録媒体。
【0058】
(付記7)前記更新ステップにおいて、更新すべき小行列ブロックの三角行列部分を前記複数のプロセッサモジュールの数の2倍の数のブロックに分割し、該分割された三角行列部分のブロックを2つずつ組み合わせて、各プロセッサモジュールのローカルメモリ域に格納し、演算をプロセッサモジュールに行わせることを特徴とする付記5に記載の記録媒体。
【0059】
(付記8)複数のプロセッサモジュールを持つ共有メモリ型スカラ並列計算機の行列演算において、
行列を小行列ブロックに分けるブロック化ステップと、
該小行列ブロックの内、対角ブロックと対角ブロックでない小行列ブロックとを該複数のプロセッサモジュールのローカルメモリ領域に格納する格納ステップと、
該複数のプロセッサモジュールが並列に、それぞれ有するブロックを演算することにより、複数のプロセッサモジュールで、対角ブロックを冗長に演算する演算ステップと、
該演算ステップで得られた小行列ブロックの演算結果を使って、該行列を更新する更新ステップと、
を備えることを特徴とする並列行列処理方法。
【0060】
(付記9)複数のプロセッサモジュールを持つ共有メモリ型スカラ並列計算機において、
行列を小行列ブロックに分けるブロック化手段と、
該小行列ブロックの内、対角ブロックと対角ブロックでない小行列ブロックとを該複数のプロセッサモジュールのローカルメモリ領域に格納する格納手段と、該複数のプロセッサモジュールが並列に、それぞれ有するブロックを演算することにより、複数のプロセッサモジュールで、対角ブロックを冗長に演算する演算手段と、
該演算ステップで得られた小行列ブロックの演算結果を使って、該行列を更新する更新手段と、
を備えることを特徴とする並列行列処理装置。
【0061】
【発明の効果】
本発明によれば、高性能かつスケーラビリティのある行列の処理方法が得られる。
【図面の簡単な説明】
【図1】共有メモリ型スカラ並列計算機のハードウェア構成例を示す図である。
【図2】本発明の実施形態に従ったLU分解の並列処理の概念を説明する図(その1)である。
【図3】本発明の実施形態に従ったLU分解の並列処理の概念を説明する図(その2)である。
【図4】本実施形態のLU分解の処理の流れを示す概略フローチャートである。
【図5】本実施形態のLU分解の方法をより詳細に説明する図(その1)である。
【図6】本実施形態のLU分解の方法をより詳細に説明する図(その2)である。
【図7】本実施形態のLU分解の方法をより詳細に説明する図(その3)である。
【図8】本実施形態のLU分解の方法をより詳細に説明する図(その4)である。
【図9】本実施形態のLU分解の方法をより詳細に説明する図(その5)である。
【図10】本実施形態のLU分解の方法をより詳細に説明する図(その6)である。
【図11】正値対称行列の場合にコレスキー分解を行う処理の概念を説明する図(その1)である。
【図12】正値対称行列の場合にコレスキー分解を行う処理の概念を説明する図(その2)である。
【図13】正値対称行列の場合にコレスキー分解を行う処理の概念を説明する図(その3)である。
【図14】変形コレスキー分解のアルゴリズムをより詳細に説明する図(その1)である。
【図15】変形コレスキー分解のアルゴリズムをより詳細に説明する図(その2)である。
【図16】変形コレスキー分解のアルゴリズムをより詳細に説明する図(その3)である。
【符号の説明】
10−1〜10−n プロセッサ
11−1〜11−n メモリモジュール
12 相互結合網
13−1〜13−n 2次キャッシュメモリ
Claims (6)
- 複数のプロセッサモジュールと、各プロセッサモジュールに対応して設けられた2次キャッシュと、各プロセッサモジュールに内蔵された1次キャッシュと、各プロセッサモジュールと該2次キャッシュを介して接続された相互結合網と、該相互結合網を介して各プロセッサモジュールがアクセス可能な複数のメモリモジュールとを持つ共有メモリ型スカラ並列計算機のLU分解を行う行列演算において、
行列を対角部分に設定される対角ブロックと、対角ブロックに隣接する行方向に長いブロックと、対角ブロックに隣接する列方向に長いブロックと、その他の正方ブロックからなる小行列ブロックに分けるブロック化ステップと、
該小行列ブロックの内、対角ブロックと対角ブロックの下に位置する列方向に長いブロックを均等に行方向(1次元目を)分割した小行列ブロックとを該複数のプロセッサモジュールのローカルメモリ領域に格納する格納ステップと、
該複数のプロセッサモジュールが並列に、該ローカルメモリ領域にそれぞれ有するブロックを対角ブロックと列方向に長い、あるいは、行方向に長いブロックとをいっしょに、正方行列をLU分解したときの対応するブロック部分の演算結果と同じになるように演算する演算ステップと、
該演算ステップで得られた小行列ブロックの演算結果を使って、その他の正方ブロックから該演算された行方向に長いブロックと列方向に長いブロックの積を減算することにより、該正方ブロックを更新する更新ステップと、
を備えることを特徴とする並列行列処理方法を情報装置に実現させるプログラムを格納した、情報装置読み取り可能な記録媒体。 - 前記複数のプロセッサモジュールが有する小行列ブロックのデータからそれぞれが、該小行列ブロックのデータの中で最も大きな値を持つ行列要素であるピボットの候補を抽出する抽出ステップと、
該ピボットの候補から該複数のプロセッサモジュールに共通のメモリ領域において、データ値が最大値を示すピボット候補を最終的なピボットと決定するピボット決定ステップと、
を更に備え、該決定されたピボットを用いてLU分解を行うことを特徴とする請求項1に記載の記録媒体。 - 前記LU分解は、前記行列の外側から再帰的なアルゴリズムによって順次行列の更新を行い、前記行列の内、最後に更新し残った部分を、1つのプロセッサモジュールでLU分解することにより、該行列全体についてLU分解を完了することを特徴とする請求項1に記載の記録媒体。
- 該行列演算は、行列のコレスキー分解あるいは変形コレスキー分解であることを特徴とする請求項1に記載の記録媒体。
- 前記コレスキー分解あるいは変形コレスキー分解は、前記行列の外側から再帰的なアルゴリズムによって順次行列の更新を行い、前記行列の内、最後に更新し残った部分を、1つのプロセッサモジュールでLU分解することにより、該行列全体についてLU分解を完了することを特徴とする請求項1に記載の記録媒体。
- 前記更新ステップにおいて、更新すべき小行列ブロックの三角行列部分を前記複数のプロセッサモジュールの数の2倍の数のブロックに分割し、該分割された三角行列部分のブロックを2つずつ組み合わせて、各プロセッサモジュールのローカルメモリ域に格納し、演算をプロセッサモジュールに行わせることを特徴とする請求項4に記載の記録媒体。
Priority Applications (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2000358232A JP3639206B2 (ja) | 2000-11-24 | 2000-11-24 | 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 |
US09/811,484 US6907513B2 (en) | 2000-11-24 | 2001-03-20 | Matrix processing method of shared-memory scalar parallel-processing computer and recording medium |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2000358232A JP3639206B2 (ja) | 2000-11-24 | 2000-11-24 | 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2002163246A JP2002163246A (ja) | 2002-06-07 |
JP3639206B2 true JP3639206B2 (ja) | 2005-04-20 |
Family
ID=18830173
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2000358232A Expired - Fee Related JP3639206B2 (ja) | 2000-11-24 | 2000-11-24 | 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 |
Country Status (2)
Country | Link |
---|---|
US (1) | US6907513B2 (ja) |
JP (1) | JP3639206B2 (ja) |
Families Citing this family (30)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3809062B2 (ja) * | 2000-11-21 | 2006-08-16 | 富士通株式会社 | マルチレベル不完全ブロック分解による前処理を行う処理装置 |
ATE283580T1 (de) * | 2001-07-05 | 2004-12-15 | Mitsubishi Electric Inf Tech | Mehrbenutzerdetektion in einem mc-cdma telekommunikationssystem |
EP1300977A1 (en) * | 2001-10-04 | 2003-04-09 | Mitsubishi Electric Information Technology Centre Europe B.V. | Parallel interference cancellation in an MC-CDMA telecommunication system |
US20030182518A1 (en) | 2002-03-22 | 2003-09-25 | Fujitsu Limited | Parallel processing method for inverse matrix for shared memory type scalar parallel computer |
EP1376380A1 (en) * | 2002-06-14 | 2004-01-02 | EADS Deutschland GmbH | Procedure for computing the Choleski decomposition in a parallel multiprocessor system |
JP3983193B2 (ja) * | 2003-03-31 | 2007-09-26 | 富士通株式会社 | 行列処理方法及び装置 |
US7318071B2 (en) * | 2003-05-27 | 2008-01-08 | Emc Corporation | System and method for transfering data from a source machine to a target machine |
US7363200B2 (en) * | 2004-02-05 | 2008-04-22 | Honeywell International Inc. | Apparatus and method for isolating noise effects in a signal |
US7574333B2 (en) * | 2004-02-05 | 2009-08-11 | Honeywell International Inc. | Apparatus and method for modeling relationships between signals |
JP2006085619A (ja) | 2004-09-17 | 2006-03-30 | Fujitsu Ltd | 帯係数行列を持つ連立1次方程式の解法プログラム |
US7937709B2 (en) | 2004-12-29 | 2011-05-03 | Intel Corporation | Synchronizing multiple threads efficiently |
US20060265445A1 (en) * | 2005-05-20 | 2006-11-23 | International Business Machines Corporation | Method and structure for improving processing efficiency in parallel processing machines for rectangular and triangular matrix routines |
US7849126B1 (en) * | 2006-03-06 | 2010-12-07 | Intellectual Property Systems, LLC | System and method for fast matrix factorization |
US8074210B1 (en) * | 2006-06-29 | 2011-12-06 | Xilinx, Inc. | Method and apparatus for producing optimized matrix triangulation routines |
US7813948B2 (en) * | 2006-08-25 | 2010-10-12 | Sas Institute Inc. | Computer-implemented systems and methods for reducing cost flow models |
US7882061B1 (en) | 2006-12-21 | 2011-02-01 | Emc Corporation | Multi-thread replication across a network |
JP4823928B2 (ja) * | 2007-01-22 | 2011-11-24 | 三菱電機株式会社 | 連立一次方程式の並列求解装置 |
US20080208553A1 (en) * | 2007-02-27 | 2008-08-28 | Fastrack Design, Inc. | Parallel circuit simulation techniques |
WO2008136045A1 (ja) * | 2007-04-19 | 2008-11-13 | Fujitsu Limited | 共有メモリ型スカラ並列計算機向け、実対称行列の三重対角化の並列処理方法 |
US8024241B2 (en) * | 2007-07-13 | 2011-09-20 | Sas Institute Inc. | Computer-implemented systems and methods for cost flow analysis |
US8380778B1 (en) * | 2007-10-25 | 2013-02-19 | Nvidia Corporation | System, method, and computer program product for assigning elements of a matrix to processing threads with increased contiguousness |
US8200518B2 (en) | 2008-02-25 | 2012-06-12 | Sas Institute Inc. | Computer-implemented systems and methods for partial contribution computation in ABC/M models |
US20090254319A1 (en) * | 2008-04-03 | 2009-10-08 | Siemens Aktiengesellschaft | Method and system for numerical simulation of a multiple-equation system of equations on a multi-processor core system |
US8417755B1 (en) | 2008-05-28 | 2013-04-09 | Michael F. Zimmer | Systems and methods for reducing memory traffic and power consumption in a processing environment by solving a system of linear equations |
US20120304042A1 (en) * | 2011-05-28 | 2012-11-29 | Jose Bento Ayres Pereira | Parallel automated document composition |
KR101585980B1 (ko) * | 2014-04-11 | 2016-01-19 | 전자부품연구원 | 멀티-프로세서의 공유 메모리를 적극 활용한 cr 알고리즘 처리 방법 및 이를 적용한 프로세서 |
US9805001B2 (en) | 2016-02-05 | 2017-10-31 | Google Inc. | Matrix processing apparatus |
US10915075B2 (en) * | 2017-10-02 | 2021-02-09 | University Of Dayton | Reconfigurable hardware-accelerated model predictive controller |
GB2598996A (en) * | 2020-04-09 | 2022-03-23 | Nvidia Corp | Fifth Generation (5G) New Radio Channel Equalization |
US20210320825A1 (en) * | 2020-04-09 | 2021-10-14 | Nvidia Corporation | Fifth generation (5g) new radio channel equalization |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5157778A (en) * | 1986-08-20 | 1992-10-20 | Digital Equipment Corporation | Method and apparatus for circuit simulation using parallel processors including memory arrangements and matrix decomposition synchronization |
EP0523544B1 (en) * | 1991-07-12 | 2002-02-27 | Matsushita Electric Industrial Co., Ltd. | Apparatus to solve a system of linear equations |
JP2956800B2 (ja) * | 1991-09-19 | 1999-10-04 | 株式会社日立製作所 | 連立一次方程式に関する計算装置 |
US5644517A (en) * | 1992-10-22 | 1997-07-01 | International Business Machines Corporation | Method for performing matrix transposition on a mesh multiprocessor architecture having multiple processor with concurrent execution of the multiple processors |
JPH08227405A (ja) | 1995-02-21 | 1996-09-03 | Hitachi Ltd | 反復処理の並列実行方法 |
US6182270B1 (en) * | 1996-12-04 | 2001-01-30 | Lucent Technologies Inc. | Low-displacement rank preconditioners for simplified non-linear analysis of circuits and other devices |
JP2959525B2 (ja) * | 1997-06-02 | 1999-10-06 | 日本電気株式会社 | データ処理装置および方法、情報記憶媒体 |
TW388921B (en) * | 1997-11-28 | 2000-05-01 | Nippon Electric Co | Semiconductor process device simulation method and storage medium storing simulation program |
-
2000
- 2000-11-24 JP JP2000358232A patent/JP3639206B2/ja not_active Expired - Fee Related
-
2001
- 2001-03-20 US US09/811,484 patent/US6907513B2/en not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
JP2002163246A (ja) | 2002-06-07 |
US20020091909A1 (en) | 2002-07-11 |
US6907513B2 (en) | 2005-06-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP3639206B2 (ja) | 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 | |
JP6977239B2 (ja) | 行列乗算器 | |
Pavel et al. | Algorithms for efficient computation of convolution | |
US7814297B2 (en) | Algebraic single instruction multiple data processing | |
Kiss et al. | Parallel realization of the element-by-element FEM technique by CUDA | |
Nehab et al. | GPU-efficient recursive filtering and summed-area tables | |
US4821224A (en) | Method and apparatus for processing multi-dimensional data to obtain a Fourier transform | |
JPH07271760A (ja) | メモリ分散型並列計算機による連立1次方程式計算処理方法および計算機 | |
Membarth et al. | Towards domain-specific computing for stencil codes in HPC | |
Shrivastava et al. | A survey of hardware architectures for generative adversarial networks | |
JP3639207B2 (ja) | 共有メモリ型スカラ並列計算機における多次元フーリエ変換の並列処理方法 | |
Nehab et al. | Parallel recursive filtering of infinite input extensions | |
Roberts et al. | Multithreaded implicitly dealiased convolutions | |
Espinoza-Valverde et al. | Coarsest-level improvements in multigrid for lattice QCD on large-scale computers | |
Kanakia et al. | Espresso-gpu: blazingly fast two-level logic minimization | |
Nobile et al. | Efficient implementation of multidimensional fast Fourier transforms on a cray X-MP | |
JP2021005242A (ja) | 情報処理装置、情報処理プログラム、及び情報処理方法 | |
Vishwanath | Time-frequency distributions: Complexity, algorithms and architectures | |
Graham et al. | Parallel algorithms and architectures for optimal state estimation | |
López et al. | Unified framework for the parallelization of divide and conquer based tridiagonal systems | |
Botros et al. | Hardware realization of Krawtchouk transform using VHDL modeling and FPGAs | |
US20040078412A1 (en) | Parallel processing method of an eigenvalue problem for a shared-memory type scalar parallel computer | |
JP2862388B2 (ja) | 超高速画像処理システムのフィルタリング処理方式 | |
WO2003034269A1 (en) | Method of performing a fft transform on parallel processors | |
Pic et al. | Wavelet transform on parallel SIMD architectures |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20040330 |
|
A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040527 |
|
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: 20050111 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20050113 |
|
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: 20080121 Year of fee payment: 3 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090121 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100121 Year of fee payment: 5 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110121 Year of fee payment: 6 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110121 Year of fee payment: 6 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120121 Year of fee payment: 7 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130121 Year of fee payment: 8 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130121 Year of fee payment: 8 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140121 Year of fee payment: 9 |
|
LAPS | Cancellation because of no payment of annual fees |