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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : fft_triadic_threads.hpp
* DESCRIPTION: multi-threaded triadic FFT multiplication
* 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_TRIADIC_THREADS__HPP
#define __MMX__FFT_TRIADIC_THREADS__HPP
#include <algebramix/fft_triadic_naive.hpp>
#include <basix/threads.hpp>

namespace mmx {

/******************************************************************************
* Multi-threaded and triadic FFT transformation
******************************************************************************/

template<typename C, typename FFTER= fft_triadic_naive_transformer<C>,
         nat thr= 9>
class fft_triadic_threads_transformer {
public:
  typedef implementation<vector_linear,vector_naive> NVec;
  typedef typename FFTER::V V;
  typedef typename V::template helper<C>::roots_type R;
  typedef typename R::U U;
  static const nat min_reps= 27; // Must be a power of 3

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

public:
  inline fft_triadic_threads_transformer (nat n):
    ffter (new FFTER (n)),
    depth (ffter->depth), len (ffter->len),
    roots (ffter->roots), stoor (ffter->stoor) {}
  inline ~fft_triadic_threads_transformer () { delete ffter; }

  template<typename CC> static inline void
  copy (CC* d, nat drs, nat dcs, CC* s, nat srs, nat scs, nat r, nat c) {
    for (nat j=0; j<c; j++, d+=dcs, s+=scs) {
      CC* dd= d; CC* ss= s;
      for (nat i=0; i<r; i++, dd+=drs, ss+=srs)
        *dd= *ss;
    }
  }

  template<typename CC>
  struct outer_fft_triadic_task_rep: public task_rep {
    FFTER* ffter;
    bool direct;
    CC *buf;    
    nat thr_reps, tot_reps;
    CC *c;
    nat stride, shift, steps;
    nat l;
  public:
    inline outer_fft_triadic_task_rep
     (FFTER* ffter2, bool direct2,
      CC* buf2, nat thr_reps2, nat tot_reps2,
      CC* c2, nat stride2, nat shift2, nat steps2):
        ffter (ffter2), direct (direct2),
        buf (buf2), thr_reps (thr_reps2), tot_reps (tot_reps2),
        c (c2), stride (stride2), shift (shift2), steps (steps2) {
      l= binpow ((nat) 3, steps); }
    inline ~outer_fft_triadic_task_rep () {}
    void execute () {
      for (nat i=0; i*min_reps < thr_reps; i++, c += stride * min_reps) {
        CC* aux= buf;
        copy (buf, 1, min_reps, c, stride, stride * tot_reps,
              min_reps, l);
        if (direct)
          for (nat j=0; j<min_reps; j++, aux++)
            ffter->dfft_triadic (aux, min_reps, shift, steps);
        else
          for (nat j=0; j<min_reps; j++, aux++)
            ffter->ifft_triadic (aux, min_reps, shift, steps);
        copy (c, stride, stride * tot_reps, buf, 1, min_reps,
              min_reps, l);
      }
    }
  };

  template<typename CC>
  struct inner_fft_triadic_task_rep: public task_rep {
    FFTER* ffter;
    bool direct;
    nat start, inc, total;
    CC *c;
    nat stride, shift, steps;
    nat l;
  public:
    inline inner_fft_triadic_task_rep
     (FFTER* ffter2, bool direct2, nat start2, nat inc2, nat total2,
      CC* c2, nat stride2, nat shift2, nat steps2):
        ffter (ffter2), direct (direct2),
        start (start2), inc (inc2), total (total2),
        c (c2), stride (stride2), shift (shift2), steps (steps2) {
      l= binpow ((nat) 3, steps); }
    inline ~inner_fft_triadic_task_rep () {}
    void execute () {
      if (direct)
        for (nat i=start; i<total; i+=inc)
          ffter->dfft_triadic (c + i * l * stride, stride,
                               shift + i * l, steps);
      else
        for (nat i=start; i<total; i+=inc)
          ffter->ifft_triadic (c + i * l * stride, stride,
                               shift + i * l, steps);   
    }
  };

  template<typename CC> void
  fft_triadic (bool direct, CC* c, nat stride, nat shift, nat steps) {
    nat nt= threads_number;
    nat half1= steps >> 1;
    nat half2= steps - half1;
    nat len1= binpow ((nat) 3, half1);
    nat len2= binpow ((nat) 3, half2);
    if (steps <= thr || len2 <= min_reps || nt == 1) {
      if (direct) ffter->dfft_triadic (c, stride, shift, steps);
      else ffter->ifft_triadic (c, stride, shift, steps);
    }
    else {
      for (nat stage=0; stage<2; stage++) {
        if ((stage == 0) ^ (!direct)) {
          // Below is multi-threaded the following loop:
          // C*  cc1  = c;
          // for (nat i=0; i< len2; i++, cc1 += stride)
          //   if (direct)
          //     ffter->dfft_triadic (cc1, stride * len2, shift, half1);
          //   else
          //     ffter->ifft_triadic (cc1, stride * len2, shift, half1);
          nat bsz= min_reps * len1; // Recall min_reps < len2
          CC* buf= mmx_new<CC> (nt * bsz);
          task tasks[nt];
          for (nat i=0; i<nt; i++) {
            nat sta= min_reps * ((  i   * len2) / (min_reps * nt));
            nat end= min_reps * (((i+1) * len2) / (min_reps * nt));
            tasks[i]= new outer_fft_triadic_task_rep<CC>
              (ffter, direct, buf + i*bsz, end-sta, len2,
               c + sta*stride, stride, shift / len2, half1);
          }
          threads_execute (tasks, nt);
          mmx_delete<CC> (buf, nt * bsz);
        }
        else {
          // Below is multi-threaded the following loop:
          // CC* cc2= c;
          // for (nat i=0; i<len1; i++, cc2 += stride*len2)
          //   if (direct)
          //     ffter->dfft_triadic (cc2, stride, shift + i*len2, half2);
          //   else
          //     ffter->ifft_triadic (cc2, stride, shift + i*len2, half2);
          task tasks[nt];
          for (nat i=0; i<nt; i++)
            tasks[i]= new inner_fft_triadic_task_rep<CC>
              (ffter, direct, i, nt, len1, c, stride, shift, half2);
          threads_execute (tasks, nt);
        }
      }
    }
  }

  template<typename CC> inline void
  dfft_triadic (CC* c, nat stride, nat shift, nat steps) {
    if (steps <= thr) ffter->dfft_triadic (c, stride, shift, steps);
    else fft_triadic (true, c, stride, shift, steps);
  }

  template<typename CC> inline void
  ifft_triadic (CC* c, nat stride, nat shift, nat steps) {
    if (steps <= thr) ffter->ifft_triadic (c, stride, shift, steps);
    else fft_triadic (false, c, stride, shift, steps);
  }

  template<typename CC> inline void
  dfft_triadic (CC* c, nat stride, nat shift, nat steps,
                nat step1, nat step2) {
    if (step1 == 0 && step2 == steps && steps > thr)
      fft_triadic (true, c, stride, shift, steps);
    else ffter->dfft_triadic (c, stride, shift, steps, step1, step2);
  }

  template<typename CC> inline void
  ifft_triadic (CC* c, nat stride, nat shift, nat steps,
                nat step1, nat step2) {
    if (step1 == 0 && step2 == steps && steps > thr)
      fft_triadic (false, c, stride, shift, steps);
    else ffter->ifft_triadic (c, stride, shift, steps, step1, step2);
  }

  inline void
  direct_transform_triadic (C* c) {
    dfft_triadic (c, 1, 0, depth);
  }

  inline void
  inverse_transform_triadic (C* c, bool shift=true) {
    ifft_triadic (c, 1, 0, depth);
    if (shift) NVec::mul (c, invert (C (len)), len);
  }
};

} // namespace mmx
#endif //__MMX__FFT_TRIADIC_THREADS__HPP