//// spectral toolkit // copyright (c) 2005 university corporation for atmospheric research// licensed under the gnu general public license//#include "rfft2.h"#include "decompose.h"#include "alloc.h"#include "alias.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 rfft2::rfft2(int nx,int ny,int nt){  if((nx<2)||(ny<2)||(nt<1))    throw bad_parameter();  n1=nx;  N2=ny;  n2=ny/2+1;  threads=nt;  t=new rfft2::thread*[threads];  g=new group(threads);  for(int i=0;i<threads;i++)    t[i]=new rfft2::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<2rfft2::rfft2(int nx,int ny){  if((nx<2)||(ny<2))    throw bad_parameter();  threads=1;  g=0;  t=new rfft2::thread*[1];  t[0]=new rfft2::thread(nx,ny,0,g);}/// Destructor for 2-d real FFT.rfft2::~rfft2(){  if(g)    delete g;  for(int i=0;i<threads;i++)    delete t[i];  delete [] t;}/// 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 real input array of size [m][n1][n2]/// \param b 3-d complex output array of size [m][n2/2+1][n1]/// \param m number of sequences to transform/// \param p offset (from the default zero) of first sequence in array void rfft2::analysis(real ***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 complex input array of size [m][n2/2+1][n1]/// \param a 3-d real 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 rfft2::synthesis(complex ***b,real ***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();}/// Forward transform.  Computes the 2-d FFT of a real 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 real input array is a[N1][N2] and the complex output array is normalized and transposed /// as b[N2/2+1][N1], for \f$ k_1 = 0, \ldots, N_1-1 \f$ and /// \f$ k_2 = 0, \ldots, \lfloor N_2/2 \rfloor\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 rfft2::analysis(real **a,complex **b){  analysis(&a,&b,1);}/// Backward transform.  Computes the 2-d inverse FFT for a complex array/// of spectral coefficients of a real symmetric 2-d array/// \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 n_2/N_2} \f]/// using the property \f[ \hat{a}_{k_2,k_1} = \overline{\hat{a}}_{N_2-k_2,N_1-k_1} \f]/// where the complex input array is b[N2/2+1][N1] and the real 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 rfft2::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 index thread index which is from 0 to size(G)-1/// \param G group object for cfft2rfft2::thread::thread(int nx,int ny,int index,group *G){  n1=nx;  N2=ny;  n2=ny/2+1;  g=G;  FFT1=new cfft(n1);  FFT2=new rfft(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.rfft2::thread::~thread(){  delete FFT2;  delete FFT1;}/// Internal thread object join method.void rfft2::thread::wait(){  join();}/// Internal thread object forward transform driver.void rfft2::thread::analysis(real ***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 rfft2::thread::synthesis(complex ***b,real ***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 rfft2::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]);	  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);	  sync(g);	  transpose(B[i],A[i]);	  FFT2->synthesis(A[i],m1,p1);	}      break;    }}void rfft2::thread::transpose(real **a,complex **b){  const int block=32;  const int n2p=(N2+1)/2;  int I1,I2;  int i1,i2;  for(I1=p1;I1<p1+m1;I1+=block)    for(I2=1;I2<n2p;I2+=block)      for(i1=I1;i1<min(I1+block,p1+m1);i1++)	for(i2=I2;i2<min(I2+block,n2p);i2++)	  {	    b[i2][i1].re=fac*a[i1][2*i2-1];	    b[i2][i1].im=fac*a[i1][2*i2  ];	  }                      for(i1=p1;i1<p1+m1;i1++)    {      b[0][i1].re=fac*a[i1][0];      b[0][i1].im=0.0;    }                      if(N2%2==0)    for(i1=p1;i1<p1+m1;i1++)      {	b[n2-1][i1].re=fac*a[i1][N2-1];	b[n2-1][i1].im=0.0;      }}void rfft2::thread::transpose(complex **b,real **a){  const int block=32;  const int n2p=(N2+1)/2;  int I1,I2;  int i1,i2;      for(I2=1;I2<n2p;I2+=block)    for(I1=p1;I1<p1+m1;I1+=block)      for(i2=I2;i2<min(I2+block,n2p);i2++)	for(i1=I1;i1<min(I1+block,p1+m1);i1++)	  {	    a[i1][2*i2-1]=b[i2][i1].re;	    a[i1][2*i2  ]=b[i2][i1].im;	  }                      for(i1=p1;i1<p1+m1;i1++)    a[i1][0]=b[0][i1].re;      if(N2%2==0)    for(i1=p1;i1<p1+m1;i1++)      a[i1][N2-1]=b[n2-1][i1].re;}