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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_simd.hpp
* DESCRIPTION: FFT using SIMD SSE operations
* COPYRIGHT  : (C) 2007  Joris van der Hoeven
*******************************************************************************
* 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_SIMD__HPP
#define __MMX__FFT_SIMD__HPP
#include <numerix/simd.hpp>
#include <numerix/sse.hpp>
#include <algebramix/fft_naive.hpp>

namespace mmx {

// Fallback for when no simd type is available
template<typename C,
         typename FFTER= fft_naive_transformer<C>,
         typename FFTER_SIMD= fft_naive_transformer<typename Simd_type(C)>,
         nat thr= 2>
class fft_simd_transformer {
public:
  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_simd_transformer (nat n):
    ffter (new FFTER (n)),
    depth (ffter->depth), len (ffter->len), roots (ffter->roots)  {}

  inline ~fft_simd_transformer () { delete ffter; }

  template<typename CC> inline void
  dfft (CC* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    ffter->dfft (c, stride, shift, steps, step1, step2); }

  template<typename CC> inline void
  ifft (CC* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    ffter->ifft (c, stride, shift, steps, step1, step2); }

  template<typename CC> inline void
  dfft (CC* c, nat stride, nat shift, nat steps) {
    dfft (c, stride, shift, steps, 0, steps); }

  template<typename CC> inline void
  ifft (CC* 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) {
    typedef implementation<vector_linear,vector_naive> NVec;
    ifft (c, 1, 0, depth);
    if (divide) NVec::mul (c, invert (binpow (S (2), depth)), len); }
};
  
#ifdef NUMERIX_ENABLE_SIMD

/******************************************************************************
* Support for SSE instructions
******************************************************************************/

#ifdef __SSE2__

inline void
simd_encode (complex<double>* c, nat n) {
  double  temp;
  double* r= (double*) ((void*) c);
  for (; n!=0; r+=8, n-=4) {
    temp= r[1];
    r[1]= r[2];
    r[2]= temp;
    temp= r[5];
    r[5]= r[6];
    r[6]= temp;
  }
}

inline void
simd_decode (complex<double>* c, nat n) {
  double  temp;
  double* r= (double*) ((void*) c);
  for (; n!=0; r+=8, n-=4) {
    temp= r[1];
    r[1]= r[2];
    r[2]= temp;
    temp= r[5];
    r[5]= r[6];
    r[6]= temp;
  }
}

STMPL
struct roots_helper<complex<sse_double> >:
  public roots_helper<complex<double> > {

  typedef complex<sse_double> C;
  typedef complex<double>     U;
  typedef double              S;

  static inline void
  fft_cross (C* c1, C* c2) {
    double* z1= (double*) ((void*) c1);
    double* z2= (double*) ((void*) c2);
    sse_double re_z1= simd_load_aligned (z1);
    sse_double im_z1= simd_load_aligned (z1+2);
    sse_double re_z2= simd_load_aligned (z2);
    sse_double im_z2= simd_load_aligned (z2+2);
    simd_save_aligned (z2  , re_z1 - re_z2);
    simd_save_aligned (z2+2, im_z1 - im_z2);
    simd_save_aligned (z1  , re_z1 + re_z2);
    simd_save_aligned (z1+2, im_z1 + im_z2);
  }

  static inline void
  dfft_cross (C* c1, C* c2, const U* u) {
    double* z1= (double*) ((void*) c1);
    double* z2= (double*) ((void*) c2);
    double* u1= (double*) ((void*) u );
    sse_double re_z1= simd_load_aligned (z1);
    sse_double im_z1= simd_load_aligned (z1+2);
    sse_double re_z2= simd_load_aligned (z2);
    sse_double im_z2= simd_load_aligned (z2+2);
    sse_double re_u1= simd_load_duplicate (u1);
    sse_double im_u1= simd_load_duplicate (u1+1);
    sse_double re_u2= re_u1 * re_z2 - im_u1 * im_z2;
    sse_double im_u2= re_u1 * im_z2 + im_u1 * re_z2;
    simd_save_aligned (z2  , re_z1 - re_u2);
    simd_save_aligned (z2+2, im_z1 - im_u2);
    simd_save_aligned (z1  , re_z1 + re_u2);
    simd_save_aligned (z1+2, im_z1 + im_u2);
  }

  static inline void
  ifft_cross (C* c1, C* c2, const U* u) {
    double* z1= (double*) ((void*) c1);
    double* z2= (double*) ((void*) c2);
    double* u1= (double*) ((void*) u );
    sse_double re_z1= simd_load_aligned (z1);
    sse_double im_z1= simd_load_aligned (z1+2);
    sse_double re_z2= simd_load_aligned (z2);
    sse_double im_z2= simd_load_aligned (z2+2);
    sse_double re_u1= simd_load_duplicate (u1);
    sse_double im_u1= simd_load_duplicate (u1+1);
    sse_double re_u2= re_z1 - re_z2;
    sse_double im_u2= im_z1 - im_z2;
    simd_save_aligned (z2  , re_u1 * re_u2 - im_u1 * im_u2);
    simd_save_aligned (z2+2, re_u1 * im_u2 + im_u1 * re_u2);
    simd_save_aligned (z1  , re_z1 + re_z2);
    simd_save_aligned (z1+2, im_z1 + im_z2);
  }
};

STMPL
struct std_roots<complex<sse_double> > {
  typedef complex<sse_double> C;
  typedef cached_roots_helper<roots_helper<C> > roots_type;
};

/******************************************************************************
* The FFT transformer class
******************************************************************************/

template<typename FFTER, typename FFTER_SIMD, nat thr>
class fft_simd_transformer<complex<double>, FFTER, FFTER_SIMD, thr> {
public:
  typedef complex<double> C;
  typedef typename FFTER::R R;
  typedef typename R::U U;
  typedef typename R::S S;
  
  FFTER* ffter;
  nat    depth;
  nat    len;
  U*     roots;

  typedef complex<sse_double> C_SIMD;
  FFTER_SIMD* ffter_simd;

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

  inline ~fft_simd_transformer () { delete ffter; delete ffter_simd; }

  inline void
  dfft (C* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    if (steps <= thr || stride != 1)
      ffter->dfft (c, stride, shift, steps, step1, step2);
    else {
      if (steps == step2) {
        simd_encode (c, 1 << steps);
        ffter_simd->dfft ((C_SIMD*) ((void*) c), 1, shift>>1,
                          steps-1, step1, steps-1);
        simd_decode (c, 1 << steps);
        ffter->dfft (c, 1, shift, steps, steps-1, steps);
      }
      else {
        simd_encode (c, 1 << steps);
        ffter_simd->dfft ((C_SIMD*) ((void*) c), 1,
                          shift>>1, steps-1, step1, step2);
        simd_decode (c, 1 << steps);
      }
    }
  }

  inline void
  ifft (C* c, nat stride, nat shift, nat steps, nat step1, nat step2) {
    if (steps <= thr || stride != 1)
      ffter->ifft (c, stride, shift, steps, step1, step2);
    else {
      if (steps == step2) {
        ffter->ifft (c, 1, shift, steps, steps-1, steps);
        simd_encode (c, 1 << steps);
        ffter_simd->ifft ((C_SIMD*) ((void*) c), 1, shift>>1,
                          steps-1, step1, steps-1);
        simd_decode (c, 1 << steps);
      }
      else {
        simd_encode (c, 1 << steps);
        ffter_simd->ifft ((C_SIMD*) ((void*) c), 1,
                          shift>>1, steps-1, step1, step2);
        simd_decode (c, 1 << steps);
      }
    }
  }

  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, ffter->depth); }

  inline void
  inverse_transform (C* c, bool divide=true) {
    typedef implementation<vector_linear,vector_naive> NVec;
    ifft (c, 1, 0, depth);
    if (divide) NVec::mul (c, invert (binpow (S (2), depth)), len); }
};

#endif //__SSE2__
#endif // NUMERIX_ENABLE_SIMD
} // namespace mmx
#endif //__MMX__FFT_SIMD__HPP