Magnonics

/* 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 */