cfft3.cpp

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
//

#include "cfft3.h"
#include "transpose.h"
#include "alloc.h"
#include "alias.h"
#include "decompose.h"

using namespace spectral;

/// Multithreaded constructor.
/// \param nx first dimension of transform
/// \param ny second dimension of transform
/// \param nz third dimension of transform
/// \param nt number of threads used in transform
/// \throws bad_parameter if nx<2, ny<2, nz<2  or nt<1 

cfft3::cfft3(int nx,int ny,int nz,int nt)
{
    if((nx<2)||(ny<2)||(nz<2)||(nt<1))
        throw bad_parameter();
    threads=nt;
    t=new thread*[threads];
    C=alloc<complex>(nx,nz,ny);
    g=new group(threads);
    for(int i=0;i<threads;i++)
        t[i]=new thread(nx,ny,nz,i,g,C);
}

/// Single thread constructor.
/// \param nx first dimension of transform
/// \param ny second dimension of transform
/// \param nz third dimension of transform
/// \throws bad_parameter if nx<2, ny<2 or nz<1 

cfft3::cfft3(int nx,int ny,int nz)
{
    if((nx<2)||(ny<2)||(nz<2))
        throw bad_parameter();
    threads=1;
    g=0;
    C=alloc<complex>(nx,nz,ny);
    t=new cfft3::thread*[1];
    t[0]=new cfft3::thread(nx,ny,nz,0,g,C);
}

/// Destructor for 3-d complex FFT.

cfft3::~cfft3()
{
    if(g)
        delete g;
    dealloc<complex>(C);
    for(int i=0;i<threads;i++)
        delete t[i];
    delete [] t;
}

/// Multiple instance forward transform.
/// \param a 4-d input array of size [m][n1][n2][n3]
/// \param b 4-d output array of size [m][n3][n2][n1]
/// \param m number of sequences to transform

void cfft3::analysis(complex ****a,complex ****b,int m)
{
    if(m<1)
        throw bad_parameter();
    int i;
    for(i=0;i<threads;i++)
        t[i]->analysis(a,b,m);
    for(i=0;i<threads;i++)
        t[i]->wait();
}

/// Multiple instance backward transform.
/// \param b 4-d input array of size [m][n3][n2][n1]
/// \param a 4-d output array of size [m][n1][n2][n3]
/// \param m number of sequences to transform

void cfft3::synthesis(complex ****b,complex ****a,int m)
{
    if(m<1)
        throw bad_parameter();
    int i;
    for(i=0;i<threads;i++)
        t[i]->synthesis(b,a,m);
    for(i=0;i<threads;i++)
        t[i]->wait();
}

/// Forward transform.  Computes the 3-d FFT of a complex 3-d array
/// \f[ \hat{a}_{k_3,k_2,k_1} = \frac{1}{N_1 N_2 N_3} \sum_{n_1=0}^{N_1-1} 
/// \sum_{n_2=0}^{N_2-1} \sum_{n_3=0}^{N_3-1} a_{n_1,n_2,n_3} 
/// e^{-2 \pi i k_1 n_1/N_1} e^{-2 \pi i k_2 n_2/N_2} e^{-2 \pi i k_3 n_3/N_3} \f]
/// where the input array is a[N1][N2][N3] and the output array is normalized and transposed 
/// as b[N3][N2][N1].  The input and output arrays can reference the same memory (using 
/// the alias function template) or can be distinct.
/// Note that the input array is used in the computation and is overwritten.  
/// \param a 3-d input array in physical space
/// \param b 3-d output array containing Fourier coefficients in transposed order

void cfft3::analysis(complex ***a,complex ***b)
{
    analysis(&a,&b,1);
}

/// Backward transform.  Computes the 3-d inverse FFT of a complex 3-d array
/// of spectral coefficients
/// \f[ a_{n_1,n_2,n_3} = \sum_{k_1=0}^{N_1-1} \sum_{k_2=0}^{N_2-1} \sum_{k_3=0}^{N_3-1} 
/// \hat{a}_{k_3,k_2,k_1} e^{2 \pi i k_1 n_1/N_1} e^{2 \pi i k_2 n_2/N_2} e^{2 \pi i k_3 n_3/N_3} \f]
/// where the input array is b[N3][N2][N1] and the output array is a[N1][N2][N3].
/// The input and output arrays may refer to the same location.
/// Note that the input array is used in the computation and is overwritten.
/// \param a 3-d input array containing Fourier coefficients in transposed order
/// \param b 3-d output array in physical space

void cfft3::synthesis(complex ***b,complex ***a)
{
    synthesis(&b,&a,1);
}

/// Internal thread object constructor.
/// \param nx first dimension of transform
/// \param ny second dimension of transform
/// \param nz second dimension of transform
/// \param index thread index which is from 0 to size(G)-1
/// \param G group object
/// \param W internal work array C used for transpose

cfft3::thread::thread(int nx,int ny,int nz,int index,group *G,complex ***W)
{
    n1=nx;
    n2=ny;
    n3=nz;
    C=W;
    g=G;
    FFT1=new cfft(n1);
    FFT2=new cfft2(n2,n3);
    fac=1.0/(real)n1;
    int nt=max(size(g),1);
    decompose(n1,nt,index,m1,p1);
    decompose(n3*n2,nt,index,m32,p32);
}

/// Internal thread destructor.

cfft3::thread::~thread()
{
    delete FFT2;
    delete FFT1;
}

/// Internal thread join.

void cfft3::thread::wait()
{
    join();
}

/// Internal thread analysis driver.

void cfft3::thread::analysis(complex ****a,complex ****b,int m)
{
    A=a;
    B=b;
    M=m;
    Method=ANALYSIS;
    if(g)
        spawn();
    else
        start();
}

/// Internal thread synthesis driver.

void cfft3::thread::synthesis(complex ****b,complex ****a,int m)
{
    A=a;
    B=b;
    M=m;
    Method=SYNTHESIS;
    if(g)
        spawn();
    else
        start();
}

/// Internal thread object start method that executes upon thread spawn.

void cfft3::thread::start()
{
    complex ***b=alias<complex,complex>(M,n2*n3,n1,B);
    complex **c=alias<complex,complex>(n1,n2*n3,C);
    switch(Method)
    {
        case ANALYSIS:
            for(int i=0;i<M;i++)
            {
                FFT2->analysis(A[i],C,m1,p1);
                sync(g);
                transpose(c,b[i],m1,p1,n2*n3,fac);
                sync(g);
                FFT1->analysis(b[i],m32,p32);
            }
            break;
        case SYNTHESIS:
            for(int i=0;i<M;i++)
            {
                FFT1->synthesis(b[i],m32,p32);
                sync(g);
                transpose(b[i],c,m32,p32,n1);
                sync(g);
                FFT2->synthesis(C,A[i],m1,p1);
            }
            break;
    }
    dealias<complex>(b);
    dealias<complex>(c);
}