JP2004348493A - 通信隠蔽型の並列高速フーリエ変換方法 - Google Patents

Abstract

【課題】並列計算機による３次元の高速フーリエ変換において、プロセッサ間でのデータ転送によるオーバーヘッドを削減し、並列化効率を向上させることを目的とする。
【解決手段】データをＸ座標が偶数の要素と奇数の要素とに分割し、前者に対してＹ方向の変換を行うと同時に後者に対してはプロセッサ間での転置処理を行い（処理３４）、後者に対してＸ方向の変換を行うと同時に前者に対してはプロセッサ間での転置処理を行う（処理３５）。
【効果】プロセッサ間での転置処理のための通信時間を演算時間に隠蔽することが可能となり、従来法に比べて並列化効率を向上させることができる。
【選択図】 図１

Description

【０００１】
【発明の属する技術分野】
本発明は分散メモリ型の並列計算機を用いて高速フーリエ変換を行うための方式に関する。
【０００２】
【従来の技術】
半導体デバイスの特性計算、電子状態計算、気象予測のための計算などの科学技術計算では、数万から数百万に上る変数を扱う大規模シミュレーションが必要である。このような大規模問題を扱う手段としては、並列計算機、特に分散メモリ型と呼ばれるタイプの並列計算機が有力である。分散メモリ型並列計算機は、各々がメモリを持つプロセッサをネットワークで複数個結合したシステムであり、従来の逐次型計算機に比べ、プロセッサ台数を増やすことでピーク性能をいくらでも高めることができるという利点を持つ。
分散メモリ型並列計算機では、計算対象のデータを各プロセッサの持つメモリに分散させて格納し、そのそれぞれに対して各プロセッサが並列に計算を行う。もし計算途中で他のプロセッサの持つデータが必要となる場合は、必要なデータをそのプロセッサから転送してもらい、その後に計算を続行する。したがって分散メモリ型並列計算機では、一般に、演算時間に加えてプロセッサ間でデータを転送するための時間が必要である。そこで、計算効率を上げるには、プロセッサ間での通信がなるべく少なくて済むような並列性の高い計算手法が必要である。また、多くの分散メモリ型並列計算機においては、あるデータに対して演算を行っている間に、別のデータを他プロセッサに転送する機構が用意されている。この場合は、演算とデータ転送が同時に行えるように計算手法を工夫すれば、転送のための時間が演算時間に隠蔽でき、計算効率を向上させることができる。
【０００３】
科学技術計算において頻繁に利用される処理の一つに、フーリエ変換がある。フーリエ変換は、ある実数区間で定義された複素数値をとる関数ｆ（ｘ）を複素指数関数ｅｘｐ（ｉｋｘ）の重ね合わせとして表す処理であり、計算機上で実現する場合には、扱いうる点の数が有限であることから、複素数の点列ｆ_０、ｆ_１，… ，ｆ_Ｎ−１をＮ個の複素指数関数
ｅｘｐ（２πｉｋｊ／Ｎ）（ただしｋ＝０，１，… ，Ｎ−１で、ｉは虚数単位、πは円周率）
の重ね合わせとして
ｆ_ｊ＝Σ_ｋ＝０ ^Ｎ−１ ｃ_ｋ ｅｘｐ（２πｉｋｊ／Ｎ） （ただしｊ＝０，１，… ，Ｎ−１）
と表す処理となる。すなわち、 ｆ_０，ｆ_１，… ，ｆ_Ｎ−１が与えられたときに、重ね合わせの係数ｃ_０，ｃ_１，… ，ｃ_Ｎ−１を求める処理がフーリエ変換である。この係数ｃ_０，ｃ_１，… ，ｃ_Ｎ−１は、式
ｃ_ｋ＝（１／Ｎ）Σ_ｊ＝０ ^Ｎ−１ ｆ_ｊ ｅｘｐ（−２πｉｋｊ／Ｎ） （ただしｋ＝０，１，… ，Ｎ−１）
により計算できることが知られているが、この定義に基づいて計算を行うと、式の数がＮ本あり、各式がＮ個の項から成るため、複素指数関数 ｅｘｐ（−２πｉｋｊ／Ｎ）の計算に加えて、複素数の加算と乗算が約Ｎ^２回必要である。そこで実際には、アルゴリズム上の工夫により計算量を約ＮｌｏｇＮのオーダーに減少させた高速フーリエ変換という手法が広く使われている。高速フーリエ変換の計算法については、たとえばＧ． Ｇｏｌｕｂ ａｎｄ Ｃ． Ｆ． ｖａｎ Ｌｏａｎ： ”Ｍａｔｒｉｘ Ｃｏｍｐｕｔａｔｉｏｎｓ”， ３^ｒｄ ｅｄｉｔｉｏｎ， Ｔｈｅ Ｊｏｈｎｓ Ｈｏｐｋｉｎｓ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ， １９９６， ｐｐ．１８９〜１９２（非特許文献１）に詳しく記述されている。
以上で述べたフーリエ変換は１次元のフーリエ変換と呼ばれるが、半導体デバイスの特性計算、電子状態計算、気象予測のための計算などでは、３次元のフーリエ変換が用いられる。これは、３つの添字ｊ_ｘ，ｊ_ｙ，ｊ_ｚを持つ複素数データ
｛ｆ_{ｊｘ，ｊｙ，ｊｚ}｝（ただしｊ_ｘ＝０，１，… ，Ｎ_ｘ−１，ｊ_ｙ＝０，１，… ，Ｎ_ｙ−１，ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）をＮ_ｘ×Ｎ_ｙ×Ｎ_ｚ個の複素指数関数
ｅｘｐ（２πｉｋ_ｘｊ_ｘ／Ｎ_ｘ） ｅｘｐ（２πｉｋ_ｙｊ_ｙ／Ｎ_ｙ） ｅｘｐ（２πｉｋ_ｚｊ_ｚ／Ｎ_ｚ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｋ_ｚ＝０，１，… ，Ｎ_ｚ−１）
の重ね合わせとして
ｆ_{ｊｘ，ｊｙ，ｊｚ}＝Σ_ｋｘ＝０ ^Ｎｘ−１Σ_ｋｙ＝０ ^Ｎｙ−１Σ_ｋｚ＝０ ^Ｎｚ−１
ｃ_{ｋｘ，ｋｙ，ｋｚ} ｅｘｐ（２πｉｋ_ｘｊ_ｘ／Ｎ_ｘ） ｅｘｐ（２πｉｋ_ｙｊ_ｙ／Ｎ_ｙ） ｅｘｐ（２πｉｋ_ｚｊ_ｚ／Ｎ_ｚ）
（ただしｊ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｊ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
と表す処理となる。すなわち、｛ｆ_{ｊｘ，ｊｙ，ｊｚ}｝が与えられたときに、重ね合わせの係数｛ｃ_{ｋｘ，ｋｙ，ｋｚ}｝を求める処理が３次元フーリエ変換である。この係数｛ｃ_{ｋｘ，ｋｙ，ｋｚ}｝は、次の式によって計算できることが知られている。
ｃ_{ｋｘ，ｋｙ，ｋｚ}＝Σ_ｊｘ＝０ ^Ｎｘ−１Σ_ｊｙ＝０ ^Ｎｙ−１Σ_ｊｚ＝０ ^Ｎｚ−１
ｆ_{ｊｘ，ｊｙ，ｊｚ} ｅｘｐ（−２πｉｋ_ｘｊ_ｘ／Ｎ_ｘ） ｅｘｐ（−２πｉｋ_ｙｊ_ｙ／Ｎ_ｙ） ｅｘｐ（−２πｉｋ_ｚｊ_ｚ／Ｎ_ｚ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｋ_ｚ＝０，１，… ，Ｎ_ｚ−１）
さらに、この式は次の３つの変換を順に行うことによって計算できることが容易に示せる。
＜Ｙ方向の変換＞
ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１） ＝Σ_ｊｙ＝０ ^Ｎｙ−１ ｆ_{ｊｘ，ｊｙ，ｊｚ} ｅｘｐ（−２πｉｋ_ｙｊ_ｙ／Ｎ_ｙ）
（ただしｊ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
＜Ｘ方向の変換＞
ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） ＝Σ_ｊｘ＝０ ^Ｎｘ−１ ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１） ｅｘｐ（−２πｉｋ_ｘｊ_ｘ／Ｎ_ｘ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
＜Ｚ方向の変換＞
ｃ_{ｋｘ，ｋｙ，ｋｚ} ＝Σ_ｊｚ＝０ ^Ｎｚ−１ ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） ｅｘｐ（−２πｉｋ_ｚｊ_ｚ／Ｎ_ｚ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｋ_ｚ＝０，１，… ，Ｎ_ｚ−１）
これらの式から明らかなように、Ｙ方向の変換では、添字ｊ_ｘ，ｊ_ｚが同じＮ_ｙ個のデータに対し、１次元のフーリエ変換を行う。そして、ｊ_ｘとｊ_ｚとを動かしてこのような変換をＮ_ｘ×Ｎ_ｚ回行うことにより、Ｙ方向の変換が完了する。Ｘ方向、Ｚ方向の変換についても同様である。したがって、図２の１に示すように、３次元のデータ｛ｆ_{ｊｘ，ｊｙ，ｊｚ}｝を各辺の長さがＮ_ｘ×Ｎ_ｙ×Ｎ_ｚの直方体状に並べた場合、 Ｙ方向の変換ではＹ軸に平行なＮ_ｙ個のデータ２に対して１次元のフーリエ変換を行い、Ｘ方向の変換ではＸ軸に平行なＮ_ｘ個のデータ３、Ｚ方向の変換ではＺ軸に平行なＮ_ｚ個のデータ４に対してそれぞれ１次元のフーリエ変換を行うことになる。この計算方式より、Ｙ方向の変換においてはＸ座標あるいはＺ座標が異なる組に対する計算は並列に行えることが明らかである。また、Ｘ方向の変換においてはＹ座標あるいはＺ座標が異なる組に対する計算が並列に行えること、Ｚ方向の変換においてはＸ座標あるいはＹ座標が異なる組に対する計算が並列に行えることも明らかである。
【０００４】
従来、３次元の高速フーリエ変換を分散メモリ型並列計算機で行うに当たっては、この並列性を利用する方式が一般的であった。その中でも最もプロセッサ間でのデータ転送量が少なく効率的な方式は転置アルゴリズムと呼ばれ、たとえばＶ． Ｋｕｍａｒ， Ａ． Ｇｒａｍａ， Ａ． Ｇｕｐｔａ ａｎｄ Ｇ． Ｋａｒｙｐｉｓ： ”Ｉｎｔｒｏｄｕｃｔｉｏｎ ｔｏ Ｐａｒａｌｌｅｌ Ｃｏｍｐｕｔｉｎｇ”， Ｔｈｅ Ｂｅｎｊａｍｉｎ／Ｃｕｍｍｉｎｇｓ Ｐｕｂｌｉｓｈｉｎｇ Ｃｏｍｐａｎｙ， １９９４， ｐｐ．３７７〜４０６（非特許文献２）に詳しく解説されている。この方式では、図３に示すように、まず３次元のデータをＺ軸に垂直な平面でプロセッサ台数に等しい数の部分データ５に分割し、このデータ５をそれぞれ１台のプロセッサのメモリに格納する。そして、この状態でＹ方向の変換を行う。この状態では、Ｙ方向の変換に必要なＮ_ｙ個のデータ２は１台のプロセッサが持っているため、プロセッサ間でデータ転送を行うことなく変換が実行できることが明らかである。Ｙ方向の変換の終了後、データの分割方式を変更し、３次元のデータをＹ軸に垂直な平面でプロセッサ台数に等しい数の部分データ６に分割し、このデータ６をそれぞれ１台のプロセッサのメモリに格納する。このためには、全プロセッサが自分以外の全プロセッサとデータの転送を行う処理が必要である。この処理を転置と呼ぶ。転置の終了後には、各プロセッサはＸ方向の変換のために必要なＮ_ｘ個のデータ３をすべて自分で持っているため、プロセッサ間でデータ転送を行うことなくＸ方向の変換を行うことができる。また、Ｚ方向の変換においても、各プロセッサは変換のために必要なＮ_ｚ個のデータ４をすべて自分で持っているため、プロセッサ間でデータ転送を行うことなく実行でき、これによって３次元のフーリエ変換が終了する。以上が分散メモリ型並列計算機による従来の３次元高速フーリエ変換の実行方法である。
【非特許文献１】
Ｇ． Ｇｏｌｕｂ ａｎｄ Ｃ． Ｆ． ｖａｎ Ｌｏａｎ： ”Ｍａｔｒｉｘ Ｃｏｍｐｕｔａｔｉｏｎｓ”， ３^ｒｄ ｅｄｉｔｉｏｎ， Ｔｈｅ Ｊｏｈｎｓ Ｈｏｐｋｉｎｓ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ， １９９６， ｐｐ．１８９〜１９２
【非特許文献２】
Ｖ． Ｋｕｍａｒ， Ａ． Ｇｒａｍａ， Ａ． Ｇｕｐｔａ ａｎｄ Ｇ． Ｋａｒｙｐｉｓ： ”Ｉｎｔｒｏｄｕｃｔｉｏｎ ｔｏ Ｐａｒａｌｌｅｌ Ｃｏｍｐｕｔｉｎｇ”， Ｔｈｅ Ｂｅｎｊａｍｉｎ／Ｃｕｍｍｉｎｇｓ Ｐｕｂｌｉｓｈｉｎｇ Ｃｏｍｐａｎｙ， １９９４， ｐｐ．３７７〜４０６
【０００５】
【発明が解決しようとする課題】
上記の転置アルゴリズムによる並列計算方法では、Ｙ方向、Ｘ方向、Ｚ方向の各変換が各プロセッサで完全に独立に実行できる。しかし、途中の転置処理においては、全プロセッサ対全プロセッサのデータ転送が必要となる。一般に、分散メモリ型並列計算機では、データ転送は演算に比べて多くの時間を必要とし、この傾向は近年のプロセッサの高速化に伴ってますます強まりつつある。また最近では、分散メモリ型並列計算機の一種として、多数のパーソナルコンピュータ（ＰＣ）をイーサネット（登録商標）などのネットワークで結合したＰＣクラスタが広く使われるようになってきている。この場合は、専用の並列計算機に比べてプロセッサ間でのデータ転送能力が低いため、特にデータ転送処理が処理時間のネックとなりやすい。このような背景から、分散メモリ型並列計算機で３次元の高速フーリエ変換を行う場合、従来の転置アルゴリズムによる方式では十分な並列性能が得られない場合が多数生じてきている。本発明では、この問題を解決することを目的とする。
【０００６】
【課題を解決するための手段】
上記目的を達成するため、本発明では、転置アルゴリズムにおけるＹ方向とＸ方向の変換をそれぞれ２つの部分処理に分割して演算の順序を変更し、通信と演算のオーバーラップを可能にする。そのため、変換に関する次の性質に着目する。
▲１▼ Ｙ方向の変換では、Ｘ方向に関して並列性がある。したがって、各プロセッサは、自分の行うべきＮ_Ｘ組の変換を、任意の順序で行うことができる。
▲２▼ Ｘ方向の変換はｌｏｇ_２ Ｎ_Ｘ個のステップからなるが、このうちで、最初のｌｏｇ_２Ｎ_Ｘ−１個のステップでは、Ｘ座標が偶数の要素は偶数の要素同士、奇数の要素は奇数の要素同士でしか演算を行わない。
Ｎ_Ｘ＝８の場合を例に取り、Ｘ方向の変換における１組のデータに対する計算の流れを図４に示す。ここで、縦に並ぶ８個の丸７は各ステップでのデータの１要素を表し、丸と丸を結ぶ線８は、右側の丸の値を計算するために左側の丸の値を使うことを示す。図より明らかなように、変換のうち最初のｌｏｇ_２Ｎ_Ｘ１＝２個のステップでは、Ｘ座標が偶数の要素は偶数の要素同士、奇数の要素は奇数の要素同士でしか演算を行わない。したがって、上記で述べた性質▲２▼が成り立っていることがわかる。この性質はＮ_Ｘ＝８の場合に限らず、一般の場合について成り立つ。
上記▲１▼、▲２▼の２つの性質を利用すると、転置アルゴリズムにおけるＹ方向と方向の変換を次の５つの処理に分割して行うことができる。
＜処理１＞
Ｘ座標が偶数の要素のみに対し、Ｙ方向の変換を行う。
＜処理２＞
Ｘ座標が奇数の要素のみに対し、Ｙ方向の変換を行う。
Ｘ座標が偶数の要素に対しては、転置処理を行う。
＜処理３＞
Ｘ座標が偶数の要素のみに対し、Ｘ方向の変換の最初のｌｏｇ_２Ｎ_Ｘ１ステップを行う。
Ｘ座標が奇数の要素に対しては、転置処理を行う。
＜処理４＞
Ｘ座標が奇数の要素のみに対し、Ｘ方向の変換の最初のｌｏｇ_２Ｎ_Ｘ１ステップを行う。
＜処理５＞
Ｘ方向の変換の最後の１ステップを行う。
本発明による処理１〜処理４の様子を図５に示す。図において、斜線部分９および１０は変換の計算を行っている要素、灰色の部分１１はデータ転送を行っている要素、白の部分１２は何も行わない要素を示す。また、直方体中の線はプロセッサへのデータ分割を示す。本発明の計算方式では、処理２においてＸ座標が奇数の要素に対して演算処理を行っている間に、Ｘ座標が偶数の要素に対して転置処理を行うことができる。また、処理３においては、Ｘ座標が偶数の要素に対して演算処理を行っている間に、Ｘ座標が奇数の要素に対して転置処理を行うことができる。これにより、転送のための時間を演算時間に隠蔽することが可能となり、計算効率を向上させることができる。
【０００７】
【発明の実施の形態】
≪実施例１≫
（１） 装置の概略構成
以下、本発明の原理および実施例を、図面により詳細に説明する。ここで実施例として挙げるのは、分散メモリ型並列計算機上で、３次元のフーリエ変換を用いて偏微分方程式を解き、これを用いてシミュレーションを行う場合である。ここでは特に、シミュレーションの例として気象予測のための計算を取り上げる。
本方法における並列プログラムを実行する並列計算機システムを図６に示す。システムは計算領域の形状、初期条件、物質定数などのパラメータを入力するための入力装置１３、それぞれがメモリ１４を備えたＰ台のプロセッサ１５と、各プロセッサの持つメモリの間でデータ転送を行うためのネットワーク１６から成る処理装置１７、計算結果を出力するための出力装置１８、プログラム１９およびデータ２０を格納するための外部記憶装置２１から構成される。また、異なるプロセッサの持つメモリの間でのネットワークを介したデータ転送は、プロセッサの演算処理と同時に実行可能であるとする。
（２） 従来の方法による３次元並列高速フーリエ変換
以下、分散メモリ型並列計算機上での３次元高速フーリエ変換の原理を、数式を用いて説明する。３次元フーリエ変換は、Ｎ_ｘ×Ｎ_ｙ×Ｎ_ｚ個の複素入力データ
｛ｆ_{ｊｘ，ｊｙ，ｊｚ}｝からＮ_ｘ×Ｎ_ｙ×Ｎ_ｚ個の複素出力データ｛ｃ_{ｋｘ，ｋｙ，ｋｚ}｝を、（数１）を用いて計算する処理である。
ｃ_{ｋｘ，ｋｙ，ｋｚ}＝Σ_ｊｘ＝０ ^Ｎｘ−１Σ_ｊｙ＝０ ^Ｎｙ−１Σ_ｊｚ＝０ ^Ｎｚ−１
ｆ_{ｊｘ，ｊｙ，ｊｚ} ｅｘｐ（−２πｉｋ_ｘｊ_ｘ／Ｎ_ｘ） ｅｘｐ（−２πｉｋ_ｙｊ_ｙ／Ｎ_ｙ） ｅｘｐ（−２πｉｋ_ｚｊ_ｚ／Ｎ_ｚ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｋ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数１）
容易にわかるように，この式は次の（数２）（数３）（数４）でそれぞれ示すＹ方向の変換、Ｘ方向の変換、Ｚ方向の変換を順に行うことによって計算することができる。
ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１） ＝Σ_ｊｙ＝０ ^Ｎｙ−１ ｆ_{ｊｘ，ｊｙ，ｊｚ} ｅｘｐ（−２πｉｋ_ｙｊ_ｙ／Ｎ_ｙ）
（ただしｊ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数２）
ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） ＝Σ_ｊｘ＝０ ^Ｎｘ−１ ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１） ｅｘｐ（−２πｉｋ_ｘｊ_ｘ／Ｎ_ｘ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数２）
ｃ_{ｋｘ，ｋｙ，ｋｚ} ＝Σ_ｊｚ＝０ ^Ｎｚ−１ ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） ｅｘｐ（−２πｉｋ_ｚｊ_ｚ／Ｎ_ｚ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｋ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数４）
これらの変換を分散メモリ型並列計算機上で行うには、（数２）がＮ_Ｘ×Ｎ_Ｚ組の独立な変換、（数３）がＮ_Ｙ×Ｎ_Ｚ組の独立な変換、（数４）がＮ_Ｘ×Ｎ_Ｙ組の独立な変換からなることを利用すればよい。この考え方に基づく計算方式が、「従来の技術」の項で説明した転置アルゴリズムである。転置アルゴリズムによる計算のフローチャートを図７に示す。この計算は、主に次の４つの処理からなる。
＜Ｙ方向の変換＞
入力装置から３次元のデータｆ_{ｊｘ，ｊｙ，ｊｚ}を入力し（図７の処理２３）、直方体状に並べ（処理２４）、Ｚ軸に垂直な平面でデータをプロセッサ台数と等しい数に分割し、各部分データをそれぞれ１台のプロセッサのメモリに格納する（処理２５）。データの分割の仕方には、ブロック分割、サイクリック分割など、いろいろな方式がありうるが（データ分割方式についての詳細はＶ． Ｋｕｍａｒ， Ａ． Ｇｒａｍａ， Ａ． Ｇｕｐｔａ ａｎｄ Ｇ． Ｋａｒｙｐｉｓ： ”Ｉｎｔｒｏｄｕｃｔｉｏｎ ｔｏ Ｐａｒａｌｌｅｌ Ｃｏｍｐｕｔｉｎｇ”， Ｔｈｅ Ｂｅｎｊａｍｉｎ／Ｃｕｍｍｉｎｇｓ Ｐｕｂｌｉｓｈｉｎｇ Ｃｏｍｐａｎｙ， １９９４参照。）、たとえばブロック分割の場合、プロセッサｐ（０≦ｐ≦Ｐ−１）はＺ座標ｊ_ｚが（Ｎ_ｚ／Ｐ）×ｐ≦ｊ_ｚ≦（Ｎ_ｚ／Ｐ）×（ｐ＋１）−１をみたすＮ_ｘ×Ｎ_ｙ×Ｎ_ｚ／Ｐ個のデータを持つ。
この状態で、各プロセッサは自分のメモリに格納された部分データに対し、（数２）によるＹ方向の変換を行う（処理２６）。
＜転置処理＞
次に、データを図３の６のようにＹ軸に垂直な面でプロセッサ台数と等しい数に分割し、各部分データをそれぞれ１台のプロセッサのメモリに格納する（処理２７）。この場合もデータ分割の仕方はいろいろな方式がありうるが、たとえばブロック分割の場合、プロセッサｐ（０≦ｐ≦Ｐ−１）はＹ座標ｋ_ｙが（Ｎ_ｙ／Ｐ）×ｐ≦ｋ_ｙ≦（Ｎ_ｙ／Ｐ）×（ｐ＋１）−１をみたすＮ_ｘ×Ｎ_ｙ／Ｐ×Ｎ_ｚ個のデータを持つ。このように分割方式を変更するため、プロセッサｐはＹ方向の変換時に自分が持っていたデータのうち、Ｙ座標ｋ_ｙが（Ｎ_ｙ／Ｐ）×ｐ’≦ｋ_ｙ≦（Ｎ_ｙ／Ｐ）×（ｐ’＋１）−１をみたすＮ_ｘ×Ｎ_ｙ／Ｐ×Ｎ_ｚ／Ｐ個のデータを、プロセッサｐ’に転送する。この処理を転置処理と呼ぶ。
＜Ｘ方向の変換＞
各プロセッサは自分のメモリに格納された部分データに対し、（３）式によるＸ方向の変換を行う（処理２８）。
＜Ｚ方向の変換＞
各プロセッサは自分のメモリに格納された部分データに対し、（４）式によるＺ方向の変換を行う（処理２９）。変換結果のデータを出力装置へ出力する（処理３０）。
【０００８】
以上が従来の転置アルゴリズムによる３次元高速フーリエ変換である。しかしこの計算方式では、「従来の技術」の項で説明したように、途中の転置処理２７において全プロセッサ対全プロセッサのデータ転送が必要となり、この部分が性能上のネックになる。
（３） 本発明の方法による３次元並列高速フーリエ変換
そこで本発明では、上記の転置アルゴリズムにおけるＹ方向とＸ方向の変換をそれぞれ２つの部分処理に分割して演算の順序を変更し、通信と演算のオーバーラップを可能にする。本発明による計算のフローチャートを図１に示す。ただし、本発明による計算は、図７に示した転置アルゴリズムによる計算と、処理２５以前の部分、および処理２９以後の部分は同じであるので、図１のフローチャートでは処理２５より後で、かつ処理２９より前の部分のみを示してある。この部分は、次の５つの処理からなる。
＜Ｙ方向の変換処理（１）＞
第７図の分割ステップ２５で複数プロセッサのメモリに分散格納された部分データのうち、Ｘ座標ｊ_ｘが偶数の要素のみに対し、（２）式によるＹ方向の変換を行う（図１の処理３３）。部分データの分割について再度述べると、変換対象データを直方体に配列した複素データ列についてＺ軸と垂直な面でプロセッサ台数と同じ数に分割する。したがって、上記複素データ列のＸ軸方向及びＹ軸方向は分割されていないので、プロセッサ間のデータ転送なしにＸ方向及びＹ方向での変換処理を完了できる分散配置となる。しかし、Ｙ方向の変換処理（１）では、このうちＸ座標ｊ_ｘが偶数の要素に対するＹ方向の変換のみを各プロセッサにて行う。
＜Ｙ方向の変換処理（２）＞
各プロセッサは、自分のメモリに格納された部分データのうち、Ｘ座標ｊ_ｘが奇数の要素のみに対し、（２）式によるＹ方向の変換を行う。
これと並行して、Ｘ座標ｊ_ｘが偶数の要素に対しては、転置処理を行う。転置処理においてプロセッサｐは、従来の方法での転置処理と同様、自分が持っているデータのうち、Ｙ座標ｋ_ｙが（Ｎ_ｙ／Ｐ）×ｐ’≦ｋ_ｙ≦（Ｎ_ｙ／Ｐ）×（ｐ’＋１）−１をみたすＮ_ｘ／２×Ｎ_ｙ／Ｐ×Ｎ_ｚ個のデータをプロセッサｐ’に転送する（図１の処理３４）。つまりこの転地処理の結果、Ｘ座標ｊ_ｘのが偶数の要素についてのみ、プロセッサ間で再配置される。その結果、プロセッサ間のデータ転送なしに各プロセッサでＺ方向の変換が完了できる状態となる。またＸ方向については、隣接する要素同士の演算であるＸ方向の変換の最終ステップのみを除いたｌｏｇ_２Ｎ_Ｘ１ステップを各プロセッサで実行できる状態となる。
＜Ｘ方向の変換処理（１）＞
各プロセッサは、自分のメモリに格納された部分データのうち、Ｘ座標ｊ_ｘが偶数の要素のみに対し、Ｘ方向の変換の最終ステップのみを除いた最初のｌｏｇ_２Ｎ_Ｘ１ステップを行う。このときの計算式は、（数３）の代わりに
ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２’） ＝Σ_{ｊｘ’＝０} ^{Ｎｘ／２−１} ｃ_{２ｊｘ’，ｋｙ，ｊｚ} ^（１） ｅｘｐ（−２πｉｋ_ｘ２ｊ_ｘ’／Ｎ_ｘ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数５）
とする。
これと並行して、Ｘ座標ｊ_ｘが奇数の要素に対しては、転置処理を行う。転置処理においてプロセッサｐは、従来の方法での転置処理と同様、自分が持っているデータのうち、Ｙ座標ｋ_ｙが（Ｎ_ｙ／Ｐ）×ｐ’≦ｋ_ｙ≦（Ｎ_ｙ／Ｐ）×（ｐ’＋１）−１をみたすＮ_ｘ／２×Ｎ_ｙ／Ｐ×Ｎ_ｚ個のデータをプロセッサｐ’に転送する（図１の処理３５）。これにより、Ｘ座標ｊ_ｘのが奇数の要素についても、各プロセッサでＺ方向の変換が完了でき、またＸ方向については、隣接する要素同士の演算であるＸ方向の変換の最終ステップのみを除いたｌｏｇ_２Ｎ_Ｘ１ステップを各プロセッサで実行できる状態に再配置される。
【０００９】
＜Ｘ方向の変換処理（２）＞
各プロセッサは、自分のメモリに格納された部分データのうち、Ｘ座標ｊ_ｘが奇数の要素のみに対し、Ｘ方向の変換の最終ステップのみを除いた最初のｌｏｇ_２Ｎ_Ｘ１ステップを行う。このときの計算式は、（数３）の代わりに
ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２”） ＝Σ_{ｊｘ’＝０} ^{Ｎｘ／２−１} ｃ_{２ｊｘ’＋１，ｋｙ，ｊｚ} ^（１） ｅｘｐ（−２πｉｋ_ｘ（２ｊ_ｘ’＋１）／Ｎ_ｘ）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数６）
とする（図１の処理３６）。
＜Ｘ方向の変換処理（３）＞
各プロセッサは、自分のメモリに格納されたデータを用いて、Ｘ方向の変換の最後の１ステップを行う。このときの計算式は、（数５）で計算したｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２’）、および（６）式で計算したｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２”）を用いて、
ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） ＝ ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２’） ＋ ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２”）
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−− （数７）
とする（図１の処理３７）。
以上の計算方式では、各プロセッサはのＹ方向の変換処理（２）においてＸ座標が奇数の要素に対してＹ方向の変換処理を行うと同時に、Ｘ座標が偶数の要素に対してデータの転送処理を行うことができる。また、各プロセッサはＸ方向の変換処理（１）において、Ｘ座標が偶数の要素に対してＸ方向の変換処理を行うと同時に、Ｘ座標が奇数の要素に対してデータの転送処理を行うことができる。これにより、転送のための時間の一部または全部を演算時間に隠蔽することが可能となり、計算効率を向上させることができる。
（４） 並列高速フーリエ変換ライブラリ
本発明を並列計算機上で３次元高速フーリエ変換を行うライブラリに適用した例を示す。本ライブラリは、サブルーチン名称をＦＦＴ３Ｄとし、実行するには
ＣＡＬＬ ＦＦＴ３Ｄ（Ｎ_Ｘ， Ｎ_Ｙ， Ｎ_Ｚ， Ｐ， Ｆ， ＴＢ， ＩＯＰＴ， ＩＥＲ）
のようにしてすべてのプロセッサで同時にコールする。ここで、Ｎ_Ｘ，Ｎ_Ｙ，Ｎ_Ｚはフーリエ変換を行うべき３次元データのそれぞれＸ、Ｙ、Ｚ方向の個数、Ｐはプロセッサ台数、Ｆは入力時はフーリエ変換すべき３次元データ｛ｆ_{ｊｘ，ｊｙ，ｊｚ}｝、出力時はフーリエ変換結果｛ｃ_{ｋｘ，ｋｙ，ｋｚ}｝を格納する配列、ＴＢは変換に用いる複素指数関数の値を格納するテーブル、ＩＯＰＴはサブルーチンの機能を指定する入力、ＩＥＲは実行時エラーが生じたか否かを示す出力である。ここで、配列Ｆは各プロセッサがそれぞれ持つ部分配列であり、Ｚ座標に対してブロック分割を行うので、第ｐ番目のプロセッサｐ（０≦ｐ≦Ｐ−１）はＺ座標ｊ_ｚが（Ｎ_ｚ／Ｐ）×ｐ≦ｊ_ｚ≦（Ｎ_ｚ／Ｐ）×（ｐ＋１）−１をみたすＮ_ｘ×Ｎ_ｙ×Ｎ_ｚ／Ｐ個のデータのみを持つ。すなわち、第ｐ番目のプロセッサの配列Ｆには、
Ｆ（ｊ_ｘ， ｊ_ｙ， ｊ_ｚ’） ＝ ｆ_{ｊｘ，ｊｙ，（Ｎｚ／Ｐ）＊ｐ＋ｊｚ’}
（ただしｊ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｊ_ｙ＝０，１，… ，Ｎ_ｙ−１；ｊ_ｚ’＝０，１，… ，Ｎ_ｚ／Ｐ−１）
−−−（数８）
を格納する。したがって各プロセッサの持つ配列Ｆの大きさはＮ_Ｘ×Ｎ_Ｙ×Ｎ_Ｚ／Ｐである。また、ＴＢは、第１回目のコールで計算した複素指数関数の値を格納しておくテーブルであり、２回目のコールからはここに格納した値を再利用することにより、新たな計算が不要となる。また、第１回目のコールではＩＯＰＴ＝１を指定し、このときは複素指数関数のテーブルを作成する。ＩＯＰＴ＝２は２回目以降のコールを意味し、このときは既にＴＢに格納されている値を用いる。
【００１０】
本ライブラリのフローチャートを図８に示す。ライブラリは、コールされるとまず入力データをチェックし、Ｎ_Ｘ，Ｎ_Ｙ，Ｎ_ＺとＰとが１以上の整数であるかどうか、ＩＯＰＴが１または２の値であるかどうかなど、入力の有効性を調べる（処理４０）。入力データに無効な値が入っていた場合は、ＩＥＲ＝１０００と設定して（処理４１）リターンする。次に、ライブラリが引数で指定した通りにＰ台のプロセッサでコールされているかどうかをチェックする（処理４３）。この条件が満たされていない場合は、ＩＥＲ＝２０００と設定して（処理４４）リターンする。次にＩＯＰＴの値をチェックし（処理４５）、ＩＯＰＴ＝１の場合は各方向のフーリエ変換で用いる複素指数関数の値を前もって計算し、配列ＴＢに格納する（処理４６）。
【００１１】
次に、Ｙ方向およびＸ方向の変換を行う（処理４７）。この変換は、本実施例における「（３） 本発明の方法による３次元並列高速フーリエ変換」の項で述べた方法に従い、Ｙ方向の変換処理（１）（図１の処理３３）からＸ方向の変換処理（３）（図１の処理３７）の５つを順に行うことによって終了する。次に、本実施例における「（２） 従来の方法による３次元並列高速フーリエ変換」の項で述べた方法に従い、Ｚ方向の変換を行う（図８の処理４８）。これにより、３次元高速フーリエ変換の処理が終了する。終了時には、従来の方法による３次元並列高速フーリエ変換と同様、配列ＦにはＹ座標に対してブロック分割されたデータが格納される。すなわち、第ｐ番目のプロセッサｐの配列Ｆには、
Ｆ（ｋ_ｘ， ｋ_ｙ’， ｋ_ｚ） ＝ ｃ_{ｋｘ，（Ｎｙ／Ｐ）＊ｐ＋ｋｙ’，ｋｚ}
（ただしｋ_ｘ＝０，１，… ，Ｎ_ｘ−１；ｋ_ｙ’＝０，１，…，Ｎ_ｙ／Ｐ−１；ｋ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−−（数９）
が格納される（図８の処理４９）。
（５） シミュレーションプログラム
本実施例において実行すべき気象計算のための並列プログラムを図９に示す。ここでは、サイズＮ_ｘ×Ｎ_ｙ×Ｎ_ｚの３次元メッシュで計算を行う場合を例にとって説明する。
【００１２】
本プログラムでは、まず計算領域のサイズＮ_ｘ，Ｎ_ｙ，Ｎ_ｚ、温度・風速・圧力などの初期条件、空気の熱伝導率をはじめとする物質定数などのパラメータを入力し（処理５１）、計算に必要な前処理を行う。ここで前処理とは、観測によって得られた温度・風速・圧力などのデータに対して補間を行い、計算に必要なメッシュポイントでのデータを得ることである。次に、温度・風速・圧力などのデータを並列計算機の各プロセッサに分割する（処理５２）。分割は、前項で述べた３次元高速フーリエ変換ライブラリＦＦＴ３Ｄを利用できるように、Ｚ方向をブロック分割する形で行う。
【００１３】
これらの処理が終わった後、ループにより各時間ステップでの温度・風速・圧力などの量を順々に求めていく。基礎となる方程式は、風速に対する運動方程式
ｄｕ／ｄｔ ＝ −２Ω×ｕ − （１／ρ）∇ｐ ＋ Ｆ_ｕ −−−（数１０）
質量保存の式
ｄρ／ｄｔ ＝ −ρ∇・ｕ −−−（数１１）
温度変化を表す式
ｄＴ／ｄｔ ＝ −κ∇^２Ｔ ＋ ｕ・∇Ｔ −−−（数１２）
の３本である。ここで、ｕは風速、ｐは圧力、Ｔは温度を表し、Ωはコリオリ力と呼ばれる地球の自転による力、 Ｆ_ｕはそれ以外の外力、ρは空気の密度、κは空気の熱伝導率を表す。これらの式から次の時刻でのデータの値を求めるには、まずフーリエ変換により格子点上の温度Ｔ、圧力ｐおよび風速ｕを波数空間でのデータに変換し（処理５３）、波数空間でこれらのデータを微分し（処理５４）、再び波数空間でのデータを逆フーリエ変換して格子点上での温度勾配∇Ｔ 、２次微分∇^２Ｔ 、圧力勾配∇ｐ 、速度の発散∇・ｕ を求める（処理５５）。この後、これらの量を（８）−（１０）式の右辺に代入し、次の時間ステップでの温度・風速・圧力を求める（処理５６）。なお、フーリエ変換により格子点上のデータを波数空間上のデータに変換して処理を行うのは、その方が微分が精度良く計算できるからであり、本プログラムではこの部分で２次元フーリエ変換ライブラリＦＦＴ３Ｄを用いる。
上記のループでは、各時間ステップ毎に、求める時刻までの計算が終了したかどうかを判定し（処理５７）、終了したら、後処理を行い（処理５８）、結果を出力する（処理５９）。後処理では主に、計算を行うメッシュポイントと結果データが必要な点とがずれている場合に、結果を補間して必要な点での値を計算するなどの処理を行う。
上記の例では気象予測計算を行う場合を例にとって説明したが、本発明の手法は、これ以外の応用例についても、並列計算機上で３次元高速フーリエ変換を用いてシミュレーションを行う場合に適用できることは明らかである。また、上記の例ではフーリエ変換前には３次元データをＺ方向にブロック分割し、変換後にはデータはＹ方向にブロック分割されているとしたが、ブロック分割の代わりに、サイクリック分割やブロックサイクリック分割を用いた場合でも、本発明の手法が適用できることは明らかである。更に、上記の例ではＹ方向、Ｘ方向、Ｚ方向の変換をこの順に行うとしたが、Ｘ、Ｙ、Ｚなどは便宜的に定めた座標軸の名前に過ぎないため、上記の計算方法においてＸ→Ｙ、Ｙ→Ｚ、Ｚ→Ｘなどの名前の付け替えを行った計算方法も、本発明の手法と同じものであることは明らかである。
≪実施例２≫
上記実施例１では３次元フーリエ変換のための計算方法を説明したが、本方式は１次元フーリエ変換の場合へも応用できる。Ｎ点のデータに対して１次元フーリエ変換を行うには、 Ｎ_ｘ×Ｎ_ｙ×Ｎ_ｚ＝Ｎを満たす任意の整数Ｎ_ｘ，Ｎ_ｙ，Ｎ_ｚを定めてデータをＮ_ｘ×Ｎ_ｙ×Ｎ_ｚの直方体状に並べ、これに対して「ひねり係数の乗算」という処理を加えた３次元フーリエ変換を行えば良いことが知られている。ここで、ひねり係数の乗算とは、（２）式によってＹ方向の変換を行った後、中間結果のｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１） に対して
ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１） ：＝ ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１）×ｅｘｐ（２πｉｋ_ｙｊ_ｘ／（Ｎ_ｘＮ_ｙ）
（ただしｊ_ｘ＝０，１，… ，Ｎ_ｘ−１，ｋ_ｙ＝０，１，… ，Ｎ_ｙ−１，ｊ_ｚ＝０，１，… ，Ｎ_ｚ−１）
−−−（数１３）
という処理を行い、かつ（数３）によってＸ方向の変換を行った後に中間結果のｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） に対して
＜Ｘ方向の変換＞
ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２） ＝ ｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２）×ｅｘｐ（２πｉ（Ｎ_ｙｋ_ｘ＋ｋ_ｙ）ｊ_ｘ／（Ｎ_ｘＮ_ｙＮ_ｚ））
−−−（数１４）
という処理を行うことである。このひねり係数の乗算は、配列ｃ_{ｊｘ，ｋｙ，ｊｚ} ^（１）またはｃ_{ｋｘ，ｋｙ，ｊｚ} ^（２）に対する要素毎に独立の処理であるので、実施例１の方式においても、Ｙ方向の変換後とＸ方向の変換後にそれぞれ組み入れることができる。これにより、本発明の方式を１次元高速フーリエ変換に対しても適用することが可能となる。
≪実施例３≫
本発明は更に、２次元の高速フーリエ変換に対しても適用できる。この場合は、実施例１の方式による３次元高速フーリエ変換において、単にＮ_ｚ＝１とすればよい。
【００１４】
補足すると、変換対象のデータはＸ軸とＹ軸の２次元に配列された複素データアレイである。実施例１で説明した「Ｙ方向の変換処理（１）」の前のデータ分割の段階では、上記複素データアレイをＸ軸に垂直な直線でプロセッサ台数と等しい数の部分データに分割し、各プロセッサのメモリに分散配置することになる。つまり、プロセッサ間のデータ転送なしにＹ方向の変換が完了できる状態である。しかし「Ｙ方向の変換処理（１）」では、Ｘ座標ｊ_ｘが偶数の要素に対するＹ方向の変換のみを各プロセッサにて行う。次に、Ｙ方向の変換処理（２）でＸ座標ｊ_ｘが奇数の要素に対するＹ方向の変換を各プロセッサにて行うのと並行して、Ｘ座標ｊ_ｘが偶数の要素はプロセッサ間で再配分を行う。この再配分とは、Ｘ座標ｊ_ｘが偶数の要素をＹ軸と垂直な直線でプロセッサ台数と等しい数の部分データに分割して各プロセッサに分散配置する処理であり、この結果、各プロセッサではＸ方向の変換の最終ステップを除いたｌｏｇ_２Ｎ_Ｘ１ステップが実行可能となる。さらにＸ方向の変換処理（１）では、再配分されたＸ座標ｊ_ｘが偶数の要素に対するＸ方向の変換（ｌｏｇ_２Ｎ_Ｘ１ステップ）を各プロセッサにて行うのと並行して、Ｘ座標ｊ_ｘが奇数の要素のプロセッサ間での再配分を行う。この再配分も上記と同様である。更に、Ｘ方向の変換処理（２）、 Ｘ方向の変換処理（３）と進むのは実施例１と全く同様である。変換対象のデータが２次元データであるので、Ｘ方向の変換処理（３）で処理が完了する。
≪実施例４≫
本発明による高速フーリエ変換を用いてシミュレーションを行う他の例として、半導体デバイス等における電子構造計算を説明する。電子構造計算では、３次元のメッシュで定義された電子の波動関数ｕ（ｒ）を、シュレディンガー方程式
ｄｕ（ｒ）／ｄｔ ＝ −（ｈ^２／２ｍ）∇^２ｕ（ｒ） ＋ （Ｅ−Ｖ（ｒ））ｕ（ｒ）−−−（数１５）
に従って計算することにより、半導体の性質を決定するバンドギャップの大きさや、結晶の構造安定性などを求める。ただし、上式でｈはプランク定数、ｍは電子の質量、Ｅは対象とする波動関数のエネルギーレベル、Ｖは結晶中の原子や他の電子によるポテンシャルエネルギーを表す。
【００１５】
（１５）式の計算では、波動関数の２次微分∇^２ｕ（ｒ）が必要であるが、気象計算の例において述べたのと同様な理由により、この部分はｕ（ｒ）をフーリエ変換により波数空間に移してから計算し、結果を逆フーリエ変換で再び実空間に戻す。したがって、並列計算機上で電子構造計算を行う場合には、この部分で本発明の３次元高速フーリエ変換方法が適用できる。
【００１６】
【発明の効果】
以上説明したように、本発明によれば、分散メモリ型並列計算機上での高速フーリエ変換において、演算処理とデータ転送処理を並行して行い、後者の時間の一部または全部を前者の時間に隠蔽することが可能となる。したがって本発明の方法では、従来法に比較して並列化効率を向上させることができる。なお、この向上効果は、本発明を適用する分散メモリ型並列計算機のプロセッサ間通信の性能に依存するが、たとえば１６台のプロセッサでＮ＝２５６^３の問題を処理する場合、２０％〜３０％程度の実行時間短縮が期待できる。
【図面の簡単な説明】
【図１】本発明による並列３次元高速フーリエ変換方法におけるＹ方向及びＸ方向の変換の計算方法を示すフローチャート。
【図２】従来法による３次元高速フーリエ変換方法を示す図。
【図３】従来法による並列３次元高速フーリエ変換方法を示す図。
【図４】Ｘ方向の変換におけるデータの依存関係を示す図。
【図５】本発明による並列３次元高速フーリエ変換方法におけるＹ方向及びＸ方向の変換の計算方法を示す図。
【図６】本発明を適用すべき分散メモリ型並列計算機の構成を示す図。
【図７】従来法による並列３次元高速フーリエ変換方法を示すフローチャート。
【図８】本発明による並列３次元高速フーリエ変換ライブラリのフローチャート。
【図９】本発明を利用して分散メモリ型並列計算機で気象計算を行う場合のフローチャート。
【符号の説明】
１：ＮＸ×ＮＹ×ＮＺの直方体状に並べたデータ、２：Ｙ軸に平行なＮＹ個のデータ、３：Ｘ軸に平行なＮＸ個のデータ、４：Ｚ軸に平行なＮＺ個のデータ、５：Ｚ軸に垂直な平面で分割した部分データ、６：Ｙ軸に垂直な平面で分割した部分データ、７：Ｘ方向の変換における各ステップでのデータの１要素、８：右側の要素の値を計算するために左側の要素の値を使うことを示す線、９：Ｙ方向の変換を行っている要素、１０：Ｘ方向の変換を行っている要素、１１：データ転送を行っている要素、１２：何も行わない要素、１３：入力装置、１４：メモリ、１５：プロセッサ、１６：ネットワーク、１７：処理装置、１８：出力装置、１９：プログラム、２０：データ、２１：外部記憶装置、２２：スタート、２３：変換すべきデータの入力、２４：入力データを直方体状に並べる、２５：データをＺ軸に垂直な面でプロセッサ間に分割、２６：Ｙ方向の変換、２７：データをＹ軸に垂直な面で再分割、２８：Ｘ方向の変換、２９：Ｚ方向の変換、３０：変換結果の出力、３１：終了、３２：図７の処理２５より継続、３３：Ｙ方向の変換処理（１）、３４：Ｙ方向の変換処理（２）、３５：Ｘ方向の変換処理（１）、３６：Ｘ方向の変換処理（２）、３７：Ｘ方向の変換処理（３）、３８：図７の処理２９へ継続、３９：ライブラリコール、４０：入力データの有効性チェック、４１：エラーコードの設定、４２：他のプロセッサにコールを通知、４３：Ｐ台のプロセッサでコールしていることをチェック、４４：エラーコードの設定、４５：ＩＯＰＴの値のチェック、４６：複素指数関数テーブルの設定、４７：Ｙ方向及びＸ方向の変換、４８：Ｚ方向の変換、４９：変換結果の格納、５０：リターン、５１：計算領域の形状・初期条件・パラメータの入力、５２：前処理とデータの分割、５３：温度及び風速データを端数空間でのデータに変換、５４：端数空間で温度勾配及び圧力データを計算、５５：温度勾配及び圧力データを格子点上のデータに変換、５６：次の時間ステップでの温度と風速データの計算、５７：判定処理、５８：物理量の計算、５９：結果の出力。

Claims (10)

1. 入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とを含む並列計算機上で高速フーリエ変換を行う方法であって、
   変換対象となるデータを複数の部分データに分割し、各部分データを前記複数のプロセッサのメモリに分散して格納するステップと、
   各プロセッサで前記複数の部分データをさらに分割した第１部分、第２部分のうちの各第１部分に対して演算処理を行うステップと、
   各プロセッサで前記複数の部分データの各第２部分に対して演算処理を行い、これと並行して前記各第１部分の演算処理結果を前記複数プロセッサ間で再配置するするステップを含むことを特徴とする並列計算機上での高速フーリエ変換方法。
2. 前記複数の部分データの各第１部分に対する演算処理、及び各第２部分に対する演算処理は、いずれもデータを第１軸方向に関してフーリエ変換する処理であり、
   前記各第１部分の演算処理結果の再配分のステップに引き続き、各プロセッサに再配分されたデータを第２軸方向に関してフーリエ変換するステップを更に有する請求項１の高速フーリエ変換方法。
3. 前記複数の部分データの各第１部分とはそれぞれ前記変換対象のデータ配列のうち第２軸方向で偶数番目のデータであり、各第２部分とはそれぞれ前記変換対象のデータ配列のうち第２軸方向で奇数番目のデータであることを特徴とする請求項２の高速フーリエ変換方法。
4. 前記再配分されたデータを第２軸方向に関してフーリエ変換するステップに並行して、前記複数の部分データの各第２部分の第１の方向に関するフーリエ変換の結果を各プロセッサに再配分することを特徴とする請求項３の高速フーリエ変換方法。
5. 前記各第２部分の第１軸方向に関するフーリエ変換の結果の再配分のステップの終了後に、再配分された各第２部分のデータを第２軸方向に関してフーリエ変換するステップと、
   該各第２部分の第２軸方向に関するフーリエ変換の結果と、先に得られた各第１部分の第２軸方向に関するフーリエ変換の結果との双方を用いて、前記変換対象のデータの第２軸方向に関するフーリエ変換の最終演算を行うステップとを更に有する請求項４の高速フーリエ変換方法。
6. 入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とを含む並列計算機上で３次元の高速フーリエ変換を行う方法であって、
   データをＸ軸方向、Ｙ軸方向、Ｚ軸方向の長さがそれぞれＮＸ、ＮＹ、ＮＺであるような直方体として並べ、Ｚ方向に垂直な面で分割して各部分データを前記複数のプロセッサのメモリに分散して格納し、
   各プロセッサにてＸ座標が偶数の要素のみに対し、Ｙ方向の変換処理を行い、
   各プロセッサにてＸ座標が奇数の要素のみに対し、Ｙ方向の変換処理を行い、これと並行して、Ｘ座標が偶数の要素に対しては、データがＹ軸に垂直な面で分割されて各部分データが前記複数のプロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
   各プロセッサにてＸ座標が偶数の要素のみに対し、Ｘ方向の変換の最初のｌｏｇ_２ ＮＸ １ ステップを行い、これと並行してＸ座標が偶数の要素に対しては、データがＹ軸に垂直な面で分割されて各部分データが前記複数のプロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
   各プロセッサにてＸ座標が奇数の要素のみに対し、Ｘ方向の変換の最初のｌｏｇ_２ ＮＸ １ ステップを行い、
   各プロセッサにてＸ方向の変換の最後の１ステップを行い、
   各プロセッサにてＺ方向の変換処理を行うことを特徴とする並列計算機上での高速フーリエ変換方法。
7. 入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とから構成される並列計算機上でＮ点のデータに対する１次元の高速フーリエ変換を行う方法であって、
   ＮＸ、ＮＹ、ＮＺをＮＸ×ＮＹ×ＮＺ＝Ｎを満たす整数とするとき、変換対象のデータをＸ軸方向、Ｙ軸方向、Ｚ軸方向の長さがそれぞれＮＸ、ＮＹ、ＮＺであるような直方体として並べ、Ｚ方向に垂直な面でデータを分割し、各部分データを前記複数のプロセッサのメモリに分散して格納し、
   各プロセッサにてＸ座標が偶数の要素のみに対し、Ｙ方向の変換処理とひねり係数の乗算処理を行い、
   各プロセッサにてＸ座標が奇数の要素のみに対し、Ｙ方向の変換処理とひねり係数の乗算処理を行い、これと並行してＸ座標が偶数の要素に対しては、データがＹ軸に垂直な面で分割されて各部分データが該並列計算機の各プロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
   各プロセッサにてＸ座標が偶数の要素のみに対し、Ｘ方向の変換の最初のｌｏｇ_２ ＮＸ １ ステップを行い、これと並行してＸ座標が偶数の要素に対しては、データがＹ軸に垂直な面で分割されて各部分データが該並列計算機の各プロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
   各プロセッサにてＸ座標が奇数の要素のみに対し、Ｘ方向の変換の最初のｌｏｇ_２ ＮＸ １ ステップを行い、
   各プロセッサにおいてＸ方向の変換の最後の１ステップとひねり係数の乗算処理を行い、
   各プロセッサにおいてＺ方向の変換処理を行うことを特徴とする並列計算機上での高速フーリエ変換方法。
8. 入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とを含む並列計算機上で高速フーリエ変換を行う方法であって、
   該並列計算機の各プロセッサの演算の状態及びプロセッサ間の通信の状態を表示するモニタで観察した場合に、処理が４つのフェーズから構成され、
   第１のフェーズでは各プロセッサは演算のみを行い、
   第２のフェーズでは各プロセッサは演算を行うと同時に、自分以外の全てのプロセッサに対してデータの転送処理を行い、
   第３のフェーズでは各プロセッサは演算を行うと同時に、自分以外の全てのプロセッサに対してデータの転送処理を行い、
   第４のフェーズでは各プロセッサは演算のみを行うことを特徴とする並列計算機上での高速フーリエ変換方法。
9. 並列計算機を用いて高速フーリエ変換を利用して気象予測計算を行うシステムであって、該高速フーリエ変換の部分に請求項１、請求項６、請求項７もしくは請求項８のいずれか一つに記載の方法を使うことを特徴とする並列計算機上での気象予測計算システム。
10. 並列計算機を用いて高速フーリエ変換を利用して電子構造計算を行うシステムであって、該高速フーリエ変換の部分に請求項１、請求項６、請求項７もしくは請求項８のいずれか一つに記載の方法を使うことを特徴とする並列計算機上での電子構造計算システム。

Images