rfft2.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 "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<2

rfft2::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 order

void 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 space

void 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 cfft2

rfft2::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;
}