/Users/mourrain/Devel/mmx/algebramix/include/algebramix/fft_blocks.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_blocks.hpp
* DESCRIPTION: Cache friendly FFT
* COPYRIGHT  : (C) 2008  Gregoire Lecerf
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX__FFT_BLOCKS__HPP
#define __MMX__FFT_BLOCKS__HPP
#include <basix/vector.hpp>
#include <basix/threads.hpp>
#include <algebramix/fft_naive.hpp>

namespace mmx {

/******************************************************************************
* Multi-block FFT transformation
******************************************************************************/

template<typename C, typename FFTER= fft_naive_transformer<C>,
         nat log2_outer_block_size= 9,
         nat log2_block_number= 7, // number of outer blocks
                                   // to be handled in sequence
         nat log2_inner_block_size= 11,
         nat threshold= 12>
class fft_blocks_transformer {
public:
  typedef implementation<vector_linear,vector_naive> NVec;
  typedef typename FFTER::R R;
  typedef typename R::U U;
  typedef typename R::S S;

  FFTER* ffter;
  nat    depth;
  nat    len;
  U*     roots;

public:
  inline fft_blocks_transformer (nat n):
    ffter (new FFTER (n)),
    depth (ffter->depth), len (ffter->len), roots (ffter->roots) {}

  inline ~fft_blocks_transformer () { delete ffter; }

  static inline void
  delocate (C* c1, nat s1, C* c0, nat s0, nat n) {
    for (nat i= 0; i < s1; i++, c0 += s0) {
      C* cc0= c0, * cc1= c1 + i;
      for (nat k= 0; k < n; k++, cc0++, cc1+= s1)
        *cc1= *cc0; } }

  static inline void
  relocate (C* c1, nat s1, C* c0, nat s0, nat n) {
    for (nat i= 0; i < s1; i++, c0 += s0) {
      C* cc0= c0, * cc1= c1 + i;
      for (nat k= 0; k < n; k++, cc0++, cc1+= s1)
        *cc0= *cc1; } }

  inline void
  dfft (C* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    if (steps <= max (threshold, log2_inner_block_size))
      ffter->dfft (c, stride, shift, steps, step1, step2);
    else {
      nat temp_len= default_aligned_size<C>
        ((nat) 1<<(log2_outer_block_size+log2_block_number));
      C* c0, * temp= mmx_new<C> (temp_len);
      nat step, block, blocks, todo, k_beg, k_end;

      for (step= step1;
           step < min (steps - log2_inner_block_size, step2);
           step += block) {
        block = min (min (steps - log2_inner_block_size - step, step2 - step),
                     log2_outer_block_size);
        todo  = steps - 1 - step;
        blocks= todo - block + 1;
        c0= c;
        for (nat j= 0; j < ((nat) 1<<step); j++) {
          for (k_beg= 0; k_beg < ((nat) 1<<blocks); k_beg= k_end) {
            k_end= min (((nat) 1<<blocks),
                        k_beg + ((nat) 1<< log2_block_number));
            delocate (temp, (nat) 1<<block,
                      c0 + k_beg, stride<<blocks, k_end - k_beg);
            for (nat k= 0; k < k_end - k_beg; k++)
              ffter->dfft (temp+(k<<block), 1,
                           ((shift>>todo) + j)<<(block-1), block);
            relocate (temp, (nat) 1<<block,
                      c0+k_beg, stride<<blocks, k_end - k_beg);
          }
          c0 += (stride<<(todo+1));
        }
      }
      mmx_delete<C> (temp, temp_len);
      if (step < step2) {
        c0= c;
        todo= steps - 1 - step;
        for (nat j= 0; j < ((nat) 1<<step); j++) {
          ffter->dfft (c0, stride, ((shift>>todo) + j) << todo,
                       todo+1, 0, step2-step);
          c0 += (stride<<(todo+1));
        }
      }
    }
  }

  inline void
  ifft (C* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    if (steps <= max (threshold, log2_inner_block_size))
      ffter->ifft (c, stride, shift, steps, step1, step2);
    else {
      nat temp_len= default_aligned_size<C>
        ((nat) 1<<(log2_outer_block_size+log2_block_number));
      C* c0, * temp= mmx_new<C> (temp_len);
      nat step, block, blocks, todo, k_beg, k_end;

      if (log2_inner_block_size > steps - step2) {
        step= steps - log2_inner_block_size;
        c0= c;
        todo= steps - 1 - step;
        for (nat j= 0; j < ((nat) 1<<step); j++) {
          ffter->ifft (c0, stride, ((shift>>todo) + j) <<todo,
                       todo+1, max(0, (int) step1- (int) step), step2-step);
          c0 += (stride<<(todo+1));
        }
      }
      else
        step= step2;

      while (step > step1) {
        block = min (step - step1, log2_outer_block_size);
        step -= block;
        todo  = steps - 1 - step;
        blocks= todo - block + 1;
        c0= c;
        for (nat j= 0; j < ((nat) 1<<step); j++) {
          for (k_beg= 0; k_beg < ((nat) 1<<blocks); k_beg= k_end) {
            k_end= min (((nat) 1<<blocks),
                        k_beg + ((nat) 1<< log2_block_number));
            delocate (temp, (nat) 1<<block,
                      c0 + k_beg, stride<<blocks, k_end - k_beg);
            for (nat k= 0; k < k_end - k_beg; k++)
              ffter->ifft (temp+(k<<block), 1,
                           ((shift>>todo) + j)<<(block-1), block);
            relocate (temp, (nat) 1<<block,
                      c0+k_beg, stride<<blocks, k_end - k_beg);
          }
          c0 += (stride<<(todo+1));
        }
      }
      mmx_delete<C> (temp, temp_len);
    }
  }

  inline void
  dfft (C* c, nat stride, nat shift, nat steps) {
    dfft (c, stride, shift, steps, 0, steps); }

  inline void
  ifft (C* c, nat stride, nat shift, nat steps) {
    ifft (c, stride, shift, steps, 0, steps); }

  inline void
  direct_transform (C* c) {
    dfft (c, 1, 0, depth); }

  inline void
  inverse_transform (C* c, bool divide=true) {
    ifft (c, 1, 0, depth);
    if (divide) {
      S x= binpow (S (2), depth);
      x= invert (x);
      NVec::template vec_unary_scalar<typename R::fft_mul_sc_op> (c, x, len);
    }
  }
};

} // namespace mmx
#endif //__MMX__FFT_BLOCKS__HPP