//// spectral toolkit // copyright (c) 2005 university corporation for atmospheric research// licensed under the gnu general public license//#include "cfft2.h"#include "transpose.h"#include "decompose.h"#include "alloc.h"using namespace spectral;/// Multithreaded constructor./// \param nx first dimension of transform/// \param ny second dimension of transform/// \param nt number of threads used in transform/// \throws bad_parameter if nx<2, ny<2 or nt<1 cfft2::cfft2(int nx,int ny,int nt){    if((nx<2)||(ny<2)||(nt<1))        throw bad_parameter();      threads=nt;    t=new cfft2::thread*[threads];    g=new group(threads);    for(int i=0;i<threads;i++)        t[i]=new thread(nx,ny,i,g);}/// Single thread constructor./// \param nx first dimension of transform/// \param ny second dimension of thransform/// \throws bad_parameter if nx<2 or ny<2cfft2::cfft2(int nx,int ny){    if((nx<2)||(ny<2))        throw bad_parameter();    threads=1;    g=0;    t=new cfft2::thread*[1];    t[0]=new thread(nx,ny,0,g);}/// Destructor for 2-d complex FFT.cfft2::~cfft2(){  if(g)    delete g;  for(int i=0;i<threads;i++)    delete t[i];  delete [] t;}/// Forward transform.  Computes the 2-d FFT of a complex 2-d array/// \f[ \hat{a}_{k_2,k_1} = \frac{1}{N_1 N_2} \sum_{n_1=0}^{N_1-1} /// \sum_{n_2=0}^{N_2-1} a_{n_1,n_2} e^{-2 \pi i k_1 n_1/N_1} e^{-2 \pi i k_2 n_2/N_2} \f]/// where the input array is a[N1][N2] and the output array is normalized and transposed /// as b[N2][N1], for \f$ k_1 = 0, \ldots, N_1-1 \f$ and \f$ k_2 = 0, \ldots, N_2-1 \f$./// Note that the input array is used in the computation and is overwritten.  /// Also, the input and output arrays must not reference the same locations./// \param a 2-d input array in physical space/// \param b 2-d output array containing Fourier coefficients in transposed ordervoid cfft2::analysis(complex **a,complex **b){    analysis(&a,&b,1);}/// Backward transform.  Computes the 2-d inverse FFT of a complex 2-d array/// of spectral coefficients/// \f[ a_{n_1,n_2} = \sum_{k_1=0}^{N_1-1} \sum_{k_2=0}^{N_2-1} \hat{a}_{k_2,k_1} /// e^{2 \pi i k_1 n_1/N_1} e^{2 \pi i k_2 n_2/N_2} \f]/// where the input array is b[N2][N1] and the output array is a[N1][N2]./// Note that the input array is used in the computation and is overwritten./// Also, the input and output arrays must not reference the same locations./// \param a 2-d input array containing Fourier coefficients in transposed order/// \param b 2-d output array in physical spacevoid cfft2::synthesis(complex **b,complex **a){    synthesis(&b,&a,1);}/// Multiple instance forward transform./// Note that the input array is used in the computation and is overwritten./// Also, the input and output arrays must not reference the same locations./// \param a 3-d input array of size [m][n1][n2]/// \param b 3-d output array of size [m][n2][n1]/// \param m number of sequences to transform/// \param p offset (from the default zero) of first sequence in array void cfft2::analysis(complex ***a,complex ***b,int m,int p){  if((m<1)||(p<0))    throw bad_parameter();  int i;  for(i=0;i<threads;i++)    t[i]->analysis(a,b,m,p);  for(i=0;i<threads;i++)    t[i]->wait();}/// Multiple instance backward transform./// Note that the input array is used in the computation and is overwritten./// Also, the input and output arrays must not reference the same locations./// \param b 3-d input array of size [m][n2][n1]/// \param a 3-d output array of size [m][n1][n2]/// \param m number of sequences to transform/// \param p offset (from the default zero) of first sequence in array void cfft2::synthesis(complex ***b,complex ***a,int m,int p){    if((m<1)||(p<0))        throw bad_parameter();    int i;    for(i=0;i<threads;i++)        t[i]->synthesis(b,a,m,p);    for(i=0;i<threads;i++)        t[i]->wait();}/// Internal thread object constructor./// \param nx first dimension of transform/// \param ny second dimension of transform/// \param index thread index which is from 0 to size(G)-1/// \param G group object for cfft2cfft2::thread::thread(int nx,int ny,int index,group *G){    n1=nx;    n2=ny;    g=G;    FFT1=new cfft(n1);    FFT2=new cfft(n2);    fac=1.0/(((real)n1)*((real)n2));    int nt=max(size(g),1);    decompose(n1,nt,index,m1,p1);    decompose(n2,nt,index,m2,p2);}/// Internal thread object destructor.cfft2::thread::~thread(){    delete FFT2;    delete FFT1;}/// Internal thread object join method.void cfft2::thread::wait(){    join();}/// Internal thread object forward transform driver.void cfft2::thread::analysis(complex ***a,complex ***b,int m,int p){    A=a;    B=b;    M=m;    P=p;    Method=ANALYSIS;    if(g)        spawn();    else        start();}/// Internal thread object backward transform driver.void cfft2::thread::synthesis(complex ***b,complex ***a,int m,int p){    A=a;    B=b;    M=m;    P=p;    Method=SYNTHESIS;    if(g)        spawn();    else        start();}/// Internal thread object start method that executes upon thread spawn.void cfft2::thread::start(){    switch(Method)    {        case ANALYSIS:            for(int i=P;i<P+M;i++)            {                FFT2->analysis(A[i],m1,p1);                transpose(A[i],B[i],m1,p1,n2,fac);                sync(g);                FFT1->analysis(B[i],m2,p2);            }            break;        case SYNTHESIS:            for(int i=P;i<P+M;i++)            {                FFT1->synthesis(B[i],m2,p2);                transpose(B[i],A[i],m2,p2,n1);                sync(g);                FFT2->synthesis(A[i],m1,p1);            }            break;    }}