JP2000200261A - フ―リエ変換方法、シミュレ―ション方法およびプログラム記録媒体 - Google Patents

Abstract

(57)【要約】 【課題】並列計算機上でフーリエ変換を、変換対象デー
タと変換結果データがデータ分割形式が同じに保ち、か
つ、高速に実行する。 【解決手段】変換対象一次元データを直方体状に並べて
３次元の変換対象データに写像し、この直方体をＺ方向
に垂直な面で分割して各面のデータを一つのプロセッサ
に割り当て、各プロセッサで、この割り当てにしたがっ
て３次元の変換対象データに対してＹ方向のフーリエ変
換を実行して第１変換結果データを得る(56)。各プロセ
ッサで、第１変換結果データに対してＸ方向のフーリエ
変換に類似の変換を実行して第２変換結果データを得る
(57)。第２変換結果データをＹ軸に垂直な面で分割し直
し、その分割に従いそのデータをプロセッサ間で並び替
え(58)、その後に、各プロセッサで、Ｚ方向のフーリエ
変換に類似の変換を行い最終的な三次元フーリエ変換結
果データを得る(59)。

Description

【発明の詳細な説明】

【発明の属する技術分野】本発明は、複数のプロセッサ
を有する計算機で実行するのに適したフーリエ変換方法
に係り、とくに、ベクトル演算器を内蔵する複数のプロ
セッサからなるベクトル並列計算機で実行するのに適し
たフーリエ変換方法に関する。

【従来の技術】科学技術計算において頻繁に利用される
処理の一つに、フーリエ変換がある。フーリエ変換は、
物理現象のシミュレーションその他に使用される。フー
リエ変換は、ある実数区間で定義された複素数値をとる
関数ｆ（ｘ）を複素指数関数ｅｘｐ（ｉｋｘ）の重ね合
わせとして表す処理であり、計算機上で実現する場合に
は、扱いうる点の数が有限であることから、複素数の点
列ｆ₀，ｆ₁，．．．，ｆ_N-1をＮ個の複素指数関数ｅｘ
ｐ（２πｉｋｊ／Ｎ）（ただし、ｋ＝０，１，．．．，
Ｎ−１で、ｉは虚数単位、πは円周率）の重ね合わせと
して表す処理となる。すなわち、ｆ₀，ｆ₁，．．．，ｆ
_N-1が与えられたときに、式１ａにより重ね合わせの係
数ｃ₀，ｃ₁，．．．，ｃ_N-1を求めるのがフーリエ変換
である。各点ｆ_jの値は、これらの係数を用いると、式
１ｂによりあらわされる。

【数１】 c_k＝(1/N)Σ_j=0 ^N-1 f_j exp(-2πikj/N) （ただし、k=0,1,...,N-1） （１ａ） f_j＝Σ_k=0 ^N-1 c_k exp(2πikj/N) （ただしj=0,1,...,N-1） （１ｂ） しかし、この定義に基づいて計算を行うと、式の数がＮ
本あり、各式がＮ個の項から成るため、複素指数関数ｅ
ｘｐ（−２πｉｋｊ／Ｎ）の計算に加えて、複素数の加
算と乗算が約Ｎ²回必要である。そこで実際には、アル
ゴリズム上の工夫により計算量を約ＮｌｏｇＮのオーダ
ーに減少させた高速フーリエ変換という手法が広く使わ
れている。高速フーリエ変換を並列計算機上で効率的に
行うための手法として、従来、転置アルゴリズムとバイ
ナリ・エクスチェンジと呼ばれる２つの手法が提案され
ている（たとえばV. Kumar, A. Grama, A. Gupta and
G. Karypis: "Introductionto Parallel Computing, Th
e Benjamin/Cummings Publishing Company, 1994参
照）。前者はプロセッサ間の通信を計算途中の一箇所に
まとめて行う方式、後者はプロセッサ間で通信を行いな
がら計算を進める方式であり、プロセッサの台数をｐと
すると、通信の回数は前者がｐ−１、後者がｌｏｇ₂ｐ
で、通信１回あたりに送るデータ量は、前者がＮ／
ｐ²、後者がＮ／２ｐである。後者は前者に比べて通信
の回数が少なくて済むため、通信のセットアップ時間が
支配的となる小規模問題では通信時間が少なくて済むと
いう利点があるが、通信すべきデータの総量は多くなる
ため、大規模データの場合には前者が有利となる。半導
体デバイスの特性計算、電子状態計算、気象予測のため
の計算などの科学技術計算では、数万から数百万に上る
変数を扱う大規模シミュレーションが必要である。この
ような大規模問題を扱う手段としては、並列計算機が有
力である。並列計算機は数十個から数万個に上る多数の
高速プロセッサをネットワークで結んだシステムであ
り、従来の逐次型計算機に比べ、プロセッサ台数を増や
すことでピーク性能をいくらでも高めることができると
いう利点を持つ。さらに、最近の並列計算機では、各プ
ロセッサで、一連のデータに対して同じ演算を高速に実
行するできるように演算器が構成されていることも多
い。とくに、各プロセッサに、そのような演算器として
同じ演算を複数のデータに対してパイプライン的に実行
するベクトル演算器を有するベクトル並列計算機も開発
されている。ベクトル並列計算機の中には、このベクト
ル演算器による演算を指定するベクトル命令を実行でき
るものもある。さらに、メモリとベクトル演算器の間に
複数のベクトルレジスタが設けられている並列計算機も
ある。これらのベクトルレジスタはメモリと演算器のデ
ータの転送時間が処理時間に及ぼす時間を軽減してい
る。より高速にシミュレーションを実行可能になってい
る。また、厳密にはベクトル並列計算機ではないが、ベ
クトル並列計算機に類似の並列計算機として、ある演算
を実行する演算器がベクトル演算器でなくても、一連の
データに対してその演算を高速に実行できるように構成
されている演算器を使用する並列計算機も多い。フーリ
エ変換は科学技術計算でもっともよく使われる処理の一
つであり、最近では並列計算機用のライブラリとして提
供されることも多い。たとえば日立製作所編「プログラ
ムプロダクトＨＩ−ＵＸ／ＭＰＰ行列計算副プログラム
ライブラリＭＡＴＲＩＸ／ＭＰＰ」参照。並列計算機で
実行する大規模のシミュレーションがフーリエ変換を実
行する場合には前述の転置アルゴリズムが使用されるこ
とが多い。上に記載したように、ベクトル型並列計算機
あるいはそれらに類似の並列計算機で転置アルゴリズム
を実行する場合には、変換すべき一次元空間の点列デー
タを３次元空間に直方体状に並べ、これに対してたとえ
ばＹ方向の変換、Ｘ方向の変換、Ｚ方向の変換を順次行
うことによって、全データ点列に対して高速フーリエ変
換を行ったのと同一の結果を得る。より具体的には、フ
ーリエ変換の対象となる一次元のデータｆ₀，
ｆ₁，．．．，ｆ_N-1を入力し、各辺の長さがＮＸ，Ｎ
Ｙ，ＮＺの直方体状に並べる。ここで、ＮＸ，ＮＹ，Ｎ
ＺはＮＸ＊ＮＹ＊ＮＺ＝Ｎを満たす整数である。データ
を直方体状に並べるに当たっては、原点からたとえばＺ
方向にデータを並べていき、ＮＺ個のデータを並べ終わ
ったら次はＸ座標を１だけ増やしてデータを並べ、これ
を繰り返してＮＸ＊ＮＺ個のデータを並べ終わったら次
はＹ座標を１だけ増やしてデータを並べる、という操作
を行う。このようにデータを並べた後、直方体をＺ軸に
垂直にスライスし、こうしてできる各面を並列計算機の
一つのプロセッサに割り当てる。次に、Ｙ方向の変換を
行う。プロセッサへの入力データｆ_jの割り当て方式よ
り、各ＸＹ平面は１台のプロセッサに担当されているか
ら、この変換処理は通信なしに各プロセッサで独立に行
える。次に、同様にして各プロセッサで独立にＹ方向の
変換の結果データに対してＸ方向の変換を行う。Ｘ方向
の変換の終了後、プロセッサ間でのＸ方向の変換の結果
データの入れ替えを行い、今度はその結果データが構成
する直方体をＸ軸に垂直にスライスし、こうしてできる
各面を一つのプロセッサに割り当てる。この処理を転置
と呼び、各プロセッサが自分以外の全プロセッサとデー
タの交換を行う必要がある。転置の終了後、今度は各プ
ロセッサで独立にＺ方向の変換を行う。以上で、直方体
状に並べられた一次元入力データｆ_jのフーリエ変換が
終了し、直方体状に並べられた、重ね合わせの係数を表
す出力データｃ_kが求まる。出力データｃ_kの並び方は、
原点からまずＹ方向に、Ｙ方向にＮＹ個行ったら次はＸ
座標が１だけ増え、ＸＹ平面上にＮＸ＊ＮＹ個のデータ
が並んだら次はＺ座標が１だけ増えるという順で並ぶ。
上記の転置アルゴリズムでは、入力データｆ_jの分割
は、ｆ_jを第ＭＯＤ（j，ｐ）番のプロセッサが担当する
という形でデータがプロセッサ間で分割されている。こ
のデータ分割形式はサイクリック分割と呼ばれる。デー
タ分割形式はデータのプロセッサへの割り当ての順序を
表すものでもあり、本明細書ではデータ分割形式のこと
を割り当て順序あるいは割り当て態様とも呼ぶことがあ
る。一方、出力データｃ_kの分割は、ＮＹ個の連続する
データを１台のプロセッサが担当するブロックサイクリ
ック分割となり、入力データｆ_jとはプロセッサ間のデ
ータ分割形式が異なる。上記の転置アルゴリズムでは、
入力データｆ_jの並べ方およびそのデータのプロセッサ
への分割の仕方より分かるように、入力データの分割形
式は、ｆ_jを第ＭＯＤ（j，ｐ）番のプロセッサが担当す
るというサイクリック分割となる。一方、転置後のデー
タのプロセッサへの分割の仕方、および変換で得られた
出力データｃ_kの並び方より分かるように、出力データ
の分割形式は、ＮＹ個の連続するデータを１台のプロセ
ッサが担当するというブロックサイクリック分割とな
る。しかし多くの応用では、フーリエ変換と逆フーリエ
変換とを対にして用い、しかも逆フーリエ変換はフーリ
エ変換プログラムを流用して行うため、フーリエ変換の
入力データと出力データが同じデータ分割形式（データ
割り当て順序）になっている方が都合がよい。そのた
め、従来の高速フーリエ変換方法では、以上の処理に従
ってブロックサイクリック分割の出力データｃ_kを得た
後、再びプロセッサ間でデータの転送を行い、データｃ
_kをサイクリック分割に直して出力する必要がある。

【発明が解決しようとする課題】本発明者の検討の結
果、以上の従来のフーリエ変換方法では、フーリエ変換
係数の計算後に行うデータ分割形式（データ割り当て順
序）の変更のためのプロセッサ間でのデータ転送が、フ
ーリエ変換時間の短縮を妨げていることが分かった。し
たがって、本発明の目的は、フーリエ変換係数の計算後
にデータ分割形式の変更のためにデータ転送を行わなく
ても、フーリエ変換結果データがフーリエ変換対象デー
タと同一のデータ分割形式（データ割り当て順序）を持
ち得るフーリエ変換方法を提供することである。

【課題を解決するための手段】上記目的を達成するた
め、本発明によるフーリエ変換方法は、各プロセッサに
より、第１の変換処理、第２の変換処理、第３の変換処
理を順次実行し、上記複数のプロセッサの各々による、
上記第１、第２の変換処理のいずれか一方の変換処理の
実行後に、上記複数のプロセッサでのその一方の変換処
理を実行した結果得られた一群の結果データを構成する
複数の結果データ部分群が異なるプロセッサに割り当て
られるように、上記一群の結果データを上記複数のプロ
セッサの間で交換するステップを有する。上記第１から
第３の変換処理は、一群の順序づけられた変換対象デー
タに対する一群の順序づけられたフーリエ変換係数デー
タを構成する複数のフーリエ変換係数データ部分群をそ
れぞれ異なるプロセッサにより生成するように定めら
れ、各プロセッサには、上記一群の変換対象データを構
成する複数の変換対象データ部分群の一つの変換対象デ
ータ部分群がそのプロセッサに対して予め割り当てら
れ、上記一群のフーリエ変換係数データのそれぞれを生
成したプロセッサの順序が、上記一群の変換対象データ
のそれぞれが割り当てられたプロセッサの順序と同一と
なるように、上記交換するステップで各プロセッサに割
り当てられる結果データ部分群が定められているもの。
より具体的には、上記第１、第２、第３の変換処理は、
それぞれ３次元データ空間の第１、第２、第３の座標軸
に関する変換処理であり、上記変換対象データ群の各々
は、上記３次元データ空間の直方体形状に位置する格子
点群の一つの座標をそれぞれ有し、上記複数の変換対象
データ部分群は、上記３次元データ空間の第３の座標軸
の座標値が同じであり、上記３次元データ空間の第１、
第２の座標軸の座標値が異なる全ての変換対象データが
同一の変換対象データ部分群に含まれるように定めら
れ、上記フーリエ変換係数データ群の各々は、３次元係
数空間の直方体形状に位置する格子点群の一つの座標を
それぞれ有し、上記複数のフーリエ変換係数データ部分
群は、上記３次元係数空間の第１の座標軸の座標値が同
じであり、上記３次元波数空間の第２、第３の座標軸の
座標値が異なる全てのフーリエ変換係数データが同一の
フーリエ変換係数データ部分群に含まれるように定めら
れている。本発明の具体的な態様によるフーリエ変換方
法では、上記変換対象データ群の各々は、上記３次元デ
ータ空間の直方体形状に位置する格子点群の一つの座標
をそれぞれ有し、上記複数の変換対象データ部分群は、
上記３次元データ空間の第３の座標軸の座標値が同じで
あり、上記３次元データ空間の第１、第２の座標軸の座
標値が異なる全ての変換対象データが同一の変換対象デ
ータ部分群に含まれるように定められ、上記フーリエ変
換係数データ群の各々は、３次元係数空間の直方体形状
に位置する格子点群の一つの座標をそれぞれ有し、上記
複数のフーリエ変換係数データ部分群は、上記３次元係
数空間の第１の座標軸の座標値が同じであり、上記３次
元波数空間の第２、第３の座標軸の座標値が異なる全て
のフーリエ変換係数データが同一のフーリエ変換係数デ
ータ部分群に含まれるように定められる。更に具体的に
は、上記変換対象データ群が上記３次元データ空間に直
方体形状に位置する格子点群に上記３次元空間に第３の
座標軸、第２の座標軸、第１の座標軸の順に順次割り当
てられ、上記第１から第３の変換処理は、上記複数のフ
ーリエ変換係数データが、３次元係数空間に直方体形状
に位置する格子点群に、当該３次元係数空間の第１、第
２、第３の座標軸の順序で割り当てられるように定めら
れている。更に具体的な態様では、各プロセッサが上記
第１の変換処理により生成する上記一つの第１の結果デ
ータ部分群は、上記３次元データ空間の第３の座標軸の
座標値が所定の同じ値であり、上記３次元データ空間の
第２の座標軸の座標値と上記３次元係数空間の第１の座
標軸の座標値が異なる値を有する全ての複数の第１の結
果データを含み、上記交換ステップが上記第１の変換処
理が上記複数のプロセッサにより実行された後に実行さ
れ、上記複数のプロセッサは、この交換ステップで、上
記３次元係数空間の第１の座標軸の座標値が所定の同じ
値であり、上記３次元データ空間の第２、第３の座標軸
の座標値が異なる値を有する全ての複数の第１の結果デ
ータを含む第１の結果データ部分群が同一のプロセッサ
に割り当てられるように、上記複数のプロセッサが生成
した一群の第１の結果データを上記複数のプロセッサの
間で交換し、各プロセッサが上記第２の変換処理により
生成する上記一つの第２の結果データ部分群は、上記３
次元係数空間の第１の座標軸の座標値が所定の同じ値で
あり、上記３次元波数空間の第２の座標軸の座標値と上
記３次元データ空間の第３の座標軸の座標値が異なる値
を有する全ての複数の第２の結果データを含み、各プロ
セッサが上記第３の変換処理により生成する上記一つの
フーリエ変換係数部分群は、上記３次元係数空間の第１
の座標軸の座標値が所定の値であり、上記３次元波数空
間の第２、第３の座標軸の座標値が異なる値を有する全
ての複数のフーリエ変換係数を含む。更に具体的な他の
態様では、各プロセッサが上記第１の変換処理により生
成する上記一つの第１の結果データ部分群は、上記３次
元データ空間の第３の座標軸の座標値が所定の同じ値で
あり、上記３次元データ空間の第２の座標軸の座標値と
上記３次元係数空間の第１の座標軸の座標値が異なる値
を有する全ての複数の第１の結果データを含み、上記交
換ステップが上記第２の変換処理が上記複数のプロセッ
サにより実行された後に実行され、各プロセッサが上記
第２の変換処理により生成する上記一つの第２の結果デ
ータ部分群は、上記３次元データ空間の第３の座標軸の
座標値が所定の同じ値であり、上記３次元係数空間の第
１、第２の座標軸の座標値が異なる値を有する全ての複
数の第２の結果データを含み、上記複数のプロセッサ
は、上記交換ステップにより、上記３次元係数空間の第
１の座標軸の座標値が所定の同じ値であり、上記３次元
係数空間の第１の座標軸の座標値と上記３次元データ空
間の第３の座標軸の座標値が異なる値を有する全ての複
数の第１の結果データを含む第１の結果データ部分群が
同一のプロセッサに割り当てられるように、上記複数の
プロセッサが生成した一群の第１の結果データを上記複
数のプロセッサの間で交換し、各プロセッサが上記第３
の変換処理により生成する上記一つのフーリエ変換係数
部分群は、上記３次元係数空間の第１の座標軸の座標値
が所定の値であり、上記３次元係数空間の第２、第３の
座標軸の座標値が異なる値を有する全ての複数のフーリ
エ変換係数を含む。本発明のより具体的な態様では、各
プロセッサにより、３次元空間の第１、第２、第３の座
標軸の座標にそれぞれ関する第１、第２、第３の変換処
理を順次かつ他のプロセッサと並行して実行し、各プロ
セッサが上記第１、第２の変換処理のいずれか一方を実
行した後に、その一方の変換処理の結果それぞれのプロ
セッサで得られた複数の結果データを上記複数のプロセ
ッサの間で交換するステップを有する。ここで、一群の
順序づけられた変換対象データが上記３次元空間に直方
体の形に並べられ、上記第１から第３の変換処理は、上
記一群の変換対象データに対する一群の順序づけられた
３次元空間の座標を有する複数のフーリエ変換係数デー
タを生成するように定められ、上記複数の変換対象デー
タが構成する上記直方体を分割する上記３次元空間の上
記第１の座標軸に垂直な複数の面の各々に含まれる複数
の変換対象データが同一のプロセッサに割り当てられ、
上記交換ステップは、上記一方の変換処理の結果得られ
た上記複数の結果データが構成する３次元空間の直方体
を、その３次元空間の第１の座標軸に垂直な複数の面で
分割し直して、各面に属する複数の結果データを同一の
プロセッサに割り当てるように、上記一方の変換処理の
結果得られた上記複数の結果データを上記複数のプロセ
ッサ間で交換するステップを有する。とくに、望ましく
は、上記一群の順序づけられた変換対象データを上記３
次元空間に直方体の形に並べられる順序は、第３の座標
軸、第２の座標軸、第１の座標軸の順であり、上記第１
から第３の変換処理は、上記複数のフーリエ変換係数デ
ータが３次元空間に第１、第２、第３の座標軸の順序で
並べられるように定められている。さらに望ましくは、
本発明によるフーリエ変換方法は、各プロセッサがパイ
プライン演算器を含み、その演算器での演算の対象とす
るループ長がＬのときのその各プロセッサの演算性能を
求めるための性能データを上記複数のプロセッサに共通
に記憶し、その性能データを用いて、上記直方体の上記
第１、第２、第３の座標軸方向の長さを決定し、その決
定された上記第１、第２、第３の座標軸方向の長さを有
する直方体に、上記順序づけられた複数の変換対象デー
タを並べるステップをさらに有する。本発明によるプロ
グラム記憶媒体は、上記いろいろのフーリエ変換方法の
いずれかを実行するようにプログラムされたプログラム
が記憶する。さらに、本発明によるシミュレーション方
法は、上記いろいろのフーリエ変換方法のいずれかを使
用してシミュレーションを実行する。本発明による他の
プログラム記憶媒体は、上記シミュレーション方法を実
行するようにプログラムされたプログラムを記憶する。

【発明の実施の形態】以下、本発明に係るフーリエ変換
方法、それを用いるシミュレーション方法およびプログ
ラムを図面に示したいくつかの実施の形態を参照してさ
らに詳細に説明する。なお、以下においては、同じ参照
番号は同じものもしくは類似のものを表わすものとす
る。また、第２の実施の形態以降では、第１の実施の形
態との相違点を主に説明するに止める。 ＜発明の実施の形態１＞ （１）装置の概略構成 本発明によるフーリエ変換方法を実行するための並列計
算機システムの一例を図１に示す。並列計算機２８は、
それぞれがメモリ２６を備えた複数のプロセッサ２７
と、プログラムおよびデータを格納するための複数の外
部記憶装置３１から構成され、これらの装置は、内部デ
ータ転送ネットワーク２９を介して相互にデータを交換
可能なように構成されている。外部記憶装置３１には、
たとえば、多くのユーザの利用に供するために並列計算
機２８に予め準備された複数のプログラムライブラリ４
４とそれらが使用するデータ３０等が記憶される。各プ
ロセッサのメモリ２６は、いわゆるローカルメモリであ
り、このメモリに記憶されたデータに割り当てられるア
ドレスは、そのプロセッサで定められたローカルなアド
レス空間に属するアドレスであり、この種のメモリは一
般に分散メモリと呼ばれ、この種のメモリを有する計算
機は分散メモリ型の並列計算機と呼ばれる。並列計算機
２８は、各プロセッサ２９が、一連のデータ要素からな
るベクトルデータに対して同じ演算をパイプライン的に
連続して実行できるベクトル演算器（図示せず）を備え
るベクトル並列計算機であると仮定する。これらのプロ
セッサ２７内の特定の一つのプロセッサには、ユーザが
操作可能な計算機、たとえばワークステーション１がＬ
ＡＮ等のネットワーク２を介して接続されている。この
計算機は他の計算機たとえばパーソナルコンピュータで
もよい。このワークステーション１には、並列計算機２
８に対する指示あるいはデータを入力するための入力装
置３（典型的には、キーボードとマウス）と、並列計算
機２８からの計算結果を出力するための出力装置２９
（典型的には、表示装置と印刷装置）が接続されてい
る。なお、ワークステーション１内には、並列計算機２
８に送るべきプログラムおよびそのプログラムで使用す
るデータを記憶する記憶装置（図示せず）も設けられて
いる。上記特定のプロセッサは、並列計算機２８内で計
算を司るプロセッサの役目とユーザ用のワークステーシ
ョン１との通信の役目とを兼ねる。すなわち、このプロ
セッサは、ワークステーション１から送付されるプログ
ラムとデータを受信し、それらを外部記憶装置３１の一
つに記憶し、その後、並列計算機２８の内部に記憶され
た適当なプログラムにより、ユーザ指定のプログラムを
複数のプロセッサ（具体的には全プロセッサ）にロード
し、ユーザ指定のデータの異なる部分を、それぞれそれ
らのプロセッサの異なるものに割り当て、そのユーザ指
定のプログラムを起動する。しかしながら、本発明によ
るフーリエ変換方法を実施するためには、並列計算機２
８は、複数のプロセッサを有することが必要であるが、
それ以外の点では特に限定した構造を有しなくてもよい
ことは言うまでもない。並列計算機２８は、ベクトル並
列計算機であると仮定したが、このベクトル演算器はご
く一部の演算のみを実行でき、他の演算はベクトル演算
器ではないスカラ演算器で実行されてもよい。さらに、
並列計算機２８は、対してこのような演算器を有しなく
てもよい。もちろん、一連のデータに対する同じ演算を
高速に実行できるように構成されている演算器を有する
ことが望ましい。また、並列計算機２８は、メモリ２９
と演算器（図示せず）の間に複数のベクトルレジスタを
有しないと仮定するが、これらのレジスタが使用するこ
とはより望ましいことである。さらに、それらのプロセ
ッサの具体的な構造あるいはそれらの間のデータ転送ネ
ットワークの構造、あるいはそれらのプロセッサと入力
装置あるいは出力装置との接続形態がいろいろであって
も、本発明はそれらの並列計算機に適用可能である。た
とえば、ワークステーションと通信可能な複数のプロセ
ッサが設けられていてもよく、また、ワークステーショ
ンと通信可能な少なくとも一つのプロセッサが計算用の
プロセッサとは別に設けられていてもよい。また、実行
すべきプログラムとデータを並列計算機２８に送付する
方法は他の方法に依ってもよいことは明らかである。ユ
ーザは上記複数のプロセッサ２７を使用して種々の計算
を実行できる。最も典型的な計算は、物理現象などのシ
ミュレーションであり、たとえば、地球の気象の予測も
シミュレーションにより行われる。半導体デバイスの設
計も、半導体デバイスの物理的な動作をシミュレーショ
ンして行われている。このようなシミュレーションを並
列計算機を使用して実行する場合、シミュレーション対
象の物理領域を複数の部分領域に区分し、各部分領域を
一つのプロセッサに割り当て、そのプロセッサにおい
て、その部分領域についてのシミュレーションを、たと
えば一つまたは複数の物理量に関する偏微分方程式を解
いて実行することが多い。この場合、シミュレーション
に使用される複数のプロセッサは、同じシミュレーショ
ンプログラムを互いに並列に実行する。したがって、こ
のようなプログラムは並列プログラムとも呼ばれる。各
プロセッサが実行するシミュレーションプログラムが使
用するデータは異なる。たとえばシミュレーション領域
の位置と形状を表すデータ、シミュレーションすべき物
理量の初期値、シミュレーション領域の物質に関する物
質定数、あるいは各部分領域に関する境界条件など異な
る。各プロセッサは、計算の途中で得られた結果データ
を他の適当なプロセッサに転送し、あるいは他のプロセ
ッサから計算結果データを受け取り、さらにシミュレー
ションを続けていく。このシミュレーションプログラム
の中にはフーリエ変換を使用するものもある。本実施の
形態では、いろいろのシミュレーションの利用に供する
ために、本発明によるフーリエ変換方法にしたがってフ
ーリエ変換を実行するようにプログラムされたフーリエ
変換ライブラリがいずれかの外部記憶装置３１に記憶さ
れる。さらにプロセッサ間の通信を実行するための通信
ライブラリも外部記憶装置３１に記憶される。シミュレ
ーションプログラムは、上記フーリエ変換ライブラリあ
るいは通信ライブラリを必要な時点でコールするように
プログラムされる。並列計算機２８は、ワークステーシ
ョン１から送信されたユーザ指定のシミュレーションプ
ログラムと、そのシミュレーションプログラムが使用す
るライブラリ（今の場合には上記フーリエ変換ライブラ
リと上記通信ライブラリ）を各プロセッサにロードす
る。さらに、並列計算機２８は、それぞれのプロセッサ
でシミュレーションプログラムが使用する、ワークステ
ーション１から送信されたユーザ指定のデータをそれぞ
れのプロセッサにロードする。なお、シミュレーション
プログラムは全プロセッサにロードされてもよく、一部
のプロセッサにロードされてもよいが、以下では簡単化
のために、シミュレーションプログラムは、全てのプロ
セッサにロードされると仮定する。上記ライブラリある
いは上記シミュレーションプログラムは、並列計算機２
８の命令セットあるいはハード構造の特徴、ソフトウエ
ア上の制約等を反映するコンパイラによりコンパイルさ
れたものである。本発明によりフーリエ変換方法を実行
する上記ライブラリあるいは上記シミュレーションプロ
グラムに上記ライブラリに含まれたフーリエ変換のため
のプログラム部分を組み込んだプログラムを磁気記憶装
置のようなプログラム記録媒体に記憶して販売できる。 （２）並列高速フーリエ変換の原理 すでに述べたごとく、フーリエ変換は、Ｎ個の入力デー
タｆ₀，ｆ₁，．．．，ｆ_N-1からＮ個の出力データｃ₀，
ｃ₁，．．．，ｃ_N-1を、式１ａを用いて計算する処理で
ある。入力データｆ_j、出力データｃ_kは、実数データで
あっても複素データであってもよい。入力データｆ_j、
出力データｃ_kはそれぞれ実空間のデータ、波数空間の
データと呼ばれることがある。すなわち、入力データｆ
_jの添え字jは、一次元の実空間の格子点の座標を表し、
出力データｃ_kの添え字ｋは、一次元の波数空間の格子
点の座標を表すと考えることができる。言い換えると、
上記の式によるフーリエ変換は、一次元の実空間のデー
タを一次元の波数空間のデータに変換する処理である。
したがって、本明細書では、入力データｆ_jの添え字jを
一次元実空間の格子点座標あるいは単に座標と呼び、出
力データｃ_kの添え字ｋを一次元波数空間の格子点の座
標あるいは単に座標と呼ぶことがある。あるいは、それ
らのデータはその座標を有すると呼ぶことがある。しか
し、入力データｆ_j、出力データｃ_kが実際にはそのよう
な実空間、波数空間に属するデータでなくてもよい。一
般に、並列計算機を使用して演算を行う場合、できるだ
け多くのプロセッサが互いに並列に動作する時間を増大
するとともに、プロセッサ間のデータ通信の総回数を少
なくすることが望ましいことが知られている。データ通
信は、プロセッサ内部の計算時間に比べて時間が掛かる
上に、通信はあるプロセッサからのデータの送信と他の
プロセッサでの受信となからなり、受信側のプロセッサ
では、ある処理を実行する前に他のプロセッサからそこ
での演算結果データを受信するようにプログラムされた
場合、そのプロセッサは、受信すべき演算結果データが
受信されるまで、その処理を開始することができない。
したがって、各プロセッサでは、通信の発生に伴い、受
信待ち時間が増大し、他のプロセッサと並列に動作する
時間が減少する。したがって、並列計算機で演算を高速
に行うには、プロセッサ間の通信の総回数を減らすこと
が望ましいことが知られている。このことは演算として
フーリエ変換を並列計算機で実行する場合も同じであ
る。このためには、演算で使用するデータをどのプロセ
ッサに割り当てるか、いつプロセッサ間で演算結果デー
タを交換するかが重要な問題である。並列計算機でフー
リエ変換を行うには、従来の転置アルゴリズムでは、変
換対象データｆ_jを以下のようにして３次元の実空間の
データに写像し、それを用いて変換対象データを割り当
てるプロセッサを決定することになる。いま、ＮＸ，Ｎ
Ｙ，ＮＺをＮＸ＊ＮＹ＊ＮＺ＝Ｎを満たす正の整数と
し、１次元の添字ｊ，ｋを３次元の添字（ｊ_x，ｊ_y，ｊ
_z）、（ｋ_x，ｋ_y，ｋ_z）に次の式２、３によって置換す
る。

【数２】 ｊ=ｊ_y*NX*NZ+ｊ_x*NZ+ｊ_z （ただし、ｊ_x=0,1,...,NX-1， ｊ_y=0,1,...,NY-1， ｊ_z=0,1,...,NZ-1） ．．．（２）

【数３】 ｋ=ｋ_z*NX*NY+ｋ_x*NY+ｋ_y （ただし、ｋ_x=0,1,...,NX-1， ｋ_y=0,1,...,NY-1， ｋ_z=0,1,...,NZ-1） ．．．（３） ここで、記号＊は乗算を表す。この置換は、１次元実空
間の格子点座標ｊ、１次元の波数空間の格子点座標ｋ
を、それぞれ３次元の実空間の格子点座標（ｊ_x，ｊ_y，
ｊ_z）、３次元の波数空間の格子点座標（ｋ_x，ｋ_y，
ｋ_z）に写像することであるとも言える。３次元の実空
間の座標（ｊ_x，ｊ_y，ｊ_z）は、１次元の実空間の座標
ｊから次式により計算される。 ｊ_z=MOD(j/NZ,NX) ｊ_y=(j/(NX*NZ))↓ ｊ_z=MOD(j,NZ) ここで、（ ）↓は、括弧内の数値の整数部分のみを
表し、小数点以下を切り捨てることを表す。したがっ
て、ｊが０，１，２，３，，，（Ｎ−１）と変化したと
きに、（ｊ_x，ｊ_y，ｊ_z）は、（０，０，０），（０，
０，１），（０，０，２），（０，０，３），，，
（０，０，ＮＺ−１）と変化し、さらに、（１，０，
０），（１，０，１），（１，０，２），（１，０，
３），，，（１，０，ＮＺ−１）と変化し、この変化を
（ＮＸ−１，０，ＮＺ−１）まで変化した後に、ｊ_yを
１に変えて上記変化を座標（ＮＸ−１，ＮＹ−１，ＮＺ
−１）に達するまで繰り返す。すなわち、一次元の順次
異なる座標点ｊに対応する３次元の座標点（ｊ_x，ｊ_y，
ｊ_z）は、Ｚ方向、Ｘ方向、Ｙ方向の順に変化する。本
明細書では、このようにフーリエ変換の対象となる一次
元実空間のデータｆ₀，ｆ₁，．，ｆ_N-1を３次元実空間
のデータに写像することを、簡単化のために各辺の長さ
がＮＸ，ＮＹ，ＮＺの直方体状に並べるともいう。以上
の置換は、言い換えると、図２に示すように、原点から
まずＺ方向にデータを並べていき、ＮＺ個のデータを並
べ終わったら次はＸ座標を１だけ増やしてデータを並
べ、これを繰り返してＮＸ＊ＮＺ個のデータを並べ終わ
ったら次はＹ座標を１だけ増やしてデータを並べるとい
う操作を行うことと等価である。但し、図２では、Ｎは
５１２であり、ＮＸ，ＮＹ，ＮＺはともに８に等しいと
仮定した。３次元の波数空間の座標（ｋ_x，ｋ_y，ｋ_z）
は、１次元の波数空間の座標ｋから次式により計算され
る。 ｋ_x=MOD(k/NY,NX) ｋ_y=MOD(k,NY) ｋ_z=(k/(NX*NY))↓ したがって、ｋが０，１，２，３，，，（ＮＸ＊ＮＹ＊
ＮＺ−１）と変化したときに、（ｋ_x，ｋ_y，ｋ_z）は、
（０，０，０），（０，１，０），（０，２，０），
（０，３，０），，，（０，ＮＹ−１，０）と変化し、
さらに、（１，０，０），（１，１，０），（１，２，
０），（１，３，０），，，（１，ＮＹ−１，０）と変
化し、この変化を（ＮＸ−１，ＮＹ−１，０）まで変化
した後に、ｋ_zを１に変えて上記変化を座標（ＮＸ−
１，ＮＹ−１，ＮＺ−１）に達するまで繰り返す。すな
わち、一次元の順次異なる座標点ｋに対応する３次元の
座標点（ｋ_x，ｋ_y，ｋ_z）は、Ｙ方向、Ｘ方向、Ｚ方向
の順に変化する。したがって、求めるべきフーリエ変換
係数ｃ_kとそれに対応する３次元の波数空間の座標
（ｋ_x，ｋ_y，ｋ_z）との関係は図３に示すとおりにな
る。但し、図３では、Ｎは５１２であり、ＮＸ，ＮＹ，
ＮＺはともに８に等しいと仮定した。なお，転置アルゴ
リズム自体は、１次元空間のデータ列ｆ_j、それに対す
るフーリエ変換結果データｃ_kを２次元空間のデータｆ
_jx,jy,ckx,kyに変換して行うこともできる。この場合に
は、１次元のフーリエ変換を２次元のフーリエ変換に置
き直すことになる。しかし、ここで記載するように、１
次元空間のデータ列ｆ_j、それに対するフーリエ変換結
果データｃ_kを３次元空間のデータｆ_jx,jy,jzとｃ
_kx,ky,kzに変換してフーリエ変換を行うのは，並列計算
機の個々のプロセッサがベクトル演算器を持つ場合に、
そのベクトル演算器をうまく利用するためである。この
場合には、１次元のフーリエ変換を３次元のフーリエ変
換に置き直すことになる。すなわち，以下の実施の形態
でも述べるように，変換の各ステップでは，Ｘ方向、Ｙ
方向、Ｚ方向のうちのどれかの方向について変換を行
い，残りの２方向のうちの１方向を用いて並列化を行
い，更に残りの１方向を用いてベクトル化を行う。この
ため，データを３つの方向を持つ直方体状に並べる必要
がある。原理的には，式２，３と同様の変換を行って，
データを２次元空間あるいは４次元以上の空間に並べ直
すこともできるが，２次元では並列化とベクトル化の両
方を行うには次元が足りず，また，４次元以上では不要
な次元ができ，その分だけベクトル化対象のループ長が
短くなってしまうので性能的に不利である。そのため，
ベクトル演算器を持つ並列計算機上で高速フーリエ変換
を行うには，データを式２，３のように３次元に並べ直
すことがことが望ましい。このようなデータの変換は、
各プロセッサがベクトル演算器を持たない場合にも演算
の高速化に有効な場合が多い。置換式２，３を使用する
と、式１ａは次のように書き換えられる。

【数４】 ｃ_k ＝ｃ_kx,ky,kz ＝(1/N)Σ_jz=0 ^Nz-1Σ_jx=0 ^Nx-1Σ_jy=0 ^Ny-1ｆ_jx,jy,jz *exp(-2πi(ｋ_z*NX*NY+ｋ_x*NY+ｋ_y） (ｊ_y*NX*NZ+ｊ_x*NZ+ｊ_z)/(NX*NY*NZ)) ＝(1/N)Σ_jz=0 ^Nz-1exp(-2πi((ｋ_z*NX*NY+ｋ_x*NY+ｋ_y)ｊ_z) /(NX*NY*NZ)) *(Σ_jx=0 ^Nx-1exp(-2πi((ｋ_x*NY+ｋ_y)ｊ_x)/(NX*NY)) *(Σ_jy=0 ^Ny-1ｆ_jx,jy,jz*exp(-2πiｋ_yｊ_y/NY))) （ただし、ｋ_x=0,1,...,NX-1， ｋ_y=0,1,...,NY-1， ｋ_z=0,1,...,NZ-1） ．．．（４） さらに、この変換式は次の３式で表される３ステップの
変換をそれらの式の順に順次実行することにより実現さ
れる変換であることが分かる。

【数５】 ｃ'_jx,ky,jz＝Σ_jy=0 ^Ny-1ｆ_jx,jy,jz*exp(-2πiｋ_yｊ_y/NY) （５）

【数６】 ｃ''_kx,ky,jz＝Σ_jx=0 ^Nx-1ｃ'_jx,ky,jz *exp(-2πi(ｋ_x+ｋ_y/NY)ｊ_x/NX) ．．．（６）

【数７】 ｃ_kx,ky,kz＝Σ_jz=0 ^Nz-1ｃ''_kx,ky,jz *exp(-2πi(ｋ_z+ｋ_x/NY+ｋ_y/(NX*NY))ｊ_z／NZ) ．．（７） 式５は、ｊ_x、ｊ_zが特定の値であり、ｊ_yの値が異なる
ＮＹ個の入力データｆ_jx,jy,jzに対するフーリエ変換
を、ｊ_x、ｊ_zが採りうる値の組合わせの数（ＮＸ＊ＮＺ
組）だけ行い、それにより上記二つの実空間座標
（ｊ_x、ｊ_z）の組のひとつにそれぞれ対応する複数（Ｎ
Ｘ＊ＮＺ）組の、３次元の波数空間の一つの座標
（ｋ_y）に関する一次変換結果データ（ｃ'_jx,ky,jz）
（ｋ_y＝０〜ＮＹ−１、ｊ_x＝０〜ＮＸ−１、ｊ_z＝０〜
ＮＺ−１）を得る処理を表す。式６も、複素指数関数の
中にｋ_y／ＮＹという余分な項が入る以外はフーリエ変
換と同じ変換を表し、具体的には、この式は、ｊ_z、ｋ_y
が特定の値であり、ｊ_xの値が異なるＮＸ個の一次変換
結果データｃ'_jx,ky,jzに対してフーリエ変換と類似の
変換を、ｊ_z、ｋ_yが採りうる値の組合わせの数（ＮＹ＊
ＮＺ組）だけ行い、それにより、座標ｊ_zの異なる値に
対する、３次元の波数空間の二つの座標（ｋ_x、ｋ_y）に
関する２次変換結果データ（ｃ''_kx,ky,jz）（ｋ_x＝０
〜ＮＸ−１、ｋ_y＝０〜ＮＹ−１、ｊ_z＝０〜ＮＺ−１）
を得る処理を表す。式７も、複素指数関数の中にｋ_x／
ＮＹ＋ｋ_y／（ＮＸ＊ＮＹ）という余分な項が入る以外
はフーリエ変換と同じ変換を表し、具体的には、この式
は、ｋ_x、ｋ_yが特定の値であり、ｊ_zの値が異なるＮＺ
個のデータｃ''_kx,ky,jzに対してフーリエ変換と類似の
変換を、ｋ_x、ｋ_yが採りうる値の組合わせの数（ＮＸ＊
ＮＹ組）だけ行い、それにより３次元の波数空間の３つ
の座標（ｋ_x,ｋ_y,ｋ_z）に関する最終的なフーリエ変換
結果データ（ｃ_kx,ky,kz）（ｋ_x＝０〜ＮＸ−１、ｋ_y＝
０〜ＮＹ−１、ｋ_z＝０〜ＮＺ−１）を得る処理を表
す。したがって、これらの３つの変換は、すべて高速フ
ーリエ変換のアルゴリズムを用いて実行することができ
る。以下では、これらの変換をそれぞれＹ方向の変換、
Ｘ方向の変換、Ｚ方向の変換と呼ぶ。本明細書では、こ
れらの変換を簡単化のためにそれぞれＹ方向のフーリエ
変換、Ｘ方向のフーリエ変換、Ｚ方向のフーリエ変換と
呼ぶこともある。ここで、Ｘ方向等は、式２、３で定め
た座標変換できまる方向である。すなわち、一次元の順
次異なる座標点ｊに対応して最初に順次変化する座標が
Ｚ座標であり、その後に変化する座標がＸ座標であり、
最後に変化する座標がＹ座標である。座標変換式を式
２、３から変更すれば、Ｘ方向等の変換の内容が変わる
のは言うまでもない。したがって、本明細書では、より
一般的には、これらの変換は以下の変換を指す。Ｙ方向
の変換とは、式５により例示されたように、実空間の第
１、第３の座標軸の座標（ｊ_x，ｊ_z）が特定の値であ
り、第２の座標軸の座標（ｊ_y）の値が異なる複数（Ｎ
Ｙ）個の入力データｆ_jx,jy,jzに対してフーリエ変換を
行い、３次元実空間の第１、第３の座標軸の座標
（ｊ_x，ｊ_z）の組のひとつにそれぞれ対応する複数（Ｎ
Ｘ＊ＮＺ）組の、３次元波数空間の第２の座標軸の座標
（ｋ_y）に関連する一次変換結果データ（ｃ'_jx,ky,jz）
（ｊ_x＝０〜ＮＸ−１、ｋ_y＝０〜ＮＹ−１、ｊ_z＝０〜
ＮＺ−１）を得る処理を指す。あるいは、言い換える
と、Ｙ方向の変換は、入力データｆ_jx,jy,jzに対して、
実空間の第２の座標軸に関してフーリエ変換を行い、３
次元波数空間の第２の座標軸の座標と、３次元実空間の
第１の座標軸の座標と、第３の座標軸の座標との関数で
ある一次変換結果データを得る処理を指すとも言える。
さらに、Ｘ方向の変換とは、式６により例示されたよう
に、上記第３の実空間座標系の第１の座標系の座標（ｊ
_z）と３次元波数空間の第２の座標系の座標（ｋ_y）とが
特定の値であり、上記第１の実空間座標系の座標
（ｊ_x）の値が異なる複数（ＮＸ）個の一次変換結果デ
ータ（ｃ'_jx,ky,jz）に対してフーリエ変換に類似の変
換を行い、上記第３の実空間座標軸の座標（ｊ_z）の異
なる値の一つにそれぞれ対応する複数（Ｎ_z）個の、上
記３次元の波数空間の第１、第２の座標軸の座標
（ｋ_x、ｋ_y）に関連する２次変換結果データ（ｃ''
_kx,ky,jz）（ｋ_x＝０〜ＮＸ−１、ｋ_y＝０〜ＮＹ−１、
ｊ_z＝０〜ＮＺ−１）を得る処理を指す。あるいは、言
い換えると、Ｘ方向の変換は、一次変換結果データに対
して、３次元実空間の第１の座標軸に関してフーリエ変
換に類似の変換を行い、３次元実空間の第３の座標軸の
座標と、３次元波数空間の第１、第２の座標軸の座標と
の関数である２次変換結果データを得る処理を指すとも
言える。さらに、Ｚ方向の変換とは、式７により例示さ
れたように、３次元波数空間の第１、第２の座標系の座
標（ｋ_x，ｋ_y）とが特定の値であり、３次元実空間の第
３の座標系の座標（ｊ_z）の値が異なる複数（ＮＺ）個
の２次変換結果データ（ｃ''_kx,ky,jz）に対してフーリ
エ変換に類似な変換を行い、３次元波数空間の第１、第
２、第３の座標軸の座標（ｋ_x，ｋ_y，ｋ_Z）に関連す
る、入力データに対する最終的なフーリエ変換結果デー
タ（ｃ_kx,ky,kz）（ｋ_x＝０〜ＮＸ−１、ｋ_y＝０〜ＮＹ
−１、ｋ_z＝０〜ＮＺ−１）を得る処理を指すとも言え
る。あるいは、言い換えると、Ｚ方向の変換は、２次変
換結果データに対して、３次元実空間の第３の座標軸に
関してフーリエ変換に類似の変換を行い、３次元波数空
間の第１、第２、第３の座標軸の座標の関数である最終
的なフーリエ変換結果データを得る処理を指すとも言え
る。Ｙ方向の変換は、式５にて示されるように、ＮＸ＊
ＮＺ組のデータに対する互いに独立な変換からなる。同
様に、Ｘ方向の変換は、式６にて示されるように、ＮＹ
＊ＮＺ組のデータに対する互いに独立な変換からなる。
同様に、Ｚ方向の変換は、式７にて示されるように、Ｎ
Ｘ＊ＮＹ組のデータに対する互いに独立な変換からな
る。従来の転置アルゴリズムによるフーリエ変換方法で
は、この特徴を利用してプロセッサ間の通信を少なくす
るように、変換対象データを並列計算機の異なるプロセ
ッサに割り当てている。すなわち、式２に従って、ま
た、図２に例示されるように、変換対象データｆ_j（j＝
０〜Ｎ）を直方体状に並べ、３次元実空間のＺ軸に並行
な平面でこのデータを分割し、図５に例示するように、
j_z＝０，１，，，７をそれぞれ有する複数のデータはプ
ロセッサ０，１，，７に割り当てられている。すなわ
ち、特定の値のＺ座標j_zを有する全ての変換対象データ
は、それらのＸ座標j_x、Ｙ座標j_yの値に依らないで同一
のプロセッサに割り当てる。図２では、Ｎ＝５１２，Ｎ
Ｘ＝ＮＹ＝ＮＺ＝８と仮定し、図４ではプロセッサの総
数ＮＰＵ＝８と仮定したが、これらの数値がここに仮定
の数値と異なる場合でも、特定の値のＺ座標j_zを有する
全ての変換対象データは、それらのＸ座標、Ｙ座標の値
に依らないで同一のプロセッサに割り当てればよい。た
とえば、プロセッサの総数ＮＰＵ＝ＮＸ＝ＮＺとし、Ｎ
Ｙ＝（Ｎ／（ＮＸ＊ＮＺ））↑とすればよい。ここで、
（ ）↑は、括弧内の数値の小数点以下を切り上げた後
の整数を示す。たとえば、Ｎ＝５１２、ＮＰＵ＝４のと
きには、ＮＸ＝ＮＺ＝４、ＮＹ＝３２であればよい。こ
のようなデータの割り当てを行った後、フーリエ変換を
以下のように実行する。この方法では式５によるＹ方向
の変換と式６によるＸ方向の変換とは、プロセッサ間の
通信を使用しないで行うことができる。 （ステップ１）Ｙ方向の変換 まず、Ｙ方向の変換を実行する。式５から分かるよう
に、Ｙ方向の変換では、Ｘ座標j_xとＺ座標j_zが特定の値
であり、Ｙ座標j_yが異なる複数の変換対象データに対し
てフーリエ変換が実行される。しかし、これらの複数の
データは同じプロセッサに割り当てられている。こうし
て、各プロセッサでは、式５にしたがって、プロセッサ
間の通信を使用しないで、 Ｘ座標j_xとＺ座標j_zが特定
の値であり、波数空間のＹ座標ｋ_yがいろいろの値を有
する複数の一次変換結果データが得られる。 （ステップ２）Ｘ方向の変換 次に、Ｘ方向の変換を実行する。式６から分かるよう
に、Ｘ方向の変換では、Ｚ座標j_zと波数空間のＹ座標ｋ
_yが特定の値であり、Ｘ座標j_xが異なる複数の一次変換
結果データに対してフーリエ変換に類似の変換が実行さ
れる。しかし、これらの複数のデータは同じプロセッサ
でのＹ方向の変換によりすでに得られている。こうし
て、各プロセッサは、式６にしたがって、他のプロセッ
サとの通信をしないで、そのプロセッサに割り当てられ
た、３次元実空間のＺ座標j_zの特定の値に対する、３次
元波数空間の座標（ｋ_x，ｋ_y）に関連するＮＸ＊ＮＹ個
の２次変換結果データ（ｃ''_kx,ky,jz）（ｋ_x＝０〜Ｎ
Ｘ−１、ｋ_y＝０〜ＮＹ−１）を得る。 （ステップ３）データの転置 式７によるＺ方向の変換を行うには、３次元波数空間の
座標（ｋ_x，ｋ_y）の特定の値と、３次元実空間のＺ座標
j_zの全ての値に対して得られた２次変換結果データ
（ｃ''_kx,ky,jz）が必要である。そこで、各プロセッサ
に、３次元波数空間の座標ｋ_xの特定の値を割り当て
て、Ｚ方向の変換を実行するに必要なデータをプロセッ
サ間で転送する処理が実行される。すなわち、各プロセ
ッサへの座標ｋ_xの値の割り当てでは、プロセッサ０，
１，２，，，７には、ｋ_x＝０，１，２，３，，，を順
次割り当てる。このことは、図５に示すように、３次元
波数空間を、そのｋ_x軸に垂直な平面で分割して、分割
後の部分空間の各々を一つのプロセッサに割り当てるこ
とを意味する。各プロセッサが、そのプロセッサに割り
当てられた３次元波数空間の座標ｋ_xの特定の値と、３
次元波数空間の座標ｋ_yの全ての値と、３次元実空間の
Ｚ座標j _zの全ての値とに対して得られた全ての２次変換
結果データ（ｃ''_kx,ky,jz）を使用してＺ方向の変換を
実行できるように、全プロセッサ間で２次変換結果デー
タの転送が行われる。このデータ転送は、データの転置
あるいはデータの並び替えとも言われる。すなわち、各
プロセッサは、３次元実空間のＺ座標j_zの全ての値と、
３次元波数空間の座標（ｋ_x，ｋ_y）の全ての値との組に
対して得られた２次変換結果データ（ｃ''_kx,ky,jz）
（但し、ｋ_x＝特定値、ｋ_y＝０〜ＮＹ―１、ｋ_z＝０〜
（ＮＺ―１））の内、自プロセッサが生成しなかったデ
ータを他のプロセッサから受信するように、全プロセッ
サの間で２次変換結果データを転送する。 （ステップ４）Ｚ方向の変換 各プロセッサは、そのプロセッサに割り当てられた３次
元波数空間の座標ｋ_xの特定の値と、３次元波数空間の
他の二つの座標ｋ_y，ｋ_zの全ての値を有する最終的なフ
ーリエ変換結果データｃ_kx,ky,kz（但し、ｋ_x＝特定
値、ｋ_y＝０〜ＮＹ―１、ｋｚ＝０〜ＮＺ―１）を計算
する。各プロセッサはこの計算を互いに並列に実行でき
る。 （ステップ５）データの転置 しかし以上の処理だけでは、変換対象データｆ_jと変換
結果データｃ_kのデータ分割形式が異なり、実用上不便
であるという問題がある。すなわち、変換対象データｆ
_jは、図２に示すように、３次元実空間のデータに写像
され、後者のデータは、図４にしたがって分割されて複
数のプロセッサに割り当てられた。したがって、図６に
示すように、変換対象データｆ_jの分割は、ｆ_jを第ＭＯ
Ｄ（j，ｐ）番のプロセッサが担当するサイクリック分
割となる。一方、変換結果データｃ_kは、図３に示すよ
うに、３次元波数空間のデータに写像され、３次元波数
空間は図５にしたがって分割されて複数のプロセッサに
割り当てられた。したがって、図７に示すように、変換
結果データｃ_kのデータ分割は、ＮＹ個の連続するデー
タを１台のプロセッサが担当するブロックサイクリック
分割となる。多くの応用では、変換対象データをフーリ
エ変換して得られる変換結果データに対してある処理を
施し、その処理結果に対して再び逆フーリエ変換を行
う。逆フーリエ変換はフーリエ変換のプログラムを流用
して行われることが多い。すなわち、次の逆フーリエ変
換の式

【数８】ｆ_j=Σ_k=0 ^N-1ｃ_kexp(2πikj/N) （ただし、ｊ=0,1,...,N-1） ．．．（８） を変形すると次式が得られる。

【数９】ｆ_j=(Σ_k=0 ^N-1ｃ_k ^*exp(-2πikj/N))^* （ただし、ｊ=0,1,...,N-1） ．．．（９） ここで、＊印は複素共役を示す。したがって、式１ａと
９の比較より明らかなように、逆フーリエ変換は、フー
リエ変換結果データの複素共役をフーリエ変換し、得ら
れた結果データの複素共役を取ることに等価である。し
たがって、原理的には、逆フーリエ変換はフーリエ変換
のプログラムを流用して実行できることが分かる。しか
し、並列計算機でフーリエ変換を実行するときには、変
換対象データと変換結果データとのデータ分割形式が異
なると、いずれかのプロセッサに割り当てられた変換対
象データに対する変換結果データがそのプロセッサに割
り当てられていないことになり、そのプロセッサは、そ
の変換結果データを他のプロセッサから受信しないとフ
ーリエ変換のプログラムを流用して逆フーリエ変換を行
うことができなくなる。このような不便を避けるため、
従来の転置アルゴリズムを使ったフーリエ変換プログラ
ムでは、上記Ｚ方向の変換を実行した後に、３次元フー
リエ変換結果データｃ_kx,ky,kzのプロセッサへの割り当
てを変更し、再びプロセッサ間でフーリエ変換結果デー
タｃ_kx,ky,kzの転置（入れ替え）を行い、フーリエ変換
結果データｃ_kx,ky,kzのデータ分割をサイクリック分割
に直すのが一般的であった。すなわち、図８に示すよう
に、３次元波数空間を、Ｙ座標軸に垂直な平面で切断
し、同じＹ座標値ｋ_yを有するフーリエ変換結果データ
ｃ_kx,ky,kzを同一のプロセッサに割り当てる。具体的に
は、ｋ_y＝０，１，２，，，を有するフーリエ変換結果
データｃ_kx,ky,kzを順次プロセッサ０，１，２，，，に
割り当てる。各プロセッサが、この割り当てにしたがっ
てそのプロセッサに割り当てられたフーリエ変換結果デ
ータｃ_kx,ky,kzの内、自プロセッサが生成しなかったデ
ータを他のプロセッサから受信するように、全プロセッ
サの間でフーリエ変換結果データｃ_kx,ky,kzを転送す
る。こうして、サイクリック分割された３次元フーリエ
変換結果データが得られる。こうして、得られた最終的
３次元フーリエ変換結果データｃ_kx,ky,kzから目的とす
る１次元フーリエ変換結果データｃ_kは式３より得るこ
とができる。データｃ_kとその三次元座標ｋ_x，ｋ_y，ｋ_z
との関係は図３に示されたとおりである。しかし、本発
明者による検討によれば、従来のデータのデータ入れ替
えのためのプロセッサ間での余分な通信は、並列化効率
を低下し、フーリエ変換に必要な処理時間を増大する原
因であることが判明した。そこで本発明では、従来の転
置アルゴリズムにおけるデータの分割方式を見直し、プ
ロセッサ間のデータ転送量の削減の側面から最適なデー
タ分割方式を以下のようにして決定した。従来の転置ア
ルゴリズムは、Ｙ方向、Ｘ方向、Ｚ方向の各変換と、プ
ロセッサ間でのデータの入れ替えを行う転置操作から構
成される。そのアルゴリズムでは、Ｙ方向の変換を行う
際に、変換対象データの空間をＺ軸に垂直な複数の平面
で切り、各面を一つのプロセッサに割り当てて、次にＸ
方向の変換を行う際には、その割り当てをそのまま使用
し、さらにＺ方向の変換を行う際には、Ｘ軸に垂直な複
数の平面で変換対象データの空間を切り、各面を一つの
プロセッサに割り当てていた。これにより、その変換そ
のものは、プロセッサ間のデータ転送なしに実行でき
た。しかし、たとえば最初のＹ方向の変換を行う際に
は、変換の対象となる同一のＸ座標とＺ座標を持つＮＹ
個のデータが１台のプロセッサ上にありさえすれば、そ
の変換そのものは、プロセッサ間のデータ転送なしに実
行できる。したがって、この変換対象データの分割は、
Ｚ軸に垂直な平面によってではなく、Ｘ軸に垂直な平面
によって行ってもよい。このことは、Ｘ方向、Ｚ方向の
変換についても言える。したがって、望ましいデータ分
割方式が満たすべき第一の条件は、「ある方向の変換を
行うときには、その変換の変換対象データをその方向以
外の方向に垂直な複数の平面で分割してプロセッサへの
データ割り当てをする」というデータ分割形式が、Ｙ方
向、Ｘ方向、Ｚ方向のすべてに採用されていることであ
る。今ひとつ考慮すべきことは転置のためのデータ転送
回数である。たとえば、Ｙ方向、Ｘ方向、Ｚ方向の変換
を行うとき、変換対象データをそれぞれＺ軸、Ｙ軸、Ｘ
軸に垂直な複数の平面で切って分割するというデータ分
割方式は、上記の第１の条件を満たすが、このデータ分
割方式では、データ分割形式がＸ方向、Ｚ方向の変換を
行うときという２回にわたって変更され、その変更の度
に転置のためのデータ転送が必要になる。したがって、
望ましいデータ分割方式が満たすべき条件として、上記
の第１の条件に加えて、「データ分割形式の変更のため
の転置処理は一回に限る」という第二の条件を付加す
る。これら２つの条件を満たすデータ分割方式を数え上
げた結果を図９に示す。上記第１、第２の条件を満たす
データ分割方式は４通りあり、これらの内で、フーリエ
変換対象データのデータ分割形式とフーリエ変換結果デ
ータのデータ分割形式が同一のデータ分割方式が求める
ものである。方式１が従来の転置アルゴリズムで採用さ
れているものである。方式４は入力データがＸ方向に分
割、出力データがＹ方向に分割だから、図３および図６
と照らし合わせてみると、フーリエ変換対象データがブ
ロックサイクリック分割、フーリエ変換結果データがサ
イクリック分割であり、方式１とはデータ分割形式がち
ょうど逆になってはいるものの、これらの二つの種類の
データの間でデータ分割形式が異なるという方式１と同
様の欠点を抱えていることがわかる。一方、方式２で
は、従来の転置アルゴリズムと同じく、フーリエ変換対
象データがＺ方向に沿って分割され、Ｙ方向の変換、Ｘ
方向の変換も従来の転置アルゴリズムと同じように実行
されるが、Ｚ方向の変換は、従来の転置アルゴリズムと
異なり、Ｘ方向の変換の結果データがＹ方向に沿って分
割された後に実行される。Ｘ方向の変換とＺ方向の変換
の間では、データ転置が必要である。図２および図３と
照らし合わせてみると、フーリエ変換対象データもフー
リエ変換結果データもサイクリック分割になることが分
かる。したがって、方式２のデータ分割方式を採用する
と、フーリエ変換係数の計算後にプロセッサ間でデータ
転送を行わなくても、フーリエ変換対象データとフーリ
エ変換結果データとが同じデータ分割形式を保つ。この
方式２では、フーリエ変換は具体的には以下のようにし
て実行される。以下に記載するＹ方向、Ｘ方向、Ｚ方向
の変換はＦＦＴのアルゴリズムにより計算される。 （ステップａ）Ｙ方向の変換 Ｙ方向の変換は、従来の転置アルゴリズムについて既に
述べたステップ１の要領で実行される。既に述べたごと
く、フーリエ変換対象データのデータ分割は、サイクリ
ック分割である。 （ステップｂ）Ｘ方向の変換 さらに、Ｘ方向の変換は、Ｙ方向の変換の結果データに
対して、従来の転置アルゴリズムについて既に述べたス
テップ２の要領でなされる。 （ステップｃ）データ転置 式７によるＺ方向の変換を行うには、３次元波数空間の
座標（ｋ_x，ｋ_y）の特定の値と、３次元実空間のＺ座標
j_zの全ての値に対して得られた２次変換結果データ
（ｃ''_kx,ky,jz）が必要である。方式２では、従来の転
置アルゴリズムと異なり、Ｘ方向の変換で得られた２次
変換結果データ（ｃ''_kx,ky,jz）は、Ｙ軸に垂直な複数
の平面で分割される。このことは、図８に示すように、
３次元波数空間を、そのｋ_y軸に垂直な複数の平面で分
割して、分割後の部分空間（上記複数の平面）の各々を
一つのプロセッサに割り当てることを意味する。すなわ
ち、各プロセッサへ座標ｋ_yの特定の値を割り当てる。
具体的には、ｋ_x＝０，１，２，３，，，の２次変換結
果データ（ｃ''_kx,ky,jz）を順次プロセッサ０，１，
２，，，７に割り当てる。この割り当てに従い、Ｚ方向
の変換を実行するに必要なデータをプロセッサ間で転送
する処理が実行される。各プロセッサが、そのプロセッ
サに割り当てられた３次元波数空間の座標ｋ_yの特定の
値と、３次元波数空間の座標ｋ_xの全ての値と、３次元
実空間のＺ座標j_zの全ての値とに対して得られた全ての
２次変換結果データ（ｃ''_kx,ky,jz）を使用してＺ方向
の変換を実行できるように、全プロセッサ間で２次変換
結果データの転送が行われる。すなわち、各プロセッサ
は、３次元実空間のＺ座標j_zの全ての値と、３次元波数
空間の座標ｋ_xの全ての値と、座標ｋ_yの特定の値との組
に対して得られた２次変換結果データ（ｃ''_kx,ky,jz）
（但し、ｋ_x＝０〜（ＮＸ―１）、ｋ_y＝特定値、ｋ_z＝
０〜（ＮＺ―１））の内、自プロセッサが生成しなかっ
たデータを他のプロセッサから受信するように、全プロ
セッサの間で２次変換結果データを転送する。 （ステップｄ）Ｚ方向の変換 各プロセッサは、そのプロセッサに割り当てられた３次
元波数空間の座標ｋ_yの特定の値と、３次元波数空間の
他の二つの座標ｋ_x、ｋ_zの全ての値を有する最終的なフ
ーリエ変換結果データｃ_kx,ky,kz（但し、ｋ_x＝０〜Ｎ
Ｘ―１、ｋ_y＝特定値、ｋ_z＝０〜ＮＺ―１）を計算す
る。各プロセッサはこの計算を互いに並列に実行でき
る。この結果、座標ｋ_y＝０，１，２，３，，，に対応
するフーリエ変換係数ｃ₀，ｃ₁，ｃ₂，，，が、順次プ
ロセッサ０，１，２，，，で生成され、全フーリエ変換
結果データｃ_kx,ky,kzは、プロセッサ間でサイクリック
に分割されていることが分かる。本実施の形態では、上
記データ分割方式２を使用する。なお、データ分割方式
３も後に詳細に説明するように、フーリエ変換対象デー
タもフーリエ変換結果データもサイクリック分割になっ
ている。したがって、この方式３も使用することができ
る。後に述べるように、実際にフーリエ変換ライブラリ
を構成する場合に、方式３は、方式２に比べて、ライブ
ラリが生成する複素指数関数の値のテーブルのサイズが
小さくてよいという利点を有する。なお、計算機により
入力データｆ_jに対して以上の変換を実施するときに
は、一般に配列データが使用される。すなわち、同じプ
ロセッサでＹ方向の変換を施すべき一群のデータは、３
次元配列に格納され、その配列に対してＹ方向の変換が
実行される。その結果得られた１次変換結果データは、
同じ配列あるいは他の３次元配列に格納されてもよい。
他の方向の変換も先に実行された変換の結果データを格
納する配列に対してなされる。また、プロセッサ間での
データの転置も各プロセッサが生成した配列の内容を交
換するようになされる。したがって、このようにそれら
の変換において同じ３次元配列を使用するときには、そ
の３次元配列の各次元のインデックスは、あるときには
３次元実空間の各座標軸に対応し、他の時にはある方向
の変換後の結果データが属する３次元空間の各座標軸に
対応し、最終的にはフーリエ変換係数が属する３次元波
数空間の各座標軸に対応することになる。しかし、この
ように同じ配列を異なる種類の一群のデータの格納に使
用された場合でも、ある時点でその配列に格納されてい
る一群のデータは、その一群のデータが属する３次元空
間に属し、その配列の各インデックスは、その一群のデ
ータが属する３次元空間のいずれかの座標軸を表すこと
には変わりはない。したがって、本発明を実施するにあ
たって一群のデータを格納するのに使用する配列の具体
的な構造は特定のもの限定されない。さらに、以上の原
理で説明したいくつかの３次元空間のいずれか一つに属
する一群のデータを格納する配列の構造が、その３次元
空間に直接対応しないものであっても、その配列に含ま
れた各データは、その３次元空間の座標を有すると見な
すことができ、以上の原理説明がその配列に対してもあ
てはまるのは言うまでもない。 (3) 多次元フーリエ変換への応用 以上、１次元フーリエ変換のためのアルゴリズムを説明
したが、本方式は多次元フーリエ変換の場合へも簡単に
拡張できる。次式の２次元フーリエ変換を例に採って説
明する。

【数１０】 ｃ_k1,k2＝(1/N₁N₂)Σ_j1=0 ^N1-1Σ_j2=0 ^N2-1ｆ_j1,j2 *exp(-2πi(k₁ｊ₁/N₁+ｋ₂ｊ₂/N₂)) （ただし、ｋ₁=0,1,... ,N₁-1,ｋ₂=0,1,...,N₂-1) ．．．（１０） この式は、次式１１に変形できる。この式１１は、更に
次の２ステップからなる変換式１１ａ、１１ｂとして書
くことができる。

【数１１】 ｃ_k1,k2＝(1/Ｎ₂)Σ_j2=0 ^N2-1exp(-2πiｋ₂ｊ₂/N₂) *(1/Ｎ₁)Σ_j1=0 ^N1-1ｆ_j1,j2 *exp(-2πiｋ₁ｊ₁/Ｎ₁) ．．．（１１） ｃ'_k1,j2＝(1/Ｎ₁)Σ_j1=0 ^N1-1ｆ_j1,j2 *exp(-2πiｋ₁ｊ₁/Ｎ₁) ．．（１１ａ） ｃ_k1,k2＝(1/Ｎ₂)Σ_j2=0 ^N2-1ｃ'_k1,j2 *exp(-2πiｋ₂ｊ₂/Ｎ₂) ．．．（１１ｂ） したがって、２次元フーリエ変換は、まず式１１ａのよ
うにＮ₁個のデータに対する１次元フーリエ変換をＮ₂組
行い、次に式１１ｂのようにＮ₂個のデータに対する１
次元フーリエ変換をＮ₁組行うことに帰着する。したが
って、これらの１次元フーリエ変換において、本発明の
方式を適用できる。本発明の方式を用いて並列計算機上
で２次元フーリエ変換を行うには、２次元データｆ
_j1,j2に対し、添え字ｊ₁の方向（以下、これを第１方向
と呼ぶ）にサイクリック分割を行う。すなわち、第i番
目のプロセッサにｆ_m*NPU+i,j2（但し、m=0,1,...,((N1
/NPU)-1)、j2=0,1,...,N2)番目の要素を割り当てる。こ
こで、ＮＰＵはプロセッサの台数である。すると、式１
１ａのステップは、ｊ₂が同じＮ₁個の要素の間での１次
元フーリエ変換をＮ₂組行うことであり、このＮ₁個の要
素はプロセッサ間にサイクリック分割されているから、
このステップは本発明の方式による１次元フーリエ変換
をＮ₂組行うことに帰着する。変換後のデータｃ'_k1,j2
は、第１方向にサイクリック分割されている。次に式１
１ｂのステップでは、ｊ₁が同じＮ₂個の要素の間での１
次元フーリエ変換をＮ₁組行うが、これらＮ₂個の要素は
同一プロセッサ上にあるため、この変換は通信なしに各
プロセッサごとに独立に行える。以上により２次元フー
リエ変換が完了し、変換後のデータｃ_k1,k2は第１方向
にサイクリック分割される。なお、以上では２次元の場
合を示したが、より次元の大きい場合も、第１方向にサ
イクリック分割を行い、第１方向の変換のみを本発明の
方式を用いて行い、以下の変換はプロセッサごとに独立
に行うことにより、本発明のフーリエ変換方式を適用可
能である。 （４）並列高速フーリエ変換ライブラリ 図１に戻り、並列計算機２８上で使用される高速フーリ
エ変換ライブラリは、具体的には、たとえば以下のよう
に構成される。但し、本発明を適用したフーリエ変換ラ
イブラリは、これに限定されないことは言うまでもな
い。本ライブラリは、全てのプロセッサにロードされ、
そのプロセッサ内のシミュレーションプログラムから必
要に応じてサブルーチンとしてコールされる。サブルー
チン名称をＦＦＴ１Ｄとし、これを実行するには CALL FFT1D (NX, NY, NZ, NPU, F, TB, IOPT, IER) のように所定の引数を指定して、いずれかのプロセッサ
にロードされたすべてのシミュレーションプログラムか
ら同時にコールする必要がある。ここで、Ｎ＝ＮＸ＊Ｎ
Ｙ＊ＮＺはフーリエ変換対象データｆ_jの個数、ＮＰＵ
はプロセッサ台数、Ｆはライブラリのコール時はフーリ
エ変換対象データｆ_j、ライブラリからのリターン時は
フーリエ変換結果データｃ_kを格納する配列、ＴＢは複
素指数関数の値を格納するテーブル、ＩＯＰＴはサブル
ーチンの機能を指定する入力、ＩＥＲは実行時エラーが
生じたか否かを示す出力である。ここで、配列Ｆは各プ
ロセッサがそれぞれ持つ部分配列である。フーリエ変換
の原理説明で説明したように、全入力データ（フーリエ
変換対象データ）ｆ_jは、図２のように３次元実空間に
直方体状に配置され、各プロセッサには、この直方体の
内、一つまたは複数の特定のＺ座標を有するＺ軸に垂直
な平面に属する入力データが割り当てられる。この割り
当てられた入力データが、上記引数Ｆで指定される部分
配列に格納されている。すなわち、フーリエ変換対象デ
ータとフーリエ変換結果データは、ともにサイクリック
分割されるので、第ｉ番目のプロセッサは、ｍ＊ＮＰＵ
＋ｉ（ｍ＝０，１，．．．，Ｎ／ＮＰＵ−１）番目の要
素のみを持つ。すなわち、第ｉ番目のプロセッサの配列
Ｆには、Ｎ個の入力データ列ｆ_j（ｊ＝０〜Ｎ−１）の
内、次式で示されるように、一群の入力データｆ
_m*NPU+iを格納する。 Ｆ（ｍ）＝ｆ_m*NPU+i(m=0,1,...,N/NPU- 1) したがって各プロセッサの持つ配列Ｆの大きさはＮ／Ｎ
ＰＵである。また、ＴＢは、第１回目のコールで計算し
た複素指数関数の値を格納しておくテーブルであり、２
回目のコールからはここに格納した値を再利用すること
により、新たな計算が不要となる。また、第１回目のコ
ールではＩＯＰＴ＝１を指定し、このときは複素指数関
数のテーブルを作成する。ＩＯＰＴ＝２は２回目以降の
コールを意味し、このときは既にＴＢに格納されている
値を用いる。本ライブラリのフローチャートを図１０に
示す。本ライブラリは、コールされると (ステップ４
５) 、まず引数をチェックする（ステップ４６）。す
なわち、Ｎ＝ＮＸ＊ＮＹ＊ＮＺとＮＰＵとが１以上の整
数であるかどうか、ＩＯＰＴが１または２の値であるか
どうかなど、引数の有効性を調べる。入力データに無効
な値が入っていた場合は、ＩＥＲ＝１０００と設定して
（処理４７）リターンする。次に、他のプロセッサに本
ライブラリがコールされたことを通知する（ステップ４
８）。この通知は、実際には、そのライブラリがロード
されている通信ライブラリに、他の全てのプロセッサに
当該プロセッサでのライブラリコールの発生を通知する
ことを要求し、その通信ライブラリが、その発生を他の
全てのプロセッサに通知するメッセージを送信し、それ
ぞれの他のプロセッサでは、そこにロードされた通信ラ
イブラリが、このメッセージを受信して、そのプロセッ
サでロードされた本ライブラリに、送信元のプロセッサ
での本ライブラリのコールを通知する。次に、ライブラ
リが引数で指定した通りにＮＰＵ台のプロセッサでコー
ルされているかどうかをチェックする（ステップ４
９）。このチェックは、上に述べた他の全てのプロセッ
サから本ライブラリに対するライブラリコールが発生し
たとの通知を受信したか否かに基づいて行われる。この
条件が満たされていない場合は、ＩＥＲ＝２０００と設
定して（ステップ５０）リターンする。次にＩＯＰＴの
値をチェックし（ステップ５１）、ＩＯＰＴ＝１の場合
は、現在のコールが、最初のコールである。したがっ
て、そのコール元のプロセッサでのフーリエ変換を実行
するための準備を行う。具体的には、Ｘ、Ｙ、Ｚの方向
の変換のために、そのプロセッサで式５，６，７で使用
する複素指数関数の値を前もって計算し、複素指数関数
のテーブルを生成し、配列ＴＢとして格納する（ステッ
プ５５）。計算すべき複素指数関数の値は、そのプロセ
ッサに対するデータの割り当てにより定まる。すなわ
ち、Ｘ、Ｙ、Ｚの方向の変換の各々においてそのプロセ
ッサが処理すべきデータの３次元実空間の座標ｊ_x，
ｊ_y，ｊ_zと３次元実空間の座標ｋ_x，ｋ_y，ｋ_zとを決定
し、この結果により、式５から７の複素指数の偏角が採
りうるいろいろの値を決定し、それぞれの偏角に対する
余弦関数の値と正弦関数の値を計算し、配列ＴＢに格納
する。上記決定では、各方向での変換に使用されるデー
タ分割形式とそのコール元のプロセッサに予め割り当て
られたプロセッサ番号と、式２、３が使用される。この
プロセッサ番号は、シミュレーションプログラムのロー
ド時に予め各プロセッサに並列計算機２８により指定さ
れるものである。各方向での変換に使用されるデータ分
割形式は、使用されるデータ分割方式、本実施の形態で
は前述の方式２、により定まる。なお、ＩＯＰＴ＝１で
ない場合は、現在のコールが、２回目以降のコールであ
る。このようなコールは、ライブラリのコール元のシミ
ュレーションプログラムが、異なる物理量に対するフー
リエ変換を行うようにプログラムされている場合におい
て生じる。たとえば、シミュレーションプログラムが、
第１の物理量に対するフーリエ変換のために本ライブラ
リをコールした後に、第２の物理量に対するフーリエ変
換のために本ライブラリを再度コールした場合である。
この場合、第２の物理量を表すフーリエ変換対象データ
も第１の物理量を表すフーリエ変換対象データと同じ添
え字を有することが多い。この場合には、第２の物理量
に対するフーリエ変換の実行時に、先に配列ＢＴに格納
した複素指数の値が使用できる。したがって、ＩＯＰＴ
＝１でない場合は、ステップ５５を実行しない。次に、
Ｙ方向の変換を行う（ステップ５６）。本実施の形態で
は、全プロセッサが持つフーリエ変換対象データを仮想
的に図２のような各辺の長さが引数ＮＸ，ＮＹ，ＮＺの
直方体状に並べ、図４に示すように、特定の座標を有す
るフーリエ変換対象データを同一のプロセッサに割り当
てられる。このＹ方向の変換では、既にステップ１ある
いはステップａとして述べたように、各プロセッサは、
同じＸ座標とＺ座標とを持つＮＹ個のデータについて、
高速フーリエ変換が式５に従い行う。このようなデータ
の組は全部でＮＸ＊ＮＺ組あるため、結局、ＮＸ＊ＮＺ
個の独立なＮＹ次の高速フーリエ変換を行うことにな
る。プロセッサへのデータの割り当て方式より、各ＸＹ
平面は１台のプロセッサに担当されているから、この変
換処理は通信なしに各プロセッサで独立に行える。本ラ
イブラリの場合、本ライブラリにより各プロセッサが処
理すべき変換対象データは、そのプロセッサで実行され
るシミュレーションプログラムにより、引数Ｆで指定さ
れる配列として、そのプロセッサのメモリ（２６（図
１））にコール前に格納されている。その配列Ｆの添え
字と３次元実空間の座標ｊ_x，ｊ_y，ｊ_zとの関係は、デ
ータ分割形式により定まる。したがって、このＹ方向の
変換では、この関係を使用して配列Ｆ内の変換対象デー
タに対して式５で指定される変換を実行する。変換で得
られた一次変換結果データはそのプロセッサのメモリに
記憶される。具体的には、各プロセッサでは、本ライブ
ラリがコールされると、適当なタイミングで（たとえ
ば、ステップ４６で入力データに無効な値が入っていな
いと判定された時）、各プロセッサは、データ格納用の
３次元配列及び第１、第２、第３の３次元の作業配列を
メモリ上に確保する。ここでは、データ格納用の３次元
配列は、３次元のインデックスの長さが引数ＮＸ、Ｎ
Ｙ、ＮＺに等しい。以下ではこれらの３次元の作業配列
もデータ格納用の３次元配列と同じ大きさを有すると仮
定する。しかし、これらの３次元の作業配列は、以下に
説明するデータを格納できる大きさを有すればよく、し
たがって、これらの作業配列の大きさは適宜変更可能で
ある。さらに、これらの３次元の作業配列の構造も、そ
の使用目的に合致する限り、変更することができる。デ
ータ格納用の３次元配列には、上記引数が指定する配列
Ｆに含まれるデータ点列を以下のようにして格納でき
る。そのプロセッサに、図２のＺ軸に垂直な一つの平面
が割り当てられているときには、その一つの平面に属す
るＮＸ＊ＮＹ個の入力データが、それぞれのデータの
Ｘ、Ｙ、Ｚ座標に対応する、上記データ用３次元配列の
インデックスを有する位置に格納される。そのプロセッ
サにＺ軸に垂直な複数の平面が割り当てられたときに
は、各面のデータは同様にして、上記データ用３次元配
列の対応するインデックスの位置に格納される。各プロ
セッサでは、Ｙ方向の変換はこのデータ格納用の配列を
使用して、Ｚ座標が特定の値を有し、Ｙ座標とＸ座標が
いろいろの値を有する一群の入力データに対して、式５
により行なわれる。このとき、Ｚ座標が特定の値を有
し、Ｘ座標が異なる一群の入力データに対してＹ方向の
変換が高速フーリエ変換アルゴリズム（ＦＦＴ）を用い
て実行される。この変換の実行にあっては、Ｘ座標が異
なる一群のデータに対して、プロセッサ内のベクトル演
算器（図示せず）が使用され、パイプライン的に計算が
実行される。その結果得られる１次変換結果データｃ’
_jx,ky,jz（但し、ｊ_x＝０〜ＮＸ−１，ｋ_y＝０〜ＮＹ−
１，ｊ_z＝特定値）は、第１の３次元の作業配列の、こ
れらの座標値ｊ_x，ｋ_y，ｊ_zに対応するインデックスの
ところに格納される。したがって、図２の場合、各プロ
セッサでは、図２の一つの平面上の一群の入力データに
対する１次変換結果データｃ’_jx,ky,jzが、第１の３次
元作業配列の、特定の座標ｊ_zを有する平面上に格納さ
れる。もし、図２において、Ｚ軸に垂直な複数の平面に
属する入力データがそのプロセッサに割り当てられてい
るときには、それぞれのＺ面に対応する、上記第１の３
次元作業配列内の、Ｚ軸に垂直な複数の平面のそれぞれ
に対応する１次変換結果データｃ'_jx,ky,jzが格納され
る。Ｙ方向の変換の終了後、同様にしてＸ方向の変換を
行う（ステップ５７）。すなわち、各プロセッサは、す
でにステップ２あるいはステップｂとして述べたよう
に、各プロセッサは、Ｙ方向の変換で得られた一次変換
結果データに対して式６で指定される変換を実行する。
変換で得られた２次変換結果データはそのプロセッサの
メモリに記憶される。この変換処理も通信なしに各プロ
セッサで独立に行える。具体的には、各プロセッサで
は、Ｘ方向の変換は、Ｚ座標が特定の値を有し、Ｘ座標
とｋ_y座標とがいろいろの値を有する一群の１次変換結
果データｃ’_jx,ky,jzに対して、式６により行なわれ
る。このとき、Ｚ座標が特定の値を有し、ｋ_y座標が異
なる一群の入力データに対してＸ方向の変換が高速フー
リエ変換アルゴリズム（ＦＦＴ）を用いて実行される。
この変換は、上記第１の３次元作業配列を使用して実行
される。この変換の実行にあっては、ｋ_y座標が異なる
一群のデータに対して、プロセッサ内のベクトル演算器
が使用され、計算がパイプライン的に実行される。その
結果得られる２次変換結果データｃ''_kx,ky,jz（但し、
k_x＝０〜ＮＸ−１，ｋ_y＝０〜ＮＹ−１，ｊ_z＝特定値）
は、第２の３次元作業配列の、これらの座標値ｋ_x，
ｋ_y，ｊ_zに対応するインデックスのところに格納され
る。したがって、図２の場合、各プロセッサでは、２次
変換結果データｃ''_kx,ky,jzは、第２の３次元作業配列
の、特定の座標ｊ_zを有する一つの平面に格納される。
もし、図２において、Ｚ軸に垂直な複数の平面に属する
入力データがそのプロセッサに割り当てられているとき
には、それぞれのＺ面に対応する、上記第２の３次元作
業配列内の、Ｚ軸に垂直な複数の平面のそれぞれに対応
する２次変換結果データｃ''_kx,ky,jzが格納される。Ｘ
方向の変換の終了後、プロセッサ間でのデータの転置
（入れ替え）を行う。すなわち、今度は既にステップｃ
で述べたように、２次変換結果データの直方体を図８の
ようにＹ軸に垂直にスライスし、こうしてできる各面を
一つのプロセッサに割り当てる（ステップ５８）。既に
ステップｃで述べたように、この割り当てに従い、各プ
ロセッサが自分以外の全プロセッサとの間でそれぞれの
プロセッサが生成した２次変換結果データの交換を行
う。具体的には、この転置時には、各プロセッサは、上
記第２の作業配列に、そのプロセッサに割り当てられた
座標ｋ_yの値を有するＹ軸に垂直な平面に属すべき、ｋ_y
が特定値で、k_x，ｊ_zが種々の値を持つ２次変換結果デ
ータｃ''_kx,ky,jz（但し、k_x＝０〜ＮＸ−１，ｋ_y＝特
定値，ｊ_z＝０〜ＮＺ−１）を受信するように、プロセ
ッサ間で２次変換結果データｃ''_kx,ky,jzを交換する。
転置の終了後、Ｚ方向の変換を行う（ステップ５９）。
すなわち、各プロセッサは、既にステップｄで記載した
ように、そのプロセッサに新たに割り当てられた２次変
換結果データに対して、式７により指定される変換を実
行し、最終的な３次元のフーリエ変換結果データを生成
する。転置により各ＸＺ平面は１台のプロセッサに担当
されているから、この変換処理も通信なしに各プロセッ
サで独立に行える。具体的には、各プロセッサでは、Ｚ
方向の変換は、ｋ_y座標が特定の値を有し、ｋ_x座標とＺ
座標とがいろいろの値を有する一群の２次変換結果デー
タｃ''_kx,ky,jzに対して、式７により行なわれる。この
とき、ｋ_y座標が特定の値を有し、ｋ_x座標が異なる一群
の入力データに対してＺ方向の変換が高速フーリエ変換
アルゴリズム（ＦＦＴ）を用いて実行される。この変換
は、上記第２の３次元作業配列を使用して実行される。
この変換の実行にあっては、ｋ_z座標が異なる一群のデ
ータに対して、プロセッサ内のベクトル演算器が使用さ
れ、計算がパイプライン的に実行される。その結果得ら
れる最終フーリエ変換結果データｃ_kx,ky,kz（但し、k_x
＝０〜ＮＸ−１，ｋ_y＝特定値，ｋ_z＝０〜ＮＺ−１）
は、第３の３次元作業配列の、これらの座標値ｋ_x，
ｋ_y，ｋ_zに対応するインデックスのところに格納され
る。したがって、図８のように、一つのプロセッサに一
つの座標ｋ_yを有す一つの平面が割り当てられた場合、
各プロセッサでは、最終フーリエ変換結果データｃ
_kx,ky,kzは、上記第３の３次元作業配列の、特定の座標
ｋ_yを有する一つの平面に格納される。もし、図８にお
いて、ｋ_y軸に垂直な複数の平面がそのプロセッサに割
り当てられているときには、それぞれの平面に対応す
る、上記第３の３次元作業配列内の、ｋ_y軸に垂直な複
数の平面のそれぞれに対応する最終フーリエ変換結果デ
ータｃ_kx,ky,kzが格納される。Ｚ方向の変換が終了する
と、一次元の変換対象データｆ_jのフーリエ変換が終了
し、重ね合わせの係数ｃ_kが求まる。ｃ_kの並び方は、原
点からまずＹ方向に、Ｙ方向にＮＹ個行ったら次はＸ座
標が１だけ増え、ＸＹ平面上にＮＸ＊ＮＹ個のデータが
並んだら次はＺ座標が１だけ増える、という順で並ぶ
（図３）。このデータの並び方と図８のデータの分割形
式とを照らし合わせることにより、本実施の形態では、
出力データｃ_kもサイクリック分割になっていることが
わかる。上記第３の３次元作業配列内でも、最終フーリ
エ変換結果データｃ_kx,ky,kzはこの並びに対応する並び
を有する。ライブラリはこのデータｃ_kx,ky,kzを一次元
座標ｋの順に並び替えて一次元配列Ｆに格納し（ステッ
プ６１）、リターンする（ステップ６２）。本ライブラ
リでは、従来法で必要であった変換後のデータ分割形式
の変更が不要となり、通信の削減により従来法を上回る
並列化効率を得ることが可能となり、フーリエ変換時間
を低減できる。なお、ＮＸ，ＮＹ，ＮＺの決め方として
は、プロセッサ台数をｐとすると、Ｙ方向、Ｘ方向の変
換でＺ方向に垂直な面でデータを分割することから、Ｎ
Ｚ≧ｐが成り立つ必要がある。また、Ｚ方向の変換では
Ｙ方向に垂直な面でデータを分割することから、ＮＹ≧
ｐも成り立つ必要がある。また、並列計算機２８の各プ
ロセッサ２９がベクトル演算器（図示せず）を備えると
仮定した。このような並列計算機では、このベクトル演
算器を効率的に使うためには、ベクトル化の対象となる
ループの長さ（すなわち、同じ演算を受けるデータ群
（ベクトルデータ）の要素数であり、ベクトル長とも言
われる）をできるだけ長く取る必要があることが知られ
ている。本アルゴリズムで式５から７を計算するときに
は、このベクトル演算器が使用される。ベクトル化の対
象となるループは、フーリエ変換にも並列化にも使わな
い座標軸の方向で複数のデータに対して同じ演算を実行
する計算であり、Ｙ方向、Ｘ方向、Ｚ方向の変換におい
て、それぞれＸ方向、Ｙ方向、Ｘ方向での演算となる。
したがって、ベクトル演算器の性能を引き出すには、Ｎ
Ｙ≧ｐ，ＮＺ≧ｐの２つの条件を満たしつつＮＸとＮＹ
をできるだけ大きく取るようにＮＸ，ＮＹ，ＮＺを決め
ることが望ましい。なお、並列計算機２８は、ベクトル
演算器を有すると仮定したいが、この演算器がフーリエ
変換において必要な全ての演算の一部の演算をパイプラ
イン的に実行できるものでもよい。さらに、並列計算機
２８がベクトル演算器を有しない並列計算機であって
も、ループ長を大きくすることが高速化に有効である場
合が多い。また、以上の動作の説明では、並列計算機２
８がメモリ２９と演算器（図示せず）の間に複数のベク
トルレジスタを有しないと仮定し、各方向の変換で利用
される配列はメモリ２９から直接演算器に読み出され、
あるいはその変換で生成される配列はメモリ２９に直接
演算器から書き込まれるかのように説明した。しかし、
複数のベクトルレジスタを有する並列計算機では、メモ
リ２９上の配列に対する演算あるいはその演算の結果得
られた配列のメモリ２９への格納は、これらのレジスタ
を介して実行させればよいことは当業者に明らかであ
る。本実施の形態では本ライブラリは逆フーリエ変換を
実行するためのプログラムを有しない。後で説明するよ
うに、シミュレーションプログラムが逆フーリエ変換を
必要とするときには、シミュレーションプログラムの方
で、逆フーリエ変換の対象のデータの複素共役データを
生成し、そのデータに対してフーリエ変換を本ライブラ
リに要求する。この複素共役データに対して得られたフ
ーリエ変換データの複素共役をシミュレーションプログ
ラムが生成する。しかし、本ライブラリにフーリエ変換
の機能を持たせることもできる。すなわち、本ライブラ
リの引数としてフーリエ変換か逆フーリエ変換かを指定
する引数を追加し、シミュレーションプログラムが逆フ
ーリエ変換を要求したときには、本ライブラリで、変換
対象データの複素共役を求め、これにフーリエ変換を上
記のようにして実行し、得られた結果データの複素共役
を求め、それを逆フーリエ変換結果データとしてシミュ
レーションプログラムに戻せばよい。 （５）シミュレーションプログラム 本実施の形態において使用するシミュレーションプログ
ラムの例として気象計算のための並列プログラムを図１
１に示す。気象計算は本来３次元の計算であるが、現在
は計算機能力の制約から２次元で行うことも多い。そこ
で本実施の形態では、２次元の気象予測対象とする領域
（これが計算対象領域となる）の場合を例にとってシミ
ュレーションプログラムを説明する。ユーザは、予め全
計算対象領域を複数の部分計算領域に区分し、それぞれ
をいずれか一つのプロセッサに割り当てる。さらに、各
部分計算領域のサイズＮ１，Ｎ２、フーリエ変換で用い
るＮＸ，ＮＹ，ＮＺなどのパラメータを指定する。ユー
ザが本プログラムを並列計算機２８で使用するときに
は、まず、ワークステーション１が、並列計算機２８内
の特定のプロセッサ（たとえばプロセッサ０）と交信し
て、このプログラムと上記ユーザ指定の情報と、空気の
熱伝導率などの計算に使用する物質定数と、全計算対象
領域に内の観測によって得られた温度・風速・圧力など
の初期値データとを、その特定のプロセッサ０を介して
外部記憶装置３１に記憶する。その後、その特定のプロ
セッサ０が、各プロセッサに本プログラムをロードし、
全プロセッサで本プログラムを起動する。本プログラム
は、並列計算機２８内の全プロセッサで全く同じように
して並行して実行される。本プログラムでは、起動され
ると、まず計算領域のサイズＮ１，Ｎ２、フーリエ変換
で用いるＮＸ，ＮＹ，ＮＺ、空気の熱伝導率などの物質
定数などのパラメータと、観測によって得られた温度・
風速・圧力などの初期値データとを外部記憶装置３１か
ら入力する（ステップ３２）。本プログラムは、それが
ロードされたプロセッサがどの部分計算領域に関する計
算を実行するかを判断するようにプログラムされている
と仮定する。このステップでは、各プロセッサは、プロ
セッサに依らないで使用される上記パラメータを入力す
るとともに、外部記憶装置３１に記憶された全計算対象
領域に対する初期値データの内、そのプロセッサに割り
当てられた部分計算領域に関する初期値データを選択し
て外部記憶装置３１から入力する。なお、上記（３）の
「多次元フーリエ変換への応用」の項で述べたように、
本発明の方式による２次元高速フーリエ変換では、入力
データが第１の座標方向にサイクリック分割されている
必要がある。すなわち、サイズＮ１×Ｎ２のメッシュ上
で定義されたある物理量Ａ_j1,j2（ただしｊ₁＝０, １，
... ，Ｎ１−１，ｊ₂＝０, １， ... ，Ｎ２−１）の
うち，要素Ａ_m*NPU+i,j2（ｍ=０,１,...,Ｎ１／ＮＰＵ
−１，ｉ＝０，１，...，ＮＰＵ−１， ｊ₂＝０, １，
... ，Ｎ２−１)は第ｉ番目のプロセッサに割り当てら
れている必要がある。そこで，本実施の形態のシミュレ
ーションプログラムでも，２次元高速フーリエ変換での
入力形式に合わせて，このように第１の座標方向をサイ
クリック分割することによって得られる部分計算領域を
用いる。その後計算に必要な前処理を行う（ステップ３
３）。ここで前処理とは、観測によって得られた温度・
風速・圧力などのデータに対して補間を行い、計算に必
要なメッシュポイントでの温度・風速・圧力などのデー
タを得ることである。これらの処理が終わった後、以下
に説明する繰り返しループにより各時間ステップでの温
度・風速・圧力などの量を順々に求めていく。基礎とな
る方程式は、以下に示す風速に対する運動方程式、質量
保存の式、温度変化を表す式の３本である。

【数１２】 ｄｕ／ｄｔ＝−２Ω×ｕ−（１／ρ）∇ｐ＋Ｆ_u ，．．．（１２）

【数１３】 ｄρ／ｄｔ＝−ρ∇・ｕ ．．．（１３）

【数１４】 ｄＴ／ｄｔ＝−κ∇²Ｔ＋ｕ・∇Ｔ ．．．（１４） ここで、ｕは風速、ｐは圧力、Ｔは温度を表し、Ωはコ
リオリ力と呼ばれる地球の自転による力、Ｆ_uはそれ以
外の外力、ρは空気の密度、κは空気の熱伝導率を表
す。これらの式から次の時刻でのデータの値を求めるに
は、まずフーリエ変換により格子点上の温度Ｔ、圧力ｐ
および風速ｕをそれぞれ波数空間でのデータに変換す
る。そのために、それぞれの物理量についてのデータに
ついて２次元高速フーリエ変換ライブラリＦＦＴ２Ｄを
順次コールする（ステップ３４）。ライブラリＦＦＴ２
Ｄのコール時には、既に述べた引数を指定する。波数空
間でそれぞれの物理量のデータを微分する（ステップ３
５）。すなわち、ライブラリＦＦＴ２Ｄから与えられ
る、各物理量に関するフーリエ変換係数データを波数空
間で微分し、その物理量についての、波数空間の格子点
上での温度勾配∇Ｔ、２次微分∇²Ｔ、圧力勾配∇ｐ、
風速の発散∇・ｕ等の微分に関連するデータを求める。
各物理量についての上記微分に関連するデータを逆フー
リエ変換して、実空間の格子点上での温度勾配∇Ｔ、２
次微分∇²Ｔ、圧力勾配∇ｐ、風速の発散∇・ｕ等の微
分に関連するデータを求める（ステップ３６）。逆フー
リエ変換するには、すでに述べた式８、９を使用する。
すなわち、各物理量についての上記微分後のデータの複
素共役なデータを生成し、ライブラリＦＦＴ２Ｄをコー
ルしてこの複素共役なデータに対するフーリエ変換を要
求する。さらに、得られたフーリエ変換結果データの複
素共役なデータを生成し、この生成された複素共役なフ
ーリエ変換結果データを逆フーリエ変換結果データとし
て使用する。この後、上記逆フーリエ変換で得られた微
分に関連するデータを、式１２−１４の右辺に代入し、
風速ｕ、空気密度ρ、温度Ｔのそれぞれに関する時間微
分を決定し、得られたそれらの時間微分を用いて、次の
時間ステップでの温度・風速・圧力を求める（ステップ
３７）。なお、ステップ３４，３５，３６で、フーリエ
変換により実空間上の格子点のデータを波数空間上のデ
ータに変換してそこで微分関連データを求め、得られた
微分関連データを逆フーリエ変換して実空間に関する微
分関連データを得るのは、その方が微分が精度良く計算
できるからであり、本シミュレーションプログラムで
は、この計算部分で２次元フーリエ変換ライブラリＦＦ
Ｔ２Ｄを用いる。上記のループでは、各時間ステップ毎
に、求める時刻までの計算が終了したかどうかを判定し
（ステップ３８）、終了したら、後処理を行い（ステッ
プ３９）、予測結果データとして出力する（ステップ４
０）。後処理では、主に計算を行うメッシュポイントと
予測結果データが必要な点とがずれている場合に、計算
結果データを補間して必要な点での予測値を計算するな
どの処理を行う。出力処理４０では、各プロセッサは、
生成したデータを外部記憶装置３１に書き戻し、上記特
定のプロセッサはシミュレーションプログラムの実行終
了時に、このデータをワークステーション１に一つの結
果データとして転送する。なお、以上ではシミュレーシ
ョンプログラムは、並列計算機２８内の全プロセッサで
全く同じようにして並行して実行されると仮定した。さ
らに、それがロードされたプロセッサがどの部分計算領
域に関する計算を実行するかを判断するようにプログラ
ムされていると仮定する。しかし、本発明に依るフーリ
エ変換を利用するプログラムはこのようなプログラムに
限定されないことはいうまでもない。上記フーリエ変換
ライブラリＦＦＴ２Ｄを使用するには、それぞれのライ
ブラリが要求する上記複数の引数を指定することが必要
であり、それらの引数の生成あるいは獲得は他の方法で
も良い。たとえば、本プログラムは、並列計算機２８内
の特定の一つのプロセッサで実行される単独処理部分と
全プロセッサで並行して実行される並列処理部分とから
構成されてもよい。たとえば、本プログラムがいずれか
のプロセッサで起動されたときに、そのプロセッサが上
記特定の一つのプロセッサであるときにその単独処理部
分が実行され、そうでないときには上記並列処理部分の
みが実行される。上記単独処理部分では、各プロセッサ
が担当する部分計算領域を判断し、その結果を他のプロ
セッサに通知するように構成できる。上記の例では気象
予測計算を行う場合を例にとって説明したが、本発明の
手法は、これ以外の応用例についても、並列計算機上で
高速フーリエ変換を用いてシミュレーションを行う場合
に適用できることは明らかである。一次元フーリエ変換
ライブラリＦＦＴ１Ｄについても全く同様である。 ＜発明の実施の形態２＞上記の実施の形態１では、フー
リエ変換ライブラリは、図９の方式２を採用した。しか
し、本実施の形態では、フーリエ変換ライブラリは、図
９の方式３を採用する。この方式３では、Ｙ方向の変換
を行った後で転置を行い、その後Ｘ方向とＺ方向の変換
を行う。この方式２では、フーリエ変換は具体的には以
下のようにして実行される。 （ステップａ’）Ｙ方向の変換 Ｙ方向の変換は、従来の転置アルゴリズムについて既に
述べたステップ１あるいはａ’の要領で実行される。既
に述べたごとく、フーリエ変換対象データのデータ分割
は、サイクリック分割である。 （ステップｂ’）データ転置 式６によるＸ方向の変換を行うには、３次元波数空間の
座標ｋ_yの特定の値と、３次元実空間の座標（j_x，j_z）
の全ての値に対して得られた１次変換結果データｃ'
_jx,ky,jzが必要である。方式３では、方式２と異なり、
Ｙ方向の変換で得られた１次変換結果データｃ'
_jx,ky,jzは、Ｙ軸に垂直な複数の平面で分割される。こ
のことは、図８に示すように、３次元波数空間を、その
ｋ_y軸に垂直な複数の平面で分割して、分割後の部分空
間（上記複数の平面）の各々を一つのプロセッサに割り
当てることを意味する。すなわち、各プロセッサへの座
標ｋ_yの特定の値を割り当てる。具体的には、ｋ_y＝０，
１，２，３，，，の１変換結果データｃ'_jx,ky,jzを順
次プロセッサ０，１，２，，，７に割り当てる。この割
り当てに従い、後にＸ方向の変換を実行するに必要なデ
ータをプロセッサ間で転送する処理がここのステップ
ｂ’で実行される。各プロセッサが、そのプロセッサに
割り当てられた３次元波数空間の座標ｋ_yの特定の値
と、３次元実空間の座標（j_x，j_z）の全ての値とに対し
て得られた全ての１次変換結果データｃ'_jx,ky,jzを使
用してＸ方向の変換を実行できるように、全プロセッサ
間で１次変換結果データの転送が行われる。すなわち、
各プロセッサは、３次元実空間の座標（j_x，j_z）の全て
値と、３次元波数空間の座標ｋ_yの特定の値との組に対
して得られた１次変換結果データｃ'_jx,ky,jz（但し、
ｊ_x＝０〜（ＮＸ―１）、ｊ_z＝０〜（ＮＺ―１）、ｋ_y
＝特定値）の内、自プロセッサが生成しなかったデータ
を他のプロセッサから受信するように、全プロセッサの
間で２次変換結果データを転送する。 （ステップｃ’）Ｘ方向の変換 次に、Ｘ方向の変換を実行する。各プロセッサは、式６
により、他のプロセッサとの通信をしないで、そのプロ
セッサに割り当てられた３次元波数空間の座標ｋ_yの特
定の値と、３次元波数空間の座標ｋ_xの全ての値と、３
次元実空間の座標ｊ_zの全ての値に関連するＮＸ＊ＮＺ
個の２次変換結果データ（ｃ''_kx,ky,jz）（ｋ_x＝０〜
ＮＸ−１、ｊ_z＝０〜ＮＺ−１）を得る。 （ステップｄ’）Ｚ方向の変換 各プロセッサは、そのプロセッサに割り当てられた３次
元波数空間の座標ｋ_yの特定の値と、３次元波数空間の
他の二つの座標ｋ_x、ｋ_zの全ての値を有する最終的なフ
ーリエ変換結果データｃ_kx,ky,kz（但し、ｋ_x＝０〜Ｎ
Ｘ―１、ｋ_y＝特定値、ｋ_z＝０〜ＮＺ―１）を計算す
る。各プロセッサはこの計算を互いに並列に実行でき
る。この結果、座標ｋ_y＝０，１，２，３，，，に対応
するフーリエ変換係数ｃ₀，ｃ₁，ｃ₂，，，が、順次プ
ロセッサ０，１，２，，，で生成され、全フーリエ変換
結果データｃ_kx,ky,kzは、プロセッサ間でサイクリック
に分割されていることが分かる。本方式は、実施の形態
１で使用した方式２に比べ、複素指数関数テーブルを格
納する配列ＢＴの容量の点で有利となる。実際、式６に
より、Ｘ方向の変換における複素指数関数の値はＸ方
向、Ｙ方向のインデックスのみに依存し、Ｚ方向のイン
デックスには依存しない。したがって、実施の形態１の
ようにＸ方向の変換においてＺ軸に垂直な分割を採用し
た場合には、各プロセッサが同じテーブルを重複して持
つことになる。それに対して本方式では、分割をＹ軸に
垂直な面で行うので、各プロセッサが自分の計算に必要
なテーブルの一部分のみを持つことになり、重複はな
い。これにより、本方式ではＸ方向の変換に必要なテー
ブルの大きさが１／（プロセッサ台数）に削減できる。
本実施の形態では、ベクトル化の対象となるループはＹ
方向、Ｘ方向、Ｚ方向の変換において、それぞれＸ方
向、Ｚ方向、Ｘ方向となるので、ベクトル並列計算機の
性能を引き出すには、ＮＹ≧ｐ，ＮＺ≧ｐの２つの条件
を満たしつつＮＸとＮＺをできるだけ大きく取るのがよ
い。 ＜発明の実施の形態３＞本実施の形態において対象とな
る並列計算機は、実施の形態１で説明した図１の並列計
算機システムとほぼ同様であるが、各プロセッサは同一
のベクトル演算器を内蔵し、かつ、外部記憶装置中３２
に、そのベクトル演算器の性能に関するデータベースを
持つ。ベクトル演算器の演算性能とはたとえば単位時間
あたりに実行可能な演算回数である。ベクトル演算器性
能データベース中には、ベクトル演算器の性能データが
ループ長Ｌの関数（ｇ（Ｌ））の形で格納されている。
実施の形態１では、フーリエ変換のためのパラメータＮ
Ｘ，ＮＹ，ＮＺはプログラムへの入力として決定してい
たが、これを並列計算機２８の各々のプロセッサを構成
するベクトル演算器の特性に応じて最適化することによ
り、さらに効率的な計算が可能となる。一般に、ベクト
ル演算器の演算性能は、ループ長Ｌの関数ｇ（Ｌ）であ
る。ｇ（Ｌ）は通常、Ｌに対して単調に増加する関数で
ある。いま、実施の形態１のフーリエ変換方式でＹ方向
の変換を計算するステップの演算量を考えると、ＮＹ次
のフーリエ変換を一回行うための演算量はＮＹｌｏｇＮ
Ｙであり、これをＮＸ＊ＮＺ組計算するから、全体での
演算量はＮＸ＊ＮＹ＊ＮＺｌｏｇＮＹ＝ＮｌｏｇＮＹで
ある。同様にして、ＮＸ方向、ＮＺ方向での演算量は、
それぞれＮｌｏｇＮＸ、ＮｌｏｇＮＺである。一方、そ
れぞれの演算におけるベクトル化のループ長は、実施の
形態１で述べたようにＮＸ，ＮＹ，ＮＸであるから、ベ
クトル演算器の演算性能はそれぞれｇ（ＮＸ），ｇ（Ｎ
Ｙ），ｇ（ＮＸ）となる。演算時間ｔは演算量を演算性
能で割ることによって得られ、合計で t=NlogNY/g(NX)+NlogNX/g(NY)+NlogNZ/g(NX) となる。したがって、プロセッサ台数をｐとするとき、
ＮＸ≧ｐ，ＮＺ≧ｐという条件の下でｔを最小化するよ
うにＮＸ，ＮＹ，ＮＺを決めることにより、ベクトル演
算器の性能を最大限に発揮できる高速フーリエ変換が実
現できる。本実施の形態でのライブラリのフローチャー
トを図１２に示す。処理は、ＮＸ，ＮＹ，ＮＺの決定部
分（ステップ４３）を除いては、実施の形態１（図１
０）と同様である。ステップ４３では、上記ベクトル演
算器性能データベースを用いて、ＮＸ≧ｐ，ＮＺ≧ｐと
いう条件の下で上記演算時間ｔを最小化するようにＮ
Ｘ，ＮＹ，ＮＺを決める。その後の処理は、実施の形態
１と同様である。このライブラリへのコール文ではシミ
ュレーションプログラムはこれらのパラメータＮＸ，Ｎ
Ｙ，ＮＺを指定する必要はない。本実施の形態の方式に
よれば、ユーザは自分でＮＸ，ＮＹ，ＮＺを計算するこ
となく、ベクトル演算器を内蔵する並列計算機の性能を
最大限に引き出すことが可能となる。なお、ＮとＮＰＵ
とが一般の整数の場合には、ＮＸ，ＮＹ，ＮＺを変える
ことにより、入力データのプロセッサへの分割形式も変
更する必要があるが、フーリエ変換でもっともよく利用
される、ＮおよびＮＰＵが共に２のべき乗の場合には、
ＮＸ，ＮＹ，ＮＺを変えても、分割形式を変更する必要
がない場合がある。実際、２つの組（ＮＸ，ＮＹ，Ｎ
Ｚ）＝（ＮＸ１，ＮＹ１，ＮＺ１）、（ＮＸ２，ＮＹ
２，ＮＺ２）が共にＮＹ≧ＮＰＵ，ＮＺ≧ＮＰＵの２つ
の条件を満たしているとき、入力データｆ_jを図２に示
す順番で直方体状に並べ、これをＺ軸に垂直な面でスラ
イスして、各面をサイクリックにプロセッサ０，
１，．．．，ＮＰＵ−１に割り当てたとする。すると、
（ＮＸ，ＮＹ，ＮＺ）＝（ＮＸ１，ＮＹ１，ＮＺ１）の
場合は、ｆ_jの属する面は上からＭＯＤ（ｆ_j，ＮＺ１）
＋１番目であり、この面を担当するプロセッサの番号は MOD(MOD(ｆ_j,NZ1),NPU) である。一方、（ＮＸ，ＮＹ，ＮＺ）＝（ＮＸ２，ＮＹ
２，ＮＺ２）の場合も同様にして、ｆ_jを担当するプロ
セッサの番号は MOD(MOD(ｆ_j,NZ2),NPU) となる。ところが、いまＮＺ１≧ＮＰＵ，ＮＺ２≧ＮＰ
Ｕであり、ＮＺ１，ＮＺ２，ＮＰＵはすべて２のべき乗
であるから、ＮＺ１，ＮＺ２は共にＮＰＵの倍数であ
る、したがって、 MOD(MOD(ｆ_j,NZ1),NPU) =MOD(MOD(ｆ_j,NZ2),NPU) =MOD(ｆ_j,NPU)すなわち、ｆ_jを担当するプロセッサの番
号は、どちらの場合も同じである。以上の考察より、Ｎ
およびＮＰＵが共に２のべき乗で、ＮＹ≧ＮＰＵ，ＮＺ
≧ＮＰＵの２つの条件が成り立っている限り、ＮＸ，Ｎ
Ｙ，ＮＺを変えても、入力データｆ_jのプロセッサへの
分割形式は変更する必要がないことがわかる。このこと
を利用し、分割形式を変えずに済む範囲でＮＸ，ＮＹ，
ＮＺの最適化を行えば、分割形式変更に伴う新たな通信
オーバーヘッドを生じることなく、ベクトル演算器を含
む並列計算機での処理速度、具体的には、フーリエ変換
速度を向上させることができる。 ＜発明の実施の形態４＞本発明による高速フーリエ変換
を用いてシミュレーションを行う他の例として、半導体
デバイス等における電子構造計算を説明する。電子構造
計算は、その結果を利用して半導体デバイスの設計、と
くにデバイス構造の決定に使用されている。電子構造計
算では、２次元または３次元のメッシュで定義された電
子の波動関数ｕ（ｒ）を、次のシュレディンガー方程式

【数１５】 ｄｕ（ｒ）／ｄｔ＝−（ｈ²／２ｍ）∇²ｕ（ｒ） ＋（Ｅ−Ｖ（ｒ））ｕ（ｒ） ．．． （１５） に従って計算することにより、半導体の性質を決定する
バンドギャップの大きさや、結晶の構造安定性などを求
める。ただし、上式で、ｈはプランク定数、ｍは電子の
質量、Ｅは対象とする波動関数のエネルギーレベル、Ｖ
は結晶中の原子や他の電子によるポテンシャルエネルギ
ーを表す。式１５の計算では、波動関数ｕ（ｒ）の２次
微分∇²ｕ（ｒ）が必要であるが、気象計算の例におい
て述べたのと同様な理由により、この部分はｕ（ｒ）を
フーリエ変換により波数空間に移してから計算し、結果
を逆フーリエ変換で再び実空間に戻す。したがって、並
列計算機上で電子構造計算を行う場合には、この部分で
本発明の高速フーリエ変換方法が適用できる。 ＜変形例＞本発明は、以上の実施の形態に限定されるの
ではなく、以下に例示する変更あるいは変形以外のいろ
いろの変更あるいは変形により実現可能である。 （１）本発明によるフーリエ変換方法は、シミュレーシ
ョンに限らず他の用途にも使用できるのは言うまでもな
い。たとえば、伝送される信号あるいは地震波等の波動
の解析に利用でき、解析の結果を用いて、信号伝送に関
係する装置、例えば伝送装置あるいは伝送線路の設計を
行うことができ、あるいは地震を利用した応用、例えば
資源開発等にも利用できる。 （２）以上の実施の形態では、フーリエ変換変換はその
ために用意されたフーリエ変換ライブラリにより実行さ
れた。しかし、本発明は、フーリエ変換を使用するアプ
リケーションプログラム自身にこのフーリエ変換手順を
実行するプログラムを組み込んでもよいことは明らかで
ある。このようなシミュレーションプログラムは。プロ
グラムを磁気記憶装置のようなプログラム記録媒体に記
憶して販売できる。 （３）本発明は、フーリエ変換対象データｆ_jが実デー
タであるときにも適用できる。その場合に、Ｙ方向等の
変換のときに、係数の計算においては、虚数部の計算を
省略することができる。 以上から明らかなように、本発明によれば、並列計算機
を使用してフーリエ変換を従来より高速に実行できる。
たとえば本出願人により開発された並列計算機ＳＲ２２
０１を用いて本発明の効果を評価した結果では以下の通
りである。３次元フーリエ変換を実行する場合、従来法
では、２５６×２５６×２５６のサイズのデータを２５
６台のプロセッサを用いて変換するのに、約０．２６秒
の時間が必要である。この内訳は、計算に０．１４秒、
途中でのデータの転置に０．０６秒、最後のデータの並
べ替えに０．０６秒の時間がかかる。実施の形態１ある
いは２に記載の方法によれば、計算と転置の時間は従来
法と同じであり、最後のデータの並べ替えが省略できる
ので、０．２０秒でフーリエ変換を行うことができ、約
２４％の高速化が達成できる。とくに、実施の形態１で
記載した気象計算を、３次元フーリエ変換を用いて行う
場合、気象計算では全計算時間の約５０％がフーリエ変
換で占められるため、約１２％の高速化が得られる。ま
た、実施の形態４で記載した電子構造計算を３次元フー
リエ変換を用いて行う場合、通常全実行時間の３０％程
度がフーリエ変換で占められるため、約７％の高速化が
達成できる。

【発明の効果】以上説明したように、本発明によれば、
並列計算機上でのフーリエ変換を高速に実行できる。

【図面の簡単な説明】

【図１】本発明の第１の実施の形態で使用する並列計算
機の概略構成図。

【図２】本発明の第１の実施の形態で使用する一次元変
換対象データの３次元データへの写像を説明する図。

【図３】本発明の第１の実施の形態で使用する一次元変
換結果データの３次元データへの写像を説明する図。

【図４】本発明の第１の実施の形態で使用する、プロセ
ッサへデータを割り当てる第１の方法を示す図。

【図５】従来技術で使用する、プロセッサへデータを割
り当てる他の方法を示す図。

【図６】本発明の第１の実施の形態で使用する、一次元
変換対象データのプロセッサ間データ分割形式を説明す
る図。

【図７】従来技術による、一次元変換結果データのプロ
セッサ間データ分割形式を説明する図。

【図８】本発明の第１の実施の形態で使用する、プロセ
ッサへデータを割り当てる第２の方法を示す図。

【図９】本発明に至る前に比較検討した、複数のフーリ
エ変換変換手順を示す図。

【図１０】本発明の実施の形態１で使用するフーリエ変
換ライブラリのフローチャート。

【図１１】本発明の実施の形態１で使用するシミュレー
ションプログラムのフローチャート。

【図１２】本発明の実施の形態３で使用するフーリエ変
換ライブラリのフローチャート。

Claims (28)

【特許請求の範囲】
1. 【請求項１】複数のプロセッサを有する計算機で実行す
   るためのフーリエ変換方法であって、 各プロセッサにより、第１の変換処理、第２の変換処
   理、第３の変換処理を順次かつ他のプロセッサと並行し
   て実行し、 上記複数のプロセッサの各々による、上記第１、第２の
   変換処理のいずれか一方の変換処理の実行後に、上記複
   数のプロセッサでのその一方の変換処理を実行した結果
   得られた一群の結果データを構成する複数の結果データ
   部分群が異なるプロセッサに割り当てられるように、上
   記一群の結果データを上記複数のプロセッサの間で交換
   するステップを有し、 上記第１から第３の変換処理は、一群の順序づけられた
   変換対象データに対する一群の順序づけられたフーリエ
   変換係数データを構成する複数のフーリエ変換係数デー
   タ部分群をそれぞれ異なるプロセッサにより生成するよ
   うに定められ、 各プロセッサには、上記一群の変換対象データを構成す
   る複数の変換対象データ部分群の一つの変換対象データ
   部分群がそのプロセッサに対して予め割り当てられ、 上記一群のフーリエ変換係数データのそれぞれを生成し
   たプロセッサの順序が、上記一群の変換対象データのそ
   れぞれが割り当てられたプロセッサの順序と同一となる
   ように、上記交換するステップで各プロセッサに割り当
   てられる結果データ部分群が定められているもの。
2. 【請求項２】各プロセッサは、そのプロセッサに対して
   予め割り当てられた一つの変換対象データ部分群に対し
   て、上記第１の変換処理を他のプロセッサと並行して実
   行し、第１の結果データの部分群を生成し、 各プロセッサは、上記交換するステップが上記第２の変
   換処理の実行後に実行されるときには、その各プロセッ
   サで実行した上記第１の変換処理の結果生成した上記一
   つの第１の結果データ部分群に対して、上記交換するス
   テップが上記第１の変換処理の実行後に実行されたとき
   には、その交換ステップでそのプロセッサに割り当てら
   れた一群の第１の部分結果データに対して、上記第２の
   変換処理を他のプロセッサと並行して実行して、一つの
   第２の結果データ部分群を生成し、 各プロセッサは、上記交換するステップが上記第１の変
   換処理の実行後に実行されたときには、その各プロセッ
   サで実行した上記第２の変換処理の結果生成した一つの
   第２の結果データ部分群に対して、上記交換するステッ
   プが上記第２の変換処理の実行後に実行されたときに
   は、その交換ステップでそのプロセッサに割り当てられ
   た一つの第２の結果データ部分群に対して、上記第３の
   変換処理を他のプロセッサと並行して実行して一つのフ
   ーリエ変換係数データ部分群を生成する請求項１記載の
   フーリエ変換方法。
3. 【請求項３】上記第１、第２、第３の変換処理は、それ
   ぞれ３次元データ空間の第１、第２、第３の座標軸に関
   する変換処理であり、 上記変換対象データ群の各々は、上記３次元データ空間
   の直方体形状に位置する格子点群の一つの座標をそれぞ
   れ有し、 上記複数の変換対象データ部分群は、上記３次元データ
   空間の第３の座標軸の座標値が同じであり、上記３次元
   データ空間の第１、第２の座標軸の座標値が異なる全て
   の変換対象データが同一の変換対象データ部分群に含ま
   れるように定められ、 上記フーリエ変換係数データ群の各々は、３次元係数空
   間の直方体形状に位置する格子点群の一つの座標をそれ
   ぞれ有し、 上記複数のフーリエ変換係数データ部分群は、上記３次
   元係数空間の第１の座標軸の座標値が同じであり、上記
   ３次元波数空間の第２、第３の座標軸の座標値が異なる
   全てのフーリエ変換係数データが同一のフーリエ変換係
   数データ部分群に含まれるように定められている請求項
   １記載のフーリエ変換方法。
4. 【請求項４】複数のプロセッサを有する計算機で実行す
   るためのフーリエ変換方法であって、 各プロセッサにより、３次元データ空間の第１の座標軸
   に関する第１の変換処理、上記３次元データ空間の第２
   の座標軸に関する第２の変換処理、上記３次元データ空
   間の第３の座標軸に関する第３の変換処理を順次実行
   し、 上記複数のプロセッサの各々による、上記第１、第２の
   変換処理のいずれか一方の変換処理の実行後に、上記複
   数のプロセッサでのその一方の変換処理を実行した結果
   得られた一群の結果データを構成する複数の結果データ
   部分群がそれぞれ異なるプロセッサに割り当てられるよ
   うに、上記一群の結果データを上記複数のプロセッサの
   間で交換するステップを有し、 上記第１から第３の変換処理は、一群の順序づけられた
   変換対象データに対する一群の順序づけられたフーリエ
   変換係数データを構成する複数のフーリエ変換係数デー
   タ部分群をそれぞれ異なるプロセッサにより生成するよ
   うに定められ、 上記変換対象データ群の各々は、上記３次元データ空間
   の直方体形状に位置する格子点群の一つの座標をそれぞ
   れ有し、 上記複数の変換対象データ部分群は、上記３次元データ
   空間の第３の座標軸の座標値が同じであり、上記３次元
   データ空間の第１、第２の座標軸の座標値が異なる全て
   の変換対象データが同一の変換対象データ部分群に含ま
   れるように定められ、 上記フーリエ変換係数データ群の各々は、３次元係数空
   間の直方体形状に位置する格子点群の一つの座標をそれ
   ぞれ有し、 上記複数のフーリエ変換係数データ部分群は、上記３次
   元係数空間の第１の座標軸の座標値が同じであり、上記
   ３次元波数空間の第２、第３の座標軸の座標値が異なる
   全てのフーリエ変換係数データが同一のフーリエ変換係
   数データ部分群に含まれるように定められ、 各プロセッサには、上記一群の変換対象データを構成す
   る複数の変換対象データ部分群の一つの変換対象データ
   部分群が予め割り当てられ、 上記一群のフーリエ変換係数データのそれぞれを生成し
   たプロセッサの順序が、上記一群の変換対象データのそ
   れぞれが割り当てられたプロセッサの順序と同一となる
   ように、上記交換するステップで各プロセッサに割り当
   てられる結果データ部分群が定められているもの。
5. 【請求項５】各プロセッサは、そのプロセッサに対して
   予め割り当てられた一つの変換対象データ部分群に対し
   て、上記第１の変換処理を他のプロセッサと並行して実
   行して、第１の結果データ部分群を生成し、 各プロセッサは、上記交換するステップが上記第２の変
   換処理の実行後に実行されるときには、その各プロセッ
   サで実行した上記第１の変換処理の結果生成した第１の
   結果データ部分群に対して、上記交換するステップが上
   記第１の変換処理の実行後に実行されたときには、その
   交換ステップでそのプロセッサに割り当てられた一つの
   第１の結果データ部分群に対して、上記第２の変換処理
   を他のプロセッサと並行して実行して、一つの第２の結
   果データ部分群を生成し、 各プロセッサは、上記交換するステップが上記第１の変
   換処理の実行後に実行されるときには、その各プロセッ
   サで実行した上記第２の変換処理の結果生成した一つの
   第２の結果データ部分群に対して、上記交換するステッ
   プが上記第２の変換処理の実行後に実行されたときに
   は、その交換ステップでそのプロセッサに割り当てられ
   た一つの第２の結果データ部分群に対して、上記第３の
   変換処理を他のプロセッサと並行して実行して一つのフ
   ーリエ変換係数データ部分群を生成する請求項４に記載
   のフーリエ変換方法。
6. 【請求項６】上記変換対象データ群が上記３次元データ
   空間に直方体形状に位置する格子点群に上記３次元空間
   に第３の座標軸、第２の座標軸、第１の座標軸の順に順
   次割り当てられ、 上記第１から第３の変換処理は、上記複数のフーリエ変
   換係数データが、３次元係数空間に直方体形状に位置す
   る格子点群に、当該３次元係数空間の第１、第２、第３
   の座標軸の順序で割り当てられるように定められている
   請求項４記載のフーリエ変換方法。
7. 【請求項７】各プロセッサが上記第１の変換処理により
   生成する上記一つの第１の結果データ部分群は、上記３
   次元データ空間の第３の座標軸の座標値が所定の同じ値
   であり、上記３次元データ空間の第２の座標軸の座標値
   と上記３次元係数空間の第１の座標軸の座標値が異なる
   値を有する全ての複数の第１の結果データを含み、 上記交換ステップが上記第１の変換処理が上記複数のプ
   ロセッサにより実行された後に実行され、 上記複数のプロセッサは、この交換ステップで、上記３
   次元係数空間の第１の座標軸の座標値が所定の同じ値で
   あり、上記３次元データ空間の第２、第３の座標軸の座
   標値が異なる値を有する全ての複数の第１の結果データ
   を含む第１の結果データ部分群が同一のプロセッサに割
   り当てられるように、上記複数のプロセッサが生成した
   一群の第１の結果データを上記複数のプロセッサの間で
   交換し、 各プロセッサが上記第２の変換処理により生成する上記
   一つの第２の結果データ部分群は、上記３次元係数空間
   の第１の座標軸の座標値が所定の同じ値であり、上記３
   次元波数空間の第２の座標軸の座標値と上記３次元デー
   タ空間の第３の座標軸の座標値が異なる値を有する全て
   の複数の第２の結果データを含み、 各プロセッサが上記第３の変換処理により生成する上記
   一つのフーリエ変換係数部分群は、上記３次元係数空間
   の第１の座標軸の座標値が所定の値であり、上記３次元
   波数空間の第２、第３の座標軸の座標値が異なる値を有
   する全ての複数のフーリエ変換係数を含む請求項４記載
   のフーリエ変換方法。
8. 【請求項８】各プロセッサが上記第１の変換処理により
   生成する上記一つの第１の結果データ部分群は、上記３
   次元データ空間の第３の座標軸の座標値が所定の同じ値
   であり、上記３次元データ空間の第２の座標軸の座標値
   と上記３次元係数空間の第１の座標軸の座標値が異なる
   値を有する全ての複数の第１の結果データを含み、 上記交換ステップが上記第２の変換処理が上記複数のプ
   ロセッサにより実行された後に実行され、 各プロセッサが上記第２の変換処理により生成する上記
   一つの第２の結果データ部分群は、上記３次元データ空
   間の第３の座標軸の座標値が所定の同じ値であり、上記
   ３次元係数空間の第１、第２の座標軸の座標値が異なる
   値を有する全ての複数の第２の結果データを含み、 上記複数のプロセッサは、上記交換ステップにより、上
   記３次元係数空間の第１の座標軸の座標値が所定の同じ
   値であり、上記３次元係数空間の第１の座標軸の座標値
   と上記３次元データ空間の第３の座標軸の座標値が異な
   る値を有する全ての複数の第１の結果データを含む第１
   の結果データ部分群が同一のプロセッサに割り当てられ
   るように、上記複数のプロセッサが生成した一群の第１
   の結果データを上記複数のプロセッサの間で交換し、 各プロセッサが上記第３の変換処理により生成する上記
   一つのフーリエ変換係数部分群は、上記３次元係数空間
   の第１の座標軸の座標値が所定の値であり、上記３次元
   係数空間の第２、第３の座標軸の座標値が異なる値を有
   する全ての複数のフーリエ変換係数を含む請求項４記載
   のフーリエ変換方法。
9. 【請求項９】複数のプロセッサを有する計算機で実行す
   るためのフーリエ変換方法であって、 各プロセッサにより、３次元空間の第１、第２、第３の
   座標軸の座標にそれぞれ関する第１、第２、第３の変換
   処理を順次かつ他のプロセッサと並行して実行し、 各プロセッサが上記第１、第２の変換処理のいずれか一
   方を実行した後に、その一方の変換処理の結果それぞれ
   のプロセッサで得られた複数の結果データを上記複数の
   プロセッサの間で交換するステップを有し、 ここで、一群の順序づけられた変換対象データが上記３
   次元空間に直方体の形に並べられ、 上記第１から第３の変換処理は、上記一群の変換対象デ
   ータに対する一群の順序づけられた３次元空間の座標を
   有する複数のフーリエ変換係数データを生成するように
   定められ、 上記複数の変換対象データが構成する上記直方体を分割
   する上記３次元空間の上記第１の座標軸に垂直な複数の
   面の各々に含まれる複数の変換対象データが同一のプロ
   セッサに割り当てられ、 上記交換ステップは、上記一方の変換処理の結果得られ
   た上記複数の結果データが構成する３次元空間の直方体
   を、その３次元空間の第１の座標軸に垂直な複数の面で
   分割し直して、各面に属する複数の結果データを同一の
   プロセッサに割り当てるように、上記一方の変換処理の
   結果得られた上記複数の結果データを上記複数のプロセ
   ッサ間で交換するステップを有するもの。
10. 【請求項１０】上記一群の順序づけられた変換対象デー
    タを上記３次元空間に直方体の形に並べられる順序は、
    第３の座標軸、第２の座標軸、第１の座標軸の順であ
    り、 上記第１から第３の変換処理は、上記複数のフーリエ変
    換係数データが３次元空間に第１、第２、第３の座標軸
    の順序で並べられるように定められている請求項９記載
    のフーリエ変換方法。
11. 【請求項１１】各プロセッサがパイプライン演算器を含
    み、その演算器での演算の対象とするループ長がＬのと
    きのその各プロセッサの演算性能を求めるための性能デ
    ータを上記複数のプロセッサに共通に記憶し、 その性能データを用いて、上記直方体の上記第１、第
    ２、第３の座標軸方向の長さを決定し、 その決定された上記第１、第２、第３の座標軸方向の長
    さを有する直方体に、上記順序づけられた複数の変換対
    象データを並べるステップをさらに有する請求項９記載
    のフーリエ変換方法。
12. 【請求項１２】計算機により読みとり可能なプログラム
    記録媒体であって、複数のプロセッサを有する計算機で
    フーリエ変換を実行するためのプログラムを記憶し、そ
    のプログラムは、 各プロセッサにより、第１の変換処理、第２の変換処
    理、第３の変換処理を順次かつ他のプロセッサと並行し
    て実行し、 上記複数のプロセッサの各々による、上記第１、第２の
    変換処理のいずれか一方の変換処理の実行後に、上記複
    数のプロセッサでのその一方の変換処理を実行した結果
    得られた一群の結果データを構成する複数の結果データ
    部分群が異なるプロセッサに割り当てられるように、上
    記一群の結果データを上記複数のプロセッサの間で交換
    するステップを実行するようにプログラムされ、 上記第１から第３の変換処理は、一群の順序づけられた
    変換対象データに対する一群の順序づけられたフーリエ
    変換係数データを構成する複数のフーリエ変換係数デー
    タ部分群をそれぞれ異なるプロセッサにより生成するよ
    うに定められ、 各プロセッサには、上記一群の変換対象データを構成す
    る複数の変換対象データ部分群の一つの変換対象データ
    部分群がそのプロセッサに対して予め割り当てられ、 上記一群のフーリエ変換係数データのそれぞれを生成し
    たプロセッサの順序が、上記一群の変換対象データのそ
    れぞれが割り当てられたプロセッサの順序と同一となる
    ように、上記交換するステップで各プロセッサに割り当
    てられる結果データ部分群が定められているもの。
13. 【請求項１３】各プロセッサは、そのプロセッサに対し
    て予め割り当てられた一つの変換対象データ部分群に対
    して、上記第１の変換処理を他のプロセッサと並行して
    実行し、第１の結果データの部分群を生成し、 各プロセッサは、上記交換するステップが上記第２の変
    換処理の実行後に実行されるときには、その各プロセッ
    サで実行した上記第１の変換処理の結果生成した上記一
    つの第１の結果データ部分群に対して、上記交換するス
    テップが上記第１の変換処理の実行後に実行されたとき
    には、その交換ステップでそのプロセッサに割り当てら
    れた一群の第１の部分結果データに対して、上記第２の
    変換処理を他のプロセッサと並行して実行して、一つの
    第２の結果データ部分群を生成し、 各プロセッサは、上記交換するステップが上記第１の変
    換処理の実行後に実行されたときには、その各プロセッ
    サで実行した上記第２の変換処理の結果生成した一つの
    第２の結果データ部分群に対して、上記交換するステッ
    プが上記第２の変換処理の実行後に実行されるときに
    は、その交換ステップでそのプロセッサに割り当てられ
    た一つの第２の結果データ部分群に対して、上記第３の
    変換処理を他のプロセッサと並行して実行して一つのフ
    ーリエ変換係数データ部分群を生成する請求項１２記載
    のプログラム記録媒体。
14. 【請求項１４】上記第１、第２、第３の変換処理は、そ
    れぞれ３次元データ空間の第１、第２、第３の座標軸に
    関する変換処理であり、 上記変換対象データ群の各々は、上記３次元データ空間
    の直方体形状に位置する格子点群の一つの座標をそれぞ
    れ有し、 上記複数の変換対象データ部分群は、上記３次元データ
    空間の第３の座標軸の座標値が同じであり、上記３次元
    データ空間の第１、第２の座標軸の座標値が異なる全て
    の変換対象データが同一の変換対象データ部分群に含ま
    れるように定められ、 上記フーリエ変換係数データ群の各々は、３次元係数空
    間の直方体形状に位置する格子点群の一つの座標をそれ
    ぞれ有し、 上記複数のフーリエ変換係数データ部分群は、上記３次
    元係数空間の第１の座標軸の座標値が同じであり、上記
    ３次元波数空間の第２、第３の座標軸の座標値が異なる
    全てのフーリエ変換係数データが同一のフーリエ変換係
    数データ部分群に含まれるように定められている請求項
    １３記載のプログラム記録媒体。
15. 【請求項１５】計算機により読みとり可能なプログラム
    記録媒体であって、複数のプロセッサを有する計算機で
    フーリエ変換を実行するためのプログラムを記憶し、 そのプログラムは、 各プロセッサにより、３次元データ空間の第１の座標軸
    に関する第１の変換処理、上記３次元データ空間の第２
    の座標軸に関する第２の変換処理、上記３次元データ空
    間の第３の座標軸に関する第３の変換処理を順次実行
    し、 上記複数のプロセッサの各々による、上記第１、第２の
    変換処理のいずれか一方の変換処理の実行後に、上記複
    数のプロセッサでのその一方の変換処理を実行した結果
    得られた一群の結果データを構成する複数の結果データ
    部分群がそれぞれ異なるプロセッサに割り当てられるよ
    うに、上記一群の結果データを上記複数のプロセッサの
    間で交換するステップを実行するようにプログラムさ
    れ、 上記第１から第３の変換処理は、一群の順序づけられた
    変換対象データに対する一群の順序づけられたフーリエ
    変換係数データを構成する複数のフーリエ変換係数デー
    タ部分群をそれぞれ異なるプロセッサにより生成するよ
    うに定められ、 上記変換対象データ群の各々は、上記３次元データ空間
    の直方体形状に位置する格子点群の一つの座標をそれぞ
    れ有し、 上記複数の変換対象データ部分群は、上記３次元データ
    空間の第３の座標軸の座標値が同じであり、上記３次元
    データ空間の第１、第２の座標軸の座標値が異なる全て
    の変換対象データが同一の変換対象データ部分群に含ま
    れるように定められ、 上記フーリエ変換係数データ群の各々は、３次元係数空
    間の直方体形状に位置する格子点群の一つの座標をそれ
    ぞれ有し、 上記複数のフーリエ変換係数データ部分群は、上記３次
    元係数空間の第１の座標軸の座標値が同じであり、上記
    ３次元波数空間の第２、第３の座標軸の座標値が異なる
    全てのフーリエ変換係数データが同一のフーリエ変換係
    数データ部分群に含まれるように定められ、 各プロセッサには、上記一群の変換対象データを構成す
    る複数の変換対象データ部分群の一つの変換対象データ
    部分群が予め割り当てられ、 上記一群のフーリエ変換係数データのそれぞれを生成し
    たプロセッサの順序が、上記一群の変換対象データのそ
    れぞれが割り当てられたプロセッサの順序と同一となる
    ように、上記交換するステップで各プロセッサに割り当
    てられる結果データ部分群が定められているもの。
16. 【請求項１６】各プロセッサは、上記一群の変換対象デ
    ータを構成する複数の変換対象データ部分群の内のその
    プロセッサに対して予め割り当てられた一つの変換対象
    データ部分群に対して、上記第１の変換処理を他のプロ
    セッサと並行して実行して、第１の結果データ部分群を
    生成し、 各プロセッサは、上記交換するステップが上記第２の変
    換処理の実行後に実行されるときには、その各プロセッ
    サで実行した上記第１の変換処理の結果生成した第１の
    結果データ部分群に対して、上記交換するステップが上
    記第１の変換処理の実行後に実行されたときには、その
    交換ステップでそのプロセッサに割り当てられた一つの
    第１の結果データ部分群に対して、上記第２の変換処理
    を他のプロセッサと並行して実行して、一つの第２の結
    果データ部分群を生成し、 各プロセッサは、上記交換するステップが上記第１の変
    換処理の実行後に実行されるときには、その各プロセッ
    サで実行した上記第２の変換処理の結果生成した一つの
    第２の結果データ部分群に対して、上記交換するステッ
    プが上記第２の変換処理の実行後に実行されたときに
    は、その交換ステップでそのプロセッサに割り当てられ
    た一つの第２の結果データ部分群に対して、上記第３の
    変換処理を他のプロセッサと並行して実行して一つのフ
    ーリエ変換係数データ部分群を生成する請求項１５記載
    のプログラム記録媒体。
17. 【請求項１７】上記変換対象データ群が上記３次元デー
    タ空間に直方体形状に位置する格子点群に上記３次元空
    間に第３の座標軸、第２の座標軸、第１の座標軸の順に
    順次割り当てられ、 上記第１から第３の変換処理は、上記複数のフーリエ変
    換係数データが、３次元係数空間に直方体形状に位置す
    る格子点群に、当該３次元係数空間の第１、第２、第３
    の座標軸の順序で割り当てられるように定められている
    請求項１５記載のプログラム記録媒体。
18. 【請求項１８】各プロセッサが上記第１の変換処理によ
    り生成する上記一つの第１の結果データ部分群は、上記
    ３次元データ空間の第３の座標軸の座標値が所定の同じ
    値であり、上記３次元データ空間の第２の座標軸の座標
    値と上記３次元係数空間の第１の座標軸の座標値が異な
    る値を有する全ての複数の第１の結果データを含み、 上記交換ステップが上記第１の変換処理が上記複数のプ
    ロセッサにより実行された後に実行され、 上記複数のプロセッサは、この交換ステップで、上記３
    次元係数空間の第１の座標軸の座標値が所定の同じ値で
    あり、上記３次元データ空間の第２、第３の座標軸の座
    標値が異なる値を有する全ての複数の第１の結果データ
    を含む第１の結果データ部分群が同一のプロセッサに割
    り当てられるように、上記複数のプロセッサが生成した
    一群の第１の結果データを上記複数のプロセッサの間で
    交換し、 各プロセッサが上記第２の変換処理により生成する上記
    一つの第２の結果データ部分群は、上記３次元係数空間
    の第１の座標軸の座標値が所定の同じ値であり、上記３
    次元波数空間の第２の座標軸の座標値と上記３次元デー
    タ空間の第３の座標軸の座標値が異なる値を有する全て
    の複数の第２の結果データを含み、 各プロセッサが上記第３の変換処理により生成する上記
    一つのフーリエ変換係数部分群は、上記３次元係数空間
    の第１の座標軸の座標値が所定の値であり、上記３次元
    波数空間の第２、第３の座標軸の座標値が異なる値を有
    する全ての複数のフーリエ変換係数を含む請求項１５記
    載のプログラム記録媒体。
19. 【請求項１９】各プロセッサが上記第１の変換処理によ
    り生成する上記一つの第１の結果データ部分群は、上記
    ３次元データ空間の第３の座標軸の座標値が所定の同じ
    値であり、上記３次元データ空間の第２の座標軸の座標
    値と上記３次元係数空間の第１の座標軸の座標値が異な
    る値を有する全ての複数の第１の結果データを含み、 上記交換ステップが上記第２の変換処理が上記複数のプ
    ロセッサにより実行された後に実行され、 各プロセッサが上記第２の変換処理により生成する上記
    一つの第２の結果データ部分群は、上記３次元データ空
    間の第３の座標軸の座標値が所定の同じ値であり、上記
    ３次元係数空間の第１、第２の座標軸の座標値が異なる
    値を有する全ての複数の第２の結果データを含み、 上記複数のプロセッサは、上記交換ステップにより、上
    記３次元係数空間の第１の座標軸の座標値が所定の同じ
    値であり、上記３次元係数空間の第１の座標軸の座標値
    と上記３次元データ空間の第３の座標軸の座標値が異な
    る値を有する全ての複数の第１の結果データを含む第１
    の結果データ部分群が同一のプロセッサに割り当てられ
    るように、上記複数のプロセッサが生成した一群の第１
    の結果データを上記複数のプロセッサの間で交換し、 各プロセッサが上記第３の変換処理により生成する上記
    一つのフーリエ変換係数部分群は、上記３次元係数空間
    の第１の座標軸の座標値が所定の値であり、上記３次元
    係数空間の第２、第３の座標軸の座標値が異なる値を有
    する全ての複数のフーリエ変換係数を含む請求項１５記
    載のプログラム記録媒体。
20. 【請求項２０】計算機により読みとり可能なプログラム
    記録媒体であって、複数のプロセッサを有する計算機で
    フーリエ変換を実行するためのプログラムを記憶し、 そのプログラムは、 各プロセッサにより、３次元空間の第１、第２、第３の
    座標軸の座標にそれぞれ関する第１、第２、第３の変換
    処理を順次かつ他のプロセッサと並行して実行し、 各プロセッサが上記第１、第２の変換処理のいずれか一
    方を実行した後に、その一方の変換処理の結果それぞれ
    のプロセッサで得られた複数の結果データを上記複数の
    プロセッサの間で交換するステップを実行するようにプ
    ログラムされ、 ここで、一群の順序づけられた変換対象データが上記３
    次元空間に直方体の形に並べられ、 上記第１から第３の変換処理は、上記一群の変換対象デ
    ータに対する一群の順序づけられた３次元空間の座標を
    有する複数のフーリエ変換係数データを生成するように
    定められ、 上記複数の変換対象データが構成する上記直方体を分割
    する上記３次元空間の上記第１の座標軸に垂直な複数の
    面の各々に含まれる複数の変換対象データが同一のプロ
    セッサに割り当てられ、 上記交換ステップは、上記一方の変換処理の結果得られ
    た上記複数の結果データが構成する３次元空間の直方体
    を、その３次元空間の第１の座標軸に垂直な複数の面で
    分割し直して、各面に属する複数の結果データを同一の
    プロセッサに割り当てるように、上記一方の変換処理の
    結果得られた上記複数の結果データを上記複数のプロセ
    ッサ間で交換するステップを有するもの。
21. 【請求項２１】上記一群の順序づけられた変換対象デー
    タが上記３次元空間に直方体の形に並べられる順序は、
    第３の座標軸、第２の座標軸、第１の座標軸の順であ
    り、 上記第１から第３の変換処理は、上記複数のフーリエ変
    換係数データが、３次元空間に、第１、第２、第３の座
    標軸の順序で並べられるように定められ請求項２０記載
    のプログラム記録媒体。
22. 【請求項２２】各プロセッサがパイプライン演算器を含
    み、その演算器での演算の対象とするループ長がＬのと
    きのその各プロセッサの演算性能を求めるための性能デ
    ータを上記複数のプロセッサに共通に記憶し、 その性能データを用いて、上記直方体の上記第１、第
    ２、第３の座標軸方向の長さを決定し、 その決定された上記第１、第２、第３の座標軸方向の長
    さを有する直方体に、上記順序づけられた複数の変換対
    象データを並べるステップをさらに有する請求項２０記
    載のプログラム記録媒体。
23. 【請求項２３】上記プログラムは、各プロセッサ上で実
    行されるアプリケーションプログラムから呼び出され、
    そのアプリケーションプログラムが指定するフーリエ変
    換対象データに対してフーリエ変換を実行し、生成した
    フーリエ変換係数データをそのアプリケーションプログ
    ラムに戻すライブラリである請求項１２から２２のいず
    れか一つに記載のプログラム記録媒体。
24. 【請求項２４】複数のプロセッサを有する計算機で実行
    するためのシミュレーション方法であって、 各プロセッサにより、シミュレーションすべき物理現象
    を支配する方程式に基づいて、シミュレーション対象の
    物理空間の異なる点での少なくとも一つの物理量の値を
    計算し、 各プロセッサにより、その計算に当たり、算出された複
    数の値を表すデータにし対してフーリエ変換を実行する
    ステップとを有し、 上記フーリエ変換を実行するステップは、請求項１から
    １１のいずれか一つにより実行されるもの。
25. 【請求項２５】上記物理現象は物理的な装置の動作に関
    連する物理現象であり、 上記算出された複数の値に基づいて、上記物理的な装置
    を設計するためのデータを生成するステップをさらに有
    する請求項２４記載のシミュレーション方法。
26. 【請求項２６】上記物理的な装置は、半導体装置である
    請求項２５記載のシミュレーション方法。
27. 【請求項２７】上記物理現象は気象であり、 上記算出された複数の値に基づいて、気象予測として報
    じるためのデータを生成するステップをさらに有する請
    求項２４記載のシミュレーション方法。
28. 【請求項２８】計算機により読みとり可能なプログラム
    記録媒体であって、複数のプロセッサを有する計算機で
    シミュレーションを実行するためのプログラムを記憶
    し、そのプログラムは、請求項２４から２７のいずれか
    一つによりシミュレーションを実行するようにプログラ
    ムされているもの。