//// 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 transformvoid 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 transformvoid 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 ordervoid 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 spacevoid 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 transposecfft3::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);}