JP2016119084A - 効率的な疎行列表現及び処理のためのコンピュータ実装システム及び方法 - Google Patents

Abstract

【課題】元の構造を乱すことなく疎行列処理の速度を向上させるコンピュータ実装システム及び方法を提供する。【解決手段】疎行列が処理される速度は、行列の改善された圧縮表現を使用することによって増加させることができる。構造化圧縮表現は、キャッシュがランダムにアクセスされなければならない回数を削減することによって行列処理中に経験するキャッシュミス数を低減する。さらに、非ゼロエントリのそれらの数に基づいて行列の行及び列を分割して再グループ化する行列の表現は、行列のこれらの部分を処理するための最も適切なカーネル関数を割り当てるのを可能とし、ＧＰＵベースのハードウェアの限界を克服する。【選択図】図６

Description

本特許出願は、一般に、行列データを処理することに関し、特に、効率的な疎行列表現及び処理のためのコンピュータ実装システム及び方法に関する。

疎行列は、要素の大部分がゼロである行列である。そのような行列を使用する操作は、様々なアプリケーションを有し、通常、そのようなアプリケーションの最も計算集約的な部分である。例えば、疎行列−ベクトル乗算（ＳｐＭＶ）及び疎行列転置ベクトル乗算（ＳｐＭＴＶ）、疎線形代数（ＳＬＡ）における基本的な操作は、検索結果を提供するときにウェブページをランク付けするためのグーグル（登録商標）によって使用されているページランク（登録商標）アルゴリズムなどのランク付けアルゴリズムを実行するために使用されている。ＳｐＭＶ及びＳｐＭＴＶは、そのようなアプリケーションの最も計算集約的な部分であり、行列が使用可能な速度は、ＳｐＭＶ及びＳｐＭＴＶによって制限される。

疎行列処理の速度を向上させるための試みが行われてきたが、そのような努力は、まだ改善のための重要な余地を残している。例えば、速度を高めるために、行列は、非ゼロエントリの行列内の値及び位置に関する情報の複数の配列を含む圧縮フォーマットで符号化されており、ゼロエントリに関する情報を省略している。例えば、圧縮疎行列の行フォーマットは、非ゼロエントリの値を有する配列、非ゼロエントリが配置された列、及び、各行における第１の非ゼロエントリの第１の配列内のインデックスを保持する配列を含む。圧縮疎列フォーマットは、同様の配列を含む。そのような配列は、配列データに対する高速アクセスを可能にするために計算を実行するプロセッサのキャッシュに記憶される。しかしながら、より大きな行列の場合、圧縮フォーマット配列は、キャッシュに収まらないことがあり、プロセッサが計算の単一ステップを実行するためにメインメモリにおける行列を表す異なる配列にアクセスすることを必要とする。そのような構成において、中央処理装置（ＣＰＵ）及びグラフィックス処理ユニット（ＧＰＵ）を含む現代のコンピュータプロセッサは、計算中にキャッシュミスやキャッシュから必要なデータを取得するためにプロセッサによる障害を経験する可能性がある。キャッシュミスの後に計算を終了することは、プロセッサがメインメモリから失われたデータを取得することを必要とし、はるかに遅くなることがある。

さらに、ＧＰＵがＳｐＭＶ及びＳｐＭＴＶなどの行列演算を実行するために使用される場合に追加の懸念が存在する。ＧＰＵは、そのほとんどの要素が非ゼロエントリである密行列の処理などの密計算用により良好に設計されて最適化されている。そのようなハードウェアは、一般に、行列データを処理するために単一のカーネル関数を実行する。その結果、ハードウェアは、異なる行又は列などの行列の異なる部分における非ゼロエントリの数の大幅な変動に応答することができない。例えば、行列の単一の行又は列を処理するために単一のスレッドを割り当てるカーネルは、最も密な行又は列を処理するために割り当てられたスレッドに依存する全体の処理時間により、負荷不均衡を被ることがある。一方、単一の行又は列を処理するために複数のスレッドを割り当てるカーネルは、処理に関与しないいくつかの割り当てられたスレッドにより、割り当てられたスレッドの数が行又は列における非ゼロエントリの数未満である場合にハードウェアリソースの浪費を被る。

したがって、キャッシュミスの可能性を低減させ且つ行列の異なる部分における非ゼロエントリの数の変化に応答することができるような方法で疎行列を表現する必要がある。

疎行列が処理される速度は、行列の改善された圧縮表現を使用することによって増加させることができる。構造化圧縮表現は、キャッシュがランダムにアクセスされなければならない回数を削減することにより、行列処理中に経験したキャッシュミスの数を低減する。さらに、それらの非ゼロエントリ数に基づいて行列の行及び列を分割し且つ再グループ化するという行列の表現は、行列のこれらの部分を処理するための最も適切なカーネル関数を割り当てることができ、ＧＰＵベースハードウェアの限界を克服する。その結果、処理速度は、行列の元の構造を乱すことなく増加させることができる。

１つの実施形態は、構造化疎行列表現を取得するためのコンピュータ実装システム及び方法を提供する。各部分が行及び列の一方を含む行列の部分において１つ以上の順序で配置された１つ以上の非ゼロエントリを含む行列の構造化圧縮表現が取得され、各要素が非ゼロエントリのいずれかとその非ゼロエントリを含む部分のいずれかのインデックスとを含む１つ以上の要素を含む複合配列を取得することと、その順序のうちの１つ以上における最初である非ゼロエントリを含む各要素の複合配列におけるインデックスを含み、さらに行列における複数の非ゼロエントリを含むインデックス配列を取得することとを備える。

さらなる実施形態は、効率的な疎行列表現及び処理のためのコンピュータ実装システム及び方法を提供する。１つ以上の非ゼロエントリを含む疎行列の圧縮表現が取得され、各部分が行列の行及び列の一方を含む行列の部分が行列におけるそれらの部分に基づいてインデックス付けされる。行列の部分についての複数のパーティションが定義される。各部分における複数の非ゼロエントリは、圧縮表現を使用して取得される。各部分は、その部分における複数の非ゼロエントリに基づいてパーティションの１つに関連付けられる。各パーティションに関連する全ての部分のリストが作成され、部分はそれらのインデックスの順序でリスト化される。リストを含むマッピング配列を含む行列のパーティション化された圧縮表現が作成される。

ＧＰＵベースのハードウェアについて、ＧＰＵは、通常、ＣＰＵよりもはるかに柔軟性が低い単一命令複数スレッドの実行モデルを想定していることから、パーティション化された（構造化）ＣＳＲ／ＣＳＣ符号化は、非パーティション化のものよりも好ましい。このモデルは、ＧＰＵをＣＰＵよりも負荷不均衡の問題に対してより脆弱にする。ＰＳＣＳＲ及びＰＳＣＳＣなどのパーティション化された疎行列符号化は、非ゼロのそれらの数に基づいて同様の計算特性を有する行又は列をともにグループ化することによってＧＰＵに対する負荷不均衡を効果的に低減することができる。

本発明のさらに他の実施形態は、以下の詳細な説明から当業者にとって容易に明らかになるであろう。本発明を実施するための最良の形態を説明することによって本発明の実施形態が記載される。理解されるように、本発明は、他の及び異なる実施形態が可能であり、そのいくつかの詳細は、全て本発明の精神及び範囲から逸脱することなく様々な明白な点において変更が可能である。したがって、図面及び詳細な説明は、本質的に例示であり、限定されるものではないとみなされるべきである。

図１は、１つの実施形態（従来技術）にかかる、圧縮された疎行フォーマットの符号化において符号化された行列に対してＳｐＭＶを実行するための方法を示すフロー図である。 図２は、１つの実施形態（従来技術）にかかる、圧縮された疎列フォーマットの符号化において符号化された行列に対してＳｐＭＴＶを実行するための方法を示すフロー図である。 図３は、１つの実施形態にかかる、疎行列の効率的な表現及び処理のためのコンピュータ実装システムを示すブロック図である。 図４は、一例として隣接行列を示すグラフである。 図５は、図３のシステムの様々なハードウェアセットアップのための様々な符号化の強みをまとめたグラフである。 図６は、１つの実施形態にかかる構造化疎行列表現を取得するためのコンピュータ実装方法を示すフロー図である。 図７は、１つの実施形態にかかる図６の方法において使用するための構造化ＣＳＲ符号化を使用して符号化された行列に対してＳｐＭＶを実行するためのルーチンである。 図８は、１つの実施形態にかかる図６の方法において使用するための構造化ＣＳＣ符号化を使用して符号化された行列に対してＳｐＭＴＶを実行するためのルーチンである。 図９は、１つの実施形態にかかる効率的な疎行列パーティション及び処理のためのコンピュータ実装方法を示すフロー図である。 図１０は、１つの実施形態にかかる図９の方法において使用するための圧縮表現で表される行列の行の順序維持パーティションを実行するためのルーチンを示すフロー図である。 図１１は、１つの実施形態にかかる図９の方法において使用するための圧縮表現で表される行列の列の順序維持パーティションを実行するためのルーチンを示すフロー図である。 図１２は、１つの実施形態にかかる図９の方法において使用するためのマッピング配列にパーティションをマージするためのルーチンを示すフロー図である。 図１３は、１つの実施形態にかかる図９の方法において使用するためのパーティション化された圧縮表現で符号化された行列に対してＳｐＭＶを実行するためのルーチンを示すフロー図である。 図１４は、１つの実施形態にかかる図９の方法において使用するためのパーティション化された圧縮表現で符号化された行列に対してＳｐＭＴＶを実行するためのルーチンを示すフロー図である。 図１５は、１つの実施形態にかかる図９の方法において使用するためのパーティション化された圧縮表現処理のためのカーネル関数を起動するためのルーチンである。 図１６は、１つの実施形態にかかる図１５のルーチンにおいて使用するための選択されたカーネル関数の起動引数を設定するためのルーチンである。 図１７は、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｔ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチンを示すフロー図である。 図１８は、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｔ１Ｒカーネル関数によってＳｐＭＴＶを実行するためのルーチンを示すフロー図である。 図１９Ａは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｗ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチンを示すフロー図である。 図１９Ｂは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｗ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチンを示すフロー図である。 図２０Ａは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｗ１ＲカーネルによってＳｐＭＴＶを実行するためのルーチンを示すフロー図である。 図２０Ｂは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｗ１ＲカーネルによってＳｐＭＴＶを実行するためのルーチンを示すフロー図である。 図２１Ａは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｂ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチンを示すフロー図である。 図２１Ｂは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｂ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチンを示すフロー図である。 図２２Ａは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｂ１Ｒカーネル関数によってＳｐＭＴＶを実行するためのルーチンを示すフロー図である。 図２２Ｂは、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｂ１Ｒカーネル関数によってＳｐＭＴＶを実行するためのルーチンを示すフロー図である。 図２３は、１つの実施形態にかかる図６及び図９の方法において使用するためのべき乗法を実行するためのルーチンを示すフロー図である。

上述したように、疎行列は、圧縮された疎行（ＣＳＲ）フォーマットで符号化されることができる。Ａをｅ個の非ゼロエントリを有するｍ×ｎの疎行列とする。Ａ_ｖ、Ａ_ｃ及びＡ_ｒを、それぞれ、長さｅ、ｅ及び（ｍ＋１）の３つの１次元配列とする。ＣＳＲフォーマットにおいて、Ａは、＜Ａ_ｖ、Ａ_ｃ、Ａ_ｒ＞として符号化され、
・値配列Ａ_ｖは、行優先順序でＡにおける全ての非ゼロエントリの値を保持し、
・列配列Ａ_ｃは、行優先順序でＡにおける全ての非ゼロエントリの列を保持し、
・行インデックス配列Ａ_ｒは、Ａにおける各行の第１の非ゼロエントリのＡ_ｖにおけるインデックスを保持し、Ａ_ｖ配列におけるエントリの総数であるＡ_ｒ配列における最後のエントリを有する。

例えば、以下の４×４行列

は、＜Ａ_ｖ＝［１，０．５，０．５，０．５，０．５，１］、Ａ_ｃ＝［３，０，３，０，１，２］、Ａ_ｒ＝［０，１，３，５，６］＞としてＣＳＲにおいて符号化されることができる。

以下に記載される本例及び他の例において、行及び列は、それらのｉｄとして機能するインデックスによってインデックス付けされる。インデックスは、ゼロから始まり、行列の左から右へと行くと列が増加し、上から下へと行くと行が増加する。それゆえに、上記行列Ａは、０から３まで行及び列のｉｄを有する。

ｘ及びｙを、それぞれ、サイズｎ及びｍの２つの密ベクトルとする。ＳｐＭＶのタスクは、ｙ＝Ａｘを計算することであり、Ａは疎行列である。図１は、１つの実施形態（従来技術）にかかる、圧縮された疎行フォーマット符号化において符号化された行列に対してＳｐＭＶを実行するための方法１０を示すフロー図である。ｉ＝０からｍ−１のｉｄを有する行である行列の全ての行を処理する反復処理ループ（ステップ１１〜１９）が開始される（ステップ１１）。Ａ_ｒ配列におけるｉ行目のエントリが配置され、行列におけるその行についての第１の非ゼロエントリのＡ_ｖ配列におけるインデックスを識別し、識別されたインデックスは、変数ｊとして設定される（ステップ１２）。Ａ_ｒ配列における次の（ｉ＋１）行目のエントリが配置され、変数ｊ_ｍａｘとして設定される（ステップ１３）。ｉ行目が行列における最後の行でない限り、Ａ_ｒ配列における次のエントリは、（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖ配列におけるインデックスである。ｉ行目が行列における最後の行である場合、Ａ_ｒ配列における次のエントリは、Ａ_ｖ配列におけるエントリの総数である。以下のステップ１６において記載される非ゼロ配列の値の乗算結果をともに加算する関数である総和計算部（以下の擬似コードにおける総和積算器とも称される）は、ゼロにおいて合計値を設定することによって初期化される（ステップ１４）。ｊがｊ_ｍａｘ未満である場合（ステップ１５）、以下の式による計算
ｓｕｍ←ｓｕｍ＋Ａ_ｖ［ｊ］×ｘ［Ａ_ｃ［ｊ］］
が実行される（ステップ１６）。計算において、インデックスｊを有するＡ_ｖ配列における値はｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｃ配列の数である。乗算結果は、ステップ１６の前の反復中に実行された乗算結果の合計に加算される。上記ステップ１４において設定されたように、本方法におけるステップ１６の第１の反復中において合計はゼロである。計算が終了すると、ｊの値に１が加算され、加算結果はｊに設定され（ステップ１７）、その行における次の列のエントリに処理を移動する。本方法は、上述したステップ１５に戻り、ｉ行目における非ゼロ値が処理されるまでステップ１５〜１７を繰り返す。ｊがｊ_ｍａｘ以上である場合（ステップ１５）、ループ１５〜１７における反復中に乗算結果を加算した合計は、密ベクトルｙに格納される（ステップ１８）。反復処理ループは、次の行に移動し（ステップ１９）、全ての行が処理されるまでループ（ステップ１１〜１９）を介して処理を継続し、方法１０は終了する。

図１の方法１０はまた、以下の擬似コードを使用して表すことができる。
ｆｏｒ ｉ＝０からｍ−１ ／＊疎行列Ａのｍ行に対するループ＊／
ｊ←Ａ_ｒ［ｉ］ ／＊ｊ：ｉ行目における第１の非ゼロエントリのＡ_ｖにおけるインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器を初期化＊／
ｗｈｉｌｅ （ｊ＜ｊ_ｍａｘ） ／＊行の終わりに到達したかどうかのテスト／＊
ｓｕｍ←ｓｕｍ＋Ａ_ｖ［ｊ］×ｘ［Ａ_ｃ［ｊ］］ ／＊ｙｉ＝ΣｊＡｉ，ｊ×ｘｊを計算＊／
ｊ←ｊ＋１ ／＊ｉ行目における次の非ゼロの列に移動＊／
ｗｈｉｌｅ文終了
ｙ［ｉ］←ｓｕｍ ／＊ｙに結果を格納＊／
ｆｏｒ文終了

ＳｐＭＶの最も計算集約的な部分は、上記ステップ１６において記載された総和行に起因する。
ｓｕｍ←ｓｕｍ＋Ａ_ｖ［ｊ］×ｘ［Ａ_ｃ［ｊ］］
ここで、それぞれインデックスｊ、ｊ及びＡ_ｃ［ｊ］を有するＡ_ｖ、Ａ_ｃ及びｘの３つの配列がアクセスされる。Ａ_ｖ及びＡ_ｃは、双方とも、通常はＳｐＭＶにおいてＡ_ｒ及びｘのサイズよりもはるかに大きいｅ個の要素を有することに留意されたい。ＣＰＵ及びＧＰＵの双方を含む現代のプロセッサにおけるこのサイズ差の意味は、上記総和行がこれらの配列のサイズに応じて２つ又は３つのキャッシュミスを引き起こす可能性が最も高いということである。３つの配列のうち、ｘ（入力密ベクトルを格納する配列）は、通常最小であり、それゆえに、それは、プロセッサのＬ２キャッシュに適合する最善の機会を有する。もちろん、行列Ａが十分に小さい場合、全てが適合し、自明なケースである。しかしながら、非自明のＳｐＭＶの問題について、Ａ_ｖ又はＡ_ｃのいずれもＬ２に収まることを仮定するべきではない。換言すれば、Ａ_ｖ及びＡ_ｃの双方にアクセスすることは、２つの別個のキャッシュミスをトリガすることができ、ＳｐＭＶの性能に悪影響を与えることがある。キャッシュミスは、ＳｐＭＶの速度を大幅に低減することができ、大抵の場合にメモリ参照の局所性をほとんど示さない。

疎行列はまた、大抵の場合に「転置」のＣＳＲと思われる圧縮された疎列（ＣＳＣ）フォーマットで符号化されることができる。Ａをｅ個の非ゼロエントリを有するｍ×ｎの疎行列とする。Ａ’_ｖ、Ａ’_ｃ及びＡ’_ｒを、それぞれ、長さｅ、ｅ及び（ｎ＋１）の３つの１次元配列とする。ＣＳＣフォーマットにおいて、Ａは、＜Ａ’_ｖ，Ａ’_ｃ，Ａ’_ｒ＞として符号化され、ここで、
・値配列Ａ’_ｖは、列優先順序でＡにおける全ての非ゼロエントリの値を保持し、
・列配列Ａ’_ｒは、列優先順序でＡにおける全ての非ゼロエントリの行を保持し、
・行インデックス配列Ａ’_ｃは、Ａにおける各列の第１の非ゼロエントリのＡ’_ｖにおけるインデックスを保持し、Ａ’_ｖ配列におけるエントリの総数であるＡ’_ｃ配列における最後のエントリを有する。

先に示された同じ４×４行列

は、＜Ａ’_ｖ＝［０．５，０．５，０．５，１，１，０．５］、Ａ’_ｒ＝［１，２，２，３，０，１］、Ａ’_ｃ＝［０，２，３，４，６］＞としてＣＳＣにおいて符号化されることができる。

ｘ’及びｙ’を、それぞれ、サイズｍ及びｎの２つの密ベクトルとする。ＳｐＭＴＶのタスクは、ｙ’＝ＡＴｘ’を計算することであり、Ａは元の転置されない疎行列である。図２は、１つの実施形態（従来技術）にかかる、圧縮された疎列フォーマット符号化において符号化された行列に対してＳｐＭＴＶを実行するための方法２０を示すフロー図である。ｊ＝０からｎ−１のｉｄを有する列である行列の全ての列を処理する反復処理ループ（ステップ２１〜２９）が開始される（ステップ２１）。Ａ’_ｃ配列におけるｊ列目のエントリが配置され、行列におけるその列についての第１の非ゼロエントリのＡ’_ｖ配列におけるインデックスを識別し、識別されたインデックスは、変数ｉとして設定される（ステップ２２）。Ａ’_ｃ配列における次の（ｊ＋１）列目のエントリが配置され、変数ｉ_ｍａｘとして設定される（ステップ２３）。ｊ列目が行列における最後の列でない限り、Ａ’_ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖ配列におけるインデックスである。ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖ配列におけるエントリの総数である。以下のステップ２６において記載される非ゼロ配列の値の乗算結果をともに加算する関数である総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ２４）。ｉがｉ_ｍａｘ未満である場合（ステップ２５）、以下の式による計算
ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］
が実行される（ステップ２６）。計算において、インデックスｉを有するＡ’_ｖ配列における値はｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｒ配列の数である。乗算結果は、ステップ２６の前の反復中に実行された乗算結果の合計に加算される。上記ステップ２４において設定されたように、本方法におけるステップ２６の第１の反復中において合計はゼロである。計算が終了すると、ｉの値に１が加算され、加算結果はｉに設定され（ステップ２７）、その列における次の行のエントリに処理を移動する。方法２０は、上述したステップ２５に戻り、ｊ列目における非ゼロ値が処理されるまでステップ２５〜２７を繰り返す。ｉがｉ_ｍａｘ以上である場合（ステップ２５）、ループ２５〜２７における反復中に乗算結果を加算した合計は、密ベクトルｙ’に格納される（ステップ２８）。反復処理ループは、次の列に移動し（ステップ２９）、全ての列が処理されるまでループ（ステップ２１〜２９）を介して処理を継続し、方法２０は終了する。図２の方法２０はまた、以下の擬似コードを使用して表すことができる。
ｆｏｒ ｊ＝０からｎ−１ ／＊疎行列Ａのｎ列に対するループ＊／
ｉ←Ａ’_ｃ［ｊ］ ／＊ｉ：ｊ列目における第１の非ゼロエントリのＡ’_ｖにおけるインデックス＊／
ｉ_ｍａｘ←Ａ’_ｃ［ｊ＋１］ ／＊ｉ_ｍａｘ：（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器を初期化＊／
ｗｈｉｌｅ （ｉ＜ｉ_ｍａｘ） ／＊列の終わりに到達したかどうかのテスト／＊
ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］ ／＊ｙ’ｊ＝ΣｉＡｉ，ｊ×ｘ’ｉを計算＊／
ｉ←ｉ＋１ ／＊ｊ列目における次の非ゼロの行に移動＊／
ｗｈｉｌｅ文終了
ｙ’［ｊ］←ｓｕｍ ／＊ｙ’に結果を格納＊／
ｆｏｒ文終了

上記総和行ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］は、それぞれインデックスｉ、ｉ及びＡ’_ｒ［ｉ］を有するＡ’_ｖ、Ａ’_ｒ、ｘ’の３つの配列にアクセスする。ＣＳＲの場合と同様に、３つの配列のアクセスは、図２の方法２０のステップ２５〜２７の１回の反復において３つのキャッシュミスを引き起こす可能性がある。

一般に疎行列のキャッシュミスの量及び処理速度は、疎行列の改善された符号化を使用することによって低減することができる。図３は、１つの実施形態にかかる、疎行列の効率的な表現及び処理のためのコンピュータ実装システム３０を示すブロック図である。システム３０は、１つ以上の疎行列３２を記憶するデータベースなどのコンピュータ記憶装置３１を含む。記憶装置は、行列３２の圧縮表現３５を用意する行列エンコーダ３４を実行する１つ以上のサーバ３３に接続されている。各表現は、行列における非ゼロエントリに関する情報を記憶する１つ以上のアレイを含む。以下の説明において、用語「表現」及び「符号化」は、互換的に使用される。表現３５は、記憶装置３１に記憶されることができる。各圧縮表現３５は、行列３２にゼロエントリを明示的に含まない圧縮フォーマットにおける行列３２の符号化である。例えば、そのような表現３５は、上述したＣＳＲ及びＣＳＣフォーマットにおける符号化を含むことができる。行列エンコーダ３４はまた、処理中におけるキャッシュミスの数を低減し、それゆえによりキャッシュフレンドリーである行列のより効率的な表現３５を符号化することができる。

これらのより効率的な符号化の１つは、構造化ＣＳＲ（ＳＣＳＲ）符号化と称することができる。上述したＣＳＲ符号化において、Ａ_ｖ及びＡ_ｃは、Ａ_ｖにおける要素がＡ_ｃにおける要素と同じバイトサイズを有しない場合にそれらの空間的要件が異なるものの、それらにおける同数の要素を有する。同数の要素を正確に有する値（Ａ_ｖ）及び列（Ａ_ｃ）配列は、それらが疎行列における非ゼロエントリの同じセットの異なる態様を記載することから一致しない。その結果、それらは、いかなるＣＳＲ符号化においても同数の要素を有する必要がある。

２つの別個の配列に同じ非ゼロエントリの値及び列を格納する代わりに、構造化ＣＳＲの符号化方式は、ｉ番目の要素が疎行列におけるｉ番目の非ゼロエントリの値のみならず列も記憶するように、対（値、列）を含む単一の統一された配列で値及び列を置く。換言すれば、ＣＳＲ符号化における配列Ａ_ｖ及びＡ_ｃは、以下では値−列配列と称する複合配列Ａ_ｖｃを形成するように互いに結合されている。行インデックス配列Ａ_ｒは、Ａにおける各行の第１の非ゼロエントリのＡ_ｖｃにインデックスを記憶し、Ａ_ｒ配列における最後のエントリが配列Ａ_ｖｃにおける複合エントリの総数である。それゆえに、構造化ＣＳＲ符号化は、２つの配列Ａ＝＜Ａ_ｖｃ，Ａ_ｒ＞を含む。ＳＣＳＲ符号化の名称は、値−列配列Ａ_ｖｃにおける各要素は、浮動小数又は整数などの基本的なタイプである代わりに、構造（すなわち、複合データタイプ）であることを反映する。

上記に示された例としての行列

は、以下のように構造化ＣＳＲにおいて符号化されることができる。Ａ＝＜Ａ_ｖｃ，Ａ_ｒ＞、ここで、Ａ_ｖｃ＝［＜１，３＞，＜０．５，０＞，＜０．５，３＞，＜０．５，０＞，＜０．５，１＞，＜１，２＞］及びＡ_ｒ＝［０，１，３，５，６］。

同様に、エンコーダ３４は、ＣＳＣと比較してＳｐＭＴＶなどの処理中にキャッシュミスの数を低減することができる構造化ＣＳＣ（ＳＣＳＣ）符号化を形成することができる。構造化ＣＳＣにおいて、ＣＳＣの値配列Ａ’_ｖ及び行配列Ａ’_ｒは、値−列配列と称することができる単一の複合配列Ａ’_ｖｒにマージされる。値−列配列Ａ’_ｖｒの各要素は、行列の左端列から始まる列優先順序に基づく非ゼロ行列エントリの値及び列の双方を記憶する＜値，行＞対である。列インデックス配列Ａ’_ｃは、Ａにおける各列の第１の非ゼロエントリのＡ’_ｖｒにインデックスを保持し、Ａ’_ｃ配列における最後のエントリは、Ａ’_ｖｒ配列における複合エントリの総数である。

例えば、先に示された同じ行列

は、以下のように構造化ＣＳＣにおいて符号化されることができる。Ａ＝＜Ａ’_ｖｒ，Ａ’_ｃ＞、ここで、Ａ’_ｖｒ＝［＜０．５，１＞，＜０．５，２＞，＜０．５，２＞，＜１，３＞，＜１，０＞，＜０．５，１＞］及びＡ’_ｃ＝［０，２，３，４，６］。

サーバ３３は、さらに、行列３２の圧縮表現３５に対する計算を行うことによって行列３２を処理する計算モジュール３６を含む。それゆえに、計算モジュール３６は、さらに図７及び図８を参照して後述するように、行列３２の構造化ＣＳＣ及び構造化ＣＳＲ符号化に対してＳｐＭＶ及びＳｐＭＴＶを実行することができる。また、計算モジュール３６は、多くの用途に向けてＳｐＭＶ及びＳｐＭＴＶの結果を適用することができる。例えば、ＳｐＭＶ及びＳｐＭＴＶの結果の１つの適用は、さらに図２３を参照して記載されるようにページランク（登録商標）アルゴリズムなどのランク付けアルゴリズムを実行することができる。ページランクアルゴリズムについてのグラフと疎行列との間の接続が以下に簡単に検討される。グラフは、通常は疎であり、それゆえに大抵の場合にＣＳＲ／ＣＳＣフォーマットで符号化される隣接行列として表現することができることがよく知られている。図４は、一例として、隣接行列を示すグラフ４０である。例えば、行列

は、図４におけるグラフ４０の行列表現としてみることができる。上記グラフは、以下の４つの頂点及び６つの重み付けエッジを有する。＜０，３，１．０＞，＜１，０，０．５＞，＜１，３，０．５＞，＜２，０，０．５＞，＜２，１，０．５＞，＜３，２，１．０＞、ここで、フォーム＜ｕ，ｖ，ｗ＞のエッジは、ｕがエッジの元頂点であり且つｖがエッジの先頂点であるように、重みｗを有する頂点ｕからｖまでのリンクを示す。これらの６つのエッジが＜行，列，値＞トリプルの形態で疎行列の非ゼロエントリとして表示される場合、それらが正確にＡにおける非ゼロエントリであることを確認することができる。それゆえに、重み付きグラフと疎行列との間に１対１の対応がある。グラフ４０は、エッジの重みがその元頂点から出るリンクの数の逆数であるという興味深い性質を有する。そのようなグラフは、頂点間の遷移確率をモデル化し、それは、さらに図２３を参照しながら後述するようなページランク（登録商標）アルゴリズムを計算するために使用することができる。他のアルゴリズムもまた、さらに後述するように計算モジュール３６によって実行することができる。

構造化ＣＳＣ及び構造化ＣＳＲは、任意のタイプのハードウェアを使用して行列３２の処理を高速化することを可能とする。図３を参照すると、エンコーダ３４はまた、処理のためにＧＰＵを使用するサーバ３３において行列処理を高速化するために特に適している圧縮表現３５の他のタイプを形成することが可能である。ＳｐＭＶ又はＳｐＭＴＶなどのＧＰＵベースの疎行列処理について、サーバ３３は、ＳｐＭＶ処理カーネル、後述するものなどのＧＰＵの処理アルゴリズムを実装するサーバにおける１つ以上のＧＰＵにおいて実行される処理スレッドによって実現される機能を必要とする。しかしながら、上述したように、実験は、単一のカーネルが大抵の場合には疎行列行を有することができる非ゼロエントリの数における大きな変動に応答しないことから、単一のカーネルが最適から遠いことを示している。

サーバ３３は、スレッドの並列グループ内で実行することが可能である。例えば、サーバ３３は、ワープと称される単位にグループ化された複数のＧＰＵスレッドを含むことができる。例えば、カリフォルニア州サンタクララのＮｖｉｄｉａ（登録商標）社製のＧＰＵの並列スレッドは、同じワープにおける全てのスレッドが単一のストリーミングマルチプロセッサ（ＳＭ）を共有し且つ単一命令複数スレッド（ＳＩＭＴ）実行モデルを想定するように、ワープと称される３２の単位にグループ化される。スレッドがカリフォルニア州サニーベールのアドバンスト・マイクロ・デバイス（登録商標）社製のＧＰＵにおいて実行される場合、ワープと同等の概念は、６４個のＧＰＵコアのグループにおいて実行される６４スレッドを現在含むウェーブフロントと称される。その結果、同じワープ（又はウェーブフロント）内の異なるスレッドは、それらが常に並行して異なるデータを処理するのが可能とされるものの、同時に異なる命令を実行することができない。本特許出願の目的のために、用語「ワープ」及び「ウェーブフロント」は、単一のストリーミングマルチプロセッサを共有し且つ同じ命令を実行する３２又は６４の複数の処理スレッドを並列に実行するＧＰＵコアの集合を指すように以下において互換的に使用される。さらなる実施形態において、他のスレッド数は、ワープにおいて可能である。さらに、スレッドの大きなグループ化は、ワープのグループ化によって可能である。そのような各グループは、スレッドブロックと称され、スレッドの各ブロックは、（実施形態に応じて３２又は６４以上）ワープよりも多くのスレッドを含む。各ブロック内のスレッド数は、ブロックサイズ（又は以下の説明においてはＢＬＯＣＫＳＩＺＥ）と称される一方で、ワープ内のスレッド数は、ワープサイズ（又は以下の説明においてはＷＡＲＰＳＩＺＥ）と称される。エンコーダ３４は、サーバ３３が異なる数の非ゼロエントリを有する行列３２の部分について異なるカーネルを実行するのを可能とする行列３２の符号化を形成することによってサーバ３３のハードウェアの限界を克服するのを可能とする。エンコーダ３４は、複数のパーティションを含むパーティション化された圧縮表現で行列３２のパーティション化された圧縮表現を作成する。各パーティションは、所定範囲内の複数の非ゼロエントリを有する行又は列などの行列３２の部分のグループのリストである。例えば、１つのパーティションは、１〜３１の非ゼロエントリを有する行列３２の部分を有することができ、２番目のパーティションは、３２から１０２４の非ゼロエントリを有する行列３２の部分を有することができ、３番目のパーティションは、１０２４以上の非ゼロエントリを有する行列のエントリ部分を有することができる。非ゼロエントリの他の範囲及びパーティションの他の数も可能である。エンコーダ３４は、以下のようにパーティション化された圧縮表現を定義する。以下の説明は、行であるとしてパーティション化されている行列３２の部分に言及しているが、行列の列は、同じ方法を準用してパーティション化されることができる。

ｋをゼロから始まるパーティションのインデックスとし、少数の非ゼロエントリを有する行のパーティションが多い非ゼロエントリを有するパーティションよりも低いインデックスを有する。それゆえに、３２から１０２４の非ゼロエントリを有する行のパーティションは、１０２４超の非ゼロエントリを有するパーティションよりも低いインデックスｋを有するであろう。ｐをパーティションの総数とし、Ａｓを、以下のように、各パーティションにおいて許容される行の非ゼロエントリの最小及び最大数を指定する行パーティション区切り配列と称される（ｐ＋１）要素の整数配列とする。
・Ａｓ［０］＝１
・ｋ＝０，１，２，・・・，ｐ−１について、Ａｓ［ｋ］＜Ａｓ［ｋ＋１］
・Ａｓ［ｐ］＝∞

１つの実施形態において、Ａｓ［０］及びＡｓ［ｐ］の値は、それらが行列３２の行が有することができる非ゼロエントリの最大及び最小数であることから、記憶装置に記憶されない。さらなる実施形態において、値は記憶装置３１に記憶される。

パーティションは、順序を維持する。ｅｉ〜を疎行列の行ｉにおける非ゼロエントリの数とする。Ｒｋとして示されるｋ行目の行ベースのパーティションにおける行のセットは、Ｒｋ＝｛ｉ｜Ａｓ［ｋ］≦ｅｉ〜＜Ａｓ［ｋ＋１］｝のように記載することができる。元の行列３２の内容を変更しないことが大抵の場合には望ましいことから、以下ではマッピング配列とも称される１次元の順序維持置換配列（Ａ_ｍ又はＡ’_ｍ）は、行のランク又は列のランクと称されるパーティションベースの行ｉｄからそれぞれ元の行又は列ｉｄへのマッピングを格納するために追加される。ランクのそれぞれは、行列の部分のマッピング配列におけるインデックスである。例えば、行が３のｉｄ及び０のランクを有する場合、値３は、マッピング配列にリスト化された第１の値である。通常の置換配列とは異なり、順序維持置換配列は、以下のように同じパーティションに割り当てられた行の相対的な順序を維持するために必要とされる。ｒｉ及びｒｊを置換配列における行ｉ及びｊのランクとする。ｐｉ及びｐｊをこれら２行ｉ及びｊが属するパーティションとする。ランクｒｉ及びｒｊは、以下の制約が満たされた場合に限り、順序を維持している。
・∀ｐｉ＜ｐｊ，ｒｉ＜ｒｊ−パーティションｊのインデックスがパーティションｉのインデックスよりも大きい場合、パーティションｊにおける行のランクは、パーティションｉにおける行のランクよりも大きくなければならない。
・∀ｐｉ＞ｐｊ，ｒｉ＞ｒｊ−パーティションｉのインデックスがパーティションｊのインデックスよりも大きい場合、パーティションｉにおける行のランクは、パーティションｊにおける行のランクよりも大きくなければならない。
・∀ｐｉ＝ｐｊ∧ｉ≠ｊ，ｒｉ＜ｒｊ⇔ｉ＜ｊ（又は同等に、ｒｉ＞ｒｊ⇔ｉ＞ｊ）

上記最後の制約は、小さいｉｄを有する行が同じパーティションにおいてより大きなｉｄを有する行の前にアクセスされるように、元の行列３２における同じパーティションの行の相対的な順序を尊重するように設計されている。制約は、構造化及び非構造化ＣＳＲ符号化の双方においてＡ_ｒにアクセスする際のランダム性を低減させることから、制約は、ＳｐＭＶなどのＳＬＡアルゴリズムのキャッシュ性能を向上させる。有益な副作用として、順序維持の制約を強制することもまた、ＣＳＲ行列を完全にパーティション化するために行インデックス配列Ａ_ｒの単一のキャッシュフレンドリーな線形スキャンを必要とするのみであることから、パーティション化アルゴリズムの時間的複雑性を低減する。一方、ソートベースのアルゴリズムは、かなりの数のキャッシュミスをもたらすことがあるＯ（ｍｌｏｇｍ）比較を平均的に負担する。ここで、ｍは行数である。実験は、順序維持パーティションがはるかに高速に計算できることのみならず、それらはまた、通常順序維持ではないソートベースのパーティションよりも大幅に高速なＳｐＭＶをもたらすことを示す。

行列３２の順序維持のパーティション化された圧縮表現は、ＣＳＣ、ＣＳＲ、構造化ＣＳＣ及び構造化ＣＳＲなどの行列３２の様々な既存の圧縮表現に基づいて図９から始めて以下に記載されるようなエンコーダ３４によって作成することができる。異なる圧縮表現に基づくパーティション化された圧縮表現は互いに異なる。例えば、疎行列Ａは、構造化ＣＳＲについてＡ＝＜Ａ_ｖｃ，Ａ_ｒ，Ａ_ｍ，Ａｏ，Ａｐ，Ａｓ＞として表現されることができる。ここで、Ａ_ｖｃ及びＡ_ｒは、前と同じであり、Ａ_ｍは、行パーティションマッピング配列であり、Ａｏは、パーティションのそれぞれの第１の部分のランク及びマッピング配列におけるエントリの総数を含む行パーティションオフセット配列であり、Ａｐは、行ベースのパーティションの数であり、Ａｓは、（必要に応じてパーティション化が終了した場合）行パーティション区切り配列である。行ベースのパーティションを有する構造化ＣＳＲは、この時点からパーティション化された構造化ＣＳＲ（ＰＳＣＳＲ）と称される。非構造化ＣＳＲ行列は、同様に図９を参照して以下に記載されるルーチンを使用してパーティション化されることができるため、得られた符号化は、疎行列Ａが＜Ａ_ｖ，Ａ_ｃ，Ａ_ｒ，Ａ_ｍ，Ａｏ，Ａｐ，Ａｓ＞として符号化されるパーティション化ＣＳＲ（ＰＣＳＲ）と称される。

同様に、疎行列Ａは、構造化ＣＳＣについてＡ＝＜Ａ’_ｖｒ，Ａ’ｃ，Ａ’_ｍ，Ａ’_ｏ，Ａ’_ｐ，Ａ’_ｓ＞として表現することができる。ここで、Ａ’_ｖｒ及びＡ’_ｃは、前と同じであり、Ａ’_ｍは、行パーティションマッピング配列であり、Ａ’_ｏは、行パーティションオフセット配列であり、Ａ’_ｐは、列ベースのパーティション数であり、Ａ’_ｓは、（必要に応じてパーティションが終了した場合）列パーティション区切り配列である。列ベースのパーティションを有する構造化ＣＳＣは、この時点からパーティション化された構造化ＣＳＣ（ＰＳＣＳＣ）と称される。非構造化ＣＳＣ行列は、同様に同じアルゴリズムを使用してパーティション化されることができるため、得られた符号化は、疎行列Ａが＜Ａ’_ｖ，Ａ’_ｒ，Ａ’_ｃ，Ａ’_ｍ，Ａ’_ｏ，Ａ’_ｐ，Ａ’_ｓ＞として符号化されるパーティション化ＣＳＣ（ＰＣＳＣ）と称される。

例えば、同一の例として先に示された行列

を考える。
１．第１のパーティションが単一の非ゼロエントリを有する行のみを含み、
２．第２のパーティションが複数の非ゼロエントリを有する行を含む
ように、Ａについて２つのパーティションを作成したい場合、Ａ_ｐ＝Ａ’_ｐ ＝２及びＡ_ｓ＝Ａ’_ｓ＝［１，２，∞］である。他の配列は、以下のとおりである。
・ＰＣＳＲ及びＰＳＣＳＲについて、Ａ_ｍ＝［０，３，１，２］及びＡ_ｏ＝［０，２，４］
・ＰＣＳＣ及びＰＳＣＳＣについて、Ａ’_ｍ＝［１，２，０，３］及びＡ’_ｏ＝［０，２，４］

システム３０のハードウェアリソースが表１に関連して以下に説明される式を使用して決定することができることを考えると、パーティション化された圧縮表現を作成するかどうかは実用的である。パーティション化された圧縮符号化表現は、サーバ３３が１つ以上のＧＰＵ又は１つ以上のＣＰＵのみを含むかどうかにかかわらず、図９を参照して始めて以下に記載されるものなどの計算モジュール３６によって処理することができる。しかしながら、所定数の非ゼロエントリを有する行列３２の部分をリスト化するパーティションの作成は、サーバ３３がＧＰＵを含む場合、計算モジュール３６がそれらの部分を処理するときに最も適切な処理カーネルを適用するのを可能とする。多くの処理カーネルが可能であるが、計算モジュールは、少なくとも以下の３つのカーネルを適用することができる。

１．行又は列を処理するために単一の処理スレッドを割り当てる１スレッド１行（１Ｔ１Ｒ）カーネル
２．行又は列を処理するためにワープと称されるユニット内の全てのスレッドを割り当てる１ワープ１行（１Ｗ１Ｒ）カーネル
それ自体では、いずれのカーネルも、行列３２の全ての部分について理想的ではない。上述したように、総ランタイムは、最大数の非ゼロエントリを有する行（又は列）に対応する最も遅いスレッドに依存することから、（以下ではｆ１Ｔ１Ｒカーネルとも称されることができる）１Ｔ１Ｒカーネルは負荷不均衡を被る。行（又は列）における非ゼロエントリ数が３２（又は６４）未満のワープ（又はウェーブフロント）サイズである場合、（以下ではｆ１Ｗ１Ｒカーネルとも称されることができる）１Ｗ１Ｒカーネルは、ハードウェアリソースの浪費を被る。単一の疎行列は、少数の行及び多くの非ゼロエントリを有する行の双方を有することができることに留意されたい。それゆえに、いずれのタイプのカーネルにコミットすることも、めったに問題を解決しない。さらにまた、実験は、同じ行に対して動作するように単一のワープにおけるものよりも多くのスレッドを有する利益があることを示す。これは、第３のタイプを追加する必要性を引き起こす。
３．行又は列を処理するためにスレッドのブロック（３２又は６４より多い）を割り当てる１ブロック１行（１Ｂ１Ｒ）カーネル
１Ｗ１Ｒカーネルと同様に、（以下ではｆ１Ｂ１Ｒカーネルとも称されることができる）１Ｂ１Ｒカーネルはまた、スレッドの数が行列のその部分における非ゼロエントリの数よりも大きい場合にハードウェアリソースの浪費を被ることがある。しかしながら、行列の特定のパーティションと各カーネルを相関させることにより、計算モジュールは、さらに図１５を参照して以下に記載されるように、特定の行又は列を処理するための最良のカーネルを使用することができる。

カーネル関数のさらに他の種類が可能である。

ｆ１Ｔ１Ｒ、ｆ１Ｗ１Ｒ及びｆ１Ｂ１Ｒ又は効率的に他のカーネルを混合するために、ＳｐＭＶアルゴリズムは、各カーネルが得意である行を迅速に区別する必要がある。高速な行分類について、疎行列パーティション化アルゴリズムは、図９を参照して始めて以下に記載され、同じカーネルによって最良に処理される行列の行をグループ化する。より正確には、さらに図９を参照して始めて以下に記載されるように、パーティション化アルゴリズムは、入力として、各パーティションにおける行についての非ゼロエントリの最小及び最大数をとり、各行又は各列が属するパーティションを決定するために使用することができるマッピングを生成する。さらなる実施形態において、カーネル関数は、パーティションにおける非ゼロエントリ数以外の要因に基づいてパーティションに割り当てることができる。

１つ以上のサーバ３３は、インターネット又はセルラーネットワークなどのローカルネットワーク又はインターネットワークとすることができるネットワーク３７に接続され、ネットワークを介して少なくとも１つのユーザ装置３８と通信することができる。ユーザ装置３８は、デスクトップコンピュータとして示されているが、ユーザ装置３８はまた、ラップトップコンピュータ、スマートフォン、メディアプレーヤ及びタブレットを含むことができる。ユーザ装置３８のさらに他の種類も可能である。ユーザ装置３８は、ネットワーク１５を介してサーバ３３と通信し、計算実行からコマンドを受信し、ユーザ装置に計算結果を出力し返すことができる。

ＳｐＭＶ及びＳｐＭＴＶは、（ＰＳ）ＣＳＲ及び（ＰＳ）ＣＳＣ符号化の上に構築することができる多くのＳＬＡオペレーションの２つだけである。一般に、疎行列を扱う任意のアルゴリズムは、特に、アルゴリズムが行列内の非ゼロエントリの値及び位置の双方に同時にアクセスする必要がある場合、パーティション化された（Ｐ）及び／又は構造化された（Ｓ）ＣＳＲ／ＣＳＣ符号化から利益を得ることができる。疎行列−行列乗算（ＳｐＭＭ）などの他のＳＬＡオペレーションは、同様のアクセスパターンを有し、それゆえに、ここで紹介したのと同じ符号化方式から利益を得ることができる。

しかしながら、従来のＣＳＲ／ＣＳＣ符号化が構造化することなしでさえも良好に作動することができ且つシステム３０が従来のＣＳＲ／ＣＳＣ符号化を使用可能であるタスクがある。例えば、以下のように定義される疎行列Ａのフロベニウスノルム演算

は、非ゼロエントリの位置にアクセスするためのアルゴリズムを必要とせず、それゆえに、同じ配列におけるそれらの値及び位置の混合は常に役立たないことがある。他の例は、非ゼロエントリの位置ではなくそれらの値にアクセスすることのみを必要とする疎行列が対角であるか否かをテストすることである。一般に、非ゼロ行列エントリの値又は位置（双方ではない）のいずれかにアクセスする必要があるのみである任意のアルゴリズムは、構造化ＣＳＲ／ＣＳＣから利益を受けないことがある。上記の双方の例は、データプリフェッチが有効な場合、いくつかのキャッシュミスを通常もたらすか又はもたらさない優れたメモリ参照局所性を既に有することに留意されたい。換言すれば、最悪の場合でも、それらの構造化されていない相手方を失うことは構造化ＣＳＲ／ＣＳＣについて多くはなく、一般に構造化ＣＳＲ／ＣＳＣ符号化は、これらのシナリオにおいては非構造化ＣＳＲ／ＣＳＣ符号化を実行する。それにもかかわらず、構造化ＣＳＲ／ＣＳＣは、一般に、ＣＰＵにおいて疎行列を符号化するための最も効率的且つ堅牢なデータ構造である可能性が高く且つ大部分の状況において最も有用であることから、従来のＣＳＲ／ＣＳＣを使用するために（既存のコードを除く）理由はほとんどない。図５は、図３のシステム３０の様々なハードウェアセットアップのための様々な符号化の強みをまとめたグラフ５０である。グラフ５０は、大部分のシナリオにおいて、上記符号化方式が従来の非構造化ＣＳＲ／ＣＳＣ符号化よりも優れていることを図示している。グラフ５０は、（１）ｘ軸としてハードウェアの柔軟性及び（２）ｙ軸としてメモリアクセスパターンの２つの直交次元を含む。ｙ軸について、グラフ５０は、（下方において「値

位置」としてラベル付けされた（シンボル

排他的ＯＲ関係を意味する））非ゼロエントリの値又は位置のいずれか（双方ではない）にアクセスする必要があるのみであるパターンに対して、（上方において「値＋位置」としてラベル付けされた）非ゼロ行列エントリの値及び位置の双方に同時にアクセスする必要があるメモリアクセスパターンを区別する。

グラフ５０は、４つの象限のうち３つが本特許出願において導入された疎行列符号化を好み、非ゼロエントリについてメモリアクセスパターンが値のみ又は位置のみのいずれかである場合及びハードウェアがより多くのＣＰＵのようである場合には従来の非構造化ＣＳＲ／ＣＳＣフォーマットを好む唯一のシナリオがあることを示している。先に説明したように、従来のＣＳＲ／ＣＳＣフォーマットによって達成可能な性能向上は、全ての場合に非常に限られている。一方、他の３つの象限、特に右上隅におけるその性能低下は大幅であり得て、全ての状況について上述した符号化の使用を正当化する。

図３を参照すると、１つ以上のサーバ３３は、他の要素が可能であるものの、１つ以上のＣＰＵ及びＧＰＵ及びＳＭ、メモリ、入力／出力ポート、ネットワークインターフェース並びに不揮発性記憶装置などのプログラム可能な計算装置において従来みられる要素を含むことができる。サーバ３３は、それぞれ、本願明細書に開示される実施形態を実施するための１つ以上のモジュールを含むことができる。モジュールは、従来のプログラミング言語におけるソースコードとして書かれたコンピュータプログラム又はプロシージャとして実装することができ、それは、オブジェクト又はバイトコードとして中央処理装置による実行のために提示される。あるいは、モジュールはまた、集積回路として又は読み出し専用メモリ要素に焼き付けられたものとしてハードウェアに実装されることができ、サーバ３３のそれぞれは、専用コンピュータとして機能することができる。例えば、モジュールがハードウェアとして実装される場合、その特定のハードウェアは、その目的のために使用することができないハードウェアなしでパーティション化及び他のコンピュータを実行するために特化されている。ソースコード及びオブジェクト及びバイトコードの様々な実装は、フロッピー（登録商標）ディスク、ハードドライブ、ディジタルビデオディスク（ＤＶＤ）、ランダムアクセスメモリ（ＲＡＭ）、読み出し専用メモリ（ＲＯＭ）及び同様の記憶媒体などのコンピュータ読み取り可能な記憶媒体に保持することができる。他の種類のモジュール及びモジュール機能のみならず、他の物理的なハードウェア要素も可能である。

上述したように、構造化ＣＳＲ及び構造化ＣＳＣは、ＳｐＭＶなどの行列処理中におけるキャッシュミス数を低減するのを可能とする。図６は、１つの実施形態にかかる構造化行列表現を取得するためのコンピュータ実装方法６０を示すフロー図である。他の実装が可能であるものの、方法６０は、図３のシステム３０に実装することができる。最初に、疎行列がアクセスされ、行列の部分、行及び列における非ゼロエントリについての情報が取得される（ステップ６１）。行列の構造化圧縮表現、行列の構造化ＣＳＣ又は構造化ＣＳＲのいずれかは、上述したように、複合配列（Ａ_ｖｃ又はＡ’_ｖｒ）及びインデックス配列（Ａ_ｒ又はＡ’_ｃ）に非ゼロエントリに関する情報を符号化することによって作成される（ステップ６２）。他の種類の処理も可能であるが、図７を参照して後述するＳｐＭＶ、又は、図８を参照して後述するＳｐＭＴＶを実行することなどによって構造化圧縮表現が処理される（ステップ６３）。処理結果は、処理結果の他の種類のアプリケーションも可能であるが、ページランク（登録商標）又は他のランク付けアルゴリズムを実行することなどによって適用される（ステップ６４）。方法６０は終了する。

行列の構造化圧縮表現に対して実行されるＳｐＭＶは、キャッシュミスの確率を低減する。図７は、図６の方法６０において使用するために構造化ＣＳＲ符号化を使用して符号化された行列に対してＳｐＭＶを実行するためのルーチン７０である。ａｖｃ∈Ａｖｃを値−列配列の要素とする。ａｖｃが構造であることから、以下のように配列メンバにアクセスするためにＣ＋＋及びＪａｖａ（登録商標）などの一般的なプログラミング言語におけるようなドット演算子が使用可能である。ａｖｃ・ｖは、ａｖｃの「値」フィールドを返し、ａｖｃ・ｃは、ａｖｃの「列」フィールドを返し、それらは後述するステップ７６における計算に使用される。値及び列フィールドにアクセスする他の方法も可能である。ルーチン７０及び後述する後続のＳｐＭＶルーチンを説明するために、上記例において使用された同じｍ×ｎ行列及びサイズｎ及びｍの２つの密ベクトルｘ及びｙがそれぞれ使用される。

行列の全ての行であるｉ＝０からｍ−１のｉｄを有する行を処理する反復処理ループ（ステップ７１〜７９）が開始される（ステップ７１）。Ａ_ｒ配列におけるｉ行目についてのエントリが配置され、行列のその行についての第１の非ゼロエントリのＡ_ｖｃ配列におけるインデックスを識別し、識別されたインデックスは、変数ｊとして設定される（ステップ７２）。Ａ_ｒ配列における次の（ｉ＋１）エントリが配置され、変数ｊ_ｍａｘとして設定される（ステップ７３）。ｉ行目が行列における最後の行でない限り、Ａ_ｒ配列における次のエントリは、（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃ配列におけるインデックスであり、ｉ行目が行列における最後の行である場合、Ａ_ｒ配列における次のエントリは、Ａ_ｖｃ配列におけるエントリの総数である。ステップ７６において後述する非ゼロ配列の値の乗算結果を加算する関数である総和計算部は、ゼロにおける合計の値を設定することによって初期化される（ステップ７４）。ｊがｊ_ｍａｘ未満である場合（ステップ７５）、以下の式
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］
にかかる計算が行われる（ステップ７６）。計算において、インデックスｊを有するＡ_ｖｃ配列の要素に格納された値は、そのインデックスがｊ番目のインデックスを有するＡ_ｖｃ要素に格納された列のインデックスであるｘの要素によって乗算される。乗算結果は、ステップ７６の前の反復中に実行される乗算結果の合計に加算され、ルーチンにおけるステップ７６の最初の反復中において、合計は、上記ステップ７４において設定されたようにゼロである。計算が終了すると、ｊの値に１が加算され、加算結果は、ｊとして設定され（ステップ７７）、その行における次の列のエントリに処理を移動する。方法は、上述したステップ７５に戻り、ｉ行目における非ゼロ値が処理されるまで、ステップ７５〜７７を繰り返す。ｊがｊ_ｍａｘ以上である場合（ステップ７５）、ループ７５〜７７における反復中の乗算結果を加算した合計は、密ベクトルｙに格納される（ステップ７８）。反復処理ループは、次の行に移動し（ステップ７９）、全ての行がループ（ステップ７１〜７９）を介して処理されるまで、ループを介した行の処理（ステップ７１〜７９）が継続し、その後にルーチン７０は終了する。ルーチン７０はまた、図１を参照して上記に示された擬似コードとの差異を示すテキストボックスにより、以下の擬似コードを使用して表すことができる。
ｆｏｒ ｉ＝０からｍ―１ ／＊疎行列Ａのｍ行にわたるループ＊／
ｊ←Ａ_ｒ［ｉ］ ／＊ｊ：ｉ行目における第１の非ゼロエントリのＡ_ｖｃにおけるインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器を初期化＊／
ｗｈｉｌｅ （ｊ＜ｊ_ｍａｘ） ／＊行の終わりに到達したかどうかをテスト＊／

ｊ←ｊ＋１ ／＊ｉ行目における次の非ゼロ列まで移動＊／
ｗｈｉｌｅ文終了
ｙ［ｉ］←ｓｕｍ ／＊ｙに結果を格納＊／
ｆｏｒ文終了

構造化ＣＳＲに対して実行されるＳｐＭＶと同様に、構造化ＣＳＣに対してＳｐＭＴＶを実行することは、処理が非構造化ＣＳＣに対して実行される場合よりもキャッシュミスが発生する可能性を低減する。図８は、図６の方法６０において使用するために構造化ＣＳＣ符号化を使用して符号化された行列に対してＳｐＭＴＶを実行するためのルーチン８０である。ａ’_ｖｒ∈Ａ’_ｖｒを値−行配列の要素とする。ａ’_ｖｒが構造であることから、以下のように配列メンバにアクセスするためにＣ＋＋及びＪａｖａ（登録商標）などの一般的なプログラミング言語におけるようなドット演算子が使用可能である。ａ’_ｖｒ・ｖは、ａ’_ｖｒの「値」フィールドを返し、ａ’_ｖｒ・ｒは、ａ’_ｖｒの「行」フィールドを返し、それらは後述するステップ８６における計算に使用される。値及び行フィールドにアクセスする他の方法も可能である。ルーチン８０及び以下の他のＳｐＭＴＶルーチンを説明するために、上記例において使用された同じｍ×ｎ行列及びサイズｍ及びｎの２つの密ベクトルｘ’及びｙ’がそれぞれ使用される。

行列の全ての列であるｊ＝０からｎ−１のｉｄを有する列を処理する反復処理ループ（ステップ８１〜８９）が開始される（ステップ８１）。Ａ’_ｃ配列におけるｊ列目についてのエントリが配置され、行列のその列についての第１の非ゼロエントリのＡ’_ｖｒ配列におけるインデックスを識別し、識別されたインデックスは、変数ｉとして設定される（ステップ８２）。Ａ’_ｃ配列における次の（ｊ＋１）エントリが配置され、変数ｉ_ｍａｘとして設定される（ステップ８３）。ｊ列目が行列における最後の列でない限り、Ａ’_ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒ配列におけるインデックスであり、ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖｒ配列におけるエントリの総数である。ステップ８６において後述する非ゼロ配列の値の乗算結果を加算する関数である総和計算部は、ゼロにおける合計の値を設定することによって初期化される（ステップ８４）。ｉがｉ_ｍａｘ未満である場合（ステップ８５）、以下の式
ｓｕｍ←ｓｕｍ＋Ａ’_ｖｒ［ｉ］．ｖ×ｘ’［Ａ’_ｖｒ［ｉ］．ｒ］
にかかる計算が行われる（ステップ８６）。計算において、インデックスｉを有するＡ’_ｖｒ配列の要素に格納された値は、そのインデックスがｉ番目のインデックスを有するＡ’_ｖｒ要素に格納された行のインデックスであるｘ’の要素によって乗算される。乗算結果は、ステップ８６の前の反復中に実行される乗算結果の合計に加算され、本方法におけるステップ８６の最初の反復中において、合計は、上記ステップ８４において設定されたようにゼロである。計算が終了すると、ｉの値に１が加算され、加算結果は、ｉとして設定され（ステップ８７）、その列における次の行のエントリに処理を移動する。ルーチン８０は、上述したステップ８５に戻り、ｊ列目における非ゼロ値が処理されるまで、ステップ８５〜８７を繰り返す。ｉがｉ_ｍａｘ以上である場合（ステップ８５）、ループ８５〜８７における反復中の乗算結果を加算した合計は、密ベクトルｙ’に格納される（ステップ８８）。反復処理ループは、次の列に移動し（ステップ８９）、全ての列がループ（ステップ８１〜８９）を介して処理されるまで、ループを介した列の処理（ステップ８１〜８９）が継続し、その後にルーチン８０は終了する。図８のルーチン８０はまた、図２を参照して上記に示された擬似コードとの差異を示すテキストボックスにより、以下の擬似コードを使用して表すことができる。
ｆｏｒ ｊ＝０からｎ―１ ／＊疎行列Ａのｎ列にわたるループ＊／
ｉ←Ａ’_ｃ［ｊ］ ／＊ｉ：ｊ列目における第１の非ゼロエントリのＡ’_ｖｒにおけるインデックス＊／
ｉ_ｍａｘ←Ａ’_ｃ［ｊ＋１］ ／＊ｉ_ｍａｘ：（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器を初期化＊／
ｗｈｉｌｅ （ｉ＜ｉ_ｍａｘ） ／＊列の終わりに到達したかどうかをテスト＊／

ｉ←ｉ＋１ ／＊ｊ列目における次の非ゼロ行まで移動＊／
ｗｈｉｌｅ文終了
ｙ’［ｊ］←ｓｕｍ ／＊ｙ’に結果を格納＊／
ｆｏｒ文終了

構造化された符号化は、ＧＰＵ及び排他的ＣＰＵベースのハードウェアの双方について有益であるが、処理速度のさらなる向上は、ＧＰＵベースのハードウェアについての圧縮符号化の順序維持パーティションを行うことによって得ることができる。ハードウェアシステムのハードウェアリソースが表１に関連して以下に説明する式を使用して決定することができることを考えると、パーティション化された圧縮表現を作成するかどうかは実用的である。図９は、１つの実施形態にかかる効率的な疎行列パーティション及び処理のためのコンピュータ実装方法９０を示すフロー図である。他の実装が可能であるものの、方法９０は、図３を参照して上述したシステムを使用して実装することができる。疎行列の圧縮表現が取得される（ステップ９１）。そのような圧縮表現は、ＣＳＣ符号化、ＣＳＲ符号化、構造化ＣＳＣ符号化又は構造化ＣＳＲ符号化とすることができる。表現は、サーバ３３によって記憶装置３２からアクセスされるか又は他のソースから取得されることができる。サーバ３３は、圧縮表現で表される行列の部分についてのパーティションを定義する（ステップ９２）。具体的には、サーバ３３は、作成されることになるパーティションの数（Ａｐ及びＡ’_ｐ）及び図３を参照して上述したＡｓ配列における値を定義し、パーティションにおける行列の部分が有することができる非ゼロエントリの範囲を指定する。パーティションが定義されると、さらに図１０を参照して後述するように圧縮表現の順序維持パーティションが実行される（ステップ９３）。作成されたパーティションはマージされ、さらに図１０を参照して記載されたように上述したマッピング配列（Ａ_ｍ又はＡ’_ｍ）を形成し、マッピング配列に基づいてオフセット配列（Ａｏ及びＡ’_ｏ）を定義し、それゆえに、行列の順序維持パーティション圧縮表現を完成する（ステップ９４）。そして、行列のパーティション化された圧縮表現は、図１７〜図２３を参照して後述するように、ページランク（登録商標）アルゴリズムなどの順次適用されるそのような処理の結果を用いたＳｐＭＶ及びＳｐＭＴＶなどの処理のために使用することができ（ステップ９５）、方法９０を終了する。

パーティション化は、それらの順序を維持しながら同様の数の非ゼロエントリを有する行列の部分をグループ化するのを可能とする。図１０は、１つの実施形態にかかる図９の方法９０において使用するための圧縮表現で表される行列の行の順序維持パーティションを実行するためのルーチン１００を示すフロー図である。圧縮表現は、ＣＳＲ及び構造化ＣＳＲの双方とすることができる。行列の全ての行であるｉ＝０からｍ−１のｉｄを有する行を処理する反復処理ループ（ステップ１０１〜１０７）が開始される（ステップ１０１）。Ａ_ｒ配列におけるｉ行目についてのエントリが配置され、初期の圧縮表現がＣＳＲ又は構造化ＣＳＲであったかどうかに応じて、行列におけるその行についての第１の非ゼロエントリのＡｖ又はＡ_ｖｃ配列におけるインデックスを識別し、識別されたインデックスが変数ｊとして設定される（ステップ１０２）。Ａ_ｒ配列における次の（ｉ＋１）エントリが配置され、変数ｊ_ｍａｘとして設定される（ステップ１０３）。ｉ行目が行列における最後の行でない限り、Ａ_ｒ配列における次のエントリは、（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃ（又はＡ_ｖ）配列におけるインデックスであり、ｉ行目が行列における最後の行である場合、Ａ_ｒ配列における次のエントリは、Ａ_ｖｃ（又はＡｖ）配列におけるエントリの総数である。ｉ行目における非ゼロエントリの数は、ｊ_ｍａｘからｊを差し引くことによって決定され（ステップ１０４）、非ゼロエントリの数は、ｅｉ〜として示される。ｅｉ〜がゼロよりも大きい場合（ステップ１０５）、ｅｉ〜の値に基づいて及びパーティションにおいて許容される非ゼロエントリの最大及び最小数に基づいて定義されたパーティションのいずれかに行ｉが割り当てられ、パーティションｋにおける行のリストの最後に行ｉの行ｉｄが追加される（ステップ１０６）。それゆえに、行ｉが属するインデックスｋを有するパーティションは、Ａｓ［ｋ］≦ｅｉ〜＜Ａｓ［ｋ＋１］のようにみられる。１つの実施形態において、パーティションｋは、ｂと設定されたｅｉ〜を有し、ｂ以上である昇順配列Ａｓの第１の要素のインデックスを返すｌｏｗｅｒ＿ｂｏｕｎｄ（Ａｓ，ｂ）として示される関数を使用して求めることができる。パーティションを求めるための他の方法も可能である。反復処理ループは、次の行に移動し（ステップ１０７）、全ての行が処理されるまでループ（ステップ１０１〜１０７）を介した行の処理が継続する。ｅｉ〜がゼロよりも大きくない場合（ステップ１０５）、ルーチン１００は、ステップ１０７に直接移動する。行列の全ての行がループ（ステップ１０１〜１０７）を介して処理された後、ルーチン１００は終了する。図９を参照して示されるルーチン１００はまた、以下の擬似コードを使用して表現することができる。
ｆｏｒ ｉ＝０からｍ―１ ／＊疎行列Ａのｍ行にわたるループ＊／
ｊ←Ａ_ｒ［ｉ］ ／＊ｊ：ｉ行目における第１の非ゼロエントリのインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのインデックス＊／
ｅｉ〜←ｊ_ｍａｘ−ｊ ／＊ｅｉ〜：ｉ行目における非ゼロエントリ数＊／
（ｅｉ〜＞０）の場合
ｋ←ｌｏｗｅｒ＿ｂｏｕｎｄ（Ａｓ，ｅｉ〜） ／＊パーティションｋ ｓ．ｔ．Ａｓ［ｋ］≦ｅｉ〜＜Ａｓ［ｋ＋１］をみつける＊／
ｐａｒｔｉｔｉｏｎ［ｋ］．ｐｕｓｈ＿ｂａｃｋ（ｉ） ／＊パーティションｋの終わりに行ｉｄｉを追加（すなわち、Ｒｋ）＊／
ｉｆ文終了
ｆｏｒ文終了

同様に、ＣＳＣ又は構造化ＣＳＣについての順序維持行列パーティション化アルゴリズムを設計することができる。ｅ〜ｊを疎行列の列ｊにおける非ゼロエントリ数とする。Ｃｋと表されるｋ番目の列ベースのパーティションにおける列の設定は、Ｃｋ＝｛ｊ｜Ａｓ［ｋ］≦ｅ〜ｊ＜Ａｓ［ｋ＋１］｝と記載することができる。相対的順序を尊重することは、列インデックス配列Ａ’_ｃにアクセスする際のランダム性を低減することから、パーティション化アルゴリズムは、元のＣＳＣ行列における同じパーティション列の相対的順序を尊重し、構造化及び非構造化ＣＳＣ符号化の双方のキャッシュ性能を向上させる。図１１は、１つの実施形態にかかる図９の方法９０において使用するための圧縮表現で表される行列の列の順序維持パーティションを実行するためのルーチン１１０を示すフロー図である。圧縮表現は、ＣＳＣ及び構造化ＣＳＣの双方とすることができる。

行列の全ての列であるｊ＝０からｎ−１のｉｄを有する列を処理する反復処理ループ（ステップ１１１〜１１７）が開始される（ステップ１１１）。Ａ’_ｃ配列におけるｊ列目についてのエントリが配置され、行列におけるその列についての第１の非ゼロエントリについて、元の圧縮符号化が構造化ＣＳＣ又はＣＳＣであるかどうかに応じて、Ａ’_ｖｒ又はＡ’_ｖ配列におけるインデックスを識別し、識別されたインデックスは、変数ｉと設定される（ステップ１１２）。Ａ’_ｃ配列における次の（ｊ＋１）エントリが配置され、変数ｉ_ｍａｘと設定される（ステップ１１３）。ｊ列目が行列における最後の列でない限り、Ａ’ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒ（又はＡ’_ｖ）配列におけるインデックスであり、ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖｒ（又はＡ’_ｖ）配列におけるエントリの総数である。ｊ列目における非ゼロエントリの数ｅ〜ｊは、ｉ_ｍａｘからｉを差し引くことによって決定される（ステップ１１４）。ｅ〜ｊがゼロよりも大きい場合（ステップ１１５）、ｅ〜ｊの値に基づいて及びパーティションにおいて許容される非ゼロエントリの最大及び最小数に基づいて定義されたパーティションのいずれかに列ｊが割り当てられ、パーティションｋにおける列のリストの最後に列ｊの列ｉｄが追加される（ステップ１１６）。それゆえに、列ｊが属するインデックスｋを有するパーティションは、Ａ’ｓ［ｋ］≦ｅ〜ｊ＜Ａ’_ｓ［ｋ＋１］のようにみられる。１つの実施形態において、パーティションｋは、ｂと設定されたｅ〜ｊを有し、ｂ以上である昇順配列Ａ’_ｓの第１の要素のインデックスを返すｌｏｗｅｒ＿ｂｏｕｎｄ（Ａ’_ｓ，_ｂ）として示される関数を使用して求めることができる。パーティションを求めるための他の方法も可能である。反復処理ループは、次の列に移動し（ステップ１１７）、全ての列が処理されるまでループ（ステップ１１１〜１１７）を介した列の処理が継続する。ｅ〜ｊがゼロよりも大きくない場合（ステップ１１５）、ルーチン１１０は、ステップ１１７に直接移動する。行列の全ての列が処理された後、ルーチン１１０は終了する。図１１を参照して示されるルーチン１１０はまた、以下の擬似コードを使用して表現することができる。
ｆｏｒ ｊ＝０からｎ―１ ／＊疎行列Ａのｎ列にわたるループ＊／
ｉ←Ａ’_ｃ［ｊ］ ／＊ｉ：ｊ列目における第１の非ゼロエントリのインデックス＊／
ｉ_ｍａｘ←Ａ’_ｃ［ｊ＋１］ ／＊ｉ_ｍａｘ：（ｊ＋１）列目における第１の非ゼロエントリのインデックス＊／
ｅ〜ｊ←ｉ_ｍａｘ−ｉ ／＊ｅ〜ｊ：ｊ列目における非ゼロエントリ数＊／
（ｅ〜ｊ＞０）の場合
ｋ←ｌｏｗｅｒ＿ｂｏｕｎｄ（Ａ’_ｓ，_ｅ〜ｊ） ／＊パーティションｋ ｓ．ｔ．Ａ’_ｓ［ｋ］≦ｅ〜ｊ＜Ａ’_ｓ［ｋ＋１］をみつける＊／
ｐａｒｔｉｔｉｏｎ［ｋ］．ｐｕｓｈ＿ｂａｃｋ（ｊ） ／＊パーティションｋの終わりに列ｉｄｊを追加（すなわち、Ｃｋ）＊／
ｉｆ文終了
ｆｏｒ文終了

パーティションが作成されると、パーティションは、マッピング配列にマージされることができる。図１２は、１つの実施形態にかかる図９の方法９０において使用するためのマッピング配列にパーティションをマージするためのルーチン１２０を示すフロー図である。最初に、マッピング配列（Ａ_ｍ又はＡ’_ｍ）は、元の圧縮表現に応じて、配列のサイズ、配列における行列の部分の数をゼロに設定することによって初期化される（ステップ１２１）。作成された全てのパーティションであるインデックスｋ＝０からｐ−１を有するパーティションを処理する反復処理ループ（ステップ１２２〜１２６）が開始される（ステップ１２２）。パーティションｋについて、ｋよりも小さいインデックスを有する全てのパーティションの積算サイズが計算され、ｋが０である場合、積算サイズもまた０である（ステップ１２３）。ｋ未満のインデックスを有するパーティションは、配列におけるパーティションｋに先行することから、前のパーティションの積算サイズは、マッピング配列の最後がどこにあるかを示し、積算サイズを決定すると、パーティションｋは、マッピング配列の終わりに挿入される（ステップ１２４）。必要に応じて、パーティションｋによって占有される記憶装置３１内の任意のメモリは（１２５）まで解放される。ルーチン１２０は、次のパーティションに移動し（ステップ１２６）、全てのパーティションが処理されるまでループ（ステップ１２２〜１２６）を介したパーティションの処理が継続する。全ての作成されたパーティションを処理すると、全てのパーティションの積算サイズが計算されてＡｏ［ｐ］と設定され、ルーチン１２０を終了する。図１２のルーチン１２０はまた、以下の擬似コードを使用して示すことができる。
Ａ_ｍ．ｓｉｚｅ←０ ／＊Ａ_ｍ：（順序維持置換）マッピング配列＊／
ｆｏｒ ｋ＝０からｐ―１ ／＊ｐ個のパーティションにわたるループ＊／
Ａｏ［ｋ］←Ａ_ｍ．ｓｉｚｅ（） ／＊Ａｏ［ｋ］：全ての前のパーティションの積算サイズ＜ｋ＊／
Ａ_ｍｉｎｓｅｒｔ（ｐａｒｔｉｔｉｏｎ［ｋ］） ／＊Ａ_ｍの終わりにパーティションｋを挿入＊／
ｐａｒｔｉｔｉｏｎ［ｋ］を削除 ／＊パーティションｋによって使用される空きメモリ（必要に応じて）＊／
ｆｏｒ文終了
Ａｏ［ｐ］←Ａ_ｍ．ｓｉｚｅ（） ／＊Ａｏ［ｐ］：全てのパーティションの総サイズ＊／

図１２を参照して示されたルーチン１２０に戻ると、Ａｏ［ｋ］は、パーティションｋの第１の行（又は列）のランクを格納する。パーティションｋの最後の行（又は列）のランクは、ｋ＝０，１，・・・，ｐ−１について、Ａｏ［ｋ＋１］−１によって与えられる。配列オフセットＡｏ［ｐ］の最後の要素は、ＣＳＲ（又はＣＳＣ）についての行（又は列）の数である置換配列の要素数に常に等しい。

行列のパーティション化された圧縮表現は、ＣＰＵのみを含むものと同様にＣＰＵ及びＧＰＵの双方を含むサーバ３３による処理のために使用することができる。図１３は、１つの実施形態にかかる図９の方法９０において使用するためのパーティション化された圧縮表現で符号化された行列に対してＳｐＭＶを実行するためのルーチン１３０を示すフロー図である。ルーチンは、ＰＳＣＳＲ及びＰＣＳＲ表現の双方に適用することができる。ルーチン１３０は、ＣＰＵのみを含むものと同様にＣＰＵ及びＧＰＵの双方を含むサーバによって使用することができる。作成された全てのパーティションであるインデックスｋ＝０からＡｐ−１を有するパーティションを処理する外側の反復処理ループ（ステップ１３１〜１４２）が開始される（ステップ１３１）。ｋ番目のパーティションにおける全ての行を処理する内側の反復処理ループ（ステップ１３２〜１４１）が開始され、マッピング配列における行のランクは、ｒ＝Ａｏ［ｋ］からＡｏ［ｋ＋１］−１である（ステップ１３２）。ｒ番目にランク付けされた行のｉｄは、ｉと識別されて設定される（ステップ１３３）。Ａ_ｒ配列におけるｉ行目についてのエントリが配置され、初期の圧縮表現がＣＳＲ又は構造化ＣＳＲであるかどうかに応じて、行列におけるその行についての第１の非ゼロエントリのＡｖ又はＡ_ｖｃ配列におけるインデックスを識別し、識別されたインデックスが変数ｊとして設定される（ステップ１３４）。Ａ_ｒ配列における次の（ｉ＋１）エントリが配置され、変数ｊ_ｍａｘと設定される（ステップ１３５）。ｉ行目が行列における最後の行でない限り、Ａ_ｒ配列における次のエントリは、（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃ（又はＡｖ）配列におけるインデックスであり、ｉ行目が行列における最後の行である場合、Ａ_ｒ配列における次のエントリは、Ａ_ｖｃ（又はＡｖ）配列におけるエントリの総数である。ステップ１３８において後述する非ゼロ配列の値の乗算結果を加算する関数である総和計算部は、ゼロにおける合計の値を設定することによって初期化される（ステップ１３６）。ｊがｊ_ｍａｘ未満である場合（ステップ１３７）、ＳｐＭＶが実行される符号化がＰＳＣＳＲ又はＰＣＳＲであるかどうかに応じた動作で、ｊの値に対して乗算及び加算計算が実行される（ステップ１３８）。符号化がＰＳＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］

計算において、インデックスｊを有するＡ_ｖｃ配列の要素に格納された値はｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｖｃ要素に格納された列のインデックスであり、乗算結果は、ステップ１３８の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａｖ［ｊ］×ｘ［Ａ_ｃ［ｊ］］
ここで、インデックスｊを有するＡｖ配列における値は、ｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｃ配列における数であり、乗算結果は、ステップ１３８の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｊの値に１が加算され、加算結果がｊと設定され（ステップ１３９）、その行における次の列のエントリに処理を移動する。ルーチン１３０は、上述したステップ１３７に戻り、ｉ行目における非ゼロ値が処理されるまでステップ１３７〜１３９を繰り返す。ｊがｊ_ｍａｘ以上である場合（ステップ１３７）、ループ１３７〜１３９における反復中に乗算結果を加算した合計は、密ベクトルｙに格納される（ステップ１４０）。反復処理ループは、次の行に移動し（ステップ１４１）、パーティションの全ての行がステップ１３２〜１４１を介して処理されるまで内側ループ（ステップ１３２〜１４１）を介した行の処理が継続する。ｋ番目のパーティションの全ての行が処理されると、ルーチン１３０は、次のパーティションに移動し（ステップ１４２）、全てのパーティションが処理されるまで外側処理ループ（ステップ１３１〜１４２）を介したパーティションの処理が継続する。全てのパーティションがステップ１３１〜１４２において処理されると、ルーチン１３０は終了する。ルーチン１３０はまた、ＰＳＣＳＲ符号化に対してＳｐＭＶを実行することを含む以下の擬似コードを使用して表現することができる。
ｆｏｒ ｋ＝０からＡｐ―１ ／＊Ａｐ個の行ベースのパーティションにわたるループ＊／
ｆｏｒ ｒ＝Ａｏ［ｋ］からＡｏ［ｋ＋１］−１ ／＊ｋ番目のパーティションにおける行にわたるループ＊／
ｉ←Ａ_ｍ［ｒ］ ／＊ｉ：ｒ番目にランク付けされた行のｉｄ＊／
ｊ←Ａ_ｒ［ｉ］ ／＊ｊ：ｉ行目における第１の非ゼロエントリのＡ_ｖｃにおけるインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器の初期化＊／
ｗｈｉｌｅ（ｊ＜ｊ_ｍａｘ） ／＊行の終わりに到達したかどうかのテスト＊／
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］ ／＊ｙｉ＝ΣｊＡｉ，ｊ×ｘｊの計算
ｊ←ｊ＋１ ／＊ｉ行目における次の非ゼロ列に移動＊／
ｗｈｉｌｅ文終了
ｙ［ｉ］←ｓｕｍ ／＊ｙに結果を格納＊／
ｆｏｒ文終了
ｆｏｒ文終了
ルーチン１３０はまた、準用するＰＣＳＲ符号化に対してＳｐＭＶを実行するための擬似コードを使用して表現することができる。

ＳｐＭＴＶは、ＣＰＵのみを含むものと同様にＣＰＵ及びＧＰＵの双方を含むサーバを使用してパーティション化された圧縮表現に対して行うことができる。図１４は、図９の方法９０において使用するための行列のパーティション化された圧縮表現に対してＳｐＭＴＶを実行するためのルーチン１５０を示すフロー図である。ルーチンは、ＰＳＣＳＣ及びＰＣＳＣ表現の双方に適用することができる。ルーチン１５０は、ＣＰＵのみを含むものと同様にＣＰＵ及びＧＰＵの双方を含むサーバを使用することができる。作成された全てのパーティションであるインデックスｋ＝０からＡ’_ｐ−１のパーティションを処理する外側の反復処理ループ（ステップ１５１〜１６２）が開始される（ステップ１５１）。ｋ番目のパーティションにおける全ての列を処理する内側の反復処理ループ（ステップ１５２〜１６１）が開始され、処理されるマッピング配列における列のランクは、ｒ＝Ａ’_ｏ［ｋ］からＡ’_ｏ［ｋ＋１］−１である（ステップ１５２）。ｒ番目にランク付けされた列のｉｄは、ｊと識別されて設定される（ステップ１５３）。Ａ’_ｒ配列におけるｊ列目についてのエントリが配置され、圧縮表現がＣＳＣ又は構造化ＣＳＣであるかどうかに応じて、行列におけるその列についての第１の非ゼロエントリのＡ’_ｖ又はＡ’_ｖｒ配列におけるインデックスを識別し、識別されたインデックスが変数ｉとして設定される（ステップ１５４）。Ａ’_ｃ配列における次の（ｊ＋１）エントリが配置され、変数ｉ_ｍａｘと設定される（ステップ１５５）。ｊ列目が行列における最後の列でない限り、Ａ’_ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒ（又はＡ’_ｖ）配列におけるインデックスであり、ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖｒ（又はＡ’_ｖ）配列におけるエントリの総数である。ステップ１５８において後述する非ゼロ配列の値の乗算結果を加算する関数である総和計算部は、ゼロにおける合計の値を設定することによって初期化される（ステップ１５６）。ｉがｉ_ｍａｘ未満である場合（ステップ１５７）、ＳｐＭＴＶが実行される符号化がＰＳＣＳＣ又はＰＣＳＣであるかどうかに応じた動作で、ｉの値に対して乗算及び加算計算が実行される（ステップ１５８）。符号化がＰＳＣＳＣである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖｒ［ｉ］．ｖ×ｘ’［Ａ’_ｖｒ［ｉ］．ｒ］
ここで、インデックスｉを有するＡ’_ｖｒ配列の要素に格納された値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｖｒ要素に格納された行のインデックスであり、乗算結果は、ステップ１５８の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］
ここで、インデックスｉを有するＡ’_ｖ配列における値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｒ配列における数であり、乗算結果は、ステップ１５８の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｉの値に１が加算され、加算結果がｉと設定され（ステップ１５９）、その列における次の行のエントリに処理を移動する。ルーチン１５０は、上述したステップ１５７に戻り、ｊ列目における非ゼロ値が処理されるまでステップ１５８〜１５９を繰り返す。ｉがｉ_ｍａｘ以上である場合（ステップ１５７）、ループ１５７〜１５９における反復中に乗算結果を加算した合計は、密ベクトルｙ’に格納される（ステップ１６０）。反復処理ループは、次の列に移動し（ステップ１６１）、全ての列がステップ１５２〜１６１を介して処理されるまで内側処理ループ（ステップ１５２〜１６１）を介した列の処理が継続する。ｋ番目のパーティションの全ての列が処理されると、ルーチン１５０は、次のパーティションに移動し（ステップ１６２）、全てのパーティションが処理されるまで外側処理ループ（ステップ１５１〜１６２）を介したパーティションの処理が継続する。全てのパーティションがステップ１５１〜１６２において処理されると、ルーチン１５０は終了する。ルーチン１５０はまた、ＰＳＣＳＲ符号化に対してＳｐＭＴＶを実行することを含む以下の擬似コードを使用して表現することができる。
ｆｏｒ ｋ＝０からＡ’_ｐ―１ ／＊Ａ’_ｐ個の列ベースのパーティションにわたるループ＊／
ｆｏｒ ｒ＝Ａ’_ｏ［ｋ］からＡ’_ｏ［ｋ＋１］−１ ／＊ｋ番目のパーティションにおける列にわたるループ＊／
ｊ←Ａ’_ｍ［ｒ］ ／＊ｊ：ｒ番目にランク付けされた列のｉｄ＊／
ｉ←Ａ’_ｃ［ｊ］ ／＊ｉ：ｊ列目における第１の非ゼロエントリのＡ’_ｖｒにおけるインデックス＊／
ｉ_ｍａｘ←Ａ’_ｃ［ｊ＋１］ ／＊ｉ_ｍａｘ：（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器の初期化＊／
ｗｈｉｌｅ（ｉ＜ｉ_ｍａｘ） ／＊列の終わりに到達したかどうかのテスト＊／
ｓｕｍ←ｓｕｍ＋Ａ’_ｖｒ［ｉ］．ｖ×ｘ’［Ａ’_ｖｒ［ｉ］．ｒ］ ／＊ｙ’ｊ＝ΣｉＡｉ，ｊ×ｘ’ｉの計算
ｉ←ｉ＋１ ／＊ｊ列目における次の非ゼロ行に移動＊／
ｗｈｉｌｅ文終了
ｙ’［ｊ］←ｓｕｍ ／＊ｙ’に結果を格納＊／
ｆｏｒ文終了
ｆｏｒ文終了

上記図１３及び図１４を参照して上述したルーチンは、それらを実行するための特定のハードウェアセットアップに固有のものではなかった。ＧＰＵを含むハードウェアは、特定のパーティションに対してＳｐＭＶ及びＳｐＭＴＶを実行するために特定のカーネルを割り当てることによってより高速になる可能性がある。図１５は、１つの実施形態にかかる図９の方法９０において使用するためのパーティション化された圧縮表現処理のためのカーネル関数を起動するためのルーチンである。カーネル関数Ｆの配列は、図３を参照して上述したカーネルのリストを含むＦ＝［ｆ１Ｔ１Ｒ，ｆ１Ｗ１Ｒ，ｆ１Ｂ１Ｒ］のように作成され、カーネルは、０から２にインデックス付けされ、ｆ１Ｔ１Ｒはインデックス０を有し、ｆ１Ｂ１Ｒはインデックス２を有する（ステップ１７１）。さらなる実施形態において、他のカーネルは、配列の一部とすることができる。作成された全てのパーティションであるインデックスｋ＝０からＡｐ−１を有するパーティションについての反復処理ループ（ステップ１７２〜１７６）が開始される（ステップ１７２）。カーネル関数の１つは、ｋの値及びＦ配列におけるカーネルのインデックスに基づいてパーティションｋについて選択される（ステップ１７３）。例えば、ｋ＝０の場合、０番目のインデックスを有するカーネルｆ１Ｔ１Ｒが選択され、ｋ＝１の場合、１のインデックスを有するカーネルｆ１Ｗ１Ｒが選択され、ｋ＝２の場合、２のインデックスを有するカーネルｆ１Ｂ１Ｒが選択される。３超のパーティションがある場合、２よりも大きいインデックスｋを有する全てのパーティションについてｆ１Ｂ１Ｒが選択される。パーティション及びカーネルについての他のインデックスもまた使用可能であり、カーネル及びパーティションを一致させる他の方法も可能である。起動引数は、さらに図１６を参照して記載されるように、選択されたカーネルについて設定される（ステップ１７４）。選択された関数は、各引数について起動され、さらに図１７〜図２２を参照して後述するように、パーティション化された圧縮符号化の処理に使用される（ステップ１７５）。反復処理ループは、次のパーティションに移動し（ステップ１７６）、全てのパーティションが処理されるまでループ（ステップ１７２〜１７６）を介したパーティションの処理が継続する。全てのパーティションが処理されると、ルーチン１７０は終了する。ルーチン１７０はまた、以下の擬似コードを使用して表現することができる−擬似コードは、行ベースのパーティションについて記載されているが、列ベースのパーティションについての擬似コードを準用して記載することができる。
ｆｏｒ ｋ＝０からＡｐ−１ ／＊Ａｐ個の行ベースのパーティションにわたるループ＊／
ａｒｇｓ←ｋｅｒｎｅｌ＿ｌａｕｎｃｈ＿ａｒｇｓ（Ｆ［ｋ］，Ａｏ［ｋ＋１］−Ａｏ［ｋ］） ／＊Ｆ［ｋ］：ｋ番目のカーネル関数＊／
Ｆ［ｋ］＜＜＜ａｒｇｓ＞＞＞（ｙ，ｘ，Ａ_ｖｃ，Ａ_ｒ，Ａ_ｍ，Ａｏ［ｋ］，Ａｏ［ｋ＋１］） ／＊ｋ番目のパーティションについてのｋ番目のカーネルを起動＊／
ｆｏｒ文終了
擬似コードは、ＰＳＣＳＲを参照して記載されているが、他のパーティション化された圧縮表現についての擬似コードを準用して表すことができる。図１５を参照して上述したもの以外の特定のパーティションを処理するための特定のカーネル関数を選択する方法もまた可能である。

カーネルが選択された後に設定された起動引数は、パーティション化された圧縮表現の処理を実行するサーバ３３などのシステムの実行時に入力される制約を提供する。ＧＰＵ上でカーネルを起動するために、スレッドブロックの寸法及び可能であればそのようなスレッドブロックのグリッドの寸法などの起動引数を指定する必要がある。同じベンダからの異なるＧＰＵは、ブロック又は各寸法に沿ったグリッドの最大サイズについての制約を課すことがあり、異なるＣＰＵベンダは、異なる数のスレッドブロック及び／又はグリッド寸法をサポートすることがある。図１６は、１つの実施形態にかかる図１５のルーチン１７０において使用するための選択されたカーネル関数の起動引数を設定するためのルーチン１８０である。ｎを選択された行列に割り当てられた行列の部分、行又は列の数とする。選択されたカーネル関数がｆ１Ｔ１Ｒである場合（ステップ１８１）、ｎが起動するスレッドの最小数であるという制約が設定される。残りの引数が設定されて制約に基づいて返され（ステップ１８３）、ルーチン１８０を終了する。選択された関数がｆ１Ｔ１Ｒでない場合（ステップ１８１）、選択された関数がｆ１Ｗ１Ｒであるかどうかが判定される（ステップ１８４）。選択された関数がｆ１Ｗ１Ｒである場合（ステップ１８４）、ブロックサイズがワープサイズに等しく設定され且つｎが起動するブロックの最小数であるという制約が設定される（１８５）。ルーチン１８０は、後述するステップ１８３に移動する。選択された関数がｆ１Ｗ１Ｒでない場合（ステップ１８１）、選択された関数がｆ１Ｂ１Ｒであるかどうかが判定される（ステップ１８６）。選択された関数がｆ１Ｂ１Ｒである場合（ステップ１８７）、ブロックサイズは、ブロック内のスレッド数に設定され、ｎが起動するブロックの最小数に設定される。ルーチン１８０は、後述するステップ１８３に移動する。選択された関数がｆ１Ｂ１Ｒでない場合（ステップ１８６）、カーネル関数は未知であり、ルーチン１８０は終了する。ルーチン１８０はまた、以下の擬似コードを使用して表現することができる−擬似コードは、ＳｐＭＶを実行するためのカーネルの起動について記載されているが、ＳｐＭＴＶを実行するためのカーネルの起動についての擬似コードを準用して記載することができる。
関数 ｋｅｒｎｅｌ＿ｌａｕｎｃｈ＿ａｒｇｓ（ｆ，ｎ） ／＊ｆ：ＳｐＭＶカーネル；ｎ：ｆに割り当てられた行数＊／
（ｆ＝ｆ１Ｔ１Ｒ）である場合 ／＊ｆ１Ｔ１Ｒ：１スレッド１行カーネル＊／
ａｒｇｓ．ｓｅｔ＿_ｍｉｎ＿ｔｈｒｅａｄｓ（ｎ） ／＊ｎ＝起動されるスレッドの最小数＊／
（ｆ＝ｆ１Ｗ１Ｒ）である場合 ／＊ｆ１Ｗ１Ｒ：１ワープ１行カーネル＊／
ａｒｇｓ．ｓｅｔ＿ｂｌｏｃｋ＿ｓｉｚｅ（ＷＡＲＰＳＩＺＥ） ／＊ＢＬＯＣＫＳＩＺＥ＝ＷＡＲＰＳＩＺＥに設定＊／
ａｒｇｓ．ｓｅｔ＿_ｍｉｎ＿ｂｌｏｃｋｓ（ｎ） ／＊ｎ＝起動されるブロックの最小数＊／
（ｆ＝ｆ１Ｂ１Ｒ）である場合 ／＊ｆ１Ｂ１Ｒ：１ブロック１行カーネル＊／
ａｒｇｓ．ｓｅｔ＿ｂｌｏｃｋ＿ｓｉｚｅ（ＢＬＯＣＫＳＩＺＥ） ／＊ＢＬＯＣＫＳＩＺＥ：ブロックにおけるスレッド数＊／
ａｒｇｓ．ｓｅｔ＿_ｍｉｎ＿ｂｌｏｃｋｓ（ｎ） ／＊ｎ＝起動されるブロックの最小数＊／
それ以外の場合
エラー「中断：未知のカーネル関数」
ｉｆ文終了
ａｒｇｓ．ｃｏｍｐｕｔｅ＿ｓａｔｉｓｆｙ＿ａｒｇｓ（） ／＊上記制約に基づいて残りの引数を設定＊／
リターン ａｒｇｓ

固定されたスレッド−ブロック及びグリッド寸法を使用する代わりに、上記擬似コードは、３つのＳｐＭＶカーネルのそれぞれによって課される明示的な制約を形成することにより、制約を満たす問題を設定する起動引数に接近する。例えば、ｆ１Ｔ１Ｒカーネルは、スレッドの総数がカーネルに割り当てられた行数に少なくとも等しくなければならないという単一の制約の充足を必要とするのみである。一方、ｆ１Ｗ１Ｒカーネルは、（１）ブロック内のスレッド数がＷＡＲＰＳＩＺＥと同じでなければならない、及び（２）カーネルに割り当てられた行があるのと少なくとも同数のスレッドブロックがなければならないという２つの制約を同時に充足することを要求する。ハードウェアによって課される制約を尊重しつつ、それらのカーネル起動制約を満たす方法は、ＧＰＵ又はベンダ依存とすることができ、具体的な説明は提供されない。通常、最良の起動引数は、２の整数累乗である値を想定し、それゆえに、そのような制約充足問題の探索空間は、通常、非常に小さく且つ扱いやすい。例えば、Ｎｖｉｄｉａ（登録商標）のＦｅｒｍｉ ＧＰＵ上の実装は、ｆ１Ｔ１Ｒカーネルについては８×８のスレッドブロックサイズを使用し、ｆ１Ｗ１Ｒカーネルについては３２×１のブロックサイズを使用し、ｆ１Ｂ１Ｒカーネルについては５１２×１のブロックサイズを使用する。

選択されると、カーネル関数は、行列の割り当てられた部分を処理することができる。行列パーティション化の目的は、同じパーティションにおける行がカーネルが最適化された非ゼロエントリ数などのいくつかの共通の特徴を共有するように、単一のＧＰＵカーネル関数について１つのパーティションを作成することであることを思い出されたい。上述したＧＰＵベースのＳｐＭＶアルゴリズムは、任意数のパーティション及びカーネルを扱うことができる一方で、実験は、ＧＴＸ４８０及び５８０を含むＮｖｉｄｉａ（登録商標）のＦｅｒｍｉクラスのＧＰＵにおいて最良のＳｐＭＶ結果を生成する３つの異なるカーネル（ｆ１Ｔ１Ｒ，ｆ１Ｗ１Ｒ，ｆ１Ｂ１Ｒ）を有することを示している。

ＷＡＲＰＳＩＺＥ（本特許出願において「ワープサイズ」とも称される）及びＢＬＯＣＫＳＩＺＥ（本特許出願において「ブロックサイズ」とも称される）を、それぞれ、ワープ及びスレッドブロックにおけるスレッド数とする。ＧＰＵコアに対する行列の行（又は列）の分布を容易とするために、図１７〜図２２Ｂを参照して以下に記載されるルーチンにおいて以下のヘルパー関数を使用可能である。
・ｔｈｒｅａｄ＿ｉｄ（）は、現在のスレッドについてのグローバルに固有のスレッドｉｄを返す
・ｗａｒｐ＿ｉｄ（）は、現在のスレッドについてのグローバルに固有のワープｉｄを返す
・ｗａｒｐ＿ｔｈｒｅａｄ＿ｉｄ（）は、現在のスレッドについての（ワープ内でのみ固有の）ローカルスレッドｉｄを返す
・ｂｌｏｃｋ＿ｔｈｒｅａｄ＿ｉｄ（）は、現在のスレッドについての（ブロック内でのみ固有の）ローカルスレッドｉｄを返す
・ｓｙｎｃ＿ｗａｒｐ＿ｔｈｒｅａｄｓ（）は、ワープ内の全てのスレッド間で同期する
・ｓｙｎｃ＿ｂｌｏｃｋ＿ｔｈｒｅａｄｓ（）は、ブロック内の全てのスレッド間で同期する

全てのグローバル（又はローカル）スレッド（又はワープ）ｉｄは、ゼロで始まるスレッドのインデックスを有するゼロベースである。ＷＡＲＰＳＩＺＥは、ハードウェアに依存する：Ｎｖｉｄｉａ（登録商標）については３２個又はアドバンスト・マイクロ・デバイス（登録商標）については６４個のＧＰＵ。一方、ＢＬＯＣＫＳＩＺＥは、通常はＷＡＲＰＳＩＺＥの整数倍であるプログラム可能なパラメータとすることができる。共有されるキーワードは、同じブロック（又はＢＬＯＣＫＳＩＺＥ＝ＷＡＲＰＳＩＺＥである場合にはワープ）内の全てのスレッド間で共有されるスレッド変数を宣言する。

図１７は、１つの実施形態にかかる図１５のルーチンにおいて使用するためのｆ１Ｔ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチン１９０を示すフロー図である。ルーチン１９０は、ＰＳＣＳＲ又はＰＣＳＲ符号化のいずれかでＳｐＭＶを実行するために使用することができる。カーネルにおける全ての起動されるスレッドについて反復処理ループが開始され、スレッドは並列に動作し、スレッドは並列に動作し、それゆえに、以下のステップは、１つのスレッドを参照して説明され、任意の他の起動スレッドは、並列にステップを行う（ステップ１９１）。ｔｈｒｅａｄ＿ｉｄ（）関数を使用することによってスレッドの固有ｉｄが取得され、Ａｏ［ｋ］に格納されたパーティションｋ（スレッドが割り当てられたパーティション）の第１の行のランクに等しい変数ｒ_ｍｉｎの値に追加され、加算結果は、そのスレッドに割り当てられた行のランクを識別する変数ｒによって表される（１９２）。ｒの値は、Ａｏ［ｋ＋１］によって与えられる次のパーティションｋ＋１の第１の行のランクに等しい変数ｒ_ｍａｘと比較することができ、起動されるスレッド数がパーティションｋにおける行数（ｒ_ｍａｘ−ｒ_ｍｉｎ）に等しい場合には比較は任意である（１９３）。ｒがｒ_ｍａｘ未満である場合（ステップ１９３）、ｒ番目のランク付けされた行のｉｄｉは、マッピング配列Ａ_ｍにおいて識別される（ステップ１９４）。ｒがｒ_ｍａｘ以上である場合、ルーチン１９０は、後述するステップ２０２に移動する。Ａｒ配列におけるｉ行目についてのエントリが配置され、初期の圧縮表現がＣＳＲ又は構造化ＣＳＲであったかどうかに応じて、行列におけるその行についての第１の非ゼロエントリのＡｖ又はＡ_ｖｃ配列におけるインデックスを識別し、識別されたインデックスは、変数ｊとして設定される（ステップ１９５）。Ａ_ｒ配列における次の（ｉ＋１）エントリが配置され、変数ｊ_ｍａｘとして設定される（ステップ１９６）。ｉ行目が行列における最後の行でない限り、Ａ_ｒ配列における次のエントリは、（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃ（又はＡｖ）配列におけるインデックスである。ｉ行目が行列における最後の行である場合、Ａ_ｒ配列における次のエントリは、Ａ_ｖｃ（又はＡｖ）配列におけるエントリの総数である。以下のステップ１９９において記載される非ゼロ配列の値の乗算結果をともに加算する関数である総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ１９７）。ｊがｊ_ｍａｘ未満である場合（ステップ１９８）、ＳｐＭＶが行われた符号化がＰＳＣＳＲ又はＰＣＳＲであるかどうかに応じて、乗算及び加算計算がｊの値に対して実行される（ステップ１９９）。符号化がＰＳＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］

計算において、インデックスｊを有するＡ_ｖｃ配列の要素に格納された値はｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｖｃ要素に格納された列のインデックスであり、乗算結果は、ステップ１９９の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａｖ［ｊ］×ｘ［Ａ_ｃ［ｊ］］
ここで、インデックスｊを有するＡｖ配列における値は、ｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｃ配列における数であり、乗算結果は、ステップ１９９の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｊの値に１が加算され、加算結果がｊと設定され（ステップ２００）、その行における次の列のエントリに処理を移動する。ルーチン１９０は、上述したステップ１９８に戻り、ｉ行目における非ゼロ値が処理されるまでステップ１９９〜２０１を繰り返す。ｊがｊ_ｍａｘ以上である場合（ステップ１９８）、ループ１９８〜２００における反復中に乗算結果を加算した合計は、密ベクトルｙに格納される（ステップ２０１）。スレッドの処理の実行は停止され（ステップ２０２）、ルーチン１９０は、全ての起動したスレッドの処理の停止によって終了する。ルーチン１９０はまた、以下の擬似コードを使用して表現することができ−擬似コードは、ＰＳＣＳＲ符号化に関して記載されるが、ＰＣＳＲ符号化についての擬似コードを準用して記載することができる。
カーネル ｆ１Ｔ１Ｒ（ｙ，ｘ，Ａ_ｖｃ，Ａ_ｒ，Ａ_ｍ，ｒ_ｍｉｎ，ｒ_ｍａｘ） ／＊１Ｔ１Ｒ：１スレッド１行ＳｐＭＶカーネル＊／
ｒ←ｒ_ｍｉｎ＋ｔｈｒｅａｄ＿ｉｄ（） ／＊このスレッドに割り当てられた頂点のランクを計算＊／
（ｒ＜ｒ_ｍａｘ）である場合 ／＊任意の正確に（ｒ_ｍａｘ−ｒ_ｍｉｎ）スレッドが作成された場合＊／
ｉ←Ａ_ｍ［ｒ］ ／＊ｉ：ｒ番目にランク付けされた行のｉｄ＊／
ｊ←Ａ_ｒ［ｉ］ ／＊ｊ：ｉ行目における第１の非ゼロエントリのＡ_ｖｃにおけるインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのＡ_ｖｃにおけるインデックス＊／
ｓｕｍ←０ ／＊総和積算器の初期化＊／
ｗｈｉｌｅ（ｊ＜ｊ_ｍａｘ） ／＊行の終わりに到達したかどうかのテスト＊／
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］ ／＊ｙｉ＝ΣｊＡｉ，ｊ×ｘｊの計算
ｊ←ｊ＋１ ／＊ｉ行目における次の非ゼロ列に移動＊／
ｗｈｉｌｅ文終了
ｙ［ｉ］←ｓｕｍ ／＊ｙに結果を格納＊／
ｉｆ文終了

図１８は、１つの実施形態にかかる図１５のルーチン１７０において使用するためのｆ１Ｔ１Ｒカーネル関数によってＳｐＭＴＶを実行するためのルーチン２１０を示すフロー図である。ルーチン２１０は、図１７を参照して上述したように、ＰＳＳＣＲ又はＰＣＳＲ符号化のいずれかにおいてＳｐＭＴＶを実行するために使用することができ、全ての起動したスレッドは並列に動作する。反復処理ループは、カーネル内の全ての起動したスレッドについて開始され、スレッドは並列に動作する（ステップ２１１）。スレッドの固有ｉｄがｔｈｒｅａｄ＿ｉｄ（）関数を使用して取得され、Ａ’_ｏ［ｋ］に格納されたパーティションｋ（スレッドが割り当てられたパーティション）の第１の列のランクに等しい変数ｒ_ｍｉｎの値に加算され、加算結果は、変数ｒによって表される。ｒの値は、Ａ’_ｏ［ｋ＋１］によって与えられる次のパーティションｋ＋１の第１の列のランクに等しい変数ｒ_ｍａｘと比較することができ、起動されるスレッド数がパーティションｋにおける列数（ｒ_ｍａｘ−ｒ_ｍｉｎ）に等しい場合には比較は任意である。ｒがｒ_ｍａｘ未満である場合（ステップ２１３）、ｒ番目のランク付けされた列のｉｄｊは、マッピング配列Ａ’_ｍにおいて識別される（ステップ２１４）。ｒがｒ_ｍａｘ以上である場合、ルーチン２１０は、後述するステップ２２２に移動する。Ａ’_ｃ配列におけるｊ列目についてのエントリが配置され、圧縮表現がＣＳＣ又は構造化ＣＳＣであったかどうかに応じて、行列におけるその列についての第１の非ゼロエントリのＡ’_ｖ又はＡ’_ｖｒ配列におけるインデックスを識別し、識別されたインデックスは、変数ｉとして設定される（ステップ２１５）。Ａ’_ｃ配列における次の（ｊ＋１）エントリが配置され、変数ｉ_ｍａｘとして設定される（ステップ２１６）。ｊ列目が行列における最後の列でない限り、Ａ’_ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒ（又はＡ’_ｖ）配列におけるインデックスである。ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖｒ（又はＡ’_ｖ）配列におけるエントリの総数である。以下のステップ２１９において記載される非ゼロ配列の値の乗算結果をともに加算する関数である総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ２１７）。ｉがｉ_ｍａｘ未満である場合（ステップ２１８）、ＳｐＭＴＶが行われた符号化がＰＳＣＳＣ又はＰＣＳＣであるかどうかに応じて、乗算及び加算計算がｉの値に対して実行される（ステップ２１９）。符号化がＰＳＣＳＣである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖｒ［ｉ］．ｖ×ｘ’［Ａ’_ｖｒ［ｉ］．ｒ］
ここで、インデックスｉを有するＡ’_ｖｒ配列の要素に格納された値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｖｒ要素に格納された行のインデックスであり、乗算結果は、ステップ２１９の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］
ここで、インデックスｉを有するＡ’_ｖ配列における値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｒ配列における数であり、乗算結果は、ステップ２１９の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｉの値に１が加算され、加算結果がｉと設定され（ステップ２２０）、その列における次の行のエントリに処理を移動する。ルーチン２１０は、上述したステップ２１８に戻り、ｊ列目における非ゼロ値が処理されるまでステップ２１８〜２２０を繰り返す。ｉがｉ_ｍａｘ以上である場合（ステップ２１８）、ループ２１８〜２２０における反復中に乗算結果を加算した合計は、密ベクトルｙ’に格納される（ステップ２２１）。スレッドの実行が停止される（ステップ２２２）。ルーチン２１０は、全ての起動した処理スレッドの実行が停止されると終了する。ルーチン２１０についての擬似コードは、準用する図１７を参照して上記に示された擬似コードと同様に記載することができる。

非ゼロエントリ数が多い行列の部分は、各行などの各部分を処理するためにスレッドのワープを割り当てるｆ１Ｗ１Ｒカーネルによって処理されることから利益を得ることができる。図１９Ａ〜図１９Ｂは、１つの実施形態にかかる図１５のルーチン１７０において使用するためのｆ１Ｗ１Ｒカーネル関数によってＳｐＭＶを実行するためのルーチン２３０を示すフロー図である。ルーチン２３０は、ＰＳＳＣＲ又はＰＣＳＲ符号化のいずれかにおいてＳｐＭＶを行うために使用することができる。カーネルにおける全ての起動ワープについて反復処理ループが開始される（ステップ２３１）。起動ワープ及び起動ワープにおけるスレッドは、互いに並列に実行される。それゆえに、以下の説明がワープのいずれか又はスレッドのいずれかを参照する場合には、全ての他の起動ワープ又はスレッドは、スレッド又はワープのいずれかの実行が停止されるまで（以下のステップ２５１など）参照されるワープ又はスレッドと並列にルーチン２３０の同じステップを実行する。

ローカル共有は、同じワープにおける全てのスレッド間で開始され、スレッドがステップ２４０〜２４２において後述する計算結果を共有するのを可能とする（ステップ２３２）。ワープのいずれかの固有のグローバルｉｄを取得するために関数ｗａｒｐ＿＿ｉｄ（）を使用し且つｒ_ｍｉｎの値に固有のグローバルｉｄを追加することによってワープのいずれかに割り当てられた行のランク（図１７を参照して上述したようにランクが割り当てられるパーティションの第１の行のランク）が取得され、加算結果は、ワープに割り当てられた行のランクであり、変数ｒを使用して表される（ステップ２３３）。値ｒは、Ａｏ［ｋ＋１］によって与えられる次のパーティションｋ＋１の第１の行のランクに等しい変数であるｒ_ｍａｘと比較することができ、起動されるワープ数がパーティションｋにおける行数（ｒ_ｍａｘ−ｒ_ｍｉｎ）に等しい場合には比較は任意である。ｒがｒ_ｍａｘ未満である場合（ステップ２３４）、ｒ番目のランク付けされた行のｉｄｉは、マッピング配列Ａ_ｍにおいて識別される（ステップ２３５）。ｒがｒ_ｍａｘ以上である場合（ステップ２３４）、ルーチン２３０は、後述するステップ２５１に移動する。関数ｗａｒｐ＿ｔｈｒｅａｄ＿ｉｄ（）を使用することなどによってワープにおけるスレッドのいずれかのローカルｉｄ（ワープ内のｉｄ）が取得され、変数ｔを使用して表される（ステップ２３６）。そのスレッドに割り当てられた第１の非ゼロエントリのインデックスは、Ａ_ｒ配列におけるｉ行目についてのエントリの値をｔに追加することによって取得され、加算結果は、変数ｊを使用して表される（ステップ２３７）。Ａ_ｒ配列における次の（ｉ＋１）エントリが配置され、変数ｊ_ｍａｘとして設定される（ステップ２３８）。以下のステップ２４１において記載される非ゼロ配列の値の乗算結果をともに加算する関数であるローカル総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ２３９）。ｊがｊ_ｍａｘ未満である場合（ステップ２４０）、ＳｐＭＶが行われた符号化がＰＳＣＳＲ又はＰＣＳＲであるかに応じた動作で、値ｊに対して乗算及び加算計算が実行される（ステップ２４１）。符号化がＰＳＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］

計算において、インデックスｊを有するＡ_ｖｃ配列の要素に格納された値はｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｖｃ要素に格納された列のインデックスであり、乗算結果は、ステップ２４１の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａｖ［ｊ］×ｘ［Ａ_ｃ［ｊ］］
ここで、インデックスｊを有するＡｖ配列における値は、ｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｃ配列における数であり、乗算結果は、ステップ２４１の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｊの値にワープサイズが加算され、加算結果がｊと設定され（ステップ２４２）、スレッドが割り当てられたその行における次の列のエントリに処理を移動する。ルーチン２３０は、上述したステップ２４０に戻り、そのスレッドが割り当てられるｉ行目における非ゼロ値が処理されるまでステップ２４０〜２４２を繰り返す。ｊがｊ_ｍａｘ以上である場合（ステップ２４０）、ワープにおける全てのスレッドは、必要に応じて同期され、ワープにおける全てのスレッドがステップ２４０〜２４２のループの実行を終了するのを可能とし、同期は、並列に動作するワープにおけるスレッドに起因して１つの実施形態では必要なく、同期から出ることがない一方で、さらなる実施形態においては同期が行われる。

ワープのサイズ、ワープにおけるスレッド数は、整数除算を使用して２で割られ、除算結果は、変数ｔ_ｍａｘを使用して表される（ステップ２４４）。ｔの場合、ワープにおけるスレッドのｉｄは、ｔ_ｍａｘ未満であり（ステップ２４５）、スレッドは、ともに組み合わせて、そのスレッドｔ及びスレッドｉｄがｔ＋ｔ_ｍａｘである他のスレッドによって上記実行されるステップ２４０〜２４２における計算から得られる合計の削減を行う（ステップ２４６）。ステップ２４６は、半分の合計数を削減し、組み合わせた合計は、スレッドｔ（組み合わせを行ったスレッド）についての合計として設定され、ワープにおけるスレッドｉｄがｔ＋ｔ_ｍａｘであるスレッドと前に関連付けられた合計は破棄される（ステップ２４６）。例えば、ワープにおいて３２のスレッドがあり且つスレッドについてのｔが０である場合、スレッドｔは、ステップ２４０〜２４２のスレッドｔ自身の性能から得られた合計と、スレッドｉｄが１６であるスレッドによってステップ２４０〜２４２の性能から得られた合計とを組み合わせ、組み合わせた合計は、ステップ２４６のその後の反復についてのスレッドｔの合計として設定され、１６のｉｄを有するスレッドは、もはや合計に関連付けられない。組み合わせに続いて、値ｔ_ｍａｘは、整数除算を使用して半分にカットされ、ｔ_ｍａｘとして設定される（ステップ２４７）。必要に応じて、合計の組み合わせを行ったスレッドは再度同期され（ステップ２４８）、ルーチン２３０は、上記ステップ２４５に戻る。ステップ２４５〜２４７のループの各後続反復中において、ループに参加しているスレッド数は、値ｔ_ｍａｘの削減に起因して半分に削減される。

ｔがｔ_ｍａｘ以上である場合（ステップ２４５）、スレッドがワープにおける第１のスレッドであるかどうか（ｔ＝０及びｒ＜ｒ_ｍａｘ）が判定される（ステップ２４９）。スレッドが第１のスレッドでない場合、スレッドの実行は終了する（ステップ２５１）。スレッドが第１のスレッドである場合、ステップ２４５〜２４７において記載された削減から生じる総合計は、密ベクトルｙに格納され（ステップ２５０）、そのスレッドの実行が終了する（ステップ２５１）。ルーチン２３０は、全ての起動ワープについての全てのスレッドの実行の終了によって終了する。ルーチン２３０はまた、以下の擬似コードを使用して表現することができる−擬似コードは、ＰＳＣＳＲ符号化に関連して記載されるが、ＰＣＳＲに関連する擬似コードを準用して記載することができる。
カーネル ｆ１Ｗ１Ｒ（ｙ，ｘ，Ａ_ｖｃ，Ａ_ｒ，Ａ_ｍ，ｒ_ｍｉｎ，ｒ_ｍａｘ） ／＊１Ｗ１Ｒ：１ワープ１行ＳｐＭＶカーネル＊／
共有 ｓｕｍ［ＷＡＲＰＳＩＺＥ］ ／＊ｓｕｍ：ワープにおけるスレッドによって共有されるローカル合計＊／
ｒ←ｒ_ｍｉｎ＋ｗａｒｐ＿ｉｄ（） ／＊このワープに割り当てられた頂点のランクの計算＊／
（ｒ＜ｒ_ｍａｘ）である場合 ／＊任意の正確に（ｒ_ｍａｘ−ｒ_ｍｉｎ）ワープが作成された場合＊／
ｉ←Ａ_ｍ［ｒ］ ／＊ｉ：ｒ番目にランク付けされた行のｉｄ＊／
ｔ←ｗａｒｐ＿ｔｈｒｅａｄ＿ｉｄ（） ／＊ｔ：ワープにおけるローカルスレッドｉｄ＊／
ｊ←Ａ_ｒ［ｉ］＋ｔ ／＊ｊ：このスレッドに割り当てられた第１の非ゼロエントリのインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのインデックス＊／
ｓｕｍ［ｔ］←０ ／＊ローカル総和積算器の初期化＊／
ｗｈｉｌｅ（ｊ＜ｊ_ｍａｘ） ／＊行の終わりに到達したかどうかのテスト＊／
ｓｕｍ［ｔ］←ｓｕｍ［ｔ］＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］ ／＊ｙｉ＝ΣｊＡｉ，ｊ×ｘｊの計算
ｊ←ｊ＋ＷＡＲＰＳＩＺＥ ／＊このスレッドについての次の非ゼロエントリに移動＊／
ｗｈｉｌｅ文終了
ｉｆ文終了
ｓｙｎｃ＿ｗａｒｐ＿ｔｈｒｅａｄｓ（） ／＊任意にワープにおけるスレッドが常に同期している場合＊／
ｔ_ｍａｘ←ＷＡＲＰＳＩＺＥ／２ ／＊ｔ_ｍａｘ：ローカル合計まで追加したスレッド数＊／
ｗｈｉｌｅ（ｔ＜ｔ_ｍａｘ） ／＊このスレッドが参加すべきであるかどうかのテスト＊／
ｓｕｍ［ｔ］←ｓｕｍ［ｔ］＋ｓｕｍ［ｔ＋ｔ_ｍａｘ］ ／＊２つのローカル合計を１つに削減＊／
ｔ_ｍａｘ←ｔ_ｍａｘ／２ ／＊合計追加スレッド数を半分にカット＊／
ｓｙｎｃ＿ｗａｒｐ＿ｔｈｒｅａｄｓ（） ／＊任意にワープにおけるスレッドが常に同期している場合＊／
ｗｈｉｌｅ文終了
（ｔ＝０及びｒ＜ｒ_ｍａｘ）である場合 ／＊これはワープにおける第１のスレッド？＊／
ｙ［ｉ］←ｓｕｍ［０］ ／＊ｙに総合計を格納＊／
ｉｆ文終了

ＳｐＭＴＶは、ブロックを処理するためにスレッドのワープを割り当てるｆ１Ｗ１Ｒカーネルを使用して同様に行うことができる。図２０Ａ〜図２０Ｂは、１つの実施形態にかかる図１５のルーチン１７０において使用するためのｆ１Ｗ１ＲカーネルによってＳｐＭＴＶを実行するためのルーチン２６０を示すフロー図である。ルーチン２６０は、ＰＳＣＣＲ又はＰＣＳＲ符号化のいずれかにおいてＳｐＭＶを行うために使用することができる。カーネルにおける全ての起動ワープについて反復処理ループが開始される（ステップ２６１）。起動ワープ及び起動ワープにおけるスレッドは、互いに並列に実行され、それゆえに、以下の説明がワープのいずれか又はスレッドのいずれかを参照する場合には、全ての他の起動ワープ又はスレッドは、スレッド又はワープのいずれかの実行が停止されるまで（以下のステップ２８１など）同時に並列にルーチン２６０の同じステップを実行する。

ローカル共有は、同じワープにおける全てのスレッド間で開始され、スレッドがステップ２７０〜２７２において後述する計算結果を共有するのを可能とする（ステップ２６２）。ワープのいずれかの固有のグローバルｉｄを取得するために関数ｗａｒｐ＿＿ｉｄ（）を使用し且つｒ_ｍｉｎの値に固有のグローバルｉｄを追加することによってワープのいずれかに割り当てられた列のランク（図１８を参照して上述したようにランクが割り当てられるパーティションの第１の列のランク）が取得され、加算結果は、ワープに割り当てられた列のランクであり、変数ｒを使用して表される（ステップ２６３）。値ｒは、Ａ’_ｏ［ｋ＋１］によって与えられる次のパーティションｋ＋１の第１の列のランクに等しい変数であるｒ_ｍａｘと比較することができ、起動されるワープ数がパーティションｋにおける列数（ｒ_ｍａｘ−ｒ_ｍｉｎ）に等しい場合には比較は任意である。ｒがｒ_ｍａｘ未満である場合（ステップ２６４）、ｒ番目のランク付けされた列のｉｄｊは、マッピング配列Ａ’_ｍにおいて識別される（ステップ２６５）。ｒがｒ_ｍａｘ以上である場合（ステップ２６４）、ルーチン２６０は、後述するステップ２８１に移動する。関数ｗａｒｐ＿ｔｈｒｅａｄ＿ｉｄ（）を使用することなどによってワープにおけるスレッドのいずれかのローカルｉｄ（ワープ内のｉｄ）が取得され、変数ｔを使用して表される（ステップ２６６）。そのスレッドに割り当てられた第１の非ゼロエントリのインデックスは、Ａ’_ｃ配列におけるｊ列目についてのエントリの値をｔに追加することによって取得され、加算結果は、変数ｉを使用して表される（ステップ２６７）。Ａ’_ｃ配列における次の（ｊ＋１）エントリが配置され、変数ｉ_ｍａｘとして設定される（ステップ２６８）。ｊ列目が行列における最後の列でない限り、Ａ’_ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒ（又はＡ’_ｖ）配列におけるインデックスであり、ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖｒ（又はＡ’_ｖ）配列におけるエントリの総数である。以下のステップにおいて記載される非ゼロ配列の値の乗算結果をともに加算する関数であるローカル総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ２６９）。ｊがｊ_ｍａｘ未満である場合（ステップ２７０）、ＳｐＭＴＶが行われた符号化がＰＳＣＳＣ又はＰＣＳＣであるかに応じた動作で、値ｉに対して乗算及び加算計算が実行される（ステップ２７１）。符号化がＰＳＣＳＣである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖｒ［ｉ］．ｖ×ｘ’［Ａ’_ｖｒ［ｉ］．ｒ］
ここで、インデックスｉを有するＡ’_ｖｒ配列の要素に格納された値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｖｒ要素に格納された行のインデックスであり、乗算結果は、ステップ２７１の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＣである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］
ここで、インデックスｉを有するＡ’_ｖ配列における値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｒ配列における数であり、乗算結果は、ステップ２７１の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、値ｉにワープサイズが加算され、加算結果がｉと設定され（ステップ２７２）、そのスレッドに割り当てられたその列における次のエントリに処理を移動する。ルーチン２６０は、上述したステップ２７０に戻り、ｊ列目における非ゼロ値の全てが処理されるまでステップ２７０〜２７２を繰り返す。ｉがｉ_ｍａｘ以上である場合（ステップ２７０）、ワープにおける全てのスレッドは、必要に応じて同期し、ワープにおける全てのスレッドがステップ２７０〜２７２のループの実行を終了するのを可能とし、同期は、並列に動作するワープにおけるスレッドに起因して１つの実施形態では必要なく、同期から出ることがない一方で、さらなる実施形態においては同期が行われる。

ワープのサイズ、ワープにおけるスレッド数は、整数除算を使用して２で割られ、除算結果は、変数ｔ_ｍａｘを使用して表される（ステップ２７４）。ｔの場合、ワープにおけるスレッドのｉｄは、ｔ_ｍａｘ未満であり、スレッドは、ともに組み合わせて、そのスレッドｔ及びスレッドｉｄがｔ＋ｔ_ｍａｘである他のスレッドによって上記実行されるステップ２７０〜２７２における計算から得られる合計の削減を行う（ステップ２７６）。ステップ２７６は、上述したステップ２４６と同様に、半分の合計数を削減し、組み合わせた合計は、スレッドｔ（組み合わせを行ったスレッド）についての合計として設定され、ワープにおけるスレッドｉｄがｔ＋ｔ_ｍａｘであるスレッドと前に関連付けられた合計は破棄される（ステップ２７６）。組み合わせに続いて、値ｔ_ｍａｘは、整数除算を使用して半分にカットされ、ｔ_ｍａｘとして設定される（ステップ２７７）。必要に応じて、合計の組み合わせを行ったスレッドは再度同期され（ステップ２７８）、ルーチン２６０は、上記ステップ２７５に戻る。ステップ２７５〜２７７のループの各後続反復中において、ループに参加しているスレッド数は、値ｔ_ｍａｘの削減に起因して半分に削減される。

ｔがｔ_ｍａｘ以上である場合、スレッドがワープにおける第１のスレッドであるかどうか（ｔ＝０及びｒ＜ｒ_ｍａｘ）が判定される（ステップ２７９）。スレッドが第１のスレッドでない場合、スレッドの実行は終了する（ステップ２８１）。スレッドが第１のスレッドである場合、ステップ２７５〜２７７において記載された削減から生じる総合計は、密ベクトルｙ’に格納され（ステップ２８０）、そのスレッドの実行が終了する（ステップ２８１）。ルーチン２６０は、全ての起動ワープについての全てのスレッドの実行の終了によって終了する。ルーチン２６０についての擬似コードは、準用する図１９Ａ〜図１９Ｂを参照して上記に示された擬似コードと同様である。

３２又は６４超の非ゼロエントリを有する行及び列の処理は、単一の行又は列を処理するためにスレッドのブロックを割り当てるｆ１Ｂ１Ｒカーネルを使用して最速に処理することができる。図２１Ａ〜図２１Ｂは、１つの実施形態にかかる図１５のルーチン１７０において使用するためのｆ１Ｂ１ＲカーネルによってＳｐＭＶを実行するためのルーチン２９０を示すフロー図である。ルーチン２９０は、ＰＳＣＳＲ又はＰＣＳＲ符号化のいずれかにおいてＳｐＭＶを行うために使用することができる。カーネルにおける全ての起動ブロックについて反復処理ループが開始される（ステップ２９１）。起動ブロック及び起動ブロックにおけるスレッドは、互いに並列に実行され、それゆえに、以下の説明がブロックのいずれか又はスレッドのいずれかを参照する場合には、全ての他の起動ブロック又はスレッドは、スレッド又はブロックのいずれかの実行が停止されるまで（以下のステップ３１１など）同時に並列にルーチン２９０のステップを実行する。

ローカル共有は、同じブロックにおける全てのスレッド間で開始され、スレッドがステップ３００〜３０２において後述する計算結果を共有するのを可能とする（ステップ２９２）。ブロックのいずれかの固有のグローバルｉｄを取得するために関数ｂｌｏｃｋ＿＿ｉｄ（）を使用し且つｒ_ｍｉｎの値に固有のグローバルｉｄを追加することによってブロックのいずれかに割り当てられた行のランク（図１７を参照して上述したようにランクが割り当てられるパーティションの第１の行のランク）が取得され、加算結果は、ブロックに割り当てられた行のランクであり、変数ｒを使用して表される（ステップ２９３）。値ｒは、Ａｏ［ｋ＋１］によって与えられる次のパーティションｋ＋１の第１の行のランクに等しい変数であるｒ_ｍａｘと比較することができ、起動されるブロック数がパーティションｋにおける行数（ｒ_ｍａｘ−ｒ_ｍｉｎ）に等しい場合には比較は任意である。ｒがｒ_ｍａｘ未満である場合（ステップ２９４）、ｒ番目のランク付けされた行のｉｄｉは、マッピング配列Ａ_ｍにおいて識別される（ステップ２９５）。ｒがｒ_ｍａｘ以上である場合（ステップ２９４）、ルーチン２９０は、後述するステップ３１１に移動し、スレッドの処理を停止する。関数ｂｌｏｃｋ＿ｔｈｒｅａｄ＿ｉｄ（）を使用することなどによってブロックにおけるスレッドのいずれかのローカルｉｄ（ブロック内のｉｄ）が取得され、変数ｔを使用して表される（ステップ２９６）。そのスレッドに割り当てられた第１の非ゼロエントリのインデックスは、Ａ_ｒ配列におけるｉ行目についてのエントリの値をｔに追加することによって取得され、加算結果は、変数ｊを使用して表される（ステップ２９７）。Ａ_ｒ配列における次の（ｉ＋１）エントリが配置され、変数ｊ_ｍａｘとして設定される（ステップ２９８）。以下のステップ３０１において記載される非ゼロ配列の値の乗算結果をともに加算する関数であるローカル総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ２９９）。ｊがｊ_ｍａｘ未満である場合（ステップ３００）、ＳｐＭＶが行われた符号化がＰＳＣＳＲ又はＰＣＳＲであるかに応じた動作で、値ｊに対して乗算及び加算計算が実行される（ステップ３０１）。符号化がＰＳＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］

計算において、インデックスｊを有するＡ_ｖｃ配列の要素に格納された値はｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｖｃ要素に格納された列のインデックスであり、乗算結果は、ステップ３０１の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａｖ［ｊ］×ｘ［ｖ［ｊ］］
ここで、インデックスｊを有するＡｖ配列における値は、ｘの要素によって乗算され、そのインデックスは、ｊ番目のインデックスを有するＡ_ｃ配列における数であり、乗算結果は、ステップ３０１の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｊの値にブロックサイズが加算され、加算結果がｊと設定され（ステップ３０２）、そのスレッドに割り当てられたその行における次のエントリに処理を移動する。ルーチン２９０は、上述したステップ３００に戻り、そのスレッドに割り当てられるｉ行目における非ゼロ値が処理されるまでステップ３００〜３０２を繰り返す。ｊがｊ_ｍａｘ以上である場合（ステップ３００）、ブロックにおける全てのスレッドは同期し、ブロックにおける全てのスレッドがステップ３００〜３０２のループの実行を終了するのを可能とする。

ブロックのサイズ、ブロックにおけるスレッド数は、整数除算を使用して２で割られ、除算結果は、変数ｔ_ｍａｘを使用して表される（ステップ３０４）。ｔの場合、ブロックにおけるスレッドのｉｄは、ｔ_ｍａｘ未満であり、スレッドは、ともに組み合わせて、そのスレッドｔ及びスレッドｉｄがｔ＋ｔ_ｍａｘである他のスレッドによって上記実行されるステップ３００〜３０２における計算から得られる合計の削減を行う（ステップ３０６）。ステップ３０６は、半分の合計数を削減し、組み合わせた合計は、スレッドｔ（組み合わせを行ったスレッド）についての合計として設定され、ブロックにおけるスレッドｉｄがｔ＋ｔ_ｍａｘであるスレッドと前に関連付けられた合計は破棄される（ステップ３０６）。組み合わせに続いて、値ｔ_ｍａｘは、整数除算を使用して半分にカットされ、ｔ_ｍａｘとして設定される（ステップ３０７）。合計の組み合わせを行ったスレッドは再度同期され（ステップ３０８）、ルーチン２９０は、上記ステップ３０５に戻る。ステップ３０５〜３０７のループの各後続反復中において、ループに参加しているスレッド数は、値ｔ_ｍａｘの削減に起因して半分に削減される。

ｔがｔ_ｍａｘ以上である場合、スレッドがブロックにおける第１のスレッドであるかどうか（ｔ＝０及びｒ＜ｒ_ｍａｘ）が判定される（ステップ３０９）。スレッドが第１のスレッドでない場合、スレッドの実行は終了する（ステップ３１１）。スレッドが第１のスレッドである場合、ステップ３０５〜３０７において記載された削減から生じる総合計は、密ベクトルｙに格納され（ステップ３１０）、そのスレッドの実行が終了する（ステップ３１１）。ルーチン２９０は、全てのスレッドの実行の終了によって終了する。ルーチン２９０はまた、以下の擬似コードを使用して表現することができる−擬似コードは、ＰＳＣＳＲ符号化に関連して記載されるが、ＰＣＳＲに関連する擬似コードを準用して記載することができる。
カーネル ｆ１Ｂ１Ｒ（ｙ，ｘ，Ａ_ｖｃ，Ａ_ｒ，Ａ_ｍ，ｒ_ｍｉｎ，ｒ_ｍａｘ） ／＊１Ｂ１Ｒ：１ブロック１行ＳｐＭＶカーネル＊／
共有 ｓｕｍ［ＢＬＯＣＫＳＩＺＥ］ ／＊ｓｕｍ：ブロックにおけるスレッドによって共有されるローカル合計＊／
ｒ←ｒ_ｍｉｎ＋ｂｌｏｃｋ＿ｉｄ（） ／＊このブロックに割り当てられた頂点のランクの計算＊／
（ｒ＜ｒ_ｍａｘ）である場合 ／＊任意の正確に（ｒ_ｍａｘ−ｒ_ｍｉｎ）ブロックが作成された場合＊／
ｉ←Ａ_ｍ［ｒ］ ／＊ｉ：ｒ番目にランク付けされた行のｉｄ＊／
ｔ←ｂｌｏｃｋ＿ｔｈｒｅａｄ＿ｉｄ（） ／＊ｔ：ブロックにおけるローカルスレッドｉｄ＊／
ｊ←Ａ_ｒ［ｉ］＋ｔ ／＊ｊ：このスレッドに割り当てられた第１の非ゼロエントリのインデックス＊／
ｊ_ｍａｘ←Ａ_ｒ［ｉ＋１］ ／＊ｊ_ｍａｘ：（ｉ＋１）行目における第１の非ゼロエントリのインデックス＊／
ｓｕｍ［ｔ］←０ ／＊ローカル総和積算器の初期化＊／
ｗｈｉｌｅ（ｊ＜ｊ_ｍａｘ） ／＊行の終わりに到達したかどうかのテスト＊／
ｓｕｍ［ｔ］←ｓｕｍ［ｔ］＋Ａ_ｖｃ［ｊ］．ｖ×ｘ［Ａ_ｖｃ［ｊ］．ｃ］ ／＊ｙｉ＝ΣｊＡｉ，ｊ×ｘｊの計算
ｊ←ｊ＋ＢＬＯＣＫＳＩＺＥ ／＊このスレッドについての次の非ゼロエントリに移動＊／
ｗｈｉｌｅ文終了
ｉｆ文終了
ｓｙｎｃ＿ｂｌｏｃｋ＿ｔｈｒｅａｄｓ（） ／＊ブロックにおける全てのスレッドが同期している＊／
ｔ_ｍａｘ←ＢＬＯＣＫＳＩＺＥ／２ ／＊ｔ_ｍａｘ：ローカル合計まで追加したスレッド数＊／
ｗｈｉｌｅ（ｔ＜ｔ_ｍａｘ） ／＊このスレッドが参加すべきであるかどうかのテスト＊／
ｓｕｍ［ｔ］←ｓｕｍ［ｔ］＋ｓｕｍ［ｔ＋ｔ_ｍａｘ］ ／＊２つのローカル合計を１つに削減＊／
ｔ_ｍａｘ←ｔ_ｍａｘ／２ ／＊合計追加スレッド数を半分にカット＊／
ｓｙｎｃ＿ｂｌｏｃｋ＿ｔｈｒｅａｄｓ（） ／＊ブロックにおける全てのスレッドが同期している＊／
ｗｈｉｌｅ文終了
（ｔ＝０及びｒ＜ｒ_ｍａｘ）である場合 ／＊これはブロックにおける第１のスレッド？＊／
ｙ［ｉ］←ｓｕｍ［０］ ／＊ｙに総合計を格納＊／
ｉｆ文終了

同様に、ＳｐＭＴＶは、１Ｂ１Ｒカーネルを使用して計算することができる。図２２Ａ〜図２２Ｂは、１つの実施形態にかかる図１５のルーチン１７０において使用するためのｆ１Ｂ１Ｒカーネル関数によってＳｐＭＴＶを実行するためのルーチン３２０を示すフロー図である。ルーチン３２０は、ＰＳＣＳＣ又はＰＣＳＣ符号化のいずれかにおいてＳｐＭＴＶを行うために使用することができる。カーネルにおける全ての起動ブロックについて反復処理ループが開始される（ステップ３２１）。起動ブロック及び起動ブロックにおけるスレッドは、互いに並列に実行され、それゆえに、以下の説明がブロックのいずれか又はスレッドのいずれかを参照する場合には、全ての他の起動ブロック及びスレッドは、スレッド又はブロックのいずれかの実行が停止されるまで（以下のステップ３４１など）同時に並列にルーチン３２０のステップを実行する。

ローカル共有は、同じブロックにおける全てのスレッド間で開始され、スレッドがステップ３３０〜３３２において後述する計算結果を共有するのを可能とする（ステップ３２２）。ブロックのいずれかの固有のグローバルｉｄを取得するために関数ｂｌｏｃｋ＿＿ｉｄ（）を使用し且つｒ_ｍｉｎの値に固有のグローバルｉｄを追加することによってブロックのいずれかに割り当てられた列のランク（図１８を参照して上述したようにランクが割り当てられるパーティションの第１の列のランク）が取得され、加算結果は、ブロックに割り当てられた列のランクであり、変数ｒを使用して表される（ステップ３２３）。値ｒは、Ａ’_ｏ［ｋ＋１］によって与えられる次のパーティションｋ＋１の第１の列のランクに等しい変数であるｒ_ｍａｘと比較することができ、起動されるブロック数がパーティションｋにおける列数（ｒ_ｍａｘ−ｒ_ｍｉｎ）に等しい場合には比較は任意である。ｒがｒ_ｍａｘ未満である場合（ステップ３２４）、ｒ番目のランク付けされた列のｉｄｊは、マッピング配列Ａ’_ｍにおいて識別される（ステップ３２５）。ｒがｒ_ｍａｘ以上である場合（ステップ３２４）、ルーチン３２０は、後述するステップ３４１に移動する。関数ｂｌｏｃｋ＿ｔｈｒｅａｄ＿ｉｄ（）を使用することなどによってブロックにおけるスレッドのいずれかのローカルｉｄ（ブロック内のｉｄ）が取得され、変数ｔを使用して表される（ステップ３２６）。そのスレッドに割り当てられた第１の非ゼロエントリのインデックスは、Ａ’_ｃ配列におけるｊ列目についてのエントリの値をｔに追加することによって取得され、加算結果は、変数ｉを使用して表される（ステップ３２７）。Ａ’_ｃ配列における次の（ｊ＋１）エントリが配置され、変数ｉ_ｍａｘとして設定される（ステップ３２８）。ｊ列目が行列における最後の列でない限り、Ａ’_ｃ配列における次のエントリは、（ｊ＋１）列目における第１の非ゼロエントリのＡ’_ｖｒ（又はＡ’_ｖ）配列におけるインデックスであり、ｊ列目が行列における最後の列である場合、Ａ’_ｃ配列における次のエントリは、Ａ’_ｖｒ（又はＡ’_ｖ）配列におけるエントリの総数である。以下のステップ３３１において記載される非ゼロ配列の値の乗算結果をともに加算する関数であるローカル総和計算部は、ゼロにおいて合計値を設定することによって初期化される（ステップ３２９）。ｉがｉ_ｍａｘ未満である場合（ステップ３３０）、ＳｐＭＴＶが行われた符号化がＰＳＣＳＣ又はＰＣＳＣであるかに応じた動作で、値ｉに対して乗算及び加算計算が実行される（ステップ３３１）。符号化がＰＳＣＳＣである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖｒ［ｉ］．ｖ×ｘ’［Ａ’_ｖｒ［ｉ］．ｒ］
ここで、インデックスｉを有するＡ’_ｖｒ配列の要素に格納された値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｖｒ要素に格納された行のインデックスであり、乗算結果は、ステップ３３１の前の反復中に実行された乗算結果の合計に加算される。あるいは、符号化がＰＣＳＲである場合、計算は、以下の式にしたがって行われる。
ｓｕｍ←ｓｕｍ＋Ａ’_ｖ［ｉ］×ｘ’［Ａ’_ｒ［ｉ］］
ここで、インデックスｉを有するＡ’_ｖ配列における値は、ｘ’の要素によって乗算され、そのインデックスは、ｉ番目のインデックスを有するＡ’_ｒ配列における数であり、乗算結果は、ステップ３３１の前の反復中に実行された乗算結果の合計に加算される。いずれかの式に基づいて計算が終了すると、ｉの値にブロックサイズが加算され、加算結果がｉと設定され（ステップ３３２）、その列における次のエントリに処理を移動する。ルーチン３２０は、上述したステップ３３０に戻り、ｊ列目における非ゼロ値が処理されるまでステップ３３０〜３３２を繰り返す。ｉがｉ_ｍａｘ以上である場合（ステップ３３０）、ブロックにおける全てのスレッドは同期し、ブロックにおける全てのスレッドがステップ３３０〜３３２のループの実行を終了するのを可能とする。

ブロックのサイズ、ブロックにおけるスレッド数は、整数除算を使用して２で割られ、除算結果は、変数ｔ_ｍａｘを使用して表される（ステップ３３４）。ｔの場合、ブロックにおけるスレッドのｉｄは、ｔ_ｍａｘ未満であり、スレッドは、ともに組み合わせて、そのスレッドｔ及びスレッドｉｄがｔ＋ｔ_ｍａｘである他のスレッドによって上記実行されるステップ３３０〜３３２における計算から得られる合計の削減を行う（ステップ３３６）。ステップ３３６は、半分の合計数を削減し、組み合わせた合計は、スレッドｔ（組み合わせを行ったスレッド）についての合計として設定され、上記のステップ２４６と同様に、ブロックにおけるスレッドｉｄがｔ＋ｔ_ｍａｘであるスレッドと前に関連付けられた合計は破棄される（ステップ３３６）。組み合わせに続いて、値ｔ_ｍａｘは、整数除算を使用して半分にカットされ、ｔ_ｍａｘとして設定される（ステップ３３７）。合計の組み合わせを行ったスレッドは再度同期され（ステップ３３８）、ルーチン３２０は、上記ステップ３３５に戻る。ステップ３３５〜３３７のループの各後続反復中において、ループに参加しているスレッド数は、値ｔ_ｍａｘの削減に起因して半分に削減される。

ｔがｔ_ｍａｘ以上である場合、スレッドがブロックにおける第１のスレッドであるかどうか（ｔ＝０及びｒ＜ｒ_ｍａｘ）が判定される（ステップ３３９）。スレッドが第１のスレッドでない場合、スレッドの実行は終了する（ステップ３４１）。スレッドが第１のスレッドである場合、ステップ３３５〜３３７において記載された削減から生じる総合計は、密ベクトルｙ’に格納され（ステップ３４０）、そのスレッドの実行が終了する（ステップ３４１）。ルーチン３２０は、全てのスレッドの実行の終了によって終了する。ルーチン３２０についての擬似コードは、準用する図２１Ａ〜図２１Ｂを参照して上記に示された擬似コードと同様である。

図３及び図４を参照して上述したように、ＳｐＭＶ及びＳｐＭＴＶの結果の一般的な用途の１つは、べき乗法を使用して実行されることができるページランク（登録商標）アルゴリズムなどのランク付けアルゴリズムである。図２３は、１つの実施形態にかかる図６及び図９の方法６０及び９０において使用するためのべき乗法を実行するためのルーチン３５０を示すフロー図である。ルーチン３５０は、図６及び図９の方法において使用されるルーチンを参照して上述したＳｐＭＶ及びＳｐＭＴＶの全てのバリエーションの結果に適用するために使用することができる。ｄ∈（０，１）を減衰係数とし、ｎをウェブページ数とする。Ｐをｎ×ｎの正方行列とする。

ページｊからページｉへのリンクがある場合
それ以外
ここで、Ｌｊは、ページｊから出るリンクの数である。ｘ及びｙをサイズｎの２つの密ベクトルとし、ε＊を停止閾値とする。初期のページランクの確率分布がベクトルｘについて設定される（ステップ３５１）。εの値がε＊の値未満である場合に反復処理ループ（ステップ３５２〜３５６）が実行される（ステップ３５２）。値ｙは、以下の式に基づいて設定される。

ここで、１は、（ｎ×１）列ベクトルである（ステップ３５３）。値εは、以下の式にしたがって決定される。

ここで、値εは、ｙ及びｘの差異に等しい（ステップ５３）。最後に、ｘは、値ｙに等しくなるように設定される。ルーチン３５０は、ループの次の反復に移動し（ステップ３５６）、εが値ε＊未満であるまでループ（３５２〜３５６）を介した処理を継続し、その後にルーチン３５０は終了する。ルーチン３５０はまた、以下の擬似コードを使用して表現することができる。
ｘ←ｘ０ ／＊初期のページランクの確率分布＊／
ループ

ε←｜ｙ−ｘ｜
ｘ←ｙ
ε＜ε＊まで

上述したシステム３０及び方法６０、９０は、疎行列、すなわち、ＳＣＳＲ、ＳＣＳＣ、ＰＣＳＲ、ＰＣＳＣ、ＰＳＣＳＲ及びＰＳＣＳＣについての６つの新たな符号化の利点を活用し、従来技術よりも優れている。６つのうち、ＳＣＳＲ及びＳＣＳＣは、それらの構造化されていない相手方のＣＳＲ及びＣＳＣと全く同じ空間の複雑性を有する。一方、残りの４つの新たな符号化は、僅かに高い空間的要件を有する。ＰＣＳＲ及びＰＳＣＳＲについて、過剰な空間的オーバーヘッドは、非ゼロ行（すなわち、少なくとも１つの非ゼロエントリを有する行）の数において線形である。ＰＣＳＣ及びＰＳＣＳＣについて、オーバーヘッドは、非ゼロ列（すなわち、少なくとも１つの非ゼロエントリを有する列）の数において線形である。ＰＳＣＳＲ符号化の有効性を評価するために、フロリダ大学の疎行列コレクションにおいてみられる様々な疎行列に対してＳｐＭＶを実行する実験を行った。行列のパーティションを作成するために、実験における全てのデータセットについてＡｓ＝［１，３２，１０２４，∞］が使用された。その結果、ＰＳＣＳＲの第１のパーティションは、少なくとも１つの非ゼロを有し且つ３２未満の非ゼロエントリを有する行を含み、第２のパーティションは、３２以上１０２４未満の非ゼロエントリを有する行を含み、第３のパーティションは、１０２４以上の非ゼロエントリを有する行を含んでいた。

表１は、（「＃行」とラベル付けされている）行数、非ゼロ行の数（「＃非ゼロ行」）、行列の非ゼロエントリの総数（＃非ゼロ）、行あたりの非ゼロエントリの平均数（「平均＃非ゼロ／行」）、（ＳＣＳＲと同じである）ＣＳＲ符号化のサイズ、ＰＳＣＳＲ符号化のサイズ、及び、元のＣＳＲ符号化の割合としてのＰＳＣＳＲの過剰な空間的オーバーヘッドを含む、実験において使用された疎行列の統計を示している。観察されるように、ＰＳＣＳＲの過剰な空間的オーバーヘッドは、行の非ゼロエントリの平均数が増加するのにともない低減する。オーバーヘッドは、行あたり２．１個の非ゼロエントリのみを有する行列「ｗｉｋｉ−Ｔａｌｋ」について最高である（１９．０６％）一方で、オーバーヘッドは、行あたり２２．３個の非ゼロエントリの平均を有する行列「ｅｕ−２００５」について２．１９％まで低下する。表１における全ての行列にわたって平均化したＰＳＣＳＲの過剰な空間的オーバーヘッドは５．４５％である。

λを通常は実装依存定数であるＰＳＣＳＲ（又はＰＳＣＳＣ）のＡ_ｒ（又はＡ’_ｃ）における要素のサイズに対するＡ_ｖｃ（又はＡ’_ｖｃ）における要素のサイズの比率とする。ｂをＰＳＣＳＲ（又はＰＳＣＳＣ）における行（又は列）の数に対する非ゼロエントリの数の比率とし、γをＰＳＣＳＲ（又はＰＳＣＳＣ）における行（又は列）の総数に対する非ゼロ行（又は列）の数の比率とする。δを通常のＣＳＲ（又はＣＳＣ）に対するＰＳＣＳＲ／ＰＣＳＲ（又はＰＳＣＳＣ／ＰＣＳＣ）の過剰な空間的オーバーヘッドとする。そのδは、

によって与えられて示すことができる。実施形態の１つにおいて、Ａ_ｖｃにおける要素は、Ａ_ｒにおける要素の２倍の大きさであり、それゆえにλ＝２である。行列「ｗｉｋｉ−Ｔａｌｋ」について、γ＝２，３６９，１８１／２，３９４，３８５＝０．９８９５、ｂ＝５，０２１，４１０／２，３９４，３８５＝２．０９７を有する。それゆえに、過剰な空間的オーバーヘッドは、δ＝γ（λｂ＋１）−１＝０．９８９５×（２×２．０９７＋１）−１＝１９．０５％であり、実験において観察された実際のオーバーヘッドに近い。行列「ｅｕ−２００５」について、γ＝８６２，６６４／８６２，６６４＝１、ｂ＝１９，２３５，１４０／８６２，６６４＝２２．３０である。それゆえに、δ＝γ（λｂ＋１）−１＝１×（２×２２．３０＋１）−１＝２．１９％である（経験的な数と同じである）。δ＝γ（λｂ＋１）−１という解析式は、表１に記録されたＰＳＣＳＲの実際の過剰なオーバーヘッドと非常によく一致する（大抵の場合、検出可能なエラーなし）ことを確認することができる。この式は、（１＋δ）の係数を乗じた通常のＣＳＲ（又はＣＳＣ）符号化のサイズとしてＰＳＣＳＲ／ＰＣＳＲ（又はＰＳＣＳＣ／ＰＣＳＣ）符号化のサイズを正確に予測することを可能とする。予測された過剰な空間的オーバーヘッドは、ＰＳＣＳＲ／ＰＣＳＲ符号化を形成することが、特定のハードウェアセットアップについて有用であるかどうか、又は、不十分なハードウェアリソースがそのような符号化の形成又は使用を不可能とするかどうかを判定するために使用することができる。

図１３を参照して始めて記載されたＧＰＵベースのＳｐＭＶルーチンの性能が表１の全ての行列についてテストされた。比較のため、ＣＰＵベースのＳｐＭＶ実装の結果が表３に含まれる。使用されたテストマシンは、それぞれが３．４６ＧＨｚで動作する６個のコアを有する２つのＩｎｔｅｌ Ｘｅｏｎ Ｘ５６９０プロセッサを有する。単一のＸ５６９０プロセッサのキャッシュサイズは１２ＭＢであり、同じボックス内のＮｖｉｄｉａ（登録商標）のＧＴＸ５８０ ＧＰＵのＬ２キャッシュサイズよりもかなり大きい。単一のＣＰＵコアに対するＳｐＭＶにおけるＧＰＵの高速化を測定するために、ＣＰＵ実装は、単一スレッドを使用する。

ＣＰＵベースの実装についての語は、実装が大規模グラフにおける最先端のＳｐＭＶにあるということである。テストマシンの単一コアを使用して、ＣＰＵは、Ｋｗａｋら、「Ｗｈａｔ ｉｓ Ｔｗｉｔｔｅｒ， ａ Ｓｏｃｉａｌ Ｎｅｔｗｏｒｋ ｏｒ ａ Ｎｅｗｓ Ｍｅｄｉａ？」、Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ｔｈｅ １９ｔｈ ｉｎｔｅｒｎａｔｉｏｎａｌ ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ｗｏｒｌｄ Ｗｉｄｅ Ｗｅｂ、２０１０年、５９１〜６００ページによって導入された、ある４１７０万人のユーザ（すなわち、行）及び１４．７億の接続（すなわち、非ゼロ）のＴｗｉｔｔｅｒ（登録商標）フォロワーネットワークにおいて３６秒でＳｐＭＶの反復を実行する。同じＴｗｉｔｔｅｒ（登録商標）フォロワーネットワークは、ベンチマークデータセットとしての様々なＳｐＭＶ実装によって使用されている。表２は、文献において見出された最良のＣＰＵベースの大規模ＳｐＭＶシステムのいくつかの実行時性能を比較する。最初の４つのシステムについての番号は、Ｇｏｎｚａｌｅｚら、ＰｏｗｅｒＧｒａｐｈ： Ｄｉｓｔｒｉｂｕｔｅｄ Ｇｒａｐｈ−Ｐａｒａｌｌｅｌ Ｃｏｍｐｕｔａｔｉｏｎ ｏｎ Ｎａｔｕｒａｌ Ｇｒａｐｈｓ、ＯＳＤＩ、２０１２年、第１２巻、第１号、２ページによって報告され、そのＰｏｗｅｒＧｒａｐｈシステムは、それぞれが１０ギガイーサネット（登録商標）を介して接続された２３ＧＢのＲＡＭを有する２つのクワッドコアのＩｎｔｅｌ Ｘｅｏｎ Ｘ５５７０プロセッサを有するＡｍａｚｏｎ ＥＣ２ ｃｃ１．４ｘｌａｒｇｅ Ｌｉｎｕｘ（登録商標）インスタンスの６４のノードクラスタにおいて３．６秒でＳｐＭＶの反復を実行することができる。ＨｉｐｅｒＧｒａｐｈと称される単一ノードシステムはまた、全１２個のコアが使用された場合に３．６秒でＳｐＭＶの反復を終了することができる。コアあたりの高速化を測定するために、ＨｉｐｅｒＧｒａｐｈは、同じＴｗｉｔｔｅｒ（登録商標）のフォロワーネットワーク上の完全なパスを終了するために３６秒かかった単一スレッドモードで実行された。換言すれば、ＰｏｗｅｒＧｒａｐｈは、３．４６ＧＨｚ＠の単一コアでＨｉｐｅｒＧｒａｐｈに対して１０倍の高速化を得るために２．９３ＧＨｚ＠の５１２個のコアを使用する。ＳｐＭＶに対するＨｉｐｅｒＧｒａｐｈの単一コアの性能は競合的であり、それゆえに、後述する実験についてのベースラインとして使用された。

表３は、ＣＰＵ（すなわち、１コアを有するＨｉｐｅｒＧｒａｐｈ）及び様々な疎行列に対するＧＰＵベースのＳｐＭＶの性能を示している。ＧＰＵベースのＳｐＭＶの実行は、以下の３つの部分にわけられる：（ａ）ＰＳＣＳＲ符号化によって必要とされるパーティションを生成するパーティション、（ｂ）ＣＰＵからＣＰＵに対してパーティション化された行列を転送する負荷、及び、（ｃ）ＧＰＵにおいてＳｐＭＶの単一反復を実行するＳｐＭＶ。（ａ）及び（ｂ）の双方の部分は、そのコストがＳｐＭＶの複数の反復にわたって償却されることができるワンタイム動作である。表３の最後の列は、シーケンシャルＣＰＵベースのＳｐＭＶに対するＧＰＵベースの高速化を示しており、１００回の反復にわたって積算される。行列が既にＰＳＣＳＲフォーマットで符号化されている場合、パーティション時間は、表３において単にゼロとしなければならないことに留意されたい。高速化は、２３．５×から３８．３×まで及び、平均は３１倍速である。

実験は、さらに、上述した方法を実装するＧＰＵベースのハードウェアが上述した従来技術を越えて実現することができるという符号化のさらなる利点を示す。大幅に高速化するのみならず、ＧＰＵはまた、ＣＰＵよりもはるかに安価である：ＧＴＸ ５８０のＧＰＵのコストは、現在＄４００未満であるが、単一のＸｅｏｎ Ｘ５６９０プロセッサは、約＄１８００のコストがかかる。Ｘ５６９０の６個全てのコアが使用された場合に我々のシーケンシャルなＣＰＵ ＳｐＭＶの線形高速化を仮定すると、ＧＰＵは、Ｘ５６９０プロセッサにおける完全に並列化されたＣＰＵの実装よりも３１／６＝５．２×高速化するであろう。ドルあたりの性能の観点から、ＧＰＵベースのＳｐＭＶは、上記使用されたテストマシンについて約２３×以上の費用対効果がある。

ＧＰＵベースの解決策はまた、既存のシステムのＧＰＵの追加又は更新が同じボックス内にＣＰＵを追加又は更新するよりもはるかに容易であることから、より拡張性がある。コモディティハードウェアについて、ＣＰＵの最大数は、通常は２又は４である。一方、単一のコモディティＧＰＵサーバは、８個のＧＰＵまで保持することができ、それぞれが数千コアまで有することができる（例えば、Ｎｖｉｄｉａ（登録商標）のＴｅｓｌａ Ｋ４０ ＧＰＵは２８８０のコアを有する）。それゆえに、ドルあたりの性能又はワットあたりのＦＬＯＰＳのみならず、データセンタのために特に重要である１Ｕラック空間あたりの性能においても、ＧＰＵはＣＰＵに優れている。

本発明は、その実施形態を参照して具体的に示されて説明されたが、当業者は、形態及び詳細における上述した及び他の変更が本発明の精神及び範囲から逸脱することなく行われ得ることを理解するであろう。

Claims (10)

1. 構造化疎行列表現を取得するためのコンピュータ実装方法において、
   各部分が行及び列のうちの１つを含む行列の部分における１つ以上の順序で配置された１つ以上の非ゼロエントリを含む前記行列の構造化圧縮表現を取得することと、
   前記表現に含まれる複合配列であり、各要素が非ゼロエントリのいずれかとその非ゼロエントリを含む前記部分のいずれかのインデックスとを含む１つ以上の要素を含む前記複合配列を取得することと、
   前記順序のうちの１つ以上における第１である前記非ゼロエントリを含む各要素の前記複合配列におけるインデックスを含むインデックス配列であり、さらに前記行列における複数の前記非ゼロエントリを含む前記インデックス配列を取得することとを備える、方法。
2. 各複合配列の要素インデックスが、列優先順序でその非ゼロエントリの行のインデックスを含み、前記順序が各列における前記非ゼロエントリの順序を含む、請求項１に記載の方法。
3. さらに、
   各列を処理することによって複数の要素を含む密ベクトルによって前記行列の疎行列転置ベクトル乗算を実行することと、
   前記インデックス配列を使用してその列における前記非ゼロエントリを含む全ての前記複合配列要素からなる前記複合配列におけるインデックスを識別することと、
   そのエントリの前記インデックスを使用して前記密ベクトルの要素のいずれかによってその列における前記非ゼロエントリのそれぞれを乗算することと、
   その列における前記非ゼロエントリのそれぞれについての乗算結果を加算し、異なる密ベクトルの加算結果を格納することとを備える、請求項２に記載の方法。
4. 各複合配列の要素インデックスが、行優先順序でその非ゼロエントリの列のインデックスを含み、前記順序が各行における前記非ゼロエントリの順序を含む、請求項１に記載の方法。
5. さらに、
   各行を処理することによって複数の要素を含む密ベクトルによって前記行列の疎行列ベクトルの乗算を実行することと、
   前記インデックス配列を使用してその行における前記非ゼロエントリを含む全ての前記複合配列要素からなる前記複合配列におけるインデックスを識別することと、
   そのエントリの前記インデックスを使用して前記密ベクトルの要素のいずれかによってその行における前記非ゼロエントリのそれぞれを乗算することと、
   その行における前記非ゼロエントリのそれぞれについての乗算結果を加算し、異なる密ベクトルの加算結果を格納することとを備える、請求項４に記載の方法。
6. さらに、
   前記構造化圧縮表現を処理することと、
   前記処理の結果を使用してランク付け分析を実行することとを備える、請求項１に記載の方法。
7. 効率的な疎行列表現のためのコンピュータ実装方法において、
   １つ以上の非ゼロエントリを含む疎行列の圧縮表現を取得することであり、前記行列の行及び列のいずれかを含む前記行列の部分が前記行列におけるそれらの位置に基づいてインデックス付けされることと、
   前記行列の部分について複数のパーティションを定義することと、
   前記圧縮表現を使用して前記部分のそれぞれにおける複数の前記非ゼロエントリを取得することと、
   その部分における複数の前記非ゼロエントリに基づいて前記パーティションのいずれかを有する部分のそれぞれを関連付けることと、
   前記パーティションのそれぞれに関連する全ての部分であり、それらのインデックスの順序でリスト化された前記部分のリストを作成することと、
   前記リストを含むマッピング配列を含む前記行列のパーティション化された圧縮表現を作成することとを備える、方法。
8. さらに、
   前記マッピング配列を処理することと、
   前記処理の結果に基づいてランク付け分析を実行することとを備える、請求項７に記載の方法。
9. さらに、
   前記パーティションのそれぞれに関連する前記行列の部分における非ゼロエントリの数に基づいて前記パーティションをインデックス付けることと、
   前記パーティションのそれぞれについて、前記パーティションのインデックスに基づいて前記配列におけるそのパーティションに先行する全てのパーティションのサイズであり、前記パーティションに関連する複数の部分の数を含む前記サイズを識別し、そのパーティションについての前記リストを前記サイズに基づく位置において前記配列に挿入することとを備える、請求項７に記載の方法。
10. 前記行列の行を含む前記各部分が、さらに、
    前記マッピング配列にリスト化されたパーティションのそれぞれの行のそれぞれをシーケンシャルに処理することによって複数の要素を含む密ベクトルによって前記行列の疎行列ベクトル乗算を実行することと、
    その行における前記非ゼロエントリのそれぞれに関連付けられた前記圧縮表現におけるインデックスを識別することと、
    その非ゼロエントリの前記インデックスを使用して前記密ベクトルの要素のいずれかによってその行における前記非ゼロエントリのそれぞれを乗算することと、
    その行における前記非ゼロエントリのそれぞれについての乗算結果を加算し、異なる密ベクトルにおける加算結果を格納することとを備える、請求項７に記載の方法。
    

Images