rfftwnd.c

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

				  in + i * n_after * istride, istride,
			   out + i * n_after_r * ostride, ostride, work);
     }
}

/*
 * alternate version of rfftwnd_aux -- this version pushes the howmany
 * loop down to the leaves of the computation, for greater locality
 * in cases where dist < stride.  It is also required for correctness
 * if in==out, and we must call a special version of the executor.
 * Note that work must point to 'howmany' copies of its data
 * if in == out. 
 */

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)
{
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
     int k;

     if (cur_dim == p->rank - 2) {
	  /* just do the last dimension directly: */
	  if (p->is_in_place)
	       for (k = 0; k < n; ++k)
		    rfftw_real2c_overlap_aux(p->plans[p->rank - 1], howmany,
					in + (k * n_after * istride) * 2,
					     istride, idist,
					   out + (k * n_after * ostride),
					     ostride, odist,
					     (fftw_real *) work);
	  else {
	       int nlast = p->plans[p->rank - 1]->n;
	       for (k = 0; k < n; ++k)
		    rfftw_real2c_aux(p->plans[p->rank - 1], howmany,
				     in + k * nlast * istride,
				     istride, idist,
				     out + k * n_after * ostride,
				     ostride, odist,
				     (fftw_real *) 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_howmany(p, cur_dim + 1, howmany,
			    in + i * n_after_r * istride, istride, idist,
			     out + i * n_after * ostride, ostride, odist,
					  work);
     }

     /* do the current dimension (in-place): */
     for (k = 0; k < n_after; ++k)
	  fftw(p->plans[cur_dim], howmany,
	       out + k * ostride, n_after * ostride, odist,
	       work, 1, 0);
}

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)
{
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
     int k;

     /* do the current dimension (in-place): */
     for (k = 0; k < n_after; ++k)
	  fftw(p->plans[cur_dim], howmany,
	       in + k * istride, n_after * istride, idist,
	       work, 1, 0);

     if (cur_dim == p->rank - 2) {
	  /* just do the last dimension directly: */
	  if (p->is_in_place)
	       for (k = 0; k < n; ++k)
		    rfftw_c2real_overlap_aux(p->plans[p->rank - 1], howmany,
					     in + (k * n_after * istride),
					     istride, idist,
				       out + (k * n_after * ostride) * 2,
					     ostride, odist,
					     (fftw_real *) work);
	  else {
	       int nlast = p->plans[p->rank - 1]->n;
	       for (k = 0; k < n; ++k)
		    rfftw_c2real_aux(p->plans[p->rank - 1], howmany,
				     in + k * n_after * istride,
				     istride, idist,
				     out + k * nlast * ostride,
				     ostride, odist,
				     (fftw_real *) 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_howmany(p, cur_dim + 1, howmany,
			      in + i * n_after * istride, istride, idist,
			   out + i * n_after_r * ostride, ostride, odist,
					  work);
     }
}

/********** Computing the N-Dimensional FFT: User-Visible Routines **********/

void rfftwnd_real_to_complex(fftwnd_plan p, int howmany,
			     fftw_real *in, int istride, int idist,
			     fftw_complex *out, int ostride, int odist)
{
     fftw_complex *work = p->work;
     int rank = p->rank;
     int free_work = 0;

     if (p->dir != FFTW_REAL_TO_COMPLEX)
	  fftw_die("rfftwnd_real_to_complex with complex-to-real plan");

#ifdef FFTW_DEBUG
     if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
	 && p->nwork && p->work)
	  fftw_die("bug with FFTW_THREADSAFE flag");
#endif

     if (p->is_in_place) {
	  ostride = istride;
	  odist = (idist == 1 && idist < istride) ? 1 : (idist / 2);  /* ugh */
	  out = (fftw_complex *) in;
	  if (howmany > 1 && istride > idist && rank > 0) {
	       int new_nwork;

	       new_nwork = p->n[rank - 1] * howmany;
	       if (new_nwork > p->nwork) {
		    work = (fftw_complex *)
			fftw_malloc(sizeof(fftw_complex) * new_nwork);
		    if (!work)
			 fftw_die("error allocating work array");
		    free_work = 1;
	       }
	  }
     }
     if (p->nwork && !work) {
	  work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork);
	  free_work = 1;
     }
     switch (rank) {
	 case 0:
	      break;
	 case 1:
	      if (p->is_in_place && howmany > 1 && istride > idist)
		   rfftw_real2c_overlap_aux(p->plans[0], howmany,
					    in, istride, idist,
					    out, ostride, odist,
					    (fftw_real *) work);
	      else
		   rfftw_real2c_aux(p->plans[0], howmany,
				    in, istride, idist,
				    out, ostride, odist,
				    (fftw_real *) work);
	      break;
	 default:		/* rank >= 2 */
	      {
		   if (howmany > 1 && ostride > odist)
			rfftwnd_real2c_aux_howmany(p, 0, howmany,
						   in, istride, idist,
						   out, ostride, odist,
						   work);
		   else {
			int i;

			for (i = 0; i < howmany; ++i)
			     rfftwnd_real2c_aux(p, 0,
						in + i * idist, istride,
						out + i * odist, ostride,
						(fftw_real *) work);
		   }
	      }
     }

     if (free_work)
	  fftw_free(work);
}

void rfftwnd_complex_to_real(fftwnd_plan p, int howmany,
			     fftw_complex *in, int istride, int idist,
			     fftw_real *out, int ostride, int odist)
{
     fftw_complex *work = p->work;
     int rank = p->rank;
     int free_work = 0;

     if (p->dir != FFTW_COMPLEX_TO_REAL)
	  fftw_die("rfftwnd_complex_to_real with real-to-complex plan");

#ifdef FFTW_DEBUG
     if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
	 && p->nwork && p->work)
	  fftw_die("bug with FFTW_THREADSAFE flag");
#endif

     if (p->is_in_place) {
	  ostride = istride;
	  odist = idist;
	  odist = (idist == 1 && idist < istride) ? 1 : (idist * 2);  /* ugh */
	  out = (fftw_real *) in;
	  if (howmany > 1 && istride > idist && rank > 0) {
	       int new_nwork = p->n[rank - 1] * howmany;
	       if (new_nwork > p->nwork) {
		    work = (fftw_complex *)
			fftw_malloc(sizeof(fftw_complex) * new_nwork);
		    if (!work)
			 fftw_die("error allocating work array");
		    free_work = 1;
	       }
	  }
     }
     if (p->nwork && !work) {
	  work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork);
	  free_work = 1;
     }
     switch (rank) {
	 case 0:
	      break;
	 case 1:
	      if (p->is_in_place && howmany > 1 && istride > idist)
		   rfftw_c2real_overlap_aux(p->plans[0], howmany,
					    in, istride, idist,
					    out, ostride, odist,
					    (fftw_real *) work);
	      else
		   rfftw_c2real_aux(p->plans[0], howmany,
				    in, istride, idist,
				    out, ostride, odist,
				    (fftw_real *) work);
	      break;
	 default:		/* rank >= 2 */
	      {
		   if (howmany > 1 && ostride > odist)
			rfftwnd_c2real_aux_howmany(p, 0, howmany,
						   in, istride, idist,
						   out, ostride, odist,
						   work);
		   else {
			int i;

			for (i = 0; i < howmany; ++i)
			     rfftwnd_c2real_aux(p, 0,
						in + i * idist, istride,
						out + i * odist, ostride,
						(fftw_real *) work);
		   }
	      }
     }

     if (free_work)
	  fftw_free(work);
}

void rfftwnd_one_real_to_complex(fftwnd_plan p,
				 fftw_real *in, fftw_complex *out)
{
     rfftwnd_real_to_complex(p, 1, in, 1, 1, out, 1, 1);
}

void rfftwnd_one_complex_to_real(fftwnd_plan p,
				 fftw_complex *in, fftw_real *out)
{
     rfftwnd_complex_to_real(p, 1, in, 1, 1, out, 1, 1);
}