memo.txt

# dev note for otfft

-------------------------------------------------------------------------------
  Memo
-------------------------------------------------------------------------------

2016-03-30

* ToDo: -Wunused-function の警告がうるさい。
  しかしヘッダファイル中の static 関数は必ずしも呼び出さなくて良い筈なので、
  この警告を出さないようにするべきである。一つの方法は単純に static にするのではなくて、

  static inline

  にする事である。もし inlining をしないという風にするのであれば、

  __attribute__ ((unused,noinline)) static

  という様に指定しないといけない。

* ToDo: 特定の環境で simd_malloc が適切に align されていない。

* ToDo: otfft_misc.h は fft 部分の実装でしか使わない物と、
  一般のユーザも使う物が混ざり合っている。分離して欲しい。

* ToDo: 名前空間の整理

* ToDo: 自動生成されているっぽいのをもっと何とか整理できないのか

-------------------------------------------------------------------------------
  ToDo
-------------------------------------------------------------------------------

2015-11-12

* 多次元FFTを並び替えなしで実行できるか:

  多次元FFTは各行各列のFFTの繰り返しである。
  しかし各列の変換に関してはデータが格納されている位置が飛び飛びである。
  これを既存のFFTのコードに渡す為にはデータを並べ替える必要がある。

  一方で元々FFTは複数の小さなFFTに分解して実行するという性質があることから、
  既存のFFTのコードを工夫して用いる事によって
  飛び飛びになっているデータに対しても変更せずに
  計算を実行できるのではないかという気がする。

  本当にそのようにできるか確かめる為には既存の実装について理解する必要がある。

  fwd0fft<n,s,eo>

    N = n1(global,target) * s(local,done)

  fwdcore<n,s>

    N = n1(global) * 4(target) * s(local)

  現在の実装を眺めてみると 1 FFT を 2 FFT に分割した際、
  メモリ的に、より non-local(s小) な FFT を先に実行してから
  local (s大) な FFT を実行している?

  % A できるだけローカルに計算を行うという可能性?
  %
  %   FFT では分割統治で計算を行う。
  %   通常は一番小さい単位の計算を全て完了してから、次の単位の計算を全て完了して、
  %   という具合に計算を進めていく。しかし、これだと全体のデータを何度も舐めることになる。
  %   それよりは「その単位の結果が必要になった時に、その単位が依存しているより小さな単位の計算を要求する」
  %   という形にした方が特に比較的小さい単位の変換においてキャッシュが働くのではないだろうか。
  %
  %   x 現在の実装で sequential なアクセスになっている。
  %     現在の実装を覗いてみたが途中状態のメモリ配置を工夫する事によって、
  %     どのステップでも十分に sequential なアクセスで計算できている。
  %
  %   x 無駄に複雑な制御。メモリアクセスの予測を阻害するのではないか。
  %     たしかに制御を複雑にすることでアクセスの範囲をできるだけ小さくする事はできる。
  %     しかし sequential なアクセスにはならないし寧ろメモリアクセスの予測を阻害する事になるのでは。
  %
  %   x 更にこの方法を用いると特に大きな単位の計算において
  %     アクセス範囲を絞る事ができない。
  %     寧ろ、この様な複雑な制御を採用していることにより、
  %     メモリ配置を工夫するという事ができなくなり、
  %     sequential なアクセスを工夫することが不可能である。
  %
  %   →この案は寧ろ遅くなる。全然駄目。却下。

  よく見ると fwdcore ではデータの並び替えを行って、
  次のステップでも同様に sequential なアクセスになる様に工夫をしている。
  この様な場合、一括多次元FFTのために fwdcore を使うのは得策でないかもしれない。

  B 多次元FFTを一つの fwdcore で計算出来ないか?

    x W はどうするのか? 巨大な物を用意するのか?

      →巨大な W をそのまま使っても問題は生じない。
        何しろ FFT の各 step (より低次元のFFT) に際しても
        同じ W を用いて計算ができているからである。

        具体的にはより低次元の FFT においては
        W は適当に間引かれて使用されることになるからである。

      ただ、やはり巨大な W を用意しなければならない
      ということに変わりはない。

      というか W の添字の倍率もテンプレート引数で
      指定できる様にすれば良いという気もする。
      或いは元より W として様々な物を用意しておいて、
      FFT の次数に応じて適切な W を中に渡すという事にするか。

      % →W の配列を工夫してどのサイズでも同様に使える様にする?
      %   つまり  Wx = [ 0, 1/2, 1/4, 3/4, 1/8, 5/8, 3/8, 7/8, ... ]
      %   として例えば N = 16 の時は Wx の初めの 16 要素を取ってきて使えば良い。
      %   通常の W の要素は W[i] = Wx[ ビット逆順(i) ] で取得する事ができる。
      %
      %   と思ったがビット逆順の操作はかなり面倒なようだ。
      %   調べたがすっきりと実行する方法・またはCPU命令は存在しないようである。
      %   やはり需要のない変な機能ということのようだ。
      %
      % このような事ならばやはり W は普通に複数保持したほうが良いだろう。
      % 全て保持したとしても大きさは2倍になるだけである。

    x x -> y-> x の様にして交互に計算を実行する事になっているが、
      全体を一気に変換しようとすると buffer y として巨大な物を用意しなければならない。
      しかしながら sequential なアクセスにする為なのだと考えればこれは仕方ない様にも思う。

    o メモリアクセスはシーケンシャルに全体を走る物になる。
      気分よく計算する事ができるだろう。

    - 3次元になるともっと複雑な制御をしなければならないのではないか?

      3次元の一番内側のFFTは通常通りにすれば良い。
      一番外側のFFTも2次元FFTと思って計算すれば良い。
      真ん中にある物はどうするのか?→と思ったがこれも簡単だ。
      一番外側については普通のループにして、内側の2段に関して2次元列FFTを実行をすれば良い。
      これで十分に seq なアクセスにする事ができる。

    ? s の大きさをコンパイル時に決定するのではなく、動的に指定できる様にできないか?

      以下の部分の +N0, etc の部分が変数になるということである。
      const ymm a = getpz2(xq_sp+N0);
      const ymm b = getpz2(xq_sp+N1);
      const ymm c = getpz2(xq_sp+N2);
      const ymm d = getpz2(xq_sp+N3);
      もしかするとこの部分を定数にする事によって
      x86 命令の addressing 部分に埋め込まれる事を想定している?
      →逆アセンブルしてみるとやはりそうなっていた。

      そうだとするとこの部分を変数にすることは余り得策ではない。
      実際のアドレス計算は [xq_sp + 16 * N1] という事になると予想される。
      sizeof(complex_t) == 16 なので。しかし、16 という倍率の index 指定はない。
      従ってこの部分は complex_t によるアドレス計算ではなくて
      double* によるアドレス計算に変えて取り扱う必要があるかも。
      しかしそうだとしても N1, N2, N3 をそれぞれレジスタに置いておく必要が出るのでやはり微妙?
      →でも逆アセンブルしてみたところレジスタはかなり空いている様にも思われる。
        レジスタ退避はループの外側に成るだろうからそんなに問題もなさそうだし。

      うーん命令列に埋め込まれている即値を動的に書き換える事ができればとても楽なのだけれど。

      というかアドレスを計算するために命令がいくつか増えるぐらいならば
      最終的な速度に余り影響を与えないのではないかという可能性もある?
      特に _mm256 系統の命令とアドレスの計算の命令が混ざっている場合、
      適当に並列に実行されたりしないのだろうか。でもここは最も内側のループだし…。

      或いは初めから xq0, xq1, xq2, xq3 という4つのポインタを用意するのも手かも。
      complex_vector xq0_sp = xq + N0 + sp;
      complex_vector xq1_sp = xq + N1 + sp;
      complex_vector xq2_sp = xq + N2 + sp;
      complex_vector xq3_sp = xq + N3 + sp;
      しかしこれだとレジスタを 4 つも専有する事になり駄目な気がする。
      yq の方もある事を考慮にいれると合計で 8 のレジスタが必要だ。
      やはり現在の様に定数として N1 を埋め込んでいるのは速度面で重要な気がする。

  結論としては技術的には可能であるが fwdcore に手をいれる必要がある。
  実装も面倒であるし otfft の製作者に連絡が取れない現状で version 管理も面倒である。
  これらのことから現状では愚直にデータを並び替えてから1次元FFTを実行する方法で良いと思う。
  一応多次元FFTを並び替え無しで行う方法に関してまとめておく:

  - fwdcore は W として次数 n に対応する物を受け取る事にする。
    (使用する際は今までの W の添字を s で割ったものを用いれば良い)
  - fwdcore の引数の一部を動的にした場合最終的な performance に影響が出るかもしれない。
    特に otfft の優位な点が丁度この部分から来ているのではないかという可能性…。

2015-11-11

* 更新されている。

  またもや大幅な書き換えになっている。
  ソースコードを幾らか読んでみたが冗長で繰り返しが多い。
  どうも何らかの別のプログラムから自動生成されている様に思われる。
  このままだと差分が本質的でなく使いづらい。

  もし original の otfft の更新が止まっているのであれば
  勝手に冗長な物を消して自動生成するプログラムを自分の手で書くのであるが、
  リアルタイムに更新されている物に対して大幅な書き換えを実行するのは困難である。
  これは最早作者に相談した方が良いような気がする…が連絡先が何処にも書かれていない。

  取り敢えず要望について整理しておいた方が良いように思う。

  1 自動生成のプログラムは公開不能か?

    或いは、公開可能な別の方式で自動生成する様にできないか。

  2 構造の変更?

    比較的短い同じ長さの FFT を何度も実行する時、
    長さによる switch (table jump) を何度も実行する事になるが、
    これを何とかループの外側に持って来られないか。

    これは多次元FFTを実行する際に絶対不可欠である。
    書き換えとしては複数の FFT (並列) を同時に実行できる様にすれば良い。

    或いはこちらで勝手に書き換えても良いか。

  3 実FFT

    実→複素で、複素数の後半も生成する様になっているが、
    後半は前半と等価な情報を持っているので、実際の数値計算で使用する事は滅多にない。
    メモリサイズを節約したいので後半を使用しないような version/flag を用意することはできないか。

  4 実FFT (2並列)

    多次元FFTの場合には同じ長さの実FFTを何度も実行する事になる。
    "1実FFTを半分の長さの複素FFTで実行するもの" だけでなく
    "2実FFTを同じ長さの複素FFTで実行するもの" も欲しい。
    (前者を二回繰り返した方が performance が高いのなら無意味かもしれないが)

  5 提案?

    勘違いしているだけで実はナンセンスな事を書いているかもしれないが。

    複数の方法のうち最速の物を選択する事になっているが、
    どの方法も一番長い変換から短い変換まで同一の方法を用いると想定している。
    複数の方法を混合した物も試してみて良いのではないか。
    例えば FFT(2^m) から FFT(2^{m-2}) を呼び出す時に既に分かっている最速の方法を呼ぶなど。


-------------------------------------------------------------------------------
  Done
-------------------------------------------------------------------------------

2016-09-20

* 実は otfft の作者への連絡先が v8.1 (20160322) から renraku.png
  という画像に含まれていたという事が判明した。

  内容は少々攻撃的な文章である。
  もしかするとこの身勝手なことを書いているメモを
  見られてだいぶ警戒させてしまったのかもしれない。

  さて、当初の目的が一体どんなことだったか忘れてしまったので、
  また調べてみる。コードの変更についての提案である。
  それとは別に、

  1 ライセンスについての詳細
  2 #include <pmmintrin.h> が namespace の中で行われていること
  3 complex_t が std::complex<double> に変換できないこと
  4 OPENMP を使うか使わないかのスイッチ

  などを聞きたかったが、時間が立つうちに (1), (2), (4) は対応がなされた。
  (3) については単純なことだし明らかなミスに思われるので対応をお願いする。
  (1) に関連して勝手に github で公開していることについてどう考えているか確認できれば幸いである。
  他に 2016-03-30 にある様な細かいことがあるが、
  最新版でもこれらの問題が残っているかどうかは確認していないので質問は控える。


2015-12-26

* 久しぶりに見たら丁度更新された所だった

  追加・削除されたファイルの一覧

    +LICENSE.txt
    +otfft/LICENSE.txt
    +otfft/msleep.h
    +otfft/otfft_avxdif16.h
    +otfft/otfft_avxdif16omp.h
    -otfft/otfft_avxdifx.h
    +otfft/otfft_avxdit16.h
    +otfft/otfft_avxdit16omp.h
    +otfft/README.txt

  変更を追跡していないファイル

    Makefile
    Makefile.linux
    Makefile.osx
    Makefile.windows
    src/Makefile
    src/Makefile.linux
    src/Makefile.osx
    src/Makefile.windows

  ./update.sh modify-source で対応できない物 (手動で変更するべき物)

    | diff -bwur html/d/5.4/src/otfft/otfft_misc.h src/otfft/otfft_misc.h
    | --- html/d/5.4/src/otfft/otfft_misc.h	2015-12-27 08:26:23.918671452 +0900
    | +++ src/otfft/otfft_misc.h	2015-11-14 00:21:28.886062221 +0900
    | @@ -94,7 +94,7 @@
    |      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<double>& z) : Re(z.real()), Im(z.imag()) {}
    | -    operator std::complex<double>() { return std::complex<double>(Re, Im); }
    | +    operator std::complex<double>() const { return std::complex<double>(Re, Im); }
    |
    |      complex_t& operator+=(const complex_t& z) {
    |          Re += z.Re;
    | @@ -322,10 +322,14 @@
    |
    |  #if defined(__SSE3__) && defined(USE_INTRINSIC)
    |
    | +} // namespace OTFFT_MISC
    | +
    |  extern "C" {
    |  #include <pmmintrin.h>
    |  }
    |
    | +namespace OTFFT_MISC {
    | +
    |  static inline xmm haddpz(const xmm ab, const xmm xy) force_inline;
    |  static inline xmm haddpz(const xmm ab, const xmm xy)
    |  {


2015-11-14

* original otfft.h に対する修正の一覧

  以下の様にしないと pmmintrin.h の中身が OTFFT_MISC に展開されてしまうので、
  他のヘッダで pmmintrin.h を使う物があった時にコンパイルできない。

  | +} // namespace OTFFT_MISC
  |
  |  extern "C" {
  |  #include <pmmintrin.h>
  |  }
  |
  | +namespace OTFFT_MISC {

  以下のようにしないと complex_t const& や const_complex_vector
  で受け取った値を std::complex<double> にできない。

  | -operator std::complex<double>() { return std::complex<double>(Re, Im); }
  | +operator std::complex<double>() const { return std::complex<double>(Re, Im); }

2015-10-24

* 改めてページを確認した所、更に更新がされている様だったので再度取得を試みる。

  $ mkdir html2; cd html2; wget -r http://www.moon.sannet.ne.jp/okahisa/stockham/stockham.html

  ソースコードの差分を見たら、check/cpp_fftw3.h の
  FFTW_MEASURE が FFTW_ESTIMATE に変更されただけであった。適用する。

* 更新otfft-4.0 → otfft-5.3

  元の otfft の version が上がっている様なので更新を行う。

  ファイルの増減は以下の通りである。

  --- a.find       2015-10-24 19:25:02.636561285 +0900
  +++ b.find       2015-10-24 19:25:02.643560487 +0900
  @@ -1,7 +1,6 @@

   /bstcheck.cpp
   /cpp_fftw3.h
  -/ctavx_fft.h
   /dctcheck.cpp
   /fftbench1.cpp
   /fftbench1.txt
  @@ -17,16 +16,16 @@
   /otfft/Makefile
   /otfft/otfft.cpp
   /otfft/otfft.h
  -/otfft/otfft_ctavx.h
  -/otfft/otfft_ctavxn.h
  -/otfft/otfft_difavx.h
  -/otfft/otfft_difavx8.h
  -/otfft/otfft_difavx8n.h
  -/otfft/otfft_difavxn.h
  -/otfft/otfft_ditavx.h
  -/otfft/otfft_ditavx8.h
  -/otfft/otfft_ditavx8n.h
  -/otfft/otfft_ditavxn.h
  +/otfft/otfft_avxdif4.h
  +/otfft/otfft_avxdif4omp.h
  +/otfft/otfft_avxdif8.h
  +/otfft/otfft_avxdif8omp.h
  +/otfft/otfft_avxdifx.h
  +/otfft/otfft_avxdit4.h
  +/otfft/otfft_avxdit4omp.h
  +/otfft/otfft_avxdit8.h
  +/otfft/otfft_avxdit8omp.h
  +/otfft/otfft_eightstep.h
   /otfft/otfft_fwd.h
   /otfft/otfft_fwd0.h
   /otfft/otfft_inv.h
  @@ -34,10 +33,14 @@
   /otfft/otfft_misc.h
   /otfft/otfft_setup.h
   /otfft/otfft_sixstep.h
  -/otfft/otfft_sixstepn.h
  +/otfft/otfft_sixstep0r.h
  +/otfft/otfft_sixstep0s.h
  +/otfft/otfft_sixstepnr.h
  +/otfft/otfft_sixstepns.h
   /otfft/stopwatch.h
  +/otfft/towin.sh
   /otfft/version.sh
   /rfftcheck.cpp
   /simple_fft.h
  +/towin.sh
   /unify.sh

  ot_fft_dif* から otfft_avxdif* に名称変更が行われている様である。
  一方でこちらで勝手に修正した項目も各ファイルに存在している。
  取り敢えず手許で行われた変更について確認を行う。

  ヘッダ名の変更
  -#include "otfft_misc.h"
  +#include "otfft/otfft_misc.h"
  と、#pragma omp を #ifdef _OPENMP で囲むという事だけの様だ。

  他に、以下のヘッダ名の変更がある。
  -#include "stopwatch.h"
  +#include "otfft/stopwatch.h"
  -#include "otfft.h"
  +#include "otfft/otfft.h"

  1 src ディレクトリ

    1.1 先ずファイルを新しいものに全て交換する。
      otfft.h stopwatch.h otfft_misc.h は otfft ディレクトリ以下に移動する。

    1.2 次に #ifdef _OPENMP の付加を以下のコマンドで実行する。
      $ sed -i 's/^[[:space:]]*#pragma[[:space:]]\{1,\}omp[[:space:]]\{1,\}.*$/#ifdef _OPENMP\n&\n#endif/g' *.h *.cpp otfft/*.h

    1.3 ヘッダファイル名の変更は以下で行う。
      $ refact '(otfft|otfft_misc|stopwatch)\.h' 'otfft/&' *.h *.cpp

    1.4 b/otfft/towin.sh をコピーする。

  2 check ディレクトリ以下についてもファイルを更新する。
    $ cd check; cp ../html/diff/b/*.{cpp,h} ./; rm ctavx_fft.h

    2.1 USE_FFTW_THREADS で囲む
      $ refact '^.*\bfftw_[a-z_]+threads\b.*$' '#ifdef USE_FFTW_THREADS\n&\n#endif'
      2行連続で fft_...threads の呼出をしている箇所が cpp_fftw.h にあったので微修正する。

      ※check ディレクトリ以下の *.cpp で otfft/otfft.h, otfft/stopwatch.h 等を include する時に、
        #include "..." ではなくて #include <...> を使う様にしていたが面倒なので、
        #include "..." のまま保持する事にした。

    2.2 /check/fftbench?.txt は repository から除外する。新しい物をコピーする。
      $ git rm fftbench?.txt

    2.3 b/towin.sh をコピーする。

  3 src/version.sh 差分を適用