|
/* FILE: fft.cc -*-Mode: c++-*-
|
|
*
|
|
* C++ code to do 1 and 2 dimensional FFT's.
|
|
*
|
|
*/
|
|
|
|
#include <string.h>
|
|
|
|
#include "nb.h"
|
|
#include "fft.h"
|
|
|
|
#ifdef USE_MPI
|
|
#include "mmsmpi.h"
|
|
#endif /* USE_MPI */
|
|
|
|
/* End includes */
|
|
|
|
#undef USE_COMPLEX
|
|
|
|
#ifndef OLD_CODE
|
|
#define CMULT(xr,xi,yr,yi,zr,zi) (zr) = (xr)*(yr)-(xi)*(yi), \
|
|
(zi) = (xr)*(yi)+(xi)*(yr)
|
|
#else
|
|
extern inline void CMULT(const MY_COMPLEX_REAL_TYPE &xr,
|
|
const MY_COMPLEX_REAL_TYPE &xi,
|
|
const MY_COMPLEX_REAL_TYPE &yr,
|
|
const MY_COMPLEX_REAL_TYPE &yi,
|
|
MY_COMPLEX_REAL_TYPE &zr,
|
|
MY_COMPLEX_REAL_TYPE &zi)
|
|
{
|
|
zr = xr*yr-xi*yi;
|
|
zi = xr*yi+xi*yr;
|
|
}
|
|
#endif // OLD_CODE
|
|
|
|
void FFT::ReleaseMemory(void)
|
|
{
|
|
if(vecsize>0) {
|
|
if(Uforward!=(MyComplex *)NULL) delete[] Uforward;
|
|
if(Uinverse!=(MyComplex *)NULL) delete[] Uinverse;
|
|
if(permindex!=(int *)NULL) delete[] permindex;
|
|
}
|
|
vecsize=0;
|
|
Uforward=Uinverse=(MyComplex *)NULL;
|
|
permindex=(int *)NULL;
|
|
}
|
|
|
|
FFT::~FFT(void)
|
|
{
|
|
ReleaseMemory();
|
|
}
|
|
|
|
void FFT::Setup(int size)
|
|
{
|
|
if(size==vecsize) return;
|
|
if(size<1) PlainError(1,"Error in FFT::Setup(int): "
|
|
"Requested length (%d) must be >0",size);
|
|
int k;
|
|
// Check that size is power of 2
|
|
for(k=size;k>2;k/=2)
|
|
if(k%2!=0) PlainError(1,"Error in FFT::Setup(int): "
|
|
"Requested length (%d) is not a power of 2",size);
|
|
ReleaseMemory();
|
|
vecsize=size;
|
|
|
|
// Allocate and setup MyComplex arrays
|
|
if((Uforward=new MyComplex[size])==0)
|
|
PlainError(1,"Error in FFT::Setup %s",ErrNoMem);
|
|
if((Uinverse=new MyComplex[size])==0)
|
|
PlainError(1,"Error in FFT::Setup %s",ErrNoMem);
|
|
#ifdef ORIG_CODE
|
|
double baseang= -2*PI/double(size);
|
|
Uforward[0]=Uinverse[0]=MyComplex(1,0);
|
|
double x,y;
|
|
for(k=1;k<size/2;k++) {
|
|
x = cos(baseang*k);
|
|
y = -sqrt(1-x*x);
|
|
y += (1-(x*x+y*y))/(2*y); // Tiny error correction
|
|
x += (1-(x*x+y*y))/(2*x);
|
|
Uforward[k]=Uinverse[size-k]=MyComplex(x,y);
|
|
Uforward[size-k]=Uinverse[k]=MyComplex(x,-y);
|
|
}
|
|
if(size>1) {
|
|
Uforward[size/2]=Uinverse[size/2]=MyComplex(-1,0);
|
|
}
|
|
#else
|
|
double baseang = -2*PI/double(size);
|
|
for(k=1;k<size/8;k++) {
|
|
double angle=k*baseang;
|
|
double y=sin(angle);
|
|
double x=cos(angle);
|
|
|
|
|
|
Uforward[k]=Uinverse[size-k]=MyComplex(x,y);
|
|
Uforward[size/2-k]=Uinverse[size/2+k]=MyComplex(-x,y);
|
|
Uforward[size/2+k]=Uinverse[size/2-k]=MyComplex(-x,-y);
|
|
Uforward[size-k]=Uinverse[k]=MyComplex(x,-y);
|
|
|
|
Uforward[size/4-k]=Uinverse[3*size/4+k]=MyComplex(-y,-x);
|
|
Uforward[size/4+k]=Uinverse[3*size/4-k]=MyComplex(y,-x);
|
|
Uforward[3*size/4-k]=Uinverse[size/4+k]=MyComplex(y,x);
|
|
Uforward[3*size/4+k]=Uinverse[size/4-k]=MyComplex(-y,x);
|
|
}
|
|
Uforward[0]=Uinverse[0]=MyComplex(1,0);
|
|
if(size>1) {
|
|
Uforward[size/2]=Uinverse[size/2]=MyComplex(-1,0);
|
|
}
|
|
if(size>3) {
|
|
Uforward[size/4]=Uinverse[3*size/4]=MyComplex(0,-1);
|
|
Uforward[3*size/4]=Uinverse[size/4]=MyComplex(0,1);
|
|
}
|
|
if(size>7) {
|
|
double x=SQRT1_2; // 1/sqrt(2)
|
|
double y=-x;
|
|
Uforward[size/8]=Uinverse[7*size/8]=MyComplex(x,y);
|
|
Uforward[3*size/8]=Uinverse[5*size/8]=MyComplex(-x,y);
|
|
Uforward[5*size/8]=Uinverse[3*size/8]=MyComplex(-x,-y);
|
|
Uforward[7*size/8]=Uinverse[size/8]=MyComplex(x,-y);
|
|
}
|
|
#endif // ORIG_CODE
|
|
|
|
// Allocate and setup (bit-reversal) permutation index
|
|
if((permindex=new int[size])==0)
|
|
PlainError(1,"Error in FFT::Setup %s",ErrNoMem);
|
|
permindex[0]=0;
|
|
int m,n; // The following code relies heavily on size==2^log2vecsize
|
|
for(k=1,n=size>>1;k<size;k++) {
|
|
// At each step, n is bit-reversed pattern of k
|
|
if(n>k) permindex[k]=n; // Swap index
|
|
else permindex[k]=0; // Do nothing: Index already swapped or the same
|
|
// Calculate next n
|
|
m=size>>1;
|
|
while(m>0 && n&m) { n-=m; m>>=1; }
|
|
n+=m;
|
|
}
|
|
}
|
|
|
|
void FFT::ForwardDecFreq(int size,MyComplex *vec,FFT_REAL_TYPE divisor)
|
|
{
|
|
if(divisor==0) divisor=1.0; // Default is no normalization on forward FFT
|
|
Setup(size);
|
|
BaseDecFreqForward(vec);
|
|
Permute(vec);
|
|
if(divisor!=0. && divisor!=1.) {
|
|
MY_COMPLEX_REAL_TYPE mult=1./divisor;
|
|
for(int k=0;k<size;k++) vec[k]*=mult;
|
|
}
|
|
}
|
|
|
|
void FFT::InverseDecTime(int size,MyComplex *vec,FFT_REAL_TYPE divisor)
|
|
{
|
|
if(divisor==0) divisor=(FFT_REAL_TYPE)size;
|
|
/// Default divisor on iFFT is 'size'
|
|
Setup(size);
|
|
Permute(vec);
|
|
BaseDecTimeInverse(vec);
|
|
if(divisor!=0. && divisor!=1.) {
|
|
MY_COMPLEX_REAL_TYPE mult=1./divisor;
|
|
for(int k=0;k<size;k++) vec[k]*=mult;
|
|
}
|
|
|
|
}
|
|
|
|
inline void Swap(MyComplex &a,MyComplex &b)
|
|
{ MyComplex c(a); a=b; b=c; }
|
|
|
|
void FFT::Permute(MyComplex *vec)
|
|
{ /* Bit reversal permutation */
|
|
int i,j;
|
|
for(i=0;i<vecsize;i++) {
|
|
if((j=permindex[i])!=0) Swap(vec[i],vec[j]);
|
|
}
|
|
}
|
|
|
|
void FFT::BaseDecFreqForward(MyComplex *vec)
|
|
{ // In-place forward FFT using Decimation in Frequency technique,
|
|
// *WITHOUT* resuffling of indices.
|
|
// NOTE 1: This code does not use the Complex class, because
|
|
// some compilers do not effectively optimize around
|
|
// Complex operations such as multiplication. So this
|
|
// routine just makes use of ordinary type "FFT_REAL_TYPE"
|
|
// variables, and assumes each Complex variable is
|
|
// actually two consecutive "MY_COMPLEX_REAL_TYPE" variables.
|
|
// NOTE 2: See notes in MJD's micromagnetics notebook, 11-Sep-96
|
|
// (p. 62) and 29-Sep-96 (p. 69).
|
|
// NOTE 3: This code has been optimized for performance on cascade.cam,
|
|
// a PentiumPro 200 MHz machine using a stock gcc 2.7.2
|
|
// compiler. In particular, the x86 chips suffer from a
|
|
// shortage of registers.
|
|
// NOTE 4: Some compromise made to RISC architectures on 27-May-1997,
|
|
// by moving all loads before any stores in main loop. As
|
|
// done, it hurts performance on PentiumPro-200 only a couple
|
|
// of percent. (mjd)
|
|
|
|
#define datat float /*Added by Guru on 11/10/2010: created data type datat; specified datat as required*/
|
|
typedef datat FFT_REAL_TYPE; /*Added by Guru on 11/10/2010: converted FFT_REAL_TYPE to datat from double*/
|
|
|
|
|
|
if(vecsize==1) return; // Nothing to do
|
|
|
|
MY_COMPLEX_REAL_TYPE *v;
|
|
|
|
MY_COMPLEX_REAL_TYPE const *const U=(MY_COMPLEX_REAL_TYPE *)Uforward;
|
|
MY_COMPLEX_REAL_TYPE *const dvec=(MY_COMPLEX_REAL_TYPE *)vec;
|
|
|
|
int block,blocksize,blockcount,offset,uoff1;
|
|
int halfbs,threehalfbs; // Half blocksize,3/2 blocksize
|
|
FFT_REAL_TYPE m1x,m1y,m2x,m2y,m3x,m3y;
|
|
FFT_REAL_TYPE x0,y0,x1,y1,x2,y2,x3,y3;
|
|
FFT_REAL_TYPE xs02,ys02,xd02,yd02,xs13,ys13,xd13,yd13;
|
|
FFT_REAL_TYPE t1x,t1y,t2x,t2y,t3x,t3y;
|
|
|
|
|
|
// Blocksize>4
|
|
for(blocksize=vecsize,blockcount=1;blocksize>4;
|
|
blocksize/=4,blockcount*=4) {
|
|
// Loop through double-step matrix multiplications
|
|
halfbs=blocksize/2; threehalfbs=blocksize+halfbs;
|
|
for(block=0,v=dvec;block<blockcount;block++,v+=2*blocksize) {
|
|
for(offset=0;offset<halfbs;offset+=2) {
|
|
uoff1=offset*blockcount;
|
|
m1x=U[uoff1]; m1y=U[uoff1+1];
|
|
m2x=U[2*uoff1]; m2y=U[2*uoff1+1];
|
|
m3x=U[3*uoff1]; m3y=U[3*uoff1+1];
|
|
|
|
x0=v[offset];
|
|
y0=v[offset+1];
|
|
x1=v[halfbs+offset];
|
|
y1=v[halfbs+offset+1];
|
|
x2=v[blocksize+offset];
|
|
y2=v[blocksize+offset+1];
|
|
x3=v[threehalfbs+offset];
|
|
y3=v[threehalfbs+offset+1];
|
|
|
|
xs02=x0+x2; xs13=x1+x3;
|
|
v[offset]=xs02+xs13;
|
|
t1x=xs02-xs13;
|
|
|
|
ys02=y0+y2; ys13=y1+y3;
|
|
v[offset+1]=ys02+ys13;
|
|
t1y=ys02-ys13;
|
|
|
|
v[halfbs+offset] =t1x*m2x-t1y*m2y;
|
|
v[halfbs+offset+1] =t1y*m2x+t1x*m2y;
|
|
|
|
xd02=x0-x2; yd13=y1-y3;
|
|
t3x=xd02-yd13; t2x=xd02+yd13;
|
|
yd02=y0-y2; xd13=x1-x3;
|
|
t3y=yd02+xd13; t2y=yd02-xd13;
|
|
v[blocksize+offset] =t2x*m1x-t2y*m1y;
|
|
v[blocksize+offset+1]=t2y*m1x+t2x*m1y;
|
|
v[threehalfbs+offset] =t3x*m3x-t3y*m3y;
|
|
v[threehalfbs+offset+1]=t3y*m3x+t3x*m3y;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Do smallest blocks; size is either 4 or 2
|
|
if(blocksize==4) {
|
|
blockcount=vecsize/4;
|
|
for(block=0,v=dvec;block<blockcount;block++,v+=8) {
|
|
x0=v[0]; y0=v[1]; x1=v[2]; y1=v[3];
|
|
x2=v[4]; y2=v[5]; x3=v[6]; y3=v[7];
|
|
|
|
xs02=x0+x2;
|
|
xs13=x1+x3;
|
|
v[0]=xs02+xs13;
|
|
v[2]=xs02-xs13;
|
|
|
|
ys02=y0+y2;
|
|
ys13=y1+y3;
|
|
v[1]=ys02+ys13;
|
|
v[3]=ys02-ys13;
|
|
|
|
xd02=x0-x2;
|
|
yd13=y1-y3;
|
|
v[4]=xd02+yd13;
|
|
v[6]=xd02-yd13;
|
|
|
|
yd02=y0-y2;
|
|
xd13=x1-x3;
|
|
v[5]=yd02-xd13;
|
|
v[7]=yd02+xd13;
|
|
}
|
|
}
|
|
else { // blocksize==2
|
|
blockcount=vecsize/2;
|
|
for(block=0,v=dvec;block<blockcount;block++,v+=4) {
|
|
x0=v[0]; y0=v[1];
|
|
x1=v[2]; y1=v[3];
|
|
v[0]=x0+x1; v[2]=x0-x1;
|
|
v[1]=y0+y1; v[3]=y0-y1;
|
|
}
|
|
}
|
|
}
|
|
|
|
void FFT::BaseDecTimeInverse(MyComplex *vec)
|
|
{ // In-place inverse FFT using Decimation in Time technique,
|
|
// *WITHOUT* resuffling of indices.
|
|
// NOTE 1: This code does not use the Complex class, because
|
|
// some compilers do not effectively optimize around
|
|
// Complex operations such as multiplication. So this
|
|
// routine just makes use of ordinary type "FFT_REAL_TYPE"
|
|
// variables, and assumes each Complex variable is
|
|
// actually two consecutive "MY_COMPLEX_REAL_TYPE" variables.
|
|
// NOTE 2: See notes in MJD's micromagnetics notebook, 11-Sep-96
|
|
// (p. 62) and 29-Sep-96 (p. 69).
|
|
// NOTE 3: This code has been optimized for performance on cascade.cam,
|
|
// a PentiumPro 200 MHz machine using a stock gcc 2.7.2
|
|
// compiler. In particular, the x86 chips suffer from a
|
|
// shortage of registers.
|
|
// NOTE 4: Some compromise made to RISC architectures on 27-May-1997,
|
|
// by moving all loads before any stores in main loop. As
|
|
// done, it hurts performance on PentiumPro-200 only a couple
|
|
// of percent. (mjd)
|
|
|
|
#define datat float /*Added by Guru on 11/10/2010: created data type datat; specified datat as required*/
|
|
typedef datat FFT_REAL_TYPE; /*Added by Guru on 11/10/2010: converted FFT_REAL_TYPE to datat from double*/
|
|
|
|
(vecsize==1) return; // Nothing to do
|
|
|
|
MY_COMPLEX_REAL_TYPE *v;
|
|
|
|
MY_COMPLEX_REAL_TYPE const *const U=(MY_COMPLEX_REAL_TYPE *)Uinverse;
|
|
MY_COMPLEX_REAL_TYPE *const dvec=(MY_COMPLEX_REAL_TYPE *)vec;
|
|
|
|
int block,blocksize,blockcount,offset,uoff1;
|
|
int halfbs,threehalfbs; // Half blocksize,3/2 blocksize
|
|
FFT_REAL_TYPE m1x,m1y,m2x,m2y,m3x,m3y;
|
|
FFT_REAL_TYPE x0,y0,x1,y1,x2,y2,x3,y3;
|
|
FFT_REAL_TYPE xs01,ys01,xd01,yd01,xs23,ys23,xd23,yd23;
|
|
FFT_REAL_TYPE t1x,t1y,t2x,t2y,t3x,t3y;
|
|
|
|
// Do smallest blocks; size is either 4 or 2
|
|
if(vecsize>2) {
|
|
blockcount=vecsize/4;
|
|
for(block=0,v=dvec;block<blockcount;block++,v+=8) {
|
|
x0=v[0]; y0=v[1];
|
|
x1=v[2]; y1=v[3];
|
|
x2=v[4]; y2=v[5];
|
|
x3=v[6]; y3=v[7];
|
|
xs01=x0+x1; xs23=x2+x3; // See NOTE 3 above
|
|
v[0]=xs01+xs23; v[4]=xs01-xs23;
|
|
ys01=y0+y1; ys23=y2+y3;
|
|
v[1]=ys01+ys23; v[5]=ys01-ys23;
|
|
xd01=x0-x1; yd23=y2-y3;
|
|
v[2]=xd01-yd23; v[6]=xd01+yd23;
|
|
yd01=y0-y1; xd23=x2-x3;
|
|
v[3]=yd01+xd23; v[7]=yd01-xd23;
|
|
}
|
|
}
|
|
else { // vecsize==2
|
|
x0=dvec[0]; y0=dvec[1];
|
|
x1=dvec[2]; y1=dvec[3];
|
|
dvec[0]=x0+x1;
|
|
dvec[2]=x0-x1;
|
|
dvec[1]=y0+y1;
|
|
dvec[3]=y0-y1;
|
|
return;
|
|
}
|
|
|
|
// Blocksize>4
|
|
for(blocksize=16,blockcount=vecsize/16;blocksize<=vecsize;
|
|
blocksize*=4,blockcount/=4) {
|
|
// Loop through double-step matric multiplications
|
|
halfbs=blocksize/2; threehalfbs=blocksize+halfbs;
|
|
for(block=0,v=dvec;block<blockcount;block++,v+=2*blocksize) {
|
|
for(offset=0;offset<blocksize/2;offset+=2) {
|
|
x0=v[offset];
|
|
y0=v[offset+1];
|
|
t2x=v[blocksize+offset];
|
|
t2y=v[blocksize+offset+1];
|
|
uoff1=offset*blockcount;
|
|
m1x=U[uoff1]; m1y=U[uoff1+1];
|
|
x2=t2x*m1x-t2y*m1y;
|
|
y2=t2y*m1x+t2x*m1y;
|
|
|
|
m2x=U[2*uoff1]; m2y=U[2*uoff1+1];
|
|
t1x=v[halfbs+offset];
|
|
t1y=v[halfbs+offset+1];
|
|
x1=t1x*m2x-t1y*m2y;
|
|
y1=t1y*m2x+t1x*m2y;
|
|
|
|
t3x=v[threehalfbs+offset];
|
|
t3y=v[threehalfbs+offset+1];
|
|
m3x=U[3*uoff1]; m3y=U[3*uoff1+1];
|
|
x3=t3x*m3x-t3y*m3y;
|
|
y3=t3y*m3x+t3x*m3y;
|
|
|
|
xs01=x0+x1; xs23=x2+x3;
|
|
v[ offset ] = xs01+xs23;
|
|
v[ blocksize+offset ] = xs01-xs23;
|
|
|
|
ys01=y0+y1; ys23=y2+y3;
|
|
v[ offset+1] = ys01+ys23;
|
|
v[ blocksize+offset+1] = ys01-ys23;
|
|
|
|
yd01=y0-y1; xd23=x2-x3;
|
|
v[ halfbs +offset+1] = yd01+xd23;
|
|
v[threehalfbs+offset+1] = yd01-xd23;
|
|
|
|
xd01=x0-x1; yd23=y2-y3;
|
|
v[ halfbs +offset ] = xd01-yd23;
|
|
v[threehalfbs+offset ] = xd01+yd23;
|
|
}
|
|
}
|
|
}
|
|
|
|
if(blocksize==2*vecsize) {
|
|
// We still have to do one single-step matrix multiplication
|
|
blocksize=vecsize; v=dvec;
|
|
for(offset=0;offset<blocksize;offset+=2,v+=2) {
|
|
x0=v[0]; y0=v[1];
|
|
t1x=v[vecsize]; t1y=v[vecsize+1];
|
|
m1x=U[offset]; m1y=U[offset+1];
|
|
x1=t1x*m1x-t1y*m1y; y1=t1y*m1x+t1x*m1y;
|
|
v[0] = x0+x1;
|
|
v[vecsize] = x0-x1;
|
|
v[1] = y0+y1;
|
|
v[vecsize+1] = y0-y1;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
const REALWIDE FFTReal2D::CRRCspeedratio=1.10;
|
|
/// Relative speed of ForwardCR as compared to ForwardRC.
|
|
/// If bigger than 1, then ForwardCR is faster. This will
|
|
/// be machine & compiler dependent...oh well.
|
|
|
|
void FFTReal2D::Setup(int size1,int size2)
|
|
{ // Note: This routine is also called by FFTReal2D::SetupInverse()
|
|
if(size1==vecsize1 && size2==vecsize2) return; // Nothing to do
|
|
|
|
// Check that sizes are powers of 2, and >= 1. Also extract
|
|
// base-2 log of sizes
|
|
if(size1<1) PlainError(1,"Error in FFTReal2D::Setup(int): "
|
|
"Requested size1 (%d) must be >=1",size1);
|
|
if(size2<1) PlainError(1,"Error in FFTReal2D::Setup(int): "
|
|
"Requested size2 (%d) must be >=1",size2);
|
|
|
|
int k;
|
|
if(size1==1) {
|
|
logsize1=0;
|
|
} else {
|
|
for(k=size1,logsize1=1;k>2;k/=2,logsize1++)
|
|
if(k%2!=0) PlainError(1,"Error in FFTReal2D::Setup(int): "
|
|
"Requested size1 (%d) is not a power of 2",size1);
|
|
}
|
|
if(size2==1) {
|
|
logsize2=0;
|
|
} else {
|
|
for(k=size2,logsize2=1;k>2;k/=2,logsize2++)
|
|
if(k%2!=0) PlainError(1,"Error in FFTReal2D::Setup(int): "
|
|
"Requested size2 (%d) is not a power of 2",size2);
|
|
}
|
|
|
|
// Allocate new space
|
|
ReleaseMemory();
|
|
scratch=new MyComplex[OC_MAX(size1,size2)];
|
|
scratchb=new MyComplex[OC_MAX(size1,size2)];
|
|
vecsize1=size1; vecsize2=size2;
|
|
}
|
|
|
|
void FFTReal2D::SetupInverse(int size1,int size2)
|
|
{
|
|
if(size1==vecsize1 && size2==vecsize2 && workarr!=NULL)
|
|
return; // Nothing to do
|
|
Setup(size1,size2);
|
|
if(workarr!=NULL) { delete[] workarr[0]; delete[] workarr; } // Safety
|
|
int rowcount=(vecsize1/2)+1;
|
|
workarr=new MyComplex*[rowcount];
|
|
workarr[0]=new MyComplex[rowcount*vecsize2];
|
|
for(int i=1;i<rowcount;i++) workarr[i]=workarr[i-1]+vecsize2;
|
|
}
|
|
|
|
void FFTReal2D::ReleaseMemory()
|
|
{
|
|
if(vecsize1==0 || vecsize2==0) return;
|
|
delete[] scratch; scratch=NULL;
|
|
delete[] scratchb; scratchb=NULL;
|
|
if(workarr!=NULL) {
|
|
delete[] workarr[0];
|
|
delete[] workarr;
|
|
workarr=NULL;
|
|
}
|
|
vecsize1=0; vecsize2=0;
|
|
fft1.ReleaseMemory(); fft2.ReleaseMemory();
|
|
}
|
|
|
|
void FFTReal2D::FillOut(int csize1,int csize2,MyComplex** carr)
|
|
{ // This routine assumes carr is a top-half-filled DFT of
|
|
// a real function, and fills in the bottom half using the
|
|
// relation
|
|
// carr[csize1-i][csize2-j]=conj(carr[i][j])
|
|
// for i>csize1/2, with the second indices interpreted 'mod csize2'.
|
|
int i,j;
|
|
for(i=1;i<csize1/2;i++) {
|
|
carr[csize1-i][0]=conj(carr[i][0]);
|
|
for(j=1;j<csize2;j++) {
|
|
carr[csize1-i][j]=conj(carr[i][csize2-j]);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
void FFTReal2D::ForwardCR(int rsize1,int rsize2,
|
|
const double* const* rarr,
|
|
int csize1,int csize2,MyComplex** carr)
|
|
{
|
|
Setup(OC_MAX(1,2*(csize1-1)),csize2); // Safety
|
|
if(vecsize1<2 || vecsize2<2)
|
|
PlainError(1,"Error in FFTReal2D::ForwardCR(...): "
|
|
"Full array dimensions (%dx%d) must be both >=2",
|
|
vecsize1,vecsize2);
|
|
|
|
int i,j;
|
|
FFT_REAL_TYPE x1,x2,y1,y2;
|
|
FFT_REAL_TYPE xb1,xb2,yb1,yb2;
|
|
|
|
// Do FFT on columns, 2 at a time
|
|
for(j=0;j+3<rsize2;j+=4) {
|
|
// Pack into MyComplex scratch array
|
|
for(i=0;i<rsize1;i++) {
|
|
scratch[i]=MyComplex(rarr[i][j],rarr[i][j+1]);
|
|
scratchb[i]=MyComplex(rarr[i][j+2],rarr[i][j+3]);
|
|
}
|
|
// Zero pad scratch space
|
|
for(i=rsize1;i<vecsize1;i++) scratch[i]=MyComplex(0.,0.);
|
|
for(i=rsize1;i<vecsize1;i++) scratchb[i]=MyComplex(0.,0.);
|
|
// Do complex FFT
|
|
fft1.ForwardDecFreq(vecsize1,scratchb);
|
|
fft1.ForwardDecFreq(vecsize1,scratch);
|
|
// Unpack into top half of 2D complex array
|
|
// Rows 0 & vecsize1/2 are real-valued, so pack them together
|
|
// into row 0 (row 0 as real part, row vecsize1/2 as imag. part).
|
|
carr[0][j] =MyComplex(scratch[0].real(),scratch[vecsize1/2].real());
|
|
carr[0][j+1]=MyComplex(scratch[0].imag(),scratch[vecsize1/2].imag());
|
|
carr[0][j+2]=MyComplex(scratchb[0].real(),scratchb[vecsize1/2].real());
|
|
carr[0][j+3]=MyComplex(scratchb[0].imag(),scratchb[vecsize1/2].imag());
|
|
for(i=1;i<vecsize1/2;i++) { // ASSUMES vecsize1 is even!
|
|
x1=scratch[i].real()/2; y1=scratch[i].imag()/2;
|
|
x2=scratch[vecsize1-i].real()/2; y2=scratch[vecsize1-i].imag()/2;
|
|
xb1=scratchb[i].real()/2; yb1=scratchb[i].imag()/2;
|
|
xb2=scratchb[vecsize1-i].real()/2; yb2=scratchb[vecsize1-i].imag()/2;
|
|
carr[i][j] =MyComplex(x1+x2,y1-y2);
|
|
carr[i][j+1] =MyComplex(y1+y2,x2-x1);
|
|
carr[i][j+2] =MyComplex(xb1+xb2,yb1-yb2);
|
|
carr[i][j+3] =MyComplex(yb1+yb2,xb2-xb1);
|
|
}
|
|
}
|
|
// Case rsize2 not divisible by 4
|
|
for(;j<rsize2;j+=2) {
|
|
// Pack into complex scratch array
|
|
if(j+1<rsize2) {
|
|
for(i=0;i<rsize1;i++) scratch[i]=MyComplex(rarr[i][j],rarr[i][j+1]);
|
|
}
|
|
else { // rsize2 == 1 mod 2.
|
|
for(i=0;i<rsize1;i++) scratch[i]=MyComplex(rarr[i][j],0.);
|
|
}
|
|
for(i=rsize1;i<vecsize1;i++) scratch[i]=MyComplex(0.,0.);
|
|
fft1.ForwardDecFreq(vecsize1,scratch);
|
|
carr[0][j] =MyComplex(scratch[0].real(),scratch[vecsize1/2].real());
|
|
carr[0][j+1]=MyComplex(scratch[0].imag(),scratch[vecsize1/2].imag());
|
|
for(i=1;i<vecsize1/2;i++) { // ASSUMES vecsize1 is even!
|
|
x1=scratch[i].real()/2; y1=scratch[i].imag()/2;
|
|
x2=scratch[vecsize1-i].real()/2; y2=scratch[vecsize1-i].imag()/2;
|
|
carr[i][j] =MyComplex(x1+x2,y1-y2);
|
|
carr[i][j+1] =MyComplex(y1+y2,x2-x1);
|
|
}
|
|
}
|
|
// Zero-pad remaining columns
|
|
if(rsize2<csize2) {
|
|
for(i=0;i<csize1;i++) for(j=rsize2;j<csize2;j++)
|
|
carr[i][j]=MyComplex(0.,0.);
|
|
// Note: One _may_ be able to gain a few percent speedup
|
|
// by using the 'memcpy' C-library routine.
|
|
}
|
|
|
|
// Do FFT on top half of rows
|
|
for(i=0;i<csize1-1;i++) fft2.ForwardDecFreq(csize2,carr[i]);
|
|
|
|
// Pull out row 0 & row csize1-1 from (packed) row 0
|
|
carr[csize1-1][0] = MyComplex(carr[0][0].imag(),0.);
|
|
carr[0][0] = MyComplex(carr[0][0].real(),0.);
|
|
for(j=1;j<csize2/2;j++) {
|
|
x1=carr[0][j].real()/2; y1=carr[0][j].imag()/2;
|
|
x2=carr[0][csize2-j].real()/2; y2=carr[0][csize2-j].imag()/2;
|
|
MyComplex temp1(x1+x2,y1-y2);
|
|
MyComplex temp2(y1+y2,x2-x1);
|
|
carr[0][j] = temp1;
|
|
carr[0][csize2-j] = conj(temp1);
|
|
carr[csize1-1][j] = temp2;
|
|
carr[csize1-1][csize2-j] = conj(temp2);
|
|
}
|
|
carr[csize1-1][csize2/2] = MyComplex(carr[0][csize2/2].imag(),0.);
|
|
carr[0][csize2/2] = MyComplex(carr[0][csize2/2].real(),0.);
|
|
}
|
|
|
|
void FFTReal2D::ForwardRC(int rsize1,int rsize2,
|
|
const double* const* rarr,
|
|
int csize1,int csize2,MyComplex** carr)
|
|
{
|
|
Setup(OC_MAX(1,2*(csize1-1)),csize2); // Safety
|
|
if(vecsize1<2 || vecsize2<2)
|
|
PlainError(1,"Error in FFTReal2D::ForwardRC(...): "
|
|
"Full array dimensions (%dx%d) must be both >=2",
|
|
vecsize1,vecsize2);
|
|
|
|
int i,j;
|
|
FFT_REAL_TYPE x1,x2,y1,y2;
|
|
FFT_REAL_TYPE xb1,xb2,yb1,yb2;
|
|
|
|
// Do row FFT's
|
|
for(i=0;i+1<rsize1;i+=2) {
|
|
// Pack 'MyComplex' row
|
|
for(j=0;j<rsize2;j++)
|
|
carr[i/2][j]=MyComplex(rarr[i][j],rarr[i+1][j]);
|
|
for(j=rsize2;j<csize2;j++) carr[i/2][j]=MyComplex(0.,0.); // Zero pad
|
|
// Do FFT
|
|
fft2.ForwardDecFreq(csize2,carr[i/2]);
|
|
}
|
|
for(;i<rsize1;i+=2) { // In case rsize1 == 1 mod 2
|
|
for(j=0;j<rsize2;j++)
|
|
carr[i/2][j]=MyComplex(rarr[i][j],0.);
|
|
for(j=rsize2;j<csize2;j++) carr[i/2][j]=MyComplex(0.,0.); // Zero pad
|
|
// Do FFT
|
|
fft2.ForwardDecFreq(csize2,carr[i/2]);
|
|
}
|
|
// Any remaining rows are zero padding on the fly during
|
|
// the column FFT's (see below).
|
|
|
|
// Do column FFT's
|
|
// Do column 0 and csize2/2, making use of the fact that
|
|
// these 2 columns are 'real'.
|
|
for(i=0;i<(rsize1+1)/2;i++) {
|
|
x1=carr[i][0].real(); x2=carr[i][0].imag();
|
|
y1=carr[i][csize2/2].real(); y2=carr[i][csize2/2].imag();
|
|
scratch[2*i] = MyComplex(x1,y1);
|
|
scratch[(2*i)+1] = MyComplex(x2,y2);
|
|
}
|
|
for(i*=2;i<vecsize1;i++) scratch[i]=MyComplex(0.,0.); // Zero pad
|
|
fft1.ForwardDecFreq(vecsize1,scratch);
|
|
carr[0][0] = MyComplex(scratch[0].real(),0.);
|
|
carr[0][csize2/2] = MyComplex(scratch[0].imag(),0.);
|
|
for(i=1;i<csize1-1;i++) {
|
|
x1=scratch[i].real()/2; y1=scratch[i].imag()/2;
|
|
x2=scratch[vecsize1-i].real()/2; y2=scratch[vecsize1-i].imag()/2;
|
|
carr[i][0] = MyComplex(x1+x2,y1-y2);
|
|
carr[i][csize2/2] = MyComplex(y1+y2,x2-x1);
|
|
}
|
|
carr[csize1-1][0] = MyComplex(scratch[csize1-1].real(),0.);
|
|
carr[csize1-1][csize2/2] = MyComplex(scratch[csize1-1].imag(),0.);
|
|
|
|
// Do remaining columns
|
|
for(j=1;j+1<csize2/2;j+=2) {
|
|
for(i=0;i<(rsize1+1)/2;i++) {
|
|
x1 =carr[i][j].real()/2; y1 =carr[i][j].imag()/2;
|
|
xb1=carr[i][j+1].real()/2; yb1=carr[i][j+1].imag()/2;
|
|
xb2=carr[i][csize2-1-j].real()/2; yb2=carr[i][csize2-1-j].imag()/2;
|
|
x2 =carr[i][csize2-j].real()/2; y2 =carr[i][csize2-j].imag()/2;
|
|
scratch[2*i] = MyComplex(x1+x2,y1-y2);
|
|
scratch[(2*i)+1] = MyComplex(y1+y2,x2-x1);
|
|
scratchb[2*i] = MyComplex(xb1+xb2,yb1-yb2);
|
|
scratchb[(2*i)+1] = MyComplex(yb1+yb2,xb2-xb1);
|
|
}
|
|
for(i*=2;i<vecsize1;i++)
|
|
scratch[i]= scratchb[i]=MyComplex(0.,0.); // Zero pad
|
|
fft1.ForwardDecFreq(vecsize1,scratchb);
|
|
fft1.ForwardDecFreq(vecsize1,scratch);
|
|
carr[0][j]=scratch[0]; carr[0][csize2-j]=conj(scratch[0]);
|
|
carr[0][j+1]=scratchb[0]; carr[0][csize2-1-j]=conj(scratchb[0]);
|
|
for(i=1;i<csize1;i++) {
|
|
carr[i][j]=scratch[i];
|
|
carr[i][j+1]=scratchb[i];
|
|
carr[i][csize2-1-j]=conj(scratchb[vecsize1-i]);
|
|
carr[i][csize2-j]=conj(scratch[vecsize1-i]);
|
|
}
|
|
}
|
|
// There should be 1 column left over
|
|
if(j<csize2/2) {
|
|
for(i=0;i<(rsize1+1)/2;i++) {
|
|
x1=carr[i][j].real()/2; y1=carr[i][j].imag()/2;
|
|
x2=carr[i][csize2-j].real()/2; y2=carr[i][csize2-j].imag()/2;
|
|
scratch[2*i] = MyComplex(x1+x2,y1-y2);
|
|
scratch[(2*i)+1] = MyComplex(y1+y2,x2-x1);
|
|
}
|
|
for(i*=2;i<vecsize1;i++) scratch[i]=MyComplex(0.,0.); // Zero pad
|
|
fft1.ForwardDecFreq(vecsize1,scratch);
|
|
carr[0][j]=scratch[0]; carr[0][csize2-j]=conj(scratch[0]);
|
|
for(i=1;i<csize1;i++) {
|
|
carr[i][j]=scratch[i];
|
|
carr[i][csize2-j]=conj(scratch[vecsize1-i]);
|
|
}
|
|
}
|
|
}
|
|
|
|
void FFTReal2D::InverseRC(int csize1,int csize2,
|
|
const MyComplex* const* carr,
|
|
int rsize1,int rsize2,double** rarr)
|
|
{
|
|
SetupInverse(OC_MAX(1,2*(csize1-1)),csize2); // Safety.
|
|
if(vecsize1<2 || vecsize2<2)
|
|
PlainError(1,"Error in FFTReal2D::InverseRC(...): "
|
|
"Full array dimensions (%dx%d) must be both >=2",
|
|
vecsize1,vecsize2);
|
|
|
|
int i,j;
|
|
FFT_REAL_TYPE x1,y1,x2,y2;
|
|
FFT_REAL_TYPE xb1,yb1,xb2,yb2;
|
|
|
|
// Do row inverse FFT's
|
|
// Handle the first & csize1'th row specially. These rows are
|
|
// the DFT's of real sequences, so they each satisfy the conjugate
|
|
// symmetry condition
|
|
workarr[0][0]=MyComplex(carr[0][0].real(),carr[csize1-1][0].real());
|
|
for(j=1;j<csize2/2;j++) {
|
|
x1=carr[0][j].real(); y1=carr[0][j].imag();
|
|
x2=carr[csize1-1][j].real(); y2=carr[csize1-1][j].imag();
|
|
workarr[0][j] = MyComplex(x1-y2,x2+y1);
|
|
workarr[0][csize2-j] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
workarr[0][csize2/2]=MyComplex(carr[0][csize2/2].real(),
|
|
carr[csize1-1][csize2/2].real());
|
|
fft2.InverseDecTime(csize2,workarr[0],1.);
|
|
|
|
// iFFT the remaining rows
|
|
for(i=1;i<csize1-1;i++) {
|
|
for(j=0;j<csize2;j++) workarr[i][j]=carr[i][j];
|
|
fft2.InverseDecTime(csize2,workarr[i],1.);
|
|
}
|
|
|
|
// Now do iFFT's on columns. These are conj. symmetric, so we
|
|
// process them 2 at a time. Also, recall the 1st row of workarr
|
|
// contains the iFFT's of the 1st and csize1'th row of the given carr.
|
|
for(j=0;j+3<rsize2;j+=4) {
|
|
scratch[0]=
|
|
MyComplex(workarr[0][j].real(),workarr[0][j+1].real());
|
|
scratch[csize1-1]=
|
|
MyComplex(workarr[0][j].imag(),workarr[0][j+1].imag());
|
|
scratchb[0]=
|
|
MyComplex(workarr[0][j+2].real(),workarr[0][j+3].real());
|
|
scratchb[csize1-1]=
|
|
MyComplex(workarr[0][j+2].imag(),workarr[0][j+3].imag());
|
|
for(i=1;i<csize1-1;i++) {
|
|
x1 =workarr[i][j].real(); y1 =workarr[i][j].imag();
|
|
x2 =workarr[i][j+1].real(); y2 =workarr[i][j+1].imag();
|
|
xb1=workarr[i][j+2].real(); yb1=workarr[i][j+2].imag();
|
|
xb2=workarr[i][j+3].real(); yb2=workarr[i][j+3].imag();
|
|
scratch[i] = MyComplex(x1-y2,x2+y1);
|
|
scratch[vecsize1-i] = MyComplex(x1+y2,x2-y1);
|
|
scratchb[i] = MyComplex(xb1-yb2,xb2+yb1);
|
|
scratchb[vecsize1-i] = MyComplex(xb1+yb2,xb2-yb1);
|
|
}
|
|
fft1.InverseDecTime(vecsize1,scratchb,FFT_REAL_TYPE(vecsize1*vecsize2));
|
|
fft1.InverseDecTime(vecsize1,scratch,FFT_REAL_TYPE(vecsize1*vecsize2));
|
|
for(i=0;i<rsize1;i++) {
|
|
rarr[i][j] =scratch[i].real();
|
|
rarr[i][j+1]=scratch[i].imag();
|
|
rarr[i][j+2]=scratchb[i].real();
|
|
rarr[i][j+3]=scratchb[i].imag();
|
|
}
|
|
}
|
|
// Remaining columns if rsize2 is not divisible by 4. OTOH, csize2
|
|
// *is* divisible by 2, so we can assume workarr[i][j+1] exists.
|
|
for(;j<rsize2;j+=2) {
|
|
scratch[0]=
|
|
MyComplex(workarr[0][j].real(),workarr[0][j+1].real());
|
|
scratch[csize1-1]=
|
|
MyComplex(workarr[0][j].imag(),workarr[0][j+1].imag());
|
|
for(i=1;i<csize1-1;i++) {
|
|
x1 =workarr[i][j].real(); y1 =workarr[i][j].imag();
|
|
x2 =workarr[i][j+1].real(); y2 =workarr[i][j+1].imag();
|
|
scratch[i] = MyComplex(x1-y2,x2+y1);
|
|
scratch[vecsize1-i] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
fft1.InverseDecTime(vecsize1,scratch,FFT_REAL_TYPE(vecsize1*vecsize2));
|
|
for(i=0;i<rsize1;i++) {
|
|
rarr[i][j] =scratch[i].real();
|
|
if(j+1<rsize2) rarr[i][j+1]=scratch[i].imag();
|
|
}
|
|
}
|
|
}
|
|
|
|
void FFTReal2D::InverseCR(int csize1,int csize2,
|
|
const MyComplex* const* carr,
|
|
int rsize1,int rsize2,double** rarr)
|
|
{
|
|
SetupInverse(OC_MAX(1,2*(csize1-1)),csize2); // Safety
|
|
if(vecsize1<2 || vecsize2<2)
|
|
PlainError(1,"Error in FFTReal2D::InverseCR(...): "
|
|
"Full array dimensions (%dx%d) must be both >=2",
|
|
vecsize1,vecsize2);
|
|
|
|
int i,j;
|
|
FFT_REAL_TYPE x1,y1,x2,y2,xb1,yb1,xb2,yb2;
|
|
|
|
// Column iFFT's
|
|
// Handle the first & csize2/2'th column specially. These cols are
|
|
// the DFT's of real sequences, so they each satisfy the conjugate
|
|
// symmetry condition
|
|
scratch[0]=MyComplex(carr[0][0].real(),carr[0][csize2/2].real());
|
|
for(i=1;i<csize1-1;i++) {
|
|
x1=carr[i][0].real(); y1=carr[i][0].imag();
|
|
x2=carr[i][csize2/2].real(); y2=carr[i][csize2/2].imag();
|
|
scratch[i] = MyComplex(x1-y2,x2+y1);
|
|
scratch[vecsize1-i] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
scratch[csize1-1]=MyComplex(carr[csize1-1][0].real(),
|
|
carr[csize1-1][csize2/2].real());
|
|
fft1.InverseDecTime(vecsize1,scratch,1);
|
|
for(i=0;i<vecsize1;i+=2) { // ASSUMES vecsize1 is even
|
|
// See packing note below.
|
|
workarr[i/2][0] = MyComplex(scratch[i].real(),scratch[i+1].real());
|
|
workarr[i/2][csize2/2] = MyComplex(scratch[i].imag(),scratch[i+1].imag());
|
|
}
|
|
//
|
|
// Do remaining column iFFT's, two at a time for better memory
|
|
// access locality.
|
|
for(j=1;j+1<csize2/2;j+=2) {
|
|
scratch[0]=carr[0][j];
|
|
scratchb[0]=carr[0][j+1];
|
|
for(i=1;i<csize1-1;i++) {
|
|
scratch[i]=carr[i][j];
|
|
scratchb[i]=carr[i][j+1];
|
|
scratchb[vecsize1-i]=conj(carr[i][csize2-1-j]);
|
|
scratch[vecsize1-i]=conj(carr[i][csize2-j]);
|
|
}
|
|
scratch[csize1-1]=carr[csize1-1][j];
|
|
scratchb[csize1-1]=carr[csize1-1][j+1];
|
|
fft1.InverseDecTime(vecsize1,scratchb,1.);
|
|
fft1.InverseDecTime(vecsize1,scratch,1.);
|
|
// Pack into workarr. Rows will be conjugate symmetric, so we
|
|
// can pack two rows into 1 via r[k]+i.r[k+1] -> workarr[k/2].
|
|
for(i=0;i<rsize1;i+=2) {
|
|
// CAREFUL! The above 'rsize1' bound may depend on how the
|
|
// iFFT's are calculated in the 'Row iFFT's' code section,
|
|
// and how 'i' is initialized.
|
|
x1=scratch[i].real(); y1=scratch[i].imag();
|
|
x2=scratch[i+1].real(); y2=scratch[i+1].imag();
|
|
xb1=scratchb[i].real(); yb1=scratchb[i].imag();
|
|
xb2=scratchb[i+1].real(); yb2=scratchb[i+1].imag();
|
|
workarr[i/2][j] = MyComplex(x1-y2,x2+y1);
|
|
workarr[i/2][j+1] = MyComplex(xb1-yb2,xb2+yb1);
|
|
workarr[i/2][csize2-j-1] = MyComplex(xb1+yb2,xb2-yb1);
|
|
workarr[i/2][csize2-j] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
}
|
|
// There should be 1 column left over
|
|
if((j=(csize2/2)-1)%2==1) {
|
|
// Column (csize2/2)-1 *not* processed above
|
|
scratch[0]=carr[0][j];
|
|
for(i=1;i<csize1-1;i++) {
|
|
scratch[i]=carr[i][j];
|
|
scratch[vecsize1-i]=conj(carr[i][csize2-j]);
|
|
}
|
|
scratch[csize1-1]=carr[csize1-1][j];
|
|
fft1.InverseDecTime(vecsize1,scratch,1);
|
|
for(i=0;i<rsize1;i+=2) {
|
|
// CAREFUL! The above 'rsize1' bound may depend on how the
|
|
// iFFT's are calculated in the 'Row iFFT's' code section,
|
|
// and how 'i' is initialized.
|
|
x1=scratch[i].real(); y1=scratch[i].imag();
|
|
x2=scratch[i+1].real(); y2=scratch[i+1].imag();
|
|
workarr[i/2][j] = MyComplex(x1-y2,x2+y1);
|
|
workarr[i/2][csize2-j] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
}
|
|
|
|
// Row iFFT's
|
|
for(i=0;i<rsize1;i+=2) {
|
|
fft2.InverseDecTime(vecsize2,workarr[i/2],
|
|
FFT_REAL_TYPE(vecsize1*vecsize2));
|
|
for(j=0;j<rsize2;j++) rarr[i][j] = workarr[i/2][j].real();
|
|
if(i+1<rsize1) {
|
|
for(j=0;j<rsize2;j++) rarr[i+1][j] = workarr[i/2][j].imag();
|
|
}
|
|
}
|
|
}
|
|
|
|
void FFTReal2D::Forward1D(int rsize1,int rsize2,
|
|
const double* const* rarr,
|
|
int csize1,int csize2,MyComplex** carr)
|
|
{
|
|
// The ForwardRC/CR routines assume full array dimensions >1.
|
|
// This routine handles the special case where (at least) one of
|
|
// the dimension is 1, which degenerates into a simple 1D FFT.
|
|
if(csize1==1) { // Single row FFT
|
|
int j;
|
|
for(j=0;j<rsize2;j++)
|
|
carr[0][j]=MyComplex(rarr[0][j],0.);
|
|
for(;j<csize2;j++)
|
|
carr[0][j]=MyComplex(0.,0.);
|
|
fft2.ForwardDecFreq(csize2,carr[0]);
|
|
} else if(csize2==1) { // Single column FFT
|
|
int i;
|
|
for(i=0;i<rsize1;i++)
|
|
scratch[i]=MyComplex(rarr[i][0],0.0);
|
|
for(;i<vecsize1;i++)
|
|
scratch[i]=MyComplex(0.0,0.0);
|
|
fft1.ForwardDecFreq(vecsize1,scratch);
|
|
for(i=0;i<csize1;i++)
|
|
carr[i][0]=scratch[i]; // Last half, from csize1 to vecsize1
|
|
/// is conj. sym. since input is real.
|
|
} else {
|
|
PlainError(1,"Error in FFTReal2D::Forward1D(...): "
|
|
"One array dimension (of %dx%d) must be==1",
|
|
vecsize1,vecsize2);
|
|
}
|
|
}
|
|
|
|
void FFTReal2D::Inverse1D(int csize1,int csize2,
|
|
const MyComplex* const* carr,
|
|
int rsize1,int rsize2,double** rarr)
|
|
{
|
|
// The InverseRC/CR routines assume full array dimensions >1.
|
|
// This routine handles the special case where (at least) one of
|
|
// the dimension is 1, which degenerates into a simple 1D FFT.
|
|
if(csize1==1) { // Single row iFFT
|
|
int j;
|
|
for(j=0;j<csize2;j++)
|
|
scratch[j]=carr[0][j];
|
|
fft2.InverseDecTime(csize2,scratch);
|
|
for(j=0;j<rsize2;j++)
|
|
rarr[0][j]=scratch[j].real();
|
|
} else if(csize2==1) { // Single column iFFT
|
|
int i;
|
|
scratch[0]=carr[0][0];
|
|
for(i=1;i<csize1-1;i++) {
|
|
scratch[i]=carr[i][0];
|
|
scratch[vecsize1-i]=conj(carr[i][0]); // Last half obtained
|
|
/// by using conjugate symmetry of real data FFT.
|
|
}
|
|
scratch[csize1-1]=carr[csize1-1][0];
|
|
fft2.InverseDecTime(vecsize1,scratch);
|
|
for(i=0;i<rsize1;i++)
|
|
rarr[i][0]=scratch[i].real();
|
|
} else {
|
|
PlainError(1,"Error in FFTReal2D::Inverse1D(...): "
|
|
"One array dimension (of %dx%d) must be==1",
|
|
vecsize1,vecsize2);
|
|
}
|
|
}
|
|
|
|
void FFTReal2D::Forward(int rsize1,int rsize2,
|
|
const double* const* rarr,
|
|
int csize1,int csize2,MyComplex** carr)
|
|
{
|
|
if(csize2<rsize2)
|
|
PlainError(1,"Error in FFTRealD::Forward(int,int,OC_REAL8,**,int,int,"
|
|
"MyComplex**): csize2 (=%d) *must* be >= rsize2 (=%d)\n",
|
|
csize2,rsize2);
|
|
|
|
if(csize1<(rsize1/2)+1)
|
|
PlainError(1,"Error in FFTRealD::Forward(int,int,double**,int,int,"
|
|
"MyComplex**): csize1 (=%d) *must* be >= (rsize1/2)+1"
|
|
" (=%d)\n",csize1,(rsize1/2)+1);
|
|
|
|
Setup(OC_MAX(1,2*(csize1-1)),csize2);
|
|
if(vecsize1<1 || vecsize2<1) return; // Nothing to do
|
|
|
|
// Check for 1D degenerate cases
|
|
if(vecsize1==1 || vecsize2==1) {
|
|
Forward1D(rsize1,rsize2,rarr,csize1,csize2,carr);
|
|
} else {
|
|
// Determine which Forward routine to call (ForwardCR or ForwardRC)
|
|
// (Use double arithmetic to protect against integer overflow. If
|
|
// the two times are very close, then the choice doesn't really
|
|
// matter.)
|
|
// 1) Estimated (proportional) time for ForwardCR
|
|
REALWIDE crtime=REALWIDE(vecsize1*vecsize2)*REALWIDE(logsize2)
|
|
+ REALWIDE(vecsize1*rsize2)*REALWIDE(logsize1);
|
|
// 2) Estimated (proportional) time for ForwardRC
|
|
REALWIDE rctime=REALWIDE(rsize1*vecsize2)*REALWIDE(logsize2)
|
|
+ REALWIDE(vecsize1*vecsize2)*REALWIDE(logsize1);
|
|
// Introduce empirical adjustment factor
|
|
rctime*=CRRCspeedratio;
|
|
if(crtime<=rctime) ForwardCR(rsize1,rsize2,rarr,csize1,csize2,carr);
|
|
else ForwardRC(rsize1,rsize2,rarr,csize1,csize2,carr);
|
|
}
|
|
}
|
|
|
|
void FFTReal2D::Inverse(int csize1,int csize2,
|
|
const MyComplex* const* carr,
|
|
int rsize1,int rsize2,double** rarr)
|
|
{
|
|
if(csize2<rsize2)
|
|
PlainError(1,"Error in FFTRealD::Inverse(int,int,double**,"
|
|
"int,int,MyComplex**): csize2 (=%d) *must* be >="
|
|
" rsize2 (=%d)\n",csize2,rsize2);
|
|
|
|
if(csize1<(rsize1/2)+1)
|
|
PlainError(1,"Error in FFTRealD::Inverse(int,int,double**,"
|
|
"int,int,MyComplex**): csize1 (=%d) *must* be >="
|
|
" (rsize1/2)+1 (=%d)\n",csize1,(rsize1/2)+1);
|
|
|
|
SetupInverse(OC_MAX(1,2*(csize1-1)),csize2);
|
|
if(vecsize1<1 || vecsize2<1) return; // Nothing to do
|
|
|
|
// Check for 1D degenerate cases
|
|
if(vecsize1==1 || vecsize2==1) {
|
|
Inverse1D(csize1,csize2,carr,rsize1,rsize2,rarr);
|
|
} else {
|
|
// Determine which Inverse routine to call (InverseRC or InverseCR)
|
|
// (Use double arithmetic to protect against integer overflow.
|
|
// If the two times are very close, then the choice doesn't really
|
|
// matter.)
|
|
// 1) Estimated (proportional) time for InverseRC (==ForwardCR)
|
|
REALWIDE irctime=REALWIDE(vecsize1*vecsize2)*REALWIDE(logsize2)
|
|
+ REALWIDE(vecsize1*rsize2)*REALWIDE(logsize1);
|
|
// 2) Estimated (proportional) time for InverseCR (==ForwardRC)
|
|
REALWIDE icrtime=REALWIDE(rsize1*vecsize2)*REALWIDE(logsize2)
|
|
+ REALWIDE(vecsize1*vecsize2)*REALWIDE(logsize1);
|
|
// Introduce empirical adjustment factor
|
|
icrtime*=CRRCspeedratio;
|
|
if(irctime<=icrtime) InverseRC(csize1,csize2,carr,rsize1,rsize2,rarr);
|
|
else InverseCR(csize1,csize2,carr,rsize1,rsize2,rarr);
|
|
}
|
|
}
|
|
|
|
|
|
#ifdef USE_MPI
|
|
|
|
static FFT fft1_mpi,fft2_mpi;
|
|
static int vecsize1_mpi(0),vecsize2_mpi(0);
|
|
static MyComplex* scratch_mpi(NULL);
|
|
static MyComplex** workarr_mpi(NULL);
|
|
|
|
FFTReal2D_mpi::FFTReal2D_mpi()
|
|
{
|
|
}
|
|
|
|
void SetupMemory_mpi(int size1,int size2)
|
|
{
|
|
if(size1!=vecsize1_mpi) {
|
|
if(scratch_mpi!=NULL) delete[] scratch_mpi;
|
|
vecsize1_mpi=size1;
|
|
scratch_mpi=new MyComplex[vecsize1_mpi];
|
|
}
|
|
if(size2!=vecsize2_mpi) {
|
|
if(workarr_mpi!=NULL) {
|
|
delete[] workarr_mpi[0];
|
|
delete[] workarr_mpi;
|
|
}
|
|
vecsize2_mpi=size2;
|
|
int rowcount=(vecsize1_mpi/2)+1;
|
|
workarr_mpi=new MyComplex*[rowcount];
|
|
workarr_mpi[0]=new MyComplex[rowcount*vecsize2_mpi];
|
|
for(int i=1;i<rowcount;i++)
|
|
workarr_mpi[i]=workarr_mpi[i-1]+vecsize2_mpi;
|
|
}
|
|
}
|
|
|
|
void ReleaseMemory_mpi()
|
|
{
|
|
fft1_mpi.ReleaseMemory();
|
|
fft2_mpi.ReleaseMemory();
|
|
if(scratch_mpi!=NULL) delete[] scratch_mpi;
|
|
scratch_mpi=NULL;
|
|
vecsize1_mpi=0;
|
|
if(workarr_mpi!=NULL) {
|
|
delete[] workarr_mpi[0];
|
|
delete[] workarr_mpi;
|
|
}
|
|
workarr_mpi=NULL;
|
|
vecsize2_mpi=0;
|
|
}
|
|
|
|
void FFTReal2D_mpi::ReleaseMemory()
|
|
{
|
|
Mmsolve_MpiWakeUp(ReleaseMemory_mpi); // On slaves
|
|
ReleaseMemory_mpi(); // On master
|
|
}
|
|
|
|
static int vecsize1_mpi_b(0),vecsize2_mpi_b(0);
|
|
static MyComplex **work1_mpi(NULL),**work2_mpi(NULL);
|
|
|
|
static void SetupMemory_mpi_base_b(int size1,int size2)
|
|
{
|
|
int i;
|
|
if(size1!=vecsize1_mpi_b || size2!=vecsize2_mpi_b) {
|
|
if(work1_mpi!=NULL) {
|
|
delete[] work1_mpi[0];
|
|
delete[] work1_mpi;
|
|
}
|
|
if(work2_mpi!=NULL) {
|
|
delete[] work2_mpi[0];
|
|
delete[] work2_mpi;
|
|
}
|
|
vecsize1_mpi_b=size1;
|
|
vecsize2_mpi_b=size2;
|
|
work1_mpi=new MyComplex*[vecsize1_mpi_b];
|
|
work1_mpi[0]=new MyComplex[vecsize1_mpi_b*vecsize2_mpi_b];
|
|
for(i=1;i<vecsize1_mpi_b;i++)
|
|
work1_mpi[i]=work1_mpi[i-1]+vecsize2_mpi_b;
|
|
int rowcount=vecsize2_mpi_b/2;
|
|
work2_mpi=new MyComplex*[rowcount];
|
|
work2_mpi[0]=new MyComplex[rowcount*vecsize1_mpi_b];
|
|
for(i=1;i<rowcount;i++)
|
|
work2_mpi[i]=work2_mpi[i-1]+vecsize1_mpi_b;
|
|
}
|
|
}
|
|
|
|
static void SetupMemory_mpi_slave_b()
|
|
{
|
|
int size1,size2;
|
|
MPI_Bcast(&size1,1,MPI_INT,0,MPI_COMM_WORLD);
|
|
MPI_Bcast(&size2,1,MPI_INT,0,MPI_COMM_WORLD);
|
|
if(size1<1 || size2<1)
|
|
PlainError(1,"Error propagating array size info in "
|
|
"SetupMemory_mpi_slave_b(): size1=%d, size2=%d",
|
|
size1,size2);
|
|
SetupMemory_mpi_base_b(size1,size2);
|
|
}
|
|
|
|
static void SetupMemory_mpi_master_b(int size1,int size2)
|
|
{
|
|
Mmsolve_MpiWakeUp(SetupMemory_mpi_slave_b);
|
|
MPI_Bcast(&size1,1,MPI_INT,0,MPI_COMM_WORLD);
|
|
MPI_Bcast(&size2,1,MPI_INT,0,MPI_COMM_WORLD);
|
|
SetupMemory_mpi_base_b(size1,size2);
|
|
}
|
|
|
|
void ForwardFFT1_mpi_slave_b()
|
|
{
|
|
// Get data
|
|
int rowcount,colcount;
|
|
colcount=vecsize2_mpi_b;
|
|
MPI_Request request[1];
|
|
MPI_Status status[1];
|
|
MPI_Irecv(&rowcount,1,MPI_INT,0,1,MPI_COMM_WORLD,request);
|
|
MPI_Waitall(1,request,status);
|
|
MPI_Irecv(work1_mpi[0],rowcount*colcount,MMS_COMPLEX,0,2,
|
|
MPI_COMM_WORLD,request);
|
|
MPI_Waitall(1,request,status);
|
|
|
|
// Transform rows
|
|
for(int i=0;i<rowcount;i++) {
|
|
fft1_mpi.ForwardDecFreq(colcount,work1_mpi[i]);
|
|
}
|
|
|
|
// Return results
|
|
MPI_Isend(work1_mpi[0],rowcount*colcount,MMS_COMPLEX,0,3,
|
|
MPI_COMM_WORLD,request);
|
|
MPI_Waitall(1,request,status);
|
|
}
|
|
|
|
void ForwardFFT2_mpi_slave_b()
|
|
{
|
|
// Get data
|
|
int rowcount,colcount;
|
|
colcount=vecsize1_mpi_b;
|
|
MPI_Request request[1];
|
|
MPI_Status status[1];
|
|
MPI_Irecv(&rowcount,1,MPI_INT,0,1,MPI_COMM_WORLD,request);
|
|
MPI_Waitall(1,request,status);
|
|
MPI_Irecv(work2_mpi[0],rowcount*colcount,MMS_COMPLEX,0,2,
|
|
MPI_COMM_WORLD,request);
|
|
MPI_Waitall(1,request,status);
|
|
|
|
// Transform rows
|
|
for(int i=0;i<rowcount;i++) {
|
|
fft2_mpi.ForwardDecFreq(colcount,work2_mpi[i]);
|
|
}
|
|
|
|
// Return results
|
|
MPI_Isend(work2_mpi[0],rowcount*colcount,MMS_COMPLEX,0,3,
|
|
MPI_COMM_WORLD,request);
|
|
MPI_Waitall(1,request,status);
|
|
}
|
|
|
|
void FFTReal2D_mpi::Forward(int rsize1,int rsize2,
|
|
const double* const* rarr,
|
|
int csize1,int csize2,MyComplex** carr)
|
|
{ // Computes FFT of rsize1 x rsize2 double** rarr, leaving result in
|
|
// csize1 x csize2 MyComplex** carr. rsize2 *must* be <= csize2, and
|
|
// rsize1 *must* be <= 2*(csize1-1). This routine returns only the
|
|
// top half +1 of the transform. The bottom half is given by
|
|
// carr[2*(csize1-1)-i][csize2-j]=conj(carr[i][j])
|
|
// for i>=csize1, with the second indices interpreted 'mod csize2'.
|
|
// The dimensions csize2 and 2*(csize1-1) *must* be powers of 2,
|
|
// but rsize1 and rsize2 do not. On import, the rarr will be
|
|
// implicitly zero-padded as necessary to fit into the specified
|
|
// output array carr. To get a non-periodic transform, set
|
|
//
|
|
// 2*(csize1-1) to 2*(first power of 2 >= rsize1), and
|
|
// csize2 to 2*(first power of 2 >= rsize2).
|
|
|
|
int i,j;
|
|
FFT_REAL_TYPE x1,y1,x2,y2;
|
|
int vecsize1=OC_MAX(1,2*(csize1-1));
|
|
int vecsize2=csize2;
|
|
for(i=vecsize1;i>2;i/=2)
|
|
if(i%2!=0) PlainError(1,"Error in FFTReal2D_mpi::Forward(int): "
|
|
"Requested csize1 - 1 (%d - 1) is not a power of 2",
|
|
csize1);
|
|
for(j=vecsize2;j>2;j/=2)
|
|
if(j%2!=0) PlainError(1,"Error in FFTReal2D_mpi::Forward(int): "
|
|
"Requested csize2 (%d) is not a power of 2",
|
|
csize2);
|
|
if(vecsize1==0 || vecsize2==0) return; // Nothing to do
|
|
SetupMemory_mpi_master_b(vecsize1,vecsize2);
|
|
|
|
// Copy input data into packed complex array
|
|
for(i=0;i<rsize1/2;i++) {
|
|
for(j=0;j<rsize2;j++)
|
|
work1_mpi[i][j]=MyComplex(rarr[2*i][j],rarr[2*i+1][j]);
|
|
for(;j<vecsize2;j++)
|
|
work1_mpi[i][j]=MyComplex(0,0); // Zero pad
|
|
}
|
|
if(2*i<rsize1) {
|
|
// Odd number of rows. Pack last with zeros.
|
|
for(j=0;j<rsize2;j++)
|
|
work1_mpi[i][j]=MyComplex(rarr[2*i][j],0.0);
|
|
for(;j<vecsize2;j++)
|
|
work1_mpi[i][j]=MyComplex(0,0); // Zero pad
|
|
}
|
|
|
|
// Do FFT across rows
|
|
int proc,proc_count,base_size,base_orphan,whole_size,chunk_size;
|
|
proc_count=mms_mpi_size;
|
|
Mmsolve_MpiWakeUp(ForwardFFT1_mpi_slave_b);
|
|
whole_size=(rsize1+1)/2;
|
|
base_size=whole_size/proc_count;
|
|
base_orphan=whole_size%proc_count;
|
|
i=0;
|
|
chunk_size=base_size+(0<base_orphan?1:0);
|
|
int msg_count=3*(proc_count-1);
|
|
MPI_Request *request=new MPI_Request[msg_count];
|
|
MPI_Status *status=new MPI_Status[msg_count];
|
|
for(proc=1;proc<proc_count;proc++) {
|
|
i+=chunk_size;
|
|
chunk_size=base_size+(proc<base_orphan?1:0);
|
|
MPI_Isend(&chunk_size,1,MPI_INT,proc,1,MPI_COMM_WORLD,
|
|
request+3*(proc-1));
|
|
MPI_Isend(work1_mpi[i],chunk_size*vecsize2,MMS_COMPLEX,proc,2,
|
|
MPI_COMM_WORLD,request+3*(proc-1)+1);
|
|
MPI_Irecv(work1_mpi[i],chunk_size*vecsize2,MMS_COMPLEX,proc,3,
|
|
MPI_COMM_WORLD,request+3*(proc-1)+2);
|
|
}
|
|
chunk_size=base_size+(0<base_orphan?1:0);
|
|
for(i=0;i<chunk_size;i++) {
|
|
fft1_mpi.ForwardDecFreq(vecsize2,work1_mpi[i]);
|
|
}
|
|
MPI_Waitall(msg_count,request,status);
|
|
|
|
// Unpack and transpose to prepare for FFT in cross dimension.
|
|
// We fill the first row of work2_mpi with the first and middle
|
|
// columns from unpacked work1_mpi, because these are real-valued.
|
|
for(i=0;i<(rsize1+1)/2;i++) {
|
|
work2_mpi[0][2*i]
|
|
=MyComplex(work1_mpi[i][0].real(),work1_mpi[i][vecsize2/2].real());
|
|
work2_mpi[0][2*i+1]
|
|
=MyComplex(work1_mpi[i][0].imag(),work1_mpi[i][vecsize2/2].imag());
|
|
}
|
|
for(i=2*i;i<vecsize1;i++) work2_mpi[0][i]=MyComplex(0.,0.); // Zero pad
|
|
// Process rest of the rows.
|
|
for(j=1;j<vecsize2/2;j++) {
|
|
for(i=0;i<(rsize1+1)/2;i++) {
|
|
x1=work1_mpi[i][j].real()/2.;
|
|
y1=work1_mpi[i][j].imag()/2.;
|
|
x2=work1_mpi[i][vecsize2-j].real()/2.;
|
|
y2=work1_mpi[i][vecsize2-j].imag()/2.;
|
|
work2_mpi[j][2*i]=MyComplex(x1+x2,y1-y2);
|
|
work2_mpi[j][2*i+1]=MyComplex(y1+y2,x2-x1);
|
|
}
|
|
for(i=2*i;i<vecsize1;i++) work2_mpi[j][i]=MyComplex(0.,0.); // Zero pad
|
|
}
|
|
|
|
|
|
// Do FFT's on transposed matrix
|
|
Mmsolve_MpiWakeUp(ForwardFFT2_mpi_slave_b);
|
|
whole_size=vecsize2/2;
|
|
base_size=whole_size/proc_count;
|
|
base_orphan=whole_size%proc_count;
|
|
i=0;
|
|
chunk_size=base_size+(0<base_orphan?1:0);
|
|
msg_count=3*(proc_count-1);
|
|
for(proc=1;proc<proc_count;proc++) {
|
|
i+=chunk_size;
|
|
chunk_size=base_size+(proc<base_orphan?1:0);
|
|
MPI_Isend(&chunk_size,1,MPI_INT,proc,1,MPI_COMM_WORLD,
|
|
request+3*(proc-1));
|
|
MPI_Isend(work2_mpi[i],chunk_size*vecsize1,MMS_COMPLEX,proc,2,
|
|
MPI_COMM_WORLD,request+3*(proc-1)+1);
|
|
MPI_Irecv(work2_mpi[i],chunk_size*vecsize1,MMS_COMPLEX,proc,3,
|
|
MPI_COMM_WORLD,request+3*(proc-1)+2);
|
|
}
|
|
chunk_size=base_size+(0<base_orphan?1:0);
|
|
for(i=0;i<chunk_size;i++) {
|
|
fft2_mpi.ForwardDecFreq(vecsize1,work2_mpi[i]);
|
|
}
|
|
MPI_Waitall(msg_count,request,status);
|
|
delete[] request;
|
|
delete[] status;
|
|
|
|
|
|
// Un-transpose and pack into carr.
|
|
// First row of work2_mpi needs to be handled separately, because
|
|
// of unique packing described above.
|
|
carr[0][0]=MyComplex(work2_mpi[0][0].real(),0.);
|
|
carr[0][vecsize2/2]=MyComplex(work2_mpi[0][0].imag(),0.);
|
|
carr[vecsize1/2][0]=MyComplex(work2_mpi[0][vecsize1/2].real(),0.);
|
|
carr[vecsize1/2][vecsize2/2]=MyComplex(work2_mpi[0][vecsize1/2].imag(),0.);
|
|
for(i=1;i<vecsize1/2;i++) {
|
|
x1=work2_mpi[0][i].real()/2.;
|
|
y1=work2_mpi[0][i].imag()/2.;
|
|
x2=work2_mpi[0][vecsize1-i].real()/2.;
|
|
y2=work2_mpi[0][vecsize1-i].imag()/2.;
|
|
carr[i][0]=MyComplex(x1+x2,y1-y2);
|
|
carr[i][vecsize2/2]=MyComplex(y1+y2,x2-x1);
|
|
}
|
|
// Process remaining rows
|
|
for(j=1;j<vecsize2/2;j++) {
|
|
carr[0][j]=work2_mpi[j][0];
|
|
carr[0][csize2-j]=conj(work2_mpi[j][0]);
|
|
for(i=1;i<vecsize1/2;i++)
|
|
carr[i][j]=work2_mpi[j][i];
|
|
carr[i][j]=work2_mpi[j][i];
|
|
carr[i][csize2-j]=conj(work2_mpi[j][i]);
|
|
for(i=csize1;i<vecsize1;i++)
|
|
carr[vecsize1-i][csize2-j]=conj(work2_mpi[j][i]);
|
|
}
|
|
}
|
|
|
|
void FFTReal2D_mpi::Inverse(int csize1,int csize2,
|
|
const MyComplex* const* carr,
|
|
int rsize1,int rsize2,double** rarr)
|
|
{
|
|
// Initialization
|
|
int i,j;
|
|
int vecsize1=OC_MAX(1,2*(csize1-1));
|
|
int vecsize2=csize2;
|
|
FFT_REAL_TYPE x1,y1,x2,y2; // FFT unpacking scratch vars
|
|
for(i=vecsize1;i>2;i/=2)
|
|
if(i%2!=0) PlainError(1,"Error in FFTReal2D_mpi::Inverse(int): "
|
|
"Requested csize1 - 1 (%d - 1) is not a power of 2",
|
|
csize1);
|
|
for(j=vecsize2;j>2;j/=2)
|
|
if(j%2!=0) PlainError(1,"Error in FFTReal2D_mpi::Inverse(int): "
|
|
"Requested csize2 (%d) is not a power of 2",
|
|
csize2);
|
|
if(vecsize1==0 || vecsize2==0) return; // Nothing to do
|
|
SetupMemory_mpi(vecsize1,vecsize2);
|
|
|
|
|
|
// Do row inverse FFT's
|
|
// Handle the first & csize1'th row specially. These rows are
|
|
// the DFT's of real sequences, so they each satisfy the conjugate
|
|
// symmetry condition
|
|
workarr_mpi[0][0]=MyComplex(carr[0][0].real(),carr[csize1-1][0].real());
|
|
for(j=1;j<csize2/2;j++) {
|
|
x1=carr[0][j].real(); y1=carr[0][j].imag();
|
|
x2=carr[csize1-1][j].real(); y2=carr[csize1-1][j].imag();
|
|
workarr_mpi[0][j] = MyComplex(x1-y2,x2+y1);
|
|
workarr_mpi[0][csize2-j] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
workarr_mpi[0][csize2/2]=MyComplex(carr[0][csize2/2].real(),
|
|
carr[csize1-1][csize2/2].real());
|
|
fft2_mpi.InverseDecTime(csize2,workarr_mpi[0],1.);
|
|
|
|
// iFFT the remaining rows
|
|
for(i=1;i<csize1-1;i++) {
|
|
for(j=0;j<csize2;j++) workarr_mpi[i][j]=carr[i][j];
|
|
fft2_mpi.InverseDecTime(csize2,workarr_mpi[i],1.);
|
|
}
|
|
|
|
// Now do iFFT's on columns. These are conj. symmetric, so we
|
|
// process them 2 at a time. Also, recall the 1st row of workarr
|
|
// contains the iFFT's of the 1st and csize1'th row of the given carr.
|
|
// Note that csize is guaranteed divisible by 2, so if rsize2 is odd
|
|
// then rsize2+1<=csize2.
|
|
for(j=0;j<rsize2;j+=2) {
|
|
scratch_mpi[0]=
|
|
MyComplex(workarr_mpi[0][j].real(),workarr_mpi[0][j+1].real());
|
|
scratch_mpi[csize1-1]=
|
|
MyComplex(workarr_mpi[0][j].imag(),workarr_mpi[0][j+1].imag());
|
|
for(i=1;i<csize1-1;i++) {
|
|
x1 =workarr_mpi[i][j].real();
|
|
y1 =workarr_mpi[i][j].imag();
|
|
x2 =workarr_mpi[i][j+1].real();
|
|
y2 =workarr_mpi[i][j+1].imag();
|
|
scratch_mpi[i] = MyComplex(x1-y2,x2+y1);
|
|
scratch_mpi[vecsize1-i] = MyComplex(x1+y2,x2-y1);
|
|
}
|
|
fft1_mpi.InverseDecTime(vecsize1,scratch_mpi,
|
|
FFT_REAL_TYPE(vecsize1*vecsize2));
|
|
if(j+1<rsize2) {
|
|
for(i=0;i<rsize1;i++) {
|
|
rarr[i][j]=scratch_mpi[i].real();
|
|
rarr[i][j+1]=scratch_mpi[i].imag();
|
|
}
|
|
} else {
|
|
for(i=0;i<rsize1;i++) rarr[i][j]=scratch_mpi[i].real();
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
#endif /* USE_MPI */
|
