FFTBase.h 8.56 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
// -*- C++ -*-
/***************************************************************************
 *
 * The IPPL Framework
 * 
 *
 * Visit http://people.web.psi.ch/adelmann/ for more details
 *
 ***************************************************************************/

//--------------------------------------------------------------------------
// Class FFTBase
//--------------------------------------------------------------------------

#ifndef IPPL_FFT_FFTBASE_H
#define IPPL_FFT_FFTBASE_H

// include files
#include "Utility/PAssert.h"
#include "Index/NDIndex.h"
#include "Field/GuardCellSizes.h"

#if defined(IPPL_USE_SCSL_FFT)
#include "FFT/SCSL_FFT.h"
#else
#include "FFT/fftpack_FFT.h"
#endif

#include <map>
#include <iostream>

// forward declarations
template <unsigned Dim, class T> class FFTBase;
template <unsigned Dim, class T>
std::ostream& operator<<(std::ostream&, const FFTBase<Dim,T>&);

/// character strings for transform types
inline 
std::string getTransformType(unsigned int i) 
{
adelmann's avatar
adelmann committed
41
    static const char* transformTypeString_g[4] = { "complex-to-complex FFT",
gsell's avatar
gsell committed
42
                                                    "real-to-complex FFT",
adelmann's avatar
adelmann committed
43 44
                                                    "sine transform",
                                                    "cosine transform" };
gsell's avatar
gsell committed
45

adelmann's avatar
adelmann committed
46
    return std::string(transformTypeString_g[i % 4]);
gsell's avatar
gsell committed
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
}

/*!
  The FFTBase class handles duties for the FFT class that do not involve
  the type of transform to be done.  FFTBase is templated on dimensionality
  of the Field to transform and the floating-point precision type of the
  Field (float or double).

  FFT Base Class to do stuff that is independent of transform type
*/
template <unsigned Dim, class T>
class FFTBase {

public: 
    // Some externally visible typedefs and enums.
    enum { dimensions = Dim };               // dimension
    typedef T Precision_t;                   // precision
    typedef NDIndex<Dim> Domain_t;           // domain type

    // Enumeration of transform types, used by derived FFT classes
adelmann's avatar
adelmann committed
67
    enum FFT_e { ccFFT, rcFFT, sineFFT, cosineFFT };
gsell's avatar
gsell committed
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

    // Type used for performing 1D FFTs
#if defined(IPPL_USE_SCSL_FFT)
    typedef SCSL<T> InternalFFT_t;
#else
    typedef FFTPACK<T> InternalFFT_t;
#endif

    FFTBase() {}  
  
    /** 
     * inputs are enum of transform type, domain of input Field,
     * which dimensions to transform, and whether to compress
     * temporary Fields when not in use
     * 
     * @param transform 
     * @param domain 
     * @param transformTheseDims 
     * @param compressTemps 
     */
    
    FFTBase(FFT_e transform, const Domain_t& domain,
	    const bool transformTheseDims[Dim], bool compressTemps);
    
    /** 
     * 
     * 
     * @param transform 
     * @param domain 
     * @param compressTemps 
     */
    
    FFTBase(FFT_e transform, const Domain_t& domain, bool compressTemps);
    
    // destructor
    virtual ~FFTBase(void) { delete [] activeDims_m; }
  
    /** 
     * I/O for FFT object
     * 
     * @param out 
     */    
    void write(std::ostream& out) const;

    /** 
     * Allow the user to name the transform directions, for code clarity.
     * 
     * @param direction 
     * @param directionName 
     */
    void setDirectionName(int direction, const char* directionName);

    /** 
     * Set the FFT normalization factor (to something other than the default)
     * 
     * @param nf 
     */
    void setNormFact(Precision_t nf) { normFact_m = nf; }

    /** 
     * Utility to determine the number of vnodes to use in temporary transpose
     * fields; this is either -1, or a limited number set on the command line
     * 
     * @return 
     */
    int transVnodes() const {
	if (Ippl::maxFFTNodes() > 0 && Ippl::maxFFTNodes() <= Ippl::getNodes())
	    return Ippl::maxFFTNodes();
	else
	    return (-1);
    }

protected:

    /**! 
       These members are used by the derived FFT classes
    */

    /// null GuardCellSizes object for checking BareField arguments to transform
    static GuardCellSizes<Dim> nullGC;

    /// translate direction name string into dimension number
    int getDirection(const char* directionName) const;

    /// query whether this dimension is to be transformed
    bool transformDim(unsigned d) const;

    /// query number of transform dimensions
    unsigned numTransformDims(void) const { return nTransformDims_m; }

    /// get dimension number from list of transformed dimensions
    unsigned activeDimension(unsigned d) const;

    /// access the internal FFT Engine
    InternalFFT_t& getEngine(void) { return FFTEngine_m; }

    /// get the FFT normalization factor
    Precision_t& getNormFact(void) { return normFact_m; }

    /// get our domain
    const Domain_t& getDomain(void) const { return Domain_m; }

    /// compare indexes of two domains
    bool checkDomain(const Domain_t& dom1, const Domain_t& dom2) const;

    /// do we compress temps?
    bool compressTemps(void) const { return compressTempFields_m; }

private: 

    /// Stores user-defined names for FFT directions:
    std::map<const char*,int> directions_m;

    FFT_e transformType_m;     ///< Indicates which type of transform we do
    bool transformDims_m[Dim]; ///< Indicates which dimensions are transformed.
    unsigned nTransformDims_m; ///< Stores the number of dims to be transformed
    unsigned* activeDims_m;    ///< Stores the numbers of these dims (0,1,2).

    /// Internal FFT object for performing serial FFTs.
    InternalFFT_t FFTEngine_m;

    /// Normalization factor:
    Precision_t normFact_m;

    /// Domain of the input field, mainly used to check axis sizes and ordering, former const Domain_t& Domain_m;
    Domain_t Domain_m;

    /// Switch to turn on/off compression of intermediate Fields (tempFields) as algorithm is finished with them
    bool compressTempFields_m;
};


// Inline function definitions

/// Define operator<< to invoke write() member function:
template <unsigned Dim, class T>
inline std::ostream&
operator<<(std::ostream& out, const FFTBase<Dim,T>& fft)
{
    fft.write(out);
    return out;
}

/** 
    Allow the user to name the transform directions, for code clarity.
    Typical values might be "x_to_k", "k_to_x", "t_to_omega", "omega_to_t"
*/
template <unsigned Dim, class T>
inline void
FFTBase<Dim,T>::setDirectionName(int direction,
                                 const char* directionName) {
219
    PAssert_EQ(std::abs(direction), 1);
gsell's avatar
gsell committed
220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
    directions_m[directionName] = direction;
    return;
}

/** 
 * Translate direction name string into dimension number
 * 
 * @param directionName 
 * 
 * @return 
 */
template <unsigned Dim, class T>
inline int
FFTBase<Dim,T>::getDirection(const char* directionName) const {
    return (*(directions_m.find(directionName))).second;
}

/** 
 * query whether this dimension is to be transformed
 * 
 * @param d 
 * 
 * @return 
 */
template <unsigned Dim, class T>
inline bool
FFTBase<Dim,T>::transformDim(unsigned d) const {
247
    PAssert_LT(d, Dim);
gsell's avatar
gsell committed
248 249 250 251 252 253 254 255 256 257 258 259 260
    return transformDims_m[d];
}

/** 
 * get dimension number from list of transformed dimensions
 * 
 * @param d 
 * 
 * @return 
 */
template <unsigned Dim, class T>
inline unsigned
FFTBase<Dim,T>::activeDimension(unsigned d) const {
261
    PAssert_LT(d, nTransformDims_m);
gsell's avatar
gsell committed
262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
    return activeDims_m[d];
}

/** 
 * helper function for comparing domains
 * 
 * @param Dim 
 * @param dom1 
 * @param Dim 
 * @param dom2 
 * 
 * @return 
 */
template <unsigned Dim, class T>
inline bool
FFTBase<Dim,T>::checkDomain(const FFTBase<Dim,T>::Domain_t& dom1,
                            const FFTBase<Dim,T>::Domain_t& dom2) const {
    // check whether domains are equivalent
    // we require that some permutation of the axes gives a matching domain.
    static bool matched[Dim];
    bool found;
    unsigned d, d1;
    // initialize matched array to false
    for (d=0; d<Dim; ++d) matched[d] = false;
    d=0;
    while (d<Dim) {
	d1=0;
	found = false;
	while (!found && d1<Dim) {
	    // if we have not yet found a match for this dimension,
	    // compare length and base of Index objects
	    if (!matched[d1]) {
		found = ( dom1[d].length()==dom2[d1].length() &&
			  dom1[d].sameBase(dom2[d1]) );
		// if equivalent, mark this dimension as matched
		if (found) matched[d1] = true;
	    }
	    ++d1;
	}
	if (!found) return false;
	++d;
    }
    return true;
}

307
#include "FFT/FFTBase.hpp"
gsell's avatar
gsell committed
308 309 310 311 312 313 314 315 316

#endif // IPPL_FFT_FFTBASE_H

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