JP3659307B2

JP3659307B2 - ベクトル・コンピュータにおける演算方法、及び記録媒体

Info

Publication number: JP3659307B2
Application number: JP17131599A
Authority: JP
Inventors: 巧鷲尾; 克維久保
Original assignee: NEC Corp
Current assignee: NEC Corp
Priority date: 1998-06-18
Filing date: 1999-06-17
Publication date: 2005-06-15
Anticipated expiration: 2019-06-17
Also published as: JP2000029864A; DE19827238A1; DE19827238B4; US6446105B1

Description

【０００１】
【発明の属する技術分野】
本発明は、ベクトル・コンピュータにおいて行列計算を行うための演算方法に関する。
【０００２】
【従来の技術】
楕円型偏微分方程式を差分法により、以下に示す方程式（１）のように離散化し、前処理付反復法で解く場合、前処理演算として、方程式（１）の近似方程式である、以下に示す方程式（２）、（３）が解かれる場合が多い。
【０００３】
Ａｑ＝ｒ．．．（１）
（△−Ｌ）ｐ＝ｒ．．．（２）及び、
△^−１（△−Ｕ）ｑ＝ｐ．．．（３）
ここで、Ａは、ｍ×ｍ行列、Ｌはｍ×ｍの下三角行列、Ｕはｍ×ｍの上三角行列、△はｍ×ｍの対角行列であり、ｑ、ｒ、及びｐはｍ次元のベクトルである。
【０００４】
上記境界値問題は、技術計算において多くの数値を計算する際に生じる問題であり、例えば立体的な物の温度分布を計算する場合や、所定の幾何学形状における流体の流れを計算する場合に生じる。更に簡単な例としては、Dirichlet境界条件でPoisson方程式を解く場合が挙げられる。
【０００５】
通常、この問題についての数値解法として、差分近似法が用いられる。解が計算される空間は、格子点に含まれ、そこで各格子点に対する解の数値計算が行われる。
【０００６】
各格子点に関する解を与える方程式は、例えばテーラー級数に展開され、そこで微分商が差分係数に置き換えられる。格子点における近似解を求めるためには、その級数展開の剰余項が無視して良い程度のものであれば、隣接する格子点の近似解が分かれば充分である。このような方法は、W.Tbrnigによる「Numerische Mathematik fUr Ingenieure und Physiker」、Volime 2,Springer-Verlag, Berlinに開示されている。
【０００７】
例えば、正方格子を有する２次元空間では、４つの格子点に取り囲まれた１つの格子点について近似解を求めるには、隣接するその４つの格子点の解を知るだけで充分である。
【０００８】
図１は、こうしたケースを表しており、いわゆる５点ステンシル(5-point stencil)である。一方、図２は、３次元空間に関して図１に対応するケースを示した、７点ステンシルである。
【０００９】
図２の７点ステンシルに対して、各格子点に適当な数字を付け、上述したアプローチを適用すると、前記問題は、以下の式（４）のように表される。
【００１０】
a_ijk,1 q_ijk-1 + a_ijk,2 q_ij-1k + a_ijk,3 q_i-1jk + a_ijk,4 q_ijk + a_ijk,5 q_i+1jk + a_ijk,6 q_ij+1k + a_ijk,7 q_i1jk+1 = r_ijk ．．．（４）
各格子点(ijk)は、集合Ｇの中から設定される。Ｇは、以下の式（５）で表される。
【００１１】
G = {(ijk) 1 <= i <= n_x, 1 <= j <= n_y, 1 <= k <= n_z} ．．．（５）
ここで、n_x,n_y,n_zは、それぞれｘ軸、ｙ軸、ｚ軸方向における格子点の数を表す。
【００１２】
格子点の符号を適当に選択すると、行列Ａは、その対角成分以外の要素にゼロでない要素を、常に少しだけ含む。
【００１３】
そのため、上述した式（１）を直接解くことはできず、最初に式（６）に表す関係、及び（７）に示す前処理を使って、式（８）に表す式に変換する。
【００１４】
Ａ＝Ｄ−Ｌ−Ｕ．．．（６）
Ｍ＝Ｐ₁×Ｐ₂ ．．．（７）
（Ｐ₁ ^-1ＡＰ₂ ^-1）（Ｐ₂ｑ）＝Ｐ₁ ^-1ｒ．．．（８）。
【００１５】
このとき、Ｄは、元の行列Ａの対角成分のみからなる対角行列である。
【００１６】
よく使用され、好ましい結果をもたらす前処理は、一般に「incomplete lower upper preconditioner」(ＩＬＵ前処理)と呼ばれる。また、以下の式（９）が成り立つ。
【００１７】
Ｍ＝（△−Ｌ’）△^-1（△−Ｕ’）．．．（９）
従って、△＝Ｄ−対角要素（Ｌ△^-1Ｕ）−α行の合計（対角成分でない要素（Ｌ△^-1Ｕ））であり、ここでαは、０ =< α =< 1 である。また、対角成分でない要素（Ｆ）＝Ｆ−対角成分の要素（Ｆ）である。行の合計とは、対角要素のそれぞれが、同じ行の要素の合計に対応することを意味している。
【００１８】
これはよく知られたGustavson(I. Gustafson)の変形であり、「A Class of First Order Fracturization Methods」,BIT 18(1978), pages 142-156を参照のこと。式（９）の行列Ｕ’とＬ’はＩＬＵアプローチにおいて記述される文字であるが、式（６）に示した行列Ｕ、Ｌと同一である必要はないことを示すために、アポストロフィを付けて区別した。これらの行列Ｕ’とＬ’に、適当な変更を加えることによって、収束性の改善が可能となる。
【００１９】
近似的な、ＳＳＯＲというアプローチによれば、式（９）中の△はＤと等しいとすることができる。また、前処理は、以下の１次方程式（１０）の解に対応する。
【００２０】
Ｍｑ＝ｒ．．．（１０）
ここで、ｒは既知のベクトルである。一般に、ＭがＡに近いほど、反復法の反復回数が少なくなることが知られている。この式に関する概略図が図３に示されている。解は、以下の式（１１）及び（１２）を解くことにより得られる。
【００２１】
（△−Ｌ）ｐ＝ｒ．．．（１１）
△^-1（△−Ｕ）ｑ＝ｐ．．．（１２）
最初の問題が、格子に沿ったものである場合、行列Ｕ及びＬのゼロ以外の要素は、規則的に配置されうる。例えば、各軸方向に、n_x, n_y,及びn_z個の格子点を有し、そのそれぞれの格子点を表すよう符号が付けられている、通常の３次元格子上に表された７点ステンシルの場合、ゼロ以外の要素は、行列Ｌにおいて、対角要素に隣接して配置されている第１の非対角要素（第１の非対角要素は、対角要素の集合に沿って隣接し、かつ対角要素と並行に配置される要素の集合）、n_x番目の非対角要素、及び（n_xn_y）番目の非対角要素の位置に配置される。言い換えると、Ｌ_i,i-1,Ｌ_i,i-nx,及びＬ_{i,i-nx ny}の要素のみが、ゼロ以外の要素である。
【００２２】
上記式（１１）及び（１２）は、対応する前進代入(forward substitution)又は後退代入(reverse substitution)によって解を求めることができる。
【００２３】
ｐ＝△^-1（ｒ−ＬＰ）．．．（１３）
ｑ＝ｐ−△^-1Ｕｑ．．．（１４）。
【００２４】
これらの代入アルゴリズム（１）を実行するための疑似コードの例を図４に示す。
【００２５】
この方法は、超平面(hyperplane)法として知られる方法である。超平面は、互いに独立して計算可能である格子点から構成される。超平面は、その超平面のどのような格子点における解も、同じ平面の他の格子点における解に依存しないという特性を有する。３次元空間に関して、対応する超平面が図５に示されている。以下の式（１５）は、７点ステンシルに関して適用される。
【００２６】

この方法の欠点は、特に、ベクトル・プロセッサにおけるキュー処理やベクトル処理の効果が、図４のアルゴリズム（１）の間接的なアドレッシングによって大幅に低減されてしまうということにある。この方法が、常に一定の間隔でメモリをアクセスすることができれば、より好ましいものとなる。これに関しては、「Computer Architecture」、Computational Science Education Project, Section 3.3.1 の、「Interleaved Memory」を参照されたい。
【００２７】
更に、前記方法の処理は、超平面の大きさに応じて並列処理が行われず、多くの超平面Ｈ（１）が、それだけでは複数のプロセッサに論理的に分散させるには小さすぎるために、分割された複数の並列メモリブロックと、複数のベクトル・プロセッサを有するコンピュータには向いていない。また更に、超平面を扱う既知の方法では、高次元空間又は前述した５点ステンシルや７点ステンシル以外のステンシルへの一般化は知られていない。これらの問題を解決するためのアプローチは、前述の式（１１）及び（１２）を正しく解くものではないが、その近似解を計算するものであり、前処理として同様の収束性向上をもたらす。通常、これにはvon-Neumann級数展開が使用される。
【００２８】
ｐ＝（△−Ｌ）^-1ｒ＝ △^-1（Ｉ−Ｌ△^-1）^-1ｒ＝ △^-1ｒ＋ △^-1Ｌ△^-1ｒ＋（△^-1Ｌ）²△^-1ｒ＋（△^-1Ｌ）³△^-1ｒ＋．．．＋（△^-1Ｌ）^ntr△^-1ｒ．．．（１６）。
【００２９】
この方法は、近似解を得るには、常に行列とベクトルの積が計算されるだけであるため、実際には十分にベクトル処理可能なものであり、並列に処理することができる。しかし、超平面法と同様の収束性を実現するためには、von-Neumann級数において、通常３桁又は４桁のオーダーのメンバまで計算しなければならず、計算コストは、３倍から４倍にまで膨らみ、ベクトル化及び並列処理化の利点が相殺されてしまう。
【００３０】
超平面を扱う方法を、ベクトル・コンピュータに適するように変換する他のアプローチは、1995年５月の、Takumi Washio、とKen Hayamiによる「Overlapped Multicolor MILU Preconditioning」、J. SCHI Computer, Volume 16, No. 5, ページ636-651に開示されている。
【００３１】
このアプローチでは、３次元空間及び７点ステンシルについて、最初に超平面の分割が行われ、ｉｓｔｒ個の部分集合が形成される。このことは、以下の式（１７）によって表される。
【００３２】

【００３３】
言い換えれば、部分集合Ｃ（ｌ）は、超平面Ｈ（ｌ）、Ｈ（ｌ＋ｉｓｔｒ）、Ｈ（ｌ＋２ｉｓｔｒ）．．．．．の結合である。
【００３４】
カラー(colors)とも呼ばれる、この部分集合のそれぞれは、前進代入においては、Ｃ（１）、Ｃ（２）、．．．Ｃ（ｉｓｔｒ）、Ｃ（１）、Ｃ（２）、．．．Ｃ（ｉｓｔｒ）の順で計算され、後退代入においては、Ｃ（ｉｓｔｒ）、Ｃ（ｉｓｔｒ−１）、．．．Ｃ（１）、Ｃ（ｉｓｔｒ）、Ｃ（ｉｓｔｒ−１）、．．．Ｃ（１）の順で計算される。この処理順序の例が、図６に示されている。未知の開始値は、それぞれゼロにセットされる。この方法では、各カラーは、少なくとも２回計算され、Ｈ（ｌ）に対するＨ（max(ｌ-istr,1)）の依存度が考慮される。
【００３５】
図７には、この方法のアルゴリズム（２）を実現する疑似コードの例が示されており、図８には、対応するフローチャートが示されている。
【００３６】
ｉｓｔｒの選択を適切に行うことによって、各カラーが、ベクトル処理に充分な大きさの数の格子点を有するようにすることができる。前述の７点ステンシルにおいて、例えば、n_x＝n_yでかつ、n_zが任意である場合、ｉｓｔｒ＝n_x - 1 とすれば、Ｃ（ｌ）内の格子点を計算するための全てのメモリ・アクセスが、一定の間隔ｉｓｔｒで行われ、間接アドレッシングが回避される。一方、(n_xn_yn_z)/istrは、ベクトル長として十分な大きさである。
【００３７】
例えば、ｘ軸方向がn_x = 10で、ｙ軸方向がn_y = 11で、ｚ軸方向がn_z = 10で、ｉｓｔｒが前述のとおり９と求められた、格子を有する７点ステンシルの場合、カラーへの分割が行われた後、図６に示すような計算が行われる。この図から、格子点ｉを計算するには更に６つの格子点、i+1、i+10、i+110、i-1、i-10、i-110が必要とされ、これらは同じカラーに属さないことが明らかである。従って、あるカラーの、ある格子点の計算は、同じカラーの他の格子点の計算には依存しないことが分かる。更に、間接アドレッシングは、ｉを知ることによって直接把握できるので必要がない。そして、参照すべき格子点、i+1、i+10、i+110、i-1、i-10、i-110のデータが格納されているメモリ領域を事前にロードしておくことができるので、ベクトル・プロセッサの高速処理が可能である。
【００３８】
この複数のカラーに対する処理技法を更に改良したものが、前述の文献に記載されており、それは「Overlapping Multicolor」技法と呼ばれている。ここでは、前進代入に関して、カラーがＣ（ｉｓｔｒ−ｗ）、Ｃ（ｉｓｔｒ−ｗ−１）、．．．Ｃ（ｉｓｔｒ）、Ｃ（１）、Ｃ（２）、．．．Ｃ（ｉｓｔｒ）の順で行われるという点で、複数のカラーを近似する技法の改良が行われている。オーバーラッピング、即ち、追加的に計算されるカラーの数はｗになる。対応するシステムが後退代入に適用される。この、オーバーラッピング多重カラー技法は、単一の多重カラー技法に較べて、優れた収束性を示し、一方で、ベクトル長とメモリ・アクセスに関しては同様の利点を有する。
【００３９】
【発明が解決しようとする課題】
しかし、多重カラー技法、及びオーバーラッピングを備えた多重カラー技法においては、カラーの数を求めることと、前述の７点又は５点ステンシル以外のケースに当該技法を適用することが困難である。最初の近似に関する開始値が、解から離れすぎていると、カラーの数が小さくなりすぎ、結果的に計算コストを増大させてしまう。ベクトル・コンピュータのキューのベクトル長が短くなりすぎると、カラーの数が大きくなりすぎ、結果的に計算スピードが低下してしまう。
【００４０】
前述の７点又は５点ステンシル以外のステンシルについてｉｓｔｒを求めることは、カラーの個々の格子点が他の格子点に依存しないという条件を維持しなければならないため、単純には実現できない。
【００４１】
【課題を解決するための手段】
本発明による方法は、特に、規則的な正方格子に関する楕円型偏微分方程式を解くのに適している。
【００４２】
本発明の目的は、所与の長さのキューを備えたベクトル・コンピュータにおいて、前述の式（１）の解を求める方法を提供することにある。
【００４３】
本発明の更なる目的は、特に、キューの長さやメモリ・バンクの数といった、ベクトル・コンピュータの実際の構成に適応して、他の方法より高速に演算を行う方法を提供することにある。
【００４４】
本発明の更なる目的は、計算コストを追加せずに、５点ステンシル又は７点ステンシル以外のステンシルにおける依存度を考慮する方法を提供することにある。
【００４５】
本発明の更なる目的は、好適に並列処理され、即ち、複数のプロセッサ上で並列に処理が行われ、特に、分散メモリを有するコンピュータにおいても好適に処理される方法を提供することにある。
【００４６】
【発明の実施の形態】
本発明を簡単に理解するため、まず、ベクトル・コンピュータの原理をいくつか説明する。
【００４７】
既存のコンピュータは、１つ、又はそれ以上のプロセッサ、１つ、又はそれ以上のメモリ、入／出力デバイス、及びこれらの間の通信チャネルを有している。
【００４８】
プロセッサは、コマンドを実行するために、そのコマンドをメモリにロードし、デコードし、実行する。このような処理サイクルにかかる時間は、そのコマンドの複雑性及びプロセッサの内部構造に依存するが、メモリ及び通信チャネルの構成にも依存する。
【００４９】
コンピュータ・プログラムは、１つ若しくは２，３の単純なコマンドの繰り返しであるループ処理を実行することが多い。そのようなループ処理の良い例はベクトルの和の計算及び／又はベクトルのクロス積の計算である。
【００５０】
ここでは、２進演算プロセッサにおいて２つの数の積を単純に計算する場合であっても、その処理には通常、一連の個別ステップが含まれることに注意すべきである。このことは、例えば、ビットの列をそれぞれ左に１ビット移動させる（即ち２倍する）場合や、２進数の足し算を行うような単純な場合のような、最も単純な場合にも当てはまる。
【００５１】
このような処理を行うアルゴリズムが、単純な既知のスカラ・コンピュータによって処理される場合、処理時間のほとんどは、同じコマンドのロードとデコード、及び中間の処理結果の記憶を繰り返すことに費やされる。
【００５２】
この問題を解決するために、ベクトル・コンピュータが使用される。
【００５３】
ベクトル・コンピュータの処理サイクルのタイミングが、図９に概略的に示されている。
【００５４】
図９に示す例の、最初のサイクルでは、第１のコマンドがコマンド・メモリから取り出される。第２のサイクルでは、このコマンドがプロセッサに渡され、オペランドの準備がされる。このサイクルにおいて更に、第２のコマンドがコマンド・メモリから取り出される。第３のサイクルからｍ番目のサイクルまで、第１のコマンドが実行される。また、第４のサイクルから、（ｍ＋１）番目のサイクルまで、第２のコマンド他が実行される。ｍ番目のサイクルが終了すると、プロセッサは、各サイクルの結果を出力する。この処理では、ｎ個のコマンドの処理に全部でｍ＋ｎ−１回のサイクルを必要とする。一方、この処理をスカラ・コンピュータで行うと、（ｎｍ）サイクルを要することになる。
【００５５】
言い換えれば、ベクトル・コンピュータは複数のキューを使用し、データやコマンドを、ローディング・キューを介してベクトル・レジスタに転送する。その後、算術演算キューにおいて計算が行われ、演算結果はメモリ・キューを介してメモリに記憶される。算術演算キューは、多機能キューであり、浮動小数点演算といった個別の複数ステップを並列で実行することができる。
【００５６】
１秒あたりの浮動小数点演算回数(FLOPS)で表される、このベクトル・コンピュータの性能、即ち計算能力は、いわゆるベクトル長に密接に依存している。図９から分かるように、コマンド・スタックがベクトル処理される場合、第１のコマンドの処理はｍサイクル必要であり、それ以降の各コマンドは、１サイクルだけの追加で処理される。この処理方法は、ｎ番目のコマンドが、先行するｍ−１個のコマンドのうちの１つの処理結果を知ることなく実行できることを前提としている。これは、上述のベクトルの組み合わせにおいては典型的なケースだからである。
【００５７】
従って、ベクトル長は、ベクトル化して実行可能なコマンドの数を表す数字として認識される。従来より、ベクトル・コンピュータに適したベクトル長は、単純ループによって決定されうる。アルゴリズム（３）は、そのループを実行するためのものであり、対応する疑似コードが、図１０に示されている。
【００５８】
図１１は、ベクトル長がコンピュータの性能に与える影響について、例えばＮＥＣのＳＸ４のようなコンピュータを例にして示した図である。最初は、ベクトル長の増加に伴って明らかにコンピュータの性能が向上しているが、ベクトル長ｎが４００を越えると、コンピュータの性能にほとんど変化がない。
【００５９】
この説明において、最小ベクトル長は、コンピュータの性能にほとんど影響を及ぼさないようになるベクトル長を表し、この例では約４００である。
【００６０】
最小ベクトル長は、ベクトル・プロセッサの使用の形態に関わらず、そのベクトル・プロセッサに依存する。従って、例えばFujitsu-VPXは、ベクトル処理を効果的に実現するために、１０００のベクトル長を要求する。一方、Cray EL98は、かなり短いベクトル長でも効果的に機能する。
【００６１】
それぞれの最小ベクトル長より大きなベクトル長が、コンピュータの性能にほとんど影響を与えないという事実は、ある理由によるものである。実際に、図９に示すような、充分なオーバーラッピングが行われる状況では、例えば、キュー内のステップ数が、スタック内のコマンド数と等しくなる場合には、それ以上ベクトル長を大きくしても僅かな効果しか得られない。
【００６２】
コンピュータの性能を更に劣化させる要因は、メモリの数、構成、及び構造である。ある記憶位置にアクセスするため、それに関連付けられたメモリ・バンクが常に事前ロードされ、アクセスが生じた場合だけ参照される。次にメモリ・バンクは、再び戻される（ダウンシフト）。こうした事前ロードとダウンシフトの時間は、実際のアクセス時間の何倍にもなることがある。従って、従来では、連続してアクセスされるビットが、いくつかのメモリ・バンクに記憶されており、次のメモリ・バンクは常に事前ロードされていた。こうして、メモリ・アクセスに必要な時間がかなり短縮された。
【００６３】
しかし、これは、どのメモリ・バンクが連続してアクセスされるかが常に分かっている場合を想定している。間接的なアドレッシングを用いる、前述した超平面を扱う方法においては、次にアクセスされるアドレスを常に最初に計算しなければならないので、これを想定することは不可能である。従って、間接的なアドレッシングへの適用は避けなければならない。
【００６４】
一方、既知で一定の間隔によりメモリがアクセスされる方法は、ベクトル・コンピュータへの適用に適している。既知で一定の前記間隔は、例えば、連続してアクセスされる記憶位置の間隔が、常に既知で一定であるような場合である。
【００６５】
直接的アドレッシングでないアクセスが同じメモリ・バンクに対して短時間で連続して行われる場合にも、同様のメモリ・アクセスの問題が生じる。これは、メモリ・バンクが常にダウンシフトされ、その後、再開されたアクセスの前に再び事前ロードされるので、コンピュータの性能をも低下させる。
【００６６】
従って、ベクトル化処理されるループは、同一メモリ・バンクへのアクセスの間隔が十分長いものであるよう設計されることが好ましい。
【００６７】
コマンドの処理が複数の独立した並列処理ベクトル・プロセッサに亘って分散されると、コンピュータの性能は一層向上する。このための条件は、複数の独立した部分タスクを形成できることである。
【００６８】
前述したように、本発明は、ベクトル・コンピュータ上で多重カラーを扱う方法を実行し、式（１）の解を求める方法に関するものである。
【００６９】
コンピュータの性能に悪影響を及ぼす、前述した要因を回避するために、本発明は、一定のメモリ・アクセス間隔ｉｓｔｒを以下の規則に従って決定する。
【００７０】
ａ）すべてのｋ(k = 1 ...nzr)に対して、mod(izr_k, istr) ≠ 0であること。ここで、izr_k = {j:L_i, _i-j ≠ 0 又は U_i, _i+j ≠ 0}、nzrは、L_i, _i-j ≠ 0 又は U_i, _i+j ≠ 0である全てのｊの数である。
【００７１】
ｂ）istrとnbの最大公約数(gcd)が１であること。ここで、nbはメモリ・バンクの数である。
【００７２】
ｃ）istrはできるだけ小さな数であること。その結果、実際のベクトル長が、そのベクトル・プロセッサに要求される最小ベクトル長より大きいものとなる。
【００７３】
また、istrの実際の値によって、以下のような値ｘを求め、その値がなるべく大きくなるようにする。
【００７４】
x = min{int(istr/mod (izr_k, istr)) k = 1,...nzr} ．．．（１８）。
【００７５】
規則ａ）は、全ての計算処理が、あるカラーの中では独立して実行されることを保証している。このことによって、ベクトル処理が、１つのカラーにおいて可能になる。
【００７６】
これは、ｉｓｔｒによって割り切れるizr_kがないことを仮定している。この規則を用いることによって、５点ステンシルや７点ステンシルを扱うことができるだけでなく、あらゆる依存度を考慮に入れることができる。
【００７７】
規則ｂ）は、メモリの競合が無いことを保証する。即ち、同じメモリ・バンクに対するアクセスを連続してすることができず、かつ、全てのメモリ・バンクが使用されている。
【００７８】
規則ｃ）は、ベクトル長が、ベクトル・コンピュータを最適動作させるために必要な程度に大きいことを保証する。これに関しては、ベクトル長を大きくする程、この方法による収束性が低下することに注意すべきである。これは、前進代入及び後退代入の開始時点で、すべてのｉに対して要素ｐ_i又はｑ_iの値がそれぞれゼロにセットされる（又は別の適切な開始値がセットされる）ことによるものである。ループの最初の１回が終了すると、要素ｐ_i+istr, ｐ_i+2istr, ...ｑ_mは、それぞれ第１の近似値に対応した新しい値がセットされる。次に、他の要素ｐ_i又はｑ_iの全てについて計算を更に行っていく。更なる計算に利用できる第１の近似値が迅速に求められるほど、収束性がより効果的なものになる。
【００７９】
そこで、良好な収束性を得るために、式（１８）が提供される。これは、前述のvon-Neumann級数展開に関連するものであり、本発明による方法の収束性が、式（１６）において、ntr =（nd - 1）xである限り、式（１６）による級数展開と同程度であることを保証する。
【００８０】
ｉｓｔｒを計算するために好適なアルゴリズム（４）を実現する疑似コードの例を、図１２に示す。
【００８１】
また、以下の説明に関連して、前記アルゴリズム（４）のフローチャートを図１３に示す。
【００８２】
ここでは、lvcが最小ベクトル長、即ち、ベクトル・プロセッサによって要求される最小のベクトル長であり、gcd(i, nb)が、ｉとｎｂの最大公約数である。
【００８３】
また、最小ベクトル長は、厳密な境界を有していないので、特に最後の規則ｃ）と式（１８）の適用については、ファジー理論と関連付けて考えることができる。ファジー理論を取り入れると、前述のアルゴリズム（４）では、例えば、ｉによってループするループは、ｉがistrmaxになるまでではなく、istrmaxの何倍かになるまでループするようにでき、あるｉ(i > istrmax)における解は、ｉがistrmaxより大きいという点で、他の解より好ましくないものとすることができる。
以下では、前述した本発明の適用について、np=2のベクトル・プロセッサを備えた並列ベクトル・コンピュータに関連して説明する。当業者には、プロセッサが３つ以上のベクトル・コンピュータにも、本発明を同様に適用できることは明らかである。
【００８４】
図１４には、２つのプロセッサによって、図２に示す７点ステンシルを処理する順序が示されている。ここで、n_x = 10, n_y = 11, istr = 9, 及びn_zは任意である。
【００８５】
最初に、第２のプロセッサと独立して動作する第１のプロセッサは、カラーＣ（１）及びＣ（２）について第１の近似解を計算する。同時に、第２のプロセッサは、第１のプロセッサとは独立して、カラーＣ（５）及びＣ（６）についての第１の近似解を計算する。このために、カラーＣ（８）、Ｃ（９）又はＣ（３）、Ｃ（４）の格子点の第１の近似解の値は、それぞれゼロにセットされる。
【００８６】
カラーＣ（２）又はＣ（６）の計算の後、同期がとられる。即ち、次のカラーＣ（３）又はＣ（７）の計算は、両方のプロセッサがそれぞれの先行するカラーの計算を終えた後でなければ、始めることができない。カラーＣ（３）又はＣ（７）の値は、カラーＣ（５）又はＣ（１）の値の計算に用いられるので、ここでも同期が必要である。カラーＣ（５）の計算がまだ終了していない間に、カラーＣ（３）の値を変更することは、計算の不確実性を招く。それは、カラーＣ（３）の値が、カラーＣ（５）の計算に必要だからである。結果的に、カラーＣ（５）の結果は、それぞれのプロセッサのランダムな動作条件に、多かれ少なかれ依存することになる。
【００８７】
同期をとった後で、カラーＣ（３）及びＣ（４）は、第１のプロセッサによって計算され、カラーＣ（７）、Ｃ（８）、及びＣ（９）は他のプロセッサによって計算される。次に、これらのプロセッサは再び同期をとり、処理は新たなパスに進み、カラーＣ（１）又はＣ（５）の計算が開始される。また、カラーＣ（８）、Ｃ（９）又はＣ（３）、Ｃ（４）の第１の近似解は、以前のパスで計算されている。このループは、条件id = ndになるまで、アルゴリズム（２）によって繰り返される。即ち、近似解が十分良好なものになることが保証されるまで繰り返される。通常、満足な近似解を得るには、２つのパス(nd - 2)で充分である。
【００８８】
並列ベクトル・コンピュータで上記処理を行うための好適アルゴリズム（５）を実行するための疑似コードの例が図１５、図１６に示されている。図１５に示す疑似コードと図１６に示す疑似コードとは連続しているものであり、これら２つで１つのアルゴリズム（５）を示していることに注意すべきである。
【００８９】
また、以降の説明に関連して、このアルゴリズム（５）に対応するフローチャートが、図１７に示されている。
【００９０】
２つ以上のプロセッサがあれば、ループ処理は、ｉａを使って、０からint(istr/np) - 1 まで実行される。またここで、npはプロセッサの数を表している。
【００９１】
【発明の効果】
以上に説明したように、本発明による方法によって、様々なキューの長さ、及びメモリ・バンクの数を有するようなベクトル・コンピュータに適用させることができ、結果的に既存の方法より高速な動作を実現できる。また、いかなる追加コストも必要とせずに、５点ステンシルや７点ステンシルにおける依存度以外の依存度を更に考慮することができる。
【００９２】
例えば、図１８及び図１９にそれぞれ示す、１９及び２７点ステンシルについても計算が可能である。
【００９３】
更に、本発明による方法は、複数の独立プロセッサを用いた並列処理に非常に適している。特に、分散メモリを備え、各プロセッサがそのメモリにアロケートされているメモリ領域を有し、アクセスの競合が回避されるようなコンピュータ・システムに適している。
【図面の簡単な説明】
【図１】５点ステンシルの２次元格子を示す図である。
【図２】７点ステンシルの３次元格子を示す図である。
【図３】７点ステンシルに関する行列Ｍの概略を示す図である。
【図４】１次方程式の解を前進代入及び後退代入を用いて解くためのアルゴリズム（１）を実現する疑似コードの例を示す図である。
【図５】３次元空間における超平面の位置の概略を示す図である。
【図６】 n_x = 10, n_y = 11, istr = 9で、n_zが任意である場合の計算順序を概略示す図である。
【図７】アルゴリズム（２）を実現する疑似コードの例を示す図である。
【図８】アルゴリズム（２）についてのフローチャートである。
【図９】ベクトル・コンピュータの処理タイミングを示す図である。
【図１０】ベクトル長を決定するアルゴリズム（３）を実現する疑似コードの例を示す図である。
【図１１】既存のベクトル・コンピュータのベクトル長と処理能力の関係を表すグラフである。
【図１２】ｉｓｔｒを計算するために好適なアルゴリズム（４）を実現する疑似コードの例を示す図である。
【図１３】図１２のアルゴリズム（４）についてのフローチャートである。
【図１４】 n_x = 10, n_y = 11, istr = 9で、n_zが任意であり、２つのプロセッサ上に分散された処理が並列に行われる場合の計算順序を概略示す図である。
【図１５】アルゴリズム（５）を実現するための疑似コードの例の一部を示す図である。
【図１６】アルゴリズム（５）を実現するための疑似コードの例の一部を示す図である。
【図１７】図１５及び図１６のアルゴリズム（５）についてのフローチャートである。
【図１８】３次元空間における、１９点のステンシルの配置を示す図である。
【図１９】３次元空間における、２７点のステンシルの配置を示す図である。

Claims

少なくとも１つのベクトル・プロセッサとメモリ・バンクを有するベクトル・コンピュータで、下三角行列と上三角行列を含んだ行列式の近似解を求めるための、前進代入及び後退代入を行う多重カラー技法に基づく計算方法における、前記代入を行うために、メモリをアクセスする一定の間隔ｉｓｔｒを決定する方法において、
前記代入の前に、所定の基準により、前記ｉｓｔｒを決定する方法であって、前記ｉｓｔｒが以下の規則に従って決定されることを特徴とするメモリ・アクセス間隔決定方法、
ａ）前記下三角行列及び上三角行列において対角成分の並びと平行な、非対角成分のそれぞれの並びの中で、少なくとも１つの成分がゼロでない各並びに対して、その対角成分の並びから見て、その非対角成分の並びがいくつ目の平行な要素かを表す数ｉｚｒをｉｓｔｒで割った場合に、その剰余がゼロにならないｉｓｔｒを選択し、
ｂ）ｉｓｔｒとメモリ・バンクの数ｎｂの最大公約数（ＧＣＤ）を計算し、その結果が１となるｉｓｔｒを選択し、
ｃ）可能な数の候補の中で最も小さいｉｓｔｒを求め、
ｄ）各ｉｚｒに関して、ｉｚｒをｉｓｔｒで割った前記剰余で更にｉｓｔｒを割り、その整数部の値が最小であるものをｘとする場合、前記ｘの値が、前記ベクトル・コンピュータで並列に動作するベクトル・プロセッサの数の少なくとも２倍以上となるように、ｉｓｔｒの値を決定する。
請求項１において、前記行列式が、
（△−Ｌ）△^−１（△―Ｕ）ｑ＝ｒ
であり、ここで、△はｍ×ｍの対角行列、Ｌはｍ×ｍの下三角行列、Ｕはｍ×ｍの上三角行列、ｑ及びｒはｍ次元のベクトルであることを特徴とするメモリ・アクセス間隔決定方法。
少なくとも１つのベクトル・プロセッサとメモリ・バンクを有するベクトル・コンピュータで、下三角行列と上三角行列を含んだ行列式の近似解を求めるための、前進代入及び後退代入を行う多重カラー技法に基づく計算方法における、前記代入を行うために、メモリをアクセスする一定の間隔ｉｓｔｒを決定するプログラムを記録したコンピュータ読み取り可能な記録媒体であって、前記プログラムは、
前記代入の前に、所定の基準により、前記ｉｓｔｒを決定するものであって、前記ｉｓｔｒが以下の規則に従って決定されることを特徴とするプログラムを記録したコンピュータ読み取り可能な記録媒体、
ａ）前記下三角行列及び上三角行列において対角成分の並びと平行な、非対角成分のそれぞれの並びの中で、少なくとも１つの成分がゼロでない各並びに対して、その対角成分の並びから見て、その非対角成分の並びがいくつ目の平行な要素かを表す数ｉｚｒをｉｓｔｒで割った場合に、その剰余がゼロにならないｉｓｔｒを選択し、
ｂ）ｉｓｔｒとメモリ・バンクの数ｎｂの最大公約数（ＧＣＤ）を計算し、その結果が１となるｉｓｔｒを選択し、
ｃ）可能な数の候補の中で最も小さいｉｓｔｒを求め、
ｄ）各ｉｚｒに関して、ｉｚｒをｉｓｔｒで割った前記剰余で更にｉｓｔｒを割り、その整数部の値が最小であるものをｘとする場合、前記ｘの値が、前記ベクトル・コンピュータで並列に動作するベクトル・プロセッサの数の少なくとも２倍以上となるように、ｉｓｔｒの値を決定する。