JP2000029864A - ベクトル・コンピュ―タにおける演算方法、及び記録媒体 - Google Patents
ベクトル・コンピュ―タにおける演算方法、及び記録媒体Info
- Publication number
- JP2000029864A JP2000029864A JP11171315A JP17131599A JP2000029864A JP 2000029864 A JP2000029864 A JP 2000029864A JP 11171315 A JP11171315 A JP 11171315A JP 17131599 A JP17131599 A JP 17131599A JP 2000029864 A JP2000029864 A JP 2000029864A
- Authority
- JP
- Japan
- Prior art keywords
- istr
- vector
- computer
- triangular matrix
- substitution
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Complex Calculations (AREA)
Abstract
(57)【要約】
【課題】 ベクトル・コンピュータにおいて行列計算を
行うための高速な演算方法を提供する。 【解決手段】 本発明の方法は、ベクトル・コンピュー
タにおける効率的な動作のために、既知で一定の間隔に
よってメモリがアクセスされるよう制御する。これによ
って、本発明による方法を、様々なキューの長さ、及び
メモリ・バンクの数を有するようなベクトル・コンピュ
ータに適用させることができ、結果的に、既存の方法に
較べて高速な動作を実現できる。また、いかなる追加コ
ストも必要とせずに、５点ステンシルや７点ステンシル
における依存度以外の依存度を更に考慮することができ
る。
行うための高速な演算方法を提供する。 【解決手段】 本発明の方法は、ベクトル・コンピュー
タにおける効率的な動作のために、既知で一定の間隔に
よってメモリがアクセスされるよう制御する。これによ
って、本発明による方法を、様々なキューの長さ、及び
メモリ・バンクの数を有するようなベクトル・コンピュ
ータに適用させることができ、結果的に、既存の方法に
較べて高速な動作を実現できる。また、いかなる追加コ
ストも必要とせずに、５点ステンシルや７点ステンシル
における依存度以外の依存度を更に考慮することができ
る。
Description
【０００１】
【発明の属する技術分野】本発明は、ベクトル・コンピ
ュータにおいて行列計算を行うための演算方法に関す
る。
ュータにおいて行列計算を行うための演算方法に関す
る。
【０００２】
【従来の技術】楕円型偏微分方程式を差分法により、以
下に示す方程式（１）のように離散化し、前処理付反復
法で解く場合、前処理演算として、方程式（１）の近似
方程式である、以下に示す方程式（２）、（３）が解か
れる場合が多い。
下に示す方程式（１）のように離散化し、前処理付反復
法で解く場合、前処理演算として、方程式（１）の近似
方程式である、以下に示す方程式（２）、（３）が解か
れる場合が多い。
【０００３】Ａｑ＝ｒ ．．．（１） （△−Ｌ）ｐ＝ｒ ．．．（２）及び、 △−１（△−Ｕ）ｑ＝ｐ ．．．（３） ここで、Ａは、ｍ×ｍ行列、Ｌはｍ×ｍの下三角行列、
Ｕはｍ×ｍの上三角行列、△はｍ×ｍの対角行列であ
り、ｑ、ｒ、及びｐはｍ次元のベクトルである。
Ｕはｍ×ｍの上三角行列、△はｍ×ｍの対角行列であ
り、ｑ、ｒ、及びｐはｍ次元のベクトルである。
【０００４】上記境界値問題は、技術計算において多く
の数値を計算する際に生じる問題であり、例えば立体的
な物の温度分布を計算する場合や、所定の幾何学形状に
おける流体の流れを計算する場合に生じる。更に簡単な
例としては、Dirichlet境界条件でPoisson方程式を解く
場合が挙げられる。
の数値を計算する際に生じる問題であり、例えば立体的
な物の温度分布を計算する場合や、所定の幾何学形状に
おける流体の流れを計算する場合に生じる。更に簡単な
例としては、Dirichlet境界条件でPoisson方程式を解く
場合が挙げられる。
【０００５】通常、この問題についての数値解法とし
て、差分近似法が用いられる。解が計算される空間は、
格子点に含まれ、そこで各格子点に対する解の数値計算
が行われる。
て、差分近似法が用いられる。解が計算される空間は、
格子点に含まれ、そこで各格子点に対する解の数値計算
が行われる。
【０００６】各格子点に関する解を与える方程式は、例
えばテーラー級数に展開され、そこで微分商が差分係数
に置き換えられる。格子点における近似解を求めるため
には、その級数展開の剰余項が無視して良い程度のもの
であれば、隣接する格子点の近似解が分かれば充分であ
る。このような方法は、W.Tbrnigによる「NumerischeMa
thematik fUr Ingenieure und Physiker」、Volime 2,S
pringer-Verlag, Berlinに開示されている。
えばテーラー級数に展開され、そこで微分商が差分係数
に置き換えられる。格子点における近似解を求めるため
には、その級数展開の剰余項が無視して良い程度のもの
であれば、隣接する格子点の近似解が分かれば充分であ
る。このような方法は、W.Tbrnigによる「NumerischeMa
thematik fUr Ingenieure und Physiker」、Volime 2,S
pringer-Verlag, Berlinに開示されている。
【０００７】例えば、正方格子を有する２次元空間で
は、４つの格子点に取り囲まれた１つの格子点について
近似解を求めるには、隣接するその４つの格子点の解を
知るだけで充分である。
は、４つの格子点に取り囲まれた１つの格子点について
近似解を求めるには、隣接するその４つの格子点の解を
知るだけで充分である。
【０００８】図１は、こうしたケースを表しており、い
わゆる５点ステンシル(5-point stencil)である。一
方、図２は、３次元空間に関して図１に対応するケース
を示した、７点ステンシルである。
わゆる５点ステンシル(5-point stencil)である。一
方、図２は、３次元空間に関して図１に対応するケース
を示した、７点ステンシルである。
【０００９】図２の７点ステンシルに対して、各格子点
に適当な数字を付け、上述したアプローチを適用する
と、前記問題は、以下の式（４）のように表される。
に適当な数字を付け、上述したアプローチを適用する
と、前記問題は、以下の式（４）のように表される。
【００１０】 aijk,1 qijk-1 + aijk,2 qij-1k + aijk,3 qi-1jk + aijk,4 qijk + aijk,5 q i+1jk + aijk,6 qij+1k + aijk,7 qi1jk+1 = rijk ．．．（４） 各格子点(ijk)は、集合Ｇの中から設定される。Ｇは、
以下の式（５）で表される。
以下の式（５）で表される。
【００１１】 G = {(ijk) 1 <= i <= nx, 1 <= j <= ny, 1 <= k <= nz} ．．．（５） ここで、nx,ny,nzは、それぞれｘ軸、ｙ軸、ｚ軸方向に
おける格子点の数を表す。
おける格子点の数を表す。
【００１２】格子点の符号を適当に選択すると、行列Ａ
は、その対角成分以外の要素にゼロでない要素を、常に
少しだけ含む。
は、その対角成分以外の要素にゼロでない要素を、常に
少しだけ含む。
【００１３】そのため、上述した式（１）を直接解くこ
とはできず、最初に式（６）に表す関係、及び（７）に
示す前処理を使って、式（８）に表す式に変換する。
とはできず、最初に式（６）に表す関係、及び（７）に
示す前処理を使って、式（８）に表す式に変換する。
【００１４】 Ａ＝Ｄ−Ｌ−Ｕ ．．．（６） Ｍ＝Ｐ1×Ｐ2 ．．．（７） （Ｐ1 -1ＡＰ2 -1）（Ｐ2ｑ）＝Ｐ1 -1ｒ ．．．
（８）。
（８）。
【００１５】このとき、Ｄは、元の行列Ａの対角成分の
みからなる対角行列である。
みからなる対角行列である。
【００１６】よく使用され、好ましい結果をもたらす前
処理は、一般に「incomplete lowerupper precondition
er」(ＩＬＵ前処理)と呼ばれる。また、以下の式（９）
が成り立つ。
処理は、一般に「incomplete lowerupper precondition
er」(ＩＬＵ前処理)と呼ばれる。また、以下の式（９）
が成り立つ。
【００１７】 Ｍ＝（△−Ｌ’）△-1（△−Ｕ’） ．．．（９） 従って、△＝Ｄ−対角要素（Ｌ△-1Ｕ）−α行の合計
（対角成分でない要素（Ｌ△-1Ｕ））であり、ここでα
は、０ =< α =< 1 である。また、対角成分でない要素
（Ｆ）＝Ｆ−対角成分の要素（Ｆ）である。行の合計と
は、対角要素のそれぞれが、同じ行の要素の合計に対応
することを意味している。
（対角成分でない要素（Ｌ△-1Ｕ））であり、ここでα
は、０ =< α =< 1 である。また、対角成分でない要素
（Ｆ）＝Ｆ−対角成分の要素（Ｆ）である。行の合計と
は、対角要素のそれぞれが、同じ行の要素の合計に対応
することを意味している。
【００１８】これはよく知られたGustavson(I. Gustafs
on)の変形であり、「A Class of First Order Fracturi
zation Methods」,BIT 18(1978), pages 142-156を参照
のこと。式（９）の行列Ｕ’とＬ’はＩＬＵアプローチ
において記述される文字であるが、式（６）に示した行
列Ｕ、Ｌと同一である必要はないことを示すために、ア
ポストロフィを付けて区別した。これらの行列Ｕ’と
Ｌ’に、適当な変更を加えることによって、収束性の改
善が可能となる。
on)の変形であり、「A Class of First Order Fracturi
zation Methods」,BIT 18(1978), pages 142-156を参照
のこと。式（９）の行列Ｕ’とＬ’はＩＬＵアプローチ
において記述される文字であるが、式（６）に示した行
列Ｕ、Ｌと同一である必要はないことを示すために、ア
ポストロフィを付けて区別した。これらの行列Ｕ’と
Ｌ’に、適当な変更を加えることによって、収束性の改
善が可能となる。
【００１９】近似的な、ＳＳＯＲというアプローチによ
れば、式（９）中の△はＤと等しいとすることができ
る。また、前処理は、以下の１次方程式（１０）の解に
対応する。
れば、式（９）中の△はＤと等しいとすることができ
る。また、前処理は、以下の１次方程式（１０）の解に
対応する。
【００２０】Ｍｑ＝ｒ ．．．（１０） ここで、ｒは既知のベクトルである。一般に、ＭがＡに
近いほど、反復法の反復回数が少なくなることが知られ
ている。この式に関する概略図が図３に示されている。
解は、以下の式（１１）及び（１２）を解くことにより
得られる。
近いほど、反復法の反復回数が少なくなることが知られ
ている。この式に関する概略図が図３に示されている。
解は、以下の式（１１）及び（１２）を解くことにより
得られる。
【００２１】（△−Ｌ）ｐ＝ｒ ．．．（１１） △-1（△−Ｕ）ｑ＝ｐ ．．．（１２） 最初の問題が、格子に沿ったものである場合、行列Ｕ及
びＬのゼロ以外の要素は、規則的に配置されうる。例え
ば、各軸方向に、nx, ny,及びnz個の格子点を有し、そ
のそれぞれの格子点を表すよう符号が付けられている、
通常の３次元格子上に表された７点ステンシルの場合、
ゼロ以外の要素は、行列Ｌにおいて、対角要素に隣接し
て配置されている第１の非対角要素（第１の非対角要素
は、対角要素の集合に沿って隣接し、かつ対角要素と並
行に配置される要素の集合）、nx番目の非対角要素、及
び（nxny）番目の非対角要素の位置に配置される。言い
換えると、Ｌi,i-1,Ｌi,i-nx,及びＬi,i-nx nyの要素の
みが、ゼロ以外の要素である。
びＬのゼロ以外の要素は、規則的に配置されうる。例え
ば、各軸方向に、nx, ny,及びnz個の格子点を有し、そ
のそれぞれの格子点を表すよう符号が付けられている、
通常の３次元格子上に表された７点ステンシルの場合、
ゼロ以外の要素は、行列Ｌにおいて、対角要素に隣接し
て配置されている第１の非対角要素（第１の非対角要素
は、対角要素の集合に沿って隣接し、かつ対角要素と並
行に配置される要素の集合）、nx番目の非対角要素、及
び（nxny）番目の非対角要素の位置に配置される。言い
換えると、Ｌi,i-1,Ｌi,i-nx,及びＬi,i-nx nyの要素の
みが、ゼロ以外の要素である。
【００２２】上記式（１１）及び（１２）は、対応する
前進代入(forward substitution)又は後退代入(reverse
substitution)によって解を求めることができる。
前進代入(forward substitution)又は後退代入(reverse
substitution)によって解を求めることができる。
【００２３】 ｐ＝△-1（ｒ−ＬＰ） ．．．（１３） ｑ＝ｐ−△-1Ｕｑ ．．．（１４）。
【００２４】これらの代入アルゴリズム（１）を実行す
るための疑似コードの例を図４に示す。
るための疑似コードの例を図４に示す。
【００２５】この方法は、超平面(hyperplane)法として
知られる方法である。超平面は、互いに独立して計算可
能である格子点から構成される。超平面は、その超平面
のどのような格子点における解も、同じ平面の他の格子
点における解に依存しないという特性を有する。３次元
空間に関して、対応する超平面が図５に示されている。
以下の式（１５）は、７点ステンシルに関して適用され
る。
知られる方法である。超平面は、互いに独立して計算可
能である格子点から構成される。超平面は、その超平面
のどのような格子点における解も、同じ平面の他の格子
点における解に依存しないという特性を有する。３次元
空間に関して、対応する超平面が図５に示されている。
以下の式（１５）は、７点ステンシルに関して適用され
る。
【００２６】 Ｈ（ｌ）＝｛i = ix + (iy - 1)nx + (iz - 1)nx ny | ix + iy + iz = ｌ + 2｝ ．．（１５） この方法の欠点は、特に、ベクトル・プロセッサにおけ
るキュー処理やベクトル処理の効果が、図４のアルゴリ
ズム（１）の間接的なアドレッシングによって大幅に低
減されてしまうということにある。この方法が、常に一
定の間隔でメモリをアクセスすることができれば、より
好ましいものとなる。これに関しては、「Computer Arc
hitecture」、Computational Science Education Proje
ct, Section 3.3.1 の、「Interleaved Memory」を参照
されたい。
るキュー処理やベクトル処理の効果が、図４のアルゴリ
ズム（１）の間接的なアドレッシングによって大幅に低
減されてしまうということにある。この方法が、常に一
定の間隔でメモリをアクセスすることができれば、より
好ましいものとなる。これに関しては、「Computer Arc
hitecture」、Computational Science Education Proje
ct, Section 3.3.1 の、「Interleaved Memory」を参照
されたい。
【００２７】更に、前記方法の処理は、超平面の大きさ
に応じて並列処理が行われず、多くの超平面Ｈ（１）
が、それだけでは複数のプロセッサに論理的に分散させ
るには小さすぎるために、分割された複数の並列メモリ
ブロックと、複数のベクトル・プロセッサを有するコン
ピュータには向いていない。また更に、超平面を扱う既
知の方法では、高次元空間又は前述した５点ステンシル
や７点ステンシル以外のステンシルへの一般化は知られ
ていない。これらの問題を解決するためのアプローチ
は、前述の式（１１）及び（１２）を正しく解くもので
はないが、その近似解を計算するものであり、前処理と
して同様の収束性向上をもたらす。通常、これにはvon-
Neumann級数展開が使用される。
に応じて並列処理が行われず、多くの超平面Ｈ（１）
が、それだけでは複数のプロセッサに論理的に分散させ
るには小さすぎるために、分割された複数の並列メモリ
ブロックと、複数のベクトル・プロセッサを有するコン
ピュータには向いていない。また更に、超平面を扱う既
知の方法では、高次元空間又は前述した５点ステンシル
や７点ステンシル以外のステンシルへの一般化は知られ
ていない。これらの問題を解決するためのアプローチ
は、前述の式（１１）及び（１２）を正しく解くもので
はないが、その近似解を計算するものであり、前処理と
して同様の収束性向上をもたらす。通常、これにはvon-
Neumann級数展開が使用される。
【００２８】 ｐ ＝（△−Ｌ）-1ｒ ＝ △-1（Ｉ−Ｌ△-1）-1ｒ ＝ △-1ｒ ＋ △-1 Ｌ△-1ｒ ＋（△-1Ｌ）2△-1ｒ ＋（△-1Ｌ）3△-1ｒ ＋ ．．．＋ （△-1 Ｌ）ntr△-1ｒ ．．．（１６）。
【００２９】この方法は、近似解を得るには、常に行列
とベクトルの積が計算されるだけであるため、実際には
十分にベクトル処理可能なものであり、並列に処理する
ことができる。しかし、超平面法と同様の収束性を実現
するためには、von-Neumann級数において、通常３桁又
は４桁のオーダーのメンバまで計算しなければならず、
計算コストは、３倍から４倍にまで膨らみ、ベクトル化
及び並列処理化の利点が相殺されてしまう。
とベクトルの積が計算されるだけであるため、実際には
十分にベクトル処理可能なものであり、並列に処理する
ことができる。しかし、超平面法と同様の収束性を実現
するためには、von-Neumann級数において、通常３桁又
は４桁のオーダーのメンバまで計算しなければならず、
計算コストは、３倍から４倍にまで膨らみ、ベクトル化
及び並列処理化の利点が相殺されてしまう。
【００３０】超平面を扱う方法を、ベクトル・コンピュ
ータに適するように変換する他のアプローチは、1995年
５月の、Takumi Washio、とKen Hayamiによる「Overlap
pedMulticolor MILU Preconditioning」、J. SCHI Comp
uter, Volume 16, No. 5,ページ636-651に開示されてい
る。
ータに適するように変換する他のアプローチは、1995年
５月の、Takumi Washio、とKen Hayamiによる「Overlap
pedMulticolor MILU Preconditioning」、J. SCHI Comp
uter, Volume 16, No. 5,ページ636-651に開示されてい
る。
【００３１】このアプローチでは、３次元空間及び７点
ステンシルについて、最初に超平面の分割が行われ、ｉ
ｓｔｒ個の部分集合が形成される。このことは、以下の
式（１７）によって表される。
ステンシルについて、最初に超平面の分割が行われ、ｉ
ｓｔｒ個の部分集合が形成される。このことは、以下の
式（１７）によって表される。
【００３２】 Ｃ（ｌ）＝{(i, j, r) ∈ Ｇ | mod(i + j + k - 3, istr) + 1 ＝ ｌ} ． ．．（１７）。
【００３３】言い換えれば、部分集合Ｃ（ｌ）は、超平
面Ｈ（ｌ）、Ｈ（ｌ＋ｉｓｔｒ）、Ｈ（ｌ＋２ｉｓｔ
ｒ）．．．．．の結合である。
面Ｈ（ｌ）、Ｈ（ｌ＋ｉｓｔｒ）、Ｈ（ｌ＋２ｉｓｔ
ｒ）．．．．．の結合である。
【００３４】カラー(colors)とも呼ばれる、この部分集
合のそれぞれは、前進代入においては、Ｃ（１）、Ｃ
（２）、．．．Ｃ（ｉｓｔｒ）、Ｃ（１）、Ｃ
（２）、．．．Ｃ（ｉｓｔｒ）の順で計算され、後退代
入においては、Ｃ（ｉｓｔｒ）、Ｃ（ｉｓｔｒ−
１）、．．．Ｃ（１）、Ｃ（ｉｓｔｒ）、Ｃ（ｉｓｔｒ
−１）、．．．Ｃ（１）の順で計算される。この処理順
序の例が、図６に示されている。未知の開始値は、それ
ぞれゼロにセットされる。この方法では、各カラーは、
少なくとも２回計算され、Ｈ（ｌ）に対するＨ（max(ｌ
-istr,1)）の依存度が考慮される。
合のそれぞれは、前進代入においては、Ｃ（１）、Ｃ
（２）、．．．Ｃ（ｉｓｔｒ）、Ｃ（１）、Ｃ
（２）、．．．Ｃ（ｉｓｔｒ）の順で計算され、後退代
入においては、Ｃ（ｉｓｔｒ）、Ｃ（ｉｓｔｒ−
１）、．．．Ｃ（１）、Ｃ（ｉｓｔｒ）、Ｃ（ｉｓｔｒ
−１）、．．．Ｃ（１）の順で計算される。この処理順
序の例が、図６に示されている。未知の開始値は、それ
ぞれゼロにセットされる。この方法では、各カラーは、
少なくとも２回計算され、Ｈ（ｌ）に対するＨ（max(ｌ
-istr,1)）の依存度が考慮される。
【００３５】図７には、この方法のアルゴリズム（２）
を実現する疑似コードの例が示されており、図８には、
対応するフローチャートが示されている。
を実現する疑似コードの例が示されており、図８には、
対応するフローチャートが示されている。
【００３６】ｉｓｔｒの選択を適切に行うことによっ
て、各カラーが、ベクトル処理に充分な大きさの数の格
子点を有するようにすることができる。前述の７点ステ
ンシルにおいて、例えば、nx＝nyでかつ、nzが任意であ
る場合、ｉｓｔｒ＝nx - 1 とすれば、Ｃ（ｌ）内の格
子点を計算するための全てのメモリ・アクセスが、一定
の間隔ｉｓｔｒで行われ、間接アドレッシングが回避さ
れる。一方、(nxnynz)/istrは、ベクトル長として十分
な大きさである。
て、各カラーが、ベクトル処理に充分な大きさの数の格
子点を有するようにすることができる。前述の７点ステ
ンシルにおいて、例えば、nx＝nyでかつ、nzが任意であ
る場合、ｉｓｔｒ＝nx - 1 とすれば、Ｃ（ｌ）内の格
子点を計算するための全てのメモリ・アクセスが、一定
の間隔ｉｓｔｒで行われ、間接アドレッシングが回避さ
れる。一方、(nxnynz)/istrは、ベクトル長として十分
な大きさである。
【００３７】例えば、ｘ軸方向がnx = 10で、ｙ軸方向
がny = 11で、ｚ軸方向がnz = 10で、ｉｓｔｒが前述の
とおり９と求められた、格子を有する７点ステンシルの
場合、カラーへの分割が行われた後、図６に示すような
計算が行われる。この図から、格子点ｉを計算するには
更に６つの格子点、i+1、i+10、i+110、i-1、i-10、i-1
10が必要とされ、これらは同じカラーに属さないことが
明らかである。従って、あるカラーの、ある格子点の計
算は、同じカラーの他の格子点の計算には依存しないこ
とが分かる。更に、間接アドレッシングは、ｉを知るこ
とによって直接把握できるので必要がない。そして、参
照すべき格子点、i+1、i+10、i+110、i-1、i-10、i-110
のデータが格納されているメモリ領域を事前にロードし
ておくことができるので、ベクトル・プロセッサの高速
処理が可能である。
がny = 11で、ｚ軸方向がnz = 10で、ｉｓｔｒが前述の
とおり９と求められた、格子を有する７点ステンシルの
場合、カラーへの分割が行われた後、図６に示すような
計算が行われる。この図から、格子点ｉを計算するには
更に６つの格子点、i+1、i+10、i+110、i-1、i-10、i-1
10が必要とされ、これらは同じカラーに属さないことが
明らかである。従って、あるカラーの、ある格子点の計
算は、同じカラーの他の格子点の計算には依存しないこ
とが分かる。更に、間接アドレッシングは、ｉを知るこ
とによって直接把握できるので必要がない。そして、参
照すべき格子点、i+1、i+10、i+110、i-1、i-10、i-110
のデータが格納されているメモリ領域を事前にロードし
ておくことができるので、ベクトル・プロセッサの高速
処理が可能である。
【００３８】この複数のカラーに対する処理技法を更に
改良したものが、前述の文献に記載されており、それは
「Overlapping Multicolor」技法と呼ばれている。ここ
では、前進代入に関して、カラーがＣ（ｉｓｔｒ−
ｗ）、Ｃ（ｉｓｔｒ−ｗ−１）、．．．Ｃ（ｉｓｔ
ｒ）、Ｃ（１）、Ｃ（２）、．．．Ｃ（ｉｓｔｒ）の順
で行われるという点で、複数のカラーを近似する技法の
改良が行われている。オーバーラッピング、即ち、追加
的に計算されるカラーの数はｗになる。対応するシステ
ムが後退代入に適用される。この、オーバーラッピング
多重カラー技法は、単一の多重カラー技法に較べて、優
れた収束性を示し、一方で、ベクトル長とメモリ・アク
セスに関しては同様の利点を有する。
改良したものが、前述の文献に記載されており、それは
「Overlapping Multicolor」技法と呼ばれている。ここ
では、前進代入に関して、カラーがＣ（ｉｓｔｒ−
ｗ）、Ｃ（ｉｓｔｒ−ｗ−１）、．．．Ｃ（ｉｓｔ
ｒ）、Ｃ（１）、Ｃ（２）、．．．Ｃ（ｉｓｔｒ）の順
で行われるという点で、複数のカラーを近似する技法の
改良が行われている。オーバーラッピング、即ち、追加
的に計算されるカラーの数はｗになる。対応するシステ
ムが後退代入に適用される。この、オーバーラッピング
多重カラー技法は、単一の多重カラー技法に較べて、優
れた収束性を示し、一方で、ベクトル長とメモリ・アク
セスに関しては同様の利点を有する。
【００３９】
【発明が解決しようとする課題】しかし、多重カラー技
法、及びオーバーラッピングを備えた多重カラー技法に
おいては、カラーの数を求めることと、前述の７点又は
５点ステンシル以外のケースに当該技法を適用すること
が困難である。最初の近似に関する開始値が、解から離
れすぎていると、カラーの数が小さくなりすぎ、結果的
に計算コストを増大させてしまう。ベクトル・コンピュ
ータのキューのベクトル長が短くなりすぎると、カラー
の数が大きくなりすぎ、結果的に計算スピードが低下し
てしまう。
法、及びオーバーラッピングを備えた多重カラー技法に
おいては、カラーの数を求めることと、前述の７点又は
５点ステンシル以外のケースに当該技法を適用すること
が困難である。最初の近似に関する開始値が、解から離
れすぎていると、カラーの数が小さくなりすぎ、結果的
に計算コストを増大させてしまう。ベクトル・コンピュ
ータのキューのベクトル長が短くなりすぎると、カラー
の数が大きくなりすぎ、結果的に計算スピードが低下し
てしまう。
【００４０】前述の７点又は５点ステンシル以外のステ
ンシルについてｉｓｔｒを求めることは、カラーの個々
の格子点が他の格子点に依存しないという条件を維持し
なければならないため、単純には実現できない。
ンシルについてｉｓｔｒを求めることは、カラーの個々
の格子点が他の格子点に依存しないという条件を維持し
なければならないため、単純には実現できない。
【００４１】
【課題を解決するための手段】本発明による方法は、特
に、規則的な正方格子に関する楕円型偏微分方程式を解
くのに適している。
に、規則的な正方格子に関する楕円型偏微分方程式を解
くのに適している。
【００４２】本発明の目的は、所与の長さのキューを備
えたベクトル・コンピュータにおいて、前述の式（１）
の解を求める方法を提供することにある。
えたベクトル・コンピュータにおいて、前述の式（１）
の解を求める方法を提供することにある。
【００４３】本発明の更なる目的は、特に、キューの長
さやメモリ・バンクの数といった、ベクトル・コンピュ
ータの実際の構成に適応して、他の方法より高速に演算
を行う方法を提供することにある。
さやメモリ・バンクの数といった、ベクトル・コンピュ
ータの実際の構成に適応して、他の方法より高速に演算
を行う方法を提供することにある。
【００４４】本発明の更なる目的は、計算コストを追加
せずに、５点ステンシル又は７点ステンシル以外のステ
ンシルにおける依存度を考慮する方法を提供することに
ある。
せずに、５点ステンシル又は７点ステンシル以外のステ
ンシルにおける依存度を考慮する方法を提供することに
ある。
【００４５】本発明の更なる目的は、好適に並列処理さ
れ、即ち、複数のプロセッサ上で並列に処理が行われ、
特に、分散メモリを有するコンピュータにおいても好適
に処理される方法を提供することにある。
れ、即ち、複数のプロセッサ上で並列に処理が行われ、
特に、分散メモリを有するコンピュータにおいても好適
に処理される方法を提供することにある。
【００４６】
【発明の実施の形態】本発明を簡単に理解するため、ま
ず、ベクトル・コンピュータの原理をいくつか説明す
る。
ず、ベクトル・コンピュータの原理をいくつか説明す
る。
【００４７】既存のコンピュータは、１つ、又はそれ以
上のプロセッサ、１つ、又はそれ以上のメモリ、入／出
力デバイス、及びこれらの間の通信チャネルを有してい
る。
上のプロセッサ、１つ、又はそれ以上のメモリ、入／出
力デバイス、及びこれらの間の通信チャネルを有してい
る。
【００４８】プロセッサは、コマンドを実行するため
に、そのコマンドをメモリにロードし、デコードし、実
行する。このような処理サイクルにかかる時間は、その
コマンドの複雑性及びプロセッサの内部構造に依存する
が、メモリ及び通信チャネルの構成にも依存する。
に、そのコマンドをメモリにロードし、デコードし、実
行する。このような処理サイクルにかかる時間は、その
コマンドの複雑性及びプロセッサの内部構造に依存する
が、メモリ及び通信チャネルの構成にも依存する。
【００４９】コンピュータ・プログラムは、１つ若しく
は２，３の単純なコマンドの繰り返しであるループ処理
を実行することが多い。そのようなループ処理の良い例
はベクトルの和の計算及び／又はベクトルのクロス積の
計算である。
は２，３の単純なコマンドの繰り返しであるループ処理
を実行することが多い。そのようなループ処理の良い例
はベクトルの和の計算及び／又はベクトルのクロス積の
計算である。
【００５０】ここでは、２進演算プロセッサにおいて２
つの数の積を単純に計算する場合であっても、その処理
には通常、一連の個別ステップが含まれることに注意す
べきである。このことは、例えば、ビットの列をそれぞ
れ左に１ビット移動させる（即ち２倍する）場合や、２
進数の足し算を行うような単純な場合のような、最も単
純な場合にも当てはまる。
つの数の積を単純に計算する場合であっても、その処理
には通常、一連の個別ステップが含まれることに注意す
べきである。このことは、例えば、ビットの列をそれぞ
れ左に１ビット移動させる（即ち２倍する）場合や、２
進数の足し算を行うような単純な場合のような、最も単
純な場合にも当てはまる。
【００５１】このような処理を行うアルゴリズムが、単
純な既知のスカラ・コンピュータによって処理される場
合、処理時間のほとんどは、同じコマンドのロードとデ
コード、及び中間の処理結果の記憶を繰り返すことに費
やされる。
純な既知のスカラ・コンピュータによって処理される場
合、処理時間のほとんどは、同じコマンドのロードとデ
コード、及び中間の処理結果の記憶を繰り返すことに費
やされる。
【００５２】この問題を解決するために、ベクトル・コ
ンピュータが使用される。
ンピュータが使用される。
【００５３】ベクトル・コンピュータの処理サイクルの
タイミングが、図９に概略的に示されている。
タイミングが、図９に概略的に示されている。
【００５４】図９に示す例の、最初のサイクルでは、第
１のコマンドがコマンド・メモリから取り出される。第
２のサイクルでは、このコマンドがプロセッサに渡さ
れ、オペランドの準備がされる。このサイクルにおいて
更に、第２のコマンドがコマンド・メモリから取り出さ
れる。第３のサイクルからｍ番目のサイクルまで、第１
のコマンドが実行される。また、第４のサイクルから、
（ｍ＋１）番目のサイクルまで、第２のコマンド他が実
行される。ｍ番目のサイクルが終了すると、プロセッサ
は、各サイクルの結果を出力する。この処理では、ｎ個
のコマンドの処理に全部でｍ＋ｎ−１回のサイクルを必
要とする。一方、この処理をスカラ・コンピュータで行
うと、（ｎｍ）サイクルを要することになる。
１のコマンドがコマンド・メモリから取り出される。第
２のサイクルでは、このコマンドがプロセッサに渡さ
れ、オペランドの準備がされる。このサイクルにおいて
更に、第２のコマンドがコマンド・メモリから取り出さ
れる。第３のサイクルからｍ番目のサイクルまで、第１
のコマンドが実行される。また、第４のサイクルから、
（ｍ＋１）番目のサイクルまで、第２のコマンド他が実
行される。ｍ番目のサイクルが終了すると、プロセッサ
は、各サイクルの結果を出力する。この処理では、ｎ個
のコマンドの処理に全部でｍ＋ｎ−１回のサイクルを必
要とする。一方、この処理をスカラ・コンピュータで行
うと、（ｎｍ）サイクルを要することになる。
【００５５】言い換えれば、ベクトル・コンピュータは
複数のキューを使用し、データやコマンドを、ローディ
ング・キューを介してベクトル・レジスタに転送する。
その後、算術演算キューにおいて計算が行われ、演算結
果はメモリ・キューを介してメモリに記憶される。算術
演算キューは、多機能キューであり、浮動小数点演算と
いった個別の複数ステップを並列で実行することができ
る。
複数のキューを使用し、データやコマンドを、ローディ
ング・キューを介してベクトル・レジスタに転送する。
その後、算術演算キューにおいて計算が行われ、演算結
果はメモリ・キューを介してメモリに記憶される。算術
演算キューは、多機能キューであり、浮動小数点演算と
いった個別の複数ステップを並列で実行することができ
る。
【００５６】１秒あたりの浮動小数点演算回数(FLOPS)
で表される、このベクトル・コンピュータの性能、即ち
計算能力は、いわゆるベクトル長に密接に依存してい
る。図９から分かるように、コマンド・スタックがベク
トル処理される場合、第１のコマンドの処理はｍサイク
ル必要であり、それ以降の各コマンドは、１サイクルだ
けの追加で処理される。この処理方法は、ｎ番目のコマ
ンドが、先行するｍ−１個のコマンドのうちの１つの処
理結果を知ることなく実行できることを前提としてい
る。これは、上述のベクトルの組み合わせにおいては典
型的なケースだからである。
で表される、このベクトル・コンピュータの性能、即ち
計算能力は、いわゆるベクトル長に密接に依存してい
る。図９から分かるように、コマンド・スタックがベク
トル処理される場合、第１のコマンドの処理はｍサイク
ル必要であり、それ以降の各コマンドは、１サイクルだ
けの追加で処理される。この処理方法は、ｎ番目のコマ
ンドが、先行するｍ−１個のコマンドのうちの１つの処
理結果を知ることなく実行できることを前提としてい
る。これは、上述のベクトルの組み合わせにおいては典
型的なケースだからである。
【００５７】従って、ベクトル長は、ベクトル化して実
行可能なコマンドの数を表す数字として認識される。従
来より、ベクトル・コンピュータに適したベクトル長
は、単純ループによって決定されうる。アルゴリズム
（３）は、そのループを実行するためのものであり、対
応する疑似コードが、図１０に示されている。
行可能なコマンドの数を表す数字として認識される。従
来より、ベクトル・コンピュータに適したベクトル長
は、単純ループによって決定されうる。アルゴリズム
（３）は、そのループを実行するためのものであり、対
応する疑似コードが、図１０に示されている。
【００５８】図１１は、ベクトル長がコンピュータの性
能に与える影響について、例えばＮＥＣのＳＸ４のよう
なコンピュータを例にして示した図である。最初は、ベ
クトル長の増加に伴って明らかにコンピュータの性能が
向上しているが、ベクトル長ｎが４００を越えると、コ
ンピュータの性能にほとんど変化がない。
能に与える影響について、例えばＮＥＣのＳＸ４のよう
なコンピュータを例にして示した図である。最初は、ベ
クトル長の増加に伴って明らかにコンピュータの性能が
向上しているが、ベクトル長ｎが４００を越えると、コ
ンピュータの性能にほとんど変化がない。
【００５９】この説明において、最小ベクトル長は、コ
ンピュータの性能にほとんど影響を及ぼさないようにな
るベクトル長を表し、この例では約４００である。
ンピュータの性能にほとんど影響を及ぼさないようにな
るベクトル長を表し、この例では約４００である。
【００６０】最小ベクトル長は、ベクトル・プロセッサ
の使用の形態に関わらず、そのベクトル・プロセッサに
依存する。従って、例えばFujitsu-VPXは、ベクトル処
理を効果的に実現するために、１０００のベクトル長を
要求する。一方、Cray EL98は、かなり短いベクトル長
でも効果的に機能する。
の使用の形態に関わらず、そのベクトル・プロセッサに
依存する。従って、例えばFujitsu-VPXは、ベクトル処
理を効果的に実現するために、１０００のベクトル長を
要求する。一方、Cray EL98は、かなり短いベクトル長
でも効果的に機能する。
【００６１】それぞれの最小ベクトル長より大きなベク
トル長が、コンピュータの性能にほとんど影響を与えな
いという事実は、ある理由によるものである。実際に、
図９に示すような、充分なオーバーラッピングが行われ
る状況では、例えば、キュー内のステップ数が、スタッ
ク内のコマンド数と等しくなる場合には、それ以上ベク
トル長を大きくしても僅かな効果しか得られない。
トル長が、コンピュータの性能にほとんど影響を与えな
いという事実は、ある理由によるものである。実際に、
図９に示すような、充分なオーバーラッピングが行われ
る状況では、例えば、キュー内のステップ数が、スタッ
ク内のコマンド数と等しくなる場合には、それ以上ベク
トル長を大きくしても僅かな効果しか得られない。
【００６２】コンピュータの性能を更に劣化させる要因
は、メモリの数、構成、及び構造である。ある記憶位置
にアクセスするため、それに関連付けられたメモリ・バ
ンクが常に事前ロードされ、アクセスが生じた場合だけ
参照される。次にメモリ・バンクは、再び戻される（ダ
ウンシフト）。こうした事前ロードとダウンシフトの時
間は、実際のアクセス時間の何倍にもなることがある。
従って、従来では、連続してアクセスされるビットが、
いくつかのメモリ・バンクに記憶されており、次のメモ
リ・バンクは常に事前ロードされていた。こうして、メ
モリ・アクセスに必要な時間がかなり短縮された。
は、メモリの数、構成、及び構造である。ある記憶位置
にアクセスするため、それに関連付けられたメモリ・バ
ンクが常に事前ロードされ、アクセスが生じた場合だけ
参照される。次にメモリ・バンクは、再び戻される（ダ
ウンシフト）。こうした事前ロードとダウンシフトの時
間は、実際のアクセス時間の何倍にもなることがある。
従って、従来では、連続してアクセスされるビットが、
いくつかのメモリ・バンクに記憶されており、次のメモ
リ・バンクは常に事前ロードされていた。こうして、メ
モリ・アクセスに必要な時間がかなり短縮された。
【００６３】しかし、これは、どのメモリ・バンクが連
続してアクセスされるかが常に分かっている場合を想定
している。間接的なアドレッシングを用いる、前述した
超平面を扱う方法においては、次にアクセスされるアド
レスを常に最初に計算しなければならないので、これを
想定することは不可能である。従って、間接的なアドレ
ッシングへの適用は避けなければならない。
続してアクセスされるかが常に分かっている場合を想定
している。間接的なアドレッシングを用いる、前述した
超平面を扱う方法においては、次にアクセスされるアド
レスを常に最初に計算しなければならないので、これを
想定することは不可能である。従って、間接的なアドレ
ッシングへの適用は避けなければならない。
【００６４】一方、既知で一定の間隔によりメモリがア
クセスされる方法は、ベクトル・コンピュータへの適用
に適している。既知で一定の前記間隔は、例えば、連続
してアクセスされる記憶位置の間隔が、常に既知で一定
であるような場合である。
クセスされる方法は、ベクトル・コンピュータへの適用
に適している。既知で一定の前記間隔は、例えば、連続
してアクセスされる記憶位置の間隔が、常に既知で一定
であるような場合である。
【００６５】直接的アドレッシングでないアクセスが同
じメモリ・バンクに対して短時間で連続して行われる場
合にも、同様のメモリ・アクセスの問題が生じる。これ
は、メモリ・バンクが常にダウンシフトされ、その後、
再開されたアクセスの前に再び事前ロードされるので、
コンピュータの性能をも低下させる。
じメモリ・バンクに対して短時間で連続して行われる場
合にも、同様のメモリ・アクセスの問題が生じる。これ
は、メモリ・バンクが常にダウンシフトされ、その後、
再開されたアクセスの前に再び事前ロードされるので、
コンピュータの性能をも低下させる。
【００６６】従って、ベクトル化処理されるループは、
同一メモリ・バンクへのアクセスの間隔が十分長いもの
であるよう設計されることが好ましい。
同一メモリ・バンクへのアクセスの間隔が十分長いもの
であるよう設計されることが好ましい。
【００６７】コマンドの処理が複数の独立した並列処理
ベクトル・プロセッサに亘って分散されると、コンピュ
ータの性能は一層向上する。このための条件は、複数の
独立した部分タスクを形成できることである。
ベクトル・プロセッサに亘って分散されると、コンピュ
ータの性能は一層向上する。このための条件は、複数の
独立した部分タスクを形成できることである。
【００６８】前述したように、本発明は、ベクトル・コ
ンピュータ上で多重カラーを扱う方法を実行し、式
（１）の解を求める方法に関するものである。
ンピュータ上で多重カラーを扱う方法を実行し、式
（１）の解を求める方法に関するものである。
【００６９】コンピュータの性能に悪影響を及ぼす、前
述した要因を回避するために、本発明は、一定のメモリ
・アクセス間隔ｉｓｔｒを以下の規則に従って決定す
る。
述した要因を回避するために、本発明は、一定のメモリ
・アクセス間隔ｉｓｔｒを以下の規則に従って決定す
る。
【００７０】ａ）すべてのｋ(k = 1 ...nzr)に対して、
mod(izrk, istr) ≠ 0であること。ここで、izrk = {j:
Li, i-j ≠ 0 又は Ui, i+j ≠ 0}、nzrは、Li, i-j ≠
0 又は Ui, i+j ≠ 0である全てのｊの数である。
mod(izrk, istr) ≠ 0であること。ここで、izrk = {j:
Li, i-j ≠ 0 又は Ui, i+j ≠ 0}、nzrは、Li, i-j ≠
0 又は Ui, i+j ≠ 0である全てのｊの数である。
【００７１】ｂ）istrとnbの最大公約数(gcd)が１であ
ること。ここで、nbはメモリ・バンクの数である。
ること。ここで、nbはメモリ・バンクの数である。
【００７２】ｃ）istrはできるだけ小さな数であるこ
と。その結果、実際のベクトル長が、そのベクトル・プ
ロセッサに要求される最小ベクトル長より大きいものと
なる。
と。その結果、実際のベクトル長が、そのベクトル・プ
ロセッサに要求される最小ベクトル長より大きいものと
なる。
【００７３】また、istrの実際の値によって、以下のよ
うな値ｘを求め、その値がなるべく大きくなるようにす
る。
うな値ｘを求め、その値がなるべく大きくなるようにす
る。
【００７４】 x = min{int(istr/mod (izrk, istr)) k = 1,...nzr} ．．．（１８）。
【００７５】規則ａ）は、全ての計算処理が、あるカラ
ーの中では独立して実行されることを保証している。こ
のことによって、ベクトル処理が、１つのカラーにおい
て可能になる。
ーの中では独立して実行されることを保証している。こ
のことによって、ベクトル処理が、１つのカラーにおい
て可能になる。
【００７６】これは、ｉｓｔｒによって割り切れるizrk
がないことを仮定している。この規則を用いることによ
って、５点ステンシルや７点ステンシルを扱うことがで
きるだけでなく、あらゆる依存度を考慮に入れることが
できる。
がないことを仮定している。この規則を用いることによ
って、５点ステンシルや７点ステンシルを扱うことがで
きるだけでなく、あらゆる依存度を考慮に入れることが
できる。
【００７７】規則ｂ）は、メモリの競合が無いことを保
証する。即ち、同じメモリ・バンクに対するアクセスを
連続してすることができず、かつ、全てのメモリ・バン
クが使用されている。
証する。即ち、同じメモリ・バンクに対するアクセスを
連続してすることができず、かつ、全てのメモリ・バン
クが使用されている。
【００７８】規則ｃ）は、ベクトル長が、ベクトル・コ
ンピュータを最適動作させるために必要な程度に大きい
ことを保証する。これに関しては、ベクトル長を大きく
する程、この方法による収束性が低下することに注意す
べきである。これは、前進代入及び後退代入の開始時点
で、すべてのｉに対して要素ｐi又はｑiの値がそれぞれ
ゼロにセットされる（又は別の適切な開始値がセットさ
れる）ことによるものである。ループの最初の１回が終
了すると、要素ｐi+istr, ｐi+2istr, ...ｑmは、それ
ぞれ第１の近似値に対応した新しい値がセットされる。
次に、他の要素ｐi又はｑiの全てについて計算を更に行
っていく。更なる計算に利用できる第１の近似値が迅速
に求められるほど、収束性がより効果的なものになる。
ンピュータを最適動作させるために必要な程度に大きい
ことを保証する。これに関しては、ベクトル長を大きく
する程、この方法による収束性が低下することに注意す
べきである。これは、前進代入及び後退代入の開始時点
で、すべてのｉに対して要素ｐi又はｑiの値がそれぞれ
ゼロにセットされる（又は別の適切な開始値がセットさ
れる）ことによるものである。ループの最初の１回が終
了すると、要素ｐi+istr, ｐi+2istr, ...ｑmは、それ
ぞれ第１の近似値に対応した新しい値がセットされる。
次に、他の要素ｐi又はｑiの全てについて計算を更に行
っていく。更なる計算に利用できる第１の近似値が迅速
に求められるほど、収束性がより効果的なものになる。
【００７９】そこで、良好な収束性を得るために、式
（１８）が提供される。これは、前述のvon-Neumann級
数展開に関連するものであり、本発明による方法の収束
性が、式（１６）において、ntr =（nd - 1）xである限
り、式（１６）による級数展開と同程度であることを保
証する。
（１８）が提供される。これは、前述のvon-Neumann級
数展開に関連するものであり、本発明による方法の収束
性が、式（１６）において、ntr =（nd - 1）xである限
り、式（１６）による級数展開と同程度であることを保
証する。
【００８０】ｉｓｔｒを計算するために好適なアルゴリ
ズム（４）を実現する疑似コードの例を、図１２に示
す。
ズム（４）を実現する疑似コードの例を、図１２に示
す。
【００８１】また、以下の説明に関連して、前記アルゴ
リズム（４）のフローチャートを図１３に示す。
リズム（４）のフローチャートを図１３に示す。
【００８２】ここでは、lvcが最小ベクトル長、即ち、
ベクトル・プロセッサによって要求される最小のベクト
ル長であり、gcd(i, nb)が、ｉとｎｂの最大公約数であ
る。
ベクトル・プロセッサによって要求される最小のベクト
ル長であり、gcd(i, nb)が、ｉとｎｂの最大公約数であ
る。
【００８３】また、最小ベクトル長は、厳密な境界を有
していないので、特に最後の規則ｃ）と式（１８）の適
用については、ファジー理論と関連付けて考えることが
できる。ファジー理論を取り入れると、前述のアルゴリ
ズム（４）では、例えば、ｉによってループするループ
は、ｉがistrmaxになるまでではなく、istrmaxの何倍か
になるまでループするようにでき、あるｉ(i > istrma
x)における解は、ｉがistrmaxより大きいという点で、
他の解より好ましくないものとすることができる。以下
では、前述した本発明の適用について、np=2のベクトル
・プロセッサを備えた並列ベクトル・コンピュータに関
連して説明する。当業者には、プロセッサが３つ以上の
ベクトル・コンピュータにも、本発明を同様に適用でき
ることは明らかである。
していないので、特に最後の規則ｃ）と式（１８）の適
用については、ファジー理論と関連付けて考えることが
できる。ファジー理論を取り入れると、前述のアルゴリ
ズム（４）では、例えば、ｉによってループするループ
は、ｉがistrmaxになるまでではなく、istrmaxの何倍か
になるまでループするようにでき、あるｉ(i > istrma
x)における解は、ｉがistrmaxより大きいという点で、
他の解より好ましくないものとすることができる。以下
では、前述した本発明の適用について、np=2のベクトル
・プロセッサを備えた並列ベクトル・コンピュータに関
連して説明する。当業者には、プロセッサが３つ以上の
ベクトル・コンピュータにも、本発明を同様に適用でき
ることは明らかである。
【００８４】図１４には、２つのプロセッサによって、
図２に示す７点ステンシルを処理する順序が示されてい
る。ここで、nx = 10, ny = 11, istr = 9, 及びnzは任
意である。
図２に示す７点ステンシルを処理する順序が示されてい
る。ここで、nx = 10, ny = 11, istr = 9, 及びnzは任
意である。
【００８５】最初に、第２のプロセッサと独立して動作
する第１のプロセッサは、カラーＣ（１）及びＣ（２）
について第１の近似解を計算する。同時に、第２のプロ
セッサは、第１のプロセッサとは独立して、カラーＣ
（５）及びＣ（６）についての第１の近似解を計算す
る。このために、カラーＣ（８）、Ｃ（９）又はＣ
（３）、Ｃ（４）の格子点の第１の近似解の値は、それ
ぞれゼロにセットされる。
する第１のプロセッサは、カラーＣ（１）及びＣ（２）
について第１の近似解を計算する。同時に、第２のプロ
セッサは、第１のプロセッサとは独立して、カラーＣ
（５）及びＣ（６）についての第１の近似解を計算す
る。このために、カラーＣ（８）、Ｃ（９）又はＣ
（３）、Ｃ（４）の格子点の第１の近似解の値は、それ
ぞれゼロにセットされる。
【００８６】カラーＣ（２）又はＣ（６）の計算の後、
同期がとられる。即ち、次のカラーＣ（３）又はＣ
（７）の計算は、両方のプロセッサがそれぞれの先行す
るカラーの計算を終えた後でなければ、始めることがで
きない。カラーＣ（３）又はＣ（７）の値は、カラーＣ
（５）又はＣ（１）の値の計算に用いられるので、ここ
でも同期が必要である。カラーＣ（５）の計算がまだ終
了していない間に、カラーＣ（３）の値を変更すること
は、計算の不確実性を招く。それは、カラーＣ（３）の
値が、カラーＣ（５）の計算に必要だからである。結果
的に、カラーＣ（５）の結果は、それぞれのプロセッサ
のランダムな動作条件に、多かれ少なかれ依存すること
になる。
同期がとられる。即ち、次のカラーＣ（３）又はＣ
（７）の計算は、両方のプロセッサがそれぞれの先行す
るカラーの計算を終えた後でなければ、始めることがで
きない。カラーＣ（３）又はＣ（７）の値は、カラーＣ
（５）又はＣ（１）の値の計算に用いられるので、ここ
でも同期が必要である。カラーＣ（５）の計算がまだ終
了していない間に、カラーＣ（３）の値を変更すること
は、計算の不確実性を招く。それは、カラーＣ（３）の
値が、カラーＣ（５）の計算に必要だからである。結果
的に、カラーＣ（５）の結果は、それぞれのプロセッサ
のランダムな動作条件に、多かれ少なかれ依存すること
になる。
【００８７】同期をとった後で、カラーＣ（３）及びＣ
（４）は、第１のプロセッサによって計算され、カラー
Ｃ（７）、Ｃ（８）、及びＣ（９）は他のプロセッサに
よって計算される。次に、これらのプロセッサは再び同
期をとり、処理は新たなパスに進み、カラーＣ（１）又
はＣ（５）の計算が開始される。また、カラーＣ
（８）、Ｃ（９）又はＣ（３）、Ｃ（４）の第１の近似
解は、以前のパスで計算されている。このループは、条
件id = ndになるまで、アルゴリズム（２）によって繰
り返される。即ち、近似解が十分良好なものになること
が保証されるまで繰り返される。通常、満足な近似解を
得るには、２つのパス(nd - 2)で充分である。
（４）は、第１のプロセッサによって計算され、カラー
Ｃ（７）、Ｃ（８）、及びＣ（９）は他のプロセッサに
よって計算される。次に、これらのプロセッサは再び同
期をとり、処理は新たなパスに進み、カラーＣ（１）又
はＣ（５）の計算が開始される。また、カラーＣ
（８）、Ｃ（９）又はＣ（３）、Ｃ（４）の第１の近似
解は、以前のパスで計算されている。このループは、条
件id = ndになるまで、アルゴリズム（２）によって繰
り返される。即ち、近似解が十分良好なものになること
が保証されるまで繰り返される。通常、満足な近似解を
得るには、２つのパス(nd - 2)で充分である。
【００８８】並列ベクトル・コンピュータで上記処理を
行うための好適アルゴリズム（５）を実行するための疑
似コードの例が図１５、図１６に示されている。図１５
に示す疑似コードと図１６に示す疑似コードとは連続し
ているものであり、これら２つで１つのアルゴリズム
（５）を示していることに注意すべきである。
行うための好適アルゴリズム（５）を実行するための疑
似コードの例が図１５、図１６に示されている。図１５
に示す疑似コードと図１６に示す疑似コードとは連続し
ているものであり、これら２つで１つのアルゴリズム
（５）を示していることに注意すべきである。
【００８９】また、以降の説明に関連して、このアルゴ
リズム（５）に対応するフローチャートが、図１７に示
されている。
リズム（５）に対応するフローチャートが、図１７に示
されている。
【００９０】２つ以上のプロセッサがあれば、ループ処
理は、ｉａを使って、０からint(istr/np) - 1 まで実
行される。またここで、npはプロセッサの数を表してい
る。
理は、ｉａを使って、０からint(istr/np) - 1 まで実
行される。またここで、npはプロセッサの数を表してい
る。
【００９１】
【発明の効果】以上に説明したように、本発明による方
法によって、様々なキューの長さ、及びメモリ・バンク
の数を有するようなベクトル・コンピュータに適用させ
ることができ、結果的に既存の方法より高速な動作を実
現できる。また、いかなる追加コストも必要とせずに、
５点ステンシルや７点ステンシルにおける依存度以外の
依存度を更に考慮することができる。
法によって、様々なキューの長さ、及びメモリ・バンク
の数を有するようなベクトル・コンピュータに適用させ
ることができ、結果的に既存の方法より高速な動作を実
現できる。また、いかなる追加コストも必要とせずに、
５点ステンシルや７点ステンシルにおける依存度以外の
依存度を更に考慮することができる。
【００９２】例えば、図１８及び図１９にそれぞれ示
す、１９及び２７点ステンシルについても計算が可能で
ある。
す、１９及び２７点ステンシルについても計算が可能で
ある。
【００９３】更に、本発明による方法は、複数の独立プ
ロセッサを用いた並列処理に非常に適している。特に、
分散メモリを備え、各プロセッサがそのメモリにアロケ
ートされているメモリ領域を有し、アクセスの競合が回
避されるようなコンピュータ・システムに適している。
ロセッサを用いた並列処理に非常に適している。特に、
分散メモリを備え、各プロセッサがそのメモリにアロケ
ートされているメモリ領域を有し、アクセスの競合が回
避されるようなコンピュータ・システムに適している。
【図１】５点ステンシルの２次元格子を示す図である。
【図２】７点ステンシルの３次元格子を示す図である。
【図３】７点ステンシルに関する行列Ｍの概略を示す図
である。
である。
【図４】１次方程式の解を前進代入及び後退代入を用い
て解くためのアルゴリズム（１）を実現する疑似コード
の例を示す図である。
て解くためのアルゴリズム（１）を実現する疑似コード
の例を示す図である。
【図５】３次元空間における超平面の位置の概略を示す
図である。
図である。
【図６】nx = 10, ny = 11, istr = 9で、nzが任意であ
る場合の計算順序を概略示す図である。
る場合の計算順序を概略示す図である。
【図７】アルゴリズム（２）を実現する疑似コードの例
を示す図である。
を示す図である。
【図８】アルゴリズム（２）についてのフローチャート
である。
である。
【図９】ベクトル・コンピュータの処理タイミングを示
す図である。
す図である。
【図１０】ベクトル長を決定するアルゴリズム（３）を
実現する疑似コードの例を示す図である。
実現する疑似コードの例を示す図である。
【図１１】既存のベクトル・コンピュータのベクトル長
と処理能力の関係を表すグラフである。
と処理能力の関係を表すグラフである。
【図１２】ｉｓｔｒを計算するために好適なアルゴリズ
ム（４）を実現する疑似コードの例を示す図である。
ム（４）を実現する疑似コードの例を示す図である。
【図１３】図１２のアルゴリズム（４）についてのフロ
ーチャートである。
ーチャートである。
【図１４】nx = 10, ny = 11, istr = 9で、nzが任意で
あり、２つのプロセッサ上に分散された処理が並列に行
われる場合の計算順序を概略示す図である。
あり、２つのプロセッサ上に分散された処理が並列に行
われる場合の計算順序を概略示す図である。
【図１５】アルゴリズム（５）を実現するための疑似コ
ードの例の一部を示す図である。
ードの例の一部を示す図である。
【図１６】アルゴリズム（５）を実現するための疑似コ
ードの例の一部を示す図である。
ードの例の一部を示す図である。
【図１７】図１５及び図１６のアルゴリズム（５）につ
いてのフローチャートである。
いてのフローチャートである。
【図１８】３次元空間における、１９点のステンシルの
配置を示す図である。
配置を示す図である。
【図１９】３次元空間における、２７点のステンシルの
配置を示す図である。
配置を示す図である。
Claims (6)
- 【請求項１】 少なくとも１つのベクトル・プロセッサ
とメモリ・バンクを有するベクトル・コンピュータで、
下三角行列と上三角行列を含んだ行列式の近似解を求め
るための、前進代入及び後退代入を行う多重カラー技法
に基づく計算方法において、前記計算方法が、前記代入
を行うために、メモリをアクセスする一定の間隔ｉｓｔ
ｒを有し、 前記代入の前に、所定の基準により、前記ｉｓｔｒを決
定するステップを有し、前記ｉｓｔｒが以下のステップ
により求められることを特徴とする計算方法、 ａ）前記下三角行列及び上三角行列において対角成分の
並びと平行な、非対角成分のそれぞれの並びの中で、少
なくとも１つの成分がゼロでない各並びに対して、その
対角成分の並びから見て、その非対角成分の並びがいく
つ目の平行な要素かを表す数ｉｚｒをｉｓｔｒで割った
場合に、その剰余がゼロにならないｉｓｔｒを選択する
ステップ、 ｂ）ｉｓｔｒとメモリ・バンクの数ｎｂの最大公約数
（ＧＣＤ）を計算し、その結果が１となるｉｓｔｒを選
択するステップ、及び ｃ）可能な数の候補の中で最も小さいｉｓｔｒを求める
ステップ。 - 【請求項２】 請求項１において、前記行列式が、 （△−Ｌ）△-1（△−Ｕ）ｑ＝ｒ であり、ここで、△はｍ×ｍの対角行列、Ｌはｍ×ｍの
下三角行列、Ｕはｍ×ｍの上三角行列、ｑ及びｒはｍ次
元のベクトルであることを特徴とする計算方法。 - 【請求項３】 請求項１において、前記ｉｓｔｒを求め
るステップが更に、 ｄ）各ｉｚｒに関して、ｉｚｒをｉｓｔｒで割った前記
剰余で更にｉｓｔｒを割り、その整数部の値が最小であ
るものをｘとする場合、前記ｘの値が、前記ベクトル・
コンピュータで並列に動作するベクトル・プロセッサの
数の少なくとも２倍以上となるように、ｉｓｔｒの値を
決定するステップを有する計算方法。 - 【請求項４】 請求項１において、前記ステップｃ）
が、ファジー理論によって判断されることを特徴とする
計算方法。 - 【請求項５】 請求項３において、前記ステップｄ）
が、ファジー理論によって判断されることを特徴とする
計算方法。 - 【請求項６】 少なくとも１つのベクトル・プロセッサ
とメモリ・バンクを有するベクトル・コンピュータで、
下三角行列と上三角行列を含んだ行列式の近似解を求め
るための、前進代入及び後退代入を行う多重カラー技法
に基づく計算方法を実現させるプログラムを記録したコ
ンピュータ読み取り可能な記録媒体であって、前記プロ
グラムは、 前記代入の前に、所定の基準により、前記ｉｓｔｒを決
定するステップを有し、前記ｉｓｔｒが以下のステップ
により求められることを特徴とするプログラムを記録し
たコンピュータ読み取り可能な記録媒体、 ａ）前記下三角行列及び上三角行列において対角成分の
並びと平行な、非対角成分のそれぞれの並びの中で、少
なくとも１つの成分がゼロでない各並びに対して、その
対角成分の並びから見て、その非対角成分の並びがいく
つ目の平行な要素かを表す数ｉｚｒをｉｓｔｒで割った
場合に、その剰余がゼロにならないｉｓｔｒを選択する
ステップ、 ｂ）ｉｓｔｒとメモリ・バンクの数ｎｂの最大公約数
（ＧＣＤ）を計算し、その結果が１となるｉｓｔｒを選
択するステップ、及び ｃ）可能な数の候補の中で最も小さいｉｓｔｒを求める
ステップ。
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| DE19827238.3 | 1998-06-18 | ||
| DE19827238A DE19827238B4 (de) | 1998-06-18 | 1998-06-18 | Verfahren zum Betrieb eines Vektorrechners |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2000029864A true JP2000029864A (ja) | 2000-01-28 |
| JP3659307B2 JP3659307B2 (ja) | 2005-06-15 |
Family
ID=7871327
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP17131599A Expired - Fee Related JP3659307B2 (ja) | 1998-06-18 | 1999-06-17 | ベクトル・コンピュータにおける演算方法、及び記録媒体 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US6446105B1 (ja) |
| JP (1) | JP3659307B2 (ja) |
| DE (1) | DE19827238B4 (ja) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6694343B2 (en) * | 2001-02-08 | 2004-02-17 | International Business Machines Corporation | Method for solving a large sparse triangular system of linear equations |
| US7313788B2 (en) * | 2003-10-29 | 2007-12-25 | International Business Machines Corporation | Vectorization in a SIMdD DSP architecture |
| US8306798B2 (en) * | 2008-12-24 | 2012-11-06 | Intel Corporation | Fluid dynamics simulator |
| US20100199067A1 (en) * | 2009-02-02 | 2010-08-05 | International Business Machines Corporation | Split Vector Loads and Stores with Stride Separated Words |
| CN103514109B (zh) * | 2013-09-24 | 2016-04-13 | 创新科存储技术有限公司 | 一种打开磁盘写缓存的方法和装置 |
| US9632801B2 (en) * | 2014-04-09 | 2017-04-25 | Intel Corporation | Banked memory access efficiency by a graphics processor |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5267185A (en) * | 1989-04-14 | 1993-11-30 | Sharp Kabushiki Kaisha | Apparatus for calculating matrices |
| JPH0444165A (ja) * | 1990-06-12 | 1992-02-13 | Nec Corp | 対称連立一次方程式の求解方式 |
| JP2718254B2 (ja) * | 1990-10-02 | 1998-02-25 | 日本電気株式会社 | ベクトル処理装置 |
| US5446908A (en) * | 1992-10-21 | 1995-08-29 | The United States Of America As Represented By The Secretary Of The Navy | Method and apparatus for pre-processing inputs to parallel architecture computers |
| US5680338A (en) * | 1995-01-04 | 1997-10-21 | International Business Machines Corporation | Method and system for vector processing utilizing selected vector elements |
-
1998
- 1998-06-18 DE DE19827238A patent/DE19827238B4/de not_active Expired - Fee Related
-
1999
- 1999-06-17 JP JP17131599A patent/JP3659307B2/ja not_active Expired - Fee Related
- 1999-06-18 US US09/335,705 patent/US6446105B1/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| DE19827238A1 (de) | 1999-12-23 |
| DE19827238B4 (de) | 2004-09-30 |
| US6446105B1 (en) | 2002-09-03 |
| JP3659307B2 (ja) | 2005-06-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10410112B2 (en) | Apparatus and method for performing a forward operation of artificil neural networks | |
| US7337205B2 (en) | Matrix multiplication in a vector processing system | |
| KR101202445B1 (ko) | 프로세서 | |
| US5600843A (en) | Ring systolic array system for synchronously performing matrix/neuron computation using data transferred through cyclic shift register connected in cascade of trays | |
| US11983616B2 (en) | Methods and apparatus for constructing digital circuits for performing matrix operations | |
| US20190138922A1 (en) | Apparatus and methods for forward propagation in neural networks supporting discrete data | |
| Bikov et al. | Parallel fast Walsh transform algorithm and its implementation with CUDA on GPUs | |
| US10884736B1 (en) | Method and apparatus for a low energy programmable vector processing unit for neural networks backend processing | |
| US20200226201A1 (en) | Methods and Apparatus for Constructing Digital Circuits for Performing Matrix Operations | |
| US7467053B2 (en) | Three-dimensional fourier transform processing method for shared memory scalar parallel computer | |
| JP2000029864A (ja) | ベクトル・コンピュ―タにおける演算方法、及び記録媒体 | |
| JP2023070746A (ja) | 情報処理プログラム、情報処理装置、及び情報処理方法 | |
| CN110399976B (zh) | 计算装置和计算方法 | |
| US20190130274A1 (en) | Apparatus and methods for backward propagation in neural networks supporting discrete data | |
| JP2001067206A (ja) | モジュラー乗算を実行するためのシステム並びに方法 | |
| CN110673824A (zh) | 矩阵向量乘电路以及循环神经网络硬件加速器 | |
| Ichimura et al. | Threaded accurate matrix-matrix multiplications with sparse matrix-vector multiplications | |
| US9600446B2 (en) | Parallel multicolor incomplete LU factorization preconditioning processor and method of use thereof | |
| US20190073584A1 (en) | Apparatus and methods for forward propagation in neural networks supporting discrete data | |
| JP3808925B2 (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 | |
| Myllykoski et al. | On solving separable block tridiagonal linear systems using a GPU implementation of radix-4 PSCR method | |
| JP2806262B2 (ja) | マルチプロセッサシステムのプロセス割当方法 | |
| GB2523805A (en) | Data processing apparatus and method for performing vector scan operation | |
| JP7506086B2 (ja) | データ処理 | |
| JP2006164307A (ja) | 多様な行列格納法を使用可能な連立方程式の並列処理装置および方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20040331 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040526 |
|
| TRDD | Decision of grant or rejection written | ||
| RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20050216 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20050223 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20050308 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080325 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090325 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100325 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100325 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110325 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110325 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120325 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120325 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130325 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130325 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140325 Year of fee payment: 9 |
|
| LAPS | Cancellation because of no payment of annual fees |