JP5262177B2 - ベクトル積の並列処理方法 - Google Patents

Description

本発明は、圧縮列格納法を用いたスパース行列とベクトルとの積を並列に処理する並列処理方法に関する。

数値解析を行なう場合などに使用するスパース行列のデータ記憶領域への格納方法（以下、単に「スパース行列の格納方法」という）には、いろいろな方法が提案されている。有限要素法などでは、圧縮列格納法（Ｈａｒｗｅｌｌ−ｂｏｅｉｎｇ ｆｏｒｍａｔ ｓｔｏｒａｇｅ ｍｅｔｈｏｄ）といわれる格納方法が広く使用されている。

圧縮列格納法とは、スパース行列の列ベクトルにある非ゼロ要素を圧縮し、圧縮した列ベクトルを順次１次元配列に格納する格納方法である。
例えば、ｎ行ｎ列のスパース行列（ｎは自然数）を行列ｍａｔ、行列ｍａｔの非ゼロ要素の総数をｎｚ（ｎｚは自然数）とすると、行列ｍａｔは、ｎｚ個の要素からなる１次元配列ａに格納される。

同時に、行列ｍａｔの各列における先頭の非ゼロ要素が配列ａの何番目の要素に格納されているかを示す情報が、ｎ個の要素からなる１次元配列ｎｆｃｎｚに格納される。
さらに、配列ａに格納された各要素が、行列ｍａｔの何番目の行ベクトルに属するかを示す情報が、ｎｚ個の要素からなる１次元配列ｎｒｏｗに格納される。

したがって、例えば、行列ｍａｔの第ｋ列目（ｋは自然数）の非ゼロ要素は、配列ａのｎｆｃｎｚ（ｋ）番目からｎｆｃｎｚ（ｋ＋１）−１番目の要素に格納される。そして、各要素が属する行ベクトルの行番号は、配列ｎｒｏｗのｎｆｃｎｚ（ｋ）番目からｎｆｃｎｚ（ｋ＋１）−１番目の要素に格納される。

ここで、行列ｍａｔと列ベクトルｘとの積ｙを求める場合、ｉ番目の行の要素ｙ（ｉ）は、１≦ｊ≦ｎについて（ｉ，ｊは自然数）、
ｙ（ｉ）＝ｙ（ｉ）＋ｍａｔ（ｉ，ｊ）＊ｘ（ｊ）
の計算を行なうことで得られる。したがって、１≦ｊ≦ｎの要素について上式の計算を行なうことにより、ｉ番目の行の要素ｙ（ｉ）の計算をすることができる。

上述の演算を行なう簡単な方法としては、例えば、各ＣＰＵにｙと同じ記憶領域ｙ１〜ｙｍ（ｍはＣＰＵの数を表す自然数）を割り当てるとともに、列の総数ｎを均等に分割して各ＣＰＵに割り当て、ＣＰＵ毎に当該分割した行列ｍａｔの部分行列と行の総数ｎを均等に分割して得る行列ｘの部分行列とについてベクトル積を算出し、結果を領域ｙ１〜ｙｍに格納し、最後にこれらの結果を加え合わせる方法がある。この演算方法は、簡単であるが、計算時間がかかりすぎてしまうという問題があった。

また、領域ｙｉ（ｉ＝１、２、・・・、ｍ）のバンド幅ｂａｎｄ（各列で対角要素の位置から最も離れた非ゼロ要素までの距離。例えば、対角要素を第ｉ列目の対角要素ｍａｔ（ｉ，ｉ）と同じ列にある最も離れた非ゼロ要素がｍａｔ（ｊ，ｉ）である場合のバンド幅はａｂｓ（ｊ−ｉ）となる。）を考慮し、２＊ｂａｎｄ＋ｎｎ（ｎｎはｎをＣＰＵ数ｍで均等に分割した大きさ（ｎ＋ｍ−１）／ｍ）に演算する範囲を制限することも可能であるが、バンド幅が大きく境界との間に非ゼロ要素がない場合には、大きな領域を使用することになる。そのため、上述したように、ｙｉに格納された演算結果を加算して積ｙを求める処理にコストがかかってしまい、並列化処理の効果が低下してしまうという問題があ
った。

上記技術に関連して、特許文献１には、反復解法により連立一次方程式を解くメモリ分散型並列計算機において、多様なデータ格納方法に対応して効率的な並列処理を行う反復解法について開示されている。

また、特許文献２には、プロセッサ間の同期回数を減らして、共有メモリベクトル並列計算機上でランダムスパース行列とベクトルとの積を高速に実行する行列ベクトル積演算システムについて開示されている。
特開平０９−２１２４８３号公報 特開２００１−２０９６３１号公報

本発明は、上述した問題に鑑みてなされたものであり、その解決しようとする課題は、圧縮列格納法を用いたスパース行列とベクトルとの積を効率よく並列に処理する並列処理方法を提供することである。

上記課題を解決するために、演算装置に、複数のスレッドを同時並列的に使用して行列と列ベクトルとの積を算出する演算方法であって、前記行列を所定の範囲で分割して部分行列を生成し、該部分行列についての行列のベクトル積の演算処理を、前記スレッド毎に割り当てる演算割り当て範囲を決定する演算割り当て範囲決定処理と、前記スレッドが算出した演算結果から、前記行列と列ベクトルとの積の一部を求める演算結果の更新処理を、前記スレッド毎に割り当てる更新割り当て範囲を決定する更新割り当て範囲決定処理と、前記演算割り当て範囲における前記部分行列を、前記行列を圧縮列格納法にしたがって圧縮して記憶する行列記憶手段から読み出し、該部分行列についての行列のベクトル積の演算処理を前記スレッドに実行させるベクトル演算処理と、該演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記行列と前記列ベクトルとの積を一時的に記憶する演算結果退避手段に記憶させる演算結果振り分け処理と、該演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、前記演算結果記憶手段に記憶されている演算結果を更新する演算結果更新処理と、を実行させる。

上記演算方法によると、各スレッドは、演算割り当て範囲決定処理で割り当てた演算割り当て範囲の演算処理を並列的に実行する。そして、演算結果が更新割り当て範囲の場合には行列と演算結果記憶手段に記憶し、演算結果が更新割り当て範囲でない場合には演算結果退避手段に記憶させる。そのため、各スレッドは、自身に割り当てられた更新割り当て範囲の更新処理についても並列的に実行することが可能となる。その結果、スパース行列と列ベクトルとの積を効率よく並列に処理することが可能となる。

以上に説明したように、本発明によると、圧縮列格納法を用いたスパース行列とベクトルとの積を効率よく並列に処理する並列処理方法を提供することが可能となる。

以下、本発明の実施形態について、図１〜図８に基づいて説明する。
（本実施例に係る並列処理の概要）
説明を簡単にするために、下記に示す４行４列の行列ｍａｔと４行１列の列ベクトルｘ
のベクトル積における並列処理（並列処理が可能な２つのＣＰＵ（ＣＰＵ＃１、ＣＰＵ＃２）を有する場合）について説明する。

なお、本実施例に係る並列処理を、行列ｍａｔが４行４列の場合に限定する趣旨ではない。必要に応じて、ｎ行ｎ列の行列ｍａｔについて下記の処理を行なうことで同様の効果を得ることができるのは当然である。また、並列処理を行なうＣＰＵの数についても２つに限定する趣旨ではなく、必要に応じて２以上のＣＰＵを使用してもよい。

・・・ （１）

・・・ （２）
ｙを４行１列の列ベクトルとすると、行列ｍａｔと列ベクトルｘのベクトル積は、次式によって求められる。

・・・ （３）
ここで、行列ｍａｔの列を２つに分割した下記の４行２列の部分行列ｍ１、ｍ２と、列ベクトルｘの行を２つに分割した下記に示す２行１列の部分列ベクトルｘ１、ｘ２を用いると、

・・・ （４）

・・・ （５）
式（３）は、下記のように変形することができる。

・・・ （６）
本実施例に係る行列のベクトル積では、図１に示すようにｍ１＊ｘ１’の演算処理をＣＰＵ＃１に割り当て、ｍ２＊ｘ２’の演算処理をＣＰＵ＃２に割り当てる。このように各ＣＰＵに割り当てられた演算処理を、以下「演算割り当て範囲」という。なお、演算装置１００が生成する各スレッドに当該演算処理を割り当てた場合も同様である。

また、各ＣＰＵによる演算結果ｙの更新処理において、（ｙ１，ｙ２）の更新処理をＣＰＵ＃１に割り当て、（ｙ３，ｙ４）の更新処理をＣＰＵ＃２に割り当てる。このように各ＣＰＵに割り当てられた更新処理を、以下「更新割り当て範囲」という。なお、演算装置１００が生成する各スレッドに当該更新処理を割り当てた場合も同様である。

例えば、ＣＰＵ＃１は、ＣＰＵ＃１自身による演算「ａ１，１＊ｘ１＋ａ１，２＊ｘ２」の結果と、ＣＰＵ＃２による演算「ａ１，３＊ｘ３＋ａ１，４＊ｘ４」の結果と、からｙ１を求める（ｙ１の値を更新する）。同様に、ＣＰＵ＃１は、ＣＰＵ＃１自身による演算「ａ２，１＊ｘ１＋ａ２，２＊ｘ２」の結果と、ＣＰＵ＃２による演算「ａ２，３＊ｘ３＋ａ２，４＊ｘ４」の結果と、からｙ２を求める（ｙ２の値を更新する）。

したがって、ＣＰＵ＃１自身による演算「ａ１，１＊ｘ１＋ａ１，２＊ｘ２」および「ａ２，１＊ｘ１＋ａ２，２＊ｘ２」、ＣＰＵ＃２による演算「ａ１，３＊ｘ３＋ａ１，４＊ｘ４」および「ａ２，３＊ｘ３＋ａ２，４＊ｘ４」は、ＣＰＵ＃１の更新割り当て範囲となる。

同様に、ＣＰＵ＃２は、ＣＰＵ＃１による演算「ａ３，１＊ｘ１＋ａ３，２＊ｘ２」の結果と、ＣＰＵ＃２自身による演算「ａ３，３＊ｘ３＋ａ３，４＊ｘ４」の結果と、から
ｙ３を求め、ＣＰＵ＃１による演算「ａ４，１＊ｘ１＋ａ４，２＊ｘ２」の結果と、ＣＰＵ＃２自身による演算「ａ４，３＊ｘ３＋ａ４，４＊ｘ４」の結果と、からｙ４を求めるので、
ＣＰＵ＃１による演算「ａ３，１＊ｘ１＋ａ３，２＊ｘ２」および「ａ４，１＊ｘ１＋ａ４，２＊ｘ２」、ＣＰＵ＃２自身による演算「ａ３，３＊ｘ３＋ａ３，４＊ｘ４」および「ａ４，３＊ｘ３＋ａ４，４＊ｘ４」は、ＣＰＵ＃２の更新割り当て範囲となる。

なお、並列処理を割り当て可能なＣＰＵをｍ個有する場合には、行列ｍａｔの列をｍ分割したｍ個の部分行列と、列ベクトルｘをｍ分割したｍ個の部分行列を用いればよい。この場合、各部分行列のベクトル積の演算処理を、各ＣＰＵに割り当てればよい。また、行列ｙの行をｍ分割した各部分行列の更新処理を、各ＣＰＵに割り当てればよい。なお、演算装置１００が生成する各スレッドに対して演算処理、更新処理を割り当ててもよい。
（本実施例に係る圧縮列格納法）
図２は、本実施例に係る圧縮列格納法の説明をする図である。

図２では、式（１）に示した行列ｍａｔを具体化した下記のスパース行列ｍａｔを圧縮列格納法により記憶装置に記憶する場合について説明する。なお、当該圧縮列格納法を当該スパース行列の場合に限定する趣旨でないのは当然である。

・・・ （７）
行列ｍａｔを圧縮列格納法を用いて記憶装置に記憶するためには、図２に示す１次元配列ａ、ｎｆｃｎｚおよびｎｒｏｗ（すなわち、配列ａ、ｎｆｃｎｚおよびｎｒｏｗを確保する記憶領域）を使用する。

配列ａには、行列ｍａｔの列の非ゼロ要素を圧縮して順次格納する。例えば、行列ｍａｔの第１列は（１ ０ ２ ０）であるから、０を除く要素「１」、「２」を配列ａに順次格納する。同様に、行列ｍａｔの第２列は（０ −４ ０ １）であるから、０を除く要素「−４」、「１」を配列ａに順次格納する。

配列ｎｒｏｗには、配列ａに格納された要素が属する行数を格納する。例えば、配列ａの第１、２番目の要素には、行列ｍａｔの第１列の非ゼロ要素「１」、「２」が圧縮されて格納されるので、配列ｎｒｏｗには、その要素の属する行番号「１」、「３」を格納する。

配列ｎｆｃｎｚには、行列ｍａｔにおける各列の最初の非ゼロ要素が格納されている配列ａ（または、配列ｎｒｏｗ）の要素番号を格納する。例えば、行列ｍａｔにおける第１列目の最初の非ゼロ要素は「１」である。この要素は、配列ａの第１番目の要素に格納されているので、配列ｎｆｃｎｚの第１番目の要素には、「１」を格納する。同様に、行列ｍａｔにおける第２列目の最初の非ゼロ要素は「−４」である。この要素は、配列ａの第３番目の配列に格納されているので、配列ｎｆｃｎｚの第２番目の配列には、「３」を格納する。

ここで、ｎ行ｎ列の行列ｍａｔを考える。この行列ｍａｔに非ゼロ要素がｎｚ個含まれる場合、配列ｎｆｃｎｚは大きさがｎ＋１の１次元配列となり、配列ａおよびｎｒｏｗは大きさがｎｚの１次元配列となる。

なお、配列ｎｆｃｎｚ（ｎ＋１）には、計算機の都合等を考慮して行列ｍａｔの仮想要素が格納される配列番号（例えば、図２に示す行列ｍａｔの要素ｍａｔ（１，５）が格納される配列ａの配列番号）が格納されるが、必須ではない。

また、これらの配列についてメモリ等の記憶装置に記憶領域を確保する場合、配列ａは行列ｍａｔの要素と同じデータ型（例えば、倍精度実数型など）の１次元配列として領域を確保し、配列ｎｆｃｎｚおよびｎｒｏｗは整数型の１次元配列として領域を確保すればよい。
（第１の実施例）
図３は、本実施例に係る行列ベクトル積の演算の概要を説明する図である。なお、図３に示す１次元配列ｗ、２次元配列ｉｗ、３次元配列ｎｅｎｔｒｙは、図２に示した行列ｍａｔについて演算を行なう場合の例を示しているが、これに限定する趣旨でないのは当然である。

配列ｗは、配列ａに格納されている要素について、次式による演算結果を格納する１次元配列である。以下、「演算結果退避領域」という。

・・・ （８）
したがって、配列ｗは配列ａと同じ大きさ（同じ構成）の１次元配列とすればよい。なお、本実施例では、更新割り当て範囲に属する演算結果については、配列ｗに格納せずに、演算結果を格納する配列ｙ（以下、「演算結果記憶領域」という）に格納するので、配列ａより小さい大きさであってもよい。

配列ｉｗは、配列ｗに格納されている要素について、式（８）による演算を行なった場合の要素ａｉ，ｊが属する行番号（すなわち、ｎｒｏｗ（ｉ）に格納されている値）を格納する２次元配列である。

ここで、配列ｉｗの１要素は、次の要素の位置を示す位置データを格納する領域ｉｗ（１，ｉ）（以下、「位置情報領域」という）と、上記ｎｒｏｗ（ｉ）の値を格納する領域ｉｗ（２，ｉ）（以下、「データ領域」という）と、で構成される。本実施例では、両領域ともに整数型として記憶領域を確保する。そして、複数の要素が連なってチェイン構造となる。以下、この配列ｉｗの１要素ｉｗ（１：２，ｉ）を「チェイン」という。

配列ｎｅｎｔｒｙは、更新割り当て範囲毎に、演算割り当て範囲の演算結果ｗ（ｉ）に対応するチェインｉｗ（１：２，ｉ）が格納されている先頭位置情報と終端位置情報を格納する。

本実施例では、図１で示したように更新割り当て範囲が２つ、演算割り当て範囲が２つなので、配列ｎｅｎｔｒｙは図３に示す２×２の２次元配列となる。ただし、各要素は２つの要素で構成され、１つの要素は、チェインｉｗ（１，ｉ）の先頭位置情報を格納する領域ｎｅｎｔｒｙ（１，ｉ，ｊ）（以下、「先頭位置領域」という）であり、もう１つの要素は、チェインｉｗ（ｉ）の終端位置情報を格納する領域ｎｅｎｔｒｙ（２，ｉ，ｊ）
（以下、「終端位置領域」という）である。本実施例では、両領域とも整数型として記憶領域を確保する。

したがって、配列ｎｅｎｔｒｙは、
（１）ｎｅｎｔｒｙ（１：２，１，１）に、ＣＰＵ＃１の更新割り当て範囲かつＣＰＵ＃１の演算割り当て範囲、
（２）ｎｅｎｔｒｙ（１：２，１，２）に、ＣＰＵ＃１の更新割り当て範囲かつＣＰＵ＃２の演算割り当て範囲、
（３）ｎｅｎｔｒｙ（１：２，２，１）に、ＣＰＵ＃２の更新割り当て範囲かつＣＰＵ＃１の演算割り当て範囲、
（４）ｎｅｎｔｒｙ（１：２，２，２）に、ＣＰＵ＃２の更新割り当て範囲かつＣＰＵ＃２の演算割り当て範囲、の演算結果に対応するチェインｉｗの先頭位置情報と終端位置情報を格納する。

なお、本実施例では、更新割り当て範囲かつ演算割り当て範囲の演算結果は、演算結果を格納する配列ｙに格納する。そのため、本実施例に係るｎｅｎｔｒｙ（１：２，１，１）とｎｅｎｔｒｙ（１：２，２，２）は使用しない。そこで、各先頭位置領域ｎｅｎｔｒｙ（１，１，１）とｎｅｎｔｒｙ（１，２，２）には、データが存在しない旨を示すデータ（図３の例では、「０」、以下「ターミナル情報」という）を格納する。

ここで、例えば、式（２）に示した列ベクトルｘが次式の場合について、行列ｍａｔと列ベクトルｘのベクトル積を考える。

・・・ （９）
ＣＰＵ＃１は、演算割り当て範囲における式（２）の演算を実行する。そして、当該演算が更新割り当て範囲である場合には、演算結果を格納する配列ｙに演算結果を格納し、当該演算が更新割り当て範囲でない場合には、演算結果を配列ｗに格納する。

例えば、ＣＰＵ＃１は、演算結果「１」（＝ａ（１）×ｘ（１））、「−８」（＝ａ（３）×ｘ（２））を、それぞれ配列ｙ（１）、ｙ（２）に格納し、演算結果「２」（＝ａ（２）×ｘ（１））、「２」（＝ａ（４）×ｘ（２））を、それぞれ配列ｗ（１）、ｗ（２）に格納する。

さらに、ＣＰＵ＃１は、演算結果「２」（＝ａ（２）×ｘ（１））を配列ｗ（１）に格納する時、新たなチェインのデータ領域ｉｗ（２，１）に、ａ（２）に対応するｎｒｏｗ（２）の値を格納し、当該チェインの位置情報領域ｉｗ（１，１）にターミナル情報を格納する。

同様に、ＣＰＵ＃１は、演算結果「２」（＝ａ（４）×ｘ（２））を配列ｗ（２）に格納する時、チェインの位置情報領域ｉｗ（１，１）に新たなチェインのデータ領域の位置情報３を格納する。そして、当該新たなチェインのデータ領域ｉｗ（２，２）に、ａ（４
）に対応するｎｒｏｗ（４）の値を格納し、当該チェインの位置情報領域ｉｗ（１，２）にターミナル情報を格納する。

ＣＰＵ＃２も同様の処理を行なう。
各ＣＰＵが自身に割り当てられた演算割り当て範囲について演算を完了すると、各ＣＰＵは、自身に割り当てられた更新割り当て範囲について、配列ｎｅｎｔｒｙを参照する。そして、当該ｎｅｎｔｒｙの先頭位置領域が示すチェインｉｗから、当該ｎｅｎｔｒｙの終端位置領域が示すチェインｉｗまで、たどりながらチェインｉｗのデータ領域を参照する。そして、当該データ領域に格納されている行番号に基づいて、各チェインｉｗに対応する配列ｗに格納されている演算結果を、配列ｙのいずれの要素に加算するか判断し、該当する配列ｙの要素に加算する。

例えば、ＣＰＵ＃１は、ｎｅｎｔｒｙ（１，１，２）の先頭位置領域を参照する。そして、配列ｉｗの３番目のチェインのデータ領域を参照する。すると、行番号が第１行とわかるので、当該チェインに対応する演算結果ｗ（３）は、ｙ（１）に加算するものと判断し、ｙ（１）の値に演算結果ｗ（３）を加算した値をｙ（１）に格納する。以上の処理を、チェインの終端まで行なう。ＣＰＵ＃２についても同様である。

以上の処理によって、行列ｍａｔと列ベクトルｘとのベクトル積の演算結果を配列ｙに得ることができる。
（本実施例に係る行列のベクトル積の演算の具体的な説明）
以下、ｎ行ｎ列の行列ｍａｔとｎ行１列の列ベクトルｘの行列ベクトル積ｙ＝ｍａｔ＊ｘについて説明する。

図４は、本実施例に係る行列のベクトル積の演算の処理を示すフローチャートである。
ステップＳ４０１において、演算装置１００は、行列ｍａｔの列の総数ｎを演算処理の割り当て可能なＣＰＵ数ｍで均等に分割したｎ行（ｎ／ｍ）列の部分行列ｍａｔ’１、ｍａｔ’２、ｍａｔ’３、・・・、ｍａｔ’ｍ、を生成する。

同様に、列ベクトルｘの行の総数ｎをｍで均等に分割した（ｎ／ｍ）行１列の部分列ベクトルｘ’１、ｘ’２、ｘ’３、・・・、ｘ’ｍ、を生成する。
そして、部分行列ベクトル積ｍａｔ’ｋ＊ｘ’ｋの処理をＣＰＵ＃ｋ（ｋは１以上ｍ以下の自然数）に割り当てる。この時、例えば、行列ベクトル積ｍａｔ’ｋ＊ｘ’ｋが、ＣＰＵ＃ｋの演算割り当て範囲となる。

ステップＳ４０２において、演算装置１００は、行列ベクトル積の演算結果を格納する行列ｙを演算処理を割り当て可能なＣＰＵ数ｍで均等に分割し、ｙ（ｍ＊（ｋ−１）＋１：ｍ＊ｋ）の更新処理をＣＰＵ＃ｋに割り当てる。この時、例えば、ｙ（ｍ＊（ｋ−１）＋１：ｍ＊ｋ）の更新処理がＣＰＵ＃ｋの更新割り当て範囲となる。

ここで、並列実行するスレッドの総数をｎｕｍｔｈｒｄ（＝ｍ）とし、各スレッドに割り当てられるスレッド番号（１〜ｍ）をｎｏｔｈｒｄとする。
ステップＳ４０３において、各スレッドは、自身に割り当てられた演算割り当て範囲について行列ベクトル積の計算（ｍａｔ（ｉ，ｊ）＊ｘ（ｊ））を行なう。

ステップＳ４０４において、スレッド番号ｎｏｔｈｒｄのスレッドによる演算結果が自身の更新割り当て範囲ｙ（ｍ＊（ｎｏｔｈｒｄ−１）＋１：ｍ＊ｎｏｔｈｒｄ）の要素の場合、当該スレッドは、処理をステップＳ４０５に移行し、当該演算結果を該当するｙの要素に格納する。また、既に値が格納されている場合には、当該値に演算結果を加算した値にｙの要素を更新する。

ステップＳ４０４において、スレッド番号ｎｏｔｈｒｄのスレッドによる演算結果が自身の更新割り当て範囲ｙ（ｍ＊（ｎｏｔｈｒｄ−１）＋１：ｍ＊ｎｏｔｈｒｄ）以外のｙの要素の場合、当該スレッドは、処理をステップＳ４０６に移行する。そして、当該スレッドは、配列ａと同じ大きさの作業領域の配列ｗに、自身が担当する配列ａの領域と同じ領域を対応付け、当該領域の先頭から順に演算結果を格納する。

ステップＳ４０７において、当該スレッドは、チェインを構成する１次元配列ｃｈａｉｎを複数格納する記憶領域であって、ｃｈａｉｎ（１）には次のチェインの先頭位置を示す位置情報を格納しｃｈａｉｎ（２）には当該演算結果に係る行列ｍａｔの計算要素の行数を格納する配列ｉｗに、チェインを追加し、当該演算結果に係る行列ｍａｔの計算要素の行数を格納する。

ステップＳ４０８において、当該スレッドは、更新割り当て範囲毎に、各演算割り当て範囲における演算より配列ｉｗに格納されたチェインの先頭位置と終端位置とを配列ｎｅｎｔｒｙに格納する。

例えば、スレッド番号ｎｏｔｈｒｄのスレッドによる演算結果が、スレッド番号ｋのスレッドに割り当てられた更新割り当て範囲である場合、当該スレッドは、配列ｉｗにチェインを追加して当該演算結果に係る行列ｍａｔの計算要素の行数を格納する。そして、当該チェインの先頭位置をｎｅｎｔｒｙ（１，ｋ，ｎｏｔｈｒｄ）に、当該チェインの終端位置をｎｅｎｔｒｙ（２，ｋ，ｎｏｔｈｒｄ）に、格納する。

なお、上述の配列ｗ、ｉｗおよびｎｅｎｔｒｙは、配列ａ、ｎｆｃｎｚ、ｎｒｏｗと同様に共有メモリに配置すればよい。
ステップＳ４０９において、当該スレッドは、自身に割り当てられた演算割り当て範囲について、行列のベクトル積の演算を全て実行したか否かを判別する。そして、まだ、演算割り当て範囲の全ての演算を完了していない場合、当該スレッドは、処理をステップＳ４０３に移行し、ステップＳ４０３〜Ｓ４０９の処理を繰り返す。また、演算割り当て範囲の全ての演算を完了した場合、当該スレッドは、処理をステップＳ４１０に移行する。

ステップＳ４１０において、当該スレッドは、バリア同期をとって、他のスレッドが自身に割り当てられた演算割り当て範囲の処理を完了するまで処理を停止する。
ステップＳ４１１において、スレッド番号がｎｏｔｈｒｄのスレッドは、他のスレッドに割り当てられた演算割り当て範囲、かつ自身のスレッドに割り当てられた範囲について、配列ｎｅｎｔｒｙがポイントする先頭チェインから終端チェインまでを参照し、当該チェインに格納された位置情報ｐに対応するｗ（ｐ）から演算結果を取得し、当該チェインに格納された行数ｑに対応するｙ（ｑ）に加算する。

以上の処理において、例えば、演算割り当て範囲決定処理はステップＳ４０１に対応し、更新割り当て範囲決定処理はステップＳ４０２に対応し、ベクトル演算処理はステップＳ４０３に対応し、演算結果振り分け処理はステップＳ４０４〜４０８に対応し、演算結果更新処理はステップ４１１に対応する。

また、例えば、行列記憶手段、演算結果退避手段、演算割り当て範囲決定手段、更新割り当て範囲決定手段、ベクトル演算処理手段、演算結果振り分け手段および演算結果更新手段は、図８に示す各ＣＰＵが、メモリモジュール等に配置された所定のプログラムに記載された命令を実行することによって実現される。

図５Ａ〜５Ｃは、本実施例に係る行列のベクトル積の演算処理の具体例を示すフローチ
ャートである。図５Ａ〜５Ｃに示すフローチャートは、ｎ行ｎ列のスパース行列ｍａｔとｎ行１列の列ベクトルｘとの演算処理を示す。

なお、以下の処理において、行列ｍａｔの非ゼロ要素の総数をｎｚ、並列実行するスレッドの数をｎｕｍｔｈｒｄとする。また、行列ｍａｔを圧縮列格納方式にしたがって非ゼロ要素を格納した１次元配列をａ（ｎｚ）とする。また、配列ａ（ｎｚ）に格納した各要素の行番号を、当該要素を格納したａ（ｎｚ）の配列番号と同じ配列番号の位置に格納する１次元配列をｎｒｏｗ（ｎｚ）とする。また、行列ｍａｔの各列における最初の非ゼロ要素が格納されているａ（ｎｚ）の位置を格納する１次元配列をｎｆｃｚ（ｎ）とする。

また、列ベクトルｘを格納する１次元配列をｘ（ｎ）、演算結果を格納する１次元配列をｙ（ｎ）とする。そして、図３で説明した作業領域を１次元配列ｗ（ｎｚ）、２次元配列ｉｗ（２，ｎｚ）、３次元配列ｎｅｎｔｒｙ（２，ｎｕｍｔｈｒｄ，ｎｕｍｔｈｒｄ）とする。

ステップＳ５００において、演算装置１００は、例えば、行列ｍａｔと列ベクトルｘが入力されると、行列ｍａｔを圧縮列格納方式にしたがって配列ａに格納する。そして、以下に示す行列のベクトル積ｍａｔ＊ｘの演算を開始する。

ステップＳ５０１において、演算装置１００は、並列処理を実行するスレッドをｎｕｍｔｈｒｄだけ確保する。さらに、演算装置１００は、配列ｎｅｎｔｒｙ（２，ｎｕｍｔｈｒｄ，ｎｕｍｔｈｒｄ）の領域を共有メモリに確保する。そして、演算装置１００は、ｎｅｎｔｒｙ＝０により当該ポインタをゼロクリアする。

ステップＳ５０２において、演算装置１００は、ステップＳ５０１で確保したスレッド数ｎｕｍｔｈｒｄのスレッドを生成する。そして、演算装置１００は、各スレッドにスレッド番号ｎｏｔｈｒｄ（１〜ｎｕｍｔｈｒｄ）を割り当てる。

ステップＳ５０３において、スレッドは、行列ｍａｔの次数ｎをｎｕｍｔｈｒｄで均等に分割する。そして、各スレッドが分担する区分ｊｓ〜ｊｅを決定する。この時、スレッドは、次式によりｊｓおよびｊｅを算出する。

ｎｎ ＝ （ｎ＋ｎｕｍｔｈｒｄ−１）／ｎｕｍｔｈｒｄ
ｊｓ ＝ ｎｎ＊（ｎｏｔｈｒｄ−１）＋１
ｊｅ ＝ ｍｉｎ（ｎ，ｎｎ＊ｎｏｔｈｒｄ）
例えば、スレッド番号ｋのスレッドの演算割り当て範囲は、行列ｍａｔの部分行列ｍａｔ（１：ｎ，ｎｎ＊（ｋ−１）＋１：ｍｉｎ（ｎ，ｎｎ＊ｋ））と、列ベクトルｘ（ｎ，１）の部分列ベクトルをｘ（ｎｎ＊（ｋ−１）＋１：ｍｉｎ（ｎ，ｎｎ＊ｋ），１）とすると、ｍａｔ（１：ｎ，ｎｎ＊（ｋ−１）＋１：ｍｉｎ（ｎ，ｎｎ＊ｋ））＊ｘ（ｎｎ＊（ｋ−１）＋１：ｍｉｎ（ｎ，ｎｎ＊ｋ），１） ・・・ （１０）
となる。そして、更新割り当て範囲は、行列ｙの部分行列ｙ（ｎｎ＊（ｋ−１）＋１：ｍｉｎ（ｎ，ｎｎ＊ｋ），１）となる。

ステップＳ５０４において、スレッドは、行列ベクトル積の演算結果を格納する領域ｙ（ｊｓ：ｊｅ）＝０．０ｄ０を実行して当該領域をゼロクリアする。
ステップＳ５０５において、スレッドは、バリア同期をとって、他のスレッドが自身に割り当てられた演算割り当て範囲の処理を開始できる状態になるまで処理を停止する。

ステップＳ５０６において、スレッドは、行列ベクトル積の演算に使用する各種変数ｊ、ｎｃｎｔ、ｎｂａｓｅ，ｎｓおよびｎｅに対して以下の処理を実行して初期化を行なう
。そして、スレッドは、以下に示す処理（ステップＳ５０７〜Ｓ５１９）を実行する。

ｊ ＝ ｊｓ
ｎｃｎｔ ＝ １
ｎｂａｓｅ ＝ ｎｆｃｎｚ（ｊｓ）−１
ｎｓ ＝ ｎｆｃｎｚ（ｊｓ）
ｎｅ ＝ ｎｆｃｎｚ（ｊｅ）−１
ステップＳ５０７において、スレッドは、行列ベクトル積の演算に使用する変数ｉに対して下記の処理を実行して初期化を行なう。

ｉ ＝ ｎｓ
ステップＳ５０８において、スレッドは、以降の演算（ステップＳ５１０又はＳ５１１）によって得る演算結果の更新割り当て範囲を担当するスレッド番号（以下、「インデックス」という）を以下の処理によって計算し、変数ｉｎｄｅｘに代入する。また、当該演算に使用するｘの要素を取り出して、変数ｘｘに代入する。

ｉｉ ＝ ｎｒｏｗ（ｉ）
ｉｎｄｅｘ ＝ （ｉｉ＋ｎｎ−１）／ｎｎ
ｘｘ ＝ ｘ（ｉｉ）
ステップＳ５０９において、スレッドは、ステップＳ５０８で算出したｉｎｄｅｘとｎｏｔｈｒｄとを比較する。そして、一致する場合、スレッドはステップＳ５１０に処理を移行する。また、一致しない場合、スレッドはステップＳ５１１に処理を移行する。

ステップＳ５１０において、スレッドは、以下のベクトル積の演算を実行し、演算結果を格納する行列ｙの要素を更新する。
ｙ（ｉｉ） ＝ ｙ（ｉｉ）＋ａ（ｉ）＊ｘｘ
ステップＳ５１１において、スレッドは、以下のベクトル積の演算を実行し、演算結果を演算結果退避領域に格納する。すなわち、演算結果を配列ｗに格納するとともに、配列ｉｗにチェインを追加して、配列ｉｗの位置情報領域に終端（ｔａｉｌ）を示す値０を格納し、同じくデータ領域にａ（ｉ）の行番号を格納する。

ｗ（ｎｂａｓｅ＋ｎｃｎｔ） ＝ ａ（ｉ）＊ｘｘ
ｉｗ（１，ｎｂａｓｅ＋ｎｃｎｔ） ＝ ０
ｉｗ（２，ｎｂａｓｅ＋ｎｃｎｔ） ＝ ｉｉ
ステップＳ５１２において、スレッドは、配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）が０か否かを判別する。配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）が０の場合、スレッドは、処理をステップ６１３に移行する。また、配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）が０でない場合、スレッドは、配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）に既にチェインが追加されていると判断し、処理をステップＳ５１４に移行する。

ステップＳ５１３において、スレッドは、以下の処理を行なって配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）にチェインを登録する。
ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ） ＝ ｎｂａｓｅ＋ｎｃｎｔ
ステップＳ５１４において、スレッドは、以下の処理を行なって配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）にチェインを登録する。

ｎｐｔｒ ＝ ｎｅｎｔｒｙ（２，ｉｎｄｅｘ，ｎｏｔｈｒｄ）
ｎｅｎｔｒｙ（２，ｉｎｄｅｘ，ｎｏｔｈｒｄ） ＝ ｎｂａｓｅ＊ｎｃｎｔ
ｉｗ（１，ｎｐｔｒ） ＝ ｎｂａｓｅ＋ｎｃｎｔ
ステップＳ５１５において、スレッドは、以下の処理を行なって変数ｎｃｎｔ、ｉをインクリメントする。

ｎｃｎｔ ＝ ｎｃｎｔ＋１
ｉ ＝ ｉ＋１
ステップＳ５１６において、スレッドは、変数ｉと変数ｎｅとを比較する。そして、ｉ＞ｎｅの場合、スレッドは、ステップＳ５１７に処理を移行する。また、ｉ≦ｅの場合、スレッドは、ステップＳ５０８に処理を移行する。そして、ステップＳ５０８〜Ｓ５１６の処理を繰り返し行なう。

ステップＳ５１７において、スレッドは、以下の処理を行なって変数ｊの値をインクリメントする。
ｊ ＝ ｊ＋１
ステップＳ５１８において、スレッドは、変数ｊと変数ｊｅとを比較する。そして、ｊ＞ｊｅの場合、スレッドは、ステップＳ５１９に処理を移行する。また、ｊ≦ｊｅの場合、スレッドは、ステップＳ５０７に処理を移行する。そして、ステップＳ５０７〜Ｓ５１８に処理を繰り返し行なう。

以上に示したステップＳ５０７〜Ｓ５１９の処理により、スレッドは、自身に割り当てられた演算割り当て範囲についての演算が完了する。
ステップＳ５１９において、スレッドは、バリア同期をとって、他のスレッドが自身に割り当てられた演算割り当て範囲の処理を完了するまで処理を停止する。

ステップＳ５２０において、スレッドは、変数ｉを以下の処理によって初期化する。
ｉ ＝ １
ステップＳ５２１において、スレッドは、以下の処理を行なって自身のスレッド番号以外のスレッドのスレッド番号を取得し、変数ｎｏｔｈｒｄ＿ｎｘｔに格納する。

ｎｏｔｈｒｄ＿ｎｘｔ ＝ ｍｏｄ（ｎｏｔｈｒｄ−１＋ｉ，ｎｕｍｔｈｒｄ）＋１
ステップＳ５２２において、スレッドは、配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）が０か否かを判別する。配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）が０の場合、スレッドは、他のスレッドが登録したチェインはないと判断し、ステップＳ５２６に処理を移行する。また、配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）が０でない場合、スレッドは、他のスレッドが登録したチェインが存在すると判断し、ステップＳ５２３に処理を移行する。

ステップＳ５２３において、スレッドは、以下の処理を行なって配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）に格納されている値を変数ｎｐｔｒに格納する。

ｎｐｔｒ ＝ ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）
ステップＳ５２４において、スレッドは、変数ｎｐｔｒが０か否かを判別する。そして、変数ｎｐｔｒが０の場合、スレッドは、ステップＳ５２６に処理を移行する。また、変数ｎｐｔｒが０でない場合、スレッドは、ステップＳ５２５に処理を移行する。

ステップＳ５２５において、スレッドは、以下の処理を行なって、自身の更新割り当て範囲ｙ（ｉｉ）について、他のスレッドの演算割り当て範囲の演算結果を更新する。そして、スレッドは、ステップＳ５２４に処理を移行する。

ｉｉ ＝ ｉｗ（２，ｎｐｔｒ）
ｙ（ｉｉ） ＝ ｙ（ｉｉ）＋ｗ（ｎｐｔｒ）
ｎｐｔｒ ＝ ｉｗ（１，ｎｐｔｒ）
ステップＳ５２６において、スレッドは、以下の処理を行なって変数ｉをインクリメントする。

ｉ ＝ ｉ＋１
ステップＳ５２７において、スレッドは、変数ｉとｎｕｍｔｈｒｄ−１と比較する。そして、ｉ＞ｎｕｍｔｈｒｄ−１の場合、スレッドは、ステップＳ５２８に処理を移行する。また、ｉ≦ｎｕｍｔｈｒｄ−１の場合、スレッドは、ステップＳ５２１に処理を移行する。そして、ステップＳ５０９〜Ｓ５２７の処理を繰り返す。

以上の処理によって、スレッドは、他のスレッドの演算割り当て範囲の演算結果を、自身の更新割り当て範囲に反映する処理が完了する。
ステップＳ５２８において、スレッドは、バリア同期をとって、他のスレッドの処理が完了するまで処理を停止する。そして、全スレッドの処理が完了すると、行列ｍａｔのベクトル積の演算結果が変数ｙ（１：ｎ）に得られる（ステップＳ５２９）。
（第２の実施例）
第１の実施例では、２次元配列ｉｗを行列のベクトル積演算に使用した場合について説明したが、例えば、配列ｉｗには１次元配列を使用してもよい。以下、配列ｉｗに１次元配列を使用した場合についての例を説明する。

図６は、本実施例に係る行列ベクトル積の演算の概要を説明する図である。なお、図３と同様に、図６に示す１次元配列ｗ、ｉｗ、３次元配列ｎｅｎｔｒｙは、図２に示した行列ｍａｔについて演算を行なう場合の例を示しているが、これに限定する趣旨でないのは当然である。

配列ｗは、図３と同様に、配列ａに格納されている要素について、式（８）による演算結果を格納する１次元配列である。したがって、配列ｗは配列ａと同じ大きさ（同じ構成）の１次元配列とすればよい。

配列ｉｗは、配列ｗに格納されている要素について、式（８）による演算を行なった場合の要素ａｉ，ｊが属する行番号（すなわち、ｎｒｏｗ（ｉ）に格納されている値）を格納する１次元配列である。当該行番号が格納される要素位置は、配列ｎｒｏｗ（または配列ａ）と同じ要素位置に格納される。例えば、当該行番号は、配列ｎｒｏｗ（ｉ）と同じ要素位置ｉｗ（ｉ）に格納される。

配列ｎｅｎｔｒｙは、図３と同様に、更新割り当て範囲毎に、演算割り当て範囲の演算結果ｗ（ｉ）に対応するｉｗ（ｉ）が格納されている先頭位置情報と終端位置情報を格納する。

本実施例では、図１に示したように更新割り当て範囲が２つ、演算割り当て範囲が２つなので、配列ｎｅｎｔｒｙは図３に示すように２×２の２次元配列となる。また、各要素は２つの要素で構成され、その１つの要素は先頭位置領域であり、もう１つの要素は終端位置領域である。本実施例では、両領域とも整数型として記憶領域を確保する。

例えば、式（２）に示した列ベクトルｘが式（９）の場合について、行列ｍａｔと列ベクトルｘのベクトル積を考える。
この時、ＣＰＵ＃１は、演算割り当て範囲における式（９）の演算を実行する。そして、当該演算が更新割り当て範囲である場合には、演算結果を格納する配列ｙに演算結果を格納し、当該演算が更新割り当て範囲でない場合には、演算結果を配列ｗに格納する。

例えば、ＣＰＵ＃１は、演算結果「１」（＝ａ（１）×ｘ（１））、「−８」（＝ａ（３）×ｘ（２））を、それぞれ配列ｙ（１）、ｙ（２）に格納し、演算結果「２」（＝ａ（２）×ｘ（１））、「２」（＝ａ（４）×ｘ（２））を、それぞれ配列ｗ（１）、ｗ（２）に格納する。

ここで、ＣＰＵ＃１は、演算結果「２」（＝ａ（２）×ｘ（１））を配列ｗ（１）に格納する時、新たなチェインのデータ領域ｉｗ（２）に、ａ（２）に対応するｎｒｏｗ（２）の値を格納する。

同様に、ＣＰＵ＃１は、演算結果「２」（＝ａ（４）×ｘ（２））を配列ｗ（２）に格納する時、新たなチェインのデータ領域ｉｗ（４）に、ａ（４）に対応するｎｒｏｗ（４）の値を格納する。

ＣＰＵ＃２も同様の処理を行なう。
各ＣＰＵが自身に割り当てられた演算割り当て範囲について演算を完了すると、各ＣＰＵは、自身に割り当てられた更新割り当て範囲について、配列ｎｅｎｔｒｙを参照する。そして、当該ｎｅｎｔｒｙの先頭位置領域が示すチェインｉｗから、当該ｎｅｎｔｒｙの終端位置領域が示すチェインｉｗまで、たどりながらチェインｉｗのデータ領域を参照する。そして、各チェインｉｗに対応する配列ｗに格納されている演算結果が、演算結果を格納する配列ｙのいずれの要素に加算されるものかを判断し、当該配列ｙの要素に加算する。

例えば、ＣＰＵ＃１は、ｎｅｎｔｒｙ（１，１，２）の先頭位置領域を参照する。そして、配列ｉｗの５番目のデータを参照する。すると、行番号が第１行とわかるので、当該チェインに対応する演算結果ｗ（５）は、ｙ（１）に加算するものと判断し、ｙ（１）の値に演算結果ｗ（５）を加算した値をｙ（１）に格納する。以上の処理を、チェインの終端まで行なう。ＣＰＵ＃２についても同様である。

以上の処理によって、行列ｍａｔと列ベクトルｘとのベクトル積の演算結果を配列ｙに得ることができる。
図７Ａ〜７Ｃは、本実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。図７Ａ〜７Ｃに示すフローチャートは、ｎ行ｎ列のスパース行列ｍａｔとｎ行１列の列ベクトルｘとの演算処理を示す。

なお、以下の処理において、行列ｍａｔの非ゼロ要素の総数をｎｚ、並列実行するスレッドの数をｎｕｍｔｈｒｄとする。また、行列ｍａｔを圧縮列格納方式にしたがって非ゼロ要素を格納した１次元配列をａ（ｎｚ）とする。また、配列ａ（ｎｚ）に格納した各要素の行番号を、当該要素を格納したａ（ｎｚ）の配列番号と同じ配列番号の位置に格納する１次元配列をｎｒｏｗ（ｎｚ）とする。また、行列ｍａｔの各列における最初の非ゼロ要素が格納されているａ（ｎｚ）の位置を格納する１次元配列をｎｆｃｚ（ｎ）とする。

また、列ベクトルｘを格納する１次元配列をｘ（ｎ）、演算結果を格納する１次元配列をｙ（ｎ）とする。そして、図３で説明した作業領域を１次元配列ｗ（ｎｚ）、ｉｗ（ｎｚ）、３次元配列ｎｅｎｔｒｙ（２，ｎｕｍｔｈｒｄ，ｎｕｍｔｈｒｄ）とする。

ステップＳ７００において、演算装置１００は、例えば、行列ｍａｔと列ベクトルｘが入力されると、行列ｍａｔを圧縮列格納方式にしたがって配列ａに格納する。そして、以下に示す行列のベクトル積演算ｍａｔ＊ｘを開始する。

ステップＳ７０１において、演算装置１００は、並列処理を実行するスレッドをｎｕｍｔｈｒｄだけ確保する。さらに、演算装置１００は、配列ｎｅｎｔｒｙ（２，ｎｕｍｔｈｒｄ，ｎｕｍｔｈｒｄ）の領域を共有メモリに確保する。そして、演算装置１００は、ｎｅｎｔｒｙ＝０により当該ポインタをゼロクリアする。

ステップＳ７０２において、演算装置１００は、ステップＳ７０１で確保したスレッド数ｎｕｍｔｈｒｄのスレッドを生成する。そして、演算装置１００は、各スレッドにスレッド番号ｎｏｔｈｒｄ（１〜ｎｕｍｔｈｒｄ）を割り当てる。

ステップＳ７０３において、スレッドは、行列ｍａｔの次数ｎをｎｕｍｔｈｒｄで均等に分割する。そして、各スレッドが分担する区分ｊｓ〜ｊｅを決定する。この時、スレッドは、次式によりｊｓおよびｊｅを算出し、各スレッドに対して、演算割り当て範囲と更新割り当て範囲を割り当てる。

ｎｎ ＝ （ｎ＋ｎｕｍｔｈｒｄ−１）／ｎｕｍｔｈｒｄ
ｊｓ ＝ ｎｎ＊（ｎｏｔｈｒｄ−１）
ｊｅ ＝ ｍｉｎ（ｎ，ｎｎ＊ｎｏｔｈｒｄ）
ステップＳ７０４において、スレッドは、行列ベクトル積の演算結果を格納する領域ｙ（ｊｓ：ｊｅ）＝０．０ｄ０により当該領域をゼロクリアする。

ステップＳ７０５において、スレッドは、バリア同期をとって、他のスレッドが自身に割り当てられた演算割り当て範囲の処理を開始できる状態になるまで処理を停止する。
ステップＳ７０６において、スレッドは、行列ベクトル積の演算に使用する各種変数ｊ、ｎｃｎｔ、ｎｓおよびｎｅに対して以下の処理を実行して初期化を行なう。そして、スレッドは、以下に示す処理（ステップＳ７０７〜Ｓ７１９）を各スレッドに実行させる。

ｊ ＝ ｊｓ
ｎｃｎｔ ＝ １
ｎｓ ＝ ｎｆｃｎｚ（ｊｓ）
ｎｅ ＝ ｎｆｃｎｚ（ｊｅ）−１
ステップＳ７０７において、スレッドは、行列ベクトル積の演算に使用する変数ｉに対して下記の処理を実行して初期化を行なう。

ｉ ＝ ｎｓ
ステップＳ７０８において、スレッドは、以降の演算（ステップＳ７１０又はＳ７１１）によって得る演算結果のインデックスを以下の処理によって計算し、変数ｉｎｄｅｘに代入する。また、当該演算に使用する列ベクトルｘの要素を取り出して、変数ｘｘに代入する。

ｉｉ ＝ ｎｒｏｗ（ｉ）
ｉｎｄｅｘ ＝ （ｉｉ＋ｎｎ−１）／ｎｎ
ｘｘ ＝ ｘ（ｉｉ）
ステップＳ７０９において、スレッドは、ステップＳ７０８で算出したｉｎｄｅｘとｎｏｔｈｒｄとを比較する。そして、一致する場合、スレッドはステップＳ７１０に処理を移行する。また、一致しない場合、スレッドはステップＳ７１１に処理を移行する。

ステップＳ７１０において、スレッドは、以下のベクトル積の演算を実行し、演算結果を格納する行列ｙの要素を更新する。
ｙ（ｉｉ） ＝ ｙ（ｉｉ）＋ａ（ｉ）＊ｘｘ
ステップＳ７１１において、スレッドは、以下のベクトル積の演算を実行し、演算結果
を演算結果退避領域に格納する。すなわち、演算結果を配列ｗに格納するとともに、配列ｉｗにチェインの終端（ｔａｉｌ）を示す０を格納する。

ｗ（ｉ） ＝ ａ（ｉ）＊ｘｘ
ｉｗ（ｉ） ＝ ０
ステップＳ７１２において、スレッドは、配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）が０か否かを判別する。配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）が０の場合、スレッドは、処理をステップ６１３に移行する。また、配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）が０でない場合、スレッドは、配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）に既にチェインが追加されていると判断し、処理をステップＳ７１４に移行する。

ステップＳ７１３において、スレッドは、以下の処理を行なって配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）にチェインを登録する。
ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ） ＝ ｉ
ステップＳ７１４において、スレッドは、以下の処理を行なって配列ｎｅｎｔｒｙ（１，ｉｎｄｅｘ，ｎｏｔｈｒｄ）にチェインを登録する。

ｎｐｔｒ ＝ ｎｅｎｔｒｙ（２，ｉｎｄｅｘ，ｎｏｔｈｒｄ）
ｎｅｎｔｒｙ（２，ｉｎｄｅｘ，ｎｏｔｈｒｄ） ＝ ｉ
ｉｗ（ｎｐｔｒ） ＝ ｉ
ステップＳ７１５において、スレッドは、以下の処理を行なって変数ｉをインクリメントする。

ｉ ＝ ｉ＋１
ステップＳ７１６において、スレッドは、変数ｉと変数ｎｅとを比較する。そして、ｉ＞ｎｅの場合、スレッドは、ステップＳ７１７に処理を移行する。また、ｉ≦ｅの場合、スレッドは、ステップＳ７０８に処理を移行する。そして、ステップＳ７０８〜Ｓ７１６の処理を繰り返し行なう。

ステップＳ７１７において、スレッドは、以下の処理を行なって変数ｊの値をインクリメントする。
ｊ ＝ ｊ＋１
ステップＳ７１８において、スレッドは、変数ｊと変数ｊｅとを比較する。そして、ｊ＞ｊｎの場合、スレッドは、ステップＳ７１９に処理を移行する。また、ｊ≦ｊｅの場合、スレッドは、ステップＳ７０７に処理を移行する。そして、ステップＳ７０７〜Ｓ７１８に処理を繰り返し行なう。

以上に示したステップＳ７０７〜Ｓ７１９の処理により、スレッドは、自身に割り当てられた演算割り当て範囲についての演算が完了する。
ステップＳ７１９において、スレッドは、バリア同期をとって、他のスレッドが自身に割り当てられた演算割り当て範囲の処理を完了するまで処理を停止する。

ステップＳ７２０において、スレッドは、変数ｉを以下の処理によって初期化する。
ｉ ＝ １
ステップＳ７２１において、スレッドは、以下の処理を行なって自身のスレッド番号以外のスレッドのスレッド番号を取得し、変数ｎｏｔｈｒｄ＿ｎｘｔに格納する。

ｎｏｔｈｒｄ＿ｎｘｔ ＝ ｍｏｄ（ｎｏｔｈｒｄ−１＋ｉ，ｎｕｍｔｈｒｄ）＋１
ステップＳ７２２において、スレッドは、配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏ
ｔｈｒｄ＿ｎｘｔ）が０か否かを判別する。配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）が０の場合、スレッドは、他のスレッドが登録したチェインはないと判断し、ステップＳ７２６に処理を移行する。また、配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）が０でない場合、スレッドは、他のスレッドが登録したチェインが存在すると判断し、ステップＳ７２３に処理を移行する。

ステップＳ７２３において、スレッドは、以下の処理を行なって配列ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）に格納されている値を変数ｎｐｔｒに格納する。

ｎｐｔｒ ＝ ｎｅｎｔｒｙ（１，ｎｏｔｈｒｄ，ｎｏｔｈｒｄ＿ｎｘｔ）
ステップＳ７２４において、スレッドは、変数ｎｐｔｒが０か否かを判別する。そして、変数ｎｐｔｒが０の場合、スレッドは、ステップＳ７２６に処理を移行する。また、変数ｎｐｔｒが０でない場合、スレッドは、ステップＳ７２５に処理を移行する。

ステップＳ７２５において、スレッドは、以下の処理を行なって、自身の更新割り当て範囲ｙ（ｉｉ）について、他のスレッドの演算割り当て範囲の演算結果を更新する。そして、スレッドは、ステップＳ７２４に処理を移行する。

ｉｉ ＝ ｎｒｏｗ（ｎｐｔｒ）
ｙ（ｉｉ） ＝ ｙ（ｉｉ）＋ｗ（ｎｐｔｒ）
ｎｐｔｒ ＝ ｉｗ（ｎｐｔｒ）
ステップＳ７２６において、スレッドは、以下の処理を行なって変数ｉをインクリメントする。

ｉ ＝ ｉ＋１
ステップＳ７２７において、スレッドは、変数ｉとｎｕｍｔｈｒｄ−１と比較する。そして、ｉ＞ｎｕｍｔｈｒｄ−１の場合、スレッドは、ステップＳ７２８に処理を移行する。また、ｉ≦ｎｕｍｔｈｒｄ−１の場合、スレッドは、ステップＳ７２１に処理を移行する。そして、ステップＳ７０９〜Ｓ７２７の処理を繰り返す。

以上の処理によって、スレッドは、他のスレッドの演算割り当て範囲の演算結果を、自身の更新割り当て範囲に反映する処理が完了する。
ステップＳ７２８において、スレッドは、バリア同期をとって、他のスレッドの処理が完了するまで処理を停止する。そして、全スレッドの処理が完了すると、行列ｍａｔのベクトル積の演算結果が変数ｙ（１：ｎ）に得られる（ステップＳ８２９）。

図８は、本実施例に係る行列のベクトル積演算を実行する演算装置１００の構成例を示す図である。
図８に示す演算装置１００は、複数のメモリモジュール＃１、＃２、・・・、＃ｎと、複数のＣＰＵ＃１、＃２、・・・、＃ｎと、相互結合網９０１と、を少なくとも備える共有メモリ型演算装置である。

各メモリモジュール＃１、＃２、・・・、＃ｎは、各ＣＰＵ＃１、＃２、・・・、＃ｎと、相互結合網９０１を介して接続されている。そして、例えば、本実施例に係るベクトル演算に使用する配列ａ、ｎｆｃｎｚ、ｎｒｏｗ、ｗ，ｉｗ，ｎｅｎｔｒｙ等の領域を提供する共有メモリとして使用される。

各ＣＰＵ＃１、＃２、・・・、＃ｎは、Ｌ２キャッシュおよびバスインタフェースと、Ｌ１キャッシュを有する２つのｃｐｕコアと、を備える。そして、各ｃｐｕコアは、バス
インタフェースおよび相互結合網９０１を介して他のｃｐｕコアやメモリモジュールとアクセス可能である。

なお、図８は、１つのＣＰＵユニットにｃｐｕコアが２ユニット搭載されている場合について示しているが、この構成に限定する趣旨ではないのは当然である。例えば、１つのＣＰＵユニットにｃｐｕコアが４ユニット搭載されていてもよい。

その他、図示しないが、演算装置１００には、例えば、本実施例に係る行列のベクトル積演算を実行するプログラム等を記憶するために磁気ディスク装置等で構成される記憶装置も備えてもよい。

以上に説明したように、本実施例に係る行列のベクトル積の演算方法は、各スレッドに対して、別個に演算割り当て範囲と更新割り当て範囲を割り当てる。各スレッドは、並列的に自身に割り当てた演算割り当て範囲の演算を実施し、演算結果が更新割り当て範囲である場合には演算結果記憶領域に直接記憶し、演算結果が更新割り当て範囲でない場合には演算結果退避領域に記憶する。そして、全スレッドの演算処理終了後、各スレッドは、並列的に自身に割り当てられた更新割り当て範囲の更新処理を実施する。

したがって、各スレッドは、演算割り当て範囲における演算処理と、更新割り当て範囲における更新処理と、を並列に実行することが可能となる。その結果、行列のベクトル積を効率よく並列処理することが可能となる。

一般に、圧縮列格納法で格納されたスパース行列と列ベクトルとの積は、逐次プログラムを単純にＯｐｅｎＭＰ ＦｏｒｔｒａｎのＯＣＬ（Ｏｂｊｅｃｔ Ｃｏｎｓｔｒａｉｎ
Ｌａｎｇｕａｇｅ）挿入などでは簡単に並列化することはできないが、本実施例に係る演算方法によれば、上述の理由から、簡単かつ効率的に並列化することが可能となる。

その結果、例えば、スパース行列の連立１次方程式を反復解法で解くときにも、行列のベクトル積を効率よく繰り返し計算することが可能となる。そのため、例えば、図８に示したＳＭＰ（Ｓｙｍｍｅｔｒｉｃ Ｍｕｌｔｉｐｌｅ Ｐｒｏｃｅｓｓｏｒ）システムの演算装置を用いて圧縮列格納法で格納された行列のベクトル積を効率的かつ並列化して行なうことが可能となる。

また、行列のベクトル積を使う連立１次方程式で使用する反復法（例えば、ＢＩＣＧＳＴＡＢ（Ｌ）法）などに本実施例に係る演算方法を使用することにより、ほぼ線形な台数効果を得ることができる。

さらに、複数のＣＰＵで行列のベクトル積を並列に計算する場合、行方向に行を均等に分割したそれぞれの区間ごとの、区間に属する列に存在する非ゼロの対角要素からの距離に関する分布がほぼ同じ場合に、より高い並列処理による台数効果を得ることができる。十分に大きな問題の場合にも、ほぼ線形な台数効果を引き出すことができる。

以上の実施例１〜ｎを含む実施形態に関し、さらに以下の付記を開示する。
（付記１） 複数のスレッドを同時並列的に使用して行列と列ベクトルとの積を算出するプログラムであって、
前記行列を所定の範囲で分割して部分行列を生成し、該部分行列についての行列のベクトル積の演算処理を、前記スレッド毎に割り当てる演算割り当て範囲を決定する演算割り当て範囲決定処理と、
前記スレッドが算出した演算結果から、前記行列と列ベクトルとの積の一部を求める演算結果の更新処理を、前記スレッド毎に割り当てる更新割り当て範囲を決定する更新割り
当て範囲決定処理と、
前記演算割り当て範囲における前記部分行列を、前記行列を圧縮列格納法にしたがって圧縮して記憶する行列記憶手段から読み出し、該部分行列についての行列のベクトル積の演算処理を前記スレッドに実行させるベクトル演算処理と、
該演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記行列と前記列ベクトルとの積を一時的に記憶する演算結果退避手段に記憶させる演算結果振り分け処理と、
該演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、前記演算結果記憶手段に記憶されている演算結果を更新する演算結果更新処理と、
を演算装置に実行させるプログラム。
（付記２） 前記演算割り当て範囲決定処理は、前記行列の列を均等幅に分割して部分行列を生成する、
ことを特徴とする付記１に記載のプログラム。
（付記３） 前記更新割り当て範囲決定処理は、前記行列と前記列ベクトルとの積を格納する列ベクトルの列を均等幅に分割して部分列ベクトルを生成し、該部分列ベクトルの更新処理を、前記スレッド毎に割り当てる、
ことを特徴とする付記１に記載のプログラム。
（付記４） 前記行列記憶手段は、
前記行列の非ゼロ要素を列ごとに抽出して記憶する第１の配列と、
該第１の配列の各要素が属する前記行列における行番号を記憶する第２の配列と、
前記第１および第２の配列の要素位置であって前記行列の各列の先頭要素が格納される要素位置を記憶する第３の配列と、
を備えることを特徴とする付記１に記載のプログラム。
（付記５） 前記演算結果退避手段は、
前記更新割り当て範囲外の演算結果を記憶する第１の配列と、
該演算結果の記憶する順に対応して、該演算に使用した前記行列の要素の行番号を記憶する第２の配列と、
前記更新割り当て範囲毎に、該更新割り当て範囲の演算結果が記憶されている前記第２の配列の位置を記憶する第３の配列と、
を備えることを特徴とする付記１に記載のプログラム。
（付記６） 前記演算結果更新処理は、
前記第３の配列から、他のスレッドが算出した更新割り当て範囲の演算結果を得る演算に使用した前記行列の要素の行番号が記憶されている前記第２の配列の位置を取得し、該行番号に応じて、前記演算結果記憶手段に記憶されている更新すべき演算結果を特定する、
を備えることを特徴とする付記５に記載のプログラム。
（付記７） 前記第２の配列は、前記行番号を記憶する第１の要素と次の行番号が記憶されている要素の位置を記憶する第２の要素とを１つの単位とするチェイン構造を有する、
ことを特徴とする付記５に記載のプログラム。
（付記８） 前記第１の配列には、前記行列の非ゼロ要素を列ごとに抽出して記憶する配列の要素の位置と同じ位置に、該要素から算出される演算結果が格納され、前記第２の配列には、前記行列の非ゼロ要素を列ごとに抽出して記憶する配列の要素の位置と同じ位置に、該要素の演算に使用した前記行列の要素の行番号を記憶する、
ことを特徴とする付記５に記載のプログラム。
（付記９） 複数のスレッドを同時並列的に使用して行列と列ベクトルとの積を算出する演算方法であって、
前記行列を所定の範囲で分割して部分行列を生成し、該部分行列についての行列のベクトル積の演算処理を、前記スレッド毎に割り当てる演算割り当て範囲を決定する演算割り当て範囲決定処理と、
前記スレッドが算出した演算結果から、前記行列と列ベクトルとの積の一部を求める演算結果の更新処理を、前記スレッド毎に割り当てる更新割り当て範囲を決定する更新割り当て範囲決定処理と、
前記演算割り当て範囲における前記部分行列を、前記行列を圧縮列格納法にしたがって圧縮して記憶する行列記憶手段から読み出し、該部分行列についての行列のベクトル積の演算処理を前記スレッドに実行させるベクトル演算処理と、
該演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記行列と前記列ベクトルとの積を一時的に記憶する演算結果退避手段に記憶させる演算結果振り分け処理と、
該演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、前記演算結果記憶手段に記憶されている演算結果を更新する演算結果更新処理と、
を演算装置に実行させる方法。
（付記１０） 複数のスレッドを同時並列的に使用して行列と列ベクトルとの積を算出する演算装置であって、
前記行列を圧縮列格納法にしたがって圧縮して記憶する行列記憶手段と、
前記行列と前記列ベクトルとの積を一時的に記憶する演算結果退避手段と、
前記行列を所定の範囲で分割して部分行列を生成し、該部分行列についての行列のベクトル積の演算処理を、前記スレッド毎に割り当てる演算割り当て範囲を決定する演算割り当て範囲決定手段と、
前記スレッドが算出した演算結果から、前記行列と列ベクトルとの積の一部を求める演算結果の更新処理を、前記スレッド毎に割り当てる更新割り当て範囲を決定する更新割り当て範囲決定手段と、
前記演算割り当て範囲における前記部分行列を前記行列記憶手段から読み出し、該部分行列についての行列のベクトル積の演算処理を前記スレッドに実行させるベクトル演算処理手段と、
該演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記演算結果退避手段に記憶させる演算結果振り分け手段と、
該演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、前記演算結果記憶手段に記憶されている演算結果を更新する演算結果更新手段と、
を備える演算装置。

本実施例に係る演算割り当て範囲および更新割り当て範囲を説明する図である。 本実施例に係る圧縮列格納法の説明をする図である。 第１の実施例に係る行列ベクトル積の演算の概要を説明する図である。 第１の実施例に係る行列のベクトル積の演算の処理を示すフローチャートである。 第１の実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。 第１の実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。 第１の実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。 第２の実施例に係る行列ベクトル積の演算の概要を説明する図である。 第２の実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。 第２の実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。 第２の実施例に係る行列のベクトル積の演算処理の具体例を示すフローチャートである。 本実施例に係る行列のベクトル積演算を実行する演算装置１００の構成例を示す図である。

Claims (9)

1. 複数のスレッドを使用して行列と列ベクトルとの積を算出するプログラムであって、
   行列の非ゼロ要素を列ごとに抽出して記憶する第１の配列と該第１の配列の各要素が属する行列における行番号を記憶する第２の配列と前記第１の配列および前記第２の配列の要素位置であって行列の各列における非ゼロ要素の先頭要素が格納される要素位置を記憶する第３の配列とを含む行列記憶手段に圧縮列格納法にしたがって圧縮して記憶した演算対象の行列を所定の範囲で分割した複数の部分行列それぞれについての行列のベクトル積の演算処理を複数の前記スレッドそれぞれに割り当てる演算割り当て範囲を決定し、
   複数の前記スレッドが算出した演算結果から前記演算対象の行列と列ベクトルとの積の一部を求める更新処理を、複数の前記スレッドそれぞれに割り当てる更新割り当て範囲を決定し、
   複数の前記スレッドそれぞれが、演算割り当て範囲における部分行列を前記行列記憶手段から読み出し、読み出した該部分行列についての行列のベクトル積の演算処理を並列に実行し、
   複数の前記スレッドそれぞれが、前記演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記演算対象の行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記演算対象の行列と列ベクトルとの積を記憶する演算結果退避手段に記憶し、
   複数の前記スレッドそれぞれが、前記演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、読み出した前記演算結果を前記演算結果記憶手段に記憶されている更新割り当て範囲の演算結果に加算して前記演算対象の行列と列ベクトルとの積を算出する、
   処理を演算装置に実行させるプログラム。
2. 前記演算結果退避手段は、
   更新割り当て範囲外の演算結果を記憶する第４の配列と、
   該第４の配列に演算結果を記憶する順に対応して、演算結果を得た演算に使用した行列の要素の行番号を記憶する第５の配列と、
   演算割り当て範囲毎に、演算割り当て範囲の演算結果が記憶されている前記第５の配列の範囲を記憶する第６の配列と、
   を備えることを特徴とする請求項１に記載のプログラム。
3. 複数のスレッドを使用して行列と列ベクトルとの積を算出する演算方法であって、
   行列の非ゼロ要素を列ごとに抽出して記憶する第１の配列と該第１の配列の各要素が属する行列における行番号を記憶する第２の配列と前記第１の配列および前記第２の配列の要素位置であって行列の各列における非ゼロ要素の先頭要素が格納される要素位置を記憶する第３の配列とを含む行列記憶手段に圧縮列格納法にしたがって圧縮して記憶した演算対象の行列を所定の範囲で分割した複数の部分行列それぞれについての行列のベクトル積の演算処理を複数の前記スレッドそれぞれに割り当てる演算割り当て範囲を決定し、
   複数の前記スレッドが算出した演算結果から前記演算対象の行列と列ベクトルとの積の一部を求める更新処理を、複数の前記スレッドそれぞれに割り当てる更新割り当て範囲を決定し、
   複数の前記スレッドそれぞれが、演算割り当て範囲における部分行列を前記行列記憶手段から読み出し、読み出した該部分行列についての行列のベクトル積の演算処理を並列に実行し、
   複数の前記スレッドそれぞれが、該演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記演算対象の行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記演算対象の行列と列ベクトルとの積を記憶する演算結果退避手段に記憶し、
   複数の前記スレッドそれぞれが、前記演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、読み出した前記演算結果を前記演算結果記憶手段に記憶されている更新割り当て範囲の演算結果に加算して前記演算対象の行列と列ベクトルとの積を算出する、
   処理を演算装置に実行させる方法。
4. 複数のスレッドを使用して行列と列ベクトルとの積を算出する演算装置であって、
   行列の非ゼロ要素を列ごとに抽出して記憶する第１の配列と該第１の配列の各要素が属する行列における行番号を記憶する第２の配列と前記第１の配列および前記第２の配列の要素位置であって行列の各列における非ゼロ要素の先頭が格納される要素位置を記憶する第３の配列とを含む行列記憶手段と、
   行列と列ベクトルとの積を記憶する演算結果退避手段と、
   前記圧縮列格納法にしたがって圧縮して記憶した演算対象の行列を所定の範囲で分割した複数の部分行列それぞれについての行列のベクトル積の演算処理を複数の前記スレッドそれぞれに割り当てる演算割り当て範囲を決定する演算割り当て範囲決定手段と、
   複数の前記スレッドが算出した演算結果から前記演算対象の行列と列ベクトルとの積の一部を求める更新処理を、複数の前記スレッドそれぞれに割り当てる更新割り当て範囲を決定する更新割り当て範囲決定手段と、
   複数の前記スレッドそれぞれが、前記演算割り当て範囲における部分行列を前記行列記憶手段から読み出し、読み出した該部分行列についての行列のベクトル積の演算処理を並列に実行するベクトル演算処理手段と、
   複数の前記スレッドそれぞれが、前記演算結果が前記更新割り当て範囲か否かを判別し、前記更新割り当て範囲の場合には前記演算結果を、前記演算対象の行列と列ベクトルとの積を記憶する演算結果記憶手段に記憶し、前記更新割り当て範囲でない場合には前記演算結果を、前記演算結果退避手段に記憶する演算結果振り分け手段と、
   複数の前記スレッドそれぞれが、前記演算結果退避手段から他のスレッドが算出した更新割り当て範囲の演算結果を読み出し、読み出した前記演算結果を前記演算結果記憶手段に記憶されている更新割り当て範囲の演算結果に加算して前記演算対象の行列と列ベクトルとの積を算出する演算結果更新手段と、
   を備える演算装置。
5. 前記演算割り当て範囲の決定において、前記演算対象の行列として前記行列記憶手段に記憶されたスパース行列を列単位で前記スレッドの総数に均等に分割した複数の部分行列それぞれについての行列のベクトル積の演算処理を複数の前記スレッドそれぞれに割り当てる、
   ことを特徴とする請求項１に記載のプログラム。
6. 前記部分行列についての行列のベクトル積の演算処理は、複数の前記スレッドそれぞれが、演算割り当て範囲における部分行列の列に含まれる非ゼロ要素を前記行列記憶手段から１要素ずつ読み出して、読み出した要素の列番号に対応する前記列ベクトルの要素との積を算出することによって行なわれる、
   ことを特徴とする請求項１に記載のプログラム。
7. 前記更新割り当て範囲は、前記演算対象の行列と列ベクトルの積によって得られる演算結果に含まれる要素の行の位置に基づいて決定され、
   複数の前記スレッドのうちの第１のスレッドに演算割り当て範囲として割り当てられた第１の部分行列の演算処理の演算結果に含まれる第１の要素が、前記第１のスレッドに割り当てられた第１の更新割り当て範囲に含まれる行番号のいずれかに対応する場合、前記第１の要素は、前記第１の更新割り当て範囲であると判別し、前記第１の要素が、前記第１の更新割り当て範囲に含まれる行番号のいずれでもない場合、前記第１の要素は、前記第１の更新割り当て範囲でないと判別する、
   ことを特徴とする請求項１に記載のプログラム。
8. 前記更新割り当て範囲には、前記演算対象の行列と列ベクトルの積によって得られる演算結果に含まれる要素を行単位で分割した複数の行の範囲を使用し、複数の前記行の範囲それぞれに、複数の前記スレッドそれぞれにあらかじめ割り当てられた識別番号を割り当てる、
   ことを特徴とする請求項１に記載のプログラム。
9. 複数の前記スレッドそれぞれが演算割り当て範囲の演算処理を完了したことをバリア同期を使用して確認すると、
   複数の前記スレッドそれぞれは、前記演算結果退避手段を参照し、前記第６の配列から特定される他のスレッドの演算割り当て範囲の演算結果のうち、前記第５の配列から特定される更新割り当て範囲の演算結果を前記第４の配列から読み出し、読み出した前記演算結果を前記演算結果記憶手段に記憶されている更新割り当て範囲の演算結果に加算して前記演算対象の行列と列ベクトルとの積を算出する、
   ことを特徴とする請求項２に記載のプログラム。