OTFFT の使い方

【OTFFT の使い方】

　さて、いったい何人の方が使ってくださるか分かりませんが、筆者 OK おじさんの提案する OTFFT の使い方について説明しておきましょう。

【OTFFT とは】

　OTFFT とは、Stockham FFT のアルゴリズムと SSE2/AVX/AVX2 などを組み合わせた高速な FFT ライブラリです。 C++ の template によるメタプログラミング技法も用いています。 OTFFT は、OK おじさん Template FFT の略です。 OK おじさんは、OTFFT の開発者のニックネームです。

【使用上の注意】

　OTFFT は、Clang でコンパイルする場合、コンパイルする環境に最適化されます。つまり、AVX が有効な環境でコンパイルすれば、AVX を使うようになります。このバイナリを AVX をサポートしない環境で動かせば、当然、例外を起こして落ちます。そのため、コンパイルしたバイナリを、広く一般に配布するような利用形態には、OTFFT は向きません。OTFFT は、数値計算のプログラムをソースからコンパイルするような利用形態を想定しています。

【必要なファイル】

　ダウンロードのページから otfft-11.6.tar.gz(もしくは otfft-11.5v.zip) をダウンロードし展開すると、いくつかのファイルができますが、 OTFFT を使うのに必要なファイルはその中の otfft フォルダとその中身です。このうち、

    otfft_gen_new.h
    otfft_gen_setup.h
    otfft_gen_fwd.h
    otfft_gen_fwd0.h
    otfft_gen_fwdu.h
    otfft_gen_inv.h
    otfft_gen_invu.h
    otfft_gen_invn.h
    otfft_gen_delete.h

は自動生成されるファイルで、自分の環境で最大の効果を発揮したければ ffttune を使って生成し直してやる必要があります。バージョン 10.2 までは名前に _gen_ が付きません。 Clang でコンパイルする場合、otfft フォルダで以下のように実行し、 ffttune コマンドを作ります。 Makefile は自分の環境に合わせていじってください。

    make ffttune

Mac の場合、デフォルトの Makefile は、 MacPorts の clang-3.8 以上がインストールされていると想定しています。 Visual Studio の場合、otfft フォルダで以下のように実行します。

    nmake ffttune

C++ の template によるメタプログラミング技法を用いているので、コンパイルにはかなりの時間がかかります。次に

    ./ffttune

あるいは Windows の場合、

    ffttune

と ffttune コマンドを実行します。すると先に示したファイルの自動生成が始まります。しばらく時間がかかるのでお待ちください。このとき、Core i7(4 コア 8 スレッド) の環境の場合、環境変数 OMP_NUM_THREADS を 7 に設定する必要があるかもしれません。一部の Linux では、8 スレッドで動かすと謎の速度低下を起こします。

　OTFFT は内部で幾つかのアルゴリズムを使い分けていますが、どのアルゴリズムがどのサイズで一番速いか、実際に計測して決定します。しかし、パソコンは FFT だけを黙々と実行しているわけではなく、裏で色々なプロセスが動いています。計測には常に突発的なゆれが混入します。そのため、１回の計測で最適な組み合わせが求まるとは限りません。納得いかない場合は何度か ffttune を実行し直してください。

　自動生成が終わったら、次は otfft.o あるいは otfft.obj ファイルを生成します。 Clang の場合は、otfft フォルダで次のコマンドを実行します。

    make otfft.o

Visual Studio の場合は、otfft フォルダで次のコマンドを実行します。

    nmake otfft.obj

　C++ の template によるメタプログラミング技法を用いているので、コンパイルには非常に長い時間がかかります。辛抱強くお待ちください。れで準備は終わりです。

【ご自分のバイナリのコンパイル方法】

　otfft.o あるいは otfft.obj の生成が終われば、あとは otfft/otfft.h をインクルードし、コンパイラのインクルードパスに otfft フォルダのあるフォルダを加えて、 otfft.o あるいは otfft.obj ファイルを目的のバイナリにリンクしてやれば OTFFT は使えます。

　Clang でコンパイルするには、 -lomp や -liomp5 など OpenMP のランタイムをリンクする必要があります。例えば Mac の場合、カレントフォルダに otfft フォルダがあるなら以下のようにコンパイルします。

    clang++-mp-8.0 -c -O3 hello.cpp
    clang++-mp-8.0 hello.o otfft/otfft.o -L/opt/local/lib/libomp -lomp -o hello

　clang++-mp-8.0 は MacPorts の clang++ コンパイラです。 Visual C++ の場合、 AVX を使うには、/arch:AVX を指定してください。 AVX2 も使う場合は /arch:AVX2 も指定してください。 AVX-512 も使う場合は /arch:AVX512 も指定してください。 OpenMP も使うには /openmp を指定してください。お分かりとは思いますが、otfft.obj のコンパイル時点で、これらのオプションを付けておかないと有効にはなりません。必要なら OTFFT の Makefile に記載してください (OTFFT のコンパイルが終われば、これらのオプションは必要ないかもしれません)。 Visual Studio 用のデフォルトの Makefile は、 AVX2 と OpenMP を使うよう指定されています。

【シングルスレッドモード】

　OTFFT はデフォルトでは OpenMP によるマルチスレッドで動作します。しかし、中にはマルチスレッドの実装は、自前の仕組みを使いたい人もいるかもしれません。そこで、シングルスレッドで動作するモードも用意しました。 OTFFT は DO_SINGLE_THREAD マクロを定義してコンパイルすると、シングルスレッドモードになります。シングルスレッドモードにしたい場合は、 otfft/otfft_misc.h で定義して、ffttune コマンドの生成から始めてください。

【アンアラインモード】

　OTFFT はデフォルトでは、変換する系列を保存するメモリにアラインされたメモリを要求します。しかし、既存のプログラムに OTFFT を組み込もうとすると、この制限のせいで大幅な書き換えが必要になることもあるでしょう。そこで、アラインされていないメモリでも動くモードを用意しました。OTFFT は、 USE_UNALIGNED_MEMORY マクロを定義してコンパイルすると、アンアラインモードになります。もちろん、アラインモードに比べると速度は低下します。アンアラインモードにしたい場合は、 otfft/otfft_misc.h で定義して、ffttune コマンドの生成から始めてください。

　以下で、サポートされる変換を説明します。変換する系列長 N は 2^30 より小さい必要がありますが、2^30 より小さくても、計算機のメモリが足りなくなれば変換できません。

【複素離散フーリエ変換】

　系列長 N の複素離散フーリエ変換を実行するには、

    #include "otfft/otfft.h"
    using OTFFT::complex_t;
    using OTFFT::simd_malloc;
    using OTFFT::simd_free;

    void f(int N)
    {
        complex_t* x = (complex_t*) simd_malloc(N*sizeof(complex_t));
        // 何かする...
        OTFFT::FFT fft(N); // FFT オブジェクトの作成
        fft.fwd(x);        // 複素離散フーリエ変換を実行。x が入力かつ出力
        // 何かする...
        simd_free(x);
    }

のようにします。ただし、 N の素因数が小さな素数にならない場合は、計算量が O(N^2) になります。

complex_t は

    struct complex_t
    {
        double Re, Im;

        complex_t() : Re(0), Im(0) {}
        complex_t(const double& x) : Re(x), Im(0) {}
        complex_t(const double& x, const double& y) : Re(x), Im(y) {}
        complex_t(const std::complex& z) : Re(z.real()), Im(z.imag()) {}
        operator std::complex() const { return std::complex(Re, Im); }

        // その他のメンバ関数...
    };

のように定義されています。

　フーリエ変換(fwd)には係数 1/N が掛かっています。もしそれが不都合なら係数の掛かっていない(fwd0)もあります。以下のようなメンバ関数が用意されています。

    fwd(x)  離散フーリエ変換(1/N による正規化付き)
    fwd0(x) 離散フーリエ変換(正規化無し)
    fwdu(x) 離散フーリエ変換(ユニタリ変換。つまり sqrt(1/N) で正規化)
    fwdn(x) 離散フーリエ変換(1/N による正規化付き)

    inv(x)  逆離散フーリエ変換(正規化無し)
    inv0(x) 逆離散フーリエ変換(正規化無し)
    invu(x) 逆離散フーリエ変換(ユニタリ変換)
    invn(x) 逆離散フーリエ変換(1/N による正規化付き)

FFT の系列長を変更したい場合は、

    fft.setup(2 * N);

のようにします。

　OTFFT は Stockham のアルゴリズムを使っていますので、実行には入力系列と同じ系列長の作業領域が必要です。普通は内部でその領域を確保しますが、マルチスレッドのプログラムを組む場合など、それだと都合が悪い場合があります。そんな時は外部から作業領域を渡すバージョンもあります。次のように使います。

    #include "otfft/otfft.h"
    using OTFFT::complex_t;
    using OTFFT::simd_malloc;
    using OTFFT::simd_free;

    void f(int N)
    {
        complex_t* x = (complex_t*) simd_malloc(N*sizeof(complex_t));
        complex_t* y = (complex_t*) simd_malloc(N*sizeof(complex_t));
        // 何かする...
        OTFFT::FFT0 fft(N);
        fft.fwd(x, y); // x が入力かつ出力、y が作業領域
        // 何かする...
        fft.inv(x, y); // x が入力かつ出力、y が作業領域
        // 何かする...
        simd_free(y);
        simd_free(x);
    }

　OTFFT::FFT だったところが OTFFT::FFT0 になっていることに注意してください。

【実離散フーリエ変換】

　実離散フーリエ変換では、系列長 N は偶数である必要があります。系列長 N の実離散フーリエ変換を実行するには以下のようにします。

    #include "otfft/otfft.h"
    using OTFFT::complex_t;
    using OTFFT::simd_malloc;
    using OTFFT::simd_free;

    void f(int N)
    {
        double*    x = (double*)    simd_malloc(N*sizeof(double));
        complex_t* y = (complex_t*) simd_malloc(N*sizeof(complex_t));
        // 何かする...
        OTFFT::RFFT rfft(N);
        rfft.fwd(x, y); // 実離散フーリエ変換を実行。x が入力、y が出力
        // 何かする...
        simd_free(y);
        simd_free(x);
    }

ただし、N の素因数が小さな素数にならない場合は、計算量が O(N^2) になります。

　実離散フーリエ変換には以下のようなメンバ関数が用意されています。

    fwd(x, y)  実離散フーリエ変換(1/N による正規化付き) x:入力、y:出力
    fwd0(x, y) 実離散フーリエ変換(正規化無し)           x:入力、y:出力
    fwdu(x, y) 実離散フーリエ変換(ユニタリ変換)         x:入力、y:出力
    fwdn(x, y) 実離散フーリエ変換(1/N による正規化付き) x:入力、y:出力

    inv(y, x)  逆実離散フーリエ変換(正規化無し)           y:入力、x:出力
    inv0(y, x) 逆実離散フーリエ変換(正規化無し)           y:入力、x:出力
    invu(y, x) 逆実離散フーリエ変換(ユニタリ変換)         y:入力、x:出力
    invn(y, x) 逆実離散フーリエ変換(1/N による正規化付き) y:入力、x:出力

　逆実離散フーリエ変換は複素系列 y を受け取って、実系列 x を返します。y が y[N-k] == conj(y[k]) でないと正しい結果を返しません。また、逆実離散フーリエ変換は入力 y を破壊します。保存しておきたい場合は、コピーを取っておく必要があります。内部的には系列長 N/2 の複素離散フーリエ変換で実装されています。

　実離散フーリエ変換は、作業領域として出力系列を使います。そして逆実離散フーリエ変換は、作業領域として入力系列を使います。マルチスレッドでプログラムする場合も、それぞれのスレッド専用の入力と出力を与えてやれば OK です。

【離散コサイン変換(DCT-II)】

　離散コサイン変換では、系列長 N は偶数である必要があります。系列長 N の離散コサイン変換(DCT-II)を実行するには以下のようにします。

    #include "otfft/otfft.h"
    using OTFFT::complex_t;
    using OTFFT::simd_malloc;
    using OTFFT::simd_free;

    void f(int N)
    {
        double* x = (double*) simd_malloc(N*sizeof(double));
        // 何かする...
        OTFFT::DCT dct(N);
        dct.fwd(x); // DCT-II を実行する。x が入力かつ出力
        // 何かする...
        simd_free(x);
    }

　ただし、N の素因数が小さな素数にならない場合は、計算量が O(N^2) になります。

　離散コサイン変換には以下のようなメンバ関数が用意されています。

    fwd(x)  離散コサイン変換(1/N による正規化付き)
    fwd0(x) 離散コサイン変換(正規化無し)
    fwdn(x) 離散コサイン変換(1/N による正規化付き)

    inv(x)  逆離散コサイン変換(正規化無し)
    inv0(x) 逆離散コサイン変換(正規化無し)
    invn(x) 逆離散コサイン変換(1/N による正規化付き)

　離散コサイン変換は DCT-II を採用しています。ただし、直交化はしていません。内部的には系列長 N/2 の複素離散フーリエ変換で実装されています。マルチスレッドで使う場合、作業領域を外部から与えるバージョンもあります。以下のように使います。

    #include "otfft/otfft.h"
    using OTFFT::complex_t;
    using OTFFT::simd_malloc;
    using OTFFT::simd_free;

    void f(int N)
    {
        double*    x = (double*)    simd_malloc(N*sizeof(double));
        double*    y = (double*)    simd_malloc(N*sizeof(double));
        complex_t* z = (complex_t*) simd_malloc(N*sizeof(complex_t));
        // 何かする...
        OTFFT::DCT0 dct(N);
        dct.fwd(x, y, z); // DCT を実行する。x が入力かつ出力、y,z が作業領域
        // 何かする...
        dct.inv(x, y, z); // 逆DCT を実行する。x が入力かつ出力、y,z が作業領域
        // 何かする...
        simd_free(z);
        simd_free(y);
        simd_free(x);
    }

　OTFFT::DCT だったところが OTFFT::DCT0 になっていることに注意してください。

【Bluestein's FFT】

　系列長 N の Bluestein's FFT(任意の系列長の FFT) を実行するには以下のようにします。

    #include "otfft/otfft.h"
    using OTFFT::complex_t;
    using OTFFT::simd_malloc;
    using OTFFT::simd_free;

    void f(int N)
    {
        complex_t* x = (complex_t*) simd_malloc(N*sizeof(complex_t));
        // 何かする...
        OTFFT::Bluestein bst(N); // N は 2^28 以下の任意の自然数。
        bst.fwd(x); // Bluestein's FFT を実行する。x が入力かつ出力
        // 何かする...
        simd_free(x);
    }

　Bluestein's FFT では、離散フーリエ変換の系列長が小さな素因数に分解できる必要はありません。例えば系列長が大きな素数でも計算量は O(N log N) になります。他の変換では、N が小さな素因数に分解できない場合は計算量が O(N^2) になります。指定出来る系列長の上限は 2^28 です。もちろん、メモリが足りなくなれば、それより小さい系列長でも変換できません。マルチスレッド用のメンバ関数は用意されていません。どうしてもマルチスレッドで使いたい場合は、動作確認はしていませんが、スレッドの数だけ Bluestein オブジェクトを生成し、各スレッドでそれらを使い分ければ大丈夫と思います。当然、メモリをゴージャスに使ってしまいます。

　Bluestein's FFT には以下のようなメンバ関数が用意されています。

    fwd(x)  離散フーリエ変換(1/N による正規化付き)
    fwd0(x) 離散フーリエ変換(正規化無し)
    fwdu(x) 離散フーリエ変換(ユニタリ変換)
    fwdn(x) 離散フーリエ変換(1/N による正規化付き)

    inv(x)  逆離散フーリエ変換(正規化無し)
    inv0(x) 逆離散フーリエ変換(正規化無し)
    invu(x) 逆離散フーリエ変換(ユニタリ変換)
    invn(x) 逆離散フーリエ変換(1/N による正規化付き)

【ベンチマーク】

　最後に、OTFFT にはベンチマークプログラムが付属しています。しかし、このベンチマークには FFTW3 と大浦さんのプログラムを用いているため、単独では動作しません。まずは FFTW3 をインストールし、大浦さんのページから大浦さんの FFT ライブラリをダウンロードしてください。そしてその中の fftsg.c というファイルを fftbench1.cpp と同じフォルダに入れてください。その上で以下のように fftbench1 コマンドを make します。

    make fftbench1

Visual Studio の場合 nmake します。

    nmake fftbench1

以下のように実行してください。

    ./fftbench1

あるいは Windows の場合、

    fftbench1

とします。

【追記：ベンチマークの変更点】

　FFTW の性能を十分に発揮したベンチマークを行うには、FFTW_MEASURE モードでベンチマークする必要があります。その場合、

    make fftbench5

を実行してから

    ./fftbench5

を実行してください。もちろん、make するには FFTW をインストールする必要があります。また、大浦さんの FFT パッケージから fftsg.c を持ってくる必要もあります。FFTW と OTFFT の比較のみで十分という方は

    make ibench wisdom="wisdom_name"

を実行してください。ただし、後述する wisdom ファイルは自分で作成する必要があります。

　fftbench5 の実行には FFTW の wisdom ファイルが必要です。make するときに自動的に作成されますが、とてつもなく長い時間がかかりますので覚悟してから実行してください。実行したらお茶でも飲んでお待ちください。たぶん１時間３０分ほどかかります。

　自分で wisdom ファイルを作成したいという方は、

    bash makewisdom.sh wisdom_name

を実行してください。「wisdom_name」は作成される wisdom ファイルのファイル名です。wisdom ファイルのファイル名を指定して fftbench5 を実行するには、

    ./fftbench5 1 22 wisdom_name

を実行します。22 を 24 にすれば、2^24 のサイズまでベンチマークします。