コンピュータの処理により、2つの行列以下、それぞれABとするのFMMを計算する方法であって、i1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列aiを順に計算するステップと、前記i番目の行において、j1から前記行列Aの列の数までについて、Cijに先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのij成分であるCijの値を所定の変数bestに格納し、k1から1つずつ増分しながら順次、bestminbestAiakBakjを計算し、ここで、akは前記インデックスの順列aiのk番目の要素であり、kが前記行列Aの行の数を超えるかAiakがbest2を超えることに応答して、CijbestによりCijを更新するステップと、j1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列bjを順に計算するステップと、前記j番目の列について、j1から前記行列Bの列の数までについて、前記行列Cのij成分であるCijをbestとおき、k1から1つずつ増分しながら順次、bestminbestAibkBbkjを計算し、ここで、bkは前記インデックスの順列bjのk番目の要素であり、kが前記行列Bの列の数を超えるかBbkjがbest2を超えることに応答して、CijbestによりCijを更新するステップを有する、行列計算処理方法。前記行列Aが列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、前記Aiakbest2であることに応答して、CijbestによりCijを更新するステップと、前記Bbkjbest2であることに応答して、CijbestによりCijを更新するステップが、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項1に記載の方法。前記行列A、B及びCが、n×nの正方行列である、請求項1に記載の方法。コンピュータの処理により、2つの行列以下、それぞれABとするのFMMを計算するプログラムであって、前記コンピュータに、i1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列aiを順に計算するステップと、前記i番目の行において、j1から前記行列Aの列の数までについて、Cijに先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのij成分であるCijの値を所定の変数bestに格納し、k1から1つずつ増分しながら順次、bestminbestAiakBakjを計算し、ここで、akは前記インデックスの順列aiのk番目の要素であり、kが前記行列Aの行の数を超えるかAiakがbest2を超えることに応答して、CijbestによりCijを更新するステップと、j1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列bjを順に計算するステップと、前記j番目の列について、j1から前記行列Bの列の数までについて、前記行列Cのij成分であるCijをbestとおき、k1から1つずつ増分しながら順次、bestminbestAibkBbkjを計算し、ここで、bkは前記インデックスの順列bjのk番目の要素であり、kが前記行列Bの列の数を超えるかBbkjがbest2を超えることに応答して、CijbestによりCijを更新するステップを実行させる、行列計算処理プログラム。前記行列Aが列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、前記Aiakbest2であることに応答して、CijbestによりCijを更新するステップと、前記Bbkjbest2であることに応答して、CijbestによりCijを更新するステップが、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項4に記載のプログラム。前記行列A、B及びCが、n×nの正方行列である、請求項4に記載のプログラム。コンピュータの処理により、2つの行列以下、それぞれABとするのFMMを計算するシステムであって、i1から前記行列Aの行の数までにおいて、各々の行について、値が昇順になるインデックスの順列aiを順に計算する手段と、前記i番目の行について、j1から前記行列Aの列の数までについて、Cijに先ず行列の値として想定されるよりも十分大きい値を格納し、FMMの計算結果である結果の行列Cのij成分であるCijの値を所定の変数bestに格納し、k1から1つずつ増分して順次、bestminbestAiakBakjを計算し、ここで、akは前記インデックスの順列aiのk番目の要素であり、kが前記行列Aの行の数を超えるかAiakがbest2を超えることに応答して、CijbestによりCijを更新する手段と、j1から前記行列Bの列の数までにおいて、各々の列について、値が昇順になるインデックスの順列bjを順に計算する手段と、前記j番目の列について、j1から前記行列Bの列の数までについて、前記行列Cのij成分であるCijをbestとおき、k1から1つずつ増分して順次、bestminbestAibkBbkjを計算し、ここで、bkは前記インデックスの順列bjのk番目の要素であり、kが前記行列Bの列の数を超えるかBbkjがbest2を超えることに応答して、CijbestによりCijを更新する手段を有する、行列計算処理システム。前記行列Aが列を主のレイアウトになるようにするとともに、前記行列Bが行を主のレイアウトになるようにし、前記Aiakbest2であることに応答して、CijbestによりCijを更新する手段と、前記Bbkjbest2であることに応答して、CijbestによりCijを更新する手段が、ループを融合して分割された部分をそれぞれ、SIMD命令により並列実行される、請求項7に記載のシステム。前記行列A、B及びCが、n×nの正方行列である、請求項7に記載のシステム。この発明の目的は、FMMの計算において、上記非特許文献1に開示されている技法を更に改良して高速にすることにある。本発明においても、minkvakvbkを計算する前に、vaa1≦vaa2≦...≦vaanvbb1≦vbb2≦...≦vbbnとなるような順列a1a2...anとb1b2...bnを事前に計算することは行われる。本発明の特徴は、FMMで中心となるvakvbkの最小値を計算する処理において、k1...nについて順番に処理するのではなく、best∞に初期設定してから、以下の処理Xと処理Yを一回ずつ適用した上で、bestの値をvakvbkの最小値として出力することにある。ka1a2...の順にvakvbkの値を計算していき、それまでに見つかった最小値をbestとしたときに、vakbest2となるkで処理をやめるそのようなkが無ければkanまで処理したらやめる処理Xと同様の処理をkb1b2...についても行い、vbkbest2の値となるkで処理をやめるそのようなkが無ければkbnまで処理したらやめるこのような処理を採用したことにより、FMMの計算において、SIMD命令を利用して処理を高速化することが可能となった。この場合、行列Aと行列BのFMMを計算するとすると、行列Aを、列が主のレイアウト、行列Bを、行が主のレイアウトで保存するのがポイントである。ここで、列が主のレイアウトとは、行列Aをメモリ上に保管する際に行列Aで列方向に隣り合う要素をメモリ上でもなるべく隣り合うように並べる方法のことをいう。つまり、行列Aを、A11A21A31An1A12A22A32An2A13A23Annといった順番で並べる。一方、行が主のレイアウトとは、列が主のレイアウトにおいて、列と行を入れ替えたものであり、つまり、行列Aを、A11A12A13A1nA21A22A23A2nA31A32Annといった順番で並べる。この発明によれば、bestの値を一旦計算して、その値に基づき、最小値の計算を打ち切るようにすることにより、最小値の計算を速くすることによって、FMMの計算処理を高速化することができる。また、非特許文献1の技法ではSIMD命令を利用することは困難であったが、本発明の技法では、SIMD命令を有利に適用して、更にFMMの計算処理を高速化することが可能である。以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。図1を参照すると、本発明の一実施例に係るシステム構成及び処理を実現するためのコンピュータハードウェアのブロック図が示されている。図1において、システムパス102には、CPU104と、主記憶RAM106と、ハードディスクドライブHDD108と、キーボード110と、マウス112と、ディスプレイ114が接続されている。CPU104は、好適には、32ビットまたは64ビットのアーキテクチャに基づくものであり、例えば、インテル社のPentium商標4、Core2Duo、Xeon、AMD社のAthlon商標などを使用することができる。この実施例の目的のため、CPU104は、SIMD命令をもつものである。主記憶106は、好適には、4GB以上の容量をもつものである。ハードディスクドライブ108は、行列要素の基となる大量のデータを格納できるように、例えば、500GB以上の容量をもつものであることが望ましい。ハードディスクドライブ108には、個々に図示しないが、オペレーティングシステムが、予め格納されている。オペレーティングシステムは、Linux、マイクロソフト社のWindowsXP、Windows7、アップルコンピュータのMacOS商標などの、CPU104に適合する任意のものでよい。ハードディスクドライブ108にはまた、好適には、C、C、C、Javaなどのプログラム言語処理系も格納されている。このプログラム言語処理系は、後で説明する、本発明に係るFMMの計算処理のためのモジュールを作成し、保守するために使用される。ハードディスクドライブ108はさらに、プログラム言語処理系でコンパイルするためのソースコードを書くためのテキストエディタ、及び、Eclipse商標などの開発環境を含んでいてもよい。キーボード110及びマウス112は、オペレーティングシステムまたは、ハードディスクドライブ108から主記憶106にロードされ、ディスプレイ114に表示されたプログラム図示しないを起動したり、文字を打ち込んだりするために使用される。ディスプレイ114は、好適には液晶ディスプレイであり、例えば、XGA1024×768の解像度、またはUXGA1600×1200の解像度などの任意の解像度のものを使用することができる。ディスプレイ114は、図示しないが、本発明の処理を開始するための操作ウインドウや、FMMの計算結果を表示するために使用される。次に、図2の機能ブロック図を参照して、本発明の処理を実行するための処理ルーチンについて説明する。これらの処理ルーチンはC、Cなどの、好適にはSIMD命令を利用可能なプログラム言語で作成されて、実行可能な形式でハードディスクドライブ108に保存され、オペレーティングシステムの動作で主記憶106にロードされて実行される。メインルーチン202は、本発明の全体の動作を統合するためのプログラムであり、図示しないが、ディスプレイ114に操作ウインドウを表示したり、ユーザの操作を受け付けて処理を開始したりする機能を有する。入力ルーチン204は、ハードディスクドライブ108に保存された地図データなどである処理データ206のファイルからデータを読み込んで、行列Aと行列Bの各要素を決定する機能をもつ。ソートルーチン208は、本発明に従いFMMの計算を行うためのインデックスの順列を求めるためのソートを行う機能をもつ。この際のソートのアルゴリズムは、クイックソートであるが、これには限定されず、シェルソート、ヒープソート、マージソートなどの任意の適当なソートアルゴリズムを用いることができる。更新ルーチン210は、ソートされたインデックスの順列を用いて、FMMの結果の行列Cの要素を更新するための処理を実行する機能をもつ。更新ルーチン210の処理の詳細は、図3~図5のフローチャート、あるいは図7~図9のフローチャートを参照して、後で説明する。出力ルーチン212は、行列Aと行列BとのFMM計算の結果得られた行列Cを、結果データ214のファイルとして書き出す機能をもつ。次に、図3のフローチャートを参照して、本発明による全体のFMM計算処理について説明する。この処理を開始するにあたって、メインルーチン202は予め入力ルーチン204を呼び出し、処理データ206のファイルからデータを読み込むことによって、行列Aと行列Bの各要素の値を決定する。ここでは、行列Aと行列Bは、n×nの正方行列であると仮定する。しかし、本発明の処理は、正方行列以外の行列にも適用できることを理解されたい。さて図3において、ステップ302からステップ314までは、iの1からnまでの繰り返しのループである。ステップ304では、メインルーチン202は、行列Aのi行目について、ソートルーチン208を呼び出して、値が昇順になるインデックスの順列a1anを得る。ステップ306からステップ312までは、jの1からnまでの繰り返しのループである。ステップ308では、メインルーチン202は、行列Aと行列BのFMM計算結果である行列Cについて、Cij∞を格納する。ここで∞とは、実際の処理では現れないような十分大きい値のことである。行列Cも、行列A及び行列Bと同様に、n×nの正方行列である。ステップ310では、メインルーチン202は、更新ルーチン210を呼び出して、Cijをaiを用いて更新する。ここで、aiとは、a1anの略記である。ステップ310の詳細は、図4のフローチャートを参照して後で説明する。メインルーチン202は、ステップ312でjを1つインクリメントしてステップ306に戻る。ステップ306に戻って、メインルーチン202は、jがn以下ならステップ308以下を繰り返す。jがnを超えたなら、ステップ306からステップ312までのループを抜ける。ステップ312でループを抜けると、メインルーチン202は、ステップ314でiを1つインクリメントしてステップ302に戻る。ステップ302に戻って、メインルーチン202は、iがn以下ならステップ304以下を繰り返す。iがnを超えたなら、ステップ302からステップ314までのループを抜ける。ステップ316からステップ326までは、iの1からnまでの繰り返しのループである。ステップ318では、メインルーチン202は、行列Bのi列目について、ソートルーチン208を呼び出して、値が昇順になるインデックスの順列b1bnを得る。ステップ320からステップ324までは、jの1からnまでの繰り返しのループである。ステップ322では、メインルーチン202は、更新ルーチン210を呼び出して、Cijをbjを用いて更新する。ここで、bjとは、b1bnの略記である。ステップ322の詳細は、図5のフローチャートを参照して後で説明する。メインルーチン202は、ステップ324でjを1つインクリメントしてステップ320に戻る。ステップ320に戻って、メインルーチン202は、jがn以下ならステップ322以下を繰り返す。jがnを超えたなら、ステップ320からステップ324までのループを抜ける。ステップ324でループを抜けると、メインルーチン202は、ステップ326でiを1つインクリメントしてステップ316に戻る。ステップ316に戻って、メインルーチン202は、iがn以下ならステップ318以下を繰り返す。iがnを超えたなら、ステップ316からステップ326までのループを抜ける。ここまでの処理が終了すると、Cijの要素が全て得られている。メインルーチン202は、出力ルーチン212を呼び出して、Cijの要素の値を、結果データ214を含むファイルとして、ハードディスクドライブ108に書き出す。図4は、図3のステップ310をより詳細に示すフローチャートである。図3では、この処理が、更新ルーチン210をサブルーチンとして呼び出す処理として説明されているが、図4のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。図4において、更新ルーチン210は、ステップ402で、変数kを1とおく。更新ルーチン210は、ステップ404で、bestという変数に、Cijの値を格納する。ステップ406では、更新ルーチン210は、bestminbestAiakBakjによりbestの値を更新し、ステップ408でkを1だけ増分する。ステップ410では、更新ルーチン210は、knもしくはAiakbest2であるかどうか判断し、そうでないならステップ406に戻る。更新ルーチン210がステップ410で、knもしくはAiakbest2であると判断すると、ステップ412でCijbestと格納して、ステップ310が終了する。図5は、図3のステップ322をより詳細に示すフローチャートである。図3では、この処理が、更新ルーチン210をサブルーチンとして呼び出す処理として説明されているが、図4のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。図5において、更新ルーチン210は、ステップ502で、変数kを1とおく。更新ルーチン210は、ステップ504で、bestという変数に、Cijの値を格納する。ステップ506では、更新ルーチン210は、bestminbestAibkBbkjによりbestの値を更新し、ステップ508でkを1だけ増分する。ステップ510では、更新ルーチン210は、knもしくはBbkjbest2であるかどうか判断し、そうでないならステップ506に戻る。更新ルーチン210がステップ510で、knもしくはBbkjbest2であると判断すると、ステップ512でCijbestと格納して、ステップ322が終了する。この実施例は、図4のステップ410のAiakbest2、あるいは、図5のステップ510のBbkjbest2の判断で、比較を早めに打ち切ることができるので、非特許文献1に記述されている技法よりも処理速度を向上することができる。本発明のさらなる特徴は、SIMD命令を有利に利用して、処理速度を向上することができることである。その技法を図6を参照して説明する。すなわち、本発明に従い、CijとCijここでjj1を計算するとき、-Cijについては、まず、ka1a2...の順にAikBkjを計算していく。-Cijについては、まず、ka1a2...の順にAikBkjを計算していく。この両者の計算では、行列Aについては同じ行にアクセスし、行列Bについては隣同士の値にアクセスしている。よって、CijとCijの計算を融合して、同時に行えば両者を一つのループで処理できる。そして、AikBkjとAikBkjは、2wayのSIMD命令を使えば一つのvecadd命令で処理できて高速化が可能となる。-ループを融合する際、ループの終了条件を調整する必要があるが、これも容易に実現可能-ループを融合した結果、ループ長が長くなる場合があるが、これによる計算速度の低下は一般にSIMD命令による高速化に比べて小さい。CijとCijの計算にあたっては、次にkb1b2...についても計算をする必要があるが、この両者の計算は融合できない。ただ、ii1とすると、CijとCijの計算はkb1b2...についてループの融合ができる。よって、ka1a2...の計算には列方向に隣同士の行列Cの計算を融合し、kb1b2...の計算には行方向に隣同士の行列Cの計算を融合する、という二段階で行う。非特許文献1の手法では、ka1a2...の計算とkb1b2...の計算、というように分解ができないのが問題である。SIMDの実装方法は、これには限定されないが例えば、gcc、VisualCなどの処理系を使用して、emmintrinhなどのヘッダファイルをインクルードし、m128iなどのデータ型宣言を用いる。メモリからの読み出しはmmloadusi128を使用し、レジスタの初期化にはmmsetepi32を使用し、加算にはmmaddepi32を使用する、等である。次に、SIMDを実装するのに好適な実施例の処理の全体を、図7のフローチャートを参照して説明する。この処理を開始するにあたって、メインルーチン202は予め入力ルーチン204を呼び出し、処理データ206のファイルからデータを読み込むことによって、行列Aと行列Bの各要素の値を決定する。ここでも、行列Aと行列Bは、n×nの正方行列であると仮定する。また、SIMDの多重度をsとし、nを、sで割り切れる数であると仮定する。さて図7において、ステップ702からステップ714までは、iの1からnまでの繰り返しのループである。ステップ704では、メインルーチン202は、行列Aのi行目について、ソートルーチン208を呼び出して、値が昇順になるインデックスの順列a1anを得る。ステップ706からステップ712までは、jの1からnsまでの繰り返しのループである。ステップ708では、メインルーチン202は、行列Aと行列BのFMM計算結果である行列Cについて、Cis1∞、Cis2∞、、Ciss1∞と格納する。ここで、行列Cも、行列A及び行列Bと同様に、n×nの正方行列である。ステップ710では、メインルーチン202は、更新ルーチン210を呼び出して、Cis1、Cis2、、Ciss1をaiを用いて更新する。ステップ710の詳細は、図8のフローチャートを参照して後で説明する。メインルーチン202は、ステップ712でjを1つインクリメントしてステップ706に戻る。ステップ706に戻って、メインルーチン202は、jがns以下ならステップ708以下を繰り返す。jがnsを超えたなら、ステップ706からステップ712までのループを抜ける。ステップ712でループを抜けると、メインルーチン202は、ステップ714でiを1つインクリメントしてステップ702に戻る。ステップ702に戻って、メインルーチン202は、iがn以下ならステップ704以下を繰り返す。iがnを超えたなら、ステップ702からステップ714までのループを抜ける。ステップ716からステップ726までは、iの1からnまでの繰り返しのループである。ステップ718では、メインルーチン202は、行列Bのi列目について、ソートルーチン208を呼び出して、値が昇順になるインデックスの順列b1bnを得る。ステップ720からステップ724までは、jの1からnsまでの繰り返しのループである。ステップ722では、メインルーチン202は、更新ルーチン210を呼び出して、Cs1j、Cs2j、、Css1jをbjを用いて更新する。ステップ722の詳細は、図9のフローチャートを参照して後で説明する。メインルーチン202は、ステップ724でjを1つインクリメントしてステップ720に戻る。ステップ720に戻って、メインルーチン202は、jがns以下ならステップ722以下を繰り返す。jがnsを超えたなら、ステップ720からステップ724までのループを抜ける。ステップ724でループを抜けると、メインルーチン202は、ステップ726でiを1つインクリメントしてステップ716に戻る。ステップ716に戻って、メインルーチン202は、iがn以下ならステップ718以下を繰り返す。iがnを超えたなら、ステップ716からステップ726までのループを抜ける。ここまでの処理が終了すると、Cijの要素が全て得られている。メインルーチン202は、出力ルーチン212を呼び出して、Cijの要素の値を、結果データ214を含むファイルとして、ハードディスクドライブ108に書き出す。図8は、図7のステップ710をより詳細に示すフローチャートである。図7では、この処理が、更新ルーチン210をサブルーチンとして呼び出す処理として説明されているが、図8のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。更新ルーチン210は、図8のステップ802では、変数kに、1を格納する。ステップ804からステップ808までは、ps1からss1までの繰り返しである。ステップ806では、更新ルーチン210は、tpCipと値を格納する。このループがps1からss1まで終わると、ループを抜け出てステップ810に進む。次は、ステップ810からステップ814までで、ps1からss1までの繰り返しである。ステップ812では、更新ルーチン210は、tpmintpAiakBakpを実行する。このとき、複数のAiakBakpの計算に並列的にvecaddのSIMD命令を使用して、処理速度が向上される。このループがps1からss1まで終わると、ループを抜け出てステップ816に進む。更新ルーチン210は、ステップ816で、kを1だけ増分し、ステップ818で、knもしくは、Aiakmaxptp2を判断する。この判断が否定的であるなら、処理はステップ810に戻る。ステップ816の判断が肯定的なら、ステップ820からステップ824までのループを実行する。ステップ820からステップ824は、ps1からss1までの繰り返しである。更新ルーチン210は、ステップ822で、Ciptpと格納する。このループがps1からss1まで終わると、ステップ710が終了する。図9は、図7のステップ722をより詳細に示すフローチャートである。図7では、この処理が、更新ルーチン210をサブルーチンとして呼び出す処理として説明されているが、図9のフローチャートに対応するコードがインライン的に埋め込まれて実行されるようにしてもよい。更新ルーチン210は、図9のステップ902では、変数kに、1を格納する。ステップ904からステップ908までは、ps1からss1までの繰り返しである。ステップ906では、更新ルーチン210は、tpCpjと値を格納する。このループがps1からss1まで終わると、ループを抜け出てステップ910に進む。次は、ステップ910からステップ914までで、ps1からss1までの繰り返しである。ステップ912では、更新ルーチン210は、tpmintpApbkBbkjを実行する。このとき、複数のApbkBbkjの計算に並列的にvecaddのSIMD命令を使用して、処理速度が向上される。このループがps1からss1まで終わると、ループを抜け出てステップ916に進む。更新ルーチン210は、ステップ916で、kを1だけ増分し、ステップ918で、knもしくは、Bbkjmaxptp2を判断する。この判断が否定的であるなら、処理はステップ910に戻る。ステップ916の判断が肯定的なら、ステップ920からステップ924までのループを実行する。ステップ920からステップ924は、ps1からss1までの繰り返しである。更新ルーチン210は、ステップ922で、Cpjtpと格納する。このループがps1からss1まで終わると、ステップ722が終了する。ところで、図4あるいは図5の処理を擬似コードで書くと、次のとおりである。best∞i1repeatbestminbestvaaivbaiii1untilj1repeatbestminbestvabjvbbjjj1untiloutputbestここのという判定条件については、次のような実施例も考えられる。best∞temp∞i1repeatbestminbestvaaivbaitempmintempvbaiii1untilj1repeatbestminbestvabjvbbjjj1untiloutputbestbest2≧minbest2tempなので、こうするとループの脱出が早まる。ここでは、2つのループの打ち切り条件が異なることに留意されたい。あるいは、以下のような例もありえる。best∞i1repeatbestminbestvaaivbaiii1untilj1repeatbestminbestvabjvbbjjj1untiloutputbestここでも、2つのループの打ち切り条件が異なっている。また、上記実施例では、行列A、B、Cはどれもn×nの正方行列であると想定したが、これには限定されず、通常の行列の掛け算と同様に、行列Aがm×k、行列Bがk×nであるとしてよい。この結果、行列Cは、m×nとなる。以上、特定のハードウェアおよびソフトウェアのプラットフォーム上で実施するものとして本発明を説明してきたが、本発明は示されている例に限定されず、任意のコンピュータプラットフォーム上で実施可能である。