行列計算処理方法、プログラム及びシステム
【課題】 FMM(Funny Matrix Multiplication)の計算を高速化すること。
【解決手段】 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で中心となる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(Funny Matrix Multiplication)の計算に関する。
【背景技術】
【０００２】
FMMとは、通常の行列乗算処理におけるadd（加算）とmultiply（乗算）を、それぞれ、min演算とadd演算に置き換えた行列計算処理のことである。すなわち、２つの入力行列A、Bについて、C[i,j] = mink{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の順に計算していき、そのうちの最小値をとる計算が基本となる。つまり、mink{va[k] + vb[k]}を計算する処理が中心となり、この処理をn2回行うので、
FMM全体の処理にO(n3)時間かかる。
【０００５】
ところで、近年、交通シミュレーションやカー・ナビゲーションなどの分野で、最短経路計算をより高速に行うことの要望が高まっており、そのため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について、mink{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(n2log 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の計算において、上記非特許文献１に開示されている技法を更に改良して高速にすることにある。
【課題を解決するための手段】
【０００９】
本発明においても、mink{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を超えたなら、ステップ７１６からステップ７２６までのループを抜ける。
【００６６】
ここまでの処理が終了すると、C[i,j]の要素が全て得られている。メイン・ルーチン２０２は、出力ルーチン２１２を呼び出して、C[i,j]の要素の値を、結果データ２１４を含むファイルとして、ハードディスク・ドライブ１０８に書き出す。
【００６７】
図８は、図７のステップ７１０をより詳細に示すフローチャートである。図７では、この処理が、更新ルーチン２１０をサブルーチンとして呼び出す処理として説明されているが、図８のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。
【００６８】
更新ルーチン２１０は、図８のステップ８０２では、変数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] > maxp{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] > maxp{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となる。
【００８３】
以上、特定のハードウェアおよびソフトウェアのプラットフォーム上で実施するものとして本発明を説明してきたが、本発明は示されている例に限定されず、任意のコンピュータ・プラットフォーム上で実施可能である。
【符号の説明】
【００８４】
１０２ システム・パス
１０４ ＣＰＵ
１０６ 主記憶
１０８ ハードディスク・ドライブ
１１０ キーボード
１１２ マウス
１１４ ディスプレイ
２０２ メイン・ルーチン
２０４ 入力ルーチン
２０６ 処理データ
２０８ ソート・ルーチン
２１０ 更新ルーチン
２１２ 出力ルーチン
２１４ 結果データ
【技術分野】
【０００１】
この発明は、コンピュータの処理によって、行列の計算を行う技法に関し、より詳しくは、FMM(Funny Matrix Multiplication)の計算に関する。
【背景技術】
【０００２】
FMMとは、通常の行列乗算処理におけるadd（加算）とmultiply（乗算）を、それぞれ、min演算とadd演算に置き換えた行列計算処理のことである。すなわち、２つの入力行列A、Bについて、C[i,j] = mink{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の順に計算していき、そのうちの最小値をとる計算が基本となる。つまり、mink{va[k] + vb[k]}を計算する処理が中心となり、この処理をn2回行うので、
FMM全体の処理にO(n3)時間かかる。
【０００５】
ところで、近年、交通シミュレーションやカー・ナビゲーションなどの分野で、最短経路計算をより高速に行うことの要望が高まっており、そのため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について、mink{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(n2log 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の計算において、上記非特許文献１に開示されている技法を更に改良して高速にすることにある。
【課題を解決するための手段】
【０００９】
本発明においても、mink{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を超えたなら、ステップ７１６からステップ７２６までのループを抜ける。
【００６６】
ここまでの処理が終了すると、C[i,j]の要素が全て得られている。メイン・ルーチン２０２は、出力ルーチン２１２を呼び出して、C[i,j]の要素の値を、結果データ２１４を含むファイルとして、ハードディスク・ドライブ１０８に書き出す。
【００６７】
図８は、図７のステップ７１０をより詳細に示すフローチャートである。図７では、この処理が、更新ルーチン２１０をサブルーチンとして呼び出す処理として説明されているが、図８のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。
【００６８】
更新ルーチン２１０は、図８のステップ８０２では、変数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] > maxp{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] > maxp{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となる。
【００８３】
以上、特定のハードウェアおよびソフトウェアのプラットフォーム上で実施するものとして本発明を説明してきたが、本発明は示されている例に限定されず、任意のコンピュータ・プラットフォーム上で実施可能である。
【符号の説明】
【００８４】
１０２ システム・パス
１０４ ＣＰＵ
１０６ 主記憶
１０８ ハードディスク・ドライブ
１１０ キーボード
１１２ マウス
１１４ ディスプレイ
２０２ メイン・ルーチン
２０４ 入力ルーチン
２０６ 処理データ
２０８ ソート・ルーチン
２１０ 更新ルーチン
２１２ 出力ルーチン
２１４ 結果データ
【特許請求の範囲】
【請求項１】
コンピュータの処理により、２つの行列（以下、それぞれA,Bとする）のFMM(Funny Matrix Multiplication)を計算する方法であって、
i = 1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列{ai}を順に計算するステップと、
前記i番目の行において、j = 1から前記行列Aの列の数までについて、C[i,j]に先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのi,j成分であるC[i,j]の値を所定の変数（best）に格納し、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]を更新するステップと、
j = 1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列{bj}を順に計算するステップと、
前記j番目の列について、j = 1から前記行列Bの列の数までについて、前記行列Cのi,j成分であるC[i,j]をbestとおき、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が列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、
前記A[i,ak] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップと、前記B[bk,j] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップが、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項１に記載の方法。
【請求項３】
前記行列A、B及びCが、n×nの正方行列である、請求項１に記載の方法。
【請求項４】
コンピュータの処理により、２つの行列（以下、それぞれA,Bとする）のFMM(Funny Matrix Multiplication)を計算するプログラムであって、
前記コンピュータに、
i = 1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列{ai}を順に計算するステップと、
前記i番目の行において、j = 1から前記行列Aの列の数までについて、C[i,j]に先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのi,j成分であるC[i,j]の値を所定の変数（best）に格納し、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]を更新するステップと、
j = 1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列{bj}を順に計算するステップと、
前記j番目の列について、j = 1から前記行列Bの列の数までについて、前記行列Cのi,j成分であるC[i,j]をbestとおき、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が列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、
前記A[i,ak] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップと、前記B[bk,j] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップが、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項４に記載のプログラム。
【請求項６】
前記行列A、B及びCが、n×nの正方行列である、請求項４に記載のプログラム。
【請求項７】
コンピュータの処理により、２つの行列（以下、それぞれA,Bとする）のFMM(Funny Matrix Multiplication)を計算するシステムであって、
i = 1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列{ai}を順に計算する手段と、
前記i番目の行について、j = 1から前記行列Aの列の数までについて、C[i,j]に先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのi,j成分であるC[i,j]の値を所定の変数（best）に格納し、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]を更新する手段と、
j = 1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列{bj}を順に計算する手段と、
前記j番目の列について、j = 1から前記行列Bの列の数までについて、前記行列Cのi,j成分であるC[i,j]をbestとおき、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が列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、
前記A[i,ak] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新する手段と、前記B[bk,j] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新する手段が、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項７に記載のシステム。
【請求項９】
前記行列A、B及びCが、n×nの正方行列である、請求項７に記載のシステム。
【請求項１】
コンピュータの処理により、２つの行列（以下、それぞれA,Bとする）のFMM(Funny Matrix Multiplication)を計算する方法であって、
i = 1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列{ai}を順に計算するステップと、
前記i番目の行において、j = 1から前記行列Aの列の数までについて、C[i,j]に先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのi,j成分であるC[i,j]の値を所定の変数（best）に格納し、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]を更新するステップと、
j = 1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列{bj}を順に計算するステップと、
前記j番目の列について、j = 1から前記行列Bの列の数までについて、前記行列Cのi,j成分であるC[i,j]をbestとおき、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が列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、
前記A[i,ak] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップと、前記B[bk,j] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップが、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項１に記載の方法。
【請求項３】
前記行列A、B及びCが、n×nの正方行列である、請求項１に記載の方法。
【請求項４】
コンピュータの処理により、２つの行列（以下、それぞれA,Bとする）のFMM(Funny Matrix Multiplication)を計算するプログラムであって、
前記コンピュータに、
i = 1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列{ai}を順に計算するステップと、
前記i番目の行において、j = 1から前記行列Aの列の数までについて、C[i,j]に先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのi,j成分であるC[i,j]の値を所定の変数（best）に格納し、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]を更新するステップと、
j = 1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列{bj}を順に計算するステップと、
前記j番目の列について、j = 1から前記行列Bの列の数までについて、前記行列Cのi,j成分であるC[i,j]をbestとおき、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が列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、
前記A[i,ak] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップと、前記B[bk,j] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新するステップが、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項４に記載のプログラム。
【請求項６】
前記行列A、B及びCが、n×nの正方行列である、請求項４に記載のプログラム。
【請求項７】
コンピュータの処理により、２つの行列（以下、それぞれA,Bとする）のFMM(Funny Matrix Multiplication)を計算するシステムであって、
i = 1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列{ai}を順に計算する手段と、
前記i番目の行について、j = 1から前記行列Aの列の数までについて、C[i,j]に先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのi,j成分であるC[i,j]の値を所定の変数（best）に格納し、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]を更新する手段と、
j = 1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列{bj}を順に計算する手段と、
前記j番目の列について、j = 1から前記行列Bの列の数までについて、前記行列Cのi,j成分であるC[i,j]をbestとおき、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が列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、
前記A[i,ak] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新する手段と、前記B[bk,j] > best/2であることに応答して、C[i,j] = bestによりC[i,j]を更新する手段が、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項７に記載のシステム。
【請求項９】
前記行列A、B及びCが、n×nの正方行列である、請求項７に記載のシステム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１２−１６４０３９（Ｐ２０１２−１６４０３９Ａ）
【公開日】平成２４年８月３０日（２０１２．８．３０）
【国際特許分類】
【出願番号】特願２０１１−２２３１１（Ｐ２０１１−２２３１１）
【出願日】平成２３年２月４日（２０１１．２．４）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２１年度、総務省、「自動車二酸化炭素排出量削減のための大規模モビリティ社会シミュレータの研究開発」委託事業、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（３９０００９５３１）インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】ＩＮＴＥＲＮＡＴＩＯＮＡＬ ＢＵＳＩＮＥＳＳ ＭＡＳＣＨＩＮＥＳ ＣＯＲＰＯＲＡＴＩＯＮ
【Ｆターム（参考）】
【公開日】平成２４年８月３０日（２０１２．８．３０）
【国際特許分類】
【出願日】平成２３年２月４日（２０１１．２．４）
【国等の委託研究の成果に係る記載事項】（出願人による申告）平成２１年度、総務省、「自動車二酸化炭素排出量削減のための大規模モビリティ社会シミュレータの研究開発」委託事業、産業技術力強化法第１９条の適用を受ける特許出願
【出願人】（３９０００９５３１）インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】ＩＮＴＥＲＮＡＴＩＯＮＡＬ ＢＵＳＩＮＥＳＳ ＭＡＳＣＨＩＮＥＳ ＣＯＲＰＯＲＡＴＩＯＮ
【Ｆターム（参考）】
[ Back to top ]