rfft3.h

The Spectral Toolkit is a C++ spectral transform library written by Rodney James and Chuck Panaccion

//
// spectral toolkit 
// copyright (c) 2005 university corporation for atmospheric research
// licensed under the gnu general public license
//

#ifndef __rfft3__
#define __rfft3__

#include "rfft2.h"
#include "group.h"
#include "fiber.h"

namespace spectral
{
  /// Real symmetric three-dimensional FFT object.  Computes the normalized real symmetric 3-d FFT
  /// and its inverse using multi-threaded transforms.
  ///
  /// The following example computes the forward and backward transforms on an array
  /// of real values and compares the result with the original values.  Note the use of
  /// the alias function template to map the input and output arrays to the same memory
  /// locations.  In particular, note that the overall sizes of the input and output
  /// arrays must be the same and the use of 2*(nz/2+1) as the leading dimension of
  /// the real 3-d array used to alias to the complex array storing the spectral
  /// coefficients.
  /// \code
  /// #include <iostream>
  /// #include <rfft3.h>
  /// #include <alloc.h>
  /// #include <alias.h>
  ///
  /// using namespace spectral;
  ///
  /// int main()
  /// {
  ///   int nx=50,ny=60,nz=40;
  ///   int threads=2;
  ///   real ***a=alloc<real>(nx,ny,2*(nz/2+1));
  ///   complex ***b=alias<complex,real>(nz/2+1,ny,nx,a);
  ///   real ***c=alloc<real>(nx,ny,nz);
  ///   rfft3 W(nx,ny,nz,threads);
  ///   real dx=2.0/(real)nx;
  ///   real dy=2.0/(real)ny;
  ///   real dz=2.0/(real)ny;
  ///   for(int i=0;i<nx;i++)
  ///     for(int j=0;j<ny;j++)
  ///       for(int k=0;k<nz;k++)
  ///         a[i][j][k]=c[i][j][k]=(-2.0+i*dx)*(-2.0+k*dz)*(-2.0+j*dy);
  ///   W.analysis(a,b);
  ///   W.synthesis(b,a);
  ///   real error=0.0;
  ///   for(int i=0;i<nx;i++)
  ///     for(int j=0;j<ny;j++)
  ///       for(int k=0;k<nz;k++)
  ///         error=max(error,fabs(a[i][j][k]-c[i][j][k]));
  ///   std::cout << "error=" << error << std:: endl;
  ///   dealloc<real>(a);
  ///   dealias<complex>(b);
  ///   dealloc<real>(c);
  ///   return(0);
  /// }
  /// \endcode
  ///
  /// \sa cfft3

  class rfft3
  {
  public:
    rfft3(int,int,int,int);
    rfft3(int,int,int);
    ~rfft3();          
    void analysis(real ***,complex ***);
    void synthesis(complex ***,real ***);
    void analysis(real ****,complex ****,int);
    void synthesis(complex ****,real ****,int);
  private:
    /// Internal thread object for rfft3
    class thread : public fiber
    {
    public:
      thread(int,int,int,rfft3*,int);
      ~thread();          
      void analysis(real ****,complex ****,int);
      void synthesis(complex ****,real ****,int);
      void wait();
    private:
      group *g;
      real ****A;
      complex ****B;
      complex ***C;
      int n1,n2,n3,N3;
      int m1,m32,p1,p32;
      int M,P;
      cfft *FFT1;
      rfft2 *FFT2;
      real fac;
      void start();
      enum method {ANALYSIS,SYNTHESIS} Method;
      rfft3 *Parent;
    };
    complex ***C;
    group *g;
    int threads;
    thread **t;
    friend class thread;
  };
}

#endif

// Local Variables:
// mode:C++
// End: