JP2019148969A - 行列演算装置、行列演算方法および行列演算プログラム - Google Patents

Abstract

【課題】行列積演算の並列処理化を効率的に行えるようにする。【解決手段】行列１５に含まれる複数の第１の行それぞれについて値がゼロでない非ゼロ要素の数をカウントし、非ゼロ要素の数の最大値を判定する。各第１の行から非ゼロ要素の値と列識別子とのペアを抽出し、非ゼロ要素の数が最大値より少ない第１の行に対しては値がゼロであるダミーのペアを追加することで、各第１の行に対して共通する個数のペアを含む圧縮格納データ１８を生成する。圧縮格納データ１８に含まれる各ペアに対して、行列１６から列識別子に対応する行識別子をもつ第２の行を抽出し、抽出した第２の行に対して当該ペアの値を乗算することで行ベクトルを生成する。各第１の行に対して共通する個数のスレッドを割り当て、各第１の行についてスレッドを用いて行ベクトルを合算することで、行列１５と行列１６との行列積を示す行列１７を生成する。【選択図】図１

Description

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

科学技術計算などのある種の計算分野では、ゼロ要素（値がゼロである要素）が多く非ゼロ要素（値がゼロでない要素）が少ない疎行列を扱うことがある。疎行列の内部表現形式として、全ての要素の値を列挙する通常の格納方法を用いると非効率であることから、圧縮行格納（ＣＳＲ：Compressed Sparse Row）法や圧縮列格納（ＣＳＣ：Compressed Sparse Column）法などの圧縮格納法を用いることがある。

ＣＳＲ法では、非ゼロ要素の値と当該非ゼロ要素の列番号との組を列挙した非ゼロ要素リストが生成される。非ゼロ要素リストは、疎行列の非ゼロ要素を行番号の小さい順に列挙したものであって、同じ行の中では列番号の小さい順に列挙したものである。すなわち、非ゼロ要素リストは、疎行列からゼロ要素を除去して非ゼロ要素を左詰めし、１番目の行の非ゼロ要素、２番目の行の非ゼロ要素、…と並べて一次元化したリストである。また、ＣＳＲ法では、非ゼロ要素リストのみでは行の区切りが不明であるため、各行の最初の非ゼロ要素が非ゼロ要素リストの何番目に出現するかを表す行リストが生成される。この非ゼロ要素リストと行リストによって疎行列が表現される。ＣＳＣ法はＣＳＲ法の行と列を入れ替えたものであり、非ゼロ要素リストと列リストが生成される。

疎行列を扱う計算分野では、疎行列と他の行列（密行列やベクトルなどであってもよい）との行列積をコンピュータに計算させることがある。このとき、コンピュータが保持する疎行列のデータは圧縮格納法によって表現されていることがある。

例えば、疎行列とベクトルの積の演算を高速化する行列ベクトル積演算システムが提案されている。提案の行列ベクトル積演算システムは、ＣＳＲ形式の疎行列が入力されると、疎行列のデータ構造をＣＳＲ形式からＪＡＤ（Jagged Diagonal）形式に変換し、ＪＡＤ形式の疎行列とベクトルとの積を複数のプロセッサを用いて並列に計算する。ＪＡＤ形式は、疎行列に含まれる複数の行を非ゼロ要素の多い順に並べ替え、行毎に非ゼロ要素を抽出して左詰めし、それら非ゼロ要素を列方向に辿って一次元化したデータ形式である。

また、例えば、ＣＳＲ形式の疎行列とベクトルとの積を、ＣＳＲ形式のままで複数のプロセッサを用いて並列に計算する行列演算方法が提案されている。また、例えば、疎行列のうち非ゼロ要素の数が閾値以上である列のデータをＪＡＤ形式で保持し、非ゼロ要素の数が閾値未満である列のデータをＣＳＲ形式で保持し、当該疎行列とベクトルとの積を計算する情報処理装置が提案されている。

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

コンピュータが、従来の圧縮格納法によって表された行列Ｓと他の行列Ｄとの行列積を計算して行列Ｏを生成することを考える。例えば、コンピュータは、従来のＣＳＲ法によって表された行列Ｓを用いてＳ×Ｄ＝Ｏを以下のように計算することが考えられる。

行列Ｓのｉ行目には２つの非ゼロ要素があり、ｉ行ｊ_１列の要素Ｓ［ｉ，ｊ_１］とｉ行ｊ_２列の要素Ｓ［ｉ，ｊ_２］が非ゼロ要素であるとする。コンピュータは、行列Ｄからｊ_１行目の行ベクトルを抽出し、抽出した行ベクトルの各要素に対して要素Ｓ［ｉ，ｊ_１］の値を乗算する。また、コンピュータは、行列Ｄからｊ_２行目の行ベクトルを抽出し、抽出した行ベクトルの各要素に対して要素Ｓ［ｉ，ｊ_２］の値を乗算する。この２つの行ベクトルを合算したものが行列Ｏのｉ行目に相当する。よって、コンピュータは、ＣＳＲ形式の行列Ｓに含まれる非ゼロ要素リストのレコード毎に、行列Ｄから行ベクトルを抽出して当該抽出した行ベクトルの各要素に対して乗算を行い、行列Ｓの行毎にそれら行ベクトルを合算することで行列Ｏを生成することができる。

また、例えば、コンピュータは、従来のＣＳＣ法によって表された行列Ｓを用いてＤ×Ｓ＝Ｏを以下のように計算することも考えられる。
行列Ｓのｊ列目には２つの非ゼロ要素があり、ｉ_１行ｊ列の要素Ｓ［ｉ_１，ｊ］とｉ_２行ｊ列の要素Ｓ［ｉ_２，ｊ］が非ゼロ要素であるとする。コンピュータは、行列Ｄからｉ_１列目の列ベクトルを抽出し、抽出した列ベクトルの各要素に対して要素Ｓ［ｉ_１，ｊ］の値を乗算する。また、コンピュータは、行列Ｄからｉ_２列目の列ベクトルを抽出し、抽出した列ベクトルの各要素に対して要素Ｓ［ｉ_２，ｊ］の値を乗算する。この２つの列ベクトルを合算したものが行列Ｏのｊ列目に相当する。よって、コンピュータは、ＣＳＣ形式の行列Ｓに含まれる非ゼロ要素リストのレコード毎に、行列Ｄから列ベクトルを抽出して当該抽出した列ベクトルの各要素に対して乗算を行い、行列Ｓの列毎にそれら列ベクトルを合算することで行列Ｏを生成することができる。

コンピュータは、行列Ｓと行列Ｄの行列積を、複数のスレッドを用いて並列処理化することも考えられる。しかし、従来の圧縮格納法によって表された行列Ｓをそのまま使用した場合、複数のベクトルの合算を並列処理化することが非効率になるという問題がある。

例えば、行列Ｓの行毎の非ゼロ要素数が可変であるため、上記のＳ×Ｄ＝Ｏの計算では合算すべき行ベクトルの数が行列Ｓの行によって異なり、合算する行ベクトルの範囲とスレッドとを対応付けるスレッド割り当ての制御が複雑になってしまう。また、例えば、行列Ｓの列毎の非ゼロ要素数が可変であるため、上記のＤ×Ｓ＝Ｏの計算では合算すべき列ベクトルの数が行列Ｓの列によって異なり、合算する列ベクトルの範囲とスレッドとを対応付けるスレッド割り当ての制御が複雑になってしまう。

１つの側面では、本発明は、行列積演算の並列処理化を効率的に行えるようにする行列演算装置、行列演算方法および行列演算プログラムを提供することを目的とする。

１つの態様では、記憶部と処理部とを有する行列演算装置が提供される。記憶部は、行列演算プログラムを記憶する。処理部は、行列演算プログラムに基づいて複数のスレッドを並列に実行可能である。行列演算プログラムを実行する処理部は、第１の行列に含まれる複数の第１の行それぞれについて値がゼロでない非ゼロ要素の数をカウントし、複数の第１の行の間で非ゼロ要素の数の最大値を判定する。処理部は、複数の第１の行それぞれから非ゼロ要素の値と当該非ゼロ要素が位置する列を示す列識別子とのペアを抽出し、非ゼロ要素の数が最大値より少ない第１の行に対しては値がゼロであるダミーのペアを追加することで、複数の第１の行それぞれに対して共通する個数のペアを含む圧縮格納データを生成する。処理部は、圧縮格納データに含まれるペアそれぞれに対して、第２の行列から当該ペアの列識別子に対応する行識別子をもつ第２の行を抽出し、当該抽出した第２の行に対して当該ペアの値を乗算することで、当該ペアに対応する行ベクトルを生成する。処理部は、複数の第１の行それぞれに対して共通する個数のスレッドを割り当て、複数の第１の行それぞれについて共通する個数のスレッドを用いて行ベクトルを合算することで、第１の行列と第２の行列との行列積を示す第３の行列を生成する。

また、１つの態様では、情報処理装置が実行する行列演算方法が提供される。また、１つの態様では、コンピュータに実行させる行列演算プログラムが提供される。

１つの側面では、行列積演算の並列処理化を効率的に行える。

行列演算装置の例を説明する図である。 情報処理装置のハードウェア例を示すブロック図である。 情報処理装置のソフトウェア構成例を示すブロック図である。 行列積演算の例を示す図である。 第１のＣＳＲデータの例を示す図である。 第１のＣＳＲデータを用いた行列積演算の例を示す図である。 第１の行列積演算における集計処理の例を示す図である。 第１の集計処理におけるスレッド割り当て例を示す図である。 第１のＣＳＲデータ生成の手順例を示すフローチャートである。 第１の行列積演算の手順例を示すフローチャートである。 第１の行列積演算の手順例を示すフローチャート（続き）である。 第２のＣＳＲデータの例を示す図である。 第２のＣＳＲデータを用いた行列積演算の例を示す図である。 第２の行列積演算における集計処理の例を示す図である。 第２の集計処理におけるスレッド割り当て例を示す図である。 第２のＣＳＲデータ生成の手順例を示すフローチャートである。 第２のＣＳＲデータ生成の手順例を示すフローチャート（続き）である。 第２の行列積演算の手順例を示すフローチャートである。 第２の行列積演算の手順例を示すフローチャート（続き）である。

以下、本実施の形態を図面を参照して説明する。
［第１の実施の形態］
第１の実施の形態を説明する。

図１は、行列演算装置の例を説明する図である。
第１の実施の形態の行列演算装置１０は、２つの行列の行列積を計算するコンピュータである。行列演算装置１０は、大規模な疎行列を扱う科学技術計算に用いられることがある。行列演算装置１０は、クライアントコンピュータでもよいしサーバコンピュータでもよい。行列演算装置１０は、記憶部１１および処理部１２を有する。

記憶部１１は、行列演算プログラム１３を記憶する。行列演算プログラム１３は、後述する行列積演算を処理部１２に実行させるプログラムである。行列演算プログラム１３は、ユーザが作成したユーザプログラムでもよいし、コンパイラやリンカなどの変換ソフトウェアを用いてユーザプログラムから変換されたプログラムでもよいし、ユーザプログラムから呼び出されるライブラリプログラムでもよい。記憶部１１は、ＲＡＭ（Random Access Memory）などの揮発性半導体メモリでもよいし、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性ストレージでもよい。

処理部１２は、記憶部１１に記憶された行列演算プログラム１３を実行する。処理部１２は、複数のスレッドを並列に実行することができる複数の演算部を有する。複数の演算部は、プロセッサコアでもよいしＡＬＵ（Arithmetic Logic Unit）などの比較的小さな演算回路でもよい。処理部１２は、ＣＰＵ（Central Processing Unit）でもよいし、ＧＰＵ（Graphics Processing Unit）、ＧＰＧＰＵ（General Purpose Computing on GPU）やＤＳＰ（Digital Signal Processor）でもよい。処理部１２は、数千個から数万個の多数の演算部を有してもよく、それら多数の演算部を用いて数千個から数万個の多数のスレッドを並列に実行可能であってもよい。処理部１２は、行列演算プログラム１３に基づいて、後述する行列積演算のためにスレッド１４ａ，１４ｂなどの複数のスレッドを起動して並列に実行する。

行列演算プログラム１３を実行する処理部１２は、行列１５（第１の行列、行列Ｓ）と行列１６（第２の行列、行列Ｄ）から、行列１５と行列１６の行列積を示す行列１７（第３の行列、行列Ｏ）を生成する。行列１５は、行数や列数の大きい大規模行列であって、値がゼロでない要素である非ゼロ要素が比較的少なく値がゼロの要素であるゼロ要素が比較的多い疎行列である。行列１５は正方行列であってもよい。後述するように行列演算装置１０は、行列１５のデータを、圧縮格納法（圧縮行格納法または圧縮列格納法）を拡張した方法を用いて保持する。行列１６は、行列１５との行列積を計算可能な行列、すなわち、行数と列数の少なくとも一方が行列１５と同じ行列である。行列１６は、非ゼロ要素が比較的多くゼロ要素が比較的少ない密行列であってもよい。行列演算装置１０は、行列１６のデータを圧縮格納法を用いずに保持してもよい。

以下では、行列１５が圧縮行格納法を拡張した方法で表現され、Ｓ×Ｄ＝Ｏという行列積演算を行う場合を説明する。よって以下の例では、行列１５の列数と行列１６の行数は同じである。後述するように、行列１５が圧縮列格納法を拡張した方法で表現され、Ｄ×Ｓ＝Ｏという行列積演算を行う場合には、「行」と「列」を入れ替えて解釈すればよい。また、Ｄ×Ｓ＝（Ｓ^Ｔ×Ｄ^Ｔ）^Ｔであるため、行列１５と行列１６をそれぞれ転置することで、Ｓ×Ｄ＝Ｏと同様の方法でＤ×Ｓ＝Ｏを計算することもできる。その場合、行列１６の列数と行列１５の行数は同じである。

処理部１２は、行列１５に含まれる複数の行それぞれについて非ゼロ要素の数をカウントし、それら複数の行の間で非ゼロ要素数の最大値を判定する。図１には説明を簡単にするため、４行４列の行列１５の例が記載されている。この行列１５に含まれる要素のうち、Ｓ［０，２］，Ｓ［１，０］，Ｓ［１，２］，Ｓ［２，１］，Ｓ［２，３］，Ｓ［３，０］が非ゼロ要素である。よって、行＃０の非ゼロ要素数は１、行＃１の非ゼロ要素数は２、行＃２の非ゼロ要素数は２、行＃３の非ゼロ要素数は１であり、行＃０〜＃３の非ゼロ要素数の最大値は２となる。なお、行番号および列番号はゼロから始まるものとする。

次に、処理部１２は、行列１５を表す圧縮格納データ１８を生成する。処理部１２は、圧縮格納法を用いずに表現された二次元構造データから圧縮格納データ１８を生成してもよいし、圧縮行格納法など通常の圧縮格納法を用いて表現されたデータを圧縮格納データ１８に変換してもよい。前者の場合は処理部１２は二次元構造データを受け取り、後者の場合は処理部１２は通常の圧縮格納法のデータを受け取る。

このとき、処理部１２は、行列１５の行それぞれから、非ゼロ要素の値と当該非ゼロ要素が位置する列を示す列識別子（例えば、列番号）とのペアを抽出する。処理部１２は、抽出したペアを圧縮格納データ１８に登録する。ただし、処理部１２は、非ゼロ要素数が最大値より少ない行については値がゼロであるダミーのペアを追加することで、圧縮格納データ１８において、行列１５の行毎のペア数が共通する個数になるようにする。すなわち、圧縮格納データ１８に登録された行毎のペア数を統一する。共通のペア数は、例えば、上記で判定した非ゼロ要素数の最大値とする。ダミーのペアに含まれる列識別子は任意の列識別子でよく、例えば、列＃０を示す列識別子とする。なお、非ゼロ要素の値と列識別子とは単一のテーブルで管理されなくても、両者が対応付けられていればよく、非ゼロ要素の値と列識別子とが異なるテーブルまたは異なる配列によって管理されてもよい。

上記で述べた図１の例の場合、（値，列識別子）のペアとして、行＃０からは（１，列＃２）が抽出される。また、行＃０の非ゼロ要素数は最大値未満であるため、ここでは（０，列＃０）というダミーのペアが追加されている。行＃１からは（２，列＃０）と（３，列＃２）が抽出される。行＃０の非ゼロ要素数は最大値であるため、このではダミーのペアは追加されていない。行＃２からは（１，列＃１）と（２，列＃３）が抽出される。行＃２の非ゼロ要素数は最大値であるため、ここではダミーのペアは追加されていない。行＃３からは（３，列＃０）が抽出される。行＃３の非ゼロ要素数は最大値未満であるため、ここでは（０，列＃０）というダミーのペアが追加されている。これにより、圧縮格納データ１８では行毎のペア数が２個に統一される。

次に、処理部１２は、圧縮格納データ１８に含まれる各ペアについて、行列１６から、当該ペアに含まれる列識別子に対応する行識別子（例えば、列番号と同じ行番号）をもつ行を抽出する。処理部１２は、抽出した行の各要素に対して当該ペアに含まれる値を乗算する。これにより、圧縮格納データ１８に含まれるペア毎に行ベクトルが生成される。それら行ベクトルを列挙したものがベクトルデータ１９である。処理部１２は、複数のスレッドを用いてベクトルデータ１９の生成を並列処理化してもよい。例えば、異なるペアに対応する行ベクトルの生成を異なるスレッドに実行させる。

上記で述べた図１の例の場合、行列１５の行＃０について、行列１６の行＃２を抽出して各要素に１を乗じた行ベクトルと、行列１６の行＃０を抽出して各要素にゼロを乗じた行ベクトルが生成される。行列１５の行＃１について、行列１６の行＃０を抽出して各要素に２を乗じた行ベクトルと、行列１６の行＃２を抽出して各要素に３を乗じた行ベクトルが生成される。行列１５の行＃２について、行列１６の行＃１を抽出して各要素に１を乗じた行ベクトルと、行列１６の行＃３を抽出して各要素に２を乗じた行ベクトルが生成される。行列１５の行＃３について、行列１６の行＃０を抽出して各要素に３を乗じた行ベクトルと、行列１６の行＃０を抽出して各要素に０を乗じた行ベクトルが生成される。

ダミーのペアに含まれる値はゼロであるため、ダミーのペアに対応する行ベクトルは全ての要素がゼロ要素であるゼロベクトルとなる。並列処理の制御を効率化するため、ダミーのペアに対しても他のペアと同様の手順で行ベクトルを生成することが好ましい。

次に、処理部１２は、行列１５の各行に対して共通する個数のスレッドを割り当てる。すなわち、行列１５の各行に割り当てるスレッドの数を統一する。共通のスレッド数は、例えば、圧縮格納データ１８における共通のペア数と行列１６の列数（すなわち、行ベクトルの列数）とから決定される。一例として、共通のペア数を２で割って小数点以下を切り捨てた整数に、行列１６の列数を乗じた数を、共通のスレッド数とする。上記で述べた図１の例の場合、共通のペア数が２であり行列１６の列数が２であるため、共通のスレッド数は２となり、行列１５の各行に対して２個のスレッドが割り当てられる。

次に、処理部１２は、行列１５の各行について、当該行に割り当てられたスレッドを用いて、当該行に対応する行ベクトルを合算することを並列処理化する。行列１５の各行について、共通するペア数に相当する数の行ベクトルが生成されているため、これら行ベクトルの間で同じ列の要素の値同士を足し合わせることになる。４つ以上の行ベクトルの合算は、例えば、二分木のように２つの行ベクトルを合算することを階層的に繰り返すことによって実行し得る。ある２つの行ベクトルの合算と別の２つの行ベクトルの合算とは、異なるスレッドを用いて並列に実行し得る。また、ある２つの行ベクトルの中で、ある列の値の加算と別の列の値の加算も異なるスレッドを用いて並列に実行し得る。このとき、並列処理の制御を効率化するため、ダミーのペアから生成されたゼロベクトルについても他の行ベクトルと同様の手順で合算を行うことが好ましい。

そして、処理部１２は、行列１５の各行に対応する行ベクトルの合算結果を、当該行列１５の行に対応する行列１７の行として使用する。これにより、行列１５と行列１６の行列積を示す行列１７が生成される。行ベクトルの合算は、ベクトルデータ１９の一部の行ベクトルを書き換えていくことで、ベクトルデータ１９を記憶した記憶領域の中で行うことも可能であり、２つの行ベクトルを合算する毎に新たな記憶領域を使用しなくてもよい。ベクトルデータ１９を記憶した記憶領域の中で合算を行った場合、最終的な合算結果に相当する一部の行ベクトルがベクトルデータ１９から抽出されて行列１７が生成される。ただし、ベクトルデータ１９の記憶領域とは別に行列１７の記憶領域を用意する代わりに、ベクトルデータ１９の一部の行ベクトルのみが見えるビューを定義してもよい。これにより、アプリケーションからはベクトルデータ１９のサブセットが行列１７に見える。

上記で述べた図１の例の場合、行列１５の行＃０に対応する２つの行ベクトルが合算されて行列１７の行＃０が生成される。ただし、この２つの行ベクトルの一方はゼロベクトルである。行列１５の行＃１に対応する２つの行ベクトルが合算されて行列１７の行＃１が生成される。行列１５の行＃２に対応する２つの行ベクトルが合算されて行列１７の行＃２が生成される。行列１５の行＃３に対応する２つの行ベクトルが合算されて行列１７の行＃３が生成される。ただし、この２つの行ベクトルの一方はゼロベクトルである。このように、ベクトルデータ１９には２つのゼロベクトルが存在するものの、行ベクトルを合算するスレッド毎の処理は均一にすることが可能である。

上記ではＳ×Ｄ＝Ｏという行列積演算を説明したが、Ｄ×Ｓ＝Ｏという行列積演算も「行」と「列」を入れ替えることで可能となる。
すなわち、処理部１２は、行列１５の列毎に非ゼロ要素の数をカウントし、非ゼロ要素数の最大値を判定する。処理部１２は、行列１５の各列から非ゼロ要素の値と行識別子とのペアを抽出し、非ゼロ要素数が最大値より少ない列に対しては値がゼロのダミーのペアを追加することで、行列１５の列の間で共通する個数のペアを含む圧縮格納データ１８を生成する。処理部１２は、圧縮格納データ１８に含まれるペア毎に、当該ペアの行識別子に対応する列識別子をもつ行列１６の列を抽出し、抽出した列の各要素に当該ペアの値を乗じることで列ベクトルを生成する。これら列ベクトルを列挙したものがベクトルデータ１９となる。処理部１２は、行列１５の各列について、共通する個数のスレッドを割り当て、割り当てたスレッドを用いて当該列に対応する列ベクトルを合算することで、行列１６と行列１５の行列積を示す行列１７を生成する。

第１の実施の形態の行列演算装置１０によれば、圧縮格納データ１８に値がゼロのダミーのペアを追加することで、行列１５の行毎のペア数（または、行列１５の列毎のペア数）が統一される。そのため、行列１５の行毎に生成される行ベクトルの数（または、行列１５の列毎に生成される列ベクトルの数）が統一される。よって、スレッドと合算すべき行ベクトルの範囲（または、合算すべき列ベクトルの範囲）との対応付けが容易となり、行ベクトルを合算する合算処理（または、列ベクトルを合算する合算処理）の並列処理化の制御が簡潔となる。このため、行列積演算の並列処理化が効率的となる。

［第２の実施の形態］
次に、第２の実施の形態を説明する。
図２は、情報処理装置のハードウェア例を示すブロック図である。

第２の実施の形態の情報処理装置１００は、大規模疎行列と密行列の行列積を計算する。情報処理装置１００は、ＣＰＵ１０１、ＧＰＧＰＵ１０２、ＲＡＭ１０３、ＨＤＤ１０４、画像信号処理部１０５、入力信号処理部１０６、媒体リーダ１０７および通信インタフェース１０８を有する。これらのユニットはバスに接続されている。なお、ＲＡＭ１０３またはＨＤＤ１０４は、第１の実施の形態の記憶部１１に対応する。ＧＰＧＰＵ１０２は、第１の実施の形態の処理部１２に対応する。

ＣＰＵ１０１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ１０１は、ＨＤＤ１０４に記憶されたプログラムやデータの少なくとも一部をＲＡＭ１０３にロードし、プログラムを実行する。ＣＰＵ１０１は、ＧＰＧＰＵ１０２を制御し、行列積演算をＧＰＧＰＵ１０２に実行させることがある。

ＧＰＧＰＵ１０２は、画像処理に適した演算器をもつＧＰＵを他の用途に転用したプロセッサである。ＧＰＧＰＵ１０２は、並列にスレッドを実行可能な演算器１０２ａ，１０２ｂ，１０２ｃを含む多数の演算器を有する。これらの演算器は、プロセッサコアでもよいしＡＬＵなどの比較的小さな単位回路でもよい。例えば、ＧＰＧＰＵ１０２は、数千個から数万個の多数の演算器を有し、多数のスレッドを並列に実行することができる。

ＲＡＭ１０３は、ＣＰＵ１０１やＧＰＧＰＵ１０２が実行するプログラムや演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、情報処理装置１００は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。

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

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

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

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

媒体リーダ１０７は、記録媒体１１３から読み取ったプログラムやデータを、ＲＡＭ１０３やＨＤＤ１０４などの他の記録媒体にコピーする。読み取られたプログラムは、ＣＰＵ１０１やＧＰＧＰＵ１０２によって実行され得る。なお、記録媒体１１３は可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体１１３やＨＤＤ１０４を、コンピュータ読み取り可能な記録媒体と言うことがある。

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

図３は、情報処理装置のソフトウェア構成例を示すブロック図である。
情報処理装置１００は、ユーザアプリケーション１２１、行列演算ライブラリ１２２、オペレーティングシステム１２７、スレッドプール１２８およびデータ記憶部１２９を有する。ユーザアプリケーション１２１、行列演算ライブラリ１２２およびオペレーティングシステム１２７は、プログラムを用いて実装される。スレッドプール１２８およびデータ記憶部１２９は、例えば、ＲＡＭ１０３の記憶領域を用いて実装される。

ユーザアプリケーション１２１は、ユーザが作成したユーザプログラムを用いて実装されるアプリケーションソフトウェアである。ユーザアプリケーション１２１は、疎行列と密行列を指定して行列演算ライブラリ１２２を呼び出すことで、疎行列と密行列の行列積を行列演算ライブラリ１２２に計算させて演算結果を取得する。

行列演算ライブラリ１２２は、行列積を計算するプログラムを用いて実装されるライブラリソフトウェアである。行列演算ライブラリ１２２は、行列演算ライブラリ１２２を呼び出すユーザアプリケーション１２１とリンクされる。行列演算ライブラリ１２２は、ユーザアプリケーション１２１の生成時に静的にリンクされる静的リンクライブラリでもよいし、その実行時に動的にリンクされる動的リンクライブラリでもよい。行列演算ライブラリ１２２は、ユーザアプリケーション１２１に取り込まれて配布されてもよい。

行列演算ライブラリ１２２は、ＧＰＧＰＵ１０２に複数のスレッドを並列に実行させて行列積演算を高速化する。行列演算ライブラリ１２２は、疎行列生成部１２３、データ構造変換部１２４、密行列生成部１２５および行列積演算部１２６を有する。

疎行列生成部１２３は、ユーザアプリケーション１２１からの入力に応じて、ゼロ要素も含めて疎行列の全ての要素を二次元配列として並べた二次元構造データの疎行列を生成する。データ構造変換部１２４は、疎行列生成部１２３が生成した二次元構造データを、圧縮行格納法の考え方に基づくＣＳＲデータに変換する。後述するようにＣＳＲデータとしては、通常の圧縮行格納法に従った第１のＣＳＲデータと、圧縮行格納法を拡張した方法に従った第２のＣＳＲデータを生成することが可能である。

密行列生成部１２５は、ユーザアプリケーション１２１からの入力に応じて、密行列の全ての要素を二次元配列として並べた二次元構造データの密行列を生成する。行列積演算部１２６は、データ構造変換部１２４が生成したＣＳＲデータと密行列生成部１２５が生成した密行列の二次元構造データを用いて、疎行列と密行列の行列積である出力行列を生成する。この出力行列は、出力行列の全ての要素を二次元配列として並べた二次元構造データである。後述するように、行列積演算の方法はＣＳＲデータの形式によって異なる。行列積演算部１２６は、出力行列をユーザアプリケーション１２１に出力する。

オペレーティングシステム１２７は、入力デバイス１１２を用いたユーザ入力を受け付け、ユーザ入力に応じた制御を行う。情報処理装置１００は、ユーザアプリケーション１２１を生成するために、コンパイラやリンカなどの開発用ソフトウェアを有してもよい。例えば、オペレーティングシステム１２７が受け付けたユーザ入力に応じて、コンパイラがソースコードをコンパイルしてオブジェクトコードを生成し、リンカが当該オブジェクトコードと行列演算ライブラリ１２２とをリンクする。これにより、行列演算ライブラリ１２２とリンクされたユーザアプリケーション１２１が生成される。

ただし、ユーザアプリケーション１２１に相当する実行可能プログラムは、他の情報処理装置によって生成されて情報処理装置１００に入力されてもよい。また、第２の実施の形態では行列積演算のアルゴリズムが行列演算ライブラリ１２２に記述されているが、ユーザがソースコードの中に当該アルゴリズムを直接記述することも可能である。また、ソースコードに当該アルゴリズムが記述されていない場合であっても、コンパイラが自動的にオブジェクトコードの中に当該アルゴリズムを挿入することも可能である。

スレッドプール１２８は、ＧＰＧＰＵ１０２が有する複数の演算器を用いて並列に実行可能な複数のスレッドを予め起動しておいたプールである。複数のスレッドを予め起動しておくことで、行列積演算の途中におけるスレッド起動のオーバヘッドを削減できる。これら複数のスレッドは行列演算ライブラリ１２２によって起動される。また、これら複数のスレッドは行列積演算部１２６によって行列積演算のために使用される。

データ記憶部１２９は、疎行列生成部１２３が生成した疎行列の二次元構造データと、データ構造変換部１２４が生成したＣＳＲデータと、密行列生成部１２５が生成した密行列の二次元構造データを記憶する。また、データ記憶部１２９は、行列積演算の途中で行列積演算部１２６が生成した中間データと、行列積演算部１２６が生成した出力行列の二次元構造データを記憶する。また、データ記憶部１２９は、スレッドプール１２８が保持するスレッドの割り当て状況を示す情報など各種の制御情報を記憶する。

行列演算ライブラリ１２２は、疎行列のデータ構造として第１のＣＳＲデータを用いる第１の行列積演算と、疎行列のデータ構造として第２のＣＳＲデータを用いる第２の行列積演算のうち、何れか一方または両方を実行できる。行列演算ライブラリ１２２が両方を実行できる場合、ユーザが行列積演算の方法を選択できるようにしてもよい。以下では、まず第１の行列積演算を説明し、その後に第２の行列積演算を説明する。

図４は、行列積演算の例を示す図である。
第２の実施の形態の行列積演算を説明するにあたり、図４に示す比較的コンパクトな疎行列１３１（疎行列Ｓ_２Ｄ）および密行列１３２（密行列Ｄ）を使用する。

疎行列１３１は、ゼロ要素が多い８行８列の行列である。疎行列１３１の行番号は上側の行から下側の行に向かって＃０，＃１，…，＃７と付与されており、疎行列１３１の列番号は左側の列から右側の列に向かって＃０，＃１，…，＃７と付与されている。

疎行列１３１では、［０，５］，［０，７］，［１，１］，［１，３］，［３，４］，［３，５］，［３，６］，［３，７］，［４，１］，［４，２］，［４，５］，［４，７］，［５，４］，［５，６］，［６，２］，［７，２］，［７，３］が非ゼロ要素である。それ以外の疎行列１３１の要素はゼロ要素である。すなわち、疎行列１３１に含まれる６４個の要素のうち１７個のみが非ゼロ要素である。

密行列１３２は、ゼロ要素が少ない８行３列の行列である。密行列１３２の行番号は上側の行から下側の行に向かって＃０，＃１，…，＃７と付与されており、密行列１３２の列番号は左側の列から右側の列に向かって＃０，＃１，＃２と付与されている。

なお、第２の実施の形態では、疎行列Ｓ_２Ｄを圧縮行格納法またはそれを拡張した方法によって表現してＳ_２Ｄ×Ｄを計算するが、疎行列Ｓ_２Ｄを圧縮列格納法またはそれを拡張した方法によって表現してＤ×Ｓ_２Ｄを計算することも可能である。その場合には、以下で説明するアルゴリズムの「行」と「列」を入れ替えればよい。また、Ｄ×Ｓ_２Ｄ＝（Ｓ_２Ｄ ^Ｔ×Ｄ^Ｔ）^Ｔであるため、疎行列Ｓ_２Ｄと密行列Ｄをそれぞれ転置し、生成された出力行列を転置することによっても、Ｄ×Ｓ_２Ｄを計算することができる。

図５は、第１のＣＳＲデータの例を示す図である。
第１の行列積演算では、疎行列１３１から非ゼロ要素テーブル１４１（非ゼロ要素テーブルＳ_Ｅ）および行管理テーブル１４２（行管理テーブルＳ_Ｒ）が生成される。非ゼロ要素テーブル１４１および行管理テーブル１４２は、データ記憶部１２９に記憶される。

非ゼロ要素テーブル１４１は、組番号、値および列番号の項目を含む。組番号の項目には、ゼロから始まる整数の連番であって、疎行列１３１に含まれる非ゼロ要素を識別する数字が登録される。値の項目には、非ゼロ要素の値が登録される。列番号の項目には、非ゼロ要素が位置する列を示す列番号が登録される。非ゼロ要素テーブル１４１における複数の非ゼロ要素は、行番号が小さい順に並べられており、同じ行の中では列番号が小さい順に並べられている。なお、非ゼロ要素を特定できるように格納されていれば、非ゼロ要素テーブル１４１が明示的に組番号の項目を含んでいなくてもよい。

行管理テーブル１４２は、非ゼロ要素テーブル１４１における行の区切りを示す。これは、非ゼロ要素テーブル１４１には非ゼロ要素が位置する列を特定する情報が含まれているものの行を特定する情報が含まれていないためである。行管理テーブル１４２は、行番号および組番号の項目を含む。行番号の項目には、疎行列１３１の行を示す行番号が登録される。組番号の項目には、非ゼロ要素テーブル１４１において各行の最初の非ゼロ要素を示す組番号が登録される。ただし、非ゼロ要素が１つも存在しない行に対しては、次の非ゼロ要素を示す組番号が対応付けられる。また、行管理テーブル１４２の末尾を明確にするため、行番号を「ｓｅｎｔｉｎｅｌ」とし、組番号を最大値より１だけ大きい数字としたレコードが行管理テーブル１４２に登録される。

例えば、非ゼロ要素テーブル１４１では、疎行列１３１に含まれる１７個の非ゼロ要素に対して組番号０〜１６が付与される。行管理テーブル１４２では、行＃０と当該行の最初の非ゼロ要素を示す組番号０が対応付けられる。以下同様に、行＃１と当該行の最初の非ゼロ要素を示す組番号２が対応付けられる。行＃２には非ゼロ要素が存在しないため、行＃３と当該行の最初の非ゼロ要素を示す組番号４が対応付けられ、行＃２と組番号４が対応付けられる。行＃４と当該行の最初の非ゼロ要素を示す組番号８が対応付けられる。行＃５と当該行の最初の非ゼロ要素を示す組番号１２が対応付けられる。行＃６と当該行の最初の非ゼロ要素を示す組番号１４が対応付けられる。行＃７と当該行の最初の非ゼロ要素を示す組番号１５が対応付けられる。そして、「ｓｅｎｔｉｎｅｌ」と組番号１７を対応付けたレコードが行管理テーブル１４２の末尾に登録される。

図６は、第１のＣＳＲデータを用いた行列積演算の例を示す図である。
疎行列１３１が非ゼロ要素テーブル１４１および行管理テーブル１４２として表現されている場合、行列積演算部１２６は以下のようにして行列積を計算することができる。行列積の計算は以下に説明するように、第１工程としての複製処理と、第２の工程としての乗算処理と、第３工程としての集計処理とを含む。

複製処理では、行列積演算部１２６は、非ゼロ要素テーブル１４１から中間行列１４３を生成する。中間行列１４３は、データ記憶部１２９に記憶される。中間行列１４３は、行数が非ゼロ要素テーブル１４１のレコード数と同じ、すなわち、疎行列１３１の非ゼロ要素数と同じであり、列数が密行列１３２の列数と同じである二次元行列である。

行列積演算部１２６は、非ゼロ要素テーブル１４１から列番号を抽出し、その列番号と同じ行番号をもつ行を密行列１３２から抽出して中間行列１４３に格納する。このとき、非ゼロ要素テーブル１４１の列番号の順序と中間行列１４３の行の順序とが対応しているようにする。例えば、行列積演算部１２６は、組番号０に対応する列＃５について、密行列１３２から行＃５を抽出して中間行列１４３の行＃０に複製する。また、行列積演算部１２６は、組番号１に対応する列＃７について、密行列１３２から行＃７を抽出して中間行列１４３の行＃１に複製する。このようにして、行列積演算部１２６は、組番号０〜１６の列番号に基づいて１７個の行が密行列１３２から抽出される。

乗算処理では、行列積演算部１２６は、中間行列１４３から中間行列１４４を生成する。中間行列１４４は、中間行列１４３を更新した行列であり、中間行列１４３の記憶領域を上書きすれば中間行列１４３と異なる記憶領域を使用しなくてもよい。

行列積演算部１２６は、非ゼロ要素テーブル１４１から値を抽出し、その値を中間行列１４３の対応する行の各要素に乗じる。例えば、行列積演算部１２６は、組番号０に対応する値「１」を、中間行列１４３の行＃０の各要素に乗じる（各要素の値を１倍する）。また、行列積演算部１２６は、組番号１に対応する値「２」を、中間行列１４３の行＃１の各要素に乗じる（各要素の値を２倍する）。このようにして、行列積演算部１２６は、組番号０〜１６の値と中間行列１４３の行＃０〜＃１６の値の間でそれぞれ乗算を行う。

集計処理では、行列積演算部１２６は、疎行列１３１の行毎に中間行列の行をグルーピングし、グルーピングした行を合算することで出力行列１４５を生成する。出力行列１４５は、データ記憶部１２９に記憶される。出力行列１４５は、行数が疎行列１３１の行数と同じであり、列数が密行列１３２の列数と同じ二次元行列である。

行列積演算部１２６は、疎行列１３１の行毎に、その疎行列１３１の行に存在する非ゼロ要素から生成された中間行列１４４の行範囲を特定し、特定した行範囲の値を列毎に合計する。合計値が、その疎行列１３１の行に対応する出力行列の行の値となる。例えば、行列積演算部１２６は、疎行列１３１の行＃０に対応する中間行列１４４の行＃０，＃１を特定し、２つの行を合算して出力行列１４５の行＃０を求める。また、行列積演算部１２６は、疎行列１３１の行＃１に対応する中間行列１４４の行＃２，＃３を特定し、２つの行を合算して出力行列１４５の行＃１を求める。疎行列１３１の行＃２に対応する中間行列１４４の行は存在しないため、出力行列１４５の行＃２の全ての要素はゼロとなる。

また、行列積演算部１２６は、疎行列１３１の行＃３に対応する中間行列１４４の行＃４〜＃７を特定し、４つの行を合算して出力行列１４５の行＃３を求める。また、行列積演算部１２６は、疎行列１３１の行＃４に対応する中間行列１４４の行＃８〜＃１１を特定し、４つの行を合算して出力行列１４５の行＃４を求める。また、行列積演算部１２６は、疎行列１３１の行＃５に対応する中間行列１４４の行＃１２，＃１３を特定し、２つの行を合算して出力行列１４５の行＃５を求める。疎行列１３１の行＃６に対応する中間行列１４４の行は行＃１４のみであるため、これが出力行列１４５の行＃６となる。また、行列積演算部１２６は、疎行列１３１の行＃７に対応する中間行列１４４の行＃１５，＃１６を特定し、２つの行を合算して出力行列１４５の行＃７を求める。

ここで、中間行列１４３の生成は複数のスレッドを用いて容易に並列処理化できる。例えば、非ゼロ要素テーブル１４１の異なる非ゼロ要素に対して異なるスレッドを割り当てることで、中間行列１４３の異なる行の生成を並列に実行できる。非ゼロ要素それぞれに対して密行列１３２の列数に相当する数のスレッドを割り当て、中間行列１４３の１つの要素を１つのスレッドによって生成することも可能である。

また、中間行列１４４の生成も複数のスレッドを用いて容易に並列処理化できる。例えば、非ゼロ要素テーブル１４１の異なる非ゼロ要素に対して異なるスレッドを割り当てることで、中間行列１４４の異なる行の乗算を並列に実行できる。非ゼロ要素それぞれに対して密行列１３２の列数に相当する数のスレッドを割り当て、中間行列１４４の１つの要素を１つのスレッドによって計算することも可能である。

これに対し、中間行列１４４を集計して出力行列１４５を生成することについては幾つかの並列処理化の方法が考えられる。第１の方法は、中間行列１４４の１つの要素に対して１つのスレッドを割り当てる方法である。各スレッドは出力行列１４５の何れか１つの要素に対して加算を実行する。ただし、第１の方法では出力行列１４５の同じ要素に対して複数のスレッドがアクセスする可能性があり、排他制御のオーバヘッドが生じる。

第２の方法は、出力行列１４５の１つの要素に対して１つのスレッドを割り当てる方法である。各スレッドは出力行列１４５の何れか１つの要素を独占的に計算するため、排他制御は不要である。ただし、第２の方法では疎行列１３１の中の非ゼロ要素が多い行については、並列度の不足によって計算効率が低くなることがある。疎行列１３１の特定の行にＡ個（Ａは４以上の整数）の非ゼロ要素があり、中間行列１４４のＡ個の行を合算することを考える。Ａ個の行の合算は理論上、二分木に従って２つの行の合算を階層的に繰り返すことで、Ａ／２個のスレッドを用いてｌｏｇ_２Ａステップで実行することができる。これに対して第２の方法では、スレッド数の不足によりＡ−１ステップを要する。

第３の方法は、出力行列１４５の１つの要素に対して可変個のスレッドを割り当てる方法である。例えば、出力行列１４５の１つの要素に対して、Ａ／２個のスレッドなど、中間行列１４４の中の合算すべき行数に応じた数のスレッドを割り当てる。ただし、第３の方法では、疎行列１３１の行毎の非ゼロ要素数が可変であるため、行管理テーブル１４２の先頭から順に走査して非ゼロ要素数を確認し、中間行列１４４の担当行範囲とスレッドとの対応付けを決定することになる。よって、スレッド割り当て処理が複雑となり、スレッド割り当て自体のオーバヘッドが大きくなる。

以下では、第３の方法を用いて集計処理を行うとする。
図７は、第１の行列積演算における集計処理の例を示す図である。
行列積演算部１２６は、中間行列１４４の一部要素を順次上書きすることで集計処理を進める。行列積演算部１２６は、中間行列１４４を更新して中間行列１４４ａを生成し、中間行列１４４ａを更に更新して中間行列１４４ｂを生成する。中間行列１４４ａ，１４４ｂの記憶領域としては、中間行列１４４の記憶領域をそのまま使用すればよい。

行列積演算部１２６は、行管理テーブル１４２を参照し、疎行列１３１の各行から生成された中間行列１４４の行を特定する。例えば、疎行列１３１の行＃０には中間行列１４４の行＃０，＃１が対応する。疎行列１３１の行＃１には中間行列１４４の行＃２，＃３が対応する。疎行列１３１の行＃２に対応する中間行列１４４の行は存在しない。疎行列１３１の行＃３には中間行列１４４の行＃４〜＃７が対応する。疎行列１３１の行＃４には中間行列１４４の行＃８〜＃１１が対応する。疎行列１３１の行＃５には中間行列１４４の行＃１２，＃１３が対応する。疎行列１３１の行＃６には中間行列１４４の行＃１４が対応する。疎行列１３１の行＃７には中間行列１４４の行＃１５，＃１６が対応する。

行列積演算部１２６は、上記のように区分した行範囲それぞれの中で、２つの行を合算する処理を二分木形式で階層的に繰り返す。２つの行の合算は、列毎に、行番号の大きい方の行の値を行番号の小さい方の行に加算することで行う。最初は隣接する行同士が合算され、集計処理のステップが進むにつれて離れた行同士が合算される。集計処理の第Ｉステップ（Ｉは０以上の整数）では、２^Ｉだけ離れた行同士が合算される。１行当たり非ゼロ要素数の最大値をＭとすると、集計処理の終了までのステップ数はｌｏｇ_２Ｍである。

例えば、中間行列１４４を中間行列１４４ａに更新する第０ステップでは、中間行列１４４の行＃１の値が行＃０に加算される。行＃３の値が行＃２に加算される。行＃５の値が行＃４に加算され、行＃７の値が行＃６に加算される。行＃９の値が行＃８に加算され、行＃１１の値が行＃１０に加算される。行＃１３の値が行＃１２に加算される。行＃１４は加算相手がないためそのまま維持する。行＃１６の値が行＃１５に加算される。

次に、中間行列１４４ａを中間行列１４４ｂに更新する第１ステップでは、中間行列１４４ａの行＃０は加算相手がないためそのまま維持する。行＃２は加算相手がないためそのまま維持する。行＃６の値が行＃４に加算される。行＃１０の値が行＃８に加算される。行＃１２は加算相手がないためそのまま維持する。行＃１４は加算相手がないためそのまま維持する。行＃１５は加算相手がないためそのまま維持する。中間行列１４４ａの上記以外の行は、第０ステップで合算済みであるため無視してよい。

この例では１行当たり非ゼロ要素数の最大値は４であるため、第０ステップと第１ステップで集計処理は終了する。疎行列１３１の各行に対応する中間行列１４４ｂの行は高々１個に集約されている。中間行列１４４ｂにおける集約結果の行は、行管理テーブル１４２の組番号に相当する行である。疎行列１３１の行＃０には中間行列１４４ｂの行＃０が対応する。疎行列１３１の行＃１には中間行列１４４ｂの行＃２が対応する。疎行列１３１の行＃２に対応する中間行列１４４ｂの行は存在しない。これは、疎行列１３１の行＃２には要素が全てゼロである行ベクトルが対応することを意味する。疎行列１３１の行＃３には中間行列１４４ｂの行＃４が対応する。疎行列１３１の行＃４には中間行列１４４ｂの行＃８が対応する。疎行列１３１の行＃５には中間行列１４４ｂの行＃１２が対応する。疎行列１３１の行＃６には中間行列１４４ｂの行＃１４が対応する。疎行列１３１の行＃７には中間行列１４４ｂの行＃１５が対応する。中間行列１４４ｂの上記以外の行は、第０ステップおよび第１ステップで合算済みであるため無視してよい。

行列積演算部１２６は、中間行列１４４ｂの一部の行を抽出して出力行列１４５を生成する。中間行列１４４ｂの一部の行を複製して出力行列１４５を生成してもよいし、行の複製を行わずに中間行列１４４ｂの一部の行のみが見えるビューを生成してもよい。

例えば、出力行列１４５の行＃０は中間行列１４４ｂの行＃０である。出力行列１４５の行＃１は中間行列１４４ｂの行＃２である。出力行列１４５の行＃２は要素が全てゼロの行ベクトルである。出力行列１４５の行＃３は中間行列１４４ｂの行＃４である。出力行列１４５の行＃４は中間行列１４４ｂの行＃８である。出力行列１４５の行＃５は中間行列１４４ｂの行＃１２である。出力行列１４５の行＃６は中間行列１４４ｂの行＃１４である。出力行列１４５の行＃７は中間行列１４４ｂの行＃１５である。

行列積演算部１２６は、２つの値を加算する加算演算毎に１つのスレッドを割り当てる。上記の第０ステップでは中間行列１４４を中間行列１４４ａに更新するにあたり、８行×３列＝２４個の加算演算が行われており、２４個のスレッドが並列に実行される。また、上記の第１ステップでは中間行列１４４ａを中間行列１４４ｂに更新するにあたり、２行×３列＝６個の加算演算が行われており、６個のスレッドが並列に実行される。ただし、集約処理のステップが進むにつれて加算演算の数は減少するため、前のステップで使用されたスレッドのサブセットを次のステップで使用すればよい。上記の第１ステップでは、中間行列１４４ａの行＃４を計算したスレッドが中間行列１４４ｂの行＃４を計算すればよく、中間行列１４４ａの行＃８を計算したスレッドが中間行列１４４ｂの行＃８を計算すればよい。よって、スレッド割り当ては集約処理の開始時に決定される。

図８は、第１の集計処理におけるスレッド割り当て例を示す図である。
スレッドテーブル１４６は、中間行列１４４から決定されるスレッド割り当てを示す。スレッドテーブル１４６は、データ記憶部１２９に記憶される。

スレッドテーブル１４６は、行番号およびスレッド番号の項目を有する。行番号の項目には、非ゼロ要素が存在する疎行列１３１の行の行番号が登録される。スレッド番号の項目には、疎行列１３１の行に対して割り当てられたスレッドを識別するスレッド番号が二次元配列として列挙される。この二次元配列の行数は、非ゼロ要素数をＡとするとＡ／２（小数点以下切り捨て）である。この二次元配列の列数は、密行列１３２の列数である。

上記のように、中間行列１４４に対する集計処理の第０ステップでは２４個の加算演算が行われるため、２４個のスレッドが割り当てられる。具体的には、行＃０に対してスレッド＃０〜＃２、行＃１に対してスレッド＃３〜＃５、行＃３に対してスレッド＃６〜＃１１、行＃４に対してスレッド＃１２〜＃１７、行＃５に対してスレッド＃１８〜＃２０、行＃７に対してスレッド＃２１〜＃２３が割り当てられる。

例えば、スレッド＃０は中間行列１４４の［０，０］に［１，０］の値を加算するものである。スレッド＃１は［０，１］に［１，１］の値を加算するものである。スレッド＃２は［０，２］に［１，２］の値を加算するものである。また、スレッド＃６は［４，０］に［５，０］の値を加算するものである。スレッド＃７は［４，１］に［５，１］の値を加算するものである。スレッド＃８は［４，２］に［５，２］の値を加算するものである。また、スレッド＃９は［６，０］に［７，０］の値を加算するものである。スレッド＃１０は［６，１］に［７，１］の値を加算するものである。スレッド＃１１は［６，２］に［７，２］の値を加算するものである。

集計処理の第０ステップでは、これらのスレッドが全て並列に実行される。それより後のステップでは、これらのスレッドの一部が並列に実行される。上記のように、中間行列１４４に対する集計処理の第１ステップでは６個の加算演算が行われるため、スレッド＃０〜＃２３のうちの６個のスレッドが並列に実行される。具体的には、スレッド＃６〜＃８，＃１２〜＃１４が実行され、それ以外のスレッドは実行されない。

例えば、スレッド＃６は中間行列１４４ａの［４，０］に［６，０］の値を加算する。スレッド＃７は［４，１］に［６，１］の値を加算する。スレッド＃８は［４，２］に［６，２］の値を加算する。また、スレッド＃１２は［８，０］に［１０，０］の値を加算する。スレッド＃１３は［８，１］に［１０，１］の値を加算する。スレッド＃１４は［８，２］に［１０，２］の値を加算する。これにより集計処理が終了する。

次に、第１のＣＳＲデータに関する処理手順を説明する。
図９は、第１のＣＳＲデータ生成の手順例を示すフローチャートである。
第１のＣＳＲデータ生成は、ユーザアプリケーション１２１からの入力に応じて疎行列生成部１２３が二次元構造データの疎行列１３１を生成した後に行われる。

（Ｓ１０）データ構造変換部１２４は、組番号Ｎ＝０に初期化する。
（Ｓ１１）データ構造変換部１２４は、疎行列１３１（疎行列Ｓ_２Ｄ）の行番号を小さい方から１つ選択する（行番号ｒｏｗ）。

（Ｓ１２）データ構造変換部１２４は、行管理テーブル１４２（行管理テーブルＳ_Ｒ）の末尾に、行番号ｒｏｗと組番号Ｎを含むレコード｛ｒｏｗ，Ｎ｝を追加する。
（Ｓ１３）データ構造変換部１２４は、行番号ｒｏｗについて、疎行列１３１の列番号を小さい方から１つ選択する（列番号ｃｏｌ）。

（Ｓ１４）データ構造変換部１２４は、疎行列１３１からｒｏｗ行ｃｏｌ列の値Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］（値ｖａｌ）を抽出し、ｖａｌ＝０であるか判断する。ｖａｌ＝０である場合、すなわち、Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］がゼロ要素である場合、ステップＳ１７に処理が進む。ｖａｌ＝０でない場合、すなわち、Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］が非ゼロ要素である場合、ステップＳ１５に処理が進む。

（Ｓ１５）データ構造変換部１２４は、非ゼロ要素テーブル１４１（非ゼロ要素テーブルＳ_Ｅ）の末尾に、組番号Ｎと値ｖａｌと列番号ｃｏｌを含むレコード｛Ｎ，ｖａｌ，ｃｏｌ｝を追加する。

（Ｓ１６）データ構造変換部１２４は、組番号Ｎを１だけ大きくする。
（Ｓ１７）データ構造変換部１２４は、ステップＳ１３で、行番号ｒｏｗについて全ての列番号ｃｏｌを選択したか判断する。全ての列番号ｃｏｌを選択した場合はステップＳ１８に処理が進み、未選択の列番号ｃｏｌがある場合はステップＳ１３に処理が進む。

（Ｓ１８）データ構造変換部１２４は、ステップＳ１１で、全ての行番号ｒｏｗを選択したか判断する。全ての行番号ｒｏｗを選択した場合はステップＳ１９に処理が進み、未選択の行番号ｒｏｗがある場合はステップＳ１１に処理が進む。

（Ｓ１９）データ構造変換部１２４は、行管理テーブル１４２の末尾に、「ｓｅｎｔｉｎｅｌ」と組番号Ｎを含むレコード｛ｓｅｎｔｉｎｅｌ，Ｎ｝を追加する。
図１０は、第１の行列積演算の手順例を示すフローチャートである。

第１の行列積演算は、データ構造変換部１２４が第１のＣＳＲデータを生成し、密行列生成部１２５が二次元構造データの密行列１３２を生成した後に行われる。
（Ｓ２０）行列積演算部１２６は、非ゼロ要素テーブル１４１から組番号を１つ選択し（組番号Ｎ）、組番号Ｎに対応付けられた列番号Ｓ_Ｅ［Ｎ］．ｃｏｌを選択する。

（Ｓ２１）行列積演算部１２６は、密行列１３２（密行列Ｄ）のｃｏｌ行目をコピーし、中間行列１４３（中間行列Ｔ）のＮ行目に格納する。
（Ｓ２２）行列積演算部１２６は、ステップＳ２０で、非ゼロ要素テーブル１４１の全ての組番号Ｎを選択したか判断する。全ての組番号Ｎを選択した場合はステップＳ２３に処理が進み、未選択の組番号Ｎがある場合はステップＳ２０に処理が進む。なお、ステップＳ２０，Ｓ２１の処理は複数のスレッドを用いて並列処理化することができる。

（Ｓ２３）行列積演算部１２６は、中間行列１４３から行番号と列番号を１つずつ選択し（行番号Ｎと列番号ｃｏｌ）、Ｎ行ｃｏｌ列の要素Ｔ［Ｎ］［ｃｏｌ］を選択する。
（Ｓ２４）行列積演算部１２６は、非ゼロ要素テーブル１４１から組番号Ｎに対応付けられた値Ｓ_Ｅ［Ｎ］．ｖａｌを選択する。行列積演算部１２６は、ステップＳ２３で選択したＴ［Ｎ］［ｃｏｌ］にＳ_Ｅ［Ｎ］．ｖａｌを乗じる。

（Ｓ２５）行列積演算部１２６は、ステップＳ２３で、中間行列１４３の全ての要素を選択したか判断する。全ての要素を選択した場合はステップＳ２６に処理が進み、未選択の要素がある場合はステップＳ２３に処理が進む。なお、ステップＳ２３，Ｓ２４の処理は複数のスレッドを用いて並列処理化することができる。

（Ｓ２６）行列積演算部１２６は、最大値Ｍ＝０に初期化する。
（Ｓ２７）行列積演算部１２６は、行管理テーブル１４２の行番号を小さい方から１つ選択する（行番号ｒｏｗ）。

（Ｓ２８）行列積演算部１２６は、行管理テーブル１４２から行番号ｒｏｗに対応する組番号Ｓ_Ｒ［ｒｏｗ］と行番号ｒｏｗ＋１に対応する組番号Ｓ_Ｒ［ｒｏｗ＋１］を検索する。行列積演算部１２６は、組数ｒａｎｇｅ＝Ｓ_Ｒ［ｒｏｗ＋１］−Ｓ_Ｒ［ｒｏｗ］を算出する。また、行列積演算部１２６は、現在の最大値Ｍと組数ｒａｎｇｅの何れか大きい方を最大値Ｍとする（Ｍ＝ｍａｘ（Ｍ，ｒａｎｇｅ））。

（Ｓ２９）行列積演算部１２６は、スレッドプール１２８からｆｌｏｏｒ（ｒａｎｇｅ／２）×ｃｏｌｓ個のスレッドを取得する。ｆｌｏｏｒは床関数であり、正数に対しては小数点以下切り捨てを表す。ｃｏｌｓは、中間行列１４３を更新することで得られた中間行列１４４（中間行列Ｔ）の列数を表す。行列積演算部１２６は、スレッドテーブル１４６（スレッドテーブルＨ）の行番号ｒｏｗに対応する二次元配列Ｈ［ｒｏｗ］［：］［：］に、取得したスレッドのスレッド番号を埋めていく。二次元配列Ｈ［ｒｏｗ］［：］［：］の行数はｆｌｏｏｒ（ｒａｎｇｅ／２）であり、列数はｃｏｌｓである。

（Ｓ３０）行列積演算部１２６は、ステップＳ２７で、行管理テーブル１４２の全ての行番号ｒｏｗを選択したか判断する。全ての行番号ｒｏｗを選択した場合はステップＳ３１に処理が進み、未選択の行番号ｒｏｗがある場合はステップＳ２７に処理が進む。

図１１は、第１の行列積演算の手順例を示すフローチャート（続き）である。
（Ｓ３１）行列積演算部１２６は、イテレーション数Ｉ＝０に初期化する。
（Ｓ３２）行列積演算部１２６は、スレッドテーブル１４６からスレッド番号を１つ選択する（Ｈ［ｒｏｗ］［ｍ］［ｃｏｌ］）。スレッド番号の位置は、行番号ｒｏｗと二次元配列内の行インデックスｍと二次元配列内の列番号ｃｏｌによって特定される。選択されたスレッド番号が示すスレッドによってステップＳ３３〜Ｓ３５が実行される。

（Ｓ３３）行列積演算部１２６は、行管理テーブル１４２から行番号ｒｏｗに対応する組番号Ｓ_Ｒ［ｒｏｗ］と行番号ｒｏｗ＋１に対応する組番号Ｓ_Ｒ［ｒｏｗ＋１］を検索する。行列積演算部１２６は、Ｓ_Ｒ［ｒｏｗ＋１］−Ｓ_Ｒ［ｒｏｗ］＞２^Ｉであるか、すなわち、行番号ｒｏｗに対応する組数が２^Ｉを超えるか判断する。上記条件を満たす場合はステップＳ３４に処理が進み、上記条件を満たさない場合はステップＳ３６に処理が進む。

（Ｓ３４）行列積演算部１２６は、行インデックスｍを２^Ｉで割った余り（剰余）がゼロであるか、すなわち、ｍ％２^Ｉ＝０であるか判断する。上記条件を満たす場合はステップＳ３５に処理が進み、上記条件を満たさない場合はステップＳ３６に処理が進む。

（Ｓ３５）行列積演算部１２６は、中間行列１４４から要素Ｔ［Ｓ_Ｒ［ｒｏｗ］＋２^Ｉ×２ｍ］［ｃｏｌ］と要素Ｔ［Ｓ_Ｒ［ｒｏｗ］＋２^Ｉ×（２ｍ＋１）］［ｃｏｌ］を選択し、後者の値を前者に加算することで中間行列１４４を更新する。

（Ｓ３６）行列積演算部１２６は、ステップＳ３２で、スレッドテーブル１４６の全てのスレッド番号を選択したか判断する。全てのスレッド番号を選択した場合はステップＳ３７に処理が進み、未選択のスレッド番号がある場合はステップＳ３２に処理が進む。なお、ステップＳ３３〜Ｓ３５の処理は並列処理化することができる。

（Ｓ３７）行列積演算部１２６は、イテレーション数Ｉを１だけ大きくする。
（Ｓ３８）行列積演算部１２６は、Ｉ＜ｃｅｉｌ（ｌｏｇ_２Ｍ）を満たすか判断する。ｃｅｉｌは天井関数であり、正数に対しては小数点以下切り上げを表す。ＭはステップＳ２８によって算出された非ゼロ要素数の最大値である。上記条件を満たす場合はステップＳ３２に処理が進み、上記条件を満たさない場合はステップＳ３９に処理が進む。

（Ｓ３９）行列積演算部１２６は、出力行列１４５（出力行列Ｏ）から行番号と列番号を選択する（行番号ｒｏｗと列番号ｃｏｌ）。行列積演算部１２６は、行管理テーブル１４２から行番号ｒｏｗに対応する組番号Ｓ_Ｒ［ｒｏｗ］を検索し、中間行列１４４からＴ［Ｓ_Ｒ［ｒｏｗ］］［ｃｏｌ］を抽出する。行列積演算部１２６は、Ｔ［Ｓ_Ｒ［ｒｏｗ］］［ｃｏｌ］の値をＯ［ｒｏｗ］［ｃｏｌ］の値として用いる。行列積演算部１２６は、各行と各列についてこれを繰り返すことで出力行列１４５を生成する。

以上、第１の行列積演算について説明した。第１の行列積演算では、中間行列１４４に対する集計処理の並列処理化が複雑になり、多数のスレッドを並列に実行可能なＧＰＧＰＵ１０２の演算能力が十分に活用されないおそれがある。これに対し、行列演算ライブラリ１２２は、圧縮行格納法を拡張した方法に基づく第２の行列積演算を実行することもできる。次に、第２の行列積演算について説明する。

図１２は、第２のＣＳＲデータの例を示す図である。
第２の行列積演算では、第１のＣＳＲデータに代えて、疎行列１３１から値配列１５１（値配列Ｓ_Ｖ）および列番号配列１５２（列番号配列Ｓ_Ｃ）が生成される。値配列１５１および列番号配列１５２は、データ記憶部１２９に記憶される。

値配列１５１は、疎行列１３１の非ゼロ要素の値を含む二次元配列である。列番号配列１５２は、疎行列１３１の非ゼロ要素が位置する列の列番号を含む二次元配列である。値配列１５１の行数と列番号配列１５２の行数は同じであり、値配列１５１の列数と列番号配列１５２の列数は同じである。同じ位置にある値配列１５１の要素と列番号配列１５２の要素とは対応関係にある。値配列１５１および列番号配列１５２の行数は、疎行列１３１の行数である。値配列１５１および列番号配列１５２の列数は、疎行列１３１の各行の非ゼロ要素数のうちの最大値（前述の最大値Ｍ）である。

値配列１５１では、疎行列１３１の各行の非ゼロ要素の値が左詰めで格納される。非ゼロ要素数が最大値未満の行については値がゼロであるダミー要素を挿入することで、値配列１５１の１行当り要素数を固定値に統一している。図４の疎行列１３１の場合、非ゼロ要素数の最大値は４であるため、値配列１５１は８行４列の二次元配列になる。

疎行列１３１の行＃０は非ゼロ要素数が２であるため、値配列１５１の列＃０，＃１に非ゼロ要素の値が登録され、列＃２，＃３にゼロが登録される。行＃１は非ゼロ要素数が２であるため、列＃０，＃１に非ゼロ要素の値が登録され、列＃２，＃３にゼロが登録される。行＃２は非ゼロ要素数が０であるため、列＃０〜＃３にゼロが登録される。行＃３は非ゼロ要素数が４であるため、列＃０〜＃３に非ゼロ要素の値が登録される。行＃４は非ゼロ要素数が４であるため、列＃０〜＃３に非ゼロ要素の値が登録される。行＃５は非ゼロ要素数が２であるため、列＃０，＃１に非ゼロ要素の値が登録され、列＃２，＃３にゼロが登録される。行＃６は非ゼロ要素数が１であるため、列＃０に非ゼロ要素の値が登録され、列＃１〜＃３にゼロが登録される。行＃７は非ゼロ要素数が２であるため、列＃０，＃１に非ゼロ要素の値が登録され、列＃２，＃３にゼロが登録される。

列番号配列１５２では、値配列１５１に非ゼロ要素の値が登録されている位置には、当該非ゼロ要素が存在する列の列番号が登録される。一方、値配列１５１にゼロが登録されている位置には、値配列１５１と同様にダミー要素が登録される。これにより、列番号配列１５２の１行当り要素数を固定値に統一している。列番号配列１５２のダミー要素は、ゼロなどの所定の列番号をもつ。ただし、ダミー要素の列番号は疎行列１３１に実在する列の列番号であればよく、任意の列番号でもよい。以下の説明では、ダミーの列番号の例としてゼロを用いている。図４の疎行列１３１の場合、値配列１５１と同様に列番号配列１５２は８行４列の二次元配列になる。

なお、第２の実施の形態では値配列１５１を二次元配列としたが、行番号の小さい順に非ゼロ要素およびダミー要素を並べた一次元配列とすることも可能である。同様に、列番号配列１５２を一次元配列とすることも可能である。また、第２の実施の形態では値配列１５１と列番号配列１５２を分離しているが、両者を単一のテーブルに統合してもよい。

図１３は、第２のＣＳＲデータを用いた行列積演算の例を示す図である。
疎行列１３１が値配列１５１および列番号配列１５２として表現されている場合、行列積演算部１２６は以下のようにして行列積を計算することができる。行列積の計算は第１のＣＳＲデータを用いた第１の行列積演算と同様に、第１工程としての複製処理と、第２の工程としての乗算処理と、第３工程としての集計処理とを含む。

複製処理では、行列積演算部１２６は、列番号配列１５２から中間テンソル１５３を生成する。中間テンソル１５３は、第１次元座標（Ｘ座標）と第２次元座標（Ｙ座標）と第３次元座標（Ｚ座標）により要素が特定される三次元配列である。中間テンソル１５３はデータ記憶部１２９に記憶される。Ｘ座標数は列番号配列１５２の行数と同じ、すなわち、疎行列１３１の行数と同じである。Ｙ座標数は列番号配列１５２の列数と同じ、すなわち、非ゼロ要素数の最大値と同じである。Ｚ座標数は密行列１３２の列数と同じである。

行列積演算部１２６は、列番号配列１５２から列番号を抽出し、その列番号と同じ行番号をもつ行を密行列１３２から抽出して中間テンソル１５３に格納する。このとき、列番号配列１５２の行と中間テンソル１５３のＸ座標が対応し、列番号配列１５２の列と中間テンソル１５３のＹ座標が対応するようにする。また、列番号配列１５２に含まれるダミーの列番号に対しても、他の列番号と同様の処理が行われる。

例えば、中間テンソル１５３にはＸ座標とＹ座標の組が８×４＝３２通り存在する。このうち１７通りについては、列番号配列１５２に含まれるダミーでない列番号に従って密行列１３２の何れかの行の複製が格納される。一方、残りの１５通りについては、列番号配列１５２に含まれるダミーの列番号に従って密行列１３２の特定の行の複製が格納される。第２の実施の形態ではダミーの列番号はゼロであるため、非ゼロ要素が存在しない位置には密行列１３２の行＃０が格納されることになる。

乗算処理では、行列積演算部１２６は、中間テンソル１５３から中間テンソル１５４を生成する。中間テンソル１５４は、中間テンソル１５３を更新したものであり、中間テンソル１５３の記憶領域を上書きすることで別個の記憶領域を使用しなくてもよい。

行列積演算部１２６は、値配列１５１から値を抽出し、その値を中間テンソル１５３の対応する要素に乗じる。値配列１５１に含まれるダミー要素に対しても、非ゼロ要素と同様の処理が行われる。例えば、行列積演算部１２６は、値配列１５１の０行０列に対応する値「１」を、中間テンソル１５３のＸ＝０，Ｙ＝０，Ｚ＝０〜２の各要素に乗じる（各要素の値を１倍する）。また、行列積演算部１２６は、値配列１５１の０行１列に対応する値「２」を、中間テンソル１５３のＸ＝０，Ｙ＝１，Ｚ＝０〜２の各要素に乗じる（各要素の値を２倍する）。中間テンソル１５３に存在する３２通りのＸ座標とＹ座標の組のうち１７通りについては、非ゼロ要素の値を乗じることになる。一方、残りの１５通りについては、ダミー要素の値であるゼロを乗じることになる。

集計処理では、行列積演算部１２６は、中間テンソル１５４の要素をＸ座標とＺ座標の組毎に合計することで出力行列１５５を生成する。これは、中間テンソル１５４のＸ座標毎に固定数の行ベクトルを合算することを意味する。出力行列１５５は、第１の行列積演算で生成される出力行列１５５と同じであり、データ記憶部１２９に記憶される。出力行列１５５は、行数が疎行列１３１の行数と同じであり、列数が密行列１３２の列数と同じ二次元行列である。出力行列１５５の行は中間テンソル１５４のＸ座標に対応し、出力行列１５５の列は中間テンソル１５４のＺ座標に対応する。値配列１５１のダミー要素から生成された中間テンソル１５４の要素に対しても、他の要素と同様の処理が行われる。

例えば、中間テンソル１５４のＸ＝０，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃０になる。中間テンソル１５４のＸ＝１，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃１になる。中間テンソル１５４のＸ＝２，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃２になる。中間テンソル１５４のＸ＝３，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃３になる。

また、中間テンソル１５４のＸ＝４，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃４になる。中間テンソル１５４のＸ＝５，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃５になる。中間テンソル１５４のＸ＝６，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃６になる。中間テンソル１５４のＸ＝７，Ｙ＝０〜３の４つの行ベクトルを合算して出力行列１５５の行＃７になる。このように、第１の行列積演算と異なり第２の行列積演算では、出力行列１５５の１行当りに合算する中間テンソル１５４の行ベクトルの数が固定になる。

ここで、中間テンソル１５３の生成は複数のスレッドを用いて容易に並列処理化できる。例えば、列番号配列１５２の異なる要素に対して異なるスレッドを割り当てることで、中間テンソル１５３の異なるＸ座標とＹ座標の組に対する複製処理を並列に実行できる。また、中間テンソル１５４の生成も複数のスレッドを用いて容易に並列処理化できる。例えば、値配列１５１の異なる要素に対して異なるスレッドを割り当てることで、中間テンソル１５４の異なるＸ座標とＹ座標の組に対する乗算処理を並列に実行できる。また、出力行列１５５の生成も複数のスレッドを用いて並列処理化できる。

図１４は、第２の行列積演算における集計処理の例を示す図である。
行列積演算部１２６は、中間テンソル１５４の一部要素を順次上書きすることで集計処理を進める。行列積演算部１２６は、中間テンソル１５４を更新して中間テンソル１５４ａを生成し、中間テンソル１５４ａを更に更新して中間テンソル１５４ｂを生成する。中間テンソル１５４ａ，１５４ｂのために別途記憶領域を用意しなくてもよい。

行列積演算部１２６は、中間テンソル１５４の各Ｘ座標の中で、Ｙ座標が異なる２つの行ベクトルの合算を二分木形式で階層的に繰り返す。２つの行ベクトルの合算は、Ｚ座標毎に、Ｙ座標が大きい方の要素の値をＹ座標が小さい方の要素に加算することで行う。最初はＹ座標が１だけ離れた行ベクトル同士が合算され、集計処理のステップが進むにつれてＹ座標が離れた行ベクトル同士が合算される。集計処理の第Ｉステップ（Ｉは０以上の整数）では、Ｙ座標が２^Ｉだけ離れた行ベクトル同士が合算される。中間テンソル１５４のＹ座標数をＭとすると、集計処理の終了までのステップ数はｌｏｇ_２Ｍである。

例えば、中間テンソル１５４を中間テンソル１５４ａに更新する第０ステップでは、各Ｘ座標について、Ｙ＝１，Ｚ＝０の値がＹ＝０，Ｚ＝０に加算され、Ｙ＝１，Ｚ＝１の値がＹ＝０，Ｚ＝１に加算され、Ｙ＝１，Ｚ＝２の値がＹ＝０，Ｚ＝２に加算される。また、各Ｘ座標について、Ｙ＝３，Ｚ＝０の値がＹ＝２，Ｚ＝０に加算され、Ｙ＝３，Ｚ＝１の値がＹ＝２，Ｚ＝１に加算され、Ｙ＝３，Ｚ＝２の値がＹ＝２，Ｚ＝２に加算される。

次に、中間テンソル１５４ａを中間テンソル１５４ｂに更新する第１ステップでは、各Ｘ座標について、Ｙ＝２，Ｚ＝０の値がＹ＝０，Ｚ＝０に加算され、Ｙ＝２，Ｚ＝１の値がＹ＝０，Ｚ＝１に加算され、Ｙ＝２，Ｚ＝２の値がＹ＝０，Ｚ＝２に加算される。上記以外の要素は、第０ステップで合算済みであるため無視してよい。

この例では中間テンソル１５４のＹ座標数が４であるため、第０ステップと第１ステップで集計処理は終了する。集約処理の結果は、中間テンソル１５４ｂのＹ＝０に位置する行ベクトルである。行列積演算部１２６は、中間テンソル１５４ｂからＹ＝０の行ベクトルを抽出して出力行列１５５を生成する。中間テンソル１５４ｂの一部の行ベクトルを複製して出力行列１５５を生成してもよいし、行ベクトルの複製を行わずに中間テンソル１５４ｂの一部の行ベクトルのみが見えるビューを生成してもよい。

例えば、出力行列１５５の行＃０は中間テンソル１５４ｂのＸ＝０，Ｙ＝０である。出力行列１５５の行＃１は中間テンソル１５４ｂのＸ＝１，Ｙ＝０である。出力行列１５５の行＃２は中間テンソル１５４ｂのＸ＝２，Ｙ＝０である。出力行列１５５の行＃３は中間テンソル１５４ｂのＸ＝３，Ｙ＝０である。出力行列１５５の行＃４は中間テンソル１５４ｂのＸ＝４，Ｙ＝０である。出力行列１５５の行＃５は中間テンソル１５４ｂのＸ＝５，Ｙ＝０である。出力行列１５５の行＃６は中間テンソル１５４ｂのＸ＝６，Ｙ＝０である。出力行列１５５の行＃７は中間テンソル１５４ｂのＸ＝７，Ｙ＝０である。

行列積演算部１２６は、２つの値を加算する加算演算毎に１つのスレッドを割り当てる。上記の第０ステップでは８×２×３＝４８個の加算演算が行われており、４８個のスレッドが並列に実行される。また、上記の第１ステップでは８×１×３＝２４個の加算演算が行われており、２４個のスレッドが並列に実行される。ただし、集約処理が進むにつれて加算演算の数は減少するため、前のステップで使用されたスレッドのサブセットを次のステップで使用すればよい。上記の第１ステップでは、中間テンソル１５４ａのＸ＝０，Ｙ＝０，Ｚ＝０を計算したスレッドが中間テンソル１５４ｂのＸ＝０，Ｙ＝０，Ｚ＝０を計算すればよい。よって、スレッド割り当ては集約処理の開始時に決定される。

図１５は、第２の集計処理におけるスレッド割り当て例を示す図である。
スレッド配列１５６は、中間テンソル１５４から決定されるスレッド割り当てを示す。スレッド配列１５６は、データ記憶部１２９に記憶される。スレッド配列１５６は、第１次元座標（Ｘ座標）と第２次元座標（Ｙ座標）と第３次元座標（Ｚ座標）により要素が特定される三次元配列である。スレッド配列１５６のＸ座標数は、中間テンソル１５４のＸ座標数と同じである。スレッド配列１５６のＹ座標数は、中間テンソル１５４のＹ座標数をＭとするとＭ／２（小数点以下切り捨て）である。スレッド配列１５６のＺ座標数は、中間テンソル１５４のＺ座標数と同じである。

スレッド配列１５６には、スレッドを識別するスレッド番号が登録される。中間テンソル１５４の各Ｘ座標に対して固定数のスレッドが割り当てられる。上記のように、中間テンソル１５４に対する集計処理の第０ステップでは４８個の加算演算が行われるため、４８個のスレッドが割り当てられる。例えば、Ｘ＝０に対してスレッド＃０〜＃５、Ｘ＝１に対してスレッド＃６〜＃１１、Ｘ＝２に対してスレッド＃１２〜＃１７、Ｘ＝３に対してスレッド＃１８〜＃２３が割り当てられる。また、Ｘ＝４に対してスレッド＃２４〜＃２９、Ｘ＝５に対してスレッド＃３０〜＃３５、Ｘ＝６に対してスレッド＃３６〜＃４１、Ｘ＝７に対してスレッド＃４２〜＃４７が割り当てられる。

スレッド＃０は、中間テンソル１５４の［Ｘ，Ｙ，Ｚ］＝［０，０，０］に［０，１，０］の値を加算するものである。スレッド＃１は、［０，０，１］に［０，１，１］の値を加算するものである。スレッド＃２は、［０，０，２］に［０，１，２］の値を加算するものである。スレッド＃３は、［０，２，０］に［０，３，０］の値を加算するものである。スレッド＃４は、［０，２，１］に［０，３，１］の値を加算するものである。スレッド＃５は、［０，２，２］に［０，３，２］の値を加算するものである。

集計処理の第０ステップでは、これらのスレッドが全て並列に実行される。それより後のステップでは、これらのスレッドの一部が並列に実行される。上記のように、中間テンソル１５４に対する集計処理の第１ステップでは２４個の加算演算が行われるため、スレッド＃０〜＃４７のうちの２４個のスレッドが並列に実行される。例えば、スレッド＃０〜＃２，＃６〜＃８，＃１２〜＃１４，＃１８〜＃２０，＃２４〜＃２６，＃３０〜＃３２，＃３６〜＃３８，＃４２〜＃４４が実行され、それ以外のスレッドは実行されない。

スレッド＃０は、中間テンソル１５４ａの［Ｘ，Ｙ，Ｚ］＝［０，０，０］に［０，２，０］の値を加算する。スレッド＃１は、［０，０，１］に［０，２，１］の値を加算する。スレッド＃２は、［０，０，２］に［０，２，２］の値を加算する。

第１の行列積演算では集計処理に２４個のスレッドが使用されるのに対し、第２の行列積演算では４８個のスレッドが使用されておりスレッド数が増加している。また、第２の行列積演算ではゼロの乗算やゼロの加算など、ダミー要素に基づく余分な演算が追加されている。一方で、第２の行列積演算では中間テンソル１５４のデータ範囲とスレッドとの対応付けが容易であり、並列処理の制御が簡潔になる。よって、ＧＰＧＰＳ１０２が多数のスレッドを並列実行可能であれば、並列処理が効率化されて実行時間が短縮する。

次に、第２のＣＳＲデータに関する処理手順を説明する。
図１６は、第２のＣＳＲデータ生成の手順例を示すフローチャートである。
第２のＣＳＲデータ生成は、ユーザアプリケーション１２１からの入力に応じて疎行列生成部１２３が二次元構造データの疎行列１３１を生成した後に行われる。

（Ｓ４０）データ構造変換部１２４は、最大値Ｍ＝０に初期化する。
（Ｓ４１）データ構造変換部１２４は、疎行列１３１（疎行列Ｓ_２Ｄ）から行番号を１つ選択する（行番号ｒｏｗ）。

（Ｓ４２）データ構造変換部１２４は、非ゼロ要素数ｍ＝０に初期化する。
（Ｓ４３）データ構造変換部１２４は、上記の行番号ｒｏｗについて、疎行列１３１から列番号を１つ選択する（列番号ｃｏｌ）。

（Ｓ４４）データ構造変換部１２４は、疎行列１３１からｒｏｗ行ｃｏｌ列の値Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］（値ｖａｌ）を抽出し、ｖａｌ＝０であるか判断する。ｖａｌ＝０である場合、すなわち、Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］がゼロ要素である場合、ステップＳ４６に処理が進む。ｖａｌ＝０でない場合、すなわち、Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］が非ゼロ要素である場合、ステップＳ４５に処理が進む。

（Ｓ４５）データ構造変換部１２４は、非ゼロ要素数ｍを１だけ大きくする。
（Ｓ４６）データ構造変換部１２４は、ステップＳ４３で、行番号ｒｏｗについて全ての列番号ｃｏｌを選択したか判断する。全ての列番号ｃｏｌを選択した場合はステップＳ４７に処理が進み、未選択の列番号ｃｏｌがある場合はステップＳ４３に処理が進む。

（Ｓ４７）データ構造変換部１２４は、行番号ｒｏｗの非ゼロ要素数ｍと最大値Ｍとを比較し、ｍがＭより大きいか判断する。ｍがＭより大きい場合はステップＳ４８に処理が進み、ｍがＭ以下である場合はステップＳ４９に処理が進む。

（Ｓ４８）データ構造変換部１２４は、最大値Ｍに非ゼロ要素数ｍを代入する。
（Ｓ４９）データ構造変換部１２４は、ステップＳ４１で、全ての行番号ｒｏｗを選択したか判断する。全ての行番号ｒｏｗを選択した場合はステップＳ５０に処理が進み、未選択の行番号ｒｏｗがある場合はステップＳ４１に処理が進む。

（Ｓ５０）データ構造変換部１２４は、行数が疎行列１３１と同じであり列数がＭの値配列１５１（値配列Ｓ_Ｖ）を生成する。また、データ構造変換部１２４は、行数が疎行列１３１と同じであり列数がＭの列番号配列１５２（列番号配列Ｓ_Ｃ）を生成する。データ構造変換部１２４は、値配列１５１と列番号配列１５２の各要素をゼロに初期化する。

図１７は、第２のＣＳＲデータ生成の手順例を示すフローチャート（続き）である。
（Ｓ５１）データ構造変換部１２４は、疎行列１３１から行番号を１つ選択する（行番号ｒｏｗ）。

（Ｓ５２）データ構造変換部１２４は、非ゼロ要素数ｍ＝０に初期化する。
（Ｓ５３）データ構造変換部１２４は、上記の行番号ｒｏｗについて、疎行列１３１から列番号を１つ選択する（列番号ｃｏｌ）。

（Ｓ５４）データ構造変換部１２４は、疎行列１３１からｒｏｗ行ｃｏｌ列の値Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］（値ｖａｌ）を抽出し、ｖａｌ＝０であるか判断する。ｖａｌ＝０である場合、すなわち、Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］がゼロ要素である場合、ステップＳ５７に処理が進む。ｖａｌ＝０でない場合、すなわち、Ｓ_２Ｄ［ｒｏｗ］［ｃｏｌ］が非ゼロ要素である場合、ステップＳ５５に処理が進む。

（Ｓ５５）データ構造変換部１２４は、値配列１５１のｒｏｗ行ｍ列に値ｖａｌを代入する（Ｓ_Ｖ［ｒｏｗ］［ｍ］＝ｖａｌ）。また、データ構造変換部１２４は、列番号配列１５２のｒｏｗ行ｍ列に列番号ｃｏｌを代入する（Ｓ_Ｃ［ｒｏｗ］［ｍ］＝ｃｏｌ）。

（Ｓ５６）データ構造変換部１２４は、非ゼロ要素数ｍを１だけ大きくする。
（Ｓ５７）データ構造変換部１２４は、ステップＳ５３で、行番号ｒｏｗについて全ての列番号ｃｏｌを選択したか判断する。全ての列番号ｃｏｌを選択した場合はステップＳ５８に処理が進み、未選択の列番号ｃｏｌがある場合はステップＳ５３に処理が進む。

（Ｓ５８）データ構造変換部１２４は、ステップＳ５１で、全ての行番号ｒｏｗを選択したか判断する。全ての行番号ｒｏｗを選択した場合はＣＳＲデータ生成が終了し、未選択の行番号ｒｏｗがある場合はステップＳ５１に処理が進む。

図１８は、第２の行列積演算の手順例を示すフローチャートである。
第２の行列積演算は、データ構造変換部１２４が第２のＣＳＲデータを生成し、密行列生成部１２５が二次元構造データの密行列１３２を生成した後に行われる。

（Ｓ６０）行列積演算部１２６は、値配列１５１の行数ｒｏｗｓ、値配列１５１の列数Ｍおよび密行列１３２の列数ｃｏｌｓを確認する。行列積演算部１２６は、大きさがｒｏｗｓ×Ｍ×ｃｏｌｓの中間テンソル１５３（中間テンソルＴ）を生成する。

（Ｓ６１）行列積演算部１２６は、列番号配列１５２の行番号と列番号を１つずつ選択し（行番号ｒｏｗと列番号ｍ）、列番号配列１５２からｒｏｗ行ｍ列の要素Ｓ_Ｃ［ｒｏｗ］［ｍ］（列番号ｃｏｌ）を選択する。

（Ｓ６２）行列積演算部１２６は、密行列１３２（密行列Ｄ）のｃｏｌ行目をコピーし、中間テンソル１５３のＸ＝ｒｏｗ，Ｙ＝ｍに格納する（Ｔ［ｒｏｗ］［ｍ］［：］）。なお、列番号ｃｏｌがダミーであっても通常通りコピーが行われる。

（Ｓ６３）行列積演算部１２６は、ステップＳ６１で、列番号配列１５２の全ての要素を選択したか判断する。全ての要素を選択した場合はステップＳ６４に処理が進み、未選択の要素がある場合はステップＳ６１に処理が進む。なお、ステップＳ６１，Ｓ６２の処理は複数のスレッドを用いて並列処理化することができる。

（Ｓ６４）行列積演算部１２６は、値配列１５１の行番号と列番号を１つずつ選択し（行番号ｒｏｗと列番号ｍ）、値配列１５１からｒｏｗ行ｍ列の要素Ｓ_Ｖ［ｒｏｗ］［ｍ］（値ｖａｌ）を選択する。

（Ｓ６５）行列積演算部１２６は、中間テンソル１５３のＸ＝ｒｏｗ，Ｙ＝ｍの各要素（Ｔ［ｒｏｗ］［ｍ］［：］）に対して値ｖａｌを乗じる。なお、値ｖａｌがダミーでありｖａｌ＝０であっても通常通り乗算が行われる。

（Ｓ６６）行列積演算部１２６は、行列積演算部１２６は、ステップＳ６４で、値配列１５１の全ての要素を選択したか判断する。全ての要素を選択した場合はステップＳ６７に処理が進み、未選択の要素がある場合はステップＳ６４に処理が進む。なお、ステップＳ６４，Ｓ６５の処理は複数のスレッドを用いて並列処理化することができる。

（Ｓ６７）行列積演算部１２６は、スレッドプール１２８からｒｏｗｓ×ｆｌｏｏｒ（Ｍ／２）×ｃｏｌｓ個のスレッドを取得する。ｒｏｗｓは、中間テンソル１５３を更新することで得られた中間テンソル１５４（中間テンソルＴ）のＸ座標数に相当する。Ｍは、中間テンソル１５４のＹ座標数に相当する。ｃｏｌｓは、中間テンソル１５４のＺ座標数に相当する。行列積演算部１２６は、スレッド配列１５６（スレッド配列Ｈ）に、取得したスレッドのスレッド番号を登録する（Ｈ［：］［：］［：］）。スレッド配列１５６のＸ座標数はｒｏｗ、Ｙ座標数はｆｌｏｏｒ（Ｍ／２）、Ｚ座標数はｃｏｌｓである。

図１９は、第２の行列積演算の手順例を示すフローチャート（続き）である。
（Ｓ６８）行列積演算部１２６は、イテレーション数Ｉ＝０に初期化する。
（Ｓ６９）行列積演算部１２６は、スレッド配列１５６からスレッド番号を１つ選択する（Ｈ［ｒｏｗ］［ｍ］［ｃｏｌ］）。選択されたスレッド番号の位置は、Ｘ座標ｒｏｗとＹ座標ｍとＺ座標ｃｏｌによって特定される。選択されたスレッド番号が示すスレッドによって、以下のステップＳ７０，Ｓ７１が実行される。

（Ｓ７０）行列積演算部１２６は、ｍを２^Ｉで割った余り（剰余）がゼロであるか、すなわち、ｍ％２^Ｉ＝０であるか判断する。上記条件を満たす場合はステップＳ７１に処理が進み、上記条件を満たさない場合はステップＳ７２に処理が進む。

（Ｓ７１）行列積演算部１２６は、中間テンソル１５４から要素Ｔ［ｒｏｗ］［２^Ｉ×２ｍ］［ｃｏｌ］と要素Ｔ［ｒｏｗ］［２^Ｉ×（２ｍ＋１）］［ｃｏｌ］を選択し、後者の値を前者に加算することで中間テンソル１５４を更新する。なお、前者および後者の少なくとも一方がダミーであり値がゼロであっても通常通り加算が行われる。

（Ｓ７２）行列積演算部１２６は、ステップＳ６９で、スレッド配列１５６の全てのスレッド番号を選択したか判断する。全てのスレッド番号を選択した場合はステップＳ７３に処理が進み、未選択のスレッド番号がある場合はステップＳ６９に処理が進む。なお、ステップＳ７０，Ｓ７１の処理は並列処理化することができる。

（Ｓ７３）行列積演算部１２６は、イテレーション数Ｉを１だけ大きくする。
（Ｓ７４）行列積演算部１２６は、Ｉ＜ｃｅｉｌ（ｌｏｇ_２Ｍ）を満たすか判断する。Ｍは中間テンソル１５４のＹ座標数である。上記条件を満たす場合はステップＳ６９に処理が進み、上記条件を満たさない場合はステップＳ７５に処理が進む。

（Ｓ７５）行列積演算部１２６は、中間テンソル１５４からＹ＝０のデータ範囲を抽出して出力行列１５５（出力行列Ｏ）を生成する。すなわち、行列積演算部１２６は、中間テンソル１５４の全てのＸ座標ｒｏｗとＺ座標ｃｏｌについて、Ｏ［ｒｏｗ］［ｃｏｌ］＝Ｔ［ｒｏｗ］［０］［ｃｏｌ］とする。

第２の実施の形態の情報処理装置１００によれば、大規模疎行列が圧縮格納データによって表現され、大規模疎行列と密行列との行列積が圧縮格納データのまま実行される。よって、メモリ使用量が削減されると共に行列積の計算量が削減される。また、行列積演算がデータロードや乗算や加算などの単位演算に細分化されて多数のスレッドに割り振られ、多数の演算器を有するプロセッサを用いてそれら多数のスレッドが並列に実行される。よって、行列積演算を高速に実行することができる。

また、第２の行列積演算の方法を採用した場合、大規模疎行列の１行当たり非ゼロ要素数が可変であっても、ダミー要素の挿入によって圧縮格納データに含まれる１行当たり要素数が固定化される。よって、細分化した単位演算を複数のスレッドに割り振る制御が簡潔になり並列処理が効率化される。このため、プロセッサが有する多数の演算器を有効に活用することができ、行列積演算の実行時間を短縮することができる。

１０ 行列演算装置
１１ 記憶部
１２ 処理部
１３ 行列演算プログラム
１４ａ，１４ｂ スレッド
１５，１６，１７ 行列
１８ 圧縮格納データ
１９ ベクトルデータ

Claims (8)

1. 行列演算プログラムを記憶する記憶部と、
   前記行列演算プログラムに基づいて複数のスレッドを並列に実行可能な処理部と、
   を有し、
   前記行列演算プログラムを実行する前記処理部は、
   第１の行列に含まれる複数の第１の行それぞれについて値がゼロでない非ゼロ要素の数をカウントし、前記複数の第１の行の間で前記非ゼロ要素の数の最大値を判定し、
   前記複数の第１の行それぞれから非ゼロ要素の値と当該非ゼロ要素が位置する列を示す列識別子とのペアを抽出し、前記非ゼロ要素の数が前記最大値より少ない第１の行に対しては値がゼロであるダミーのペアを追加することで、前記複数の第１の行それぞれに対して共通する個数のペアを含む圧縮格納データを生成し、
   前記圧縮格納データに含まれるペアそれぞれに対して、第２の行列から当該ペアの列識別子に対応する行識別子をもつ第２の行を抽出し、当該抽出した第２の行に対して当該ペアの値を乗算することで、当該ペアに対応する行ベクトルを生成し、
   前記複数の第１の行それぞれに対して共通する個数のスレッドを割り当て、前記複数の第１の行それぞれについて前記共通する個数のスレッドを用いて行ベクトルを合算することで、前記第１の行列と前記第２の行列との行列積を示す第３の行列を生成する、
   行列演算装置。
2. 前記ダミーのペアは、ゼロの値と所定の列識別子とのペアであり、
   前記行ベクトルの生成は、前記ダミーのペアに対しても実行される、
   請求項１記載の行列演算装置。
3. 前記圧縮格納データに含まれる第１の行毎のペアの個数は、前記最大値である、
   請求項１記載の行列演算装置。
4. 前記複数の第１の行それぞれに割り当てるスレッドの個数は、前記圧縮格納データに含まれる第１の行毎のペアの個数と前記第２の行列の列数とから決定される、
   請求項１記載の行列演算装置。
5. 前記複数の第１の行それぞれに対応して４つ以上の行ベクトルが生成された場合、前記行ベクトルの合算は、２つの行ベクトルの合算を木構造で繰り返すことで行う、
   請求項１記載の行列演算装置。
6. 行列演算プログラムを記憶する記憶部と、
   前記行列演算プログラムに基づいて複数のスレッドを並列に実行可能な処理部と、
   を有し、
   前記行列演算プログラムを実行する前記処理部は、
   第１の行列に含まれる複数の第１の列それぞれについて値がゼロでない非ゼロ要素の数をカウントし、前記複数の第１の列の間で前記非ゼロ要素の数の最大値を判定し、
   前記複数の第１の列それぞれから非ゼロ要素の値と当該非ゼロ要素が位置する行を示す行識別子とのペアを抽出し、前記非ゼロ要素の数が前記最大値より少ない第１の列に対しては値がゼロであるダミーのペアを追加することで、前記複数の第１の列それぞれに対して共通する個数のペアを含む圧縮格納データを生成し、
   前記圧縮格納データに含まれるペアそれぞれに対して、第２の行列から当該ペアの行識別子に対応する列識別子をもつ第２の列を抽出し、当該抽出した第２の列に対して当該ペアの値を乗算することで、当該ペアに対応する列ベクトルを生成し、
   前記複数の第１の列それぞれに対して共通する個数のスレッドを割り当て、前記複数の第１の列それぞれについて前記共通する個数のスレッドを用いて列ベクトルを合算することで、前記第２の行列と前記第１の行列との行列積を示す第３の行列を生成する、
   行列演算装置。
7. コンピュータが実行する行列演算方法であって、
   第１の行列に含まれる複数の第１の行それぞれについて値がゼロでない非ゼロ要素の数をカウントし、前記複数の第１の行の間で前記非ゼロ要素の数の最大値を判定し、
   前記複数の第１の行それぞれから非ゼロ要素の値と当該非ゼロ要素が位置する列を示す列識別子とのペアを抽出し、前記非ゼロ要素の数が前記最大値より少ない第１の行に対しては値がゼロであるダミーのペアを追加することで、前記複数の第１の行それぞれに対して共通する個数のペアを含む圧縮格納データを生成し、
   前記圧縮格納データに含まれるペアそれぞれに対して、第２の行列から当該ペアの列識別子に対応する行識別子をもつ第２の行を抽出し、当該抽出した第２の行に対して当該ペアの値を乗算することで、当該ペアに対応する行ベクトルを生成し、
   前記複数の第１の行それぞれに対して共通する個数のスレッドを割り当て、前記複数の第１の行それぞれについて前記共通する個数のスレッドを用いて行ベクトルを合算することで、前記第１の行列と前記第２の行列との行列積を示す第３の行列を生成する、
   行列演算方法。
8. コンピュータに、
   第１の行列に含まれる複数の第１の行それぞれについて値がゼロでない非ゼロ要素の数をカウントし、前記複数の第１の行の間で前記非ゼロ要素の数の最大値を判定し、
   前記複数の第１の行それぞれから非ゼロ要素の値と当該非ゼロ要素が位置する列を示す列識別子とのペアを抽出し、前記非ゼロ要素の数が前記最大値より少ない第１の行に対しては値がゼロであるダミーのペアを追加することで、前記複数の第１の行それぞれに対して共通する個数のペアを含む圧縮格納データを生成し、
   前記圧縮格納データに含まれるペアそれぞれに対して、第２の行列から当該ペアの列識別子に対応する行識別子をもつ第２の行を抽出し、当該抽出した第２の行に対して当該ペアの値を乗算することで、当該ペアに対応する行ベクトルを生成し、
   前記複数の第１の行それぞれに対して共通する個数のスレッドを割り当て、前記複数の第１の行それぞれについて前記共通する個数のスレッドを用いて行ベクトルを合算することで、前記第１の行列と前記第２の行列との行列積を示す第３の行列を生成する、
   処理を実行させる行列演算プログラム。

Images