rfftwnd.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: rfftwnd.c,v 1.35 2003/03/16 23:43:46 stevenj Exp $ */

#include "fftw-int.h"
#include "rfftw.h"

/********************** prototypes for rexec2 routines **********************/

extern void rfftw_real2c_aux(fftw_plan plan, int howmany,
			     fftw_real *in, int istride, int idist,
			     fftw_complex *out, int ostride, int odist,
			     fftw_real *work);
extern void rfftw_c2real_aux(fftw_plan plan, int howmany,
			     fftw_complex *in, int istride, int idist,
			     fftw_real *out, int ostride, int odist,
			     fftw_real *work);
extern void rfftw_real2c_overlap_aux(fftw_plan plan, int howmany,
				   fftw_real *in, int istride, int idist,
			       fftw_complex *out, int ostride, int odist,
				     fftw_real *work);
extern void rfftw_c2real_overlap_aux(fftw_plan plan, int howmany,
				fftw_complex *in, int istride, int idist,
				  fftw_real *out, int ostride, int odist,
				     fftw_real *work);

/********************** Initializing the RFFTWND Plan ***********************/

/*
 * 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 rfftwnd_create_plan_specific(int rank, const int *n,
					 fftw_direction dir, int flags,
					 fftw_real *in, int istride,
					 fftw_real *out, int ostride)
{
     fftwnd_plan p;
     int i;
     int rflags = flags & ~FFTW_IN_PLACE;
     /* note that we always do rfftw transforms out-of-place in rexec2.c */

     if (flags & FFTW_IN_PLACE) {
	  out = NULL;
	  ostride = istride;
     }
     istride = ostride = 1;	/* 
				 * strides don't work yet, since it is not 
				 * clear whether they apply to real 
				 * or complex data 
				 */

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

     for (i = 0; i < rank - 1; ++i)
	  p->n_after[i] = (n[rank - 1]/2 + 1) * (p->n_after[i] / n[rank - 1]);
     if (rank > 0)
	  p->n[rank - 1] = n[rank - 1] / 2 + 1;

     p->plans = fftwnd_new_plan_array(rank);
     if (rank > 0 && !p->plans) {
	  rfftwnd_destroy_plan(p);
	  return 0;
     }
     if (rank > 0) {
	  p->plans[rank - 1] = rfftw_create_plan(n[rank - 1], dir, rflags);
	  if (!p->plans[rank - 1]) {
	       rfftwnd_destroy_plan(p);
	       return 0;
	  }
     }
     if (rank > 1) {
	  if (!(flags & FFTW_MEASURE) || in == 0
	      || (!p->is_in_place && out == 0)) {
	       if (!fftwnd_create_plans_generic(p->plans, rank - 1, n,
					   dir, flags | FFTW_IN_PLACE)) {
		    rfftwnd_destroy_plan(p);
		    return 0;
	       }
	  } else if (dir == FFTW_COMPLEX_TO_REAL || (flags & FFTW_IN_PLACE)) {
	       if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
						 p->n_after,
					      dir, flags | FFTW_IN_PLACE,
						 (fftw_complex *) in,
						 istride,
						 0, 0)) {
		    rfftwnd_destroy_plan(p);
		    return 0;
	       }
	  } else {
	       if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
						 p->n_after,
					      dir, flags | FFTW_IN_PLACE,
						 (fftw_complex *) out,
						 ostride,
						 0, 0)) {
		    rfftwnd_destroy_plan(p);
		    return 0;
	       }
	  }
     }
     p->nbuffers = 0;
     p->nwork = fftwnd_work_size(rank, p->n, flags | FFTW_IN_PLACE,
				 p->nbuffers + 1);
     if (p->nwork && !(flags & FFTW_THREADSAFE)) {
	  p->work = (fftw_complex *) fftw_malloc(p->nwork
						 * sizeof(fftw_complex));
	  if (!p->work) {
	       rfftwnd_destroy_plan(p);
	       return 0;
	  }
     }
     return p;
}

fftwnd_plan rfftw2d_create_plan_specific(int nx, int ny,
					 fftw_direction dir, int flags,
					 fftw_real *in, int istride,
					 fftw_real *out, int ostride)
{
     int n[2];

     n[0] = nx;
     n[1] = ny;

     return rfftwnd_create_plan_specific(2, n, dir, flags,
					 in, istride, out, ostride);
}

fftwnd_plan rfftw3d_create_plan_specific(int nx, int ny, int nz,
					 fftw_direction dir, int flags,
					 fftw_real *in, int istride,
					 fftw_real *out, int ostride)
{
     int n[3];

     n[0] = nx;
     n[1] = ny;
     n[2] = nz;

     return rfftwnd_create_plan_specific(3, n, dir, flags,
					 in, istride, out, ostride);
}

/* Create a generic fftwnd plan: */

fftwnd_plan rfftwnd_create_plan(int rank, const int *n,
				fftw_direction dir, int flags)
{
     return rfftwnd_create_plan_specific(rank, n, dir, flags, 0, 1, 0, 1);
}

fftwnd_plan rfftw2d_create_plan(int nx, int ny,
				fftw_direction dir, int flags)
{
     return rfftw2d_create_plan_specific(nx, ny, dir, flags, 0, 1, 0, 1);
}

fftwnd_plan rfftw3d_create_plan(int nx, int ny, int nz,
				fftw_direction dir, int flags)
{
     return rfftw3d_create_plan_specific(nx, ny, nz, dir, flags, 0, 1, 0, 1);
}

/************************ Freeing the RFFTWND Plan ************************/

void rfftwnd_destroy_plan(fftwnd_plan plan)
{
     fftwnd_destroy_plan(plan);
}

/************************ Printing the RFFTWND Plan ************************/

void rfftwnd_fprint_plan(FILE *f, fftwnd_plan plan)
{
     fftwnd_fprint_plan(f, plan);
}

void rfftwnd_print_plan(fftwnd_plan plan)
{
     rfftwnd_fprint_plan(stdout, plan);
}

/*********** Computing the N-Dimensional FFT: Auxiliary Routines ************/

void rfftwnd_real2c_aux(fftwnd_plan p, int cur_dim,
			fftw_real *in, int istride,
			fftw_complex *out, int ostride,
			fftw_real *work)
{
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];

     if (cur_dim == p->rank - 2) {
	  /* just do the last dimension directly: */
	  if (p->is_in_place)
	       rfftw_real2c_aux(p->plans[p->rank - 1], n,
				in, istride, (n_after * istride) * 2,
				out, istride, n_after * istride,
				work);
	  else
	       rfftw_real2c_aux(p->plans[p->rank - 1], n,
			 in, istride, p->plans[p->rank - 1]->n * istride,
				out, ostride, n_after * ostride,
				work);
     } else {			/* we have at least two dimensions to go */
	  int nr = p->plans[p->rank - 1]->n;
	  int n_after_r = p->is_in_place ? n_after * 2 
	       : nr * (n_after / (nr/2 + 1));
	  int i;

	  /* 
	   * process the subsequent dimensions recursively, in hyperslabs,
	   * to get maximum locality: 
	   */
	  for (i = 0; i < n; ++i)
	       rfftwnd_real2c_aux(p, cur_dim + 1,
				  in + i * n_after_r * istride, istride,
			     out + i * n_after * ostride, ostride, work);
     }

     /* do the current dimension (in-place): */
     fftw(p->plans[cur_dim], n_after,
	  out, n_after * ostride, ostride,
	  (fftw_complex *) work, 1, 0);
     /* I hate this cast */
}

void rfftwnd_c2real_aux(fftwnd_plan p, int cur_dim,
			fftw_complex *in, int istride,
			fftw_real *out, int ostride,
			fftw_real *work)
{
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];

     /* do the current dimension (in-place): */
     fftw(p->plans[cur_dim], n_after,
	  in, n_after * istride, istride,
	  (fftw_complex *) work, 1, 0);

     if (cur_dim == p->rank - 2) {
	  /* just do the last dimension directly: */
	  if (p->is_in_place)
	       rfftw_c2real_aux(p->plans[p->rank - 1], n,
				in, istride, n_after * istride,
				out, istride, (n_after * istride) * 2,
				work);
	  else
	       rfftw_c2real_aux(p->plans[p->rank - 1], n,
				in, istride, n_after * istride,
			out, ostride, p->plans[p->rank - 1]->n * ostride,
				work);
     } else {			/* we have at least two dimensions to go */
	  int nr = p->plans[p->rank - 1]->n;
	  int n_after_r = p->is_in_place ? n_after * 2 : 
	       nr * (n_after / (nr/2 + 1));
	  int i;

	  /* 
	   * process the subsequent dimensions recursively, in hyperslabs,
	   * to get maximum locality: 
	   */
	  for (i = 0; i < n; ++i)
	       rfftwnd_c2real_aux(p, cur_dim + 1,