--- a/liboctave/dNDArray.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/dNDArray.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -34,10 +34,453 @@
 #include "Array-util.h"
 #include "dNDArray.h"
 #include "mx-base.h"
+#include "f77-fcn.h"
 #include "lo-error.h"
 #include "lo-ieee.h"
 #include "lo-mappers.h"
 
+#if defined (HAVE_FFTW3)
+#include "oct-fftw.h"
+
+ComplexNDArray
+NDArray::fourier (int dim) const
+{
+  dim_vector dv = dims ();
+
+  if (dim > dv.length () || dim < 0)
+    return ComplexNDArray ();
+
+  int stride = 1;
+  int n = dv(dim);
+
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  int howmany = numel () / dv (dim);
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
+  int dist = (stride == 1 ? n : 1);
+
+  const double *in (fortran_vec ());
+  ComplexNDArray retval (dv);
+  Complex *out (retval.fortran_vec ());
+
+  // Need to be careful here about the distance between fft's
+  for (int k = 0; k < nloop; k++)
+    octave_fftw::fft (in + k * stride * n, out + k * stride * n, 
+		      n, howmany, stride, dist);
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::ifourier (const int dim) const
+{
+  dim_vector dv = dims ();
+
+  if (dim > dv.length () || dim < 0)
+    return ComplexNDArray ();
+
+  int stride = 1;
+  int n = dv(dim);
+
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  int howmany = numel () / dv (dim);
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
+  int dist = (stride == 1 ? n : 1);
+
+  ComplexNDArray retval (*this);
+  Complex *out (retval.fortran_vec ());
+
+  // Need to be careful here about the distance between fft's
+  for (int k = 0; k < nloop; k++)
+    octave_fftw::ifft (out + k * stride * n, out + k * stride * n, 
+		      n, howmany, stride, dist);
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::fourier2d (void) const
+{
+  dim_vector dv = dims();
+  if (dv.length () < 2)
+    return ComplexNDArray ();
+
+  dim_vector dv2(dv(0), dv(1));
+  const double *in = fortran_vec ();
+  ComplexNDArray retval (dv);
+  Complex *out = retval.fortran_vec ();
+  int howmany = numel() / dv(0) / dv(1);
+  int dist = dv(0) * dv(1);
+
+  for (int i=0; i < howmany; i++)
+    octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::ifourier2d (void) const
+{
+  dim_vector dv = dims();
+  if (dv.length () < 2)
+    return ComplexNDArray ();
+
+  dim_vector dv2(dv(0), dv(1));
+  ComplexNDArray retval (*this);
+  Complex *out = retval.fortran_vec ();
+  int howmany = numel() / dv(0) / dv(1);
+  int dist = dv(0) * dv(1);
+
+  for (int i=0; i < howmany; i++)
+    octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::fourierNd (void) const
+{
+  dim_vector dv = dims ();
+  int rank = dv.length ();
+
+  const double *in (fortran_vec ());
+  ComplexNDArray retval (dv);
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::fftNd (in, out, rank, dv);
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::ifourierNd (void) const
+{
+  dim_vector dv = dims ();
+  int rank = dv.length ();
+
+  ComplexNDArray tmp (*this);
+  Complex *in (tmp.fortran_vec ());
+  ComplexNDArray retval (dv);
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::ifftNd (in, out, rank, dv);
+
+  return retval;
+}
+
+#else
+
+extern "C"
+{
+  // Note that the original complex fft routines were not written for
+  // double complex arguments.  They have been modified by adding an
+  // implicit double precision (a-h,o-z) statement at the beginning of
+  // each subroutine.
+
+  F77_RET_T
+  F77_FUNC (cffti, CFFTI) (const int&, Complex*);
+
+  F77_RET_T
+  F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
+
+  F77_RET_T
+  F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
+}
+
+ComplexNDArray
+NDArray::fourier (int dim) const
+{
+  dim_vector dv = dims ();
+
+  if (dim > dv.length () || dim < 0)
+    return ComplexNDArray ();
+
+  ComplexNDArray retval (dv);
+  int npts = dv(dim);
+  int nn = 4*npts+15;
+  Array<Complex> wsave (nn);
+  Complex *pwsave = wsave.fortran_vec ();
+
+  OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
+
+  int stride = 1;
+
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  int howmany = numel () / npts;
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+  int dist = (stride == 1 ? npts : 1);
+
+  F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+  for (int k = 0; k < nloop; k++)
+    {
+      for (int j = 0; j < howmany; j++)
+	{
+	  OCTAVE_QUIT;
+
+	  for (int i = 0; i < npts; i++)
+	    tmp[i] = elem((i + k*npts)*stride + j*dist);
+
+	  F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave);
+
+	  for (int i = 0; i < npts; i++)
+	    retval ((i + k*npts)*stride + j*dist) = tmp[i];
+	}
+    }
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::ifourier (int dim) const
+{
+  dim_vector dv = dims ();
+
+  if (dim > dv.length () || dim < 0)
+    return ComplexNDArray ();
+
+  ComplexNDArray retval (dv);
+  int npts = dv(dim);
+  int nn = 4*npts+15;
+  Array<Complex> wsave (nn);
+  Complex *pwsave = wsave.fortran_vec ();
+
+  OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
+
+  int stride = 1;
+
+  for (int i = 0; i < dim; i++)
+    stride *= dv(i);
+
+  int howmany = numel () / npts;
+  howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
+  int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+  int dist = (stride == 1 ? npts : 1);
+
+  F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+  for (int k = 0; k < nloop; k++)
+    {
+      for (int j = 0; j < howmany; j++)
+	{
+	  OCTAVE_QUIT;
+
+	  for (int i = 0; i < npts; i++)
+	    tmp[i] = elem((i + k*npts)*stride + j*dist);
+
+	  F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave);
+
+	  for (int i = 0; i < npts; i++)
+	    retval ((i + k*npts)*stride + j*dist) = tmp[i] / 
+	      static_cast<double> (npts);
+	}
+    }
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::fourier2d (void) const
+{
+  dim_vector dv = dims();
+  dim_vector dv2 (dv(0), dv(1));
+  int rank = 2;
+  ComplexNDArray retval (*this);
+  int stride = 1;
+
+  for (int i = 0; i < rank; i++)
+    {
+      int npts = dv2(i);
+      int nn = 4*npts+15;
+      Array<Complex> wsave (nn);
+      Complex *pwsave = wsave.fortran_vec ();
+      Array<Complex> row (npts);
+      Complex *prow = row.fortran_vec ();
+
+      int howmany = numel () / npts;
+      howmany = (stride == 1 ? howmany : 
+		 (howmany > stride ? stride : howmany));
+      int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+      int dist = (stride == 1 ? npts : 1);
+
+      F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+      for (int k = 0; k < nloop; k++)
+	{
+	  for (int j = 0; j < howmany; j++)
+	    {
+	      OCTAVE_QUIT;
+
+	      for (int l = 0; l < npts; l++)
+		prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+	      F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
+
+	      for (int l = 0; l < npts; l++)
+		retval ((l + k*npts)*stride + j*dist) = prow[l];
+	    }
+	}
+
+      stride *= dv2(i);
+    }
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::ifourier2d (void) const
+{
+  dim_vector dv = dims();
+  dim_vector dv2 (dv(0), dv(1));
+  int rank = 2;
+  ComplexNDArray retval (*this);
+  int stride = 1;
+
+  for (int i = 0; i < rank; i++)
+    {
+      int npts = dv2(i);
+      int nn = 4*npts+15;
+      Array<Complex> wsave (nn);
+      Complex *pwsave = wsave.fortran_vec ();
+      Array<Complex> row (npts);
+      Complex *prow = row.fortran_vec ();
+
+      int howmany = numel () / npts;
+      howmany = (stride == 1 ? howmany : 
+		 (howmany > stride ? stride : howmany));
+      int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+      int dist = (stride == 1 ? npts : 1);
+
+      F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+      for (int k = 0; k < nloop; k++)
+	{
+	  for (int j = 0; j < howmany; j++)
+	    {
+	      OCTAVE_QUIT;
+
+	      for (int l = 0; l < npts; l++)
+		prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+	      F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
+
+	      for (int l = 0; l < npts; l++)
+		retval ((l + k*npts)*stride + j*dist) = prow[l] / 
+		  static_cast<double> (npts);
+	    }
+	}
+
+      stride *= dv2(i);
+    }
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::fourierNd (void) const
+{
+  dim_vector dv = dims ();
+  int rank = dv.length ();
+  ComplexNDArray retval (*this);
+  int stride = 1;
+
+  for (int i = 0; i < rank; i++)
+    {
+      int npts = dv(i);
+      int nn = 4*npts+15;
+      Array<Complex> wsave (nn);
+      Complex *pwsave = wsave.fortran_vec ();
+      Array<Complex> row (npts);
+      Complex *prow = row.fortran_vec ();
+
+      int howmany = numel () / npts;
+      howmany = (stride == 1 ? howmany : 
+		 (howmany > stride ? stride : howmany));
+      int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+      int dist = (stride == 1 ? npts : 1);
+
+      F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+      for (int k = 0; k < nloop; k++)
+	{
+	  for (int j = 0; j < howmany; j++)
+	    {
+	      OCTAVE_QUIT;
+
+	      for (int l = 0; l < npts; l++)
+		prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+	      F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
+
+	      for (int l = 0; l < npts; l++)
+		retval ((l + k*npts)*stride + j*dist) = prow[l];
+	    }
+	}
+
+      stride *= dv(i);
+    }
+
+  return retval;
+}
+
+ComplexNDArray
+NDArray::ifourierNd (void) const
+{
+  dim_vector dv = dims ();
+  int rank = dv.length ();
+  ComplexNDArray retval (*this);
+  int stride = 1;
+
+  for (int i = 0; i < rank; i++)
+    {
+      int npts = dv(i);
+      int nn = 4*npts+15;
+      Array<Complex> wsave (nn);
+      Complex *pwsave = wsave.fortran_vec ();
+      Array<Complex> row (npts);
+      Complex *prow = row.fortran_vec ();
+
+      int howmany = numel () / npts;
+      howmany = (stride == 1 ? howmany : 
+		 (howmany > stride ? stride : howmany));
+      int nloop = (stride == 1 ? 1 : numel () / npts / stride);
+      int dist = (stride == 1 ? npts : 1);
+
+      F77_FUNC (cffti, CFFTI) (npts, pwsave);
+
+      for (int k = 0; k < nloop; k++)
+	{
+	  for (int j = 0; j < howmany; j++)
+	    {
+	      OCTAVE_QUIT;
+
+	      for (int l = 0; l < npts; l++)
+		prow[l] = retval ((l + k*npts)*stride + j*dist);
+
+	      F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
+
+	      for (int l = 0; l < npts; l++)
+		retval ((l + k*npts)*stride + j*dist) = prow[l] /
+		  static_cast<double> (npts);
+	    }
+	}
+
+      stride *= dv(i);
+    }
+
+  return retval;
+}
+
+#endif
+
 NDArray::NDArray (const boolNDArray& a)
   : MArrayN<double> (a.dims ())
 {