*** ./doc/interpreter/signal.txi.orig	2004-02-12 15:57:00.000000000 +0100
--- ./doc/interpreter/signal.txi	2004-02-12 15:57:04.000000000 +0100
***************
*** 19,24 ****
--- 19,28 ----
  
  @DOCSTRING(ifft2)
  
+ @DOCSTRING(fftn)
+ 
+ @DOCSTRING(ifftn)
+ 
  @DOCSTRING(fftconv)
  
  @DOCSTRING(fftfilt)
*** ./liboctave/ChangeLog.orig	2004-01-30 15:24:31.000000000 +0100
--- ./liboctave/ChangeLog	2004-02-16 17:07:00.000000000 +0100
***************
*** 1,3 ****
--- 1,21 ----
+ 2004-02-16  David Bateman <address@hidden>
+ 
+ 	* oct-fftw.cc: Add support for FFTW 3.x. Include the ability to
+ 	use the real to complex transform for fft's of real matrices
+ 	* oct-fftw.h: Decls updated for the above
+ 
+ 	* dMatrix.cc: 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: 1D fft of a matrix done as single call rather than 
+ 	loop. Update for FFTW 3.x
+ 
+ 	* dNDArray.cc (fourier, ifourier, fourierNd, ifouriourNd: New
+ 	fourier transform functions for Nd arrays
+ 	* dNArray.h (fourier, ifourier, fourierNd, ifouriourNd: decls
+ 	* CNDArray.cc (fourier, ifourier, fourierNd, ifouriourNd: New
+ 	fourier transform functions for complex Nd arrays
+ 	* CNArray.h (fourier, ifourier, fourierNd, ifouriourNd: decls
+ 	
  2004-01-22  John W. Eaton  <address@hidden>
  
  	* Array.cc (Array<T>::assign2, Array<T>::assignN):
*** ./liboctave/dMatrix.cc.orig	2004-01-28 15:17:03.000000000 +0100
--- ./liboctave/dMatrix.cc	2004-02-16 17:50:45.000000000 +0100
***************
*** 52,58 ****
  #include "oct-cmplx.h"
  #include "quit.h"
  
! #ifdef HAVE_FFTW
  #include "oct-fftw.h"
  #endif
  
--- 52,58 ----
  #include "oct-cmplx.h"
  #include "quit.h"
  
! #if defined (HAVE_FFTW3)
  #include "oct-fftw.h"
  #endif
  
***************
*** 766,772 ****
      }
  }
  
! #ifdef HAVE_FFTW
  
  ComplexMatrix
  Matrix::fourier (void) const
--- 766,772 ----
      }
  }
  
! #if defined (HAVE_FFTW3)
  
  ComplexMatrix
  Matrix::fourier (void) const
***************
*** 789,804 ****
        nsamples = nc;
      }
  
!   ComplexMatrix tmp (*this);
!   Complex *in (tmp.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);
!     }
  
    return retval;
  }
--- 789,798 ----
        nsamples = nc;
      }
  
!   const double *in (fortran_vec ());
    Complex *out (retval.fortran_vec ());
  
!   octave_fftw::fft (in, out, npts, nsamples); 
  
    return retval;
  }
***************
*** 828,839 ****
    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);
!     }
  
    return retval;
  }
--- 822,828 ----
    Complex *in (tmp.fortran_vec ());
    Complex *out (retval.fortran_vec ());
  
!   octave_fftw::ifft (in, out, npts, nsamples); 
  
    return retval;
  }
***************
*** 841,853 ****
  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);
  
    return retval;
  }
--- 830,840 ----
  ComplexMatrix
  Matrix::fourier2d (void) const
  {
!   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,867 ****
  ComplexMatrix
  Matrix::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);
  
    return retval;
  }
--- 842,853 ----
  ComplexMatrix
  Matrix::ifourier2d (void) const
  {
!   dim_vector dv(rows (), cols ());
  
    ComplexMatrix retval (*this);
!   Complex *out (retval.fortran_vec ());
! 
!   octave_fftw::ifftNd (out, out, 2, dv);
  
    return retval;
  }
***************
*** 998,1005 ****
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> row (npts);
!   Complex *prow = row.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
--- 984,991 ----
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> tmp (npts);
!   Complex *prow = tmp.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
***************
*** 1067,1074 ****
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> row (npts);
!   Complex *prow = row.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
--- 1053,1060 ----
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> tmp (npts);
!   Complex *prow = tmp.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
*** ./liboctave/CMatrix.cc.orig	2004-01-28 15:17:09.000000000 +0100
--- ./liboctave/CMatrix.cc	2004-02-16 17:50:36.000000000 +0100
***************
*** 56,62 ****
  #include "mx-inlines.cc"
  #include "oct-cmplx.h"
  
! #ifdef HAVE_FFTW
  #include "oct-fftw.h"
  #endif
  
--- 56,62 ----
  #include "mx-inlines.cc"
  #include "oct-cmplx.h"
  
! #if defined (HAVE_FFTW3)
  #include "oct-fftw.h"
  #endif
  
***************
*** 1102,1108 ****
    return retval;
  }
  
! #ifdef HAVE_FFTW
  
  ComplexMatrix
  ComplexMatrix::fourier (void) const
--- 1102,1108 ----
    return retval;
  }
  
! #if defined (HAVE_FFTW3)
  
  ComplexMatrix
  ComplexMatrix::fourier (void) const
***************
*** 1128,1139 ****
    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);
!     }
  
    return retval;
  }
--- 1128,1134 ----
    const Complex *in (data ());
    Complex *out (retval.fortran_vec ());
  
!   octave_fftw::fft (in, out, npts, nsamples); 
  
    return retval;
  }
***************
*** 1162,1173 ****
    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);
!     }
  
    return retval;
  }
--- 1157,1163 ----
    const Complex *in (data ());
    Complex *out (retval.fortran_vec ());
  
!   octave_fftw::ifft (in, out, npts, nsamples); 
  
    return retval;
  }
***************
*** 1175,1187 ****
  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);
  
    return retval;
  }
--- 1165,1177 ----
  ComplexMatrix
  ComplexMatrix::fourier2d (void) const
  {
!   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,1201 ****
  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);
  
    return retval;
  }
--- 1179,1191 ----
  ComplexMatrix
  ComplexMatrix::ifourier2d (void) const
  {
!   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,1339 ****
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> row (npts);
!   Complex *prow = row.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
--- 1322,1329 ----
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> tmp (npts);
!   Complex *prow = tmp.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
***************
*** 1401,1408 ****
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> row (npts);
!   Complex *prow = row.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
--- 1391,1398 ----
    wsave.resize (nn);
    pwsave = wsave.fortran_vec ();
  
!   Array<Complex> tmp (npts);
!   Complex *prow = tmp.fortran_vec ();
  
    F77_FUNC (cffti, CFFTI) (npts, pwsave);
  
*** ./liboctave/oct-fftw.h.orig	2004-01-28 15:17:17.000000000 +0100
--- ./liboctave/oct-fftw.h	2004-02-10 11:39:44.000000000 +0100
***************
*** 22,45 ****
  #define octave_oct_fftw_h 1
  
  #include <cstddef>
! 
! #if defined (HAVE_DFFTW_H)
! #include <dfftw.h>
! #else
! #include <fftw.h>
! #endif
  
  #include "oct-cmplx.h"
  
  class
  octave_fftw
  {
  public:
!   static int fft (const Complex*, Complex *, size_t);
!   static int ifft (const Complex*, Complex *, size_t);
! 
!   static int fft2d (Complex*, size_t, size_t);
!   static int ifft2d (Complex*, size_t, size_t);
  
  private:
    octave_fftw ();
--- 22,48 ----
  #define octave_oct_fftw_h 1
  
  #include <cstddef>
! #include <fftw3.h>
  
  #include "oct-cmplx.h"
+ #include "dim-vector.h"
  
  class
  octave_fftw
  {
  public:
!   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 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 ();
*** ./liboctave/oct-fftw.cc.orig	2004-01-28 15:17:22.000000000 +0100
--- ./liboctave/oct-fftw.cc	2004-02-16 16:45:56.000000000 +0100
***************
*** 22,39 ****
  #include <config.h>
  #endif
  
! #ifdef HAVE_FFTW
  
  #include "oct-fftw.h"
  #include "lo-error.h"
! 
  
  // 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.
  
  class
  octave_fftw_planner
--- 22,47 ----
  #include <config.h>
  #endif
  
! #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.
! 
! // 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,57 ****
  public:
    octave_fftw_planner ();
  
!   fftw_plan create_plan (fftw_direction, size_t);
!   fftwnd_plan create_plan2d (fftw_direction, size_t, size_t);
  
  private:
    int plan_flags;
  
    fftw_plan plan[2];
!   fftwnd_plan plan2d[2];
! 
!   size_t n[2];
!   size_t n2d[2][2];
  };
  
  octave_fftw_planner::octave_fftw_planner ()
--- 49,83 ----
  public:
    octave_fftw_planner ();
  
!   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];
!   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];
! 
!   // 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,89 ****
    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;
  }
  
  fftw_plan
! octave_fftw_planner::create_plan (fftw_direction dir, size_t npts)
  {
!   size_t which = (dir == FFTW_FORWARD) ? 0 : 1;
    fftw_plan *cur_plan_p = &plan[which];
    bool create_new_plan = false;
  
!   if (plan[which] == 0 || n[which] != npts)
!     {
!       create_new_plan = true;
!       n[which] = npts;
!     }
  
    if (create_new_plan)
      {
        if (*cur_plan_p)
  	fftw_destroy_plan (*cur_plan_p);
  
!       *cur_plan_p = fftw_create_plan (npts, dir, plan_flags);
  
        if (*cur_plan_p == 0)
  	(*current_liboctave_error_handler) ("Error creating fftw plan");
--- 85,151 ----
    plan_flags = FFTW_ESTIMATE;
  
    plan[0] = plan[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 (int dir, const int rank,
! 				  const dim_vector dims, int howmany,
! 				  int stride, int dist, 
! 				  const Complex *in, Complex *out)
  {
!   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 || 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);
  
!       // 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,123 ****
    return *cur_plan_p;
  }
   
! fftwnd_plan
! octave_fftw_planner::create_plan2d (fftw_direction dir, 
!                                    size_t nrows, size_t ncols)
  {
!   size_t which = (dir == FFTW_FORWARD) ? 0 : 1;
!   fftwnd_plan *cur_plan_p = &plan2d[which];
    bool create_new_plan = false;
  
!   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 (create_new_plan)
      {
        if (*cur_plan_p)
! 	fftwnd_destroy_plan (*cur_plan_p);
  
!       *cur_plan_p = fftw2d_create_plan (nrows, ncols, dir, 
! 					plan_flags | FFTW_IN_PLACE);
  
        if (*cur_plan_p == 0)
! 	(*current_liboctave_error_handler) ("Error creating 2d fftw plan");
      }
  
    return *cur_plan_p;
--- 154,208 ----
    return *cur_plan_p;
  }
   
! 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)
  {
!   fftw_plan *cur_plan_p = &rplan;
    bool create_new_plan = false;
+   char in_align = ((int) in) & 0xF;
+   char out_align = ((int) out) & 0xF;
  
!   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)
! 	fftw_destroy_plan (*cur_plan_p);
  
!       // 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 fftw plan");
      }
  
    return *cur_plan_p;
***************
*** 125,175 ****
  
  static octave_fftw_planner fftw_planner;
  
  int
! octave_fftw::fft (const Complex *in, Complex *out, size_t npts)
  {
!   fftw_one (fftw_planner.create_plan (FFTW_FORWARD, npts),
!             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)
  {
!   fftw_one (fftw_planner.create_plan (FFTW_BACKWARD, npts),
!             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;
  
    return 0;
  }
  
  int
! octave_fftw::fft2d (Complex *inout, size_t nr, size_t nc)
  {
!   fftwnd_one (fftw_planner.create_plan2d (FFTW_FORWARD, nr, nc),
!               reinterpret_cast<fftw_complex *> (inout),
!               0);
  
    return 0;
  }
  
  int
! octave_fftw::ifft2d (Complex *inout, size_t nr, size_t nc)
  {
!   fftwnd_one (fftw_planner.create_plan2d (FFTW_BACKWARD, nr, nc),
!               reinterpret_cast<fftw_complex *> (inout),
!               0);
  
!   const size_t npts = nr * nc;
    const Complex scale = npts;
    for (size_t i = 0; i < npts; i++)
!     inout[i] /= scale;
  
    return 0;
  }
--- 210,393 ----
  
  static octave_fftw_planner fftw_planner;
  
+ 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)
+ {
+   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::fft (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_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 j = 0; j < nsamples; j++)
!     for (size_t i = 0; i < npts; i++)
!       out[i*stride + j*dist] /= scale;
! 
!   return 0;
! }
! 
! int
! octave_fftw::fftNd (const double *in, Complex *out, const int rank, 
! 		    const dim_vector &dv)
! {
!   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::fftNd (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_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;
  }
  
  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++)
!     out[i] /= scale;
  
    return 0;
  }
*** ./liboctave/CNDArray.cc.orig	2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/CNDArray.cc	2004-02-16 17:50:41.000000000 +0100
***************
*** 34,42 ****
--- 34,64 ----
  #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,66 ****
--- 83,511 ----
      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
*** ./liboctave/CNDArray.h.orig	2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/CNDArray.h	2004-02-10 17:24:25.000000000 +0100
***************
*** 88,93 ****
--- 88,102 ----
    ComplexNDArray sum (int dim = -1) const;
    ComplexNDArray sumsq (int dim = -1) 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;
+ 
    NDArray abs (void) const;
  
    ComplexMatrix matrix_value (void) const;
*** ./liboctave/dNDArray.cc.orig	2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/dNDArray.cc	2004-02-16 17:50:53.000000000 +0100
***************
*** 34,43 ****
--- 34,486 ----
  #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 ())
  {
*** ./liboctave/dNDArray.h.orig	2004-01-21 05:46:58.000000000 +0100
--- ./liboctave/dNDArray.h	2004-02-10 17:24:16.000000000 +0100
***************
*** 87,92 ****
--- 87,101 ----
    NDArray sum (int dim = -1) const;  
    NDArray sumsq (int dim = -1) 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;
+ 
    NDArray abs (void) const;
  
    friend NDArray real (const ComplexNDArray& a);
*** ./src/DLD-FUNCTIONS/ifft.cc.orig	2004-01-30 22:46:35.000000000 +0100
--- ./src/DLD-FUNCTIONS/ifft.cc	2004-02-11 14:08:36.000000000 +0100
***************
*** 1,130 ****
- /*
- 
- 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: ***
- */
--- 0 ----
*** ./src/DLD-FUNCTIONS/fftn.cc.orig	2004-02-11 14:09:07.000000000 +0100
--- ./src/DLD-FUNCTIONS/fftn.cc	2004-02-13 17:46:09.000000000 +0100
***************
*** 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: ***
+ */
*** ./src/DLD-FUNCTIONS/fft.cc.orig	2004-01-30 22:46:06.000000000 +0100
--- ./src/DLD-FUNCTIONS/fft.cc	2004-02-13 17:43:10.000000000 +0100
***************
*** 1,5 ****
--- 1,6 ----
  /*
  
+ Copyright (C) 2004 David Bateman
  Copyright (C) 1996, 1997 John W. Eaton
  
  This file is part of Octave.
***************
*** 32,127 ****
  #include "oct-obj.h"
  #include "utils.h"
  
! // This function should be merged with Fifft.
  
! 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")
  {
    octave_value retval;
  
    int nargin = args.length ();
  
!   if (nargin < 1 || nargin > 2)
      {
!       print_usage ("fft");
        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 ("fft: number of points must be greater than zero");
        return retval;
-     }
  
!   int arg_is_empty = empty_arg ("fft", 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.fourier ();
  	}
      }
    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.fourier ();
  	}
      }
    else
      {
!       gripe_wrong_type_arg ("fft", arg);
      }
  
    return retval;
  }
  
  /*
  ;;; Local Variables: ***
  ;;; mode: C++ ***
--- 33,201 ----
  #include "oct-obj.h"
  #include "utils.h"
  
! #if defined (HAVE_FFTW3)
! #define FFTSRC "@sc{Fftw}"
! #define WISDOM ", fft_wisdom"
! #else
! #define FFTSRC "@sc{Fftpack}"
! #define WISDOM 
! #endif
  
! 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 > 3)
      {
!       print_usage (fcn);
        return retval;
      }
  
    octave_value arg = args(0);
+   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 (error_state)
!     return retval;
  
!   if (nargin > 2)
      {
!       double dval = args(2).double_value ();
        if (xisnan (dval))
! 	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
! 	dim = NINT (dval) - 1;
      }
  
    if (error_state)
      return retval;
  
!   for (int i = 0; i < dims.length (); i++)
!     if (dims(i) < 0)
        return retval;
  
!   if (dim < 0)
!     for (int i = 0; i < dims.length (); i++)
!       if ( dims(i) > 1)
! 	{
! 	  dim = i;
! 	  break;
! 	}
  
!   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 ())
      {
!       NDArray nda = arg.array_value ();
  
        if (! error_state)
  	{
! 	  nda.resize (dims, 0.0);
! 	  retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
  	}
      }
    else if (arg.is_complex_type ())
      {
!       ComplexNDArray cnda = arg.complex_array_value ();
  
        if (! error_state)
  	{
! 	  cnda.resize (dims, 0.0);
! 	  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
  	}
      }
    else
      {
!       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, address@hidden 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, address@hidden 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++ ***
*** ./src/DLD-FUNCTIONS/fft2.cc.orig	2004-01-30 22:49:01.000000000 +0100
--- ./src/DLD-FUNCTIONS/fft2.cc	2004-02-13 17:43:22.000000000 +0100
***************
*** 1,5 ****
--- 1,6 ----
  /*
  
+ Copyright (C) 2004 David Bateman
  Copyright (C) 1996, 1997 John W. Eaton
  
  This file is part of Octave.
***************
*** 32,49 ****
  #include "oct-obj.h"
  #include "utils.h"
  
! // This function should be merged with Fifft2.
  
! 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")
  {
    octave_value retval;
  
--- 33,50 ----
  #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_fft2 (const octave_value_list &args, const char *fcn, int type)
  {
    octave_value retval;
  
***************
*** 51,129 ****
  
    if (nargin < 1 || nargin > 3)
      {
!       print_usage ("fft2");
        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 ("fft2: number of points must be greater than zero");
        return retval;
-     }
  
!   int arg_is_empty = empty_arg ("fft2", 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.fourier2d ();
  	}
      }
    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.fourier2d ();
  	}
      }
    else
      {
!       gripe_wrong_type_arg ("fft2", arg);
      }
  
    return retval;
  }
  
  /*
  ;;; Local Variables: ***
  ;;; mode: C++ ***
--- 52,184 ----
  
    if (nargin < 1 || nargin > 3)
      {
!       print_usage (fcn);
        return retval;
      }
  
    octave_value arg = args(0);
!   dim_vector dims = arg.dims ();
!   int n_rows = -1;
!   
    if (nargin > 1)
      {
        double dval = args(1).double_value ();
        if (xisnan (dval))
! 	error ("%s: NaN is invalid as the N_ROWS", fcn);
        else
! 	{
! 	  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 = -1;
    if (nargin > 2)
      {
        double dval = args(2).double_value ();
        if (xisnan (dval))
! 	error ("%s: NaN is invalid as the N_COLS", fcn);
        else
! 	{
! 	  n_cols = NINT (dval);
! 	  if (n_cols < 0)
! 	    error ("%s: number of columns must be greater than zero", fcn);
! 	}
      }
  
    if (error_state)
      return retval;
  
!   for (int i = 0; i < dims.length (); i++)
!     if (dims(i) < 0)
        return retval;
  
!   if (n_rows < 0)
!     n_rows = dims (0);
!   else
!     dims (0) = n_rows;
  
!   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 ())
      {
!       NDArray nda = arg.array_value ();
  
        if (! error_state)
  	{
! 	  nda.resize (dims, 0.0);
! 	  retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
  	}
      }
    else if (arg.is_complex_type ())
      {
!       ComplexNDArray cnda = arg.complex_array_value ();
  
        if (! error_state)
  	{
! 	  cnda.resize (dims, 0.0);
! 	  retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
  	}
      }
    else
      {
!       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++ ***
*** ./src/DLD-FUNCTIONS/ifft2.cc.orig	2004-01-30 22:50:46.000000000 +0100
--- ./src/DLD-FUNCTIONS/ifft2.cc	2004-02-11 14:08:32.000000000 +0100
***************
*** 1,131 ****
- /*
- 
- 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: ***
- */
--- 0 ----
*** ./src/DLD-FUNCTIONS/fft_wisdom.cc.orig	2004-02-11 14:09:27.000000000 +0100
--- ./src/DLD-FUNCTIONS/fft_wisdom.cc	2004-02-16 17:56:13.000000000 +0100
***************
*** 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;
+ }
*** ./src/Makefile.in.orig	2004-02-11 14:08:02.000000000 +0100
--- ./src/Makefile.in	2004-02-13 17:46:26.000000000 +0100
***************
*** 43,52 ****
  OPT_HANDLERS := DASPK-opts.cc DASRT-opts.cc DASSL-opts.cc \
  	LSODE-opts.cc NLEqn-opts.cc ODESSA-opts.cc Quad-opts.cc
  
  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 \
  	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
--- 43,56 ----
  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 \
! 	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
*** ./src/ChangeLog.orig	2004-01-30 23:14:31.000000000 +0100
--- ./src/ChangeLog	2004-02-16 17:07:09.000000000 +0100
***************
*** 1,3 ****
--- 1,21 ----
+ 2004-02-16 David Bateman <address@hidden>
+ 
+ 	* DLD-FUNCTIONS/fft.cc: Adapt for Nd arrays, combine with ifft.cc.
+ 	* DLD-FUNCTIONS/ifft.cc: remove and link to fft.cc.
+ 	* DLD-FUNCTIONS/fft2.cc: Adapt for Nd arrays, combine with ifft.cc.
+ 	* DLD-FUNCTIONS/ifft2.cc: remove and link to fft.cc.
+ 	* DLD-FUNCTIONS/fftn.cc: new function for Nd FFT, combined with 
+ 	ifftn.cc.
+ 	* DLD-FUNCTIONS/ifftn.cc: link to fftn.cc.
+ 	* 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: Add Vwisdom_prog variable and the corresponding
+ 	octave and environment variables.
+ 	
  2004-01-22  John W. Eaton  <address@hidden>
  
  	* error.cc (pr_where): New arg, print_code with default value true.
*** ./src/defaults.cc.orig	2003-08-22 18:49:46.000000000 +0200
--- ./src/defaults.cc	2004-02-13 15:34:18.000000000 +0100
***************
*** 92,97 ****
--- 92,102 ----
  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,315 ****
--- 315,333 ----
      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,405 ****
--- 418,427 ----
  
    set_default_info_prog ();
  
+ #ifdef HAVE_FFTW3
+   set_default_wisdom_prog ();
+ #endif
+ 
    set_default_editor ();
  
    set_local_site_defaults_file ();
***************
*** 484,489 ****
--- 506,531 ----
    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,592 ****
--- 629,646 ----
  @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,627 ****
  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\
--- 675,681 ----
  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\
*** ./configure.in.orig	2004-01-28 15:16:49.000000000 +0100
--- ./configure.in	2004-02-10 10:50:54.000000000 +0100
***************
*** 407,425 ****
    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
    fi
  fi
  
! if test "$with_fftw" = yes; then
    FFT_DIR=''
!   AC_DEFINE(HAVE_FFTW, 1, [Define if the FFTW library is available.])
  fi
  
  WITH_MPI=true
--- 407,423 ----
    with_fftw=$withval, with_fftw=yes)
  
  if test "$with_fftw" = "yes"; then
!   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 && test "$with_fftw3" = yes; then
    FFT_DIR=''
!   AC_DEFINE(HAVE_FFTW3, 1, [Define if the FFTW3 library is used.])
  fi
  
  WITH_MPI=true
*** ./ChangeLog.orig	2004-01-30 15:25:15.000000000 +0100
--- ./ChangeLog	2004-02-16 17:06:47.000000000 +0100
***************
*** 1,3 ****
--- 1,8 ----
+ 2004-02-16  David Bateman <address@hidden>
+ 
+ 	* configure.in: Test for the presence of FFTW 3.x and use it in
+ 	preference to FFTW 2.x. Define HAVE_FFTW3
+ 	
  2004-01-22  John W. Eaton  <address@hidden>
  
  	* octMakefile.in (maintainer-clean, distclean):