JP4607796B2

JP4607796B2 - 共用メモリ型スカラ並列計算機向け、高速３次元フーリエ変換処理方法

Info

Publication number: JP4607796B2
Application number: JP2006058875A
Authority: JP
Inventors: 誠中西
Original assignee: Fujitsu Ltd
Current assignee: Fujitsu Ltd
Priority date: 2006-03-06
Filing date: 2006-03-06
Publication date: 2011-01-05
Anticipated expiration: 2026-03-06
Also published as: JP2007241349A; US7467053B2; US20070208795A1

Description

本発明は、共用メモリ型スカラ並列計算機において、３次元フーリエ変換を高速に処理するための方法に関する。

昨今、共用メモリ型スカラ並列計算機は、科学技術計算の分野、特に、流体解析、宇宙物理、気象、衝突解析、画像解析などの広汎な分野で利用されている。共用メモリ型スカラ並列計算機をこれらに利用する場合、大規模な問題を高速に処理するために、多次元フーリエ変換の高速な性能が要求される。

また、コンピュータ技術の発達と共に、共用メモリ型スカラ並列計算機のスカラＣＰＵのアーキテクチャも変化している。
PRIMEQUESで採用されたItanium2プロセッサでは、浮動小数点レジスタが128個設けられている。また、キャッシュおよびメモリのアクセス時間も比較的高速になっている。Sparc 向けなどに開発したフーリエ変換方法は高速な１次元変換をベースに多次元フーリエ変換を構成する方法であった。Sparc システムは、キャッシュがL1,L2 と２段となっており、L1キャッシュのアクセスに比べ、L2キャッシュのアクセス速度は遅い。このため、比較的小さなL1キャッシュにデータを保持して計算する方法が優位となる。

一般的な高速フーリエ変換についての説明は、非特許文献１を参照されたい。
Charles Van Loan, "Computational Frameworks for the Fast Fourier Transform", Society for Industrial and Applied Mathematics, 1992

レジスタが豊富でキャッシュの容量が大きくかつメモリアクセス・キャッシュのアクセスが、特に、連続なアクセスで比較的高速なマシンの高性能なＣＰＵの性能を限界まで引き出す上では、Sparc向けに開発されたフーリエ変換方法は十分ではない。また、Sparcでは、レジスタ数が比較的少なかったため、大きな基数を利用して演算密度を上げたり、メモリアクセスを、各計算で必要なタイミングより幾分前のタイミングでロードする方法を利用するには、レジスタが不足して、かえって遅くなる問題があった。

Itanium2プロセッサのようにレジスタが比較的多くあり、多くのデータをレジスタに保持して高速に計算ができ、メモリ・キャッシュのデータの連続アクセスが高速で、キャッシュの容量が大きなＣＰＵからなる計算機システムで高速なフーリエ変換の方法を見出すことは重要である。

本発明の課題は、レジスタ数が多く、メモリ・キャッシュのデータへの連続アクセスが高速なプロセッサに適した３次元フーリエ変換の処理方法を提供することである。

本発明の３次元フーリエ変換処理方法は、連続したデータ領域へのアクセスが高速化された、１次元目のインデックスをｉ、２次元目をｊ、３次元目をｋとしたとき、入力データｘ（ｉ、ｊ、ｋ）のｉの連続する方向に並んだデータを連続領域に格納する共有メモリ型スカラ並列計算機における３次元フーリエ変換処理方法において、インデックスｉとｊからなる２次元配列データｘ（ｉ、ｊ、＊）（＊は、ｋの取りうる任意の値）を単位として、キャッシュにロードし、インデックスｉが変化する方向への１次元目のフーリエ変換を、ｉを連続に変化させながら、ｊ方向の多重度分行い、インデックスｊが変換する方向への２次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、インデックスｋが変化する方向への３次元目のフーリエ変換を行なうことを特徴とする。

本発明によれば、レジスタ数が多く、メモリ・キャッシュのデータへの連続アクセスが高速なプロセッサに適した３次元フーリエ変換の処理方法を提供することができる。

本発明の実施形態では、スカラ計算機向けのフーリエ変換を１次元変換ベースで構築するよりも、より多くのフーリエ変換を多重に行うことを考える。つまり、１次元と２次元を組み合わせた平面ごとの変換を行い、そのあと３次元目の変換を行うようにする。特に、３次元目の変換を行うときに、１次元目と３次元目からなる平面を考える。こうすると、１次元目の多重度方向へのデータのアクセスが連続になるため、多重フーリエ変換をデータのアクセスが連続となる方向に行うことを考える。３次元目の変換は、２次元目を固定して、１次元目と３次元目で構成される平面に関するフーリエ変換を行うことで処理する。

特に、３次元目のフーリエ変換では。フーリエ変換をいくつかの基数のフーリエ変換に分解して行なう。このとき、基数毎の変換を１次元目の多重度を利用した多重変換で行なう。
（１）３次元目と１次元目とがなす平面に対する多重変換は、キャッシュのＷＡＹ数を考慮して、比較的小さな基数で計算することで、連続アクセスの高速性能を引き出すことができる。大きな基数を使うとキャッシュのＷＡＹ数による衝突から、メモリからのアクセス性能が劣化してしまう。すなわち、フーリエ変換は、変換を行う各次元のデータ長を適当な因子に分解したときの各因子に関するフーリエ変換を行うことで計算することができるので、データ長を因数である基数に分解して、基数の大きさのフーリエ変換を繰り返し行なうことで、データ長分の計算をする。たとえば、データ長を２の階乗で表したときは、２×２の複素数からなるフーリエ変換を繰り返し行なうことで、データ長分のフーリエ変換を行なうことができる。したがって、この場合、基数×基数個の複素数は、キャッシュを用いて読み書きしながら、まとめて演算する必要がある。ところが、キャッシュにまとめて読み書きできるデータ数は、ＷＡＹ数によって限定されることから、基数がＷＡＹ数に比較して所定以上大きくなると演算しなければならないデータの量がキャッシュの読み書き可能量を超えてしまい、データのロードが発生する。データのロードは演算速度に比べて遅いため、フーリエ変換の演算速度が落ちてしまう。したがって、高速にフーリエ変換を行なうために、基数は、キャッシュのＷＡＹ数に比べて比較的小さい必要がある。
（２）特に、データ長が２の冪で表されるときが重要である。Cooley-Tukeyの方法では、各基数の変換を行った後にビットリバースのデータの並べ替えがある。

２次元目の変換では、ビットリバースの並べ替えに関して、２次元目の並び替えで、データを１次元目方向にデータが並んだベクトルと考えて並び替えを行う。このベクトルの並び替えを行うときに、１次元目のベクトル内のビットリバースの並び替えも行う。こうすることで、１次元目と２次元目のビットリバースを同時に行なうことにより、キャッシュ・レジスタ上での並び替えとして行うことができ、ビットリバースする際にデータをロードする回数を少なくできる。
（３）同じく、ビットリバースのデータの並び替えを必要とするCooley-Tukeyの方法を使う、３次元目のフーリエ変換に関しては、２次元目を固定した平面で１次元目を多重度と考えた、基数の変換を行う。すべての基数に関する変換が終わった段階で、ビットリバースを行うが、１次元目の方向のデータの並びからなるベクトルごとの入れ替えでビットリバースを行う。すべての２次元目に値に対する多重変換が終われば３次元フーリエ変換の全演算の終了である。

特に、最後の基数は、比較的大きな基数により計算を高速化するため、基数として１６、３２などを利用する。ただし、キャッシュのＷＡＹ数による衝突は起こらないようにする。各基数の計算においては、基数の変換と回転因子の演算をまとめて行う。フーリエ変換演算の最後の部分は連続領域に対する基数の変換だけとなり、回転因子の演算を行なうためのテーブルの参照は不要であり高速な計算が可能である。
（４）１次元目および２次元目の計算に関しては、この方法以外にも、かなり大きな基数を使い（ただし、キャッシュのＷＡＹ数による衝突は起こらないようにする）、ビットリバースの並び替えの不要な方法を利用する方法も有効である。この方法は、１次元目方向と２次元目方向からなる平面は、データのアクセス領域としては連続領域となっており、本実施形態で前提とするプロセッサでは連続領域へのアクセスが速い事と、十分大きなキャッシュに、演算に必要なデータ全体が格納され得ることと、レジスタが豊富にあることを利用して行う。
（ａ）１次元変換では、以下の方法を利用する。

データ長を因数分解する。各因数が基数である。基数として、８または１６を使う。基数をｎ_{ｒａｄｉｘ}，ｍ＝ｎ／ｎ_{ｒａｄｉｘ}とし（ｎはデータ長）、ｍ×ｎ_{ｒａｄｉｘ}の配列をｎ_{ｒａｄｉｘ}×ｍの配列に転置する要領で、フーリエ変換すべきデータの、長さｎ_{ｒａｄｉｘ}の行ベクトルを１本ずつ読み出し、回転因子を掛けてフーリエ変換計算したものを、これを転置した配列に格納する。計算処理されるべき残りのデータについても、同様にして、基数の値ｍ回だけ繰り返す。最後の行に対応する演算データの格納に関しては、１次元目の値ｍが、（（ｎの約数であるところの基数の総個数ｍ）−１）次元の配列を表すとみなして、これらの順番を入れ替える。

この、アクセスパターンは、参照に関しては基数の数だけある連続な１次元ベクトルを連続に順次ロードするため高速である。
更に、基数個のフーリエ変換をｍ回行うループで、演算すべき行ベクトルデータを、一回前のループでロードしておいて、レジスタに保持しておく。このことで、データのロードと計算を同時に行うことができる。これは、レジスタが豊富にあるからできることである。基数を大きくすると、１回の参照で読み出されるデータの数が増えるので、１回の参照でロードされたデータに行なうべき演算量が増えるため、データのロードを一回前のループの演算中に完了することが可能である。１回のロードで行なうべき演算量が増えるが、演算に使うレジスタと事前にロードするために使うレジスタの数がＣＰＵのレジスタ総数を超えないようにする。

２次元目の変換に関しては、１次元目の多重度を先に行うようにする。この多重部分で基数回のフーリエ変換の繰り返しが発生する。この繰り返しでも同じように先読みすることができる。

以上の意味するところは、３次元フーリエ変換すべきデータを配列ｗ（ｉ、ｊ、ｋ）に格納した場合、１次元目のインデックスｉが連続しているデータは、格納領域に連続して格納されているため、１次元目のインデックスｉについて連続してアクセスするようにすれば、本発明の実施形態で前提とするプロセッサは連続アクセスが速いので、高速に処理ができるというものである。１次元目のフーリエ変換をする場合には、通常、ｉを連続に変化させて順次計算するが、２次元目の計算をする場合には、通常だと、ｊを連続に変化させて順次計算する。しかし、本発明の実施形態では、ｉを連続に変化させるループをｊを連続に変化させるループの中に入れて、２次元目のフーリエ変換においても、１次元目のｉを連続に変化させ、データアクセスを連続アクセスとすることにより、高速化を図る。
（５）混合基底の場合
混合基底（複数の異なる基底を用いる）の場合は、各基数に関する変換と、回転因子の乗算、および、転置を組み合わせたもので、各基数に対する変換を行う。このとき、入力領域と出力領域が必要になる。そのため、１次元目と２次元目に関しては、１次元目方向と２次元目方向とで構成されるひとつの平面とみなして、変換を行う。２次元目の変換は、１次元目の変換と同じ大きさの作業域を確保して、各基数の１次元目を多重度と見なして行なう。フーリエ変換の入出力はデータの格納された領域と作業域を切り替えて行う。

以下のようなフーリエ変換を利用する。
データ長n がp ×q と因数分解できたとき次の４ステップで計算できる。
step 1: 格納域x をx(p,q)の２次元配列と考える。１次元目（pのインデックス）に関する長さq のフーリエ変換を多重度p だけ行なう。

step 2: x(p,q)の要素x(i,j)にω**(i*j) を掛ける。ここで、ω＝exp(-2π/n) である。
step 3: x(p,q)を転置してw(q,p)に格納する。

step 4: w の２次元目（qのインデックス）に関して、長さp のフーリエ変換を多重度q だけ行なう。
qが小さく、pが大きな数である場合には、２つの因子への分解をstep 4の長さp の変換に繰り返し再帰的に提供することで、長さの短いフーリエ変換の組み合わせで長さn のフーリエ変換を計算することができる。

混合基底の場合には、以下のようなフーリエ変換を利用する。
データ長n をp 個の基数の積に分解する。n=n₁×n₂×．．．×n_p
s1: m ←n/n_pとする。

s2: データの格納されている１次元領域をx(m,n_p) なる2 次元配列と見なす。
s3: 2 次元目の基数n_pのフーリエ変換を多重度m だけ行なう。
s4: x(m,n_p) に回転因子を掛ける。

s5: x を転置してw(n_p,m) に格納する。
s6: m を(m/n_p-1,n_p-1) と 2次元に分解する。そのためm ←m/n_p-1と更新する。
s7: w(n_p,m,n_p-1)と一つ次元をあげて、下記のＳＴＥＰ１の計算を適用する。
ＳＴＥＰ１：まず、３次元目のn_p-1のフーリエ変換を多重度n_p×m 行なう。

s8: 回転因子をw(*,m,n_p-1) に掛ける。＊は、０以上ｎ_ｐ以下の任意の値。
s9: 2次元目と３次元目を入れ替え（転置）する。
x(n_p,n_p-1,m) ←w(n_p,m,n_p-1)。

s10:以下配列x およびw を入出力として入れ替えながら, m を次の基数との積に分解してフーリエ変換を繰り返す。
s11:最後の基数に関してフーリエ変換を行なうとき、すべての基数の個数が偶数ならw に奇数ならx に変換すべきデータがある。最後の多重転換は、変換すべきデータがwにあるときは、出力はx に、変換すべきデータがx にあるときはx に出力する。

２次元目の変換は、各基数に関するフーリエ変換と回転因子の積を、１次元目の多重度分さらに行なう。このため、作業域w は１次元分の多重度を加えたn₁×n₂の大きさが必要になる。

３次元目のフーリエ変換も同様に、２次元目を固定して 1次元目を多重度とした多重変換として計算を行う。すべての２次元目に値に対する多重変換が終われば、３次元フーリエ変換処理の終了である。

また、各基数の計算で多重度分の繰り返しがある。このため、一つ前の繰り返し処理で次に必要なデータの先読みを行うようにすることで高速性能が引き出せる。
以上のような方法で、Ｓｐａｒｃ向けのコードをＩｔａｎｉｕｍ２に単純に移植したときの性能に比べて2 倍強の性能を引き出す効果があった。

以下に、本発明の実施形態について、より詳しく説明する。
（１）基数が２の冪となる場合。
１次元目に関して、ビットリバースの並べ替えを除いて、上記で説明したように、基数次元の行ベクトルと見て、１本ずつ、フーリエ変換と回転因子の乗算を行う。次に、２次元目に関して、１次元目を多重度として２次元目のフーリエ変換と回転因子の乗算を行う。これらは２，４，８などの基数を用いて行う。基数をｎ_{ｒａｄｉｘ}として、ｍ＝ｎ／ｎ_{ｒａｄｉｘ}として、ｍ×ｎ_{ｒａｄｉｘ}の２次元配列とみなして、各行に対して基数ｎ_{ｒａｄｉｘ}のフーリエ変換と回転因子を掛ける。この結果は、ビットリバース順に並び替えておく。

次にｎ_{ｒａｄｉｘ}を更新して、ｍ←ｍ/ ｎ_{ｒａｄｉｘ}のようにｍを更新して、一つ前のｎ_{ｒａｄｉｘ}本ある、長さｍのベクトルに関して、同様に２次元配列とみなして、フーリエ変換と回転因子を掛ける。２次元目に関して、１次元目を多重度として連続アクセスで演算を行う。

ビットリバースに関しては、入れ替えを行う要素のペアを作りリストを作る。リストは以下のように作成できる。フーリエ変換の次数をｎとしたときｍ＝ｌｏｇ２（ｎ）とする。
Ａ）フーリエ変換及び回転因子の乗算がされた後のデータｗ（ｉ、ｊ、ｋ）のインデックスｉ、ｊ、ｋのそれぞれについてビットリバースを行なうが、今、ｉの総ビット数が奇数とすると、このとき、ｉを５つの部分に分ける。左が高位ビットとする。

インデックスをビットで表して、高位ビットからＢ1 ｜Ｃ｜Ｂ２｜Ａ｜Ｂ3 と５つに分ける、Ｂ１，Ｂ２，Ｂ３は 1ビットずつ振り分けて、ＡとＣは同じビット数とする。
ａ）ＡとＣＲが異なる場合
Ａがｍビットとして、０から２＊＊ｍ−１まで変化させることを考える。

ＡをビットリバースしたものをＡＲとあらわす。ＣをビットリバースしたものをＣＲとする。
Ａ＞ＣＲを成り立たせるすべてのＣを生成する。すると、ビット列Ｃ｜Ａとビット列ＡＲ｜ＣＲは異なる値となる。

また、ビット列Ｃ｜Ａとビット列ＡＲ｜ＣＲの右半分はＡ＞ＣＲを満たしたものを生成することになり、この２つのビット列ですべてのインデックス値を生成することになる。
このように生成したペア（Ｃ｜ＡとＡＲ｜ＣＲ）はすべての組み合わせをつくしている。

もし、重複しているペアがある場合には、２つのペア（Ｃ｜ＡとＡＲ｜ＣＲ）と（Ｃ２｜Ａ２とＡ２Ｒ｜Ｃ２Ｒ）が一致するペアとして、以下のように考える。
Ａは順次生成しているので、Ａ２＞Ａとする。

両者が一致するので、ＣＲとＡ２は一致しなければならない。すると、ＣＲ＞Ａとなるが、もともと、Ａ＞ＣＲとなるようにビット列Ｃ｜ＡとＡＲ｜ＣＲを生成しているはずなので矛盾する。したがって、重複するペアは発生しないことになる。

つまり、３ビットのビット列Ｂ＝Ｂ１｜Ｂ２｜Ｂ３はなんでもよいので０から７まで生成して、Ｂ1 ｜Ｃ｜Ｂ２｜Ａ｜Ｂ3 を生成する。ビットリバースしたものは次のようになる。

Ｂ３｜ＡＲ｜Ｂ２｜ＣＲ｜Ｂ１
ｂ）ＡとＣＲが一致する場合
したがって、ＡＲとＣも一致する。よって、Ｃ｜ＡとＡＲ｜ＣＲは同じビット列となる。したがって、Ｂ1 ｜Ｃ｜Ｂ２｜Ａ｜Ｂ３とＢ３｜ＡＲ｜Ｂ２｜ＣＲ｜Ｂ１が異なるビット列となるためには、Ｂ＝Ｂ１｜Ｂ２｜Ｂ３は非対称でないといけない。つまり、ＢＲをＢのビットリバースしたものとすると、Ｂ≠ＢＲ。これは３ビットのＢの具体的なパターンを見れば決めることができる。

このａ）とｂ）のペアを生成してテーブルの保管しておき、ビットの入れ替えを行う。
Ｂ）インデックスの総ビット数が偶数のとき
中央のＢ２を総ビット数が奇数のときに対して１ビット増やして２ビットにする。
ａ）ＡとＣＲが異なる場合
Ａがｍビットとして、０から２＊＊ｍ−１まで変化させることを考える。

ＡをビットリバースしたものをＡＲとあらわす。ＣのビットリバースしたものをＣＲと表す。
Ａ＞ＣＲを保つすべてのＣを生成する。するとＣ｜ＡとＡＲ｜ＣＲは一致しない。

また、ペアの右半分はＡ＞ＣＲを満たしたものを生成することになりすべて生成することになる。
このように生成したペア（Ｃ｜ＡとＡＲ｜ＣＲ）はすべての組み合わせをつくしている。

もし、生成しているペアが重複しているなら、２つのペア（Ｃ｜ＡとＡＲ｜ＣＲ）と（Ｃ２｜Ａ２とＡ２Ｒ｜Ｃ２Ｒ）がペアとして一致することになると仮定して以下のように考える。

Ａは順次生成しているので、Ａ２＞Ａとする。
ペアが一致しているので、ＣＲとＡ２は一致しなければならない。すると、ＣＲ＞Ａとなり矛盾する。したがって、仮定が間違っていたということになるので、生成したペアには重複したペアは発生しない。

つまり、４ビットのビット列Ｂ＝Ｂ１｜Ｂ２｜Ｂ３はなんでもよいので、０から15まで生成して、Ｂ1 ｜Ｃ｜Ｂ２｜Ａ｜Ｂ3 を生成すると、ビットリバースしたものは次のようになる。

Ｂ３｜ＡＲ｜Ｂ２Ｒ｜ＣＲ｜Ｂ１
ｂ）ＡとＣＲが一致する場合
Ｂ1 ｜Ｃ｜Ｂ２｜Ａ｜Ｂ３とＢ３｜ＡＲ｜Ｂ２Ｒ｜ＣＲ｜Ｂ１が異なるビット列となるためには、Ｂ＝Ｂ１｜Ｂ２｜Ｂ３は非対称でないといけない。つまり、Ｂ≠ＢＲ。これは４ビットのＢの具体的なパターンを見れば決めることができる。

このａ）とｂ）のペアを生成してテーブルの保管しておくことにより、インデックスの入れ替えを行うことができる。
２次元目のインデックスの入れ替えに関しては、次のようにする。２次元目の入れ替えを行うとき、１次元方向のベクトルで入れかえる。そのためベクトルＶ、ベクトルＷに入れ替える列ベクトルをコピーする。コピーしたら、このベクトル上で１次元目の入れ替えを行なう。１次元目の入れ替えを行なった結果のベクトルＶとＷを入れ替えて、もとの列ベクトルに格納する。

２次元目に関して、入れ替えが不要な場合は以下の場合である。
・ＡとＣＲが一致してＢとＢＲが一致する場合。
この場合のベクトルに関しては、１次元目の入れ替えを行なう。

３次元目に関しては、２次元目と同じように１次元目を多重度と見なしたフーリエ変換を行なった後、３次元目の入れ替えとして、２次元目と同様の１次元目方向の列ベクトルの入れ替えを行なう。
（２）ビットリバースを必要としない方法で、必要なデータを先読みする方法
次数ｎとして、基数ｎ_{ｒａｄｉｘ}に対してｍ＝ｎ／ｎ_{ｒａｄｉｘ}として（ｎはデータ長）、ｍ×ｎ_{ｒａｄｉｘ}の２次元配列と見なして、各行ベクトルに対してフーリエ変換と回転因子を掛ける。フーリエ変換後のデータの格納は、この基数のフーリエ変換の格納順とする。

そして、ｎ_{ｒａｄｉｘ}×ｍの配列に転置して格納する。この繰り返し処理で、フーリエ変換するデータの参照は、１つ前の計算中にロードしておく。基数は演算量が中くらいの８および１６がＩｔａｎｉｕｍ２では妥当である。

そして、繰り返し演算の最後のフーリエ変換に関しては、値ｍを、（（基数の数）−１）次元の配列の次元を表すと見なして、次元の入れ替えを行なう。つまり、演算後のデータをＡ（Ｉ１、Ｉ２、Ｉ３）⇒Ａ（Ｉ３、Ｉ２、Ｉ１）のように、次元を逆の順に入れ替える。これで、フーリエ変換が計算できる。

２次元変換において、１次元目の変換は、上記方法を適用して、２次元目の大きさだけの多重度分、１次元目の変換を繰り返し計算する。そして、２次元目の変換を１次元目の大きさの多重度だけ計算するときには、基数がかなり大きいため、フーリエ変換を行なうよりも大きな１次元配列に格納して、１次元変換を行なう。

または、基数を小さくした方法を２次元目の変換に対応させて、フーリエ変換を行なうことができる。ビットリバースの並び替えを使わない方法を利用すると、同じ大きさの作業域を割り当てる必要がある。

３次元目のフーリエ変換に関しても、１次元目を多重度として行なう。このときは基数を比較的小さな数とする方法を選ぶ。このように選ぶことで高速なメモリアクセスとデータのキャッシュ保持を効率よく利用できる。

図１は、本発明の実施形態の２次元面でのフーリエ変換の仕方を説明する図である。
ｎをデータ長として、ｎ_{ｒａｄｉｘ}を基数とした場合、ｍ×ｎ_{ｒａｄｉｘ}の配列Ｘから行ベクトルを順次取り出して、基数に関するフーリエ変換および回転因子の乗算を行い、それを転置した配列（ｎ_{ｒａｄｉｘ}×ｍの配列）に格納する。この繰り返し処理で、参照するデータは、一つまえの繰り返しでロードしておくことで、データのロードの待ち時間中に計算を実行することができる。これは、レジスタが豊富であることにより実現できるものである。

次に、ＷとＸを交換して、同じ形状の配列として同じアクセスパターン（配列Ｗにおいては、列ベクトルを順次取り出して）で、次の基数長の列ベクトルに対する変換を行なう。繰り返し演算の最後では回転因子の乗算がない。

繰り返し演算の最後に関しては、Ｘを基数の個数の多次元配列ｘ（p₂, p₃, ..., p_n, p₁)と見なし、Ｗをw(p_n, p_n-1, ..., p₂, p₁)と見なして、このようにデータが格納されるように最高次元以外、逆順に並び替える。

繰り返し演算の最後以外は、同じアクセスパターンで、かつ連続アクセスで計算できるので、高速化が図れる。
図２は、本発明の実施形態の３次元フーリエ変換の処理を説明するフローチャートである。

図２において、ステップＳ１０で、３次元目をスレッド数で均等に分割して、各スレッドに割り当てる。ステップＳ１１で、各スレッドに割り当てられた、１次元目と２次元目からなる平面単位に２次元フーリエ変換を行なう。ステップＳ１２において、２次元目をスレッド数で均等に分割して、各スレッドに１次元目と３次元目からなる２次元平面を各スレッドに割り当てる。ステップＳ１３において、３次元目を比較的小さな基数に分解して、各基数のフーリエ変換と回転因子の積を１次元目の多重度で行なう。ビットリバースが必要な方法ならば、最後にビットリバースを１次元目の大きさのベクトルとして行なう。

図３は、図２における３次元目のフーリエ変換の処理を示すフローチャートである。
３次元ＦＦＴの３次元目のフーリエ変換に関して、各スレッドに割り当てて行なう部分は以下のようにする。

スレッド毎に２次元目の大きさを均等に分割した部分を割り当てる。この最初と最後をそれぞれn2ds,n2de とする。各基底に対するフーリエ変換を１次元の多重度を持った１次元と３次元からなる平面に関して、２次元目をn2dsからn2deまで動かしながら計算する。x(k1,n2,n3) なる配列にデータが格納されているとする。２次元目をj2∈[n2ds,n2de] に固定して計算する。ビットリバースの並び替えの不要な方法だとw(k1,n3)なる２次元の作業域が必要になる。

図３において、ステップＳ１５において、３次元目のフーリエ変換のデータ長をｎ_３として、これを基数に分解する。分解した基数の数をｍとする。すなわち、以下のように分解する。
ｎ_３＝ｒ_１×ｒ_２×・・・×ｒ_ｍ
ステップＳ１６において、ｉ＝１とする。ステップＳ１７において、ｉがｍか否かを判断する。ステップＳ１７の判断がＹｅｓの場合には、ステップＳ２２に進む。ステップＳ１７の判断がＮｏの場合には、ステップＳ１８において、ｉが偶数か、すなわち、ｉと２の剰余（mod）が０か否かを判断する。ステップＳ１８の判断がＹｅｓの場合には、ステップＳ２０に進む。ステップＳ１８の判断がＮｏの場合には、ステップＳ１９において、作業域Ｗを入力、x(k1, j2, n3)を出力にして、基数ｒ_ｘのフーリエ変換と回転因子の積演算を１次元の大きさｎ１だけの多重度分行ない、ステップＳ２１に進む。ステップＳ２０では、x(k1, j2, n3)を入力に、wを出力にして、基数ｒ_ｘのフーリエ変換と回転因子の積演算を１次元の大きさｎ１だけの多重度分行ない、ステップＳ２１に進む。ステップＳ２１では、ｉを１増加し、ステップＳ１７に戻る。

ステップＳ２２においては、ｉが偶数か否かを判断する。ステップＳ２２の判断がＹｅｓの場合には、ステップＳ２４に進む。ステップＳ２２の判断がＮｏの場合には、ステップＳ２３に進む。ステップＳ２３では、作業域wを入力、x(k1, j2, n3)を出力にして、基数ｒ_ｘのフーリエ変換の計算を１次元の大きさｎ_１だけの多重度分行ない、処理を終了する。ステップＳ２４では、x(k1, j2, n3)を入力及び出力として、基数ｒ_ｘのフーリエ変換と回転因子の積の計算を１次元の大きさｎ_１だけの多重度分行なって、処理を終了する。

図４は、本発明の実施形態に従った、アクセスパターンが連続で、かつ同じパターンで最後のループ並び替えを同時に行ない、かつ、データの先読みを行なう１次元フーリエ変換のフローチャートである。

基本的に基数は８以上を使はないとデータの先読み効果は引き出せないことが分かっている。
まず、ステップＳ３０において、フーリエ変換のデータ長を以下のように基数に分解する。

ｎ＝ｒ_１×ｒ_２×・・・×ｒ_ｍ
ここで、基数の数をｍとする。
ステップＳ３１において、必要なテーブルをｍ−１個作成する。次に、ステップＳ３２において、ｉ＝１と設定する。ステップＳ３３において、最初のテーブルで基数ｒ_１のフーリエ変換を呼び出し、ｘを入力、ｗ１を出力として実行する。この処理は、サブルーチンfft1を呼び出すことにより行なう。ステップＳ３４において、ｉを１増加する。ステップＳ３５において、ｉと２の剰余が１か否かを判断する（ｉが奇数か否かを判断する）。ステップＳ３５の判断がＹｅｓの場合には、ステップＳ３６に進む。ステップＳ３５の判断がＮｏの場合には、ステップＳ３７に進む。ステップＳ３６では、ｉ番目のテーブルで基数ｒ_ｉのフーリエ変換を呼び出す。ｗ１を入力、ｗ２を出力とする。これは、サブルーチンfft1を呼び出すことによって行なう。ステップＳ３７では、ｉ番目のテーブルで基数ｒ_ｉのフーリエ変換を呼び出す。ｗ２を入力、ｗ１を出力とする。ここでは、サブルーチンfft1を呼び出す。そして、ステップＳ３８において、ｉを１増加し、ステップＳ３９において、ｉ＜ｍか否かを判断する。ステップＳ３９の判断がＹｅｓの場合には、ステップＳ３５に戻る。ステップＳ３９の判断がＮｏの場合には、ステップＳ４０において、最後の基数ｒ_ｍについてのフーリエ変換を呼び出して、結果のデータを格納する。ｍが基数ならｗ２を偶数ならｗ１を入力として、ｘを出力とする。ここでは、サブルーチンｆｆｔ２を呼び出す。

図５は、サブルーチンfft1の処理を示すフローチャートである。
Itanium2など浮動小数点レジスタが128 個以上あるプロセッサではスカラ変数を使って、Fortranなどの高級言語でプログラムを記述するとレジスタが十分豊富である間はレジスタを使ったコードとなる。また、これらの変数は、レジスタが割り当てられているものと考える。基数をｒ_ｘ、変換データ長をnとしたとき、m=n/r_xとする。入力の配列はx(m,r_x),出力はw(r_x,m) とする。回転因子を格納したテーブルtbl(1:r_x-1,m)もサブルーチンの入力とする。

ステップＳ４５において、x(1, 1:r_x)を変数tmp1, …, tmprxにロードする。ここで、１：ｒ_ｘというのは、１からｒ_ｘまでのそれぞれの値を意味する。そして、ｉを１に設定する。ステップＳ４６において、tmp1, …, tmprxを使って、基数ｒ_ｘのフーリエ変換で、これらを参照する計算を行い、結果をスカラ変数に格納する。ステップＳ４７において、tmp1, …, tmprxの参照が終わった直後で、x(i+1, 1:r_x)をtmp1, …, tmprxにロードする命令を発行する。このロードは、プログラムでは、単なる代入文で記述することができる。ステップＳ４８において、残りの基数ｒ_ｘのフーリエ変換の計算の続きを行なう。途中、tbl(1:r_x-1, i)をtbl1, …, tblrx-1にロードする命令を適当にある間隔でプログラムに埋め込む。ステップＳ４９において、基数ｒ_ｘのフーリエ変換の出力で、最初の出力を除くｒ_ｘ−１個の出力に、順にtbl1, …, tblrx-1の内容をかける。ステップＳ５０において、結果をw(1:r_x, i)に格納する。ステップＳ５１において、ｉを１だけ増加する。ステップＳ５２において、ｉ＜ｍか否かを判断する。ステップＳ５２の判断がＹｅｓの場合には、ステップＳ４６に戻る。ステップＳ５２の判断がＮｏの場合には、ステップＳ５３において、tmp1, …, tmprxを使って、基数ｒ_ｘのフーリエ変換で、これらを参照する計算を行い、結果をスカラ変数に格納する。ステップＳ５４において、残りの基数ｒ_ｘの計算を行なう。途中、tbl1(1:r_x-1, i)をtbl1, …, tblrx-1にロードする命令を適当に、ある間隔でプログラムに埋め込んでおく。ステップＳ５５において、基数ｒ_ｘのフーリエ変換の出力で最初の出力を除く、ｒ_ｘ−１個の出力に、順に、tbl1, …, tblrx-1の内容をかける。ステップＳ５６において、結果をw(1:r_x, i)に格納して、サブルーチンから抜け出る。

図６は、本発明の実施形態に従った、サブルーチンｆｆｔ２の処理を説明するフローチャートである。
Itanium2など浮動小数点レジスタが128 個以上あるプロセッサでは、スカラ変数を使ってFortran などの高級言語でプログラムを記述するとレジスタが十分豊富である間はレジスタを使ったコードとなる。また、これらの変数はレジスタが割り当てられているものと考える。基数をr_x、変換データ長をn としたとき、m=n/r_xとする。入力の配列はx(m,r_x),出力はw(r_x,m) とする。ｘ＝１のときは、m=r₂×r₃×... ×r_pと分解され、各基数に関して変換が終わっている。因子の数はp である。基数の個数の多次元配列ｘ（r₂, r₃, ..., r_p, r_x)と見なして、Ｗをw(r_p, r_p-1, ..., r₂,r_x)と見なして、そのように格納されるように最高次元以外のインデックスを逆順に並び替える。

図６において、ステップＳ６０においては、x(1, 1:r_x)をtmp1, …, tmprxにロードし、ｉを１に設定する。ステップＳ６１において、tmp1, …, tmprxを使って、基数ｒ_ｘのフーリエ変換で、これらを参照する計算を行い、結果をスカラ変数に格納する。ステップＳ６２において、tmp1, …, tmprxの参照が終わった直後で、x(i+1, 1:r_x)をtmp1, …, tmprxにロードする命令を発行する。実際のプログラムでは、単なる代入文で十分である。ステップＳ６３において、残りの基数ｒ_ｘのフーリエ変換の計算を行なう。ステップＳ６４において、
dis1=r₂, dis2=dis1*r₃, …, disp-1=disp-2*r_p
dit1=r_p, dit2=dit1*r_p-1, ..., ditp-1=ditp-2*r₂
ｊ←ｉ−１
とおき、
j=i0+i1*dis1+i2*dis2+...+ip-1*disp-1
となるｉ０、・・・、ｉｐ−１を決める。すなわち、
i0=mod(j, dis1), i←j/dis1
,………
ix=mod(j, disx+1), j←j/disx+1を繰り返す。
そして、
k=ip-1+ip-2*dit1+...+i0*ditp-1+1と計算する。

ステップＳ６５において、結果をw(1:r_x, k)に格納する。ステップＳ６６において、ｉを１増加し、ステップＳ６７において、ｉ＜ｍか否かを判断する。ステップＳ６７の判断がＹｅｓの場合には、ステップＳ６１に戻る。ステップＳ６７の判断がＮｏの場合には、ステップＳ６８において、tmp1, …, tmprxを使って、基数ｒ_ｘのフーリエ変換で、これらを参照する計算を行い、結果をスカラ変数に格納する。ステップＳ６９において、
dis1=r₂, dis2=dis1*r₃, …, disp-1=disp-2*r_p
dit1=r_p, dit2=dit1*r_p-1, ..., ditp-1=ditp-2*r₂
ｊ←ｉ−１
とおき、
j=i0+i1*dis1+i2*dis2+...+ip-1*disp-1
となるｉ０、・・・、ｉｐ−１を決める。すなわち、
i0=mod(j, dis1), i←j/dis1
,………
ix=mod(j, disx+1), j←j/disx+1を繰り返す。
そして、
k=ip-1+ip-2*dit1+...+i0*ditp-1+1と計算する。
ステップＳ７０において、結果をw(1:r_x, k)に格納して、サブルーチンから抜け出る。

図７は、２次元変換でビットリバースの並び替えを同時に行なう方法を説明するフローチャートである。
この方法も同じように、データ長n がp ×q と因数分解できたとき次の４ステップで計算する方法から構成することができる。n は２のべき乗であらわせる数とする。

ビットリバースを別個に行なう場合は以下の通りである。
step 1: 格納域x をx(p,q)の２次元配列と考える。2 次元目に関する長さq のフーリエ変換を多重度p だけ行なう。

step 2: x(p,q)の要素x(i,j)にω**(i*j) を掛ける。ω＝exp(-2π/n) 。
step 3: x(p,q)を転置してw(q,p)に格納する。
step 4: w の２次元目に関して、長さp のフーリエ変換を多重度q だけ行なう。

２つの因子への分解をstep 4の長さp の変換でstep 1およびstep 3のフーリエ変換の結果をビットリバースの並べ替えを行なった順序で格納し、step３の転置を行なわない。この結果、最後に全体をビットリバースの並び替えを行なうことで同じ計算を行なうことになる。

ビットリバースを同時に行なう場合は以下の通りである。
step 1: 格納域x をx(p,q)の２次元配列と考える。2 次元目に関する長さq のフーリエ変換を多重度p だけ行なう。

step 2: x(p,q)の要素x(i,j)にω**(i*j) を掛ける。ω＝exp(-2π/n) 。
step 3: x(p,q)に 2次元目に関してフーリエ変換の結果をビットリバースした並びで格納する。

step 4: x の１次元目に関して長さp のフーリエ変換を多重度q だけ行なう。結果は、ビットリバースの並び替えを行なった順序で格納する。
step 5: x(p×q)のデータに対して、ビットリバースの並び替えを行なう。

step 1〜step 3を繰り返し再帰的に提供することで、長さn のフーリエ変換を計算することができる。
図７においては、ステップＳ７７において、１次元目に関して、下記のs1〜s11の計算を行なう。これを２次元目の多重度分行なう。ステップＳ７６において、２次元目に関して、s1〜s11の計算を各基底に関してのフーリエ変換及び回転因子の積に関して、１次元目の多重度分行なう。そして、ステップＳ７７で、１次元目と２次元目のビットリバースの並び替えを同時に行なう。このときは、サブルーチンbitrb2dを呼び出して行なう。

データ長n をp 個の基数の積に分解する。n=n₁×n₂×．．．×n_p
s1: m ←n/n_pとする。
s2: データの格納されている１次元領域をx(m,n_p) なる2 次元配列と見なす。

s3: 2 次元目の基数n_pのフーリエ変換を多重度m 分行なう。
s4: x(m,n_p) に回転因子を掛ける。
s5: 基数n_pのフーリエ変換結果をビットリバースの並び替えを行なった順に格納する。

これは、基数が固定されているので、結果の格納位置を変えるように基数n_pのフーリエ変換を変更することで実現できる。
s6: m を(m/n_p-1,n_p-1) と 2次元に分解する。そのためm ←m/n_p-1と更新する。

s7: w(m,n_p-1,n_p)と一つ次元をあげて、下記のＳＴＥＰ１の計算を適用する。
ＳＴＥＰ１：まず、２次元目のn_p-1のフーリエ変換を多重度m ×n_p行なう。
s8: 回転因子をx(m,n_p-1,*) に掛ける。

s9: ２次元目の基数n_p-1のフーリエ変換の結果をビットリバースの並べ替えを行なった順序に格納する。
s10:m を次の基数との積に分解してフーリエ変換を繰り返す。

s11:最後の基数に関してフーリエ変換を行ない、これも結果をビットリバースの並び替えを行なった順で格納する。
s12:x(n)のデータに関してビットリバースの並び替えを行なう。

図８〜図１１は、サブルーチンbitrb2dの処理を説明するフローチャートである。
ステップＳ８０において、２次元目の大きさのビット数ｎ_ｂ（ｎ_２＝２＊＊ｎ_ｂ）を求める。ステップＳ８１において、ｎ_ｂが偶数か否かを判断する。ステップＳ８１の判断で、奇数と判断された場合には、ステップＳ９７に進む。ステップＳ８１で、偶数と判断された場合には、ステップＳ８２に進む。ステップＳ８２において、nbitx=(n_b-3)/2, nx=2**nbitxを計算する。ステップＳ８３において、ｉを１に設定する。ステップＳ８４において、ｊを１に設定する。ステップＳ８５において、ｊ＜ｉか否かを判断する。ステップＳ８５の判断がＮｏのときは、ステップＳ８９に進む。ステップＳ８５の判断がＹｅｓの場合には、ステップＳ８６において、（ｉ−１）を１ビット左にシフトしたものに、（ｊ−１）をビットリバースしたものを、更に（nbitx+2)ビット左にシフトしたものを加えたものをｉｄｘとする。ｉとｊを入れ替えて同じ操作をしたものをｉｓｘとする。ステップＳ８７において、０から７までのビットパターンを生成して、３ビットに左から分けて、Ｂ１、Ｂ２、Ｂ３とする。ｉｄｘの2*nbit+3ビット目にＢ１を、nbit+2ビット目にＢ２を、１ビット目にＢ３を加えて、それに、１を加えてｉｄとする。同様に、ｉｓｘの2*nbit+3ビット目にＢ３を、nbit+2ビット目にＢ２を、１ビット目にＢ１を加えて、それに１を加えてｉｓとする。０から７までのビットパターンを生成する毎に、このように、ｉｄおよびｉｓをつくり、サブルーチンｐｅｒｍｘｙを呼び出す。ｐｅｒｍｘｙは、ｘ（＊、ｉｓ）とｘ（＊、ｉｄ）を入れ替える。そのときに、ｘ（＊，ｉｓ）及びｘ（＊，ｉｓ）のビットリバースの並び替えを行なう。ステップＳ８８において、ｊを１だけ増加して、ステップＳ８５に戻る。ステップＳ８９では、ｉを１だけ増加する。

ステップＳ９０では、ｉ＞ｎ_ｘか否かを判断する。ステップＳ９０の判断がＮｏのときは、ステップＳ８４に戻る。ステップＳ９０の判断がＹｅｓの場合には、ステップＳ９１で、ｉを１に設定し、ステップＳ９２において、(i-1) を１ビット左シフトしたものに(i-1) をビットリバースしたものを(nbitx+2) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ９３において、｛001, 011 ｝の 2つのビットパターンに対して、B1,B2,B3 と 3ビットを取り出して、idx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに１を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。｛001,011 ｝の 2つのビットパターン毎に、idおよびisを作り、サブルーチンpermxyを呼び出す。permxyはx(*,is) とx(*,id) を入れ替える。そのときに、x(*,is) およびx(*,id) のビットリバースの並び替えを行なう。ステップＳ９４において、2 次元目の並び替えはないが、1 次元の並び替えが必要なものである、ビットパターン［000,010,101,111]に対して、一つ前の処理と同じようにidを計算する。perm1dを呼び出して、x(*.id) のビットリバースの並び換えをビットパターン毎に行なう。ステップＳ９５では、ｉを１だけ増加し、ステップＳ９６において、ｉ＞ｎ_ｘであるか否かを判断する。ステップＳ９６の判断がＮｏの場合には、ステップＳ９２に戻る。ステップＳ９６の判断がＹｅｓの場合には、サブルーチンから抜け出る。

ステップＳ９７では、nbitx=(n_b-4)/2、nx=2**nbitxを計算する。ステップＳ９８で、ｉを１に設定し、ステップＳ９９において、ｊを１に設定する。ステップＳ１００においては、ｊ＜ｉか否かを判断する。ステップＳ１００の判断がＮｏの場合には、ステップｓ１０４に進む。ステップＳ１００の判断がＹｅｓの場合には、ステップＳ１０１において、(i-1) を１ビット左シフトしたものに(j-1) をビットリバースしたものを(nbitx+3)ビット左にシフトしたものを加えたものをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ１０２において、０から16までのビットパターンを生成して、1ビット、2ビット、１ビットに左から分けてB1,B2,B3とする。idx の2*nbit+4ビット目にB1を、nbit+3,nbit+2 ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3,nbit+2ビット目にB2を反転いたものを、1 ビット目にB1を加えてそれに1 を加えてisとする。0 から16までのビットパターンを生成する毎に、このように、idおよびisを作り、サブルーチンpermxyを呼び出す。permxyはx(*,is) とx(*,1d) を入れ替える。そのときに、x(*,is) およびx(*,id) のビットリバースの並び替えを行なう。そして、ステップＳ１０３において、ｊを１だけ増加し、ステップＳ１００に戻る。ステップＳ１０４では、ｉを１だけ増加し、ステップＳ１０５において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１０５の判断がＮｏの場合には、ステップＳ９９に戻る。ステップＳ１０５の判断がＹｅｓの場合には、ステップＳ１０６に進む。

ステップＳ１０６では、ｉに１を設定する。ステップＳ１０７において、(i-1) を１ビット左シフトしたものに(i-1) をビットリバースしたものを(nbitx+3) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ１０８において、[0001,0010,0011,0101,0111,1011｝の6 つのビットパターンに対して、B1,B2,B3として、 1ビット、2ビット、1ビットをおのおの取り出してidx の2*nbit+4ビット目にB1を、,nbit+3,nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1 を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3,nbit+2ビット目にB2を反転したものを、1 ビット目にB1を加えて、それに1 を加えてisとする。上記６つのビットパターン毎にidおよびisを作り、サブルーチンpermxyを呼び出す。permxyはx(*,is) とx(*,id) を入れ替える。そのときに、x(*,is) およびx(*,id) のビットリバースの並び替えを行なう。ステップＳ１０９において、2 次元目の並び替えはないが、1 次元の並び替えが必要なものである、ビットパターン[0000, 1001, 0110, 1111]に対して、一つ前の処理と同じようにidを計算する。perm1dを呼び出して、x(*. id) のビットリバースの並び換えをビットパターン毎に行なう。ステップＳ１１０において、ｉを１だけ増加し、ステップＳ１１１において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１１１の判断がＮｏの場合には、ステップＳ１０７に戻る。ステップＳ１１１の判断がＹｅｓの場合には、サブルーチンから抜け出る。

図１２〜図１５は、サブルーチンpermxyの処理を表すフローチャートである。
サブルーチンpermxyは、x(1:n1,id2) およびx(1:n1,is2) を呼び出し元より受け渡されて、x(1:n1,id2) およびx(1:n1,is2) の長さn1のベクトルをビットリバースの並びに並び替えて、その結果の1 次元ベクトルを入れ替える。

ステップＳ１２０において、１次元目の大きさのビット数ｎ_ｂ（ｎ_１＝２＊＊ｎ_ｂ）を求める。ステップＳ１２１において、ｎ_ｂが偶数か否かを判断する。ステップＳ１２１において、奇数と判断された場合には、ステップＳ１３７に進む。ステップＳ１２１において、偶数と判断された場合には、ステップＳ１２２において、nbitx=(n_b-3)/2、nx=2**nbitxを計算する。ステップＳ１２３において、ｉを１に設定し、ステップＳ１２４において、ｊを１に設定する。ステップＳ１２５において、ｊ＜ｉか否かを判断する。ステップＳ１２５の判断がＮｏのときは、ステップＳ１２９に進む。ステップＳ１２５の判断がＹｅｓの場合には、ステップＳ１２６において、(i-1) を１ビット左シフトしたものに(j-1) をビットリバースしたものを(nbitx+2)ビット左にシフトしたものを加えたものをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ１２７において、０から７までのビットパターンを生成して、 3ビットに左から分けてB1,B2,B3とする。idx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1 を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。0 から7 までのビットパターンを生成する毎に、このように、idおよびisを作り、x(1:n1, id2),x(1:n1, is2)のx(id,id2),x(is,id2) とx(id,is2),x(is,is2) をロードし、おのおの入れ替えたものを格納するときに、x(id, is2),x(is, is2) およびx(id, id2), x(is, id2)の２つのベクトル間の入れ替えとして行なう。ステップＳ１２８において、ｊを１だけ増加し、ステップＳ１２５に戻る。

ステップＳ１２９では、ｉを１だけ増加し、ステップＳ１３０において、ｉ＞ｎｘであるか否かを判断する。ステップＳ１３１において、ｉを１に設定する。ステップＳ１３２において、(i-1) を１ビット左シフトしたものに(i-1) をビットリバースしたものを(nbitx+2) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ１３３において、｛001,011 ｝の 2つのビットパターンに対して、B1,B2,B3 と 3ビットを取り出してidx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。｛001,011 ｝の 2つのビットパターン毎にidおよびisを作り、x(1:n1,id2),x(1:n1,is2) のx(id,id2),x(is,id2) とx(id,is2),x(is,is2) をロードし、おのおの入れ替えたものを格納するときにx(id, is2),x(is, is2) およびx(id, id2),x(is, id2) の２つのベクトル間の入れ替えとして行なう。ステップＳ１３４において、１次元目の並び替えはないが、id2,is2 次元の並び替えが必要なものである、ビットパターン［000,010,101,111]に対して、一つ前の処理と同じようにidを計算する。(id. id2) とx(id, is2) の交換を各ビットパターン毎に行なう。そして、ステップＳ１３５において、ｉを１だけ増加し、ステップＳ１３６において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１３６の判断がＮｏの場合には、ステップＳ１３２に戻る。ステップＳ１３６の判断がＹｅｓの場合には、サブルーチンから抜け出る。

ステップＳ１３７においては、nbitx=(n_b-4)/2、nx=2**nbitxを計算する。ステップＳ１３８において、ｉを１に設定し、ステップＳ１３９において、ｊを１に設定する。ステップＳ１４０において、ｊ＜ｉか否かを判断する。ステップＳ１４０の判断がＮｏの場合には、ステップＳ１４４に進む。ステップＳ１４０の判断がＹｅｓの場合には、ステップＳ１４１において、(i-1) を１ビット左シフトしたものに(j-1) をビットリバースしたものを(nbitx+3)ビット左にシフトしたものを加えたものをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ１４２において、０から16までのビットパターンを生成して、1ビット、2ビット、1 ビットに左から分けてB1,B2,B3とする。idx の2*nbit+4ビット目にB1を、nbit+3,nbit+2 ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3, nbit+2ビット目にB2を反転いたものを、1 ビット目にB1を加えてそれに1 を加えてisとする。0 から16までのビットパターンを生成する毎に、このように、idおよびisを作り、x(1:n1, id2),x(1:n1, is2) のx(id, id2), x(is, id2) とx(id, is2), x(is, is2) をロードし、おのおの入れ替えたものを格納するときにx(id, is2), x(is, is2) およびx(id, id2), x(is, id2) の２つのベクトル間の入れ替えとして行なう。ステップＳ１４３において、ｊを１だけ増加させる。ステップＳ１４４において、ｉを１だけ増加させる。ステップＳ１４５において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１４５の判断がＮｏの場合には、ステップＳ１３９に戻る。ステップＳ１４５の判断がＹｅｓの場合には、ステップＳ１４６において、ｉを１に設定する。ステップＳ１４７において、(i-1) を１ビット左シフトしたものに、(i-1) をビットリバースしたものを(nbitx+3) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ１４８において、[0001,0010,0011,0101,0111,1011｝の6 つのビットパターンに対して、B1,B2,B3と 1,2,1ビットをおのおの取り出して、idx の2*nbit+4ビット目にB1を、nbit+3,nbit+2 ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3, nbit+2ビット目にB2を反転したものを、1 ビット目にB1を加えて、それに1 を加えてisとする。上記６つのビットパターン毎にidおよびisを作り、x(1:n1, id2), x(1:n1, is2) のx(id, id2), x(is, id2) とx(id, is2), x(is, is2) をロードし、おのおの入れ替えたものを格納するときにx(id, is2), x(is, is2) およびx(id, id2), x(is, id2) の２つのベクトル間の入れ替えとして行なう。ステップＳ１４９において、2 次元目の並び替えはないが、1 次元の並び替えが必要なものである、ビットパターン［0000, 1001, 0110, 1111]に対して、一つ前の処理と同じようにidを計算する。x(id.id2) とx(id,is2) の交換をビットパターン毎に行なう。ステップＳ１５０において、ｉを１だけ増加し、ステップＳ１５１において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１５１の判断がＮｏの場合には、ステップＳ１４７に戻り、ステップＳ１５１の判断がＹｅｓの場合には、サブルーチンから抜け出る。

図１６〜図１９は、サブルーチンperm1dの処理を説明するフローチャートである。
サブルーチンperm1dは、x(1:n1,id2) を呼び出し元より受け渡されて、x(1:n1, id2) の長さn1のベクトルをビットリバースの並びに並び替えを行なう。

ステップＳ１５５において、１次元目の大きさのビット数ｎ_ｂ（ｎ_１＝２＊＊ｎ_ｂ）を求める。ステップＳ１５６において、ｎ_ｂが偶数か否かを判断する。ステップＳ１５６における判断がＹｅｓの場合には、ステップＳ１７２に進む。ステップＳ１５６の判断がＮｏの場合には、ステップＳ１５７において、nbitx=(n_b-3)/2、nx=2**nbitxを計算する。ステップＳ１５８において、ｉを１に設定する。ステップＳ１５９において、ｊを１に設定する。ステップＳ１６０において、ｊ＜ｉか否かを判断する。ステップＳ１６０の判断がＮｏの場合には、ステップＳ１６４に進む。ステップＳ１６０の判断がＹｅｓの場合には、ステップＳ１６１において、(i-1) を１ビット左シフトしたものに、(j-1) をビットリバースしたものを(nbitx+2)ビット左にシフトしたものを加えたものをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ１６２において、０から７までのビットパターンを生成して、これを3ビットに左から分けてB1,B2,B3とする。idx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1 を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。0 から7 までのビットパターンを生成する毎に、このように、idおよびisを作り、x(id, id2)と, x(is, id2) とを入れ換える。ステップＳ１６３において、ｊを１だけ増加し、ステップＳ１６４において、ｉを１だけ増加する。ステップＳ１６５において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１６５の判断がＮｏのときは、ステップＳ１５９に戻る。ステップＳ１６５の判断がＹｅｓのときは、ステップＳ１６６において、ｉを１と設定する。そして、ステップＳ１６７において、(i-1) を１ビット左シフトしたものに、(i-1) をビットリバースしたものを(nbitx+3) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ１６８において、｛001,011 ｝の 2つのビットパターンに対して、B1,B2,B3 と 3ビットを取り出して、idx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。｛001,011 ｝の 2つのビットパターン毎に、idおよびisを作り、x(id,id2),x(is,id2) を交換する。ステップＳ１６９において、１次元目の並び替えはないが、id2,is2 次元の並び替えが必要なものである、ビットパターン［000, 010, 101, 111]に対して、一つ前の処理と同じようにidを計算する。x(id.id2) とx(id,is2)の交換をビットパターン毎に行なう。ステップＳ１７０において、ｉを１だけ増加し、ステップＳ１７１において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１７１の判断がＮｏの場合には、ステップＳ１６７に戻る。ステップＳ１７１の判断がＹｅｓの場合には、サブルーチンから抜け出る。

ステップＳ１７２においては、nbitx=(n_b-4)/2、nx=2**nbitxを計算する。ステップＳ１７３において、ｉを１に設定する。ステップＳ１７４において、ｊを１に設定する。ステップＳ１７５において、ｊ＜ｉか否かを判断する。ステップＳ１７５の判断がＮｏの場合には、ステップＳ１７９に戻る。ステップＳ１７５の判断がＹｅｓの場合には、ステップＳ１７６において、(i-1) を１ビット左シフトしたものに、(j-1) をビットリバースしたものを(nbitx+3)ビット左にシフトしたものを加えたものをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ１７７において、０から16までのビットパターンを生成して、1,2,1 ビットに左から分けてB1,B2,B3とする。idx の2*nbit+4ビット目にB1を、nbit+3,nbit+2 ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3,nbit+2ビット目にB2を反転いたものを、1 ビット目にB1を加え、それに1 を加えてisとする。0 から16までのビットパターンを生成する毎に、このように、idおよびisを作り、x(1:n1,id2) のx(id,id2),x(is,id2) を入れ換える。ステップＳ１７８において、ｊを１だけ増加し、ステップＳ１７５に戻る。

ステップＳ１７９においては、ｉを１だけ増加し、ステップＳ１８０において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１８０の判断がＮｏのときは、ステップＳ１７４に戻る。ステップＳ１８０の判断がＹｅｓの場合には、ステップＳ１８１において、ｉを１に設定する。ステップＳ１８２において、(i-1) を１ビット左シフトしたものに、(i-1) をビットリバースしたものを(nbitx+3) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ１８３において、[0001, 0010, 0011, 0101, 0111, 1011｝の6 つのビットパターンに対して、B1,B2,B3と 1,2,1ビットをおのおの取り出して、idx の2*nbit+4ビット目にB1を、nbit+3,nbit+2 ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3,nbit+2ビット目にB2を反転したものを、1 ビット目にB1を加えて、それに1 を加えてisとする。上記６つのビットパターン毎にidおよびisを作り、x(1:n1,id2) のx(id,id2),x(is,id2) を入れ換える。ステップＳ１８４において、ｉを１だけ増加し、ステップＳ１８５において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ１８５の判断がＮｏの場合には、ステップＳ１８２に戻る。ステップＳ１８５の判断がＹｅｓの場合には、サブルーチンから抜け出る。

図２０は、３次元目のフーリエ変換処理のフローチャートである。
x(1:n1,j2,1:3n) に関して 1次元の多重度分、 3次元目のフーリエ変換を行なう。
ステップＳ１９０において、３次元目の処理においては、前述のs1〜s11の計算を各基数について、各基数のフーリエ変換及び回転因子の積の処理を、１次元目の多重分だけ行なう。ステップＳ１９１において、１次元目と２次元目のビットリバースの並び替えを同時に行なう。この場合、サブルーチンｂｉｔｒｂ３ｄを呼び出して実行する。

図２１〜図２４は、サブルーチンbitrb3dの処理を説明するフローチャートである。
サブルーチンbitrb3d は、x(1:n1,j2,1:n3) 1次元目と 3次元目からなる平面で３次元目に関するビットリバースの並びに並び替えを１次元のベクトル単位で行なう。

ステップＳ１９５において、１次元目の大きさのビット数ｎ_ｂ（ｎ_３＝２＊＊ｎ_ｂ）を求める。ステップＳ１９６において、ｎ_ｂが偶数か否かを判断する。ステップＳ１９６の判断がＹｅｓの場合には、ステップＳ２１１に進む。ステップＳ１９６の判断がＮｏの場合には、ステップＳ１９７に進む。

ステップＳ１９７においては、nbitx=(n_b-3)/2、nx=2**nbitxを計算する。ステップＳ１９８において、ｉを１に設定し、ステップＳ１９９において、ｊを１に設定する。ステップＳ２００において、ｊ＜ｉか否かを判断する。ステップＳ２００の判断がＮｏの場合には、ステップＳ２０４に進む。ステップＳ２００の判断がＹｅｓの場合には、ステップＳ２０１において、(i-1) を１ビット左シフトしたものに、(j-1) をビットリバースしたものを(nbitx+2)ビット左にシフトしたものを加えたものをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ２０２において、０から７までのビットパターンを生成して、 3ビットに左から分けてB1,B2,B3とする。idx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1 を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。0 から7 までのビットパターンを生成する毎に、このように、idおよびisを作り、x(1:n1,j2,id),x(1:n1,j2,is) を入れ換える。ステップＳ２０３において、ｊを１だけ増加し、ステップＳ２００に戻る。

ステップＳ２０４では、ｉを１だけ増加し、ステップＳ２０５において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ２０５の判断の結果がＮｏの場合には、ステップＳ１９９に戻る。ステップＳ２０５の判断の結果がＹｅｓの場合には、ステップＳ２０６において、ｉを１に設定する。ステップＳ２０７において、(i-1) を１ビット左シフトしたものに、(i-1) をビットリバースしたものを(nbitx+2) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ２０８において、｛001,011 ｝の 2つのビットパターンに対して、B1,B2,B3 と 3ビットを取り出して、idx の2*nbit+3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+2ビット目にB2を、1 ビット目にB1を加えそれに1 を加えてisとする。｛001,011 ｝の 2つのビットパターン毎に id およびisを作り、x(1:n1,j2,id),x(1:n1,j2,is) を交換する。ステップＳ２０９においては、１次元目の並び替えはないが、id2,is2 次元の並び替えが必要なものである、ビットパターン［000,010,101,111]に対して、一つ前の処理と同じようにidを計算する。x(1:n1,j2,id),x(1:n1,j2,is) の交換をビットパターン毎に行なう。ステップＳ２１０において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ２１０における判断がＮｏの場合には、ステップＳ２０７に戻る。ステップＳ２１０における判断がＹｅｓの場合には、ステップＳ２１１に進む。

ステップＳ２１１では、nbitx=(n_b-4)/2、nx=2**nbitxを計算する。ステップＳ２１２において、ｉを１だけ増加し、ステップＳ２１３において、ｊを１だけ増加する。ステップＳ２１４において、ｊ＜ｉか否かを判断する。ステップＳ２１４の判断がＮｏの場合には、ステップＳ２１８に進む。ステップＳ２１４の判断がＹｅｓの場合には、ステップＳ２１５において、(i-1) を１ビット左シフトしたものに、(j-1) をビットリバースしたものを(nbitx+3)ビット左にシフトしたものを加えたもをidx とする。i とｊを入れ替えて同じ操作をしたものをisx とする。ステップＳ２１６において、０から16までのビットパターンを生成して、1,2,1 ビットに左から分けてB1,B2,B3とする。idx の2*nbit+4ビット目にB1を、nbit+3,nbit+2 ビット目にB2を、1 ビット目にB3を加えそれに1を加えてidとする。同様に、isx の2*nbit+3ビット目にB3を、nbit+3, nbit+2ビット目にB2を反転いたものを、1 ビット目にB1を加えて、それに1 を加えてisとする。0 から16までのビットパターンを生成する毎に、このように、idおよびisを作り、x(1:n1,j2,id),x(1:n1,j2,is) を入れ換える。ステップＳ２１７において、ｊを１だけ増加させ、ステップＳ２１４に戻る。

ステップＳ２１８では、ｉを１だけ増加し、ステップＳ２１９で、ｉ＞ｎ_ｘか否かを判断する。ステップＳ２１９の判断がＮｏの場合には、ステップＳ２１３に戻る。ステップＳ２１９の判断がＹｅｓの場合には、ステップＳ２２０に進む。ステップＳ２２０においては、ｉを１と設定する。ステップＳ２２１においては、(i-1) を１ビット左シフトしたものに、(i-1) をビットリバースしたものを(nbitx+3) ビット左にシフトしたものを加えたものをidx とする。isx はidx と同じ値に設定する。ステップＳ２２２において、[0001, 0010, 0011, 0101, 0111, 1011｝の6 つのビットパターンに対して、B1, B2, B3と 1,2,1ビットをおのおの取り出して、idx の2*nbit*3ビット目にB1を、nbit+2ビット目にB2を、1 ビット目にB3を加えそれに1 を加えてidとする。同様に、isx の2*nbit+4ビット目にB3を、nbit+3,nbit+2ビット目にB2を反転したものを、1 ビット目にB1を加えて、それに1 を加えてisとする。上記６つのビットパターン毎に、idおよびisを作り、x(1:n1,j2,id),x(1:n1,j2,is) を入れ換える。ステップＳ２２３において、ｉを１だけ増加し、ステップＳ２２４において、ｉ＞ｎ_ｘか否かを判断する。ステップＳ２２４の判断がＮｏの場合には、ステップＳ２２１に戻り、ステップＳ２２４の判断がＹｅｓの場合には、サブルーチンから抜け出る。

（付記１）
連続したデータ領域へのアクセスが高速化された、１次元目のインデックスをｉ、２次元目をｊ、３次元目をｋとしたとき、入力データｘ（ｉ、ｊ、ｋ）のｉの連続する方向に並んだデータを連続領域に格納する共有メモリ型スカラ並列計算機における３次元フーリエ変換処理方法において、
インデックスｉとｊからなる２次元配列データｘ（ｉ、ｊ、＊）（＊は、ｋの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｉが変化する方向への１次元目のフーリエ変換を、ｉを連続に変化させながら、ｊ方向の多重度分行い、
インデックスｊが変換する方向への２次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、
インデックスｋが変化する方向への３次元目のフーリエ変換を行なう
ことを特徴とする３次元フーリエ変換処理方法。
（付記２）
前記３次元目のフーリエ変換を行うステップは、
インデックスｉとｋからなる２次元配列データｘ（ｉ、＊、ｋ）（＊は、ｊの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｋが変化する方向への３次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行なう
ことを特徴とする付記１に記載の３次元フーリエ変換処理方法。
（付記３）
前記１次元目、２次元目、３次元目のそれぞれのフーリエ変換は、それぞれの次元を表すインデックスの取りうる値の数を基数に分解し、該基数を基礎とした小さなフーリエ変換と回転因子の積の繰り返し処理によって行なうことを特徴とする付記２に記載の３次元フーリエ変換処理方法。
（付記４）
前記小さなフーリエ変換は、キャッシュにロードされた２次元配列のうちから１行ずつ読み出し、データ長を基数に分解して、データを基数×整数の２次元配列とみなして処理して、作業域に転置した配列として結果を格納することを特徴とする付記３に記載の３次元フーリエ変換処理方法。
（付記５）
１次元目方向のビットリバースと２次元目方向のビットリバースを、ビットリバースを行なうべきデータがレジスタ上にある間に、同時に行なうことを特徴とする付記１に記載の３次元フーリエ変換処理方法。
（付記６）
前記同時に行われるビットリバースは、フーリエ変換されるべきデータの長さが、２の冪で表される場合に適用されることを特徴とする付記５に記載の３次元フーリエ変換処理方法。
（付記７）
前記フーリエ変換は、Cooley-Tukeyの方法で実行されることを特徴とする付記１または２に記載の３次元フーリエ変換処理方法。
（付記８）
連続したデータ領域へのアクセスが高速化された、１次元目のインデックスをｉ、２次元目をｊ、３次元目をｋとしたとき、入力データｘ（ｉ、ｊ、ｋ）のｉの連続する方向に並んだデータを連続領域に格納する共有メモリ型スカラ並列計算機における３次元フーリエ変換処理方法を該共有メモリ型スカラ並列計算機に実現させるプログラムにおいて、
インデックスｉとｊからなる２次元配列データｘ（ｉ、ｊ、＊）（＊は、ｋの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｉが変化する方向への１次元目のフーリエ変換を、ｉを連続に変化させながら、ｊ方向の多重度分行い、
インデックスｊが変換する方向への２次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、
インデックスｋが変化する方向への３次元目のフーリエ変換を行なう
ことを特徴とする３次元フーリエ変換処理方法を該共有メモリ型スカラ並列計算機に実現させるプログラム。
（付記９）
前記３次元目のフーリエ変換を行うステップは、
インデックスｉとｋからなる２次元配列データｘ（ｉ、＊、ｋ）（＊は、ｊの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｋが変化する方向への３次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行なう
ことを特徴とする付記８に記載の３次元フーリエ変換処理方法を該共有メモリ型スカラ並列計算機に実現させるプログラム。
（付記１０）
１次元目方向のビットリバースと２次元目方向のビットリバースを、ビットリバースを行なうべきデータがレジスタ上にある間に、同時に行なうことを特徴とする付記８に記載の３次元フーリエ変換処理方法を該共有メモリ型スカラ並列計算機に実現させるプログラム。
（付記１１）
連続したデータ領域へのアクセスが高速化された、１次元目のインデックスをｉ、２次元目をｊ、３次元目をｋとしたとき、入力データｘ（ｉ、ｊ、ｋ）のｉの連続する方向に並んだデータを連続領域に格納する共有メモリ型スカラ並列計算機における３次元フーリエ変換処理装置において、
インデックスｉとｊからなる２次元配列データｘ（ｉ、ｊ、＊）（＊は、ｋの取りうる任意の値）を単位として、キャッシュにロードするロード手段と、
インデックスｉが変化する方向への１次元目のフーリエ変換を、ｉを連続に変化させながら、ｊ方向の多重度分行う１次元目変換手段と、
インデックスｊが変換する方向への２次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行う２次元目変換手段と、
インデックスｋが変化する方向への３次元目のフーリエ変換を行なう３次元目変換手段と、
を備えることを特徴とする３次元フーリエ変換処理装置。
（付記１２）
前記３次元目変換手段は、
インデックスｉとｋからなる２次元配列データｘ（ｉ、＊、ｋ）（＊は、ｊの取りうる任意の値）を単位として、キャッシュにロードする手段と、
インデックスｋが変化する方向への３次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行なう手段と、
を備えることを特徴とする付記１１に記載の３次元フーリエ変換処理装置。
（付記１３）
１次元目方向のビットリバースと２次元目方向のビットリバースを、ビットリバースを行なうべきデータがレジスタ上にある間に、同時に行なうことを特徴とする付記１１に記載の３次元フーリエ変換処理装置。

本発明の実施形態の２次元面でのフーリエ変換の仕方を説明する図である。本発明の実施形態の３次元フーリエ変換の処理を説明するフローチャートである。図２における３次元目のフーリエ変換の処理を示すフローチャートである。本発明の実施形態に従った、アクセスパターンが連続で、かつ同じパターンで最後のループ並び替えを同時に行ない、かつ、データの先読みを行なう１次元フーリエ変換のフローチャートである。サブルーチンfft1の処理を示すフローチャートである。本発明の実施形態に従った、サブルーチンｆｆｔ２の処理を説明するフローチャートである。２次元変換でビットリバースの並び替えを同時に行なう方法を説明するフローチャートである。サブルーチンbitrb2dの処理を説明するフローチャート（その１）である。サブルーチンbitrb2dの処理を説明するフローチャート（その２）である。サブルーチンbitrb2dの処理を説明するフローチャート（その３）である。サブルーチンbitrb2dの処理を説明するフローチャート（その４）である。サブルーチンpermxyの処理を表すフローチャート（その１）である。サブルーチンpermxyの処理を表すフローチャート（その２）である。サブルーチンpermxyの処理を表すフローチャート（その３）である。サブルーチンpermxyの処理を表すフローチャート（その４）である。サブルーチンperm1dの処理を説明するフローチャート（その１）である。サブルーチンperm1dの処理を説明するフローチャート（その２）である。サブルーチンperm1dの処理を説明するフローチャート（その３）である。サブルーチンperm1dの処理を説明するフローチャート（その４）である。３次元目のフーリエ変換処理のフローチャートである。サブルーチンbitrb3dの処理を説明するフローチャート（その１）である。サブルーチンbitrb3dの処理を説明するフローチャート（その２）である。サブルーチンbitrb3dの処理を説明するフローチャート（その３）である。サブルーチンbitrb3dの処理を説明するフローチャート（その４）である。

符号の説明

Ｓ１０３次元目スレッド割り当てステップ
Ｓ１１２次元フーリエ変換ステップ
Ｓ１２２次元目スレッド割り当てステップ
Ｓ１３３次元目フーリエ変換ステップ

Claims

連続したデータ領域へのアクセスが高速化された、１次元目のインデックスをｉ、２次元目をｊ、３次元目をｋとしたとき、入力データｘ（ｉ、ｊ、ｋ）のｉの連続する方向に並んだデータを連続領域に格納する共有メモリ型スカラ並列計算機における３次元フーリエ変換処理方法において、
該共有メモリ型スカラ並列計算機は、
インデックスｉとｊからなる２次元配列データｘ（ｉ、ｊ、＊）（＊は、ｋの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｉが変化する方向への１次元目のフーリエ変換を、ｉを連続に変化させながら、ｊ方向の多重度分行い、
インデックスｊが変化する方向への２次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、
インデックスｋが変化する方向への３次元目のフーリエ変換を行い、
前記３次元目のフーリエ変換を行うステップは、
インデックスｉとｋからなる２次元配列データｘ（ｉ、＊、ｋ）（＊は、ｊの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｋが変化する方向への３次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、
１次元目方向のビットリバースと２次元目方向のビットリバースを、ビットリバースを行なうべきデータがレジスタ上にある間に、同時に行なう
ことを特徴とする３次元フーリエ変換処理方法。
前記１次元目、２次元目、３次元目のそれぞれのフーリエ変換は、それぞれの次元を表すインデックスの取りうる値の数を基数に分解し、該基数を基礎とした小さなフーリエ変換と回転因子の積の繰り返し処理によって行なうことを特徴とする請求項１に記載の３次元フーリエ変換処理方法。
前記小さなフーリエ変換は、キャッシュにロードされた２次元配列のうちから１行ずつ読み出し、データ長を基数に分解して、データを基数×整数の２次元配列とみなして処理して、作業域に転置した配列として結果を格納することを特徴とする請求項２に記載の３次元フーリエ変換処理方法。
前記同時に行われるビットリバースは、フーリエ変換されるべきデータの長さが、２の冪で表される場合に適用されることを特徴とする請求項１に記載の３次元フーリエ変換処理方法。
前記フーリエ変換は、Cooley-Tukeyの方法で実行されることを特徴とする請求項１に記載の３次元フーリエ変換処理方法。
連続したデータ領域へのアクセスが高速化された、１次元目のインデックスをｉ、２次元目をｊ、３次元目をｋとしたとき、入力データｘ（ｉ、ｊ、ｋ）のｉの連続する方向に並んだデータを連続領域に格納する共有メモリ型スカラ並列計算機における３次元フーリエ変換処理方法を該共有メモリ型スカラ並列計算機に実現させるプログラムにおいて、
インデックスｉとｊからなる２次元配列データｘ（ｉ、ｊ、＊）（＊は、ｋの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｉが変化する方向への１次元目のフーリエ変換を、ｉを連続に変化させながら、ｊ方向の多重度分行い、
インデックスｊが変化する方向への２次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、
インデックスｋが変化する方向への３次元目のフーリエ変換を行い、
前記３次元目のフーリエ変換を行うステップは、
インデックスｉとｋからなる２次元配列データｘ（ｉ、＊、ｋ）（＊は、ｊの取りうる任意の値）を単位として、キャッシュにロードし、
インデックスｋが変化する方向への３次元目のフーリエ変換を、ｉが連続して変化する方向への演算を優先して、ｉ方向の多重度分行ない、
１次元目方向のビットリバースと２次元目方向のビットリバースを、ビットリバースを行なうべきデータがレジスタ上にある間に、同時に行なう
とを特徴とする３次元フーリエ変換処理方法を該共有メモリ型スカラ並列計算機に実現させるプログラム。