JP5782265B2

JP5782265B2 - 行列計算処理方法、プログラム及びシステム

Info

Publication number: JP5782265B2
Application number: JP2011022311A
Authority: JP
Inventors: 弘揮 ▲柳▼澤
Original assignee: International Business Machines Corp
Current assignee: International Business Machines Corp
Priority date: 2011-02-04
Filing date: 2011-02-04
Publication date: 2015-09-24
Anticipated expiration: 2031-02-04
Also published as: JP2012164039A; US9098460B2; US20120203815A1; US20120317160A1

Description

この発明は、コンピュータの処理によって、行列の計算を行う技法に関し、より詳しくは、FMM(Funny Matrix Multiplication)の計算に関する。

FMMとは、通常の行列乗算処理におけるadd（加算）とmultiply（乗算）を、それぞれ、min演算とadd演算に置き換えた行列計算処理のことである。すなわち、２つの入力行列A、Bについて、C[i,j] = min_k{A[i,k] + B[k,j]}を全ての(i,j)について計算することがFMMの定義である。FMMは正方行列に限定されないが、ここでは便宜上、行列Aと行列Bは、n×nの正方行列であると仮定する。

FMMは、最短経路計算や画像処理などへの応用がある。例えば、A. V. Aho, J. E. Hopcroft, and J. D. Ullman, The Design and Analysis of Computer Algorithms. Addison-Wesley, 1974には、FMMを使った最短経路計算の技法の例が記述されている。

FMMにおいては、長さnの２つのベクトルva、vbについて、va[k] + vb[k]の価を、k = 1,..,nの順に計算していき、そのうちの最小値をとる計算が基本となる。つまり、min_k{va[k] + vb[k]}を計算する処理が中心となり、この処理をn²回行うので、
FMM全体の処理にO(n³)時間かかる。

ところで、近年、交通シミュレーションやカー・ナビゲーションなどの分野で、最短経路計算をより高速に行うことの要望が高まっており、そのためFMMを高速化するアルゴリズムが、次の文献に記述されている。

J. J. McAuley and T. S. Caetano, “An expected-case sub-cubic solution to the all-pairs shortest path problem in R,” arXiv:0912.0975v1, 2009は、二つのベクトルvaとvbについて、min_k{va[k] + vb[k]}を計算する前に、va[a1] ≦ va[a2] ≦…≦ va[an], vb[b1] ≦ vb[b2] ≦…≦ vb[bn] となるような順列a1,a2,…,anとb1,b2,…,bnを事前に計算することを開示する。FMM全体の処理では、この事前計算をするには、全体で2n回のソート処理を行えば十分であり、すると、全体でO(n²log n)時間で計算を行うことが可能である。

J. J. McAuley and T. S. Caetano, "An expected-case sub-cubic solution to the all-pairs shortest path problem in R," arXiv:0912.0975v1, 2009

この発明の目的は、FMMの計算において、上記非特許文献１に開示されている技法を更に改良して高速にすることにある。

本発明においても、min_k{va[k] + vb[k]}を計算する前に、va[a1] ≦ va[a2] ≦…≦ va[an], vb[b1] ≦ vb[b2] ≦…≦ vb[bn] となるような順列a1,a2,…,anとb1,b2,…,bnを事前に計算することは行われる。

本発明の特徴は、FMMで中心となるva[k] + vb[k]の最小値を計算する処理において、k=1,…nについて順番に処理するのではなく、best = ∞に初期設定してから、以下の処理Xと処理Yを一回ずつ適用した上で、bestの値をva[k]+vb[k]の最小値として出力することにある。
(処理X) k=a1,a2,…の順にva[k]+vb[k]の値を計算していき、それまでに見つかった最小値をbestとしたときに、va[k] > best/2となるkで処理をやめる（そのようなkが無ければk = anまで処理したらやめる）
(処理Y) 処理Xと同様の処理をk = b1,b2,…についても行い、vb[k] > best/2の値となるkで処理をやめる（そのようなkが無ければk = bnまで処理したらやめる）

このような処理を採用したことにより、FMMの計算において、SIMD命令を利用して処理を高速化することが可能となった。この場合、行列Aと行列BのFMMを計算するとすると、行列Aを、列が主のレイアウト(column-majorlayout)、行列Bを、行が主のレイアウト(row-majorlayout)で保存するのがポイントである。ここで、列が主のレイアウトとは、行列 A をメモリ上に保管する際に（行列Aで）列方向に隣り合う要素を（メモリ上でも）なるべく隣り合うように並べる方法のことをいう。つまり、行列 A を、A[1,1], A[2,1], A[3,1],..., A[n,1], A[1,2], A[2,2], A[3,2], ..., A[n,2], A[1,3], A[2,3], ........., A[n,n] といった順番で並べる。一方、行が主のレイアウトとは、列が主のレイアウトにおいて、列と行を入れ替えたものであり、つまり、行列 A を、A[1,1], A[1,2], A[1,3],..., A[1,n], A[2,1], A[2,2], A[2,3], ..., A[2,n], A[3,1], A[3,2], ........., A[n,n] といった順番で並べる。

この発明によれば、bestの値を一旦計算して、その値に基づき、最小値の計算を打ち切るようにすることにより、最小値の計算を速くすることによって、FMMの計算処理を高速化することができる。

また、非特許文献１の技法ではSIMD命令を利用することは困難であったが、本発明の技法では、SIMD命令を有利に適用して、更にFMMの計算処理を高速化することが可能である。

本発明を実施するためのハードウェア構成のブロック図である。本発明に係る機能論理ブロック図である。本発明の一実施例の処理全体のフローチャートを示す図である。図３における、行列を更新する処理のフローチャートを示す図である。図３における、行列を更新する処理のフローチャートを示す図である。行列におけるSIMD命令に対応した処理を説明するための図である。本発明の、SIMD命令に対応した実施例の処理全体のフローチャートを示す図である。図７における、行列を更新する処理のフローチャートを示す図である。図７における、行列を更新する処理のフローチャートを示す図である。

以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。

図１を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータ・ハードウェアのブロック図が示されている。図１において、システム・パス１０２には、ＣＰＵ１０４と、主記憶（ＲＡＭ）１０６と、ハードディスク・ドライブ（ＨＤＤ）１０８と、キーボード１１０と、マウス１１２と、ディスプレイ１１４が接続されている。ＣＰＵ１０４は、好適には、３２ビットまたは６４ビットのアーキテクチャに基づくものであり、例えば、インテル社のPentium（商標）４、Core(商標)2 Duo、Xeon(商標)、AMD社のAthlon（商標）などを使用することができる。この実施例の目的のため、ＣＰＵ１０４は、ＳＩＭＤ(Single Instruction Multiple Data)命令をもつものである。主記憶１０６は、好適には、４ＧＢ以上の容量をもつものである。ハードディスク・ドライブ１０８は、行列要素の基となる大量のデータを格納できるように、例えば、５００ＧＢ以上の容量をもつものであることが望ましい。

ハードディスク・ドライブ１０８には、個々に図示しないが、オペレーティング・システムが、予め格納されている。オペレーティング・システムは、Linux(商標)、マイクロソフト社のWindows XP(商標)、Windows(商標)7、アップルコンピュータのMac OS（商標）などの、ＣＰＵ１０４に適合する任意のものでよい。

ハードディスク・ドライブ１０８にはまた、好適には、C、C++、C#、Java(R)などのプログラム言語処理系も格納されている。このプログラム言語処理系は、後で説明する、本発明に係るFMM(Funny Matrix Multiplication)の計算処理のためのモジュールを作成し、保守するために使用される。

ハードディスク・ドライブ１０８はさらに、プログラム言語処理系でコンパイルするためのソースコードを書くためのテキスト・エディタ、及び、Eclipse（商標）などの開発環境を含んでいてもよい。

キーボード１１０及びマウス１１２は、オペレーティング・システムまたは、ハードディスク・ドライブ１０８から主記憶１０６にロードされ、ディスプレイ１１４に表示されたプログラム（図示しない）を起動したり、文字を打ち込んだりするために使用される。

ディスプレイ１１４は、好適には液晶ディスプレイであり、例えば、ＸＧＡ（１０２４×７６８の解像度）、またはＵＸＧＡ（１６００×１２００の解像度）などの任意の解像度のものを使用することができる。ディスプレイ１１４は、図示しないが、本発明の処理を開始するための操作ウインドウや、FMMの計算結果を表示するために使用される。

次に、図２の機能ブロック図を参照して、本発明の処理を実行するための処理ルーチンについて説明する。これらの処理ルーチンはC、C++などの、好適にはSIMD命令を利用可能なプログラム言語で作成されて、実行可能な形式でハードディスク・ドライブ１０８に保存され、オペレーティング・システムの動作で主記憶１０６にロードされて実行される。

メイン・ルーチン２０２は、本発明の全体の動作を統合するためのプログラムであり、図示しないが、ディスプレイ１１４に操作ウインドウを表示したり、ユーザの操作を受け付けて処理を開始したりする機能を有する。

入力ルーチン２０４は、ハードディスク・ドライブ１０８に保存された地図データなどである処理データ２０６のファイルからデータを読み込んで、行列Aと行列Bの各要素を決定する機能をもつ。

ソート・ルーチン２０８は、本発明に従いFMMの計算を行うためのインデックスの順列を求めるためのソートを行う機能をもつ。この際のソートのアルゴリズムは、クイック・ソートであるが、これには限定されず、シェル・ソート、ヒープ・ソート、マージ・ソートなどの任意の適当なソート・アルゴリズムを用いることができる。

更新ルーチン２１０は、ソートされたインデックスの順列を用いて、FMMの結果の行列Cの要素を更新するための処理を実行する機能をもつ。更新ルーチン２１０の処理の詳細は、図３〜図５のフローチャート、あるいは図７〜図９のフローチャートを参照して、後で説明する。

出力ルーチン２１２は、行列Aと行列BとのFMM計算の結果得られた行列Cを、結果データ２１４のファイルとして書き出す機能をもつ。

次に、図３のフローチャートを参照して、本発明による全体のFMM計算処理について説明する。この処理を開始するにあたって、メイン・ルーチン２０２は予め入力ルーチン２０４を呼び出し、処理データ２０６のファイルからデータを読み込むことによって、行列Aと行列Bの各要素の値を決定する。ここでは、行列Aと行列Bは、n×nの正方行列であると仮定する。しかし、本発明の処理は、正方行列以外の行列にも適用できることを理解されたい。

さて図３において、ステップ３０２からステップ３１４までは、iの1からnまでの繰り返しのループである。ステップ３０４では、メイン・ルーチン２０２は、行列Aのi行目について、ソート・ルーチン２０８を呼び出して、値が昇順になるインデックスの順列{a1,...,an}を得る。

ステップ３０６からステップ３１２までは、jの1からnまでの繰り返しのループである。ステップ３０８では、メイン・ルーチン２０２は、行列Aと行列BのFMM計算結果である行列Cについて、C[i,j] = ∞を格納する。ここで∞とは、実際の処理では現れないような十分大きい値のことである。行列Cも、行列A及び行列Bと同様に、n×nの正方行列である。

ステップ３１０では、メイン・ルーチン２０２は、更新ルーチン２１０を呼び出して、C[i,j]を{ai}を用いて更新する。ここで、{ai}とは、{a1,...,an}の略記である。ステップ３１０の詳細は、図４のフローチャートを参照して後で説明する。

メイン・ルーチン２０２は、ステップ３１２でjを1つインクリメントしてステップ３０６に戻る。ステップ３０６に戻って、メイン・ルーチン２０２は、jがn以下ならステップ３０８以下を繰り返す。jがnを超えたなら、ステップ３０６からステップ３１２までのループを抜ける。

ステップ３１２でループを抜けると、メイン・ルーチン２０２は、ステップ３１４でiを1つインクリメントしてステップ３０２に戻る。ステップ３０２に戻って、メイン・ルーチン２０２は、iがn以下ならステップ３０４以下を繰り返す。iがnを超えたなら、ステップ３０２からステップ３１４までのループを抜ける。

ステップ３１６からステップ３２６までは、iの1からnまでの繰り返しのループである。ステップ３１８では、メイン・ルーチン２０２は、行列Bのi列目について、ソート・ルーチン２０８を呼び出して、値が昇順になるインデックスの順列{b1,...,bn}を得る。

ステップ３２０からステップ３２４までは、jの1からnまでの繰り返しのループである。

ステップ３２２では、メイン・ルーチン２０２は、更新ルーチン２１０を呼び出して、C[i,j]を{bj}を用いて更新する。ここで、{bj}とは、{b1,...,bn}の略記である。ステップ３２２の詳細は、図５のフローチャートを参照して後で説明する。

メイン・ルーチン２０２は、ステップ３２４でjを1つインクリメントしてステップ３２０に戻る。ステップ３２０に戻って、メイン・ルーチン２０２は、jがn以下ならステップ３２２以下を繰り返す。jがnを超えたなら、ステップ３２０からステップ３２４までのループを抜ける。

ステップ３２４でループを抜けると、メイン・ルーチン２０２は、ステップ３２６でiを1つインクリメントしてステップ３１６に戻る。ステップ３１６に戻って、メイン・ルーチン２０２は、iがn以下ならステップ３１８以下を繰り返す。iがnを超えたなら、ステップ３１６からステップ３２６までのループを抜ける。

ここまでの処理が終了すると、C[i,j]の要素が全て得られている。メイン・ルーチン２０２は、出力ルーチン２１２を呼び出して、C[i,j]の要素の値を、結果データ２１４を含むファイルとして、ハードディスク・ドライブ１０８に書き出す。

図４は、図３のステップ３１０をより詳細に示すフローチャートである。図３では、この処理が、更新ルーチン２１０をサブルーチンとして呼び出す処理として説明されているが、図４のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。

図４において、更新ルーチン２１０は、ステップ４０２で、変数kを1とおく。更新ルーチン２１０は、ステップ４０４で、bestという変数に、C[i,j]の値を格納する。

ステップ４０６では、更新ルーチン２１０は、best = min{best, A[i,ak] + B[ak,j]}によりbestの値を更新し、ステップ４０８でkを1だけ増分する。

ステップ４１０では、更新ルーチン２１０は、k > nもしくはA[i,ak] > best/2であるかどうか判断し、そうでないならステップ４０６に戻る。

更新ルーチン２１０がステップ４１０で、k > nもしくはA[i,ak] > best/2であると判断すると、ステップ４１２でC[i,j] = bestと格納して、ステップ３１０が終了する。

図５は、図３のステップ３２２をより詳細に示すフローチャートである。図３では、この処理が、更新ルーチン２１０をサブルーチンとして呼び出す処理として説明されているが、図４のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。

図５において、更新ルーチン２１０は、ステップ５０２で、変数kを1とおく。更新ルーチン２１０は、ステップ５０４で、bestという変数に、C[i,j]の値を格納する。

ステップ５０６では、更新ルーチン２１０は、best = min{best, A[i,bk] + B[bk,j]}によりbestの値を更新し、ステップ５０８でkを1だけ増分する。

ステップ５１０では、更新ルーチン２１０は、k > nもしくはB[bk,j] > best/2であるかどうか判断し、そうでないならステップ５０６に戻る。

更新ルーチン２１０がステップ５１０で、k > nもしくはB[bk,j] > best/2であると判断すると、ステップ５１２でC[i,j] = bestと格納して、ステップ３２２が終了する。

この実施例は、図４のステップ４１０のA[i,ak] > best/2、あるいは、図５のステップ５１０のB[bk,j] > best/2の判断で、比較を早めに打ち切ることができるので、非特許文献１に記述されている技法よりも処理速度を向上することができる。

本発明のさらなる特徴は、SIMD(Single Instruction Multiple Data)命令を有利に利用して、処理速度を向上することができることである。その技法を図６を参照して説明する。

すなわち、本発明に従い、C[i,j]とC[i,j']（ここでj'= j+1）を計算するとき、
− C[i,j]については、まず、k = a1,a2,…の順にA[i,k] + B[k,j]を計算していく。
− C[i,j'] については、まず、k = a1,a2,… の順にA[i,k] + B[k,j']を計算していく。
この両者の計算では、行列Aについては同じ行にアクセスし、行列Bについては隣同士の値にアクセスしている。

よって、C[i,j] と C[i,j'] の計算を融合して、同時に行えば両者を一つのループで処理できる。そして、A[i,k] + B[k,j] と A[i,k] + B[k,j'] は、2-wayのSIMD命令を使えば一つのvec_add命令で処理できて高速化が可能となる。
− ループを融合する際、ループの終了条件を調整する必要があるが、これも容易に実現可能
− ループを融合した結果、ループ長が長くなる場合があるが、これによる計算速度の低下は一般にSIMD命令による高速化に比べて小さい。
C[i,j]とC[i,j']の計算にあたっては、次にk = b1, b2, … についても計算をする必要があるが、この両者の計算は融合できない。ただ、i' = i+1とすると、C[i,j]とC[i',j]の計算は k = b1, b2, … についてループの融合ができる。よって、k = a1, a2, … の計算には列方向に隣同士の行列 C の計算を融合し、k = b1, b2, … の計算には行方向に隣同士の行列 C の計算を融合する、という二段階で行う。非特許文献１の手法では、k = a1, a2, … の計算と k = b1, b2, … の計算、というように分解ができないのが問題である。

SIMDの実装方法は、これには限定されないが例えば、gcc、Visual C++などの処理系を使用して、emmintrin.hなどのヘッダファイルをインクルードし、__m128iなどのデータ型宣言を用いる。メモリからの読み出しは_mm_loadu_si128()を使用し、レジスタの初期化には_mm_set_epi32()を使用し、加算には_mm_add_epi32()を使用する、等である。

次に、SIMDを実装するのに好適な実施例の処理の全体を、図７のフローチャートを参照して説明する。この処理を開始するにあたって、メイン・ルーチン２０２は予め入力ルーチン２０４を呼び出し、処理データ２０６のファイルからデータを読み込むことによって、行列Aと行列Bの各要素の値を決定する。ここでも、行列Aと行列Bは、n×nの正方行列であると仮定する。また、SIMDの多重度をsとし、nを、sで割り切れる数であると仮定する。

さて図７において、ステップ７０２からステップ７１４までは、iの1からnまでの繰り返しのループである。ステップ７０４では、メイン・ルーチン２０２は、行列Aのi行目について、ソート・ルーチン２０８を呼び出して、値が昇順になるインデックスの順列{a1,...,an}を得る。

ステップ７０６からステップ７１２までは、jの1からn/sまでの繰り返しのループである。ステップ７０８では、メイン・ルーチン２０２は、行列Aと行列BのFMM計算結果である行列Cについて、C[i,(j-1)*s+1] = ∞、C[i,(j-1)*s+2] = ∞、...、C[i,(j-1)*s+s-1] = ∞と格納する。ここで、行列Cも、行列A及び行列Bと同様に、n×nの正方行列である。

ステップ７１０では、メイン・ルーチン２０２は、更新ルーチン２１０を呼び出して、C[i,(j-1)*s+1]、C[i,(j-1)*s+2]、...、C[i,(j-1)*s+s-1]を{ai}を用いて更新する。ステップ７１０の詳細は、図８のフローチャートを参照して後で説明する。

メイン・ルーチン２０２は、ステップ７１２でjを1つインクリメントしてステップ７０６に戻る。ステップ７０６に戻って、メイン・ルーチン２０２は、jがn/s以下ならステップ７０８以下を繰り返す。jがn/sを超えたなら、ステップ７０６からステップ７１２までのループを抜ける。

ステップ７１２でループを抜けると、メイン・ルーチン２０２は、ステップ７１４でiを1つインクリメントしてステップ７０２に戻る。ステップ７０２に戻って、メイン・ルーチン２０２は、iがn以下ならステップ７０４以下を繰り返す。iがnを超えたなら、ステップ７０２からステップ７１４までのループを抜ける。

ステップ７１６からステップ７２６までは、iの1からnまでの繰り返しのループである。ステップ７１８では、メイン・ルーチン２０２は、行列Bのi列目について、ソート・ルーチン２０８を呼び出して、値が昇順になるインデックスの順列{b1,...,bn}を得る。

ステップ７２０からステップ７２４までは、jの1からn/sまでの繰り返しのループである。

ステップ７２２では、メイン・ルーチン２０２は、更新ルーチン２１０を呼び出して、C[(i-1)*s+1,j]、C[(i-1)*s+2,j]、...、C[(i-1)*s+s-1,j]を{bj}を用いて更新する。ステップ７２２の詳細は、図９のフローチャートを参照して後で説明する。

メイン・ルーチン２０２は、ステップ７２４でjを1つインクリメントしてステップ７２０に戻る。ステップ７２０に戻って、メイン・ルーチン２０２は、jがn/s以下ならステップ７２２以下を繰り返す。jがn/sを超えたなら、ステップ７２０からステップ７２４までのループを抜ける。

ステップ７２４でループを抜けると、メイン・ルーチン２０２は、ステップ７２６でiを1つインクリメントしてステップ７１６に戻る。ステップ７１６に戻って、メイン・ルーチン２０２は、iがn以下ならステップ７１８以下を繰り返す。iがnを超えたなら、ステップ７１６からステップ７２６までのループを抜ける。

図８は、図７のステップ７１０をより詳細に示すフローチャートである。図７では、この処理が、更新ルーチン２１０をサブルーチンとして呼び出す処理として説明されているが、図８のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。

更新ルーチン２１０は、図８のステップ８０２では、変数kに、1を格納する。

ステップ８０４からステップ８０８までは、p = (j-1)*s + 1から(j-1)*s + s - 1までの繰り返しである。ステップ８０６では、更新ルーチン２１０は、t[p] = C[i,p]と値を格納する。このループがp = (j-1)*s + 1から(j-1)*s + s - 1まで終わると、ループを抜け出てステップ８１０に進む。

次は、ステップ８１０からステップ８１４までで、p = (j-1)*s + 1から(j-1)*s + s - 1までの繰り返しである。ステップ８１２では、更新ルーチン２１０は、t[p] = min{t[p],A[i,ak]+B{ak,p]}を実行する。このとき、複数のA[i,ak]+B{ak,p]の計算に並列的にvec_addのSIMD命令を使用して、処理速度が向上される。このループがp = (j-1)*s + 1から(j-1)*s + s - 1まで終わると、ループを抜け出てステップ８１６に進む。

更新ルーチン２１０は、ステップ８１６で、kを1だけ増分し、ステップ８１８で、k > nもしくは、A[i,ak] > max_p{t[p]/2}を判断する。この判断が否定的であるなら、処理はステップ８１０に戻る。

ステップ８１６の判断が肯定的なら、ステップ８２０からステップ８２４までのループを実行する。ステップ８２０からステップ８２４は、p = (j-1)*s + 1から(j-1)*s + s - 1までの繰り返しである。更新ルーチン２１０は、ステップ８２２で、C[i,p] = t[p]と格納する。このループがp = (j-1)*s + 1から(j-1)*s + s - 1まで終わると、ステップ７１０が終了する。

図９は、図７のステップ７２２をより詳細に示すフローチャートである。図７では、この処理が、更新ルーチン２１０をサブルーチンとして呼び出す処理として説明されているが、図９のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。

更新ルーチン２１０は、図９のステップ９０２では、変数kに、1を格納する。

ステップ９０４からステップ９０８までは、p = (i-1)*s + 1から(i-1)*s + s - 1までの繰り返しである。ステップ９０６では、更新ルーチン２１０は、t[p] = C[p,j]と値を格納する。このループがp = (i-1)*s + 1から(i-1)*s + s - 1まで終わると、ループを抜け出てステップ９１０に進む。

次は、ステップ９１０からステップ９１４までで、p = (i-1)*s + 1から(i-1)*s + s - 1までの繰り返しである。ステップ９１２では、更新ルーチン２１０は、t[p] = min{t[p],A[p,bk]+B[bk,j]}を実行する。このとき、複数のA[p,bk]+B[bk,j]の計算に並列的にvec_addのSIMD命令を使用して、処理速度が向上される。このループがp = (i-1)*s + 1から(i-1)*s + s - 1まで終わると、ループを抜け出てステップ９１６に進む。

更新ルーチン２１０は、ステップ９１６で、kを1だけ増分し、ステップ９１８で、k > nもしくは、B[bk,j] > max_p{t[p]/2}を判断する。この判断が否定的であるなら、処理はステップ９１０に戻る。

ステップ９１６の判断が肯定的なら、ステップ９２０からステップ９２４までのループを実行する。ステップ９２０からステップ９２４は、p = (i-1)*s + 1から(i-1)*s + s - 1までの繰り返しである。更新ルーチン２１０は、ステップ９２２で、C[p,j] = t[p]と格納する。このループがp = (i-1)*s + 1から(i-1)*s + s - 1まで終わると、ステップ７２２が終了する。

ところで、図４あるいは図５の処理を擬似コードで書くと、次のとおりである。
best = ∞
i = 1
repeat
best = min { best, va[ai] + vb[ai] }
i = i+1
until (i > n or va[ai] > best/2)
j = 1
repeat
best = min { best, va[bj] + vb[bj] }
j = j+1
until (j > n or vb[bj] > best/2)
output best

ここの(j > n or vb[bj] > best/2)という判定条件については、次のような実施例も考えられる。
best = ∞, temp = ∞
i = 1
repeat
best = min { best, va[ai] + vb[ai] }
temp = min { temp, vb[ai] }
i = i+1
until (i > n or va[ai] > best/2)
j = 1
repeat
best = min { best, va[bj] + vb[bj] }
j = j+1
until (j > n or vb[bj] > min { best/2, temp })
output best
best/2 ≧ min { best/2, temp }なので、こうするとループの脱出が早まる。ここでは、２つのループの打ち切り条件が異なることに留意されたい。

あるいは、以下のような例もありえる。
best = ∞
i = 1
repeat
best = min { best, va[ai] + vb[ai] }
i = i+1
until (i > n or va[ai] > best/2)
j = 1
repeat
best = min { best, va[bj] + vb[bj] }
j = j+1
until (j > n - 1 or vb[bj] > best/2)
output best
ここでも、２つのループの打ち切り条件が異なっている。

また、上記実施例では、行列A、B、Cはどれもn×nの正方行列であると想定したが、これには限定されず、通常の行列の掛け算と同様に、行列Aがm×k、行列Bがk×nであるとしてよい。この結果、行列Cは、m×nとなる。

以上、特定のハードウェアおよびソフトウェアのプラットフォーム上で実施するものとして本発明を説明してきたが、本発明は示されている例に限定されず、任意のコンピュータ・プラットフォーム上で実施可能である。

１０２システム・パス
１０４ＣＰＵ
１０６主記憶
１０８ハードディスク・ドライブ
１１０キーボード
１１２マウス
１１４ディスプレイ
２０２メイン・ルーチン
２０４入力ルーチン
２０６処理データ
２０８ソート・ルーチン
２１０更新ルーチン
２１２出力ルーチン
２１４結果データ

Claims

CPUとメモリとを備えるコンピュータの処理により、２つの行列（以下、それぞれA,Bとする）からFMM(Funny Matrix Multiplication)の結果である１つの行列（以下、Cとする）を計算する方法であって、
前記CPUにより、前記行列Cのi,j成分であるC[i,j]の値を所定の変数（best）とし、想定されるよりも大きい値を初期値として前記メモリに格納するステップと、
前記CPUにより、前記行列Aの最初の行から最後の行までの各々の行について、前記行列Aのi番目の行に対応し、値が昇順になるインデックスの順列{ai}を計算するステップと、
前記CPUにより、前記行列Bの最初の列から最後の列までの各々の列について、前記行列Bのｊ番目の列に対応し、値が昇順になるインデックスの順列{bj}を計算するステップと、
前記CPUにより、変数ｋをk = 1から1つずつ増分しながら順次、best= min{best, A[i,ak]+B[ak,j]}を計算し、ここで、akは前記インデックスの順列{ai}のk番目の要素であり、kが前記行列Aの行の数を超えるかA[i,ak]がbest/2を超えることに応答して、C[i,j] = bestによりC[i,j]を更新する第一更新ステップと、
前記CPUにより、変数ｋをk = 1から1つずつ増分しながら順次、best= min{best, A[i,bk]+B[bk,j]}を計算し、ここで、bkは前記インデックスの順列{bj}のk番目の要素であり、kが前記行列Bの列の数を超えるかB[bk,j]がbest/2を超えることに応答して、C[i,j] = bestによりC[i,j]を更新する第二更新ステップとを有する
方法。
前記行列Aを前記メモリ上に格納する際に、前記行列Aにおいて列方向に隣り合う要素を、前記メモリ上において隣り合うように並べ、
前記行列Bを前記メモリ上に格納する際に、前記行列Bにおいて行方向に隣り合う要素を、前記メモリ上において隣り合うように並べ、
前記第一更新ステップにおいて、A[i,ak]+B[ak,j]の計算が前記CPUのSIMD命令により並列実行され、
前記第二更新ステップにおいて、A[i,bk]+B[bk,j]の計算が前記CPUのSIMD命令により並列実行される
請求項１に記載の方法。
前記CPUのSIMD命令がアクセス可能なレジスタに前記行列Aにおいて列方向に隣り合う要素及び前記行列Bにおいて行方向に隣り合う要素が読み込まれる
請求項２に記載の方法。
前記SIMD命令は、ベクトル加算命令である
請求項２又は３に記載の方法。
前記第一更新ステップにおいて、行列Cの列方向に隣り合う要素の計算が並列化され、
前記第二更新ステップにおいて、行列Cの行方向に隣り合う要素の計算が並列化される
請求項１から４のいずれか１項に記載の方法。
前記行列A、B及びCが、n×nの正方行列である
請求項１から４のいずれか１項に記載の方法。
CPUとメモリとを備えるコンピュータに、２つの行列（以下、それぞれA,Bとする）からFMM(Funny MatrixMultiplication)の結果である１つの行列（以下、Cとする）を計算させるコンピュータ・プログラムであって、
前記コンピュータに、
前記CPUにより、前記行列Cのi,j成分であるC[i,j]の値を所定の変数（best）とし、想定されるよりも大きい値を初期値として前記メモリに格納するステップと、
前記CPUにより、前記行列Aの最初の行から最後の行までの各々の行について、前記行列Aのi番目の行に対応し、値が昇順になるインデックスの順列{ai}を計算するステップと、
前記CPUにより、前記行列Bの最初の列から最後の列までの各々の列について、前記行列Bのｊ番目の列に対応し、値が昇順になるインデックスの順列{bj}を計算するステップと、
前記CPUにより、変数ｋをk = 1から1つずつ増分しながら順次、best = min{best, A[i,ak]+B[ak,j]}を計算し、ここで、akは前記インデックスの順列{ai}のk番目の要素であり、kが前記行列Aの行の数を超えるかA[i,ak]がbest/2を超えることに応答して、C[i,j] = bestによりC[i,j]を更新する第一更新ステップと、
前記CPUにより、変数ｋをk = 1から1つずつ増分しながら順次、best = min{best, A[i,bk]+B[bk,j]}を計算し、ここで、bkは前記インデックスの順列{bj}のk番目の要素であり、kが前記行列Bの列の数を超えるかB[bk,j]がbest/2を超えることに応答して、C[i,j] = bestによりC[i,j]を更新する第二更新ステップ
とを実行させるコンピュータ・プログラム。
CPUとメモリとを備え、２つの行列（以下、それぞれA,Bとする）からFMM(Funny Matrix Multiplication)の結果である１つの行列（以下、Cとする）を計算するコンピュータ・システムであって、
前記CPUにより、前記行列Cのi,j成分であるC[i,j]の値を所定の変数（best）とし、想定されるよりも大きい値を初期値として前記メモリに格納する手段と、
前記CPUにより、前記行列Aの最初の行から最後の行までの各々の行について、前記行列Aのi番目の行に対応し、値が昇順になるインデックスの順列{ai}を計算する手段と、
前記CPUにより、前記行列Bの最初の列から最後の列までの各々の列について、前記行列Bのｊ番目の列に対応し、値が昇順になるインデックスの順列{bj}を計算する手段と、
前記CPUにより、変数ｋをk = 1から1つずつ増分しながら順次、best = min{best, A[i,ak]+B[ak,j]}を計算し、ここで、akは前記インデックスの順列{ai}のk番目の要素であり、kが前記行列Aの行の数を超えるかA[i,ak]がbest/2を超えることに応答して、C[i,j] = bestによりC[i,j]を更新する第一更新手段と、
前記CPUにより、変数ｋをk = 1から1つずつ増分しながら順次、best = min{best, A[i,bk]+B[bk,j]}を計算し、ここで、bkは前記インデックスの順列{bj}のk番目の要素であり、kが前記行列Bの列の数を超えるかB[bk,j]がbest/2を超えることに応答して、C[i,j] = bestによりC[i,j]を更新する第二更新手段とを有する、
コンピュータ・システム。
前記行列Aを前記メモリ上に格納する際に、前記行列Aにおいて列方向に隣り合う要素を、前記メモリ上において隣り合うように並べ、
前記行列Bを前記メモリ上に格納する際に、前記行列Bにおいて行方向に隣り合う要素を、前記メモリ上において隣り合うように並べ、
前記第一更新手段は、A[i,ak]+B[ak,j]の計算を前記CPUのSIMD命令により並列実行し、
前記第二更新手段は、A[i,bk]+B[bk,j]の計算を前記CPUのSIMD命令により並列実行する
請求項８に記載のコンピュータ・システム。