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