src/FFT/FFT.h

Go to the documentation of this file.

// -*- C++ -*-
/***************************************************************************
 *
 * The IPPL Framework
 * 
 *
 * Visit http://people.web.psi.ch/adelmann/ for more details
 *
 ***************************************************************************/

//--------------------------------------------------------------------------
// Class FFT
//--------------------------------------------------------------------------

#ifndef IPPL_FFT_FFT_H
#define IPPL_FFT_FFT_H

// include files
#include "FFT/FFTBase.h"

// forward declarations
template <unsigned Dim> class FieldLayout;
template <class T, unsigned Dim> class BareField;
template <class T, unsigned Dim> class LField;


class CCTransform {};
class RCTransform {};
class SineTransform {};

template <class Transform, unsigned Dim, class T>
class FFT : public FFTBase<Dim,T> {};


template <unsigned Dim, class T>
class FFT<CCTransform,Dim,T> : public FFTBase<Dim,T> {

public: 
  
  // typedefs
  typedef FieldLayout<Dim> Layout_t;
#ifdef IPPL_HAS_TEMPLATED_COMPLEX
  typedef complex<T> Complex_t;
#else
  typedef complex Complex_t;
#endif
  typedef BareField<Complex_t,Dim> ComplexField_t;
  typedef LField<Complex_t,Dim> ComplexLField_t;
  typedef typename FFTBase<Dim,T>::Domain_t Domain_t;
  
  FFT(const Domain_t& cdomain, const bool transformTheseDims[Dim],
      const bool& compressTemps=false);
  
  FFT(const Domain_t& cdomain, const bool& compressTemps=false)
    : FFTBase<Dim,T>(FFT<CCTransform,Dim,T>::ccFFT, cdomain,compressTemps) { 
    
    // construct array of axis lengths
    int lengths[Dim];
    unsigned d;
    for (d=0; d<Dim; ++d)
      lengths[d] = cdomain[d].length();
        
    // construct array of transform types for FFT Engine, compute normalization
    int transformTypes[Dim];
    T& normFact = this->getNormFact();
    normFact = 1.0;
    for (d=0; d<Dim; ++d) {
      transformTypes[d] = FFTBase<Dim,T>::ccFFT;  // all transforms are complex-to-complex
      normFact /= lengths[d];
    }
        
    // set up FFT Engine
    this->getEngine().setup(Dim, transformTypes, lengths);
    // set up the temporary fields
    setup();
  }
  
  
  // Destructor
  ~FFT(void);

  void transform(int direction, ComplexField_t& f, ComplexField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, ComplexField_t& f,
                 ComplexField_t& g, const bool& constInput=false);

  void transform(int direction, ComplexField_t& f);
  void transform(const char* directionName, ComplexField_t& f) {
    // invoke in-place transform function using direction name string
    int direction = this->getDirection(directionName);
    
    // Check domain of incoming Field
    const Layout_t& in_layout = f.getLayout();
    const Domain_t& in_dom = in_layout.getDomain();
    PAssert(checkDomain(this->getDomain(),in_dom));

    // Common loop iterate and other vars:
    unsigned d;
    int idim;            // idim loops over the number of transform dims.
    int begdim, enddim;  // beginning and end of transform dim loop
    unsigned nTransformDims = this->numTransformDims();
    // Field* for temp Field management:
    ComplexField_t* temp = &f;
    // Local work array passed to FFT:
    Complex_t* localdata;
    
    // Loop over the dimensions be transformed:
    begdim = (direction == +1) ? 0 : (nTransformDims-1);
    enddim = (direction == +1) ? nTransformDims : -1;
    for (idim = begdim; idim != enddim; idim += direction) {
      TAU_PROFILE_START(timer_swap);
      // Now do the serial transforms along this dimension:
      
      bool skipTranspose = false;
      // if this is the first transform dimension, we might be able
      // to skip the transpose into the first temporary Field
      if (idim == begdim) {
        // get domain for comparison
        const Domain_t& first_dom = tempLayouts_m[idim]->getDomain();
        // check that zeroth axis is the same and is serial
        // and that there are no guard cells
        skipTranspose = ( (in_dom[0].sameBase(first_dom[0])) &&
                          (in_dom[0].length() == first_dom[0].length()) &&
                          (in_layout.getDistribution(0) == SERIAL) &&
                          (f.getGC() == FFT<CCTransform,Dim,T>::nullGC) );
      }
      
      // if this is the last transform dimension, we might be able
      // to skip the last temporary and transpose right into f
      if (idim == enddim-direction) {
        // get domain for comparison
        const Domain_t& last_dom = tempLayouts_m[idim]->getDomain();
        // check that zeroth axis is the same and is serial
        // and that there are no guard cells
        skipTranspose = ( (in_dom[0].sameBase(last_dom[0])) &&
                          (in_dom[0].length() == last_dom[0].length()) &&
                          (in_layout.getDistribution(0) == SERIAL) &&
                          (f.getGC() == FFT<CCTransform,Dim,T>::nullGC) );
      }
      
      if (!skipTranspose) {
        // transpose and permute to Field with transform dim first
        (*tempFields_m[idim])[tempLayouts_m[idim]->getDomain()] = 
          (*temp)[temp->getLayout().getDomain()];
        
        // Compress out previous iterate's storage:
        if (this->compressTemps() && temp != &f) *temp = 0;
        temp = tempFields_m[idim];  // Field* management aid
      }
      else if (idim == enddim-direction && temp != &f) {
        // last transform and we can skip the last temporary field
        // so do the transpose here using f instead
        
        // transpose and permute to Field with transform dim first
        f[in_dom] = (*temp)[temp->getLayout().getDomain()];
        
        // Compress out previous iterate's storage:
        if (this->compressTemps()) *temp = 0;
        temp = &f;  // Field* management aid
      }
      TAU_PROFILE_STOP(timer_swap);
      
      TAU_PROFILE_START(timer_fft);
      // Loop over all the Vnodes, working on the LField in each.
      typename ComplexField_t::const_iterator_if l_i, l_end = temp->end_if();
      for (l_i = temp->begin_if(); l_i != l_end; ++l_i) {

        // Get the LField
        ComplexLField_t* ldf = (*l_i).second.get();
        // make sure we are uncompressed
        ldf->Uncompress();
        // get the raw data pointer
        localdata = ldf->getP();

        // Do 1D complex-to-complex FFT's on all the strips in the LField:
        int nstrips = 1, length = ldf->size(0);
        for (d=1; d<Dim; ++d) nstrips *= ldf->size(d);
        for (int istrip=0; istrip<nstrips; ++istrip) {
          // Do the 1D FFT:
          this->getEngine().callFFT(idim, direction, localdata);
          // advance the data pointer
          localdata += length;
        } // loop over 1D strips
      } // loop over all the LFields
      TAU_PROFILE_STOP(timer_fft);
      
    } // loop over all transformed dimensions 

    // skip final assignment and compress if we used f as final temporary
    if (temp != &f) {
      TAU_PROFILE_START(timer_swap);
      // Now assign back into original Field, and compress last temp's storage:
      f[in_dom] = (*temp)[temp->getLayout().getDomain()];
      if (this->compressTemps()) *temp = 0;
      TAU_PROFILE_STOP(timer_swap);
    }

    // Normalize:
    if (direction == +1)
      f *= Complex_t(this->getNormFact(), 0.0);
    return;
  }
private: 
  
  void setup(void);

  Layout_t** tempLayouts_m;

  ComplexField_t** tempFields_m;
  
};


template <unsigned Dim, class T>
inline void
FFT<CCTransform,Dim,T>::transform(
                                  const char* directionName,
                                  typename FFT<CCTransform,Dim,T>::ComplexField_t& f,
                                  typename FFT<CCTransform,Dim,T>::ComplexField_t& g,
                                  const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}


template <class T>
class FFT<CCTransform,1U,T> : public FFTBase<1U,T> {
  
public: 
  
  // typedefs
  typedef FieldLayout<1U> Layout_t;
#ifdef IPPL_HAS_TEMPLATED_COMPLEX
  typedef complex<T> Complex_t;
#else
  typedef complex Complex_t;
#endif
  typedef BareField<Complex_t,1U> ComplexField_t;
  typedef LField<Complex_t,1U> ComplexLField_t;
  typedef typename FFTBase<1U,T>::Domain_t Domain_t;

  // Constructors:

  FFT(const Domain_t& cdomain, const bool transformTheseDims[1U],
      const bool& compressTemps=false);
  FFT(const Domain_t& cdomain, const bool& compressTemps=false);

  // Destructor
  ~FFT(void);
  
  void transform(int direction, ComplexField_t& f, ComplexField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, ComplexField_t& f,
                 ComplexField_t& g, const bool& constInput=false);
  
  void transform(int direction, ComplexField_t& f);
  void transform(const char* directionName, ComplexField_t& f);

private: 

  void setup(void);
  
  Layout_t* tempLayouts_m;

  ComplexField_t* tempFields_m;
  
};


// inline function definitions

template <class T>
inline void
FFT<CCTransform,1U,T>::transform(
                                 const char* directionName,
                                 typename FFT<CCTransform,1U,T>::ComplexField_t& f,
                                 typename FFT<CCTransform,1U,T>::ComplexField_t& g,
                                 const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <class T>
inline void
FFT<CCTransform,1U,T>::transform(
                                 const char* directionName,
                                 typename FFT<CCTransform,1U,T>::ComplexField_t& f)
{
  int dir = this->getDirection(directionName);
  transform(dir, f);
  return;
}




template <unsigned Dim, class T>
class FFT<RCTransform,Dim,T> : public FFTBase<Dim,T> {
  
public: 
  
  // typedefs
  typedef FieldLayout<Dim> Layout_t;
  typedef BareField<T,Dim> RealField_t;
  typedef LField<T,Dim> RealLField_t;
#ifdef IPPL_HAS_TEMPLATED_COMPLEX
  typedef complex<T> Complex_t;
#else
  typedef complex Complex_t;
#endif
  typedef BareField<Complex_t,Dim> ComplexField_t;
  typedef LField<Complex_t,Dim> ComplexLField_t;
  typedef typename FFTBase<Dim, T>::Domain_t Domain_t;

  // Constructors:

  FFT(const Domain_t& rdomain, const Domain_t& cdomain,
      const bool transformTheseDims[Dim], const bool& compressTemps=false);

  FFT(const Domain_t& rdomain, const Domain_t& cdomain,
      const bool& compressTemps=false, int serialAxes = 1);

  // Destructor
  ~FFT(void);

  void transform(int direction, RealField_t& f, ComplexField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, RealField_t& f,
                 ComplexField_t& g, const bool& constInput=false);

  void transform(int direction, ComplexField_t& f, RealField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, ComplexField_t& f,
                 RealField_t& g, const bool& constInput=false);
  
private: 
  
  void setup(void);
  
  Layout_t** tempLayouts_m;
  
  Layout_t* tempRLayout_m;

  ComplexField_t** tempFields_m;

  RealField_t* tempRField_m;

  Domain_t complexDomain_m;

  int serialAxes_m;
};

   // Inline function definitions

template <unsigned Dim, class T>
inline void
FFT<RCTransform,Dim,T>::transform(
                                  const char* directionName,
                                  typename FFT<RCTransform,Dim,T>::RealField_t& f,
                                  typename FFT<RCTransform,Dim,T>::ComplexField_t& g,
                                  const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <unsigned Dim, class T>
inline void
FFT<RCTransform,Dim,T>::transform(
                                  const char* directionName,
                                  typename FFT<RCTransform,Dim,T>::ComplexField_t& f,
                                  typename FFT<RCTransform,Dim,T>::RealField_t& g,
                                  const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}


template <class T>
class FFT<RCTransform,1U,T> : public FFTBase<1U,T> {
  
public: 
  
  // typedefs
  typedef FieldLayout<1U> Layout_t;
  typedef BareField<T,1U> RealField_t;
  typedef LField<T,1U> RealLField_t;
#ifdef IPPL_HAS_TEMPLATED_COMPLEX
  typedef complex<T> Complex_t;
#else
  typedef complex Complex_t;
#endif
  typedef BareField<Complex_t,1U> ComplexField_t;
  typedef LField<Complex_t,1U> ComplexLField_t;
  typedef typename FFTBase<1U,T>::Domain_t Domain_t;
  
  // Constructors:
  
  FFT(const Domain_t& rdomain, const Domain_t& cdomain,
      const bool transformTheseDims[1U], const bool& compressTemps=false);
  FFT(const Domain_t& rdomain, const Domain_t& cdomain,
      const bool& compressTemps=false);

  ~FFT(void);
  
  void transform(int direction, RealField_t& f, ComplexField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, RealField_t& f,
                 ComplexField_t& g, const bool& constInput=false);

  void transform(int direction, ComplexField_t& f, RealField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, ComplexField_t& f,
                 RealField_t& g, const bool& constInput=false);
  
private: 
  
  void setup(void);
   
  Layout_t* tempLayouts_m;

  ComplexField_t* tempFields_m;

  Layout_t* tempRLayout_m;

  //  const Domain_t& complexDomain_m;
  Domain_t complexDomain_m;
};

template <class T>
   inline void
FFT<RCTransform,1U,T>::transform(
                                 const char* directionName,
                                 typename FFT<RCTransform,1U,T>::RealField_t& f,
                                 typename FFT<RCTransform,1U,T>::ComplexField_t& g,
                                 const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <class T>
inline void
FFT<RCTransform,1U,T>::transform(
                                 const char* directionName,
                                 typename FFT<RCTransform,1U,T>::ComplexField_t& f,
                                 typename FFT<RCTransform,1U,T>::RealField_t& g,
                                 const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <unsigned Dim, class T>
class FFT<SineTransform,Dim,T> : public FFTBase<Dim,T> {
  
public: 
  
  // typedefs
  typedef FieldLayout<Dim> Layout_t;
  typedef BareField<T,Dim> RealField_t;
  typedef LField<T,Dim> RealLField_t;
#ifdef IPPL_HAS_TEMPLATED_COMPLEX
  typedef complex<T> Complex_t;
#else
  typedef complex Complex_t;
#endif
  typedef BareField<Complex_t,Dim> ComplexField_t;
  typedef LField<Complex_t,Dim> ComplexLField_t;
  typedef typename FFTBase<Dim,T>::Domain_t Domain_t;
  
  FFT(const Domain_t& rdomain, const Domain_t& cdomain,
      const bool transformTheseDims[Dim],
      const bool sineTransformDims[Dim], const bool& compressTemps=false);
  FFT(const Domain_t& rdomain, const Domain_t& cdomain,
      const bool sineTransformDims[Dim], const bool& compressTemps=false);
  FFT(const Domain_t& rdomain, const bool sineTransformDims[Dim],
      const bool& compressTemps=false);
  FFT(const Domain_t& rdomain, const bool& compressTemps=false);
  
  ~FFT(void);

  void transform(int direction, RealField_t& f, ComplexField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, RealField_t& f,
                 ComplexField_t& g, const bool& constInput=false);

  void transform(int direction, ComplexField_t& f, RealField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, ComplexField_t& f,
                 RealField_t& g, const bool& constInput=false);

  void transform(int direction, RealField_t& f, RealField_t& g,
                 const bool& constInput=false);
  void transform(const char* directionName, RealField_t& f,
                 RealField_t& g, const bool& constInput=false);
  
  void transform(int direction, RealField_t& f);
  void transform(const char* directionName, RealField_t& f);
  
private: 
  
  void setup(void);
  
  
  
  bool sineTransformDims_m[Dim]; 
  
  unsigned numSineTransforms_m;  
  
  Layout_t** tempLayouts_m;

  Layout_t** tempRLayouts_m;

  ComplexField_t** tempFields_m;
  
  RealField_t** tempRFields_m;

  const Domain_t* complexDomain_m;  
};

template <unsigned Dim, class T>
inline void
FFT<SineTransform,Dim,T>::transform(
                                    const char* directionName,
                                    typename FFT<SineTransform,Dim,T>::RealField_t& f,
                                    typename FFT<SineTransform,Dim,T>::ComplexField_t& g,
                                    const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <unsigned Dim, class T>
inline void
FFT<SineTransform,Dim,T>::transform(
                                    const char* directionName,
                                    typename FFT<SineTransform,Dim,T>::ComplexField_t& f,
                                    typename FFT<SineTransform,Dim,T>::RealField_t& g,
                                    const bool& constInput)
{
  int dir = this->getDirection(directionName);
   transform(dir, f, g, constInput);
   return;
}

template <unsigned Dim, class T>
inline void
FFT<SineTransform,Dim,T>::transform(
                                    const char* directionName,
                                    typename FFT<SineTransform,Dim,T>::RealField_t& f,
                                    typename FFT<SineTransform,Dim,T>::RealField_t& g,
                                    const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <unsigned Dim, class T>
inline void
FFT<SineTransform,Dim,T>::transform(
                                    const char* directionName,
                                    typename FFT<SineTransform,Dim,T>::RealField_t& f)
{
  int dir = this->getDirection(directionName);
  transform(dir, f);
  return;
}


 template <class T>
 class FFT<SineTransform,1U,T> : public FFTBase<1U,T> {
   
 public: 
   
   typedef FieldLayout<1U> Layout_t;
   typedef BareField<T,1U> RealField_t;
   typedef LField<T,1U> RealLField_t;
   typedef typename FFTBase<1U,T>::Domain_t Domain_t;
   
   FFT(const Domain_t& rdomain, const bool sineTransformDims[1U],
       const bool& compressTemps=false);
   FFT(const Domain_t& rdomain, const bool& compressTemps=false);

   
   ~FFT(void);

   void transform(int direction, RealField_t& f, RealField_t& g,
                  const bool& constInput=false);
   void transform(const char* directionName, RealField_t& f,
                  RealField_t& g, const bool& constInput=false);
   
   void transform(int direction, RealField_t& f);
   void transform(const char* directionName, RealField_t& f);
   
   private: 

   void setup(void);

   Layout_t* tempRLayouts_m;

   RealField_t* tempRFields_m;

 };

template <class T>
inline void
FFT<SineTransform,1U,T>::transform(
                                   const char* directionName,
                                   typename FFT<SineTransform,1U,T>::RealField_t& f,
                                   typename FFT<SineTransform,1U,T>::RealField_t& g,
                                   const bool& constInput)
{
  int dir = this->getDirection(directionName);
  transform(dir, f, g, constInput);
  return;
}

template <class T>
inline void
FFT<SineTransform,1U,T>::transform(
                                   const char* directionName,
                                   typename FFT<SineTransform,1U,T>::RealField_t& f)
{
  int dir = this->getDirection(directionName);
  transform(dir, f);
  return;
}
#include "FFT/FFT.cpp"
#endif // IPPL_FFT_FFT_H

/***************************************************************************
 * $RCSfile: FFT.h,v $   $Author: adelmann $
 * $Revision: 1.1.1.1 $   $Date: 2003/01/23 07:40:25 $
 * IPPL_VERSION_ID: $Id: FFT.h,v 1.1.1.1 2003/01/23 07:40:25 adelmann Exp $ 
 ***************************************************************************/

---