(57)【特許請求の範囲】
【請求項１】
入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とを含む並列計算機上で高速フーリエ変換を行う方法であって、
変換対象となるデータを複数の部分データに分割し、各部分データを前記複数のプロセッサのメモリに分散して格納するステップと、
各プロセッサで前記複数の部分データをさらに分割した第１部分、第２部分のうちの各第１部分に対して演算処理を行うステップと、
各プロセッサで前記複数の部分データの各第２部分に対して演算処理を行い、これと並行して前記各第１部分の演算処理結果を前記複数プロセッサ間で再配分するステップとを含み、
前記複数の部分データの各第１部分に対する演算処理、及び各第２部分に対する演算処理は、いずれもデータを第１軸方向に関してフーリエ変換する処理であり、
前記各第１部分の演算処理結果の再配分のステップに引き続き、各プロセッサに再配分されたデータを第２軸方向に関してフーリエ変換するステップを更に有し、
前記複数の部分データの各第１部分とはそれぞれ前記変換対象のデータ配列のうち第２軸方向で偶数番目のデータであり、各第２部分とはそれぞれ前記変換対象のデータ配列のうち第２軸方向で奇数番目のデータであることを特徴とする並列計算機上での高速フーリエ変換方法。
【請求項２】
前記再配分されたデータを第２軸方向に関してフーリエ変換するステップに並行して、前記複数の部分データの各第２部分の第１の方向に関するフーリエ変換の結果を各プロセッサに再配分することを特徴とする請求項１記載の並列計算機上での高速フーリエ変換方法。
【請求項３】
前記各第２部分の第１軸方向に関するフーリエ変換の結果の再配分のステップの終了後に、再配分された各第２部分のデータを第２軸方向に関してフーリエ変換するステップと、
該各第２部分の第２軸方向に関するフーリエ変換の結果と、先に得られた各第１部分の第２軸方向に関するフーリエ変換の結果との双方を用いて、前記変換対象のデータの第２軸方向に関するフーリエ変換の最終演算を行うステップとを更に有する請求項２記載の並列計算機上の高速フーリエ変換方法。
【請求項４】
入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とを含む並列計算機上で３次元の高速フーリエ変換を行う方法であって、
データをX軸方向、Y軸方向、Z軸方向の長さがそれぞれNX、NY、NZであるような直方体として並べ、Z方向に垂直な面で分割して各部分データを前記複数のプロセッサのメモリに分散して格納し、
各プロセッサにてX座標が偶数の要素のみに対し、Y方向の変換処理を行い、
各プロセッサにてX座標が奇数の要素のみに対し、Y方向の変換処理を行い、これと並行して、X座標が偶数の要素に対しては、データがY軸に垂直な面で分割されて各部分データが前記複数のプロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
各プロセッサにてX座標が偶数の要素のみに対し、X方向の変換の前半のlog2 NX-1ステップを行い、これと並行してX座標が偶数の要素に対しては、データがY軸に垂直な面で分割されて各部分データが前記複数のプロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
各プロセッサにてX座標が奇数の要素のみに対し、X方向の変換の前半のlog2 NX-1ステップを行い、
各プロセッサにてX方向の変換の後半の１ステップを行い、
各プロセッサにてZ方向の変換処理を行うことを特徴とする並列計算機上での高速フーリエ変換方法。
【請求項５】
入力装置と、それぞれがメモリを備えた複数のプロセッサ及びそれらのメモリ間でデータの転送を行うネットワークからなる処理装置と、出力装置と、外部記憶装置とから構成される並列計算機上でN点のデータに対する１次元の高速フーリエ変換を行う方法であって、
NX、NY、NZをNX×ＮＹ×NZ=Nを満たす整数とするとき、変換対象のデータをX軸方向、Y軸方向、Z軸方向の長さがそれぞれNX、NY、NZであるような直方体として並べ、Z方向に垂直な面でデータを分割し、各部分データを前記複数のプロセッサのメモリに分散して格納し、
各プロセッサにてX座標が偶数の要素のみに対し、Y方向の変換処理とひねり係数の乗算処理を行い、
各プロセッサにてX座標が奇数の要素のみに対し、Y方向の変換処理とひねり係数の乗算処理を行い、これと並行してX座標が偶数の要素に対しては、データがY軸に垂直な面で分割されて各部分データが該並列計算機の各プロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
各プロセッサにてX座標が偶数の要素のみに対し、X方向の変換の前半のlog2 NX-1ステップを行い、これと並行してX座標が偶数の要素に対しては、データがY軸に垂直な面で分割されて各部分データが該並列計算機の各プロセッサのメモリに分散して格納されるよう、データの転送処理を行い、
各プロセッサにてX座標が奇数の要素のみに対し、X方向の変換の前半のlog2 NX-1ステップを行い、
各プロセッサにおいてX方向の変換の後半の１ステップとひねり係数の乗算処理を行い、
各プロセッサにおいてZ方向の変換処理を行うことを特徴とする並列計算機上での高速フーリエ変換方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は分散メモリ型の並列計算機を用いて高速フーリエ変換を行うための方式に関する。
【０００２】
【従来の技術】
半導体デバイスの特性計算、電子状態計算、気象予測のための計算などの科学技術計算では、数万から数百万に上る変数を扱う大規模シミュレーションが必要である。このような大規模問題を扱う手段としては、並列計算機、特に分散メモリ型と呼ばれるタイプの並列計算機が有力である。分散メモリ型並列計算機は、各々がメモリを持つプロセッサをネットワークで複数個結合したシステムであり、従来の逐次型計算機に比べ、プロセッサ台数を増やすことでピーク性能をいくらでも高めることができるという利点を持つ。
分散メモリ型並列計算機では、計算対象のデータを各プロセッサの持つメモリに分散させて格納し、そのそれぞれに対して各プロセッサが並列に計算を行う。もし計算途中で他のプロセッサの持つデータが必要となる場合は、必要なデータをそのプロセッサから転送してもらい、その後に計算を続行する。したがって分散メモリ型並列計算機では、一般に、演算時間に加えてプロセッサ間でデータを転送するための時間が必要である。そこで、計算効率を上げるには、プロセッサ間での通信がなるべく少なくて済むような並列性の高い計算手法が必要である。また、多くの分散メモリ型並列計算機においては、あるデータに対して演算を行っている間に、別のデータを他プロセッサに転送する機構が用意されている。この場合は、演算とデータ転送が同時に行えるように計算手法を工夫すれば、転送のための時間が演算時間に隠蔽でき、計算効率を向上させることができる。
【０００３】
科学技術計算において頻繁に利用される処理の一つに、フーリエ変換がある。フーリエ変換は、ある実数区間で定義された複素数値をとる関数f(x)を複素指数関数exp(ikx)の重ね合わせとして表す処理であり、計算機上で実現する場合には、扱いうる点の数が有限であることから、複素数の点列f0、f1，… ，fN-1をN個の複素指数関数
exp(2πikj/N)（ただしk=0，1，… ，N-1で、iは虚数単位、πは円周率）
の重ね合わせとして
fj=Σk=0N-1 ck exp(2πikj/N)  （ただしj=0，1，… ，N-1）
と表す処理となる。すなわち、 f0，f1，… ，fN-1が与えられたときに、重ね合わせの係数c0，c1，… ，cN-1を求める処理がフーリエ変換である。この係数c0，c1，… ，cN-1は、式
ck=(1/N)Σj=0N-1 fj exp(-2πikj/N)  （ただしk=0，1，… ，N-1）
により計算できることが知られているが、この定義に基づいて計算を行うと、式の数がN本あり、各式がN個の項から成るため、複素指数関数 exp(-2πikj/N)の計算に加えて、複素数の加算と乗算が約N2回必要である。そこで実際には、アルゴリズム上の工夫により計算量を約NlogNのオーダーに減少させた高速フーリエ変換という手法が広く使われている。高速フーリエ変換の計算法については、たとえばG. Golub and C. F. van Loan: "Matrix Computations", 3rd edition, The Johns Hopkins University Press, 1996, pp.189〜192（非特許文献１）に詳しく記述されている。
以上で述べたフーリエ変換は１次元のフーリエ変換と呼ばれるが、半導体デバイスの特性計算、電子状態計算、気象予測のための計算などでは、３次元のフーリエ変換が用いられる。これは、３つの添字jx，jy，jzを持つ複素数データ
{fjx,jy,jz}（ただしjx=0，1，… ，Nx-1，jy=0，1，… ，Ny-1，jz=0，1，… ，Nz-1）をNx×Ny×Nz個の複素指数関数
exp(2πikxjx/Nx) exp(2πikyjy/Ny) exp(2πikzjz/Nz)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；kz=0，1，… ，Nz-1）
の重ね合わせとして
fjx,jy,jz=Σkx=0Nx-1Σky=0Ny-1Σkz=0Nz-1
ckx,ky,kz exp(2πikxjx/Nx) exp(2πikyjy/Ny) exp(2πikzjz/Nz)
（ただしjx=0，1，… ，Nx-1；jy=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
と表す処理となる。すなわち、{fjx,jy,jz}が与えられたときに、重ね合わせの係数{ckx,ky,kz}を求める処理が３次元フーリエ変換である。この係数{ckx,ky,kz}は、次の式によって計算できることが知られている。
ckx,ky,kz=Σjx=0Nx-1Σjy=0Ny-1Σjz=0Nz-1
fjx,jy,jz exp(-2πikxjx/Nx) exp(-2πikyjy/Ny) exp(-2πikzjz/Nz)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；kz=0，1，… ，Nz-1）
さらに、この式は次の３つの変換を順に行うことによって計算できることが容易に示せる。
＜Y方向の変換＞
cjx,ky,jz(1) =Σjy=0Ny-1 fjx,jy,jz exp(-2πikyjy/Ny)
（ただしjx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
＜X方向の変換＞
ckx,ky,jz(2) =Σjx=0Nx-1 cjx,ky,jz(1) exp(-2πikxjx/Nx)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
＜Z方向の変換＞
ckx,ky,kz =Σjz=0Nz-1 ckx,ky,jz(2) exp(-2πikzjz/Nz)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；kz=0，1，… ，Nz-1）
これらの式から明らかなように、Y方向の変換では、添字jx，jzが同じNy個のデータに対し、１次元のフーリエ変換を行う。そして、jxとjzとを動かしてこのような変換をNx×Nz回行うことにより、Y方向の変換が完了する。X方向、Z方向の変換についても同様である。したがって、図２の1に示すように、３次元のデータ{fjx,jy,jz}を各辺の長さがNx×Ny×Nzの直方体状に並べた場合、 Y方向の変換ではY軸に平行なNy個のデータ2に対して１次元のフーリエ変換を行い、X方向の変換ではX軸に平行なNx個のデータ3、Z方向の変換ではZ軸に平行なNz個のデータ4に対してそれぞれ１次元のフーリエ変換を行うことになる。この計算方式より、Y方向の変換においてはX座標あるいはZ座標が異なる組に対する計算は並列に行えることが明らかである。また、X方向の変換においてはY座標あるいはZ座標が異なる組に対する計算が並列に行えること、Z方向の変換においてはX座標あるいはY座標が異なる組に対する計算が並列に行えることも明らかである。
【０００４】
従来、３次元の高速フーリエ変換を分散メモリ型並列計算機で行うに当たっては、この並列性を利用する方式が一般的であった。その中でも最もプロセッサ間でのデータ転送量が少なく効率的な方式は転置アルゴリズムと呼ばれ、たとえばV. Kumar, A. Grama, A. Gupta and G. Karypis: "Introduction to Parallel Computing", The Benjamin/Cummings Publishing Company, 1994, pp.377〜406（非特許文献２）に詳しく解説されている。この方式では、図３に示すように、まず３次元のデータをZ軸に垂直な平面でプロセッサ台数に等しい数の部分データ5に分割し、このデータ5をそれぞれ１台のプロセッサのメモリに格納する。そして、この状態でY方向の変換を行う。この状態では、Y方向の変換に必要なNy個のデータ2は１台のプロセッサが持っているため、プロセッサ間でデータ転送を行うことなく変換が実行できることが明らかである。Y方向の変換の終了後、データの分割方式を変更し、３次元のデータをY軸に垂直な平面でプロセッサ台数に等しい数の部分データ6に分割し、このデータ6をそれぞれ１台のプロセッサのメモリに格納する。このためには、全プロセッサが自分以外の全プロセッサとデータの転送を行う処理が必要である。この処理を転置と呼ぶ。転置の終了後には、各プロセッサはX方向の変換のために必要なNx個のデータ3をすべて自分で持っているため、プロセッサ間でデータ転送を行うことなくX方向の変換を行うことができる。また、Z方向の変換においても、各プロセッサは変換のために必要なNz個のデータ4をすべて自分で持っているため、プロセッサ間でデータ転送を行うことなく実行でき、これによって３次元のフーリエ変換が終了する。以上が分散メモリ型並列計算機による従来の３次元高速フーリエ変換の実行方法である。
【非特許文献１】
G. Golub and C. F. van Loan: "Matrix Computations", 3rd edition, The Johns Hopkins University Press, 1996, pp.189〜192
【非特許文献２】
V. Kumar, A. Grama, A. Gupta and G. Karypis: "Introduction to Parallel Computing", The Benjamin/Cummings Publishing Company, 1994, pp.377〜406
【０００５】
【発明が解決しようとする課題】
上記の転置アルゴリズムによる並列計算方法では、Y方向、X方向、Z方向の各変換が各プロセッサで完全に独立に実行できる。しかし、途中の転置処理においては、全プロセッサ対全プロセッサのデータ転送が必要となる。一般に、分散メモリ型並列計算機では、データ転送は演算に比べて多くの時間を必要とし、この傾向は近年のプロセッサの高速化に伴ってますます強まりつつある。また最近では、分散メモリ型並列計算機の一種として、多数のパーソナルコンピュータ（PC）をイーサネット（登録商標）などのネットワークで結合したPCクラスタが広く使われるようになってきている。この場合は、専用の並列計算機に比べてプロセッサ間でのデータ転送能力が低いため、特にデータ転送処理が処理時間のネックとなりやすい。このような背景から、分散メモリ型並列計算機で３次元の高速フーリエ変換を行う場合、従来の転置アルゴリズムによる方式では十分な並列性能が得られない場合が多数生じてきている。本発明では、この問題を解決することを目的とする。
【０００６】
【課題を解決するための手段】
上記目的を達成するため、本発明では、転置アルゴリズムにおけるY方向とX方向の変換をそれぞれ２つの部分処理に分割して演算の順序を変更し、通信と演算のオーバーラップを可能にする。そのため、変換に関する次の性質に着目する。
▲１▼  Y方向の変換では、X方向に関して並列性がある。したがって、各プロセッサは、自分の行うべきNX組の変換を、任意の順序で行うことができる。
▲２▼  X方向の変換はlog2 NX個のステップからなるが、このうちで、最初のlog2NX-1個のステップでは、X座標が偶数の要素は偶数の要素同士、奇数の要素は奇数の要素同士でしか演算を行わない。
NX=8の場合を例に取り、X方向の変換における１組のデータに対する計算の流れを図４に示す。ここで、縦に並ぶ８個の丸7は各ステップでのデータの１要素を表し、丸と丸を結ぶ線8は、右側の丸の値を計算するために左側の丸の値を使うことを示す。図より明らかなように、変換のうち最初のlog2NX1=2個のステップでは、X座標が偶数の要素は偶数の要素同士、奇数の要素は奇数の要素同士でしか演算を行わない。したがって、上記で述べた性質▲２▼が成り立っていることがわかる。この性質はNX=8の場合に限らず、一般の場合について成り立つ。
上記▲１▼、▲２▼の２つの性質を利用すると、転置アルゴリズムにおけるY方向と方向の変換を次の５つの処理に分割して行うことができる。
＜処理１＞
X座標が偶数の要素のみに対し、Y方向の変換を行う。
＜処理２＞
X座標が奇数の要素のみに対し、Y方向の変換を行う。
X座標が偶数の要素に対しては、転置処理を行う。
＜処理３＞
X座標が偶数の要素のみに対し、X方向の変換の最初のlog2NX1ステップを行う。
X座標が奇数の要素に対しては、転置処理を行う。
＜処理４＞
X座標が奇数の要素のみに対し、X方向の変換の最初のlog2NX1ステップを行う。
＜処理５＞
X方向の変換の最後の１ステップを行う。
本発明による処理１〜処理４の様子を図５に示す。図において、斜線部分9および10は変換の計算を行っている要素、灰色の部分11はデータ転送を行っている要素、白の部分12は何も行わない要素を示す。また、直方体中の線はプロセッサへのデータ分割を示す。本発明の計算方式では、処理２においてX座標が奇数の要素に対して演算処理を行っている間に、X座標が偶数の要素に対して転置処理を行うことができる。また、処理３においては、X座標が偶数の要素に対して演算処理を行っている間に、X座標が奇数の要素に対して転置処理を行うことができる。これにより、転送のための時間を演算時間に隠蔽することが可能となり、計算効率を向上させることができる。
【０００７】
【発明の実施の形態】
≪実施例１≫
(1) 装置の概略構成
以下、本発明の原理および実施例を、図面により詳細に説明する。ここで実施例として挙げるのは、分散メモリ型並列計算機上で、３次元のフーリエ変換を用いて偏微分方程式を解き、これを用いてシミュレーションを行う場合である。ここでは特に、シミュレーションの例として気象予測のための計算を取り上げる。
本方法における並列プログラムを実行する並列計算機システムを図６に示す。システムは計算領域の形状、初期条件、物質定数などのパラメータを入力するための入力装置13、それぞれがメモリ14を備えたP台のプロセッサ15と、各プロセッサの持つメモリの間でデータ転送を行うためのネットワーク16から成る処理装置17、計算結果を出力するための出力装置18、プログラム19およびデータ20を格納するための外部記憶装置21から構成される。また、異なるプロセッサの持つメモリの間でのネットワークを介したデータ転送は、プロセッサの演算処理と同時に実行可能であるとする。
(2) 従来の方法による３次元並列高速フーリエ変換
以下、分散メモリ型並列計算機上での３次元高速フーリエ変換の原理を、数式を用いて説明する。３次元フーリエ変換は、Nx×Ny×Nz個の複素入力データ
{fjx,jy,jz}からNx×Ny×Nz個の複素出力データ{ckx,ky,kz}を、（数１）を用いて計算する処理である。
ckx,ky,kz=Σjx=0Nx-1Σjy=0Ny-1Σjz=0Nz-1
fjx,jy,jz exp(-2πikxjx/Nx) exp(-2πikyjy/Ny) exp(-2πikzjz/Nz)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；kz=0，1，… ，Nz-1）
--- （数１)
容易にわかるように，この式は次の(数２)(数３)(数４)でそれぞれ示すY方向の変換、X方向の変換、Z方向の変換を順に行うことによって計算することができる。
cjx,ky,jz(1) =Σjy=0Ny-1 fjx,jy,jz exp(-2πikyjy/Ny)
（ただしjx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
--- （数２)
ckx,ky,jz(2) =Σjx=0Nx-1 cjx,ky,jz(1) exp(-2πikxjx/Nx)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
--- （数２)
ckx,ky,kz =Σjz=0Nz-1 ckx,ky,jz(2) exp(-2πikzjz/Nz)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；kz=0，1，… ，Nz-1）
--- （数４）
これらの変換を分散メモリ型並列計算機上で行うには、（数２）がNX×NZ組の独立な変換、（数３）がNY×NZ組の独立な変換、（数４）がNX×NY組の独立な変換からなることを利用すればよい。この考え方に基づく計算方式が、「従来の技術」の項で説明した転置アルゴリズムである。転置アルゴリズムによる計算のフローチャートを図７に示す。この計算は、主に次の４つの処理からなる。
＜Y方向の変換＞
入力装置から３次元のデータfjx,jy,jzを入力し（図７の処理23）、直方体状に並べ（処理24）、Z軸に垂直な平面でデータをプロセッサ台数と等しい数に分割し、各部分データをそれぞれ１台のプロセッサのメモリに格納する（処理25）。データの分割の仕方には、ブロック分割、サイクリック分割など、いろいろな方式がありうるが（データ分割方式についての詳細はV. Kumar, A. Grama, A. Gupta and G. Karypis: "Introduction to Parallel Computing", The Benjamin/Cummings Publishing Company, 1994参照。）、たとえばブロック分割の場合、プロセッサp（0≦p≦P-1）はZ座標jzが(Nz/P)×p≦jz≦(Nz/P)×(p+1)-1をみたすNx×Ny×Nz/P個のデータを持つ。
この状態で、各プロセッサは自分のメモリに格納された部分データに対し、(数2)によるY方向の変換を行う（処理26）。
＜転置処理＞
次に、データを図３の6のようにY軸に垂直な面でプロセッサ台数と等しい数に分割し、各部分データをそれぞれ１台のプロセッサのメモリに格納する（処理27）。この場合もデータ分割の仕方はいろいろな方式がありうるが、たとえばブロック分割の場合、プロセッサp（0≦p≦P-1）はY座標kyが(Ny/P)×p≦ky≦(Ny/P)×(p+1)-1をみたすNx×Ny/P×Nz個のデータを持つ。このように分割方式を変更するため、プロセッサpはY方向の変換時に自分が持っていたデータのうち、Y座標kyが(Ny/P)×p'≦ky≦(Ny/P)×(p'+1)-1をみたすNx×Ny/P×Nz/P個のデータを、プロセッサp'に転送する。この処理を転置処理と呼ぶ。
＜X方向の変換＞
各プロセッサは自分のメモリに格納された部分データに対し、(3)式によるX方向の変換を行う（処理28）。
＜Z方向の変換＞
各プロセッサは自分のメモリに格納された部分データに対し、(4)式によるZ方向の変換を行う（処理29）。変換結果のデータを出力装置へ出力する（処理30）。
【０００８】
以上が従来の転置アルゴリズムによる３次元高速フーリエ変換である。しかしこの計算方式では、「従来の技術」の項で説明したように、途中の転置処理27において全プロセッサ対全プロセッサのデータ転送が必要となり、この部分が性能上のネックになる。
(3) 本発明の方法による３次元並列高速フーリエ変換
そこで本発明では、上記の転置アルゴリズムにおけるY方向とX方向の変換をそれぞれ２つの部分処理に分割して演算の順序を変更し、通信と演算のオーバーラップを可能にする。本発明による計算のフローチャートを図１に示す。ただし、本発明による計算は、図７に示した転置アルゴリズムによる計算と、処理25以前の部分、および処理29以後の部分は同じであるので、図１のフローチャートでは処理25より後で、かつ処理29より前の部分のみを示してある。この部分は、次の５つの処理からなる。
＜Y方向の変換処理(1)＞
第７図の分割ステップ２５で複数プロセッサのメモリに分散格納された部分データのうち、X座標jxが偶数の要素のみに対し、(2)式によるY方向の変換を行う（図１の処理33）。部分データの分割について再度述べると、変換対象データを直方体に配列した複素データ列についてＺ軸と垂直な面でプロセッサ台数と同じ数に分割する。したがって、上記複素データ列のＸ軸方向及びＹ軸方向は分割されていないので、プロセッサ間のデータ転送なしにＸ方向及びＹ方向での変換処理を完了できる分散配置となる。しかし、Y方向の変換処理(1)では、このうちX座標jxが偶数の要素に対するＹ方向の変換のみを各プロセッサにて行う。
＜Y方向の変換処理(2)＞
各プロセッサは、自分のメモリに格納された部分データのうち、X座標jxが奇数の要素のみに対し、(2)式によるY方向の変換を行う。
これと並行して、X座標jxが偶数の要素に対しては、転置処理を行う。転置処理においてプロセッサpは、従来の方法での転置処理と同様、自分が持っているデータのうち、Y座標kyが(Ny/P)×p'≦ky≦(Ny/P)×(p'+1)-1をみたすNx/2×Ny/P×Nz個のデータをプロセッサp'に転送する（図１の処理34）。つまりこの転地処理の結果、X座標jxのが偶数の要素についてのみ、プロセッサ間で再配置される。その結果、プロセッサ間のデータ転送なしに各プロセッサでＺ方向の変換が完了できる状態となる。またＸ方向については、隣接する要素同士の演算であるＸ方向の変換の最終ステップのみを除いたlog2NX1ステップを各プロセッサで実行できる状態となる。
＜X方向の変換処理(1)＞
各プロセッサは、自分のメモリに格納された部分データのうち、X座標jxが偶数の要素のみに対し、X方向の変換の最終ステップのみを除いた最初のlog2NX1ステップを行う。このときの計算式は、（数３）の代わりに
ckx,ky,jz(2') =Σjx'=0Nx/2-1 c2jx',ky,jz(1) exp(-2πikx2jx'/Nx)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
--- （数５）
とする。
これと並行して、X座標jxが奇数の要素に対しては、転置処理を行う。転置処理においてプロセッサpは、従来の方法での転置処理と同様、自分が持っているデータのうち、Y座標kyが(Ny/P)×p'≦ky≦(Ny/P)×(p'+1)-1をみたすNx/2×Ny/P×Nz個のデータをプロセッサp'に転送する（図１の処理35）。これにより、X座標jxのが奇数の要素についても、各プロセッサでＺ方向の変換が完了でき、またＸ方向については、隣接する要素同士の演算であるＸ方向の変換の最終ステップのみを除いたlog2NX1ステップを各プロセッサで実行できる状態に再配置される。
【０００９】
＜X方向の変換処理(2)＞
各プロセッサは、自分のメモリに格納された部分データのうち、X座標jxが奇数の要素のみに対し、X方向の変換の最終ステップのみを除いた最初のlog2NX1ステップを行う。このときの計算式は、(数３)の代わりに
ckx,ky,jz(2") =Σjx'=0Nx/2-1 c2jx'+1,ky,jz(1) exp(-2πikx(2jx'+1)/Nx)
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
--- （数６）
とする（図１の処理36）。
＜X方向の変換処理(3)＞
各プロセッサは、自分のメモリに格納されたデータを用いて、X方向の変換の最後の１ステップを行う。このときの計算式は、(数5)で計算したckx,ky,jz(2')、および(6)式で計算したckx,ky,jz(2")を用いて、
ckx,ky,jz(2) = ckx,ky,jz(2') + ckx,ky,jz(2")
（ただしkx=0，1，… ，Nx-1；ky=0，1，… ，Ny-1；jz=0，1，… ，Nz-1）
--- （数７）
とする（図１の処理37）。
以上の計算方式では、各プロセッサはのY方向の変換処理(2)においてX座標が奇数の要素に対してY方向の変換処理を行うと同時に、X座標が偶数の要素に対してデータの転送処理を行うことができる。また、各プロセッサはX方向の変換処理(1)において、X座標が偶数の要素に対してX方向の変換処理を行うと同時に、X座標が奇数の要素に対してデータの転送処理を行うことができる。これにより、転送のための時間の一部または全部を演算時間に隠蔽することが可能となり、計算効率を向上させることができる。
(4) 並列高速フーリエ変換ライブラリ
本発明を並列計算機上で３次元高速フーリエ変換を行うライブラリに適用した例を示す。本ライブラリは、サブルーチン名称をFFT3Dとし、実行するには
CALL FFT3D(NX, NY, NZ, P, F, TB, IOPT, IER)
のようにしてすべてのプロセッサで同時にコールする。ここで、NX，NY，NZはフーリエ変換を行うべき３次元データのそれぞれX、Y、Z方向の個数、Pはプロセッサ台数、Fは入力時はフーリエ変換すべき３次元データ{fjx,jy,jz}、出力時はフーリエ変換結果{ckx,ky,kz}を格納する配列、TBは変換に用いる複素指数関数の値を格納するテーブル、IOPTはサブルーチンの機能を指定する入力、IERは実行時エラーが生じたか否かを示す出力である。ここで、配列Fは各プロセッサがそれぞれ持つ部分配列であり、Z座標に対してブロック分割を行うので、第p番目のプロセッサp（0≦p≦P-1）はZ座標jzが(Nz/P)×p≦jz≦(Nz/P)×(p+1)-1をみたすNx×Ny×Nz/P個のデータのみを持つ。すなわち、第p番目のプロセッサの配列Fには、
F(jx, jy, jz') = fjx,jy,(Nz/P)*p+jz'
（ただしjx=0，1，… ，Nx-1；jy=0，1，… ，Ny-1；jz'=0，1，… ，Nz/P-1）
---（数８）
を格納する。したがって各プロセッサの持つ配列Fの大きさはNX×NY×NZ/Pである。また、TBは、第１回目のコールで計算した複素指数関数の値を格納しておくテーブルであり、２回目のコールからはここに格納した値を再利用することにより、新たな計算が不要となる。また、第１回目のコールではIOPT=1を指定し、このときは複素指数関数のテーブルを作成する。IOPT=2は２回目以降のコールを意味し、このときは既にTBに格納されている値を用いる。
【００１０】
本ライブラリのフローチャートを図８に示す。ライブラリは、コールされるとまず入力データをチェックし、NX，NY，NZとPとが１以上の整数であるかどうか、IOPTが１または２の値であるかどうかなど、入力の有効性を調べる（処理40）。入力データに無効な値が入っていた場合は、IER=1000と設定して（処理41）リターンする。次に、ライブラリが引数で指定した通りにP台のプロセッサでコールされているかどうかをチェックする（処理43）。この条件が満たされていない場合は、IER=2000と設定して（処理44）リターンする。次にIOPTの値をチェックし（処理45）、IOPT=1の場合は各方向のフーリエ変換で用いる複素指数関数の値を前もって計算し、配列TBに格納する（処理46）。
【００１１】
次に、Y方向およびX方向の変換を行う（処理47）。この変換は、本実施例における「(3) 本発明の方法による３次元並列高速フーリエ変換」の項で述べた方法に従い、Y方向の変換処理(1)（図１の処理33）からX方向の変換処理(3)（図１の処理37）の５つを順に行うことによって終了する。次に、本実施例における「(2) 従来の方法による３次元並列高速フーリエ変換」の項で述べた方法に従い、Z方向の変換を行う（図８の処理48）。これにより、３次元高速フーリエ変換の処理が終了する。終了時には、従来の方法による３次元並列高速フーリエ変換と同様、配列FにはY座標に対してブロック分割されたデータが格納される。すなわち、第p番目のプロセッサpの配列Fには、
F(kx, ky', kz) = ckx,(Ny/P)*p+ky',kz
（ただしkx=0，1，… ，Nx-1；ky'=0，1，…，Ny/P-1；kz=0，1，… ，Nz-1）
---（数９）
が格納される（図８の処理49）。
(5) シミュレーションプログラム
本実施例において実行すべき気象計算のための並列プログラムを図９に示す。ここでは、サイズNx×Ny×Nzの３次元メッシュで計算を行う場合を例にとって説明する。
【００１２】
本プログラムでは、まず計算領域のサイズNx，Ny，Nz、温度・風速・圧力などの初期条件、空気の熱伝導率をはじめとする物質定数などのパラメータを入力し（処理51）、計算に必要な前処理を行う。ここで前処理とは、観測によって得られた温度・風速・圧力などのデータに対して補間を行い、計算に必要なメッシュポイントでのデータを得ることである。次に、温度・風速・圧力などのデータを並列計算機の各プロセッサに分割する（処理52）。分割は、前項で述べた３次元高速フーリエ変換ライブラリFFT3Dを利用できるように、Z方向をブロック分割する形で行う。
【００１３】
これらの処理が終わった後、ループにより各時間ステップでの温度・風速・圧力などの量を順々に求めていく。基礎となる方程式は、風速に対する運動方程式
du/dt = -2Ω×u - (1/ρ)∇p + Fu ---（数１０)
質量保存の式
dρ/dt = -ρ∇・u                          ---（数１１)
温度変化を表す式
dT/dt = -κ∇2T + u・∇T            ---（数１２）
の３本である。ここで、uは風速、pは圧力、Tは温度を表し、Ωはコリオリ力と呼ばれる地球の自転による力、 Fuはそれ以外の外力、ρは空気の密度、κは空気の熱伝導率を表す。これらの式から次の時刻でのデータの値を求めるには、まずフーリエ変換により格子点上の温度T、圧力pおよび風速uを波数空間でのデータに変換し（処理53）、波数空間でこれらのデータを微分し（処理54）、再び波数空間でのデータを逆フーリエ変換して格子点上での温度勾配∇T 、２次微分∇2T 、圧力勾配∇p 、速度の発散∇・u を求める（処理55）。この後、これらの量を(8)-(10)式の右辺に代入し、次の時間ステップでの温度・風速・圧力を求める（処理56）。なお、フーリエ変換により格子点上のデータを波数空間上のデータに変換して処理を行うのは、その方が微分が精度良く計算できるからであり、本プログラムではこの部分で２次元フーリエ変換ライブラリFFT3Dを用いる。
上記のループでは、各時間ステップ毎に、求める時刻までの計算が終了したかどうかを判定し（処理57）、終了したら、後処理を行い（処理58）、結果を出力する（処理59）。後処理では主に、計算を行うメッシュポイントと結果データが必要な点とがずれている場合に、結果を補間して必要な点での値を計算するなどの処理を行う。
上記の例では気象予測計算を行う場合を例にとって説明したが、本発明の手法は、これ以外の応用例についても、並列計算機上で３次元高速フーリエ変換を用いてシミュレーションを行う場合に適用できることは明らかである。また、上記の例ではフーリエ変換前には３次元データをZ方向にブロック分割し、変換後にはデータはY方向にブロック分割されているとしたが、ブロック分割の代わりに、サイクリック分割やブロックサイクリック分割を用いた場合でも、本発明の手法が適用できることは明らかである。更に、上記の例ではY方向、X方向、Z方向の変換をこの順に行うとしたが、X、Y、Zなどは便宜的に定めた座標軸の名前に過ぎないため、上記の計算方法においてX→Y、Y→Z、Z→Xなどの名前の付け替えを行った計算方法も、本発明の手法と同じものであることは明らかである。
≪実施例２≫
上記実施例１では３次元フーリエ変換のための計算方法を説明したが、本方式は１次元フーリエ変換の場合へも応用できる。N点のデータに対して１次元フーリエ変換を行うには、 Nx×Ny×Nz=Nを満たす任意の整数Nx，Ny，Nzを定めてデータをNx×Ny×Nzの直方体状に並べ、これに対して「ひねり係数の乗算」という処理を加えた３次元フーリエ変換を行えば良いことが知られている。ここで、ひねり係数の乗算とは、(2)式によってY方向の変換を行った後、中間結果のcjx,ky,jz(1) に対して
cjx,ky,jz(1) := cjx,ky,jz(1)×exp(2πikyjx/(NxNy)
（ただしjx=0，1，… ，Nx-1，ky=0，1，… ，Ny-1，jz=0，1，… ，Nz-1）
---（数１３）
という処理を行い、かつ（数３）によってX方向の変換を行った後に中間結果のckx,ky,jz(2) に対して
＜X方向の変換＞
ckx,ky,jz(2) = ckx,ky,jz(2)×exp(2πi(Nykx+ky)jx/(NxNyNz))
---（数１４）
という処理を行うことである。このひねり係数の乗算は、配列cjx,ky,jz(1)またはckx,ky,jz(2)に対する要素毎に独立の処理であるので、実施例１の方式においても、Y方向の変換後とX方向の変換後にそれぞれ組み入れることができる。これにより、本発明の方式を１次元高速フーリエ変換に対しても適用することが可能となる。
≪実施例３≫
本発明は更に、２次元の高速フーリエ変換に対しても適用できる。この場合は、実施例１の方式による３次元高速フーリエ変換において、単にNz=1とすればよい。
【００１４】
補足すると、変換対象のデータはＸ軸とＹ軸の２次元に配列された複素データアレイである。実施例１で説明した「Y方向の変換処理(1)」の前のデータ分割の段階では、上記複素データアレイをＸ軸に垂直な直線でプロセッサ台数と等しい数の部分データに分割し、各プロセッサのメモリに分散配置することになる。つまり、プロセッサ間のデータ転送なしにＹ方向の変換が完了できる状態である。しかし「Y方向の変換処理(1)」では、X座標jxが偶数の要素に対するＹ方向の変換のみを各プロセッサにて行う。次に、Y方向の変換処理(2)でX座標jxが奇数の要素に対するＹ方向の変換を各プロセッサにて行うのと並行して、X座標jxが偶数の要素はプロセッサ間で再配分を行う。この再配分とは、X座標jxが偶数の要素をＹ軸と垂直な直線でプロセッサ台数と等しい数の部分データに分割して各プロセッサに分散配置する処理であり、この結果、各プロセッサではX方向の変換の最終ステップを除いたlog2NX1ステップが実行可能となる。さらにX方向の変換処理(1)では、再配分されたX座標jxが偶数の要素に対するX方向の変換（log2NX1ステップ）を各プロセッサにて行うのと並行して、X座標jxが奇数の要素のプロセッサ間での再配分を行う。この再配分も上記と同様である。更に、X方向の変換処理(2)、 X方向の変換処理(3)と進むのは実施例１と全く同様である。変換対象のデータが２次元データであるので、X方向の変換処理(3)で処理が完了する。
≪実施例４≫
本発明による高速フーリエ変換を用いてシミュレーションを行う他の例として、半導体デバイス等における電子構造計算を説明する。電子構造計算では、３次元のメッシュで定義された電子の波動関数u(r)を、シュレディンガー方程式
du(r)/dt = -(h2/2m)∇2u(r) + (E-V(r))u(r)---（数１５）
に従って計算することにより、半導体の性質を決定するバンドギャップの大きさや、結晶の構造安定性などを求める。ただし、上式でhはプランク定数、mは電子の質量、Eは対象とする波動関数のエネルギーレベル、Vは結晶中の原子や他の電子によるポテンシャルエネルギーを表す。
【００１５】
(15)式の計算では、波動関数の２次微分∇2u(r)が必要であるが、気象計算の例において述べたのと同様な理由により、この部分はu(r)をフーリエ変換により波数空間に移してから計算し、結果を逆フーリエ変換で再び実空間に戻す。したがって、並列計算機上で電子構造計算を行う場合には、この部分で本発明の３次元高速フーリエ変換方法が適用できる。
【００１６】
【発明の効果】
以上説明したように、本発明によれば、分散メモリ型並列計算機上での高速フーリエ変換において、演算処理とデータ転送処理を並行して行い、後者の時間の一部または全部を前者の時間に隠蔽することが可能となる。したがって本発明の方法では、従来法に比較して並列化効率を向上させることができる。なお、この向上効果は、本発明を適用する分散メモリ型並列計算機のプロセッサ間通信の性能に依存するが、たとえば16台のプロセッサでN=2563の問題を処理する場合、20%〜30%程度の実行時間短縮が期待できる。
【図面の簡単な説明】
【図１】本発明による並列３次元高速フーリエ変換方法におけるY方向及びX方向の変換の計算方法を示すフローチャート。
【図２】従来法による３次元高速フーリエ変換方法を示す図。
【図３】従来法による並列３次元高速フーリエ変換方法を示す図。
【図４】  X方向の変換におけるデータの依存関係を示す図。
【図５】本発明による並列３次元高速フーリエ変換方法におけるY方向及びX方向の変換の計算方法を示す図。
【図６】本発明を適用すべき分散メモリ型並列計算機の構成を示す図。
【図７】従来法による並列３次元高速フーリエ変換方法を示すフローチャート。
【図８】本発明による並列３次元高速フーリエ変換ライブラリのフローチャート。
【図９】本発明を利用して分散メモリ型並列計算機で気象計算を行う場合のフローチャート。
【符号の説明】
1：NX×NY×NZの直方体状に並べたデータ、2：Y軸に平行なNY個のデータ、3：X軸に平行なNX個のデータ、4：Z軸に平行なNZ個のデータ、5：Z軸に垂直な平面で分割した部分データ、6：Y軸に垂直な平面で分割した部分データ、7：X方向の変換における各ステップでのデータの１要素、8：右側の要素の値を計算するために左側の要素の値を使うことを示す線、9：Y方向の変換を行っている要素、10：X方向の変換を行っている要素、11：データ転送を行っている要素、12：何も行わない要素、13：入力装置、14：メモリ、15：プロセッサ、16：ネットワーク、17：処理装置、18：出力装置、19：プログラム、20：データ、21：外部記憶装置、22：スタート、23：変換すべきデータの入力、24：入力データを直方体状に並べる、25：データをZ軸に垂直な面でプロセッサ間に分割、26：Y方向の変換、27：データをY軸に垂直な面で再分割、28：X方向の変換、29：Z方向の変換、30：変換結果の出力、31：終了、32：図７の処理25より継続、33：Y方向の変換処理(1)、34：Y方向の変換処理(2)、35：X方向の変換処理(1)、36：X方向の変換処理(2)、37：X方向の変換処理(3)、38：図７の処理29へ継続、39：ライブラリコール、40：入力データの有効性チェック、41：エラーコードの設定、42：他のプロセッサにコールを通知、43：P台のプロセッサでコールしていることをチェック、44：エラーコードの設定、45：IOPTの値のチェック、46：複素指数関数テーブルの設定、47：Y方向及びX方向の変換、48：Z方向の変換、49：変換結果の格納、50：リターン、51：計算領域の形状・初期条件・パラメータの入力、52：前処理とデータの分割、53：温度及び風速データを端数空間でのデータに変換、54：端数空間で温度勾配及び圧力データを計算、55：温度勾配及び圧力データを格子点上のデータに変換、56：次の時間ステップでの温度と風速データの計算、57：判定処理、58：物理量の計算、59：結果の出力。
【図１】
0004052181.tif 000002
【図２】
0004052181.tif 000003
【図３】
0004052181.tif 000004
【図４】
0004052181.tif 000005
【図５】
0004052181.tif 000006
【図６】
0004052181.tif 000007
【図７】
0004052181.tif 000008
【図８】
0004052181.tif 000009
【図９】
0004052181.tif 000010