rfftwnd_threads.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
 *
 */

#include "fftw_threads-int.h"
#include "rfftw_threads.h"

/****************** prototypes for rexec2_threads routines *****************/

extern void rfftw_real2c_threads_aux(fftw_plan plan, int howmany,
				     fftw_real *in, int istride, int idist,
				     fftw_complex *out, int ostride, int odist,
				     fftw_real *work,
				     int nthreads);
extern void rfftw_c2real_threads_aux(fftw_plan plan, int howmany,
				     fftw_complex *in, int istride, int idist,
				     fftw_real *out, int ostride, int odist,
				     fftw_real *work,
				     int nthreads);
extern void rfftw_real2c_overlap_threads_aux(fftw_plan plan, int howmany,
					     fftw_real *in, int istride,
					     int idist,
					     fftw_complex *out,
					     int ostride, int odist,
					     fftw_real *work,
					     int nthreads);
extern void rfftw_c2real_overlap_threads_aux(fftw_plan plan, int howmany,
					     fftw_complex *in,
					     int istride, int idist,
					     fftw_real *out, 
					     int ostride, int odist,
					     fftw_real *work,
					     int nthreads);

/********************** 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);

/****************** prototypes for rfftwnd routines *****************/

extern void rfftwnd_real2c_aux(fftwnd_plan p, int cur_dim,
			       fftw_real *in, int istride,
			       fftw_complex *out, int ostride,
			       fftw_real *work);
extern void rfftwnd_c2real_aux(fftwnd_plan p, int cur_dim,
			       fftw_complex *in, int istride,
			       fftw_real *out, int ostride,
			       fftw_real *work);
extern void rfftwnd_real2c_aux_howmany(fftwnd_plan p, int cur_dim,
                                int howmany,
                                fftw_real *in, int istride, int idist,
                                fftw_complex *out, int ostride, int odist,
                                fftw_complex *work);
extern void rfftwnd_c2real_aux_howmany(fftwnd_plan p, int cur_dim,
                                int howmany,
                                fftw_complex *in, int istride, int idist,
                                fftw_real *out, int ostride, int odist,
                                fftw_complex *work);


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

typedef struct {
     fftwnd_plan p;
     int cur_dim;
     void *in;
     int istride, idist;
     void *out;
     int ostride, odist;
     fftw_complex *work;
} aux_data;

static void *real2c_aux_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     aux_data *d = (aux_data *) ldata->data;
     fftwnd_plan p = d->p;
     int cur_dim = d->cur_dim;
     fftw_real *in = (fftw_real *) d->in;
     int istride = d->istride, idist = d->idist;
     fftw_complex *out = (fftw_complex *) d->out;
     int ostride = d->ostride, odist = d->odist;
     fftw_real *work = (fftw_real*) (d->work + p->nwork * ldata->thread_num);

     for (; min < max; ++min)
	  rfftwnd_real2c_aux(p, cur_dim, in + idist * min, istride,
			     out + odist * min, ostride, work);
     return 0;
}

void rfftwnd_real2c_threads_aux(fftwnd_plan p, int cur_dim,
				fftw_real *in, int istride,
				fftw_complex *out, int ostride,
				fftw_complex *work,
				int nthreads)
{
     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_threads_aux(p->plans[p->rank - 1], n,
					in, istride, (n_after * istride) * 2,
					out, istride, n_after * istride,
					(fftw_real *) work, nthreads);
	  else
	       rfftw_real2c_threads_aux(p->plans[p->rank - 1], n,
					in, istride,
					p->plans[p->rank - 1]->n * istride,
					out, ostride, n_after * ostride,
					(fftw_real *) work, nthreads);
     }
     else {    /* we have at least two dimensions to go */
	  int nr = p->plans[p->rank - 1]->n;
	  aux_data d;

	  d.p = p;
	  d.cur_dim = cur_dim + 1;
	  d.in = in;
	  d.istride = istride;
	  d.idist = istride * (p->is_in_place ? n_after * 2
			       : nr * (n_after / (nr/2 + 1)));
	  d.out = out;
	  d.ostride = ostride;
	  d.odist = ostride * n_after;
	  d.work = work;

	  fftw_thread_spawn_loop(n, nthreads, real2c_aux_thread, &d);
     }

     /* do the current dimension (in-place): */
     /* (Use internal function instead of fftw_threads so that we can
	pass our workspace array.) */
     fftw_executor_many_inplace_threads(p->plans[cur_dim]->n,
					out, work, p->plans[cur_dim]->root,
					n_after * ostride, n_after, ostride,
					nthreads);
}

static void *c2real_aux_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     aux_data *d = (aux_data *) ldata->data;
     fftwnd_plan p = d->p;
     int cur_dim = d->cur_dim;
     fftw_complex *in = (fftw_complex *) d->in;
     int istride = d->istride, idist = d->idist;
     fftw_real *out = (fftw_real *) d->out;
     int ostride = d->ostride, odist = d->odist;
     fftw_real *work = (fftw_real*) (d->work + p->nwork * ldata->thread_num);

     for (; min < max; ++min)
	  rfftwnd_c2real_aux(p, cur_dim, in + idist * min, istride,
			     out + odist * min, ostride, work);
     return 0;
}

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

     /* do the current dimension (in-place): */
     /* (Use internal function instead of fftw_threads so that we can
	pass our workspace array.) */
     fftw_executor_many_inplace_threads(p->plans[cur_dim]->n,
					in, work, p->plans[cur_dim]->root,
					n_after * istride, n_after, istride,
					nthreads);

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

	  d.p = p;
	  d.cur_dim = cur_dim + 1;
	  d.in = in;
	  d.istride = istride;
	  d.odist = ostride * (p->is_in_place ? n_after * 2
			       : nr * (n_after / (nr/2 + 1)));
	  d.out = out;
	  d.ostride = ostride;
	  d.idist = istride * n_after;
	  d.work = work;

	  fftw_thread_spawn_loop(n, nthreads, c2real_aux_thread, &d);
     }
}

typedef struct {
     fftw_plan p;
     int howmany;
     void *in;
     int istride, idist, idist0;
     void *out;
     int ostride, odist, odist0;
     fftw_real *work;
     int wdist;
} howmany_aux_data;

static void *r2c_overlap_howmany_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     howmany_aux_data *d = (howmany_aux_data *) ldata->data;
     fftw_plan p = d->p;
     int howmany = d->howmany;
     fftw_real *in = (fftw_real *) d->in;
     int istride = d->istride, idist = d->idist, idist0 = d->idist0;
     fftw_complex *out = (fftw_complex *) d->out;
     int ostride = d->ostride, odist = d->odist, odist0 = d->odist0;
     fftw_real *work = d->work + d->wdist * ldata->thread_num;

     for (; min < max; ++min)
	  rfftw_real2c_overlap_aux(p, howmany,
				   in + min * idist0, istride, idist,
				   out + min * odist0, ostride, odist,
				   work);

     return 0;
}

static void *c2r_overlap_howmany_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     howmany_aux_data *d = (howmany_aux_data *) ldata->data;
     fftw_plan p = d->p;
     int howmany = d->howmany;
     fftw_real *out = (fftw_real *) d->out;
     int istride = d->istride, idist = d->idist, idist0 = d->idist0;
     fftw_complex *in = (fftw_complex *) d->in;
     int ostride = d->ostride, odist = d->odist, odist0 = d->odist0;
     fftw_real *work = d->work + d->wdist * ldata->thread_num;

     for (; min < max; ++min)
	  rfftw_c2real_overlap_aux(p, howmany,
				   in + min * idist0, istride, idist,
				   out + min * odist0, ostride, odist,
				   work);

     return 0;
}

static void *r2c_howmany_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     howmany_aux_data *d = (howmany_aux_data *) ldata->data;
     fftw_plan p = d->p;
     int howmany = d->howmany;
     fftw_real *in = (fftw_real *) d->in;
     int istride = d->istride, idist = d->idist, idist0 = d->idist0;
     fftw_complex *out = (fftw_complex *) d->out;
     int ostride = d->ostride, odist = d->odist, odist0 = d->odist0;
     fftw_real *work = d->work + d->wdist * ldata->thread_num;

     for (; min < max; ++min)
	  rfftw_real2c_aux(p, howmany,
			   in + min * idist0, istride, idist,
			   out + min * odist0, ostride, odist,
			   work);

     return 0;
}

static void *c2r_howmany_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     howmany_aux_data *d = (howmany_aux_data *) ldata->data;
     fftw_plan p = d->p;
     int howmany = d->howmany;
     fftw_real *out = (fftw_real *) d->out;
     int istride = d->istride, idist = d->idist, idist0 = d->idist0;
     fftw_complex *in = (fftw_complex *) d->in;
     int ostride = d->ostride, odist = d->odist, odist0 = d->odist0;
     fftw_real *work = d->work + d->wdist * ldata->thread_num;

     for (; min < max; ++min)
	  rfftw_c2real_aux(p, howmany,
			   in + min * idist0, istride, idist,
			   out + min * odist0, ostride, odist,
			   work);

     return 0;
}

typedef struct {
     fftwnd_plan p;
     int cur_dim;
     int howmany;
     void *in;
     int istride, idist, idist0;
     void *out;
     int ostride, odist, odist0;
     fftw_complex *work;
     int wdist;
} howmany_hyperslab_aux_data;

static void *r2c_hyperslab_howmany_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     howmany_hyperslab_aux_data *d = (howmany_hyperslab_aux_data*) ldata->data;
     fftwnd_plan p = d->p;
     int cur_dim = d->cur_dim;
     int howmany = d->howmany;
     fftw_real *in = (fftw_real *) d->in;
     int istride = d->istride, idist = d->idist, idist0 = d->idist0;
     fftw_complex *out = (fftw_complex *) d->out;
     int ostride = d->ostride, odist = d->odist, odist0 = d->odist0;
     fftw_complex *work = d->work + d->wdist * ldata->thread_num;

     for (; min < max; ++min)
	  rfftwnd_real2c_aux_howmany(p, cur_dim, howmany,
				     in + min * idist0, istride, idist,
				     out + min * odist0, ostride, odist,
				     work);

     return 0;
}

static void *c2r_hyperslab_howmany_thread(fftw_loop_data *ldata)
{
     int min = ldata->min, max = ldata->max;
     howmany_hyperslab_aux_data *d = (howmany_hyperslab_aux_data*) ldata->data;
     fftwnd_plan p = d->p;
     int cur_dim = d->cur_dim;
     int howmany = d->howmany;
     fftw_real *out = (fftw_real *) d->out;