Mercurial > octave-nkf

changeset 4773:ccfbd6047a54

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-02-16 19:02:32 by jwe]
author jwe
date Mon, 16 Feb 2004 19:02:33 +0000
parents 9eed17b2c8d1
children 0ff45249d321
files ChangeLog configure.in doc/interpreter/signal.txi liboctave/CMatrix.cc liboctave/CNDArray.cc liboctave/CNDArray.h liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dNDArray.cc liboctave/dNDArray.h liboctave/oct-fftw.cc liboctave/oct-fftw.h src/ChangeLog src/DLD-FUNCTIONS/fft.cc src/DLD-FUNCTIONS/fft2.cc src/DLD-FUNCTIONS/fft_wisdom.cc src/DLD-FUNCTIONS/fftn.cc src/DLD-FUNCTIONS/ifft.cc src/DLD-FUNCTIONS/ifft2.cc src/Makefile.in src/defaults.cc
diffstat 21 files changed, 1922 insertions(+), 486 deletions(-) [+]
line wrap: on
line diff
--- a/ChangeLog	Mon Feb 16 17:57:34 2004 +0000
+++ b/ChangeLog	Mon Feb 16 19:02:33 2004 +0000
@@ -1,3 +1,8 @@
+2004-02-16  David Bateman  <dbateman@free.fr>
+
+	* configure.in: Test for the presence of FFTW 3.x and use it in
+	preference to FFTW 2.x.  Define HAVE_FFTW3
+
 2004-02-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* mkoctfile.in (LINK_DEPS): Include $LIBS and $RLD_FLAG.
--- a/configure.in	Mon Feb 16 17:57:34 2004 +0000
+++ b/configure.in	Mon Feb 16 19:02:33 2004 +0000
@@ -29,7 +29,7 @@
 EXTERN_CXXFLAGS="$CXXFLAGS"
 
 AC_INIT
-AC_REVISION($Revision: 1.445 $)
+AC_REVISION($Revision: 1.446 $)
 AC_PREREQ(2.57)
 AC_CONFIG_SRCDIR([src/octave.cc])
 AC_CONFIG_HEADER(config.h)
@@ -407,19 +407,17 @@
   with_fftw=$withval, with_fftw=yes)
 
 if test "$with_fftw" = "yes"; then
-  have_fftw_header=no
-  AC_CHECK_HEADERS(dfftw.h fftw.h, [have_fftw_header=yes; break])
-  if test "$have_fftw_header" = yes; then
-    AC_CHECK_LIB(dfftw, fftw_create_plan, FFTW_LIBS="-ldfftw",
-      [AC_CHECK_LIB(fftw, fftw_create_plan, FFTW_LIBS="-lfftw", with_fftw=no)])
-  else
-    with_fftw=no
+  have_fftw3_header=no
+  with_fftw3=no
+  AC_CHECK_HEADER(fftw3.h, [have_fftw3_header=yes; break])
+  if test "$have_fftw3_header" = yes; then
+      AC_CHECK_LIB(fftw3, fftw_plan_dft_1d, [FFTW_LIBS="-lfftw3"; with_fftw3=yes])
   fi
 fi
 
-if test "$with_fftw" = yes; then
+if test "$with_fftw" = yes && test "$with_fftw3" = yes; then
   FFT_DIR=''
-  AC_DEFINE(HAVE_FFTW, 1, [Define if the FFTW library is available.])
+  AC_DEFINE(HAVE_FFTW3, 1, [Define if the FFTW3 library is used.])
 fi
 
 WITH_MPI=true
--- a/doc/interpreter/signal.txi	Mon Feb 16 17:57:34 2004 +0000
+++ b/doc/interpreter/signal.txi	Mon Feb 16 19:02:33 2004 +0000
@@ -19,6 +19,10 @@
 
 @DOCSTRING(ifft2)
 
+@DOCSTRING(fftn)
+
+@DOCSTRING(ifftn)
+
 @DOCSTRING(fftconv)
 
 @DOCSTRING(fftfilt)
--- a/liboctave/CMatrix.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/CMatrix.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -56,7 +56,7 @@
 #include "mx-inlines.cc"
 #include "oct-cmplx.h"
 
-#ifdef HAVE_FFTW
+#if defined (HAVE_FFTW3)
 #include "oct-fftw.h"
 #endif
 
@@ -1102,7 +1102,7 @@
   return retval;
 }
 
-#ifdef HAVE_FFTW
+#if defined (HAVE_FFTW3)
 
 ComplexMatrix
 ComplexMatrix::fourier (void) const
@@ -1128,12 +1128,7 @@
   const Complex *in (data ());
   Complex *out (retval.fortran_vec ());
 
-  for (size_t i = 0; i < nsamples; i++)
-    {
-      OCTAVE_QUIT;
-
-      octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
-    }
+  octave_fftw::fft (in, out, npts, nsamples); 
 
   return retval;
 }
@@ -1162,12 +1157,7 @@
   const Complex *in (data ());
   Complex *out (retval.fortran_vec ());
 
-  for (size_t i = 0; i < nsamples; i++)
-    {
-      OCTAVE_QUIT;
-
-      octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
-    }
+  octave_fftw::ifft (in, out, npts, nsamples); 
 
   return retval;
 }
@@ -1175,13 +1165,13 @@
 ComplexMatrix
 ComplexMatrix::fourier2d (void) const
 {
-  int nr = rows ();
-  int nc = cols ();
-
-  ComplexMatrix retval (*this);
-  // Note the order of passing the rows and columns to account for
-  // column-major storage.
-  octave_fftw::fft2d (retval.fortran_vec (), nc, nr);
+  dim_vector dv(rows (), cols ());
+
+  ComplexMatrix retval (rows (), cols ());
+  const Complex *in (data ());
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::fftNd (in, out, 2, dv);
 
   return retval;
 }
@@ -1189,13 +1179,13 @@
 ComplexMatrix
 ComplexMatrix::ifourier2d (void) const
 {
-  int nr = rows ();
-  int nc = cols ();
-
-  ComplexMatrix retval (*this);
-  // Note the order of passing the rows and columns to account for
-  // column-major storage.
-  octave_fftw::ifft2d (retval.fortran_vec (), nc, nr);
+  dim_vector dv(rows (), cols ());
+
+  ComplexMatrix retval (rows (), cols ());
+  const Complex *in (data ());
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::ifftNd (in, out, 2, dv);
 
   return retval;
 }
@@ -1332,8 +1322,8 @@
   wsave.resize (nn);
   pwsave = wsave.fortran_vec ();
 
-  Array<Complex> row (npts);
-  Complex *prow = row.fortran_vec ();
+  Array<Complex> tmp (npts);
+  Complex *prow = tmp.fortran_vec ();
 
   F77_FUNC (cffti, CFFTI) (npts, pwsave);
 
@@ -1401,8 +1391,8 @@
   wsave.resize (nn);
   pwsave = wsave.fortran_vec ();
 
-  Array<Complex> row (npts);
-  Complex *prow = row.fortran_vec ();
+  Array<Complex> tmp (npts);
+  Complex *prow = tmp.fortran_vec ();
 
   F77_FUNC (cffti, CFFTI) (npts, pwsave);
 
--- a/liboctave/CNDArray.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/CNDArray.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -34,9 +34,31 @@
 #include "Array-util.h"
 #include "CNDArray.h"
 #include "mx-base.h"
+#include "f77-fcn.h"
 #include "lo-ieee.h"
 #include "lo-mappers.h"
 
+#if defined (HAVE_FFTW3)
+#  include "oct-fftw.h"
+#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*);
+}
+#endif
+
 // XXX FIXME XXX -- could we use a templated mixed-type copy function
 // here?
 
@@ -61,6 +83,429 @@
     elem (i) = a.elem (i);
 }
 
+#if defined (HAVE_FFTW3)
+ComplexNDArray
+ComplexNDArray::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 Complex *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
+ComplexNDArray::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);
+
+  const Complex *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::ifft (in + k * stride * n, out + k * stride * n, 
+		      n, howmany, stride, dist);
+
+  return retval;
+}
+
+ComplexNDArray
+ComplexNDArray::fourier2d (void) const
+{
+  dim_vector dv = dims();
+  if (dv.length () < 2)
+    return ComplexNDArray ();
+
+  dim_vector dv2(dv(0), dv(1));
+  const Complex *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
+ComplexNDArray::ifourier2d (void) const
+{
+  dim_vector dv = dims();
+  if (dv.length () < 2)
+    return ComplexNDArray ();
+
+  dim_vector dv2(dv(0), dv(1));
+  const Complex *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::ifftNd (in + i*dist, out + i*dist, 2, dv2);
+
+  return retval;
+}
+
+ComplexNDArray
+ComplexNDArray::fourierNd (void) const
+{
+  dim_vector dv = dims ();
+  int rank = dv.length ();
+
+  const Complex *in (fortran_vec ());
+  ComplexNDArray retval (dv);
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::fftNd (in, out, rank, dv);
+
+  return retval;
+}
+
+ComplexNDArray
+ComplexNDArray::ifourierNd (void) const
+{
+  dim_vector dv = dims ();
+  int rank = dv.length ();
+
+  const Complex *in (fortran_vec ());
+  ComplexNDArray retval (dv);
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::ifftNd (in, out, rank, dv);
+
+  return retval;
+}
+
+#else
+ComplexNDArray
+ComplexNDArray::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
+ComplexNDArray::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
+ComplexNDArray::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
+ComplexNDArray::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
+ComplexNDArray::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
+ComplexNDArray::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
+
 // unary operations
 
 boolNDArray
--- a/liboctave/CNDArray.h	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/CNDArray.h	Mon Feb 16 19:02:33 2004 +0000
@@ -94,6 +94,15 @@
   
   NDArray abs (void) const;
 
+  ComplexNDArray fourier (int dim = 1) const;
+  ComplexNDArray ifourier (int dim = 1) const;
+
+  ComplexNDArray fourier2d (void) const;
+  ComplexNDArray ifourier2d (void) const;
+
+  ComplexNDArray fourierNd (void) const;
+  ComplexNDArray ifourierNd (void) const;
+
   ComplexMatrix matrix_value (void) const;
 
   ComplexNDArray squeeze (void) const { return ArrayN<Complex>::squeeze (); }
--- a/liboctave/ChangeLog	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/ChangeLog	Mon Feb 16 19:02:33 2004 +0000
@@ -1,3 +1,31 @@
+2004-02-16  David Bateman <dbateman@free.fr>
+
+	* oct-fftw.cc (octave_fftw_planner::create_plan, octave_fftw::fftNd):
+	Add support for FFTW 3.x. Include the ability to
+	use the real to complex transform for fft's of real matrices
+	(octave_fftw_planner::create_plan2d): Delete.
+	(octave_fftw::fft2d): Delete.
+	(convert_packcomplex_1d, convert_packcomplex_Nd):
+	New static functions.
+	* oct-fftw.h: Update decls.
+
+	* dMatrix.cc (Matrix::fourier, Matrix::ifourier,
+	Matrix::fourier2d, Matrix::ifourier2d): FFT's use real to complex
+	transforms.  1D FFT of a matrix done as single call rather than
+	loop.  Update for FFTW 3.x
+	* CMatrix.cc (ComplexMatrix::fourier, ComplexMatrix::ifourier,
+	ComplexMatrix::fourier2d, ComplexMatrix::ifourier2d): 1D fft of a
+	matrix done as single call rather than loop.  Update for FFTW 3.x.
+
+	* dNDArray.cc (NDArray::fourier, NDArray::ifourier,
+	NDArray::fourierNd, NDArray::ifouriourNd): New fourier transform
+	functions for Nd arrays.
+	* dNArray.h Provide decls.
+	* CNDArray.cc (ComplexNDArray::fourier, ComplexNDArray::ifourier,
+	ComplexNDArray::fourierNd, ComplexNDArray::ifouriourNd): New
+	fourier transform functions for complex Nd arrays.
+	* CNArray.h: Provide decls.
+	
 2004-02-15  Petter Risholm  <risholm@stud.ntnu.no>
 
 	* Array.cc (Array<T>::insert (const Array<T>&, int, int)):
--- a/liboctave/dMatrix.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/dMatrix.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -52,7 +52,7 @@
 #include "oct-cmplx.h"
 #include "quit.h"
 
-#ifdef HAVE_FFTW
+#if defined (HAVE_FFTW3)
 #include "oct-fftw.h"
 #endif
 
@@ -766,7 +766,7 @@
     }
 }
 
-#ifdef HAVE_FFTW
+#if defined (HAVE_FFTW3)
 
 ComplexMatrix
 Matrix::fourier (void) const
@@ -789,16 +789,10 @@
       nsamples = nc;
     }
 
-  ComplexMatrix tmp (*this);
-  Complex *in (tmp.fortran_vec ());
+  const double *in (fortran_vec ());
   Complex *out (retval.fortran_vec ());
 
-  for (size_t i = 0; i < nsamples; i++)
-    {
-      OCTAVE_QUIT;
-
-      octave_fftw::fft (&in[npts * i], &out[npts * i], npts);
-    }
+  octave_fftw::fft (in, out, npts, nsamples); 
 
   return retval;
 }
@@ -828,12 +822,7 @@
   Complex *in (tmp.fortran_vec ());
   Complex *out (retval.fortran_vec ());
 
-  for (size_t i = 0; i < nsamples; i++)
-    {
-      OCTAVE_QUIT;
-
-      octave_fftw::ifft (&in[npts * i], &out[npts * i], npts);
-    }
+  octave_fftw::ifft (in, out, npts, nsamples); 
 
   return retval;
 }
@@ -841,13 +830,11 @@
 ComplexMatrix
 Matrix::fourier2d (void) const
 {
-  int nr = rows ();
-  int nc = cols ();
-
-  ComplexMatrix retval (*this);
-  // Note the order of passing the rows and columns to account for
-  // column-major storage.
-  octave_fftw::fft2d (retval.fortran_vec (), nc, nr);
+  dim_vector dv(rows (), cols ());
+
+  const double *in = fortran_vec ();
+  ComplexMatrix retval (rows (), cols ());
+  octave_fftw::fftNd (in, retval.fortran_vec (), 2, dv);
 
   return retval;
 }
@@ -855,13 +842,12 @@
 ComplexMatrix
 Matrix::ifourier2d (void) const
 {
-  int nr = rows ();
-  int nc = cols ();
+  dim_vector dv(rows (), cols ());
 
   ComplexMatrix retval (*this);
-  // Note the order of passing the rows and columns to account for
-  // column-major storage.
-  octave_fftw::ifft2d (retval.fortran_vec (), nc, nr);
+  Complex *out (retval.fortran_vec ());
+
+  octave_fftw::ifftNd (out, out, 2, dv);
 
   return retval;
 }
@@ -998,8 +984,8 @@
   wsave.resize (nn);
   pwsave = wsave.fortran_vec ();
 
-  Array<Complex> row (npts);
-  Complex *prow = row.fortran_vec ();
+  Array<Complex> tmp (npts);
+  Complex *prow = tmp.fortran_vec ();
 
   F77_FUNC (cffti, CFFTI) (npts, pwsave);
 
@@ -1067,8 +1053,8 @@
   wsave.resize (nn);
   pwsave = wsave.fortran_vec ();
 
-  Array<Complex> row (npts);
-  Complex *prow = row.fortran_vec ();
+  Array<Complex> tmp (npts);
+  Complex *prow = tmp.fortran_vec ();
 
   F77_FUNC (cffti, CFFTI) (npts, pwsave);
 
--- 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 ())
 {
--- a/liboctave/dNDArray.h	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/dNDArray.h	Mon Feb 16 19:02:33 2004 +0000
@@ -90,6 +90,15 @@
      
   NDArray abs (void) const;
 
+  ComplexNDArray fourier (int dim = 1) const;
+  ComplexNDArray ifourier (int dim = 1) const;
+
+  ComplexNDArray fourier2d (void) const;
+  ComplexNDArray ifourier2d (void) const;
+
+  ComplexNDArray fourierNd (void) const;
+  ComplexNDArray ifourierNd (void) const;
+
   friend NDArray real (const ComplexNDArray& a);
   friend NDArray imag (const ComplexNDArray& a);
 
--- a/liboctave/oct-fftw.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/oct-fftw.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -22,18 +22,26 @@
 #include <config.h>
 #endif
 
-#ifdef HAVE_FFTW
+#if defined (HAVE_FFTW3)
 
 #include "oct-fftw.h"
 #include "lo-error.h"
-
+#include <iostream>
 
 // Helper class to create and cache fftw plans for both 1d and 2d. This
 // implementation uses FFTW_ESTIMATE to create the plans, which in theory
-// is suboptimal, but provides quite reasonable performance. Future
-// enhancement will be to add a dynamically loadable interface ("fftw")
-// to manipulate fftw wisdom so that users may choose the appropriate
-// planner.
+// is suboptimal, but provides quite reasonable performance.
+
+// Also note that if FFTW_ESTIMATE is not used the planner in FFTW3
+// destroys the input and output arrays. So with the form of the 
+// current code we definitely want FFTW_ESTIMATE!! However, we use
+// any wsidom that is available, either in a FFTW3 system wide file
+// or as supplied by the user.
+
+// XXX FIXME XXX If we can ensure 16 byte alignment in Array<T> (<T> *data)
+// the FFTW3 can use SIMD instructions for further acceleration.
+
+// Note that it is profitable to store the FFTW3 plans, for small ffts
 
 class
 octave_fftw_planner
@@ -41,17 +49,35 @@
 public:
   octave_fftw_planner ();
 
-  fftw_plan create_plan (fftw_direction, size_t);
-  fftwnd_plan create_plan2d (fftw_direction, size_t, size_t);
+  fftw_plan create_plan (int dir, const int rank, const dim_vector dims, 
+			 int howmany, int stride, int dist, 
+			 const Complex *in, Complex *out);
+  fftw_plan create_plan (const int rank, const dim_vector dims, 
+			 int howmany, int stride, int dist, 
+			 const double *in, Complex *out);
 
 private:
   int plan_flags;
 
+  // Plan for fft and ifft of complex values
   fftw_plan plan[2];
-  fftwnd_plan plan2d[2];
+  int d[2];  // dist
+  int s[2];  // stride
+  int r[2];  // rank
+  int h[2];  // howmany
+  dim_vector n[2]; // dims
+  char ialign[2];
+  char oalign[2];
 
-  size_t n[2];
-  size_t n2d[2][2];
+  // Plan for fft of real values
+  fftw_plan rplan;
+  int rd;  // dist
+  int rs;  // stride
+  int rr;  // rank
+  int rh;  // howmany
+  dim_vector rn; // dims
+  char rialign;
+  char roalign;
 };
 
 octave_fftw_planner::octave_fftw_planner ()
@@ -59,31 +85,67 @@
   plan_flags = FFTW_ESTIMATE;
 
   plan[0] = plan[1] = 0;
-  plan2d[0] = plan2d[1] = 0;
-  
-  n[0] = n[1] = 0;
-  n2d[0][0] = n2d[0][1] = n2d[1][0] = n2d[1][1] = 0;
+  d[0] = d[1] = s[0] = s[1] = r[0] = r[1] = h[0] = h[1] = 0;
+  ialign[0] = ialign[1] = oalign[0] = oalign[1] = 0;
+  n[0] = n[1] = dim_vector();
+
+  rplan = 0;
+  rd = rs = rr = rh = 0;
+  rialign = roalign = 0;
+  rn = dim_vector ();
+
+  // If we have a system wide wisdom file, import it
+  fftw_import_system_wisdom ( );
 }
 
 fftw_plan
-octave_fftw_planner::create_plan (fftw_direction dir, size_t npts)
+octave_fftw_planner::create_plan (int dir, const int rank,
+				  const dim_vector dims, int howmany,
+				  int stride, int dist, 
+				  const Complex *in, Complex *out)
 {
-  size_t which = (dir == FFTW_FORWARD) ? 0 : 1;
+  int which = (dir == FFTW_FORWARD) ? 0 : 1;
   fftw_plan *cur_plan_p = &plan[which];
   bool create_new_plan = false;
+  char in_align = ((int) in) & 0xF;
+  char out_align = ((int) out) & 0xF;
 
-  if (plan[which] == 0 || n[which] != npts)
-    {
-      create_new_plan = true;
-      n[which] = npts;
-    }
+  if (plan[which] == 0 || d[which] != dist || s[which] != stride ||
+      r[which] != rank || h[which] != howmany ||
+      ialign[which] != in_align || oalign[which] != out_align)
+    create_new_plan = true;
+  else
+    // We still might not have the same shape of array
+    for (int i = 0; i < rank; i++)
+      if (dims(i) != n[which](i))
+	{
+	  create_new_plan = true;
+	  break;
+	}
 
   if (create_new_plan)
     {
+      d[which] = dist;
+      s[which] = stride;
+      r[which] = rank;
+      h[which] = howmany;
+      ialign[which] = in_align;
+      oalign[which] = out_align;
+      n[which] = dims;
+
       if (*cur_plan_p)
 	fftw_destroy_plan (*cur_plan_p);
 
-      *cur_plan_p = fftw_create_plan (npts, dir, plan_flags);
+      // Note reversal of dimensions for column major storage in FFTW
+      OCTAVE_LOCAL_BUFFER (int, tmp, rank);
+      for (int i = 0, j = rank-1; i < rank; i++, j--)
+	tmp[i] = dims(j);
+
+      *cur_plan_p =
+	fftw_plan_many_dft (rank, tmp, howmany,
+	      reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in)),
+	      NULL, stride, dist, reinterpret_cast<fftw_complex *> (out),
+	      NULL, stride, dist, dir, plan_flags);
 
       if (*cur_plan_p == 0)
 	(*current_liboctave_error_handler) ("Error creating fftw plan");
@@ -92,32 +154,55 @@
   return *cur_plan_p;
 }
  
-fftwnd_plan
-octave_fftw_planner::create_plan2d (fftw_direction dir, 
-                                   size_t nrows, size_t ncols)
+fftw_plan
+octave_fftw_planner::create_plan (const int rank, const dim_vector dims, 
+				  int howmany, int stride, int dist, 
+				  const double *in, Complex *out)
 {
-  size_t which = (dir == FFTW_FORWARD) ? 0 : 1;
-  fftwnd_plan *cur_plan_p = &plan2d[which];
+  fftw_plan *cur_plan_p = &rplan;
   bool create_new_plan = false;
+  char in_align = ((int) in) & 0xF;
+  char out_align = ((int) out) & 0xF;
 
-  if (plan2d[which] == 0 || n2d[which][0] != nrows || n2d[which][1] != ncols)
-    {
-      create_new_plan = true;
-
-      n2d[which][0] = nrows;
-      n2d[which][1] = ncols;
-    }
+  if (rplan == 0 || rd != dist || rs != stride ||
+      rr != rank || rh != howmany ||
+      rialign != in_align || roalign != out_align)
+    create_new_plan = true;
+  else
+    // We still might not have the same shape of array
+    for (int i = 0; i < rank; i++)
+      if (dims(i) != rn(i))
+	{
+	  create_new_plan = true;
+	  break;
+	}
 
   if (create_new_plan)
     {
+      rd = dist;
+      rs = stride;
+      rr = rank;
+      rh = howmany;
+      rialign = in_align;
+      roalign = out_align;
+      rn = dims;
+
       if (*cur_plan_p)
-	fftwnd_destroy_plan (*cur_plan_p);
+	fftw_destroy_plan (*cur_plan_p);
 
-      *cur_plan_p = fftw2d_create_plan (nrows, ncols, dir, 
-					plan_flags | FFTW_IN_PLACE);
+      // Note reversal of dimensions for column major storage in FFTW
+      OCTAVE_LOCAL_BUFFER (int, tmp, rank);
+      for (int i = 0, j = rank-1; i < rank; i++, j--)
+	tmp[i] = dims(j);
+
+      *cur_plan_p =
+	fftw_plan_many_dft_r2c (rank, tmp, howmany,
+	      (const_cast<double *> (in)),
+	      NULL, stride, dist, reinterpret_cast<fftw_complex *> (out),
+	      NULL, stride, dist, plan_flags);
 
       if (*cur_plan_p == 0)
-	(*current_liboctave_error_handler) ("Error creating 2d fftw plan");
+	(*current_liboctave_error_handler) ("Error creating fftw plan");
     }
 
   return *cur_plan_p;
@@ -125,51 +210,184 @@
 
 static octave_fftw_planner fftw_planner;
 
-int
-octave_fftw::fft (const Complex *in, Complex *out, size_t npts)
+static inline void convert_packcomplex_1d (Complex *out, size_t nr, 
+					   size_t nc, int stride, int dist)
+{
+  // Fill in the missing data
+  for (size_t i = 0; i < nr; i++)
+    for (size_t j = nc/2+1; j < nc; j++)
+      out[j*stride + i*dist] = conj(out[(nc - j)*stride + i*dist]);
+}
+
+static inline void convert_packcomplex_Nd (Complex *out, 
+					   const dim_vector &dv)
 {
-  fftw_one (fftw_planner.create_plan (FFTW_FORWARD, npts),
-            reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in)),
-            reinterpret_cast<fftw_complex *> (out));
+  size_t nc = dv(0);
+  size_t nr = dv(1);
+  size_t np = (dv.length() > 2 ? dv.numel () / nc / nr : 1);
+  size_t nrp = nr * np;
+  Complex *ptr1, *ptr2;
+
+  // Create space for the missing elements
+  for (size_t i = 0; i < nrp; i++)
+    {
+      ptr1 = out + i * (nc/2 + 1) + nrp*((nc-1)/2);
+      ptr2 = out + i * nc;
+      for (size_t j = 0; j < nc/2+1; j++)
+	*ptr2++ = *ptr1++;
+    }
+
+  // Fill in the missing data for the rank = 2 case directly for speed
+  for (size_t i = 0; i < np; i++)
+    {
+      for (size_t j = 1; j < nr; j++)
+	for (size_t k = nc/2+1; k < nc; k++)
+	  out[k + (j + i*nr)*nc] = conj(out[nc - k + ((i+1)*nr - j)*nc]);
+
+      for (size_t j = nc/2+1; j < nc; j++)
+	out[j + i*nr*nc] = conj(out[(i*nr+1)*nc - j]);
+    }
+
+  // Now do the permutations needed for rank > 2 cases
+  size_t jstart = dv(0) * dv(1);
+  size_t kstep = dv(0);
+  size_t nel = dv.numel ();
+  for (int inner = 2; inner < dv.length(); inner++) 
+    {
+      size_t jmax = jstart * dv(inner);
+      for (size_t i = 0; i < nel; i+=jmax)
+	for (size_t j = jstart, jj = jmax-jstart; j < jj; 
+	     j+=jstart, jj-=jstart)
+	  for (size_t k = 0; k < jstart; k+= kstep)
+	    for (size_t l = nc/2+1; l < nc; l++)
+	      {
+		Complex tmp = out[i+ j + k + l];
+		out[i + j + k + l] =  out[i + jj + k + l];
+		out[i + jj + k + l] = tmp;
+	      }
+      jstart = jmax;
+    }
+}
+
+int
+octave_fftw::fft (const double *in, Complex *out, size_t npts, 
+		  size_t nsamples, int stride, int dist)
+{
+  dist = (dist < 0 ? npts : dist);
+
+  dim_vector dv (npts);
+  fftw_plan plan = fftw_planner.create_plan (1, dv, nsamples, stride, dist,
+					     in, out);
+
+  fftw_execute_dft_r2c (plan, (const_cast<double *>(in)),
+			reinterpret_cast<fftw_complex *> (out));
+
+  // Need to create other half of the transform
+  convert_packcomplex_1d (out, nsamples, npts, stride, dist);
 
   return 0;
 }
 
 int
-octave_fftw::ifft (const Complex *in, Complex *out, size_t npts)
+octave_fftw::fft (const Complex *in, Complex *out, size_t npts, 
+		  size_t nsamples, int stride, int dist)
 {
-  fftw_one (fftw_planner.create_plan (FFTW_BACKWARD, npts),
-            reinterpret_cast<fftw_complex *> (const_cast<Complex *> (in)),
-            reinterpret_cast<fftw_complex *> (out));
+  dist = (dist < 0 ? npts : dist);
+
+  dim_vector dv (npts);
+  fftw_plan plan = fftw_planner.create_plan (FFTW_FORWARD, 1, dv, nsamples,
+					     stride, dist, in, out);
+
+  fftw_execute_dft (plan, 
+	reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
+	reinterpret_cast<fftw_complex *> (out));
+
+  return 0;
+}
+
+int
+octave_fftw::ifft (const Complex *in, Complex *out, size_t npts, 
+		   size_t nsamples, int stride, int dist)
+{
+  dist = (dist < 0 ? npts : dist);
+
+  dim_vector dv (npts);
+  fftw_plan plan = fftw_planner.create_plan (FFTW_BACKWARD, 1, dv, nsamples,
+					     stride, dist, in, out);
+
+  fftw_execute_dft (plan, 
+	reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
+	reinterpret_cast<fftw_complex *> (out));
 
   const Complex scale = npts;
-  for (size_t i = 0; i < npts; i++)
-    out[i] /= scale;
+  for (size_t j = 0; j < nsamples; j++)
+    for (size_t i = 0; i < npts; i++)
+      out[i*stride + j*dist] /= scale;
 
   return 0;
 }
 
 int
-octave_fftw::fft2d (Complex *inout, size_t nr, size_t nc)
+octave_fftw::fftNd (const double *in, Complex *out, const int rank, 
+		    const dim_vector &dv)
 {
-  fftwnd_one (fftw_planner.create_plan2d (FFTW_FORWARD, nr, nc),
-              reinterpret_cast<fftw_complex *> (inout),
-              0);
+  int dist = 1;
+  for (int i = 0; i < rank; i++)
+    dist *= dv(i);
+
+  // Fool with the position of the start of the output matrix, so that
+  // creating other half of the matrix won't cause cache problems
+  int offset = (dv.numel () / dv(0)) * ((dv(0) - 1) / 2); 
+  
+  fftw_plan plan = fftw_planner.create_plan (rank, dv, 1, 1, dist,
+					     in, out + offset);
+
+  fftw_execute_dft_r2c (plan, (const_cast<double *>(in)),
+			reinterpret_cast<fftw_complex *> (out+ offset));
+
+  // Need to create other half of the transform
+  convert_packcomplex_Nd (out, dv);
 
   return 0;
 }
 
 int
-octave_fftw::ifft2d (Complex *inout, size_t nr, size_t nc)
+octave_fftw::fftNd (const Complex *in, Complex *out, const int rank, 
+		    const dim_vector &dv)
 {
-  fftwnd_one (fftw_planner.create_plan2d (FFTW_BACKWARD, nr, nc),
-              reinterpret_cast<fftw_complex *> (inout),
-              0);
+  int dist = 1;
+  for (int i = 0; i < rank; i++)
+    dist *= dv(i);
+
+  fftw_plan plan = fftw_planner.create_plan (FFTW_FORWARD, rank, dv, 1, 1,
+					     dist, in, out);
+
+  fftw_execute_dft (plan, 
+	reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
+	reinterpret_cast<fftw_complex *> (out));
+
+  return 0;
+}
 
-  const size_t npts = nr * nc;
+int
+octave_fftw::ifftNd (const Complex *in, Complex *out, const int rank, 
+		    const dim_vector &dv)
+{
+  int dist = 1;
+  for (int i = 0; i < rank; i++)
+    dist *= dv(i);
+
+  fftw_plan plan = fftw_planner.create_plan (FFTW_BACKWARD, rank, dv, 1, 1,
+					     dist, in, out);
+
+  fftw_execute_dft (plan, 
+	reinterpret_cast<fftw_complex *> (const_cast<Complex *>(in)),
+	reinterpret_cast<fftw_complex *> (out));
+
+  const size_t npts = dv.numel ();
   const Complex scale = npts;
   for (size_t i = 0; i < npts; i++)
-    inout[i] /= scale;
+    out[i] /= scale;
 
   return 0;
 }
--- a/liboctave/oct-fftw.h	Mon Feb 16 17:57:34 2004 +0000
+++ b/liboctave/oct-fftw.h	Mon Feb 16 19:02:33 2004 +0000
@@ -22,24 +22,27 @@
 #define octave_oct_fftw_h 1
 
 #include <cstddef>
-
-#if defined (HAVE_DFFTW_H)
-#include <dfftw.h>
-#else
-#include <fftw.h>
-#endif
+#include <fftw3.h>
 
 #include "oct-cmplx.h"
+#include "dim-vector.h"
 
 class
 octave_fftw
 {
 public:
-  static int fft (const Complex*, Complex *, size_t);
-  static int ifft (const Complex*, Complex *, size_t);
+  static int fft (const double *in, Complex *out, size_t npts, 
+		  size_t nsamples = 1, int stride = 1, int dist = -1);
+  static int fft (const Complex *in, Complex *out, size_t npts, 
+		  size_t nsamples = 1, int stride = 1, int dist = -1);
+  static int ifft (const Complex *in, Complex *out, size_t npts,
+		   size_t nsamples = 1, int stride = 1, int dist = -1);
 
-  static int fft2d (Complex*, size_t, size_t);
-  static int ifft2d (Complex*, size_t, size_t);
+  static int fftNd (const double*, Complex*, const int, const dim_vector &);
+  static int fftNd (const Complex*, Complex*, const int, 
+		    const dim_vector &);
+  static int ifftNd (const Complex*, Complex*, const int, 
+		     const dim_vector &);
 
 private:
   octave_fftw ();
--- a/src/ChangeLog	Mon Feb 16 17:57:34 2004 +0000
+++ b/src/ChangeLog	Mon Feb 16 19:02:33 2004 +0000
@@ -1,3 +1,20 @@
+2004-02-16 David Bateman <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/fft.cc: Adapt for Nd arrays, combine with ifft.cc.
+	* DLD-FUNCTIONS/ifft.cc: Delete.
+	* DLD-FUNCTIONS/fft2.cc: Adapt for Nd arrays, combine with ifft.cc.
+	* DLD-FUNCTIONS/ifft2.cc: Delete.
+	* DLD-FUNCTIONS/fftn.cc: New function for Nd FFT and inverse FFT.
+	* DLD-FUNCTIONS/fft_wisdom.cc: New function to manipulate FFTW 
+	wisdom.
+
+	* Makefile.in: Remove ifft.cc and ifft2.cc.  Add fftn.cc and 
+	ifftn.cc.  Make building of fft-wisdom.cc conditional on the 
+	value of FFTW_LIBS.
+	* defaults.cc (Vwisdom_prog): New variable
+	(set_default_wisdom_prog, wisdom_prog): New functions.
+	(symbols_of_defaults): Add DEFVAR for wisdom_prog.
+	
 2004-02-16  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* ov-list.cc (octave_list::subsasgn): Call assign for Cell objects
--- a/src/DLD-FUNCTIONS/fft.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/src/DLD-FUNCTIONS/fft.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -1,5 +1,6 @@
 /*
 
+Copyright (C) 2004 David Bateman
 Copyright (C) 1996, 1997 John W. Eaton
 
 This file is part of Octave.
@@ -32,96 +33,169 @@
 #include "oct-obj.h"
 #include "utils.h"
 
-// This function should be merged with Fifft.
+#if defined (HAVE_FFTW3)
+#define FFTSRC "@sc{Fftw}"
+#define WISDOM ", fft_wisdom"
+#else
+#define FFTSRC "@sc{Fftpack}"
+#define WISDOM 
+#endif
 
-DEFUN_DLD (fft, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} fft (@var{a}, @var{n})\n\
-Compute the FFT of @var{a} using subroutines from @sc{Fftpack}.  If @var{a}\n\
-is a matrix, @code{fft} computes the FFT for each column of @var{a}.\n\
-\n\
-If called with two arguments, @var{n} is expected to be an integer\n\
-specifying the number of elements of @var{a} to use.  If @var{a} is a\n\
-matrix, @var{n} specifies the number of rows of @var{a} to use.  If\n\
-@var{n} is larger than the size of @var{a}, @var{a} is resized and\n\
-padded with zeros.\n\
-@end deftypefn")
+static octave_value
+do_fft (const octave_value_list &args, const char *fcn, int type)
 {
   octave_value retval;
 
   int nargin = args.length ();
 
-  if (nargin < 1 || nargin > 2)
+  if (nargin < 1 || nargin > 3)
     {
-      print_usage ("fft");
+      print_usage (fcn);
       return retval;
     }
 
   octave_value arg = args(0);
-
-  int n_points = arg.rows ();
-  if (n_points == 1)
-    n_points = arg.columns ();
+  dim_vector dims = arg.dims ();
+  int n_points = -1;
+  int dim = -1;
+  
+  if (nargin > 1)
+    {
+      if (! args(1).is_empty ())
+	{
+	  double dval = args(1).double_value ();
+	  if (xisnan (dval))
+	    error ("%s: NaN is invalid as the N_POINTS", fcn);
+	  else
+	    {
+	      n_points = NINT (dval);
+	      if (n_points < 0)
+		error ("%s: number of points must be greater than zero", fcn);
+	    }
+	}
+    }
 
-  if (nargin == 2)
+  if (error_state)
+    return retval;
+
+  if (nargin > 2)
     {
-      double dval = args(1).double_value ();
+      double dval = args(2).double_value ();
       if (xisnan (dval))
-	error ("fft: NaN is invalid as the N_POINTS");
+	error ("%s: NaN is invalid as the N_POINTS", fcn);
+      else if (dval < 1 || dval > dims.length ())
+	error ("%s: invalid dimension along which to perform fft", fcn);
       else
-	n_points = NINT (dval);
+	dim = NINT (dval) - 1;
     }
 
   if (error_state)
     return retval;
 
-  if (n_points < 0)
-    {
-      error ("fft: number of points must be greater than zero");
+  for (int i = 0; i < dims.length (); i++)
+    if (dims(i) < 0)
       return retval;
-    }
 
-  int arg_is_empty = empty_arg ("fft", arg.rows (), arg.columns ());
+  if (dim < 0)
+    for (int i = 0; i < dims.length (); i++)
+      if ( dims(i) > 1)
+	{
+	  dim = i;
+	  break;
+	}
 
-  if (arg_is_empty < 0)
-    return retval;
-  else if (arg_is_empty || n_points == 0)
+  if (n_points < 0)
+    n_points = dims (dim);
+  else
+    dims (dim) = n_points;
+
+  if (dims.all_zero () || n_points == 0)
     return octave_value (Matrix ());
 
   if (arg.is_real_type ())
     {
-      Matrix m = arg.matrix_value ();
+      NDArray nda = arg.array_value ();
 
       if (! error_state)
 	{
-	  if (m.rows () == 1)
-	    m.resize (1, n_points, 0.0);
-	  else
-	    m.resize (n_points, m.columns (), 0.0);
-	  retval = m.fourier ();
+	  nda.resize (dims, 0.0);
+	  retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
 	}
     }
   else if (arg.is_complex_type ())
     {
-      ComplexMatrix m = arg.complex_matrix_value ();
+      ComplexNDArray cnda = arg.complex_array_value ();
 
       if (! error_state)
 	{
-	  if (m.rows () == 1)
-	    m.resize (1, n_points, 0.0);
-	  else
-	    m.resize (n_points, m.columns (), 0.0);
-	  retval = m.fourier ();
+	  cnda.resize (dims, 0.0);
+	  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
 	}
     }
   else
     {
-      gripe_wrong_type_arg ("fft", arg);
+      gripe_wrong_type_arg (fcn, arg);
     }
 
   return retval;
 }
 
+
+DEFUN_DLD (fft, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} fft (@var{a}, @var{n}, @var{dim})\n\
+Compute the FFT of @var{a} using subroutines from\n"
+FFTSRC
+". The FFT is calculated along the first non-singleton dimension of the\n\
+array. Thus if @var{a} is a matrix, @code{fft (@var{a})} computes the\n\
+FFT for each column of @var{a}.\n\
+\n\
+If called with two arguments, @var{n} is expected to be an integer\n\
+specifying the number of elements of @var{a} to use, or an empty\n\
+matrix to specify that its value should be ignored. If @var{n} is\n\
+larger than the dimension along which the FFT is calculated, then\n\
+@var{a} is resized and padded with zeros. Otherwise, if@var{n} is\n\
+smaller than the dimension along which the FFT is calculated, then\n\
+@var{a} is truncated.\n\
+\n\
+If called with three agruments, @var{dim} is an integer specifying the\n\
+dimension of the matrix along which the FFT is performed\n\
+@end deftypefn\n\
+@seealso {ifft, fft2, fftn"
+WISDOM
+"}")
+{
+  return do_fft(args, "fft", 0);
+}
+
+
+DEFUN_DLD (ifft, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} ifft (@var{a}, @var{n}, @var{dim})\n\
+Compute the inverse FFT of @var{a} using subroutines from\n"
+FFTSRC
+". The inverse FFT is calculated along the first non-singleton dimension\n\
+of the array. Thus if @var{a} is a matrix, @code{fft (@var{a})} computes\n\
+the inverse FFT for each column of @var{a}.\n\
+\n\
+If called with two arguments, @var{n} is expected to be an integer\n\
+specifying the number of elements of @var{a} to use, or an empty\n\
+matrix to specify that its value should be ignored. If @var{n} is\n\
+larger than the dimension along which the inverse FFT is calculated, then\n\
+@var{a} is resized and padded with zeros. Otherwise, if@var{n} is\n\
+smaller than the dimension along which the inverse FFT is calculated,\n\
+then @var{a} is truncated.\n\
+\n\
+If called with three agruments, @var{dim} is an integer specifying the\n\
+dimension of the matrix along which the inverse FFT is performed\n\
+@end deftypefn\n\
+@seealso {fft, ifft2, ifftn"
+WISDOM
+"}")
+{
+  return do_fft(args, "ifft", 1);
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/src/DLD-FUNCTIONS/fft2.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/src/DLD-FUNCTIONS/fft2.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -1,5 +1,6 @@
 /*
 
+Copyright (C) 2004 David Bateman
 Copyright (C) 1996, 1997 John W. Eaton
 
 This file is part of Octave.
@@ -32,18 +33,18 @@
 #include "oct-obj.h"
 #include "utils.h"
 
-// This function should be merged with Fifft2.
+// This function should be merged with Fifft.
 
-DEFUN_DLD (fft2, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} fft2 (@var{a}, @var{n}, @var{m})\n\
-Compute the two dimensional FFT of @var{a}.\n\
-\n\
-The optional arguments @var{n} and @var{m} may be used specify the\n\
-number of rows and columns of @var{a} to use.  If either of these is\n\
-larger than the size of @var{a}, @var{a} is resized and padded with\n\
-zeros.\n\
-@end deftypefn")
+#if defined (HAVE_FFTW3)
+#define FFTSRC "@sc{Fftw}"
+#define WISDOM ", fft_wisdom"
+#else
+#define FFTSRC "@sc{Fftpack}"
+#define WISDOM 
+#endif
+
+static octave_value
+do_fft2 (const octave_value_list &args, const char *fcn, int type)
 {
   octave_value retval;
 
@@ -51,79 +52,133 @@
 
   if (nargin < 1 || nargin > 3)
     {
-      print_usage ("fft2");
+      print_usage (fcn);
       return retval;
     }
 
   octave_value arg = args(0);
-
-  int n_rows = arg.rows ();
+  dim_vector dims = arg.dims ();
+  int n_rows = -1;
+  
   if (nargin > 1)
     {
       double dval = args(1).double_value ();
       if (xisnan (dval))
-	error ("fft2: NaN is invalid as N_ROWS");
+	error ("%s: NaN is invalid as the N_ROWS", fcn);
       else
-	n_rows = NINT (dval);
+	{
+	  n_rows = NINT (dval);
+	  if (n_rows < 0)
+	    error ("%s: number of rows must be greater than zero", fcn);
+	}
     }
 
   if (error_state)
     return retval;
 
-  int n_cols = arg.columns ();
+  int n_cols = -1;
   if (nargin > 2)
     {
       double dval = args(2).double_value ();
       if (xisnan (dval))
-	error ("fft2: NaN is invalid as N_COLS");
+	error ("%s: NaN is invalid as the N_COLS", fcn);
       else
-	n_cols = NINT (dval);
+	{
+	  n_cols = NINT (dval);
+	  if (n_cols < 0)
+	    error ("%s: number of columns must be greater than zero", fcn);
+	}
     }
 
   if (error_state)
     return retval;
 
-  if (n_rows < 0 || n_cols < 0)
-    {
-      error ("fft2: number of points must be greater than zero");
+  for (int i = 0; i < dims.length (); i++)
+    if (dims(i) < 0)
       return retval;
-    }
 
-  int arg_is_empty = empty_arg ("fft2", arg.rows (), arg.columns ());
+  if (n_rows < 0)
+    n_rows = dims (0);
+  else
+    dims (0) = n_rows;
 
-  if (arg_is_empty < 0)
-    return retval;
-  else if (arg_is_empty || n_rows == 0 || n_cols == 0)
+  if (n_cols < 0)
+    n_cols = dims (1);
+  else
+    dims (1) = n_cols;
+
+  if (dims.all_zero () || n_rows == 0 || n_cols == 0)
     return octave_value (Matrix ());
 
   if (arg.is_real_type ())
     {
-      Matrix m = arg.matrix_value ();
+      NDArray nda = arg.array_value ();
 
       if (! error_state)
 	{
-	  m.resize (n_rows, n_cols, 0.0);
-	  retval = m.fourier2d ();
+	  nda.resize (dims, 0.0);
+	  retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
 	}
     }
   else if (arg.is_complex_type ())
     {
-      ComplexMatrix m = arg.complex_matrix_value ();
+      ComplexNDArray cnda = arg.complex_array_value ();
 
       if (! error_state)
 	{
-	  m.resize (n_rows, n_cols, 0.0);
-	  retval = m.fourier2d ();
+	  cnda.resize (dims, 0.0);
+	  retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
 	}
     }
   else
     {
-      gripe_wrong_type_arg ("fft2", arg);
+      gripe_wrong_type_arg (fcn, arg);
     }
 
   return retval;
 }
 
+DEFUN_DLD (fft2, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} fft2 (@var{a}, @var{n}, @var{m})\n\
+Compute the two dimensional FFT of @var{a} using subroutines from\n"
+FFTSRC
+". The optional arguments @var{n} and @var{m} may be used specify the\n\
+number of rows and columns of @var{a} to use.  If either of these is\n\
+larger than the size of @var{a}, @var{a} is resized and padded with\n\
+zeros.\n\
+\n\
+If @var{a} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
+of @var{a} is treated seperately\n\
+@end deftypefn\n\
+@seealso {ifft2, fft, fftn"
+WISDOM
+"}")
+{
+  return do_fft2(args, "fft2", 0);
+}
+
+
+DEFUN_DLD (ifft2, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} fft2 (@var{a}, @var{n}, @var{m})\n\
+Compute the inverse two dimensional FFT of @var{a} using subroutines from\n"
+FFTSRC
+". The optional arguments @var{n} and @var{m} may be used specify the\n\
+number of rows and columns of @var{a} to use.  If either of these is\n\
+larger than the size of @var{a}, @var{a} is resized and padded with\n\
+zeros.\n\
+\n\
+If @var{a} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
+of @var{a} is treated seperately\n\
+@end deftypefn\n\
+@seealso {fft2, ifft, ifftn"
+WISDOM
+"}")
+{
+  return do_fft2(args, "ifft2", 1);
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/fft_wisdom.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -0,0 +1,190 @@
+/*
+
+Copyright (C) 2004 David Bateman
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <fftw3.h>
+
+#include "file-ops.h"
+#include "lo-sstream.h"
+#include "oct-env.h"
+#include "defaults.h"
+#include "lo-mappers.h"
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "utils.h"
+#include "sighandlers.h"
+
+// Name of the FFTW wisdom program we'd like to use.
+extern std::string Vwisdom_prog;
+
+DEFUN_DLD (fft_wisdom, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} fft_wisdom (@var{file}, @var{ow})\n\
+@deftypefnx {Loadable Function} {} fft_wisdom (@var{n})\n\
+\n\
+This function is used to manipulate the FFTW wisdom data. It can\n\
+be used to load previously stored wisdom from a file, create wisdom\n\
+needed by Octave and to save the current wisdom to a file. Wisdom\n\
+data can be used to significantly accelerate the calculation of the FFTs,\n\
+but is only profitable if the same FFT is called many times due to the\n\
+overhead in calculating the wisdom data.\n\
+\n\
+Called with a single string argument, @code{fft_wisdom (@var{file})}\n\
+will load the wisdom data found in @var{file}. If @var{file} does\n\
+not exist, then the current wisdom used by FFTW is saved to this\n\
+file. If the flag @var{ow} is non-zero, then even if @var{file}\n\
+exists, it will be overwritten.\n\
+\n\
+It is assumed that each row of @var{n} represents the size of a FFT for\n\
+which it is desired to pre-calculate the wisdom needed to accelerate it.\n\
+Any value of the matrix that is less than 1, is assumed to indicate an\n\
+absent dimension. That is\n\
+\n\
+@example\n\
+fftw_wisdom ([102, 0, 0; 103, 103, 0; 102, 103, 105]);\n\
+a = fft (rand (1,102));\n\
+b = fft (rand (103,103));\n\
+c = fftn (rand ([102, 103, 105]));\n\
+@end example\n\
+\n\
+calculates the wisdom necessary to accelerate the 103, 102x102 and\n\
+the 102x103x105 FFTs. Note that calculated wisdom will be lost when\n\
+when restarting Octave. However, if it is saved as discussed above, it\n\
+can be reloaded. Also, any system-wide wisdom file that has been found\n\
+will also be used. Saved wisdom files should not be used on different\n\
+platforms since they will not be efficient and the point of calculating\n\
+the wisdom is lost.\n\
+\n\
+Note that the program @code{fftw-wisdom} supplied with FFTW can equally\n\
+be used to create a file containing wisdom that can be imported into\n\
+Octave.\n\
+@end deftypefn\n\
+@seealso {fft, ifft, fft2, ifft2, fftn, ifftn}")
+{
+  octave_value retval;
+  int nargin = args.length();
+
+  if (nargin < 1 || nargin > 2)
+    {
+      print_usage ("fft_wisdom");
+      return retval;
+    }
+
+  if (args(0).is_string ())
+    {
+      bool overwrite = false;
+
+      if (nargin != 1)
+	{
+	  double dval = args (1).double_value ();
+	  if (NINT (dval) != 0)
+	    overwrite = true;
+	}
+
+      std::string wisdom = octave_env::make_absolute
+	(Vload_path_dir_path.find_first_of (args(0).string_value ().c_str ()),
+	 octave_env::getcwd ());
+
+      if (wisdom.empty () || overwrite)
+	{
+	  FILE *ofile = fopen (wisdom.c_str (), "wb");
+	  fftw_export_wisdom_to_file (ofile);
+	  fclose (ofile);
+	}
+      else
+	{
+	  FILE *ifile = fopen (wisdom.c_str (), "r");
+	  if (! fftw_import_wisdom_from_file (ifile))
+	    error ("fft_wisdom: can not import wisdom from file"); 
+	  fclose (ifile);
+	}
+
+    } 
+  else 
+    {
+      Matrix m = args (0).matrix_value ();
+
+      if (error_state)
+	{
+	  error ("fft_wisdom: argument must be a matrix or a string");
+	  return retval;
+	}
+
+      std::string name = file_ops::tempnam ("", "oct-");
+
+      if (name.empty ())
+	{
+	  error ("fft_wisdom: can not open temporary file");
+	  return retval;
+	}
+
+      OSSTREAM cmd_buf; 
+      cmd_buf << Vwisdom_prog << " -n -o \"" << name << "\"";
+
+      for (int k = 0; k < m.rows (); k++)
+	{
+	  bool first = true;
+
+	  cmd_buf << " ";
+
+	  // Note reversal of dimensions for column major storage in FFTW
+	  for (int j = m.columns()-1; j >= 0; j--)
+	    if (NINT(m(k,j)) > 0)
+	      {
+		if (first)
+		  first = false;
+		else
+		  cmd_buf << "x";
+		cmd_buf << NINT(m(k,j)) ;
+	      }
+	} 
+
+      cmd_buf << OSSTREAM_ENDS;
+
+      volatile octave_interrupt_handler old_interrupt_handler
+	= octave_ignore_interrupts ();
+
+      int status = system (OSSTREAM_C_STR (cmd_buf));
+
+      OSSTREAM_FREEZE (cmd_buf);
+
+      octave_set_interrupt_handler (old_interrupt_handler);
+
+      if (WIFEXITED (status))
+	{
+	  FILE *ifile = fopen (name.c_str (), "r");
+	  if (! fftw_import_wisdom_from_file (ifile))
+	    error ("fft_wisdom: can not import wisdom from temporary file"); 
+	  fclose (ifile);
+	}
+      else
+	error ("fft_wisdom: error running %s", Vwisdom_prog.c_str ());
+
+    }
+
+  return retval;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/fftn.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -0,0 +1,165 @@
+/*
+
+Copyright (C) 2004 David Bateman
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "lo-mappers.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "utils.h"
+
+// This function should be merged with Fifft.
+
+#if defined (HAVE_FFTW3)
+#define FFTSRC "@sc{Fftw}"
+#define WISDOM ", fft_wisdom"
+#else
+#define FFTSRC "@sc{Fftpack}"
+#define WISDOM 
+#endif
+
+static octave_value
+do_fftn (const octave_value_list &args, const char *fcn, int type)
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2)
+    {
+      print_usage (fcn);
+      return retval;
+    }
+
+  octave_value arg = args(0);
+  dim_vector dims = arg.dims ();
+  
+  for (int i = 0; i < dims.length (); i++)
+    if (dims(i) < 0)
+      return retval;
+
+  if (nargin > 1)
+    {
+      Matrix val = args(1).matrix_value ();
+      if (val.rows () > val.columns ())
+	val.transpose ();
+
+      if (val.columns () != dims.length () || val.rows () != 1)
+	error ("%s: second argument must be a vector of length dim", fcn);
+      else
+	{
+	  for (int i = 0; i < dims.length (); i++)
+	    {
+	      if (xisnan (val(i,0)))
+		error ("%s: NaN is invalid as a dimension", fcn);
+	      else if (NINT (val(i,0)) < 0)
+		error ("%s: all dimension must be greater than zero", fcn);
+	      else
+		{
+		  dims(i) = NINT(val(i,0));
+		}
+	    }
+	}
+    }
+
+  if (error_state)
+    return retval;
+
+  if (dims.all_zero ())
+    return octave_value (Matrix ());
+
+  if (arg.is_real_type ())
+    {
+      NDArray nda = arg.array_value ();
+
+      if (! error_state)
+	{
+	  nda.resize (dims, 0.0);
+	  retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
+	}
+    }
+  else if (arg.is_complex_type ())
+    {
+      ComplexNDArray cnda = arg.complex_array_value ();
+
+      if (! error_state)
+	{
+	  cnda.resize (dims, 0.0);
+	  retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
+	}
+    }
+  else
+    {
+      gripe_wrong_type_arg (fcn, arg);
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (fftn, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} fftn (@var{a}, @var{siz})\n\
+Compute the N dimensional FFT of @var{a} using subroutines from\n"
+FFTSRC
+". The optional vector argument @var{siz} may be used specify the\n\
+dimensions of the array to be used. If an element of @var{siz} is\n\
+smaller than the corresponding dimension, then the dimension is\n\
+truncated prior to performing the FFT. Otherwise if an element\n\
+of @var{siz} is larger than the corresponding dimension @var{a}\n\
+is resized and padded with zeros.\n\
+@end deftypefn\n\
+@seealso {ifftn, fft, fft2"
+WISDOM
+"}")
+{
+  return do_fftn(args, "fftn", 0);
+}
+
+DEFUN_DLD (ifftn, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} ifftn (@var{a}, @var{siz})\n\
+Compute the invesre N dimensional FFT of @var{a} using subroutines from\n"
+FFTSRC
+". The optional vector argument @var{siz} may be used specify the\n\
+dimensions of the array to be used. If an element of @var{siz} is\n\
+smaller than the corresponding dimension, then the dimension is\n\
+truncated prior to performing the inverse FFT. Otherwise if an element\n\
+of @var{siz} is larger than the corresponding dimension @var{a}\n\
+is resized and padded with zeros.\n\
+@end deftypefn\n\
+@seealso {fftn, ifft, ifft2"
+WISDOM
+"}")
+{
+  return do_fftn(args, "ifftn", 1);
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- a/src/DLD-FUNCTIONS/ifft.cc	Mon Feb 16 17:57:34 2004 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,130 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 2, or (at your option) any
-later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "lo-mappers.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-
-// This function should be merged with Ffft.
-
-DEFUN_DLD (ifft, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} ifft (@var{a}, @var{n})\n\
-Compute the inverse FFT of @var{a} using subroutines from @sc{Fftpack}.  If\n\
-@var{a} is a matrix, @code{fft} computes the inverse FFT for each column\n\
-of @var{a}.\n\
-\n\
-If called with two arguments, @var{n} is expected to be an integer\n\
-specifying the number of elements of @var{a} to use.  If @var{a} is a\n\
-matrix, @var{n} specifies the number of rows of @var{a} to use.  If\n\
-@var{n} is larger than the size of @var{a}, @var{a} is resized and\n\
-padded with zeros.\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  if (nargin < 1 || nargin > 2)
-    {
-      print_usage ("ifft");
-      return retval;
-    }
-
-  octave_value arg = args(0);
-    
-  int n_points = arg.rows ();
-  if (n_points == 1)
-    n_points = arg.columns ();
-
-  if (nargin == 2)
-    {
-      double dval = args(1).double_value ();
-      if (xisnan (dval))
-	error ("fft: NaN is invalid as the N_POINTS");
-      else
-	n_points = NINT (dval);
-    }
-
-  if (error_state)
-    return retval;
-
-  if (n_points < 0)
-    {
-      error ("ifft: number of points must be greater than zero");
-      return retval;
-    }
-
-  int arg_is_empty = empty_arg ("ifft", arg.rows (), arg.columns ());
-
-  if (arg_is_empty < 0)
-    return retval;
-  else if (arg_is_empty || n_points == 0)
-    return octave_value (Matrix ());
-
-  if (arg.is_real_type ())
-    {
-      Matrix m = arg.matrix_value ();
-
-      if (! error_state)
-	{
-	  if (m.rows () == 1)
-	    m.resize (1, n_points, 0.0);
-	  else
-	    m.resize (n_points, m.columns (), 0.0);
-	  retval = m.ifourier ();
-	}
-    }
-  else if (arg.is_complex_type ())
-    {
-      ComplexMatrix m = arg.complex_matrix_value ();
-
-      if (! error_state)
-	{
-	  if (m.rows () == 1)
-	    m.resize (1, n_points, 0.0);
-	  else
-	    m.resize (n_points, m.columns (), 0.0);
-	  retval = m.ifourier ();
-	}
-    }
-  else
-    {
-      gripe_wrong_type_arg ("ifft", arg);
-    }
-
-  return retval;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/src/DLD-FUNCTIONS/ifft2.cc	Mon Feb 16 17:57:34 2004 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,131 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 2, or (at your option) any
-later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "lo-mappers.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "utils.h"
-
-// This function should be merged with Ffft2.
-
-DEFUN_DLD (ifft2, args, ,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {} ifft2 (@var{a}, @var{n}, @var{m})\n\
-Compute the two dimensional inverse FFT of @var{a}.\n\
-\n\
-The optional arguments @var{n} and @var{m} may be used specify the\n\
-number of rows and columns of @var{a} to use.  If either of these is\n\
-larger than the size of @var{a}, @var{a} is resized and padded with\n\
-zeros.\n\
-@end deftypefn")
-{
-  octave_value retval;
-
-  int nargin = args.length ();
-
-  if (nargin < 1 || nargin > 3)
-    {
-      print_usage ("ifft2");
-      return retval;
-    }
-
-  octave_value arg = args(0);
-
-  int n_rows = arg.rows ();
-  if (nargin > 1)
-    {
-      double dval = args(1).double_value ();
-      if (xisnan (dval))
-	error ("fft2: NaN is invalid as N_ROWS");
-      else
-	n_rows = NINT (dval);
-    }
-
-  if (error_state)
-    return retval;
-
-  int n_cols = arg.columns ();
-  if (nargin > 2)
-    {
-      double dval = args(2).double_value ();
-      if (xisnan (dval))
-	error ("fft2: NaN is invalid as N_COLS");
-      else
-	n_cols = NINT (dval);
-    }
-
-  if (error_state)
-    return retval;
-
-  if (n_rows < 0 || n_cols < 0)
-    {
-      error ("ifft2: number of points must be greater than zero");
-      return retval;
-    }
-
-  int arg_is_empty = empty_arg ("ifft2", arg.rows (), arg.columns ());
-
-  if (arg_is_empty < 0)
-    return retval;
-  else if (arg_is_empty || n_rows == 0 || n_cols == 0)
-    return octave_value (Matrix ());
-
-  if (arg.is_real_type ())
-    {
-      Matrix m = arg.matrix_value ();
-
-      if (! error_state)
-	{
-	  m.resize (n_rows, n_cols, 0.0);
-	  retval = m.ifourier2d ();
-	}
-    }
-  else if (arg.is_complex_type ())
-    {
-      ComplexMatrix m = arg.complex_matrix_value ();
-
-      if (! error_state)
-	{
-	  m.resize (n_rows, n_cols, 0.0);
-	  retval = m.ifourier2d ();
-	}
-    }
-  else
-    {
-      gripe_wrong_type_arg ("ifft2", arg);
-    }
-
-  return retval;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/src/Makefile.in	Mon Feb 16 17:57:34 2004 +0000
+++ b/src/Makefile.in	Mon Feb 16 19:02:33 2004 +0000
@@ -43,10 +43,14 @@
 OPT_HANDLERS := DASPK-opts.cc DASRT-opts.cc DASSL-opts.cc \
 	LSODE-opts.cc NLEqn-opts.cc ODESSA-opts.cc Quad-opts.cc
 
+ifneq ($(FFTW_LIBS), )
+  DLD_WISDOM := fft_wisdom.cc
+endif
+
 DLD_XSRC := balance.cc besselj.cc betainc.cc chol.cc colloc.cc \
 	daspk.cc dasrt.cc dassl.cc det.cc eig.cc expm.cc fft.cc fft2.cc \
-	filter.cc find.cc fsolve.cc gammainc.cc getgrent.cc \
-	getpwent.cc getrusage.cc givens.cc hess.cc ifft.cc ifft2.cc \
+	fftn.cc $(DLD_WISDOM) filter.cc find.cc fsolve.cc gammainc.cc \
+	getgrent.cc getpwent.cc getrusage.cc givens.cc hess.cc \
 	inv.cc kron.cc lpsolve.cc lsode.cc lu.cc minmax.cc \
 	odessa.cc pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc \
 	sort.cc sqrtm.cc svd.cc syl.cc time.cc
--- a/src/defaults.cc	Mon Feb 16 17:57:34 2004 +0000
+++ b/src/defaults.cc	Mon Feb 16 19:02:33 2004 +0000
@@ -92,6 +92,11 @@
 std::string Vlocal_site_defaults_file;
 std::string Vsite_defaults_file;
 
+#ifdef HAVE_FFTW3
+// Name of the FFTW wisdom program
+std::string Vwisdom_prog;
+#endif
+
 // Each element of A and B should be directory names.  For each
 // element of A not in the list B, execute SCRIPT_FILE in that
 // directory if it exists.
@@ -310,6 +315,19 @@
     Vinfo_prog = std::string (oct_info_prog);
 }
 
+#ifdef HAVE_FFTW3
+static void
+set_default_wisdom_prog (void)
+{
+  std::string oct_wisdom_prog = octave_env::getenv ("OCTAVE_WISDOM_PROGRAM");
+
+  if (oct_wisdom_prog.empty ())
+    Vwisdom_prog = "fftw-wisdom";
+  else
+    Vwisdom_prog = std::string (oct_wisdom_prog);
+}
+#endif
+
 static void
 set_default_editor (void)
 {
@@ -400,6 +418,10 @@
 
   set_default_info_prog ();
 
+#ifdef HAVE_FFTW3
+  set_default_wisdom_prog ();
+#endif
+
   set_default_editor ();
 
   set_local_site_defaults_file ();
@@ -484,6 +506,26 @@
   return status;
 }
 
+#ifdef HAVE_FFTW3
+static int
+wisdom_prog (void)
+{
+  int status = 0;
+
+  std::string s = builtin_string_variable ("WISDOM_PROGRAM");
+
+  if (s.empty ())
+    {
+      gripe_invalid_value_specified ("WISDOM_PROGRAM");
+      status = -1;
+    }
+  else
+    Vwisdom_prog = s;
+
+  return status;
+}
+#endif
+
 static int
 default_exec_path (void)
 {
@@ -587,6 +629,18 @@
 @code{\"emacs\"}.\n\
 @end defvr");
 
+#ifdef HAVE_FFTW3
+  DEFVAR (WISDOM_PROGRAM, Vwisdom_prog, wisdom_prog,
+    "-*- texinfo -*-\n\
+@defvr {Built-in Variable} WISDOM_PROGRAM\n\
+A string naming the FFTW wisdom program to use to create wisdom data\n\
+to accelerate Fourier transforms. If the environment variable\n\
+@code{OCTAVE_WISDOM_PROGRAM} is set when Octave starts, its value is used\n\
+as the default. Otherwise, @code{WISDOM_PROGRAM} is set to\n\
+@code{\"fftw-wisdom\"}.\n\
+@end defvr");
+#endif
+  
   DEFVAR (EXEC_PATH, Vexec_path, exec_path,
     "-*- texinfo -*-\n\
 @defvr {Built-in Variable} EXEC_PATH\n\
@@ -621,7 +675,7 @@
 substituted for leading, trailing, or doubled colons that appear in the\n\
 built-in variable @code{EXEC_PATH}.\n\
 @end defvr");
-  
+
   DEFVAR (LOADPATH, Vload_path, loadpath,
     "-*- texinfo -*-\n\
 @defvr {Built-in Variable} LOADPATH\n\