fftwnd.c

FFTW, a collection of fast C routines to compute the Discrete Fourier Transform in one or more dime

/*
 * Copyright (c) 1997-1999, 2003 Massachusetts Institute of Technology
 *
 * This program 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 of the License, or
 * (at your option) any later version.
 *
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* $Id: fftwnd.c,v 1.44 2003/03/16 23:43:46 stevenj Exp $ */

#include "fftw-int.h"

/* the number of buffers to use for buffered transforms: */
#define FFTWND_NBUFFERS 8

/* the default number of buffers to use: */
#define FFTWND_DEFAULT_NBUFFERS 0

/* the number of "padding" elements between consecutive buffer lines */
#define FFTWND_BUFFER_PADDING 8

static void destroy_plan_array(int rank, fftw_plan *plans);

static void init_test_array(fftw_complex *arr, int stride, int n)
{
     int j;

     for (j = 0; j < n; ++j) {
	  c_re(arr[stride * j]) = 0.0;
	  c_im(arr[stride * j]) = 0.0;
     }
}

/*
 * Same as fftw_measure_runtime, except for fftwnd plan.
 */
double fftwnd_measure_runtime(fftwnd_plan plan,
			      fftw_complex *in, int istride,
			      fftw_complex *out, int ostride)
{
     fftw_time begin, end, start;
     double t, tmax, tmin;
     int i, iter;
     int n;
     int repeat;

     if (plan->rank == 0)
	  return 0.0;

     n = 1;
     for (i = 0; i < plan->rank; ++i)
	  n *= plan->n[i];

     iter = 1;

     for (;;) {
	  tmin = 1.0E10;
	  tmax = -1.0E10;
	  init_test_array(in, istride, n);

	  start = fftw_get_time();
	  /* repeat the measurement FFTW_TIME_REPEAT times */
	  for (repeat = 0; repeat < FFTW_TIME_REPEAT; ++repeat) {
	       begin = fftw_get_time();
	       for (i = 0; i < iter; ++i) {
		    fftwnd(plan, 1, in, istride, 0, out, ostride, 0);
	       }
	       end = fftw_get_time();

	       t = fftw_time_to_sec(fftw_time_diff(end, begin));
	       if (t < tmin)
		    tmin = t;
	       if (t > tmax)
		    tmax = t;

	       /* do not run for too long */
	       t = fftw_time_to_sec(fftw_time_diff(end, start));
	       if (t > FFTW_TIME_LIMIT)
		    break;
	  }

	  if (tmin >= FFTW_TIME_MIN)
	       break;

	  iter *= 2;
     }

     tmin /= (double) iter;
     tmax /= (double) iter;

     return tmin;
}

/********************** Initializing the FFTWND Plan ***********************/

/* Initialize everything except for the 1D plans and the work array: */
fftwnd_plan fftwnd_create_plan_aux(int rank, const int *n,
				   fftw_direction dir, int flags)
{
     int i;
     fftwnd_plan p;

     if (rank < 0)
	  return 0;

     for (i = 0; i < rank; ++i)
	  if (n[i] <= 0)
	       return 0;

     p = (fftwnd_plan) fftw_malloc(sizeof(fftwnd_data));
     p->n = 0;
     p->n_before = 0;
     p->n_after = 0;
     p->plans = 0;
     p->work = 0;
     p->dir = dir;

     p->rank = rank;
     p->is_in_place = flags & FFTW_IN_PLACE;

     p->nwork = 0;
     p->nbuffers = 0;

     if (rank == 0)
	  return 0;

     p->n = (int *) fftw_malloc(sizeof(int) * rank);
     p->n_before = (int *) fftw_malloc(sizeof(int) * rank);
     p->n_after = (int *) fftw_malloc(sizeof(int) * rank);
     p->n_before[0] = 1;
     p->n_after[rank - 1] = 1;

     for (i = 0; i < rank; ++i) {
	  p->n[i] = n[i];

	  if (i) {
	       p->n_before[i] = p->n_before[i - 1] * n[i - 1];
	       p->n_after[rank - 1 - i] = p->n_after[rank - i] * n[rank - i];
	  }
     }

     return p;
}

/* create an empty new array of rank 1d plans */
fftw_plan *fftwnd_new_plan_array(int rank)
{
     fftw_plan *plans;
     int i;

     plans = (fftw_plan *) fftw_malloc(rank * sizeof(fftw_plan));
     if (!plans)
	  return 0;
     for (i = 0; i < rank; ++i)
	  plans[i] = 0;
     return plans;
}

/* 
 * create an array of plans using the ordinary 1d fftw_create_plan,
 * which allocates its own array and creates plans optimized for
 * contiguous data. 
 */
fftw_plan *fftwnd_create_plans_generic(fftw_plan *plans,
				       int rank, const int *n,
				       fftw_direction dir, int flags)
{
     if (rank <= 0)
	  return 0;

     if (plans) {
	  int i, j;
	  int cur_flags;

	  for (i = 0; i < rank; ++i) {
	       if (i < rank - 1 || (flags & FFTW_IN_PLACE)) {
		    /* 
		     * fft's except the last dimension are always in-place 
		     */
		    cur_flags = flags | FFTW_IN_PLACE;
		    for (j = i - 1; j >= 0 && n[i] != n[j]; --j);
	       } else {
		    cur_flags = flags;
		    /* 
		     * we must create a separate plan for the last
		     * dimension 
		     */
		    j = -1;
	       }

	       if (j >= 0) {
		    /* 
		     * If a plan already exists for this size
		     * array, reuse it: 
		     */
		    plans[i] = plans[j];
	       } else {
		    /* generate a new plan: */
		    plans[i] = fftw_create_plan(n[i], dir, cur_flags);
		    if (!plans[i]) {
			 destroy_plan_array(rank, plans);
			 return 0;
		    }
	       }
	  }
     }
     return plans;
}

static int get_maxdim(int rank, const int *n, int flags)
{
     int i;
     int maxdim = 0;

     for (i = 0; i < rank - 1; ++i)
	  if (n[i] > maxdim)
	       maxdim = n[i];
     if (rank > 0 && flags & FFTW_IN_PLACE && n[rank - 1] > maxdim)
	  maxdim = n[rank - 1];

     return maxdim;
}

/* compute number of elements required for work array (has to
   be big enough to hold ncopies of the largest dimension in
   n that will need an in-place transform. */
int fftwnd_work_size(int rank, const int *n, int flags, int ncopies)
{
     return (ncopies * get_maxdim(rank, n, flags)
	     + (ncopies - 1) * FFTWND_BUFFER_PADDING);
}

/*
 * create plans using the fftw_create_plan_specific planner, which
 * allows us to create plans for each dimension that are specialized
 * for the strides that we are going to use. 
 */
fftw_plan *fftwnd_create_plans_specific(fftw_plan *plans,
					int rank, const int *n,
					const int *n_after,
					fftw_direction dir, int flags,
					fftw_complex *in, int istride,
					fftw_complex *out, int ostride)
{
     if (rank <= 0)
	  return 0;

     if (plans) {
	  int i, stride, cur_flags;
	  fftw_complex *work = 0;
	  int nwork;

	  nwork = fftwnd_work_size(rank, n, flags, 1);
	  if (nwork)
	       work = (fftw_complex*)fftw_malloc(nwork * sizeof(fftw_complex));

	  for (i = 0; i < rank; ++i) {
	       /* fft's except the last dimension are always in-place */
	       if (i < rank - 1)
		    cur_flags = flags | FFTW_IN_PLACE;
	       else
		    cur_flags = flags;

	       /* stride for transforming ith dimension */
	       stride = n_after[i];

	       if (cur_flags & FFTW_IN_PLACE)
		    plans[i] = fftw_create_plan_specific(n[i], dir, cur_flags,
						    in, istride * stride,
							 work, 1);
	       else
		    plans[i] = fftw_create_plan_specific(n[i], dir, cur_flags,
						    in, istride * stride,
						  out, ostride * stride);
	       if (!plans[i]) {
		    destroy_plan_array(rank, plans);
		    fftw_free(work);
		    return 0;
	       }
	  }

	  if (work)
	       fftw_free(work);
     }
     return plans;
}

/*
 * Create an fftwnd_plan specialized for specific arrays.  (These
 * arrays are ignored, however, if they are NULL or if the flags do
 * not include FFTW_MEASURE.)  The main advantage of being provided
 * arrays like this is that we can do runtime timing measurements of
 * our options, without worrying about allocating excessive scratch
 * space.
 */
fftwnd_plan fftwnd_create_plan_specific(int rank, const int *n,
					fftw_direction dir, int flags,
					fftw_complex *in, int istride,
					fftw_complex *out, int ostride)
{
     fftwnd_plan p;

     if (!(p = fftwnd_create_plan_aux(rank, n, dir, flags)))
	  return 0;

     if (!(flags & FFTW_MEASURE) || in == 0
	 || (!p->is_in_place && out == 0)) {

/**** use default plan ****/

	  p->plans = fftwnd_create_plans_generic(fftwnd_new_plan_array(rank),
						 rank, n, dir, flags);
	  if (!p->plans) {
	       fftwnd_destroy_plan(p);
	       return 0;
	  }
	  if (flags & FFTWND_FORCE_BUFFERED)
	       p->nbuffers = FFTWND_NBUFFERS;
	  else
	       p->nbuffers = FFTWND_DEFAULT_NBUFFERS;

	  p->nwork = fftwnd_work_size(rank, n, flags, p->nbuffers + 1);
	  if (p->nwork && !(flags & FFTW_THREADSAFE)) {
	       p->work = (fftw_complex*) fftw_malloc(p->nwork 
						     * sizeof(fftw_complex));
	       if (!p->work) {
		    fftwnd_destroy_plan(p);
		    return 0;
	       }
	  }
     } else {
/**** use runtime measurements to pick plan ****/

	  fftw_plan *plans_buf, *plans_nobuf;
	  double t_buf, t_nobuf;

	  p->nwork = fftwnd_work_size(rank, n, flags, FFTWND_NBUFFERS + 1);
	  if (p->nwork && !(flags & FFTW_THREADSAFE)) {
	       p->work = (fftw_complex*) fftw_malloc(p->nwork 
						     * sizeof(fftw_complex));
	       if (!p->work) {
		    fftwnd_destroy_plan(p);
		    return 0;
	       }
	  }
	  else
	       p->work = (fftw_complex*) NULL;

	  /* two possible sets of 1D plans: */
	  plans_buf = fftwnd_create_plans_generic(fftwnd_new_plan_array(rank),
						  rank, n, dir, flags);
	  plans_nobuf = 
	       fftwnd_create_plans_specific(fftwnd_new_plan_array(rank),
					    rank, n, p->n_after, dir,
					    flags, in, istride,
					    out, ostride);
	  if (!plans_buf || !plans_nobuf) {
	       destroy_plan_array(rank, plans_nobuf);
	       destroy_plan_array(rank, plans_buf);
	       fftwnd_destroy_plan(p);
	       return 0;
	  }
	  /* time the two possible plans */
	  p->plans = plans_nobuf;
	  p->nbuffers = 0;
	  p->nwork = fftwnd_work_size(rank, n, flags, p->nbuffers + 1);
	  t_nobuf = fftwnd_measure_runtime(p, in, istride, out, ostride);
	  p->plans = plans_buf;
	  p->nbuffers = FFTWND_NBUFFERS;
	  p->nwork = fftwnd_work_size(rank, n, flags, p->nbuffers + 1);
	  t_buf = fftwnd_measure_runtime(p, in, istride, out, ostride);

	  /* pick the better one: */
	  if (t_nobuf < t_buf) {	/* use unbuffered transform */
	       p->plans = plans_nobuf;
	       p->nbuffers = 0;

	       /* work array is unnecessarily large */
	       if (p->work)
		    fftw_free(p->work);
	       p->work = 0;

	       destroy_plan_array(rank, plans_buf);

	       /* allocate a work array of the correct size: */
	       p->nwork = fftwnd_work_size(rank, n, flags, p->nbuffers + 1);
	       if (p->nwork && !(flags & FFTW_THREADSAFE)) {
		    p->work = (fftw_complex*) fftw_malloc(p->nwork 
						       * sizeof(fftw_complex));
		    if (!p->work) {
			 fftwnd_destroy_plan(p);
			 return 0;
		    }