JP2019109626A - 疎行列ベクトル積演算装置及び疎行列ベクトル積演算方法 - Google Patents

Abstract

【課題】共役勾配法（ＣＧ法）を高速化する。【解決手段】疎行列ベクトル積演算装置は、共役勾配法の制御を行う共役勾配法制御装置から入力された演算データを記憶する。演算データは、疎行列との積を計算可能なベクトルと、列指向形式で圧縮された疎行列とを含む。疎行列ベクトル積演算装置は、演算データ記憶部からベクトル及び疎行列を読み込む読込部を有する。疎行列ベクトル積演算装置は、読込部によって読み込まれたベクトル及び疎行列を用いて疎行列ベクトル積の演算を実行する演算部を有する。疎行列ベクトル積演算装置は、演算部による演算の結果を演算データ記憶部に書き込む書込部を有する。【選択図】図８

Description

本発明は、疎行列ベクトル積演算装置及び疎行列ベクトル積演算方法に関する。

製品の設計に関するシミュレーションを、シミュレーションツールを使って計算機上で行う技術としてＣＡＥ（Computer Aided Engineering）が知られている。例えば、ＣＡＥは、構造解析や流体解析を目的として行われる場合がある。

ＣＡＥにおけるシミュレーションで用いられる計算手法の１つに、共役勾配法（ＣＧ法：Conjugate Gradient Method）がある。また、ＣＧ法の実行時間の大部分は、疎行列ベクトル積の演算（ＳｐＭＶ：Sparse matrix-vector multiplication）によって占められている。

また、ＳｐＭＶを高速化する技術として、メモリを共有する複数の処理装置が、ＣＳＲ（Compressed Sparse Row format）形式からＪＡＤ（Jagged Diagonal format）形式に変換された疎行列のデータを用いて、ＳｐＭＶを行う技術が知られている。

特開２００１−２０９６３１号公報 特開平８−２１２１８６号公報

ところで、ＣＧ法におけるＳｐＭＶの計算対象の疎行列のサイズは、非常に大きい場合がある。また、ＣＧ法ではＳｐＭＶが反復して行われるため、ＳｐＭＶの反復のたびにＣＰＵとメモリとの間で疎行列の転送が行われる場合がある。しかし、上記技術では、ＣＰＵとメモリとの間で疎行列の転送時間を短縮し、ＣＧ法を高速化することが困難な場合がある。

一つの側面では、共役勾配法を高速化できる疎行列ベクトル積演算装置及び疎行列ベクトル積演算方法を提供することを目的とする。

一つの態様において、疎行列ベクトル積演算装置は、共役勾配法の制御を行う共役勾配法制御装置から入力された演算データを記憶する。演算データは、疎行列との積を計算可能なベクトルと、列指向形式で圧縮された疎行列とを含む。疎行列ベクトル積演算装置は、演算データ記憶部からベクトル及び疎行列を読み込む読込部を有する。疎行列ベクトル積演算装置は、読込部によって読み込まれたベクトル及び疎行列を用いて疎行列ベクトル積の演算を実行する演算部を有する。疎行列ベクトル積演算装置は、演算部による演算の結果を演算データ記憶部に書き込む書込部を有する。

一つの態様によれば、共役勾配法を高速化できる。

図１は、ＳｐＭＶの一例を示す図である。 図２は、ＣＳＲ形式の一例を示す図である。 図３は、ＳｐＭＶの処理の一例を示す図である。 図４は、実施例１におけるシミュレーション装置の一例を示す図である。 図５は、ＥＬＬ形式及びＪＡＤ形式の一例を示す図である。 図６は、実施例１における演算データの一例を示す図である。 図７は、実施例１における前処理部の一例を示す図である。 図８は、実施例１における加速部の一例を示す図である。 図９は、実施例１におけるＳｐＭＶの処理の一例を示す図である。 図１０は、実施例１における演算部の一例を示す図である。 図１１は、実施例１におけるＣＧ法による処理の一例を示すフローチャートである。 図１２は、実施例１におけるＳｐＭＶの処理の一例を示すフローチャートである。 図１３は、実施例１における演算処理の一例を示すフローチャートである。 図１４は、ハードウェア構成例を示す図である。

以下に、本願の開示する疎行列ベクトル積演算装置及び疎行列ベクトル積演算方法の実施例を図面に基づいて詳細に説明する。なお、この実施例によりこの発明が限定されるものではない。また、以下に示す各実施例は、矛盾を起こさない範囲で適宜組み合わせても良い。

実施例１における、後に説明するシミュレーション装置１０は、ＣＧ法を使ったシミュレーションを実行する。ＣＧ法は、ＳｐＭＶの反復演算により、一次方程式の解を求める手法である。

なお、以下において、行列、ベクトル及び配列の各インデックスは０から始まるものとする。また、行列Ｍの行インデックスがｉ、列インデックスがｊである要素をＭ_ｉｊと表記する場合がある。また、ベクトルｖのインデックスがｉである要素をｖ_ｉと表記する場合がある。また、配列array［］のインデックスがｉである要素をarray［ｉ］と表記する場合がある。

また、行列の行方向は、行インデックスが増加していく方向であり、右方向と同意であることとする。行列の列方向は、列インデックスが増加していく方向であり、下方向と同意であることとする。

ここで、例として、シミュレーション装置１０がＣＧ法を用いてｎ元の一次方程式の解を求める場合の処理の流れを説明する。このとき、シミュレーション装置１０は、ｎ×ｎの行列Ａ及び要素数ｎのベクトルｙから、行列Ａとの積がベクトルｙとなるようなベクトルｘを求める。なお、以下において、行列Ａのような行列を係数行列と表記し、ベクトルｘのようなベクトルを解ベクトルと表記し、ベクトルｙのようなベクトルを定数ベクトルと表記する場合がある。また、係数行列及び定数ベクトルは既知であるが、解ベクトルの真値は未知である。

まず、シミュレーション装置１０は、ベクトルｘに任意の初期値を設定し、行列Ａとベクトルｘの積の計算、すなわちＳｐＭＶを実行する。そして、シミュレーション装置１０は、行列ＡとベクトルｘのＳｐＭＶの結果として得られるベクトルｙ´とベクトルｙとの残差から計算した勾配を基に、ベクトルｙ´がベクトルｙに近付くようにベクトルｘを更新する。

さらに、シミュレーション装置１０は、行列Ａと更新したベクトルｘとのＳｐＭＶを実行する。その後、シミュレーション装置１０は、ＳｐＭＶ及びベクトルｘの更新を所定の条件が満たされるまで繰り返し、そのときのベクトルｘを一次方程式の解又は近似解として得ることができる。

ＣＧ法を使って製品の設計に関するシミュレーションを行う場合、係数行列は疎行列である場合がある。疎行列とは、行列の要素の多くが０の行列である。ここで、図１を用いて、疎行列のＳｐＭＶについて説明する。図１は、ＳｐＭＶの一例を示す図である。図１に示すように、行列Ａは疎行列である。このとき、ベクトルｙの要素ｙ_０は（１）式で計算される。

ｙ_０＝ａ_００×ｘ_０＋０×ｘ_１＋０×ｘ_２＋ａ_０３×ｘ_３＋０×ｘ_４ ・・・（１）

ここで、（１）式の第２項、第３項及び第５項は必ず０になる。このため、プロセッサは（１）式の計算を行う場合、ベクトルｘの要素ｘ_１、ｘ_２及びｘ_４を参照する必要はない。

しかしながら、ＣＰＵやＧＰＵ（Graphics Processing Unit）等のプロセッサはメモリからのバースト転送によってデータを取得する。このため、例えばＣＰＵが１回の転送サイクルでメモリから取得可能なデータは、連続する一定量のデータに限られる。

ここで、１回の転送サイクルで取得可能なデータがベクトルｘの連続する３要素分であるＣＰＵが（１）式の計算を行う場合を考える。また、このとき、ベクトルｘの各要素はインデックスの順に連続してメモリに格納されているものとする。

この場合、ＣＰＵは、まず、ａ_００×ｘ_０を計算するために、１回目の転送サイクルでｘ_０を取得する。そして、ＣＰＵは、ａ_０３×ｘ_３を計算するために、２回目の転送サイクルでｘ_３を取得する。このように、ＣＰＵは、（１）式を計算するために、少なくとも２回の転送サイクルを要する。

仮に、ＣＰＵがベクトルｘから取得する要素がｘ_０、ｘ_１のように連続したものであれば、ＣＰＵは１回の転送サイクルで必要な要素を取得することができることになる。逆に、取得するベクトルの要素に連続性がない場合、すなわちベクトルへのアクセスがランダムな場合は、必要な転送サイクルが増加し、その結果ＳｐＭＶの演算速度が低下する。

さらに、シミュレーションで用いられる一次方程式の元の数は数百万から１千万を超える大きな数になる場合がある。そのような場合、係数行列や解ベクトルのサイズは非常に大きくなり、取得するベクトルの要素に連続性がないことによるＳｐＭＶの演算速度への影響も大きくなる。

また、疎行列は、圧縮した形式で表される場合がある。疎行列の圧縮形式の１つにＣＳＲ形式がある。図２は、ＣＳＲ形式の一例を示す図である。図２に示すようにＣＳＲ形式の疎行列は、row_off、col、valの３つの配列で表される。

配列valには、疎行列の０でない要素（非零要素）が、疎行列における行インデックスが小さい順に格納される。また、行インデックスが同じ非零要素については、列インデックスが小さいものが配列valにおいて先の順番になる。

また、配列colには、配列valに格納された各非零要素の対応する位置に、各非零要素の列インデックスが格納される。また、配列row_offには、配列colのインデックスであって、疎行列の各行の非零要素のうち最も列インデックスが小さい非零要素の列インデックスが格納されたインデックスが、当該非零要素の行インデックスが小さい順に格納される。

ここで、例えば、図１の行列Ａの要素ａ_００の行インデックス及び列インデックスは、いずれも０である。また、配列valにおける要素ａ_００のインデックスは０である。このため配列colのインデックス０の要素は０である。また、配列colのインデックスであって、疎行列の１行目（行インデックスが０の行）の非零要素のうち最も列インデックスが小さい非零要素であるａ_００の列インデックスが格納されたインデックスは０である。このため、配列row_offのインデックス０の要素は０である。

また、例えば、図１の行列Ａの要素ａ_３０の行インデックス及び列インデックスは、それぞれ３及び０である。また、配列valにおける要素ａ_３０のインデックスは６である。このため配列colのインデックス６の要素は０である。また、配列colのインデックスであって、疎行列の４行目（行インデックスが３の行）の非零要素のうち最も列インデックスが小さい非零要素であるａ_３０の列インデックスが格納されたインデックスは６である。このため、配列row_offのインデックス３の要素は６である。

図３の疑似コードを用いて、図１の行列ＡがＣＳＲ形式で表されている場合の、ＣＰＵによるＳｐＭＶの処理について説明する。図３は、ＳｐＭＶの処理の一例を示す図である。図３のｉは、ベクトルｙの要素のインデックスに対応している。また、ｊは、ベクトルｘのインデックスに対応している。また、ｘ［］及びｙ［］は、それぞれベクトルｘ及びベクトルｙを格納した配列である。

例えば、図２に示すように、row_off［０］は０であり、row_off［１］は２である。また、col［０］は０であり、col［１］は３である。このため、図３の３行目及び４行目に示すように、ｉ＝０の場合、ＣＰＵは、ｘ［０］及びｘ［３］を取得し、取得したｘ［］の各要素と対応するval［］の要素との積をｙ［０］に足していく。つまり、ｉ＝０の場合、ＣＰＵは、（１）式の計算を行うことになる。この場合、取得するｘ［］の要素が連続でないことの影響を受けて、ＳｐＭＶの演算速度が低下する。

実施例１のＳｐＭＶ装置１０ｅは、ＣＧ法によるシミュレーションを行う際の疎行列の転送を最小限に抑えつつ、ベクトルへのアクセスが連続でないことによるＳｐＭＶの演算速度の低下を抑制することを目的とするものである。

［機能ブロック］
次に、本実施例における疎行列ベクトル積演算装置１０ｂを有するシミュレーション装置１０の一例について、図４を用いて説明する。図４は、実施例１におけるシミュレーション装置の一例を示す図である。図４に示すように、本実施例におけるシミュレーション装置１０は、共役勾配法制御装置１０ａと、疎行列ベクトル積演算装置１０ｂとを有する。共役勾配法制御装置１０ａは、例えばＣＰＵによって実現される。一方、疎行列ベクトル積演算装置１０ｂは、例えばＦＰＧＡ（Field Programmable Gate Array）によって実現される。

共役勾配法制御装置１０ａは、変換部１１及び更新部１２を有する。また、疎行列ベクトル積演算装置１０ｂは、演算データ記憶部１３、前処理制御部１４、前処理部１５及び加速部１６を有する。

共役勾配法制御装置１０ａは、ＣＧ法の制御を行う。具体的には、共役勾配法制御装置１０ａは、ユーザから入力された設定等を基にＣＧ法の実行し、その際、疎行列ベクトル積演算装置１０ｂにＳｐＭＶを実行させる。なお、ＳｐＭＶにおける係数行列に対する前処理は、疎行列ベクトル積演算装置１０ｂによって行われる。

変換部１１は、ＳｐＭＶで用いられる疎行列である係数行列を列指向形式で圧縮する。例えば、変換部１１は、無圧縮又はＣＳＲ等の行指向形式で表された係数行列を、列指向形式であるＣＳＣ（Compressed Sparse Column format）形式、ＥＬＬ（Ellpack-Itpack generalized diagonal format）形式、ＪＡＤ形式等に変換する。

図５は、ＥＬＬ形式及びＪＡＤ形式の一例を示す図である。図５に示すように、変換部１１は、行列ＡをＥＬＬ形式に変換する場合、行列Ａの非零要素を左側に詰めた行列を仮定し、仮定した行列の列単位で要素を取得し、配列に格納していく。また、行列ＡをＪＡＤ形式に変換する場合、変換部１１は、仮定した行列の行を要素数の順でソートした上で、ＥＬＬ形式と同様に要素を配列に格納していく。

本実施例、変換部１１は、行列ＡをＪＡＤ形式に変換し、変換した係数行列を、演算データとして疎行列ベクトル積演算装置１０ｂの演算データ記憶部１３に格納するものとする。また、このとき、変換部１１は、係数行列とともに、初期値を設定した解ベクトルを演算データ記憶部１３に格納し、疎行列ベクトル積演算装置１０ｂにＳｐＭＶを実行させる。

図６は、実施例１における演算データの一例を示す図である。図６に示すようにＪＡＤ形式の係数行列である演算データは、colind、value、jdptr及びpermの４つの配列で表される。

ここで、変換部１１によるＪＡＤ形式の各配列の作成方法の一例を説明する。まず、変換部１１は、行列Ａの行インデックスを、各行の要素数順に配列permに格納する。次に、変換部１１は、配列permに格納した行インデックスの順序に従って行列Ａの行をフォーカスしていき、各行の配列valueに格納されていない要素のうち、最も左側にある要素を配列valueに格納する。変換部１１は、フォーカスしている行の配列valueに要素を格納した場合、又はフォーカスしている行に格納する要素が存在しない場合は次の行に進む。

さらに、変換部１１は、配列valueに格納した要素の列インデックスを配列colindの対応する位置に格納する。また、変換部１１は、配列permに格納した行インデックスの順序が最初である行の要素の配列valueにおける格納先のインデックスを配列jdptrに格納していく。

更新部１２は、ＣＧ法の手順に従い、疎行列ベクトル積演算装置１０ｂによるＳｐＭＶの結果を基に解ベクトルを更新する。ここで、疎行列ベクトル積演算装置１０ｂは、ＳｐＭＶの結果として得られる定数ベクトルの計算値を演算データ記憶部１３に格納するものとする。また、更新部１２は、更新した解ベクトルをさらに演算データ記憶部１３に格納し、疎行列ベクトル積演算装置１０ｂにＳｐＭＶを実行させることができる。

演算データ記憶部１３は、共役勾配法の制御を行う共役勾配法制御装置１０ａから入力された演算データであって、係数行列との積を計算可能な解ベクトルと、列指向形式で圧縮された係数行列とを含む演算データを記憶する。なお、係数行列は疎行列である。例えば、演算データ記憶部１３は、ＥＬＬ形式又はＪＡＤ形式で圧縮された係数行列を記憶することができる。ただし、本実施例では、演算データ記憶部１３は、ＪＡＤ形式で圧縮された係数行列を記憶するものとする。

前処理部１５は、予め設定された演算データに対する複数の前処理のそれぞれに対応した処理部を備え、指定された前処理に対応する処理部を用いて前処理を実行する。前処理制御部１４は、共役勾配法制御装置１０ａによって指定された前処理を、前処理部１５に対して指定する。なお、共役勾配法制御装置１０ａへの前処理の指定は、ユーザによって行われてもよい。

図７を用いて、前処理部１５について説明する。図７は、実施例１における前処理部の一例を示す図である。図７に示すように、前処理部１５は、対角スケーリング処理部１５１と、ＳＳＯＲ処理部１５２と、ＩＬＵ処理部１５３と、ＡＭＧ処理部１５４とを有する。

対角スケーリング処理部１５１は、対角スケーリングによる前処理を行う。ＳＳＯＲ処理部１５２は、ＳＳＯＲ（Symmetric Successive Over-relaxation）による前処理を行う。ＩＬＵ処理部１５３は、ＩＬＵ（Incomplete LU Factorization）（０）、ＩＬＵ（１）等による前処理を行う。ＡＭＧ処理部１５４は、ＡＭＧ（Algebraic Multigrid）による前処理を行う。なお、前処理部１５に備えられる処理部は、図７に示すものに限られない。前処理部１５は、行列に対する前処理を行う処理部を複数有していればよい。

前処理部１５に含まれる各処理部は、回路の接続先の変更等により、再構成可能に配置されている。例えば、前処理制御部１４によって指定された前処理が対角スケーリングであった場合、前処理部１５は、演算データ記憶部１３と対角スケーリング処理部１５１との接続を有効にする。また、そのとき、前処理部１５は、演算データ記憶部１３と、ＳＳＯＲ処理部１５２、ＩＬＵ処理部１５３及びＡＭＧ処理部１５４との接続を無効にする。

加速部１６は、ＳｐＭＶを実行する。また、加速部１６によるＳｐＭＶの実行結果は、演算データ記憶部１３に記憶される。なお、演算の結果は、例えば、定数ベクトルの計算値である。また、加速部１６は、ＣＧ法を高速化するように構成されている。

ここで、前処理部１５は、加速部１６によるＳｐＭＶが完了するたびに前処理を実行することができる。また、加速部１６は、前処理部１５による前処理が完了するたびにＳｐＭＶを実行する。これにより、疎行列ベクトル積演算装置１０ｂは、前処理及びＳｐＭＶを反復して行うことができる。

その際、共役勾配法制御装置１０ａは、所定のタイミングで疎行列ベクトル積演算装置１０ｂによる反復処理に割り込みを行い、定数ベクトルの計算値を取得する。共役勾配法制御装置１０ａは、取得した計算値と真値との残差から計算した勾配を基に、演算データ記憶部１３に記憶された解ベクトルを更新する。

図８を用いて、加速部１６について説明する。図８は、実施例１における加速部の一例を示す図である。図８に示すように、加速部１６は、ロード部１６１、割当部１６２及び複数の演算部１６３を有する。また、加速部１６は、演算データ記憶部１３と接続されている。また、ロード部１６１は、読込部１６１ａ及び書込部１６１ｂを有する。

読込部１６１ａは、演算データ記憶部１３から解ベクトル及び係数行列を読み込む。割当部１６２は、演算部１６３のそれぞれに、演算に必要なデータを割り当てる。演算部１６３は、読込部１６１ａによって読み込まれた解ベクトル及び係数行列を用いてＳｐＭＶを実行する。書込部１６１ｂは、演算部１６３による演算の結果を演算データ記憶部１３に書き込む。

具体的には、割当部１６２は、解ベクトルの要素と、解ベクトルの要素に対応する係数行列の要素とを、複数の演算部１６３のうちの当該係数行列の要素の行と予め対応付けられた演算部１６３に割り当てる。演算部１６３は、解ベクトルの要素及び係数行列の要素が割り当てられた場合、割り当てられた解ベクトルの要素と係数行列の要素との積及び積の和を求める演算を実行する。

ここで、演算データが図６に示すものである場合の例を挙げて、加速部１６の処理について説明する。演算データ記憶部１３は、解ベクトルであるベクトルｘの要素を記憶する。また、演算データ記憶部１３は、係数行列である行列Ａの要素を格納した配列valueと、行列Ａの要素の列を特定する配列colindと、行列Ａの要素の行を特定する配列jdptr及び配列permとを少なくとも記憶する。各配列の内容は、図６に示す通りである。

なお、ベクトルｘの要素は、第１の要素の一例である。また、行列Ａの要素は、第２の要素の一例である。また、配列valueは、第１のデータの一例である。また、配列colindは、第２のデータの一例である。また、配列jdptr及び配列permは、第３のデータの一例である。

また、演算データ記憶部１３は、少なくとも演算部１６３の数と同数のポートを備えているものとする。また、演算部１６３のそれぞれには、行列Ａの行のいずれかが予め対応付けられているものとする。なお、行列Ａの各行は、定数ベクトルであるベクトルｙの各要素に対応している。また、演算部１６３のそれぞれは、演算データ記憶部１３に備えられたポートのいずれかに対応しているものとする。

図９の疑似コードを用いて、加速部１６によるＳｐＭＶの処理について説明する。図９は、実施例１におけるＳｐＭＶの処理の一例を示す図である。図９のrho_maxは、配列jdptrのインデックスの最大値である。また、ｘ［］及びｙ［］は、それぞれベクトルｘ及びベクトルｙを格納した配列である。

例えば、ｊが０かつｉがjdptr［０］すなわち０である場合、図９に示すように、読込部１６１ａは、まずjdptr［０］、perm［０］、value［０］及びcolind［０］を読み込む。さらに、colind［０］は０なので、読込部１６１ａは、バースト転送により、x［０］、x［１］及びx［２］を読み込む。そして、割当部１６２は、jdptr［０］、perm［０］、value［０］、colind［０］及びx［０］を対応する演算部１６３に割り当てる。このとき、演算部１６３は、ｙ［perm［０］］にvalue［０］＊x［０］を足す。

次に、ｊが０のままｉが１になった場合、読込部１６１ａは、jdptr［０］、perm［０］、value［１］及びcolind［１］を読み込む。また、colind［１］は０であり、x［０］は読み込み済みであるため、読込部１６１ａは、ここではx［０］を読み込まない。そして、割当部１６２は、jdptr［０］、perm［０］、value［１］、colind［１］及びx［０］を対応する演算部１６３に割り当てる。ここで、perm［０］＝３なので、演算部１６３は、ｙ［３］にvalue［１］＊x［０］を足す。

さらに処理が進み、ｊが１かつｉがjdptr[１]、すなわち５になった場合、読込部１６１ａは、jdptr［１］、perm［１］、value［５］及びcolind［５］を読み込む。また、colind［５］は１であり、x［１］は読み込み済みであるため、読込部１６１ａは、ここではx［１］を読み込まない。そして、割当部１６２は、jdptr［１］、perm［１］、value［５］、colind［５］及びx［１］を対応する演算部１６３に割り当てる。ここで、perm［１］＝１なので、演算部１６３は、ｙ［１］にvalue［５］＊x［１］を足す。

さらに処理が進み、ｊが２かつｉがjdptr[２]、すなわち１０になった場合、読込部１６１ａは、jdptr［２］、perm［２］、value［１０］及びcolind［１０］を読み込む。また、colind［１０］は２であり、x［２］は読み込み済みであるため、読込部１６１ａは、ここではx［２］を読み込まない。そして、割当部１６２は、jdptr［２］、perm［２］、value［１０］、colind［１０］及びx［２］を対応する演算部１６３に割り当てる。ここで、perm［２］＝２なので、演算部１６３は、ｙ［２］にvalue［１０］＊x［２］を足す。

さらに処理が進み、ｊが２かつｉがjdptr[２]＋１、すなわち１１になった場合、読込部１６１ａは、value［１１］及びcolind［１１］を読み込む。また、colind［１１］は３であり、x［３］は読み込み済みでないため、読込部１６１ａは、ここでバースト転送により、x［３］、x［４］及びx［５］を読み込む。そして、割当部１６２は、jdptr［２］、perm［２］、value［１１］、colind［１１］及びx［３］を対応する演算部１６３に割り当てる。ここで、perm［２］＝２なので、このとき、演算部１６３は、ｙ［２］にvalue［１１］＊x［３］を足す。

このように、加速部１６は、連続的に読み込んだｘ[]の要素を使ってＳｐＭＶの演算を行うことができるため、ランダムアクセスを緩和し、ＣＧ法を高速化することができる。なお、１回の転送サイクルで読み取り取得可能な連続したｘ[]の要素を３つまでとし、演算データが図２のようなＣＲＳ形式である場合、ｘ[]の要素の読み取りには、少なくとも８回の転送サイクルが必要になる。一方、１回の転送サイクルで読み取り取得可能な連続したｘ[]の要素を３つまでとし、演算データが図６のようなＪＡＤ形式である場合、ｘ[]の要素の読み取りに必要な転送サイクルは２回である。

割当部１６２は、配列valueに格納された行列Ａの要素と、配列colindによって特定される当該行列Ａの要素の列に対応したベクトルｘの要素とを、配列jdptr及び配列permとともに演算部１６３に割り当てる。これにより、演算部１６３は、行列Ａの各要素の、ＪＡＤ形式におけるソート前の行インデックスを特定することができる。

図１０を用いて、演算部１６３について説明する。図１０は、実施例１における演算部の一例を示す図である。図１０に示すように、演算部１６３は、マルチポートメモリ１６３ａ、複数の乗算器１６３ｂ及び複数の乗算器１６３ｂのそれぞれと接続された複数の加算器１６３ｃを有する。また、演算部１６３は、マルチポートメモリ１６３ａ及び複数の加算器１６３ｃのそれぞれと接続された複数のパイプラインレジスタ１６３ｄを有する。

マルチポートメモリ１６３ａは、複数の書き込みポートと、複数の書き込みポートのそれぞれに対応した読み取りポートとを備え、部分和を係数行列の行と対応付けて記憶する。

マルチポートメモリ１６３ａは、書き込みポート１ｗ、２ｗ、３ｗ及び４ｗを有する。また、マルチポートメモリ１６３ａは、読み取りポート１ｒ、２ｒ、３ｒ及び４ｒを有する。書き込みポート１ｗ、２ｗ、３ｗ及び４ｗは、それぞれ読み取りポート１ｒ、２ｒ、３ｒ及び４ｒと対応している。なお、乗算器１６３ｂ、加算器１６３ｃ、パイプラインレジスタ１６３ｄ、書き込みポート及び読み取りポートの数は図１０に示すものに限られない。

乗算器１６３ｂは、ベクトルｘの要素と行列Ａの要素とを乗じた値を出力する。加算器１６３ｃは、複数の書き込みポートのうちのいずれかの書き込みポート及び当該書き込みポートに対応した読み取りポートと接続されている。加算器１６３ｃは、乗算器１６３ｂによって出力された値と、読み取りポートから読み取った第２の要素の行に対応する部分和とを足した値を書き込みポートからマルチポートメモリ１６３ａに書き込む。

演算部１６３は、複数の乗算器１６３ｂ、加算器１６３ｃ、書き込みポート及び読み取りポートの組み合わせを用いて、ベクトルｙの部分和の計算を並列して実行することができる。このため、割当部１６２は、各演算部１６３にベクトルｙの要素を担当要素として設定しておき、各担当要素の演算データを各演算部１６３に割り当てていくことができる。

［処理の流れ］
次に、本実施例における処理について、図１１を用いて説明する。図１１は、実施例１におけるＣＧ法による処理の一例を示すフローチャートである。まず、共役勾配法制御装置１０ａの変換部１１は、係数行列を圧縮し列指向形式に変換する（ステップＳ１０）。また、変換部１１は、変換した係数行列を疎行列ベクトル積演算装置１０ｂの演算データ記憶部１３に格納する。

次に、疎行列ベクトル積演算装置１０ｂの前処理部１５は、演算データ記憶部１３から疎行列を取得し、前処理を実行する（ステップＳ２０）。このとき、前処理部１５は、前処理制御部１４から指定された前処理を実行するために予め各処理部の再構成を行っておく。また、前処理部１５は、前処理済みの演算データを演算データ記憶部１３に格納する。

次に、加速部１６は、ＳｐＭＶを実行する（ステップＳ３０）。なお、加速部１６によるＳｐＭＶの実行については後に説明する。ここで、予め設定されたＣＧ法の終了条件が充足されていない場合（ステップＳ４０：Ｎｏ）、疎行列ベクトル積演算装置１０ｂは、Ｓ２０に戻って処理を繰り返す。一方、予め設定されたＣＧ法の終了条件が充足されている場合（ステップＳ４０：Ｙｅｓ）、共役勾配法制御装置１０ａ及び疎行列ベクトル積演算装置１０ｂは処理を終了する。例えば、終了条件終了条件は、所定の回数だけ反復が行われたことであってもよいし、定数ベクトルの残差が収束したことであってもよい。

図１２を用いて、疎行列ベクトル積演算装置１０ｂによるＳｐＭＶの処理について説明する。図１２は、実施例１におけるＳｐＭＶの処理の一例を示すフローチャートである。まず、加速部１６の読込部１６１ａは、演算データ記憶部１３から演算データを読み込む（ステップＳ３１０）。演算データには、列指向形式で表された係数行列及び解ベクトルの各要素が含まれる。

次に、割当部１６２は、演算データを各演算部１６３に割り当てる（ステップＳ３２０）。そして、演算部１６３は、割り当てられた演算データを使って演算を実行する（ステップＳ３３０）。なお、演算部１６３による演算処理については後に説明する。

ここで、演算データ記憶部１３に未読み込みの演算データがある場合（ステップＳ３４０：Ｙｅｓ）、読込部１６１ａは、Ｓ３１０に戻って処理を繰り返す。一方、演算データ記憶部１３に未読み込みの演算データがない場合（ステップＳ３４０：Ｎｏ）、疎行列ベクトル積演算装置１０ｂは処理を終了する。

図１３を用いて、演算部１６３による演算処理について説明する。図１３は、実施例１における演算処理の一例を示すフローチャートである。ここで、係数行列はＪＡＤ形式で表されているものとする。まず、演算部１６３には、演算データとして、value[ｉ]、ｘ[colind[ｉ]]、jdptr、permが入力される（ステップＳ３３１）。このとき、value[ｉ]及びｘ[colind[ｉ]]は乗算器１６３ｂへ入力される。また、jdptr及びpermは、パイプラインレジスタ１６３ｄへ入力される。

ここで、加算器１６３ｃは、パイプラインレジスタ１６３ｄに入力されたjdptr及びpermに従って、マルチポートメモリ１６３ａの読み取りポートのいずれかからｙ[perm［ｊ］]を読み取る（ステップＳ３３２）。次に、加算器１６３ｃは、ｙ[perm［ｊ］]に、value[ｉ]＊ｘ［colind[ｉ]］を足しこむ（ステップＳ３３３）。そして、加算器１６３ｃは、パイプラインレジスタ１６３ｄに入力されたjdptr及びpermに従って、ｙ[perm［ｊ］]を読み取った読み取りポートに対応する書き込みポートからｙ[perm［ｊ］]を書き込む（ステップＳ３３４）。

［効果］
以上説明したように、本実施例における疎行列ベクトル積演算装置の演算データ記憶部は、共役勾配法の制御を行う共役勾配法制御装置から入力された演算データを記憶する。演算データは、疎行列との積を計算可能なベクトルと、列指向形式で圧縮された疎行列とを含む。疎行列ベクトル積演算装置は、演算データ記憶部からベクトル及び疎行列を読み込む。疎行列ベクトル積演算装置は、読み込まれたベクトル及び疎行列を用いて疎行列ベクトル積の演算を実行する。疎行列ベクトル積演算装置は、演算の結果を演算データ記憶部に書き込む。このように、疎行列ベクトル積演算装置は、共役勾配法制御装置から演算データ記憶部に係数行列が入力された後は、当該疎行列の装置間の転送を行うことなくＳｐＭＶを反復して行うことができる。このため、疎行列ベクトル積演算装置は、ＣＧ法を高速化することができる。

また、疎行列ベクトル積演算装置は、ベクトルの要素である第１の要素と、第１の要素に対応する疎行列の非零の要素である第２の要素とを、複数の演算部のうちの第２の要素の行と予め対応付けられた演算部に割り当てる。演算部は、第１の要素及び第２の要素が割り当てられた場合、第１の要素と第２の要素との積及び積の和を求める演算を実行する。これにより、疎行列ベクトル積演算装置は、列指向形式で表された係数行列を使ったＳｐＭＶを実行することができる。

また、演算データ記憶部は、演算データとして、第２の要素を格納した第１のデータと、第２の要素の列を特定する第２のデータと、第２の要素の行を特定する第３のデータとを少なくとも記憶してもよい。疎行列ベクトル積演算装置は、第１のデータに格納された第２の要素と、第２のデータによって特定される当該第２の要素の列に対応したベクトルの要素とを、第３のデータとともに演算部に割り当てる。これにより、疎行列ベクトル積演算装置は、係数行列を配列として扱うことができる。

演算部のマルチポートメモリは、複数の書き込みポートと、複数の書き込みポートのそれぞれに対応した読み取りポートとを備え、部分和を疎行列の行と対応付けて記憶する。複数の乗算器は、第１の要素と第２の要素とを乗じた値を出力する。複数の加算器は、複数の書き込みポートのうちのいずれかの書き込みポート及び当該書き込みポートに対応した読み取りポートと接続され、乗算器によって出力された値と、読み取りポートから読み取った第２の要素の行に対応する部分和とを足した値を書き込みポートからマルチポートメモリに書き込む。これにより、疎行列ベクトル積演算装置は、各演算部の内部でさらに並列計算を行うことができる。

演算データ記憶部は、少なくとも演算部の数と同数のポートを備えていてもよい。これにより、疎行列ベクトル積演算装置は、演算部ごとの演算データの並列した読み取り、及びＳｐＭＶの並列処理を行うことができる。

演算データ記憶部は、ＥＬＬ形式又はＪＡＤ形式で圧縮された疎行列を記憶してもよい。これにより、疎行列ベクトル積演算装置は、列指向形式で表された係数行列を得ることができる。

前処理部は、予め設定された演算データに対する複数の前処理のそれぞれに対応した処理部を備え、指定された前処理に対応する処理部を用いて前処理を実行してもよい。前処理制御部は、共役勾配法制御装置によって指定された前処理を、前処理部に対して指定する。これにより、ユーザが共役勾配法制御装置に対して指定した前処理を、疎行列ベクトル積演算装置に実行させることができる。

前処理部は、演算部による疎行列ベクトル積の演算が完了するたびに前処理を実行してもよい。演算部は、前処理部による前処理が完了するたびに疎行列ベクトル積の演算を実行する。これにより、疎行列ベクトル積演算装置は、演算データの装置間での転送をともなうことなく、ＳｐＭＶの反復を行うことができる。

［システム］
また、各実施例において説明した各処理のうち、自動的に行われるものとして説明した処理の一部を手動的に行うこともできる。あるいは、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。

また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散や統合の具体的形態は図示のものに限られない。つまり、その全部又は一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。例えば、図４に示す前処理部１５と加速部１６とを統合してもよい。さらに、各装置にて行われる各処理機能は、その全部又は任意の一部が、所定のプロセッサ（ＣＰＵ、ＧＰＵ，ＦＰＧＡ等）及び当該プロセッサにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。

［ハードウェア構成］
図１４は、ハードウェア構成例を示す図である。図１４に示すように、シミュレーション装置１０は、通信インタフェース５０ａ、ＨＤＤ（Hard Disk Drive）５０ｂ、メモリ５０ｃ、ＣＰＵ１０ｄ及びＦＰＧＡ５０ｅを有する。

通信インタフェース５０ａは、他の装置の通信を制御するネットワークインタフェースカードなどである。ＨＤＤ５０ｂは、プログラムやデータなどを記憶する記憶装置の一例である。メモリ５０ｃの一例としては、ＳＤＲＡＭ（Synchronous Dynamic Random Access Memory）等のＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、フラッシュメモリ等が挙げられる。ＣＰＵ５０ｄは、共役勾配法制御装置１０ａとして機能する。また、ＦＰＧＡ５０ｅは、疎行列ベクトル積演算装置１０ｂとして機能する。

１０ シミュレーション装置
１０ａ 共役勾配法制御装置
１０ｂ 疎行列ベクトル積演算装置
１１ 変換部
１２ 更新部
１３ 演算データ記憶部
１４ 前処理制御部
１５ 前処理部
１６ 加速部
１６１ ロード部
１６１ａ 読込部
１６１ｂ 書込部
１６２ 割当部
１６３ 演算部
１６３ａ マルチポートメモリ
１６３ｂ 乗算器
１６３ｃ 加算器
１６３ｄ パイプラインレジスタ

Claims (9)

1. 共役勾配法の制御を行う共役勾配法制御装置から入力された演算データであって、疎行列との積を計算可能なベクトルと、列指向形式で圧縮された前記疎行列とを含む演算データを記憶する演算データ記憶部と、
   前記演算データ記憶部から前記ベクトル及び前記疎行列を読み込む読込部と、
   前記読込部によって読み込まれた前記ベクトル及び前記疎行列を用いて疎行列ベクトル積の演算を実行する演算部と、
   前記演算部による演算の結果を前記演算データ記憶部に書き込む書込部と、
   を有することを特徴とする疎行列ベクトル積演算装置。
2. 前記ベクトルの要素である第１の要素と、前記第１の要素に対応する前記疎行列の非零の要素である第２の要素とを、複数の前記演算部のうちの前記第２の要素の行と予め対応付けられた演算部に割り当てる割当部をさらに有し、
   前記演算部は、前記第１の要素及び前記第２の要素が割り当てられた場合、前記第１の要素と前記第２の要素との積及び前記積の和を求める演算を実行することを特徴とする請求項１に記載の疎行列ベクトル積演算装置。
3. 前記演算データ記憶部は、前記演算データとして、前記第２の要素を格納した第１のデータと、前記第２の要素の列を特定する第２のデータと、前記第２の要素の行を特定する第３のデータとを少なくとも記憶し、
   前記割当部は、前記第１のデータに格納された前記第２の要素と、前記第２のデータによって特定される当該第２の要素の列に対応した前記ベクトルの要素とを、前記第３のデータとともに前記演算部に割り当てることを特徴とする請求項２に記載の疎行列ベクトル積演算装置。
4. 前記演算部は、
   複数の書き込みポートと、前記複数の書き込みポートのそれぞれに対応した読み取りポートとを備え、部分和を前記疎行列の行と対応付けて記憶するマルチポートメモリと、
   前記第１の要素と前記第２の要素とを乗じた値を出力する複数の乗算器と、
   前記複数の書き込みポートのうちのいずれかの書き込みポート及び当該書き込みポートに対応した読み取りポートと接続され、前記乗算器によって出力された値と、前記読み取りポートから読み取った前記第２の要素の行に対応する部分和とを足した値を前記書き込みポートから前記マルチポートメモリに書き込む複数の加算器と、
   を有することを特徴とする請求項２又は３に記載の疎行列ベクトル積演算装置。
5. 前記演算データ記憶部は、少なくとも前記演算部の数と同数のポートを備えたことを特徴とする請求項１から４のいずれか１項に記載の疎行列ベクトル積演算装置。
6. 前記演算データ記憶部は、ＥＬＬ形式又はＪＡＤ形式で圧縮された前記疎行列を記憶することを特徴とする請求項１から５のいずれか１項に記載の疎行列ベクトル積演算装置。
7. 予め設定された前記演算データに対する複数の前処理のそれぞれに対応した処理部を備え、指定された前処理に対応する処理部を用いて前処理を実行する前処理部と、
   前記共役勾配法制御装置によって指定された前処理を、前記前処理部に対して指定する前処理制御部と、
   をさらに有することを特徴とする請求項１から６のいずれか１項に記載の疎行列ベクトル積演算装置。
8. 前記前処理部は、前記演算部による疎行列ベクトル積の演算が完了するたびに前処理を実行し、
   前記演算部は、前記前処理部による前処理が完了するたびに疎行列ベクトル積の演算を実行することを特徴とする請求項７に記載の疎行列ベクトル積演算装置。
9. 共役勾配法の制御を行う共役勾配法制御装置から入力された演算データであって、疎行列との積を計算可能なベクトルと、列指向形式で圧縮された前記疎行列とを含む演算データを記憶する演算データ記憶部を有する疎行列ベクトル積演算装置が、
   前記演算データ記憶部から前記ベクトル及び前記疎行列を読み込み、
   前記読み込む処理によって読み込まれた前記ベクトル及び前記疎行列を用いて疎行列ベクトル積の演算を実行し、
   演算を実行する処理による演算の結果を前記演算データ記憶部に書き込む
   処理を実行することを特徴とする疎行列ベクトル積演算方法。

Images