JP6083300B2 - プログラム、並列演算方法および情報処理装置 - Google Patents

Description

本発明はプログラム、並列演算方法および情報処理装置に関する。

スーパーコンピュータなどの演算能力の高いコンピュータを利用して、科学技術計算などの大規模数値計算が行われることがある。大規模数値計算は、次数の大きな行列の演算を含むことが多い。例えば、流体解析や構造解析などの分野では、連立方程式の係数を表した係数行列を用いて大規模な連立方程式を解くことがある。また、例えば、回路解析や振動解析などの分野では、大規模な行列の固有値を求めることがある。コンピュータを利用した行列演算では、行列とベクトルとの積を反復的に計算することで近似解を求めることができる。例えば、有限要素法に従って係数行列とベクトルとの積を反復的に計算することで、解析的に解くことが難しい微分方程式の近似解を求めることができる。

大規模数値計算で利用する行列は、値が０である要素（零要素）の割合が大きく、値が０でない要素（非零要素）の割合が小さい疎行列（スパース行列）であることがある。スパース行列は、零要素も含めて行列構造通りのデータとして表現するとデータ量が大きくなり非効率的であることから、圧縮格納法によって零要素が省略された圧縮データとして表現され得る。圧縮格納法には、圧縮列格納（ＣＣＳ：Compressed Column Storage）法と圧縮行格納（ＣＲＳ：Compressed Row Storage）法とが含まれる。

圧縮列格納法では、Ｎ×Ｍ行列に含まれる要素が、列優先の順序（１行１列，２行１列，…，Ｎ行１列，１行２列，２行２列，…，１行Ｍ列，…，Ｎ行Ｍ列の順）で探索されて、行列から非零要素のみが抽出される。そして、上記の順序に従って非零要素の値を列挙した第１の配列と、各非零要素の行番号を列挙した第２の配列と、第１の配列の中で列が変わる位置を列挙した第３の配列とが生成される。圧縮行格納法では、Ｎ×Ｍ行列に含まれる要素が、行優先の順序（１行１列，１行２列，…，１行Ｍ列，２行１列，２行２列，…，Ｎ行１列，…，Ｎ行Ｍ列の順）で探索されて、行列から非零要素のみが抽出される。そして、非零要素の値を列挙した第１の配列と、各非零要素の列番号を列挙した第２の配列と、第１の配列の中で行が変わる位置を列挙した第３の配列とが生成される。

ところで、大規模な行列演算は、行列を分割して複数のスレッドに割り振り、これら複数のスレッドを複数のプロセッサを用いて並列に実行することで、高速化することができる。このとき、行列の分割方法によっては、最終的な演算結果の中の同じ要素についての演算を、複数のスレッドが分担することになる場合がある。この場合、各スレッドに、中間の演算結果を記憶する記憶領域を割り当てる方法が考えられる。

例えば、圧縮列格納法で表現されたスパース行列と列ベクトルとの積を、複数のプロセッサを用いて計算する並列処理方法が提案されている。この並列処理方法では、行列に含まれる複数の列を均等に分割して複数のスレッドに割り振り、また、最終的な演算結果の列ベクトルと同じ大きさの記憶領域を各スレッドに割り当てる。そして、各スレッドが生成した中間の演算結果の列ベクトルを合算して、最終的な演算結果を求める。

なお、連立方程式の係数を表した係数行列では対角線付近と一部の正方領域内に非零要素が集中しやすいという特徴を利用して、係数行列を複数のブロックに分割し、複数のプロセッサを用いて並列計算を行う高速演算処理方法が提案されている。

特開２００９−１９９４３０号公報 国際公開第２００８／０２６２６１号

しかし、複数のスレッドそれぞれに対して中間の演算結果のベクトルを記憶する記憶領域を割り当てると、メモリの使用量が大きくなるという問題がある。例えば、１００万行×１００万行のスパース行列と列ベクトルとの積を、１０００スレッドを用いて並列計算する場合を考える。この場合、スパース行列のデータは圧縮格納法によって圧縮できる一方、上記の記憶領域は全体で１００万行の列ベクトル１０００個分の大きさになってしまう。この方法では、スレッド数の増加に伴ってメモリの使用量が増大する。

これに対し、複数のスレッドに共通の記憶領域を１つ用意し、各スレッドが演算結果のベクトルの中の要素に対して値を順次加算していく方法も考えられる。しかしながら、単純に記憶領域を共通化してしまうと、スレッド間で同じ要素（例えば、演算結果の列ベクトルの中の同じ行）に対する値の加算が同時に発生して、アクセスが競合する可能性がある。もし、アクセス競合に備えてスレッド間で排他制御を行うと、メモリアクセスのオーバヘッドが大きくなり並列処理の効率が低下するおそれがある。

１つの側面では、本発明は、行列演算においてメモリの記憶領域を効率的に利用できるプログラム、並列演算方法および情報処理装置を提供することを目的とする。

１つの態様では、複数のスレッドを並列に実行可能なコンピュータに、以下の処理を実行させるプログラムが提供される。零要素および非零要素を含む行列の中の第１の部分行列の演算を第１のスレッドに割り当て、行列の中の第２の部分行列の演算を第２のスレッドに割り当てる。第１の部分行列における行間または列間の非零要素の分布と、第２の部分行列における行間または列間の非零要素の分布とを比較する。比較の結果に応じて、第１および第２のスレッドが演算においてそれぞれ利用するベクトルを記憶する記憶領域の割り当てを変化させる。また、１つの態様では、複数のスレッドを並列に実行可能なコンピュータが行う並列演算方法が提供される。

また、１つの態様では、互いに並列にスレッドを実行可能な複数のプロセッサと、スレッドが演算においてそれぞれ利用するベクトルを記憶するメモリと、を有する情報処理装置が提供される。複数のプロセッサの１つは、以下の処理を実行する。零要素および非零要素を含む行列の中の第１の部分行列の演算を第１のスレッドに割り当て、行列の中の第２の部分行列の演算を第２のスレッドに割り当てる。第１の部分行列における行間または列間の非零要素の分布と、第２の部分行列における行間または列間の非零要素の分布とを比較する。比較の結果に応じて、第１および第２のスレッドに対するメモリにおけるベクトルを記憶する記憶領域の割り当てを変化させる。

１つの側面では、行列演算においてメモリの記憶領域を効率的に利用できる。

第１の実施の形態の情報処理装置を示す図である。 情報処理装置のハードウェア例を示すブロック図である。 スパース行列とベクトルの積の例を示す図である。 スパース行列の圧縮格納の例を示す図である。 スパース行列の分割とスレッド割り当ての例を示す図である。 スパース行列とベクトルの積の並列計算例を示す図である。 非零要素マップの例を示す図である。 非零要素マップの他の例を示す図である。 作業ベクトルの割り当て例を示す図である。 作業ベクトルの他の割り当て例を示す第１の図である。 作業ベクトルの他の割り当て例を示す第２の図である。 作業ベクトルの他の割り当て例を示す第３の図である。 作業ベクトルの他の割り当て例を示す第４の図である。 情報処理装置の機能例を示すブロック図である。 行列演算制御の手順例を示すフローチャートである。 非零要素チェックの手順例を示すフローチャートである。 作業ベクトル割り当ての手順例を示すフローチャートである。 並列行列演算の手順例を示すフローチャートである。

以下、本実施の形態を図面を参照して説明する。
［第１の実施の形態］
図１は、第１の実施の形態の情報処理装置を示す図である。

第１の実施の形態の情報処理装置１０は、行列演算（例えば、行列とベクトルとの積を求める演算）を、複数のスレッドを並列に動作させることで実行する。情報処理装置１０は、プロセッサ１１〜１３を含む複数のプロセッサとメモリ１４とを有する。

プロセッサ１１〜１３は、物理的に互いに並列にスレッドを実行することができる演算装置である。プロセッサ１１〜１３それぞれは、ＣＰＵ（Central Processing Unit）などの１つのプロセッサパッケージでもよいし、プロセッサパッケージに含まれるプロセッサコア（単にコアと呼ぶことがある）でもよい。例えば、プロセッサ１１〜１３は、メモリ１４に記憶されたプログラムを実行する。第１の実施の形態では、プロセッサ１２がスレッド２１を実行し、プロセッサ１３がスレッド２２を実行する。また、プロセッサ１１が、行列演算の並列化を制御するスレッドまたはプロセスを実行する。ただし、プロセッサ１２またはプロセッサ１３が、並列化制御を行うようにしてもよい。

メモリ１４は、プロセッサ１１〜１３からアクセスされる共有メモリであり、例えば、ＲＡＭ（Random Access Memory）である。メモリ１４は、スレッド２１，２２が演算において利用する１または２以上のベクトルを記憶する。このベクトルには、例えば、スレッド２１，２２それぞれの演算途中の値が書き込まれる。また、このベクトルには、例えば、スレッド２１，２２それぞれの演算が終了した時点の中間結果の値（最終結果が集計される前の値）が書き込まれる。行列演算が行列と列ベクトルの積を求める演算である場合、上記のベクトルは列ベクトルになる。また、行列演算が行ベクトルと行列の積を求める演算である場合、上記のベクトルは行ベクトルになる。第１の実施の形態では、以下に述べるように、スレッド２１，２２に対するメモリ１４のベクトルが記憶される記憶領域の割り当てを変化させる。例えば、情報処理装置１０は、スレッド２１，２２に対して共通の記憶領域（記憶領域２６）または異なる記憶領域（記憶領域２６，２７）を用意する。

ここで、プロセッサ１１は、零要素および非零要素を含む行列２３を分析し、行列２３に含まれる部分行列の演算をスレッド２１，２２に割り当てる。行列２３は、例えば、零要素の割合が大きく非零要素の割合が小さいスパース行列である。行列２３は、圧縮列格納法や圧縮行格納法などの圧縮格納法で表現されていてもよい。第１の実施の形態では、プロセッサ１１は、部分行列２４の演算をスレッド２１に割り当て、部分行列２５の演算をスレッド２２に割り当てる。例えば、部分行列２４と部分行列２５とを、列が重複しておらず少なくとも一部の行が重複したものとする。また、例えば、部分行列２４と部分行列２５とを、行が重複しておらず少なくとも一部の列が重複したものとする。

また、プロセッサ１１は、行列２３を分析し、演算において利用されるベクトルを記憶する記憶領域をスレッド２１，２２に割り当てる。このとき、プロセッサ１１は、部分行列２４の非零要素の分布と部分行列２５の非零要素の分布とを比較する。例えば、演算結果のベクトルが列ベクトルになる場合、プロセッサ１１は、列方向（複数の行の間）の非零要素の分布を比較する。また、例えば、演算結果のベクトルが行ベクトルになる場合、プロセッサ１１は、行方向（複数の列の間）の非零要素の分布を比較する。

そして、プロセッサ１１は、非零要素の分布の比較結果に応じて、スレッド２１，２２に対する記憶領域の割り当てを変化させる。例えば、プロセッサ１１は、部分行列２４と部分行列２５とが同じ行（または、同じ列）に非零要素を含まない場合、スレッド２１，２２に共通の記憶領域を割り当ててよいと判断する。これは、同じ行に非零要素が存在しなければ、スレッド２１，２２は、演算において利用される列ベクトルの中の同じ行要素にアクセスせず、アクセス競合が発生しないためである。また、同じ列に非零要素が存在しなければ、スレッド２１，２２は、演算において利用される行ベクトルの中の同じ列要素にアクセスせず、アクセス競合が発生しないためである。

例えば、スレッド２１，２２に共通の記憶領域２６が割り当てられた場合、スレッド２１は、部分行列２４の演算途中の値または演算終了時点の値（中間の演算結果）を記憶領域２６に書き込む。また、スレッド２２は、部分行列２５の演算途中の値または演算終了時点の値（中間の演算結果）を記憶領域２６に書き込む。このとき、スレッド２１とスレッド２２の間で、記憶領域２６へのアクセスの排他制御が行われなくてもよい。一方、スレッド２２にスレッド２１とは異なる記憶領域２７が割り当てられた場合、スレッド２１は、部分行列２４の演算についての値を記憶領域２６に書き込み、スレッド２２は、部分行列２５の演算についての値を記憶領域２７に書き込む。記憶領域２６，２７に記憶されたベクトルは、最終的な演算結果のベクトルを求めるときに合算される。

なお、スレッド２１，２２の間のアクセス競合は、特に次のような行列演算において問題となる。例えば、行列２３が圧縮列格納法で表現されており、行列２３と列ベクトルとの積を求める場合である。この場合、圧縮列格納法で表現された行列は列を分割する方が演算が効率的になるため、同じ行についてのスレッド２１，２２間の競合が問題となる。また、例えば、行列２３が圧縮行格納法で表現されており、行ベクトルと行列２３との積を求める場合である。この場合、圧縮行格納法で表現された行列は行を分割する方が演算が効率的になるため、同じ列についてのスレッド２１，２２間の競合が問題となる。

また、例えば、行列２３が対称行列であり、部分行列２４と対称の位置にある部分行列の演算もスレッド２１に割り当て、部分行列２５と対称の位置にある部分行列の演算もスレッド２２に割り当てる場合である。この場合、行列２３の下三角および上三角の何れか一方についてスレッド２１，２２間の競合が問題となる。行列２３と列ベクトルの積を求める場合、下三角の部分行列に含まれる同じ行についてスレッド２１，２２間の競合が問題となる。また、行ベクトルと行列２３の積を求める場合、上三角の部分行列に含まれる同じ列についてスレッド２１，２２間の競合が問題となる。

第１の実施の形態の情報処理装置１０によれば、部分行列２４の非零要素の分布と部分行列２５の非零要素の分布とが比較される。そして、比較結果に応じて、スレッド２１，２２への記憶領域の割り当てが変化する。よって、同じ記憶領域の同じ要素に対するアクセス競合を抑制でき（例えば、アクセス競合が発生しないようにでき）、スレッド２１，２２間の排他制御の負荷が軽減される。また、排他制御の負荷が重くならない範囲（例えば、排他制御が不要になる範囲）で、メモリ１４に確保する記憶領域の量を抑制できる。従って、行列演算においてメモリ１４の記憶領域を効率的に利用することができる。

［第２の実施の形態］
図２は、情報処理装置のハードウェア例を示すブロック図である。
第２の実施の形態の情報処理装置１００は、大規模行列演算が可能なコンピュータであり、例えば、ユーザからの要求に応じて行列演算を行うサーバコンピュータである。

情報処理装置１００は、ＣＰＵ１１０，１１０ａ，１１０ｂ，１１０ｃを含む複数のＣＰＵおよびＲＡＭ１２０を有する。複数のＣＰＵとＲＡＭ１２０とは、システムバス１３６に接続されている。また、情報処理装置１００は、ＨＤＤ（Hard Disk Drive）１３１、画像信号処理部１３２、入力信号処理部１３３、媒体リーダ１３４および通信インタフェース１３５を有する。ＨＤＤ１３１、画像信号処理部１３２、入力信号処理部１３３、媒体リーダ１３４および通信インタフェース１３５は、入出力バス１３７に接続されている。システムバス１３６と入出力バス１３７とは、例えば、ブリッジで接続されている。

ＣＰＵ１１０，１１０ａ，１１０ｂ，１１０ｃは、プログラムを実行するプロセッサパッケージである。ＣＰＵ１１０，１１０ａ，１１０ｂ，１１０ｃは、ＨＤＤ１３１からプログラムの命令やデータの少なくとも一部をＲＡＭ１２０にロードし、プログラムを実行する。各ＣＰＵは、複数のコアおよびキャッシュメモリを有する。

一例として、ＣＰＵ１１０は、コア１１１〜１１４を含む複数のコアおよびキャッシュメモリ１１５を有する。コア１１１〜１１４は、物理的に互いに並列にスレッドを実行することができる。キャッシュメモリ１１５は、ＲＡＭ１２０から読み込まれたプログラムの命令やデータを一時的に記憶する揮発性メモリであり、例えば、ＳＲＡＭ（Static Random Access Memory）である。キャッシュメモリはコア毎に設けてもよい。

ＲＡＭ１２０は、ＣＰＵ１１０，１１０ａ，１１０ｂ，１１０ｃから、高速なシステムバス１３６を介してアクセスされる共有メモリである。ＲＡＭ１２０は、プログラムの命令やデータを一時的に記憶する。なお、情報処理装置１００は、ＲＡＭ以外の種類の揮発性メモリを備えてもよいし、複数個のメモリを備えていてもよい。

ＨＤＤ１３１は、ＯＳ（Operating System）やアプリケーションソフトウェアなどのソフトウェアのプログラム、および、データを記憶する不揮発性の記憶装置である。なお、情報処理装置１００は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。

画像信号処理部１３２は、何れかのＣＰＵからの命令に従って、情報処理装置１００に接続されたディスプレイ１４１に画像を出力する。ディスプレイ１４１としては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ（ＰＤＰ：Plasma Display Panel）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなどを用いることができる。

入力信号処理部１３３は、情報処理装置１００に接続された入力デバイス１４２から入力信号を取得し、何れかのＣＰＵに出力する。入力デバイス１４２としては、マウスやタッチパネルやタッチパッドやトラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、情報処理装置１００に、複数の種類の入力デバイスが接続されてもよい。

媒体リーダ１３４は、記録媒体１４３に記録されたプログラムやデータを読み取る駆動装置である。記録媒体１４３として、例えば、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤなどの磁気ディスク、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）などの光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。媒体リーダ１３４は、例えば、記録媒体１４３から読み取ったプログラムの命令やデータをＲＡＭ１２０またはＨＤＤ１３１に格納する。

通信インタフェース１３５は、ネットワーク１４４に接続され、ネットワーク１４４を介して他の情報処理装置と通信するインタフェースである。通信インタフェース１３５は、ケーブルでスイッチやルータなどの通信装置と接続される有線通信インタフェースでもよいし、無線基地局と接続される無線通信インタフェースでもよい。

なお、情報処理装置１００は、媒体リーダ１３４を備えていなくてもよい。また、ユーザが操作する端末装置からネットワーク１４４経由で情報処理装置１００を制御できる場合、情報処理装置１００は、画像信号処理部１３２や入力信号処理部１３３を備えていなくてもよい。また、ディスプレイ１４１や入力デバイス１４２が、情報処理装置１００の筐体と一体に形成されていてもよい。なお、コア１１１〜１１３は前述のプロセッサ１１〜１３の一例であり、ＲＡＭ１２０は前述のメモリ１４の一例である。

次に、第２の実施の形態で行う行列演算について説明する。
図３は、スパース行列とベクトルの積の例を示す図である。
情報処理装置１００は、行列と入力ベクトルとの積を反復的に計算する。行列は、例えば、連立方程式の係数を表した係数行列である。情報処理装置１００は、１回目に、行列と初期の入力ベクトルとの積を計算する。情報処理装置１００は、１回目の積である結果ベクトルを所定のアルゴリズムに従って加工し、次の入力ベクトルとして使用する。情報処理装置１００は、２回目に、１回目と同じ行列と１回目の結果ベクトルを加工して得られた入力ベクトルとの積を計算する。以上の行列演算が、所定の終了条件（例えば、反復回数や結果ベクトルに含まれる値の精度などの条件）を満たすまで繰り返される。

第２の実施の形態では、行列として対称行列であり且つスパース行列であるもの（対称スパース行列）を想定し、入力ベクトルとして列ベクトルを想定する。対称スパース行列と列ベクトルとの積である結果ベクトルは、列ベクトルになる。例えば、６行×６列の対称スパース行列の１行目が（４．０，９．０，０，３．０，０，０）であり、入力ベクトルが（０．１，０．２，０．３，０．４，０．５，０．６）であるとする。この場合、結果ベクトルの１行目の値が３．４と算出される。なお、以下では説明を簡単にするために６行×６列の行列の例を用いることがあるが、第２の実施の形態では次数（１辺の要素数）が数万〜数千万程度の大規模な対称スパース行列が使用され得る。

図４は、スパース行列の圧縮格納の例を示す図である。
第２の実施の形態では、対称スパース行列の下三角の領域が圧縮列格納法によって表現される。対称スパース行列の上三角の領域（対角要素を除く）は、下三角の領域に基づいて再現することができるため行列データに含めず省略することができる。

行列データには、要素配列１２１（Ｖａｌ）、行番号配列１２２（Ｒｏｗ）および列ポインタ配列１２３（Ｃｐ）が含まれる。要素配列１２１は、対称スパース行列の下三角の領域に含まれる非零要素の値を列優先の順序で並べたものである。行番号配列１２２は、要素配列１２１に列挙された各非零要素の行番号を記載したものである。行番号配列１２２のｋ番目の値（Ｒｏｗ（ｋ））は、ｋ番目の非零要素の行番号を表す。要素配列１２１と行番号配列１２２の長さは、対称スパース行列に含まれる非零要素の数になる。列ポインタ配列１２３は、要素配列１２１の中で列が変わる位置を列挙したものである。列ポインタ配列１２３のｋ番目の値（Ｃｐ（ｋ））は、ｋ列目の先頭の非零要素の番号を表す。ただし、列ポインタ配列１２３の長さは、対称スパース行列の列数より１だけ大きい。列ポインタ配列１２３の末尾には、要素番号の最大値より１だけ大きい数が格納される。

例えば、６行×６列の対称スパース行列の下三角の領域に、次のような要素が含まれているとする。１列目：（４．０，９．０，０，３．０，０，０）、２列目：（１１．０，５．０，０，０，０）、３列目：（６．０，０，０，０）、４列目：（１．０，８．０，１２．０）、５列目：（２．０，１０．０）、６列目：（７．０）。

この場合、要素配列１２１は、（４．０，９．０，３．０，１１．０，５．０，６．０，１．０，８．０，１２．０，２．０，１０．０，７．０）となる。行番号配列１２２は、（１，２，４，２，３，３，４，５，６，５，６，６）となる。列ポインタ配列１２３は、（１，４，６，７，１０，１２，１３）となる。例えば、４列目の上から２番目の非零要素は、要素配列１２１、行番号配列１２２および列ポインタ配列１２３を参照して次のように特定される。まず、Ｃｐ（４）＝７から、４列目の先頭の非零要素の番号が７と特定され、４列目の２番目の非零要素の番号は８と特定される。そして、Ｖａｌ（８）＝８．０，Ｒｏｗ（８）＝５から、４列目の上から２番目の非零要素は５行目にあり、その値は８．０であると特定される。

なお、対称スパース行列の上三角の領域（対角要素を含む）は、要素配列１２１、行番号配列１２２および列ポインタ配列１２３を、値を変換せずに読み替えることで再現できる。すなわち、要素配列１２１を行優先の順序で非零要素を列挙したものと読み替え、行番号配列１２２を列番号配列に読み替え、列ポインタ配列１２３を行ポインタ配列と読み替える。このように読み替えた要素配列、列番号配列および列ポインタ配列は、上三角の領域（対角要素を含む）を圧縮行格納法で表現したものに相当する。

図５は、スパース行列の分割とスレッド割り当ての例を示す図である。
情報処理装置１００は、対称スパース行列を複数の部分行列に分割して複数のスレッドに割り振り、複数のコアを用いてそれらスレッドを並列に実行する。第２の実施の形態では、情報処理装置１００は、対称スパース行列の下三角の領域について、各スレッドに連続する１または２以上の列を割り当てる。また、情報処理装置１００は、対称スパース行列の上三角の領域（対角要素を除く）について、各スレッドに連続する１または２以上の行を割り当てる。このとき、対称の位置にある下三角の領域のｊ列目と上三角の領域（対角要素を除く）のｊ行目とを、同じスレッドに割り当てるようにする。

対称スパース行列は圧縮列格納法で表現されているため、情報処理装置１００は、まず下三角の領域の列とスレッドとの関係を決定する。これに伴い、上三角の領域（対角要素を除く）の行とスレッドとの関係も自動的に決まる。このとき、情報処理装置１００は、各スレッドで処理される非零要素の数ができる限り均等になるようにする。例えば、要素配列１２１の長さが１２であり、並列実行可能なスレッドが４つであるとする。この場合、情報処理装置１００は、各スレッドが担当する非零要素が１２／４＝３個に近くなるように、要素配列１２１を分割する。図５の例では、下三角の領域の１列目がスレッド＃１に割り当てられ、２列目および３列目がスレッド＃２に割り当てられ、４列目がスレッド＃３に割り当てられ、５列目および６列目がスレッド＃４に割り当てられている。

対称スパース行列が分割されると、情報処理装置１００は、スレッドポインタ配列１２４（Ｂｐ）を生成する。スレッドポインタ配列１２４は、各スレッドに割り当てられた列集合の先頭の列を示す列番号が列挙される。複数のスレッドには、連番のスレッド番号が付与されている。スレッドポインタ配列１２４のｋ番目の値（Ｂｐ（ｋ））は、スレッド＃ｋに割り当てられた列集合の先頭の列を表している。ただし、スレッドポインタ配列１２４の長さは、スレッド数より１だけ大きい。スレッドポインタ配列１２４の末尾には、対称スパース行列の列数より１だけ大きい数が格納される。

図６は、スパース行列とベクトルの積の並列計算例を示す図である。
ＲＡＭ１２０には、要素配列１２１、行番号配列１２２、列ポインタ配列１２３およびスレッドポインタ配列１２４が記憶される。更に、ＲＡＭ１２０には、作業領域１２７（Ｗｏｒｋ）、入力ベクトル１２８（Ｘ）および結果ベクトル１２９（Ｙ）が記憶される。

作業領域１２７は、中間の演算結果を記憶する１または２以上の列ベクトルとしての作業ベクトルを含む。下三角の領域の部分行列について各スレッドが計算した値は、作業領域１２７に書き込まれる。一方、上三角の領域（対角要素を除く）について各スレッドが計算した値は、結果ベクトル１２９に直接書き込まれる。１または２以上の作業ベクトルは、複数のスレッドの実行が完了した後、結果ベクトル１２９に合算される。ここでは、スレッド毎に１つの作業ベクトルが作業領域１２７に確保される場合を考える。

下三角の領域について、スレッド＃１は、１行１列の非零要素と入力ベクトル１２８の１行目との積を、スレッド＃１に割り当てられた作業領域１２７の作業ベクトルの１行目（１行１列）に加算する。同様に、スレッド＃１は、２行１列の非零要素と入力ベクトル１２８の１行目との積を作業領域１２７の２行１列に加算し、４行１列の非零要素と入力ベクトル１２８の１行目との積を作業領域１２７の４行１列に加算する。上三角の領域（対角要素を除く）について、スレッド＃１は、１行２列の非零要素と入力ベクトル１２８の２行目との積を結果ベクトル１２９の１行目に加算し、１行４列の非零要素と入力ベクトル１２８の４行目との積を結果ベクトル１２９の１行目に加算する。

下三角の領域について、スレッド＃２は、２行２列の非零要素と入力ベクトル１２８の２行目との積を作業領域１２７の２行２列に加算する。また、スレッド＃２は、３行２列の非零要素と入力ベクトル１２８の２行目との積を作業領域１２７の３行２列に加算し、３行３列の非零要素と入力ベクトル１２８の３行目との積を作業領域１２７の３行２列に加算する。上三角の領域（対角要素を除く）について、スレッド＃２は、２行３列の非零要素と入力ベクトル１２８の３行目との積を結果ベクトル１２９の２行目に加算する。

同様に、スレッド＃３は、作業領域１２７の４行３列、５行３列および６行３列に値を加算し、結果ベクトル１２９の４行目に値を加算する。スレッド＃４は、作業領域１２７の５行４列および６行４列に値を加算し、結果ベクトル１２９の５行目に値を加算する。スレッド＃１〜＃４の実行が完了すると、情報処理装置１００は、作業領域１２７に含まれる４つの作業ベクトルを結果ベクトル１２９に足し合わせる。これにより、結果ベクトル１２９が、対称スパース行列と入力ベクトル１２８の積を表すことになる。

上記では、異なるスレッドに対して異なる作業ベクトルを割り当てることとした。しかし、この方法では、作業領域１２７内の多くの要素が、値が加算されずに０のままになる。例えば、図６の例では、作業領域１２７の４行２列、５行２列および６行２列の値が０のままである。この４行２列、５行２列および６行２列の領域をスレッド＃３が使用したとしても、スレッド＃２とスレッド＃３の間でアクセス競合は生じない。一方、各作業ベクトルは結果ベクトル１２９に合算されるものであるため、値の加算される行が正しければ値の加算される列が変わっても、最終的な演算結果に影響しない。そこで、第２の実施の形態では、情報処理装置１００は、アクセス競合が生じない範囲で、２以上のスレッドに割り当てる作業領域１２７の作業ベクトルを共通化するようにする。

図７は、非零要素マップの例を示す図である。
情報処理装置１００は、対称スパース行列を分析して、確保する作業ベクトルの数とスレッドへの作業ベクトルの割り当てを決定する。対称スパース行列の分析において、情報処理装置１００は、対称スパース行列の行を複数の区間に分割する。好ましくは、区間１つ当たりの行数をできる限り均等にする。情報処理装置１００は、区間の数（分割数）を予め固定で決めておいてもよいし、対称スパース行列の次数に応じて分割数を変えてもよい。例えば、数万次元程度の対称スパース行列に対して分割数を１００とする。

そして、情報処理装置１００は、対称スパース行列の下三角の領域について、行方向および列方向に細分化されたブロック毎に非零要素の有無を確認し、非零要素の分布を示す非零要素マップ１２５（Ｍａｐ）を生成する。非零要素マップ１２５の行は上記の区間に対応し、非零要素マップ１２５の列はスレッド（すなわち、当該スレッドに割り当てられた対称スパース行列の列集合）に対応する。非零要素マップ１２５では、各ブロックの状態が１ビットのフラグで表現される。フラグ＝１は当該ブロック内に少なくとも１つの非零要素があることを示し、フラグ＝０は当該ブロック内に非零要素がないことを示す。

例えば、図５のように分割した６行×６列の対称スパース行列の下三角の領域から、６行×４列の非零要素マップ１２５が生成される。ここでは、非零要素マップ１２５の１つの行が対称スパース行列の１つの行に対応する。例えば、非零要素マップ１２５の１行１列、２行１列、４行１列、２行２列、３行２列、４行３列、５行３列、６行３列、５行４列および６行４列がフラグ＝１であり、他の要素がフラグ＝０になる。このような非零要素マップ１２５の列同士を比較することで、共通の作業ベクトルを割り当てることが可能なスレッドの組み合わせを探索することができる。

図８は、非零要素マップの他の例を示す図である。
ここでは、２５３４行×２５３４列の対称スパース行列を８スレッドで並列処理する場合を考える。行の分割数を８とすると、１区間当たりの行数は３１７行（２５３４／８の小数点以下を繰り上げた値）になる。ただし、端数処理の影響で末尾の区間の行数は３１５行になる。情報処理装置１００は、８行×８列の非零要素マップ１２５を生成し、各ブロックに少なくとも１つの非零要素が含まれるか確認する。大規模な対称スパース行列では、非零要素が対角線付近に集中することがある。その場合、図８に示すように、非零要素マップ１２５の対角要素および対角要素に隣接する一部の要素がフラグ＝１になり、非零要素マップ１２５の他の多くの要素がフラグ＝０になり得る。

図９は、作業ベクトルの割り当て例を示す図である。
上記の通り、情報処理装置１００は、非零要素マップ１２５に基づいて、共通の作業ベクトルを割り当て可能なスレッドの組み合わせを探索する。具体的には、情報処理装置１００は、非零要素マップ１２５から、フラグ＝１が同じ行で衝突していない列の組み合わせを探索する。図７の非零要素マップ１２５の例の場合、１列目と４列目との組み合わせはフラグ＝１が衝突せず、２列目と３列目との組み合わせはフラグ＝１が衝突しない。

この場合、例えば、情報処理装置１００は、スレッド＃１，＃４に作業ベクトル１（作業領域１２７の１列目）を割り当て、スレッド＃２，＃３に作業ベクトル２（作業領域１２７の２列目）を割り当てる。作業領域１２７には、４つのスレッド（スレッド＃１〜＃４）に対して２つの作業ベクトルを確保すればよいことになる。各スレッドに作業ベクトルを割り当てると、情報処理装置１００は、スレッドが使用する作業ベクトルの番号を列挙した作業ポインタ配列１２６（Ｕｐ）を生成する。作業ポインタ配列１２６のｋ番目の値（Ｕｐ（ｋ））は、スレッド＃ｋが使用する作業ベクトルを表す。

図１０は、作業ベクトルの他の割り当て例を示す第１の図である。
ここでは、対称スパース行列を８スレッドで並列処理し、対称スパース行列の行を８区間に分割した場合を考える。スレッド＃１の区間１，２、スレッド＃２の区間２，３、スレッド＃３の区間３，４、スレッド＃４の区間４，５、スレッド＃５の区間５，６、スレッド＃６の区間６，７、スレッド＃７の区間７，８およびスレッド＃８の区間８に非零要素が存在する。他の区間には非零要素は存在しない。

この場合、スレッド＃１，＃３，＃５，＃７に共通の作業ベクトルを割り当てても、これらのスレッド間で作業ベクトルの同じ要素に対するアクセス競合は生じない。また、スレッド＃２，＃４，＃６，＃８に共通の作業ベクトルを割り当てても、これらのスレッド間で作業ベクトルの同じ要素に対するアクセス競合は生じない。そこで、例えば、情報処理装置１００は、スレッド＃１，＃３，＃５，＃７に作業ベクトル１を割り当て、スレッド＃２，＃４，＃６，＃８に作業ベクトル２を割り当てる。８個のスレッド（スレッド＃１〜＃８）に対して２つの作業ベクトルを確保すればよいことになる。

このように、３以上のスレッドに共通の作業ベクトルを割り当てることもできる。多くの大規模な対称スパース行列では、スレッドの組み合わせを適切に判断することで、作業領域１２７に確保する作業ベクトルの数を２〜３個に抑えることが可能である。

図１１は、作業ベクトルの他の割り当て例を示す第２の図である。
ここでは、図１０と同様に、対称スパース行列を８スレッドで並列処理し、対称スパース行列の行を８区間に分割した場合を考える。スレッド＃１の区間１，２，４、スレッド＃２の区間２，３，６、スレッド＃３の区間３，４，８、スレッド＃４の区間４，５、スレッド＃５の区間５，６、スレッド＃６の区間６，７、スレッド＃７の区間７，８およびスレッド＃８の区間８に非零要素が存在する。他の区間には非零要素は存在しない。

この場合、例えば、情報処理装置１００は、スレッド＃１，＃５，＃７に作業ベクトル１を割り当て、スレッド＃２，＃４，＃８に作業ベクトル２を割り当て、スレッド＃３，＃６に作業ベクトル３を割り当てる。８個のスレッド（スレッド＃１〜＃８）に対して３つの作業ベクトルを作業領域１２７に確保すればよいことになる。

図１２は、作業ベクトルの他の割り当て例を示す第３の図である。
ここでは、図１０と同様に、対称スパース行列を８スレッドで並列処理し、対称スパース行列の行を８区間に分割した場合を考える。スレッド＃１の区間１，３，７、スレッド＃２の区間２，６、スレッド＃３の区間３，５、スレッド＃４の区間４，８、スレッド＃５の区間５，７、スレッド＃６の区間６、スレッド＃７の区間７およびスレッド＃８の区間８に非零要素が存在する。他の区間には非零要素は存在しない。

この場合、例えば、情報処理装置１００は、スレッド＃１，＃２，＃４に作業ベクトル１を割り当て、スレッド＃３，＃６，＃７，＃８に作業ベクトル２を割り当て、スレッド＃５に作業ベクトル３を割り当てる。８個のスレッド（スレッド＃１〜＃８）に対して３つの作業ベクトルを作業領域１２７に確保すればよいことになる。

ここで、アクセス競合が生じないスレッドの組み合わせは、様々なアルゴリズムを使用して探索することができる。図９〜１２の組み合わせ例は、非零要素マップ１２５の左側の列から順に、既に確保した作業ベクトルの中からフラグ＝１が衝突しないものを探し、該当する作業ベクトルがなければ新たな作業ベクトルを確保することで探索できる。

例えば、図１２の場合、情報処理装置１００は、作業ベクトル１を確保し、スレッド＃１に作業ベクトル１を割り当てる。次に、情報処理装置１００は、作業ベクトル１にスレッド＃２を追加可能であるため、スレッド＃２に作業ベクトル１を割り当てる。次に、情報処理装置１００は、作業ベクトル１にスレッド＃３を追加できないため、作業ベクトル２を確保し、スレッド＃３に作業ベクトル２を割り当てる。

次に、情報処理装置１００は、作業ベクトル１にスレッド＃４を追加可能であるため、スレッド＃４に作業ベクトル１を割り当てる。次に、情報処理装置１００は、作業ベクトル１，２の何れにもスレッド＃５を追加できないため、作業ベクトル３を確保し、スレッド＃５に作業ベクトル３を割り当てる。次に、情報処理装置１００は、作業ベクトル１にはスレッド＃６を追加できないが、作業ベクトル２にスレッド＃６を追加可能であるため、スレッド＃６に作業ベクトル２を割り当てる。同様にして、情報処理装置１００は、スレッド＃７，＃８に作業ベクトル２を割り当てる。

ただし、対称スパース行列の中の非零要素の分布が複雑である場合、使用する探索アルゴリズムによって、作業領域１２７に確保される作業ベクトルの数が変わることがある。
図１３は、作業ベクトルの他の割り当て例を示す第４の図である。

ここでは、対称スパース行列を６スレッドで並列処理し、対称スパース行列の行を８区間に分割した場合を考える。スレッド＃１の区間１，３，８、スレッド＃２の区間２，５、スレッド＃３の区間３，４，７，８、スレッド＃４の区間４，５，７、スレッド＃５の区間６およびスレッド＃６の区間６に非零要素が存在する。

ある探索アルゴリズム（例えば、前述のアルゴリズム）によれば、例えば、情報処理装置１００は、スレッド＃１，＃２，＃５に作業ベクトル１を割り当て、スレッド＃３，＃６に作業ベクトル２を割り当て、スレッド＃４に作業ベクトル３を割り当てる。すなわち、３つの作業ベクトルが作業領域１２７に確保される。一方、別の探索アルゴリズムによれば、例えば、情報処理装置１００は、スレッド＃１，＃４，＃５に作業ベクトル１を割り当て、スレッド＃２，＃３，＃６に作業ベクトル２を割り当てる。すなわち、先の探索アルゴリズムの場合よりも少ない２つの作業ベクトルが、作業領域１２７に確保される。

通常、作業ベクトルの数を最小化できる探索アルゴリズム（最適解を求めるアルゴリズム）は、その他の探索アルゴリズム（準最適解を求めるアルゴリズム）よりも計算量が大きくなる。情報処理装置１００は、解の精度と計算量とのバランスを考慮して、使用する探索アルゴリズムを選択するようにしてもよい。以下の説明では、情報処理装置１００は、少ない計算量で準最適解を求める探索アルゴリズムを使用するものとする。

次に、情報処理装置１００の機能および処理手順について説明する。
図１４は、情報処理装置の機能例を示すブロック図である。
情報処理装置１００は、データ記憶部１５０、行列演算要求部１６１、並列化制御部１６２、並列処理部１６５およびＯＳ１６８を有する。データ記憶部１５０は、ＲＡＭ１２０に確保した記憶領域として実現される。行列演算要求部１６１、並列化制御部１６２および並列処理部１６５は、ソフトウェアのモジュールとして実現される。特に、並列化制御部１６２および並列処理部１６５は、数値計算ライブラリであってもよい。

データ記憶部１５０は、行列記憶部１５１、制御データ記憶部１５２、中間データ記憶部１５３およびベクトル記憶部１５４を有する。行列記憶部１５１は、対称スパース行列を表す行列データを記憶する。行列データは、要素配列１２１、行番号配列１２２および列ポインタ配列１２３を含む。制御データ記憶部１５２は、並列化の制御に用いる制御データを記憶する。制御データは、スレッドポインタ配列１２４、非零要素マップ１２５および作業ポインタ配列１２６を含む。中間データ記憶部１５３は、作業領域１２７を含む。ベクトル記憶部１５４は、入力ベクトル１２８および結果ベクトル１２９を記憶する。

行列演算要求部１６１は、要素配列１２１、行番号配列１２２、列ポインタ配列１２３および入力ベクトル１２８をデータ記憶部１５０に格納し、対称スパース行列と入力ベクトル１２８との積の計算を並列化制御部１６２に要求する。結果ベクトル１２９が得られると、行列演算要求部１６１は、結果ベクトル１２９を加工し、次の入力ベクトル１２８として使用する。行列演算要求部１６１は、反復回数や結果ベクトル１２９に含まれる値の精度などの演算状況が所定の終了条件を満たすまで、対称スパース行列と入力ベクトル１２８との積の計算を繰り返し並列化制御部１６２に対して要求する。

並列化制御部１６２は、行列演算の並列化を制御する。並列化制御部１６２は、行列解析部１６３およびベクトル入出力部１６４を有する。
行列解析部１６３は、行列演算要求部１６１から最初に対称スパース行列が指定されたとき（反復演算の１回目のとき）、対称スパース行列を分析して並列化方法を決定する。行列解析部１６３は、対称スパース行列を分割して複数のスレッドに割り振る。行列演算を行うスレッドの数は、例えば、情報処理装置１００が備えるハードウェア資源の量、情報処理装置１００の現在の負荷、ユーザの契約内容などの条件に基づいて決定される。また、行列解析部１６３は、作業領域１２７に確保する作業ベクトルの数を決定し、複数のスレッドそれぞれに何れかの作業ベクトルを割り当てる。

ベクトル入出力部１６４は、行列演算要求部１６１から入力ベクトル１２８が指定される毎に（反復演算の１回毎に）、作業領域１２７および結果ベクトル１２９を初期化する。また、ベクトル入出力部１６４は、複数のスレッドの行列演算が完了すると、作業領域１２７に含まれる全ての作業ベクトルを結果ベクトル１２９に足し合わせて、対称スパース行列と入力ベクトル１２８の積としての最終的な解を求める。

並列処理部１６５は、並列に実行される複数のスレッドの機能を実現する。並列処理部１６５は、非零要素チェック部１６６および行列演算部１６７を有する。
非零要素チェック部１６６は、行列解析部１６３からの要求に応じて、自スレッドに割り当てられた部分行列における非零要素の分布を確認し、非零要素マップ１２５の自スレッドに対応する列のフラグを更新していく。非零要素マップ１２５の生成は、複数のスレッドを用いて並列化されることになる。なお、各スレッドは、スレッドポインタ配列１２４を参照して、割り当てられた対称スパース行列の列を特定できる。

行列演算部１６７は、ベクトル入出力部１６４からの要求に応じて、自スレッドに割り当てられた部分行列と入力ベクトル１２８との積を計算する。対称スパース行列の下三角の領域については、行列演算部１６７は、行列解析部１６３から割り当てられた作業ベクトルに演算結果を書き込む。一方、対称スパース行列の上三角の領域（対角要素を除く）については、行列演算部１６７は、結果ベクトル１２９に演算結果を書き込む。

ＯＳ１６８は、行列解析部１６３からの要求に応じて、複数のスレッドを起動し、ＣＰＵ１１０，１１０ａ，１１０ｂ，１１０ｃに含まれる複数のコアにそれら複数のスレッドを割り振る。原則として、異なるスレッドは異なるコアに割り当てられる。また、ＯＳ１６８は、行列解析部１６３からの要求に応じて、ＲＡＭ１２０に作業領域１２７を確保する。なお、並列化制御部１６２、並列処理部１６５およびＯＳ１６８に相当するスレッドまたはプロセスは、好ましくは、互いに異なるコアで実行される。

図１５は、行列演算制御の手順例を示すフローチャートである。
（Ｓ１）行列演算要求部１６１は、対称スパース行列を表す要素配列１２１、行番号配列１２２および列ポインタ配列１２３を行列記憶部１５１に格納する。行列解析部１６３は、行列記憶部１５１から行番号配列１２２および列ポインタ配列１２３を読み込む。

（Ｓ２）行列解析部１６３は、行列演算に使用するスレッドの数を決定し、決定した数のスレッドを起動するようＯＳ１６８に要求する。このとき、行列解析部１６３は、各スレッドに連番のスレッド番号を付与する。また、行列解析部１６３は、要素配列１２１を参照して対称スパース行列に含まれる非零要素の数を特定し、非零要素の数ができる限り均等になるように対称スパース行列の列を分割する。そして、行列解析部１６３は、スレッドポインタ配列１２４を生成して、制御データ記憶部１５２に格納する。

（Ｓ３）行列解析部１６３は、非零要素マップ１２５を生成して、制御データ記憶部１５２に格納する。そして、行列解析部１６３は、スレッド毎に非零要素チェック部１６６を呼び出す。各スレッドの非零要素チェック部１６６は、非零要素マップ１２５に含まれる自スレッドに対応する列のフラグを更新する。非零要素チェックの詳細は後述する。

（Ｓ４）行列解析部１６３は、ステップＳ３で生成された非零要素マップ１２５に基づいて、作業領域１２７に確保する作業ベクトルの数および各スレッドへの作業ベクトルの割り当てを決定する。そして、行列解析部１６３は、作業ポインタ配列１２６を生成し、制御データ記憶部１５２に格納する。作業ベクトル割り当ての詳細は後述する。

（Ｓ５）行列解析部１６３は、ステップＳ４で決定した数の作業ベクトルを含む作業領域１２７を中間データ記憶部１５３に確保するよう、ＯＳ１６８に要求する。また、行列解析部１６３は、空の結果ベクトル１２９を生成し、ベクトル記憶部１５４に格納する。

（Ｓ６）行列演算要求部１６１は、入力ベクトル１２８をベクトル記憶部１５４に格納する。ただし、反復演算の１回目に使用する初期の入力ベクトル１２８については、ステップＳ１でベクトル記憶部１５４に格納されてもよい。ベクトル入出力部１６４は、ベクトル記憶部１５４から入力ベクトル１２８を読み込む。

（Ｓ７）ベクトル入出力部１６４は、スレッド毎に行列演算部１６７を呼び出す。行列演算部１６７は、自スレッドに割り当てられた部分行列と入力ベクトル１２８との積を計算する。ベクトル入出力部１６４は、作業領域１２７に含まれる全ての作業ベクトルを結果ベクトル１２９に足し合わせる。これにより、結果ベクトル１２９は、対称スパース行列と入力ベクトル１２８との積を表す。並列行列演算の詳細は後述する。

（Ｓ８）行列演算要求部１６１は、結果ベクトル１２９をベクトル記憶部１５４から読み込む。そして、行列演算要求部１６１は、演算状況が所定の終了条件を満たしているか判断する。終了条件を満たす場合、行列演算要求部１６１は、反復演算を打ち切る。終了条件を満たさない場合、行列演算要求部１６１は、読み込んだ結果ベクトル１２９を用いて次の入力ベクトル１２８を生成し、ステップＳ６に処理を進める。

図１６は、非零要素チェックの手順例を示すフローチャートである。この非零要素チェックは、図１５のフローチャートのステップＳ３で実行される。
（Ｓ３０）行列解析部１６３は、対称スパース行列の行を複数の区間に分割する。例えば、区間の数（分割数）が予め決まっているとすると、行列解析部１６３は、幅（区間１つ当たりの行数）を、幅ｗ＝（対称スパース行列の行数＋分割数−１）／分割数と計算する。なお、以下の説明では除算は小数点以下を切り捨てるものとする。

（Ｓ３１）行列解析部１６３は、大きさが分割数×スレッド数の非零要素マップ１２５を生成し、制御データ記憶部１５２に格納する。このとき、行列解析部１６３は、非零要素マップ１２５の全ての要素を零に初期化する。そして、行列解析部１６３は、スレッド毎に非零要素チェック部１６６を呼び出す。以下のステップＳ３２〜Ｓ３８の処理が、複数のスレッドで並列に実行される。以下では、スレッド番号＝ｔのスレッド（スレッド＃ｔ）がステップＳ３２〜Ｓ３８の処理を行う場合を説明する。

（Ｓ３２）非零要素チェック部１６６は、スレッドポインタ配列１２４を参照して、スレッド＃ｔに割り当てられた列集合の中の先頭の列を選択する。具体的には、非零要素チェック部１６６は、列番号ｃ＝Ｂｐ（ｔ）を特定する。

（Ｓ３３）非零要素チェック部１６６は、列ポインタ配列１２３を参照して、ステップＳ３２またはステップＳ３８で選択した列に含まれる先頭の非零要素を選択する。具体的には、非零要素チェック部１６６は、要素番号ｅ＝Ｃｐ（ｃ）を特定する。

（Ｓ３４）非零要素チェック部１６６は、ステップＳ３３またはステップＳ３６で選択した非零要素に対応する非零要素マップ１２５のフラグを１に設定する。具体的には、非零要素チェック部１６６は、Ｍａｐ（（Ｒｏｗ（ｅ）−１）／ｗ＋１，ｔ）＝１とする。

（Ｓ３５）非零要素チェック部１６６は、列ポインタ配列１２３を参照して、ステップＳ３２またはステップＳ３８で選択した列に含まれる全ての非零要素を選択したか判断する。具体的には、非零要素チェック部１６６は、ステップＳ３３またはステップＳ３６で特定した要素番号ｅがＣｐ（ｃ＋１）−１に一致するか判断する。全て選択した場合はステップＳ３７に処理を進め、未選択のものがある場合はステップＳ３６に処理を進める。

（Ｓ３６）非零要素チェック部１６６は、次の非零要素を選択する。具体的には、非零要素チェック部１６６は、要素番号ｅをインクリメント（現在の要素番号ｅに１を加算）する。そして、ステップＳ３４に処理を進める。

（Ｓ３７）非零要素チェック部１６６は、スレッドポインタ配列１２４を参照して、スレッド＃ｔに割り当てられた列を全て選択したか判断する。具体的には、非零要素チェック部１６６は、列番号ｃがＢｐ（ｔ＋１）−１に一致するか判断する。全て選択した場合、非零要素チェック部１６６は、非要素チェックを終了して行列解析部１６３に完了を通知する。未選択の列がある場合、ステップＳ３８に処理を進める。

（Ｓ３８）非零要素チェック部１６６は、次の列を選択する。具体的には、非零要素チェック部１６６は、列番号ｃをインクリメント（現在の列番号ｃに１を加算）する。そして、ステップＳ３３に処理を進める。

図１７は、作業ベクトル割り当ての手順例を示すフローチャートである。この作業ベクトル割り当ては、図１５のフローチャートのステップＳ４で実行される。
（Ｓ４０）行列解析部１６３は、作業ベクトルを１つ用意すると決定し、用意した作業ベクトルをスレッド＃１に割り当てる。具体的には、行列解析部１６３は、変数としてベクトル数ｎを１に設定すると共に、Ｕｐ（１）＝１に設定する。

（Ｓ４１）行列解析部１６３は、スレッド＃２を選択する。具体的には、行列解析部１６３は、変数としてスレッド番号ｔを２に設定する。
（Ｓ４２）行列解析部１６３は、現在までに用意した作業ベクトルのうちの先頭を選択する。具体的には、行列解析部１６３は、変数としてベクトル番号ｖを１に設定する。

（Ｓ４３）行列解析部１６３は、非零要素マップ１２５のｖ列目とｔ列目とを比較し、フラグ＝１の分布に重複があるか判断する。具体的には、行列解析部１６３は、ｖ列目の列ベクトルとｔ列目の列ベクトルとの論理積を計算し、ｖ列目とｔ列目の両方でフラグ＝１になっている区間が存在するか判断する。フラグ＝１の分布に重複がある場合はステップＳ４６に処理を進め、重複がない場合はステップＳ４４に処理を進める。

（Ｓ４４）行列解析部１６３は、ステップＳ４１またはステップＳ５２で選択したスレッドに、ステップＳ４２またはステップＳ４７で選択した作業ベクトルを割り当てる。具体的には、行列解析部１６３は、Ｕｐ（ｔ）＝ｖに設定する。

（Ｓ４５）行列解析部１６３は、非零要素マップ１２５のｔ列目のフラグ＝１をｖ列目にコピーする。そして、ステップＳ５１に処理を進める。
（Ｓ４６）行列解析部１６３は、現在までに用意した作業ベクトルを全て選択したか判断する。具体的には、行列解析部１６３は、ベクトル番号ｖがベクトル数ｎに一致するか判断する。全ての作業ベクトルを選択した場合はステップＳ４８に処理を進め、未選択の作業ベクトルがある場合はステップＳ４７に処理を進める。

（Ｓ４７）行列解析部１６３は、現在用意している次の作業ベクトルを選択する。具体的には、行列解析部１６３は、ベクトル番号ｖをインクリメント（現在のベクトル番号ｖに１を加算）する。そして、ステップＳ４３に処理を進める。

（Ｓ４８）行列解析部１６３は、作業ベクトルを１つ追加する。具体的には、行列解析部１６３は、ベクトル数ｎをインクリメント（現在のベクトル数ｎに１を加算）する。
（Ｓ４９）行列解析部１６３は、ステップＳ４１またはステップＳ５２で選択したスレッドに、ステップＳ４８で新たに追加した作業ベクトルを割り当てる。具体的には、行列解析部１６３は、Ｕｐ（ｔ）＝ｎに設定する。

（Ｓ５０）行列解析部１６３は、非零要素マップ１２５のｎ列目のフラグをｔ列目のフラグで上書きする。すなわち、ｎ列目の列ベクトルをｔ列目の列ベクトルと一致させる。
（Ｓ５１）行列解析部１６３は、スレッドを全て選択したか判断する。具体的には、行列解析部１６３は、スレッド番号ｔがスレッド数に一致するか判断する。全てのスレッドを選択した場合、行列解析部１６３は、作業ベクトル割り当てを終了する。これにより、作業領域１２７に確保する作業ベクトルの数と、作業ポインタ配列１２６の内容が確定する。未選択のスレッドがある場合、ステップＳ５２に処理を進める。

（Ｓ５２）行列解析部１６３は、次のスレッドを選択する。具体的には、行列解析部１６３は、スレッド番号ｔをインクリメント（現在のスレッド番号ｔに１を加算）する。そして、ステップＳ４２に処理を進める。

図１８は、並列行列演算の手順例を示すフローチャートである。この並列行列演算は、図１５のフローチャートのステップＳ７で実行される。
（Ｓ７０）ベクトル入出力部１６４は、作業領域１２７に含まれる作業ベクトルの要素と結果ベクトル１２９の要素を全て零に初期化する。そして、ベクトル入出力部１６４は、スレッド毎に行列演算部１６７を呼び出す。以下のステップＳ７１〜Ｓ７９の処理が、複数のスレッドで並列に実行される。以下では、スレッド番号＝ｔのスレッド（スレッド＃ｔ）がステップＳ７１〜Ｓ７９の処理を行う場合を説明する。

（Ｓ７１）行列演算部１６７は、スレッドポインタ配列１２４を参照して、スレッド＃ｔに割り当てられた列集合の中の先頭の列を選択する。具体的には、行列演算部１６７は、列番号ｃ＝Ｂｐ（ｔ）を特定する。

（Ｓ７２）行列演算部１６７は、列ポインタ配列１２３を参照して、ステップＳ７１またはステップＳ７９で選択した列に含まれる先頭の非零要素を選択する。具体的には、行列演算部１６７は、要素番号ｅ＝Ｃｐ（ｃ）を特定する。

（Ｓ７３）行列演算部１６７は、下三角の領域に関して、ステップＳ７２またはステップＳ７７で選択した非零要素の値と入力ベクトル１２８のｃ行目の値との積を計算し、スレッド＃ｔに割り当てられた作業ベクトルに保存する。具体的には、行列演算部１６７は、Ｖａｌ（ｅ）Ｘ（ｃ）をＷｏｒｋ（Ｒｏｗ（ｅ），Ｕｐ（ｔ））に加算する。

（Ｓ７４）行列演算部１６７は、ステップＳ７２またはステップＳ７７で選択した非零要素が、対称スパース行列の対角要素であるか判断する。具体的には、行列演算部１６７は、Ｒｏｗ（ｅ）と列番号ｃが一致するか判断する。対角要素である場合はステップＳ７６に処理を進め、対角要素でない場合はステップＳ７５に処理を進める。

（Ｓ７５）行列演算部１６７は、上三角の領域（対角要素を除く）に関して、ステップＳ７２またはステップＳ７７で選択した非零要素の値と入力ベクトル１２８の当該非零要素に対応する行の値との積を計算し、結果ベクトル１２９に保存する。具体的には、行列演算部１６７は、Ｖａｌ（ｅ）Ｘ（Ｒｏｗ（ｅ））をＹ（ｃ）に加算する。

（Ｓ７６）行列演算部１６７は、列ポインタ配列１２３を参照して、ステップＳ７１またはステップＳ７９で選択した列に含まれる全ての非零要素を選択したか判断する。具体的には、行列演算部１６７は、ステップＳ７２またはステップＳ７７で特定した要素番号ｅがＣｐ（ｃ＋１）−１に一致するか判断する。全て選択した場合はステップＳ７８に処理を進め、未選択のものがある場合はステップＳ７７に処理を進める。

（Ｓ７７）行列演算部１６７は、次の非零要素を選択する。具体的には、行列演算部１６７は、要素番号ｅをインクリメントする。そして、ステップＳ７３に処理を進める。
（Ｓ７８）行列演算部１６７は、スレッドポインタ配列１２４を参照して、スレッド＃ｔに割り当てられた列を全て選択したか判断する。具体的には、行列演算部１６７は、列番号ｃがＢｐ（ｔ＋１）−１に一致するか判断する。全て選択した場合、行列演算部１６７は、部分行列と入力ベクトル１２８との積の計算を終了してベクトル入出力部１６４に完了を通知する。そして、ステップＳ８０に処理を進める。未選択の列がある場合、ステップＳ７９に処理を進める。

（Ｓ７９）行列演算部１６７は、次の列を選択する。具体的には、行列演算部１６７は、列番号ｃをインクリメントする。そして、ステップＳ７２に処理を進める。
（Ｓ８０）ベクトル入出力部１６４は、作業領域１２７に含まれる全ての作業ベクトルを結果ベクトル１２９に足し合わせる。具体的には、ベクトル入出力部１６４は、作業領域１２７のｉ行ｊ列の値をＹ（ｉ）に加算する。そして、ベクトル入出力部１６４は、行列演算要求部１６１に、行列演算の完了を通知する。

なお、第２の実施の形態では、対称スパース行列と列ベクトルとの積を計算した。しかし、上記の作業ベクトルの割り当て方法は、行ベクトルと対称スパース行列との積を計算する場合にも適用できる。その場合、対称スパース行列の上三角の領域の部分行列と入力ベクトル１２８との積を保存するために作業ベクトルが使用される。

また、第２の実施の形態では、対称スパース行列を圧縮列格納法で表現した。しかし、対称スパース行列を圧縮行格納法で表現してもよいし、圧縮せずに表現してもよい。対称の位置にある部分行列を同じスレッドに割り当てる限り、上三角の領域と下三角の領域の何れか一方についてスレッド間のアクセス競合の問題が生じる。

また、第２の実施の形態では、ベクトルと掛け算されるスパース行列が対称行列であるとした。しかし、スパース行列が対称行列でなくてもよい。圧縮列格納法で表現されたスパース行列の列を分割し、列ベクトルとの積を求める場合、スレッド間のアクセス競合の問題が生じる。また、圧縮行格納法で表現されたスパース行列の行を分割し、行ベクトルとの積を求める場合にも、スレッド間のアクセス競合の問題が生じる。

第２の実施の形態の情報処理装置１００によれば、複数のスレッドに割り当てられた部分行列の間で非零要素の分布が比較され、非零要素の分布が重複しないようなスレッドの組み合わせに対して共通の作業ベクトルが割り当てられる。よって、作業ベクトルの同じ要素に対するアクセス競合を防ぐことができ、複数のスレッド間で作業ベクトルへのアクセスの排他制御を行わなくてよい。また、アクセス競合が発生しない範囲で、作業領域１２７に確保する作業ベクトルの数を少なくする（例えば、最小化する）ことができる。よって、スレッド数が増加しても作業領域１２７のサイズを抑制することができる。このように、情報処理装置１００によれば、ＲＡＭ１２０の記憶領域を効率的に利用することができ、大規模なスパース行列とベクトルとの積を効率的に計算できる。

なお、前述のように、第１の実施の形態の情報処理は、情報処理装置１０にプログラムを実行させることで実現することができる。また、第２の実施の形態の情報処理は、情報処理装置１００にプログラムを実行させることで実現することができる。

プログラムは、コンピュータ読み取り可能な記録媒体（例えば、記録媒体１４３）に記録しておくことができる。記録媒体としては、例えば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリなどを使用できる。磁気ディスクには、ＦＤおよびＨＤＤが含まれる。光ディスクには、ＣＤ、ＣＤ−Ｒ（Recordable）／ＲＷ（Rewritable）、ＤＶＤおよびＤＶＤ−Ｒ／ＲＷが含まれる。プログラムは、可搬型の記録媒体に記録されて配布されることがある。その場合、可搬型の記録媒体からＨＤＤなどの他の記録媒体（例えば、ＨＤＤ１３１）にプログラムを複製して（インストールして）実行してもよい。

１０ 情報処理装置
１１，１２，１３ プロセッサ
１４ メモリ
２１，２２ スレッド
２３ 行列
２４，２５ 部分行列
２６，２７ 記憶領域

Claims (4)

1. 複数のスレッドを並列に実行可能なコンピュータに、
   零要素および非零要素を含む行列の中の２以上の行および１以上の列によって特定される第１の部分行列を、第１のスレッドに割り当て、前記第１の部分行列と少なくとも一部が重複する２以上の行および前記第１の部分行列と重複しない１以上の列によって特定される前記行列の中の第２の部分行列を、第２のスレッドに割り当て、
   前記第１の部分行列および前記第２の部分行列それぞれについて非零要素が存在する行を検出し、前記第１の部分行列と前記第２の部分行列との間で非零要素が同じ行に存在しないことを示す所定条件を満たすか判定し、
   前記行列の行数に応じた大きさの第１の記憶領域を前記第１のスレッドに割り当てることで、前記第１の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第１の記憶領域内の位置に書き込ませ、
   前記所定条件を満たす場合には、前記第１の記憶領域を更に前記第２のスレッドに割り当てることで、前記第２の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第１の記憶領域内の位置に書き込ませ、前記所定条件を満たさない場合には、前記行列の行数に応じた大きさの第２の記憶領域を前記第２のスレッドに割り当てることで、前記第２の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第２の記憶領域内の位置に書き込ませる、
   処理を実行させるプログラム。
2. 前記所定条件の判定では、前記行列の行を複数の区間に分割し、前記第１の部分行列の各区間内に非零要素が存在するか否かを示す第１のビットマップを生成し、前記第２の部分行列の各区間内に非零要素が存在するか否かを示す第２のビットマップを生成し、前記第１のビットマップと前記第２のビットマップとを比較する、
   請求項１記載のプログラム。
3. 複数のスレッドを並列に実行可能なコンピュータが行う並列演算方法であって、
   零要素および非零要素を含む行列の中の２以上の行および１以上の列によって特定される第１の部分行列を、第１のスレッドに割り当て、前記第１の部分行列と少なくとも一部が重複する２以上の行および前記第１の部分行列と重複しない１以上の列によって特定される前記行列の中の第２の部分行列を、第２のスレッドに割り当て、
   前記第１の部分行列および前記第２の部分行列それぞれについて非零要素が存在する行を検出し、前記第１の部分行列と前記第２の部分行列との間で非零要素が同じ行に存在しないことを示す所定条件を満たすか判定し、
   前記行列の行数に応じた大きさの第１の記憶領域を前記第１のスレッドに割り当てることで、前記第１の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第１の記憶領域内の位置に書き込ませ、
   前記所定条件を満たす場合には、前記第１の記憶領域を更に前記第２のスレッドに割り当てることで、前記第２の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第１の記憶領域内の位置に書き込ませ、前記所定条件を満たさない場合には、前記行列の行数に応じた大きさの第２の記憶領域を前記第２のスレッドに割り当てることで、前記第２の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第２の記憶領域内の位置に書き込ませる、
   並列演算方法。
4. 互いに並列にスレッドを実行可能な複数のプロセッサと、
   メモリと、
   を有し、前記複数のプロセッサの１つは、
   零要素および非零要素を含む行列の中の２以上の行および１以上の列によって特定される第１の部分行列を、第１のスレッドに割り当て、前記第１の部分行列と少なくとも一部が重複する２以上の行および前記第１の部分行列と重複しない１以上の列によって特定される前記行列の中の第２の部分行列を、第２のスレッドに割り当て、
   前記第１の部分行列および前記第２の部分行列それぞれについて非零要素が存在する行を検出し、前記第１の部分行列と前記第２の部分行列との間で非零要素が同じ行に存在しないことを示す所定条件を満たすか判定し、
   前記行列の行数に応じた大きさの前記メモリ内の第１の記憶領域を前記第１のスレッドに割り当てることで、前記第１の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第１の記憶領域内の位置に書き込ませ、
   前記所定条件を満たす場合には、前記第１の記憶領域を更に前記第２のスレッドに割り当てることで、前記第２の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第１の記憶領域内の位置に書き込ませ、前記所定条件を満たさない場合には、前記行列の行数に応じた大きさの前記メモリ内の第２の記憶領域を前記第２のスレッドに割り当てることで、前記第２の部分行列の非零要素に基づいて算出される値を当該非零要素が存在する行に対応する前記第２の記憶領域内の位置に書き込ませる、
   情報処理装置。
   

Images