//// spectral toolkit // copyright (c) 2005 university corporation for atmospheric research// licensed under the gnu general public license//#include "rfft3.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 rfft3::rfft3(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/2+1,ny);    g=new group(threads);    for(int i=0;i<threads;i++)        t[i]=new thread(nx,ny,nz,this,i);}/// 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 rfft3::rfft3(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/2+1,ny);    t=new rfft3::thread*[1];    t[0]=new rfft3::thread(nx,ny,nz,this,0);}/// Destructor for 3-d real FFT.rfft3::~rfft3(){    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/2+1][n2][n1]/// \param m number of sequences to transformvoid rfft3::analysis(real ****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/2+1][n2][n1]/// \param a 4-d output array of size [m][n1][n2][n3]/// \param m number of sequences to transformvoid rfft3::synthesis(complex ****b,real ****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 real symmetric FFT of a real 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/2+1][N2][N1].  The input and output arrays can reference the same memory (using /// the alias function template) or can be distinct.  In the case that the arrays are accessing/// the same memory, ensure that the length array is at least sizeof(real)*nx*ny*2(nz/2+1)./// 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 rfft3::analysis(real ***a,complex ***b){    analysis(&a,&b,1);}/// Backward transform.  Computes the 3-d inverse FFT of the real symmetric 3-d array/// of complex 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/2+1][N2][N1] and the output array is a[N1][N2][N3],/// using the symmetry relation \f$\hat{a}_{k_3,k_2,k_1} = \overline{\hat{a}}_{N_3-k_3,N_2-k_2,N_1-k_1} \f$./// 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 rfft3::synthesis(complex ***b,real ***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 parent pointer to rfft3 object/// \param index thread index which is from 0 to size(G)-1rfft3::thread::thread(int nx,int ny,int nz,rfft3 *parent,int index){    n1=nx;    n2=ny;    N3=nz;    n3=N3/2+1;    Parent=parent;    C=Parent->C;;    g=Parent->g;    FFT1=new cfft(n1);    FFT2=new rfft2(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.rfft3::thread::~thread(){    delete FFT2;    delete FFT1;}/// Internal thread join.void rfft3::thread::wait(){    join();}/// Internal thread analysis driver.void rfft3::thread::analysis(real ****a,complex ****b,int m){    A=a;    B=b;    M=m;    Method=ANALYSIS;    if(g)        spawn();    else        start();}/// Internal thread synthesis driver.void rfft3::thread::synthesis(complex ****b,real ****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 rfft3::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);}