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

/*
 * planner.c -- find the optimal plan
 */

/* $Id: rplanner.c,v 1.25 2003/03/16 23:43:46 stevenj Exp $ */
#ifdef FFTW_USING_CILK
#include <cilk.h>
#include <cilk-compat.h>
#endif

#include <stdlib.h>
#include <stdio.h>

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

extern fftw_codelet_desc *rfftw_config[];	/* global from rconfig.c */
extern fftw_rgeneric_codelet fftw_hc2hc_forward_generic;
extern fftw_rgeneric_codelet fftw_hc2hc_backward_generic;

fftw_plan_hook_ptr rfftw_plan_hook = (fftw_plan_hook_ptr) NULL;

/* timing rfftw plans: */
static double rfftw_measure_runtime(fftw_plan plan,
				    fftw_real *in, int istride,
				    fftw_real *out, int ostride)
{
     fftw_time begin, end, start;
     double t, tmin;
     int i, iter;
     int n;
     int repeat;
     int howmany = plan->vector_size;

     n = plan->n;

     iter = 1;

     for (;;) {
	  tmin = 1.0E10;
	  for (i = 0; i < n * howmany; ++i)
	       in[istride * i] = 0.0;

	  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)
		    rfftw(plan, howmany, in, istride, istride,
			  out, ostride, ostride);
	       end = fftw_get_time();

	       t = fftw_time_to_sec(fftw_time_diff(end, begin));
	       if (t < tmin)
		    tmin = 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;
     return tmin;
}

/* auxiliary functions */
static void rcompute_cost(fftw_plan plan,
			  fftw_real *in, int istride,
			  fftw_real *out, int ostride)
{
     if (plan->flags & FFTW_MEASURE)
	  plan->cost = rfftw_measure_runtime(plan, in, istride, out, ostride);
     else {
	  double c;
	  c = plan->n * fftw_estimate_node(plan->root) * plan->vector_size;
	  plan->cost = c;
     }
}

static void run_plan_hooks(fftw_plan p)
{
     if (rfftw_plan_hook && p) {
	  fftw_complete_twiddle(p->root, p->n);
	  rfftw_plan_hook(p);
     }
}

/* macrology */
#define FOR_ALL_RCODELETS(p) \
   fftw_codelet_desc **__q, *p;                         \
   for (__q = &rfftw_config[0]; (p = (*__q)); ++__q)

/******************************************
 *      Recursive planner                 *
 ******************************************/
static fftw_plan rplanner(fftw_plan *table, int n,
			  fftw_direction dir, int flags, int vector_size,
			  fftw_real *, int, fftw_real *, int);

/*
 * the planner consists of two parts: one that tries to
 * use accumulated wisdom, and one that does not.
 * A small driver invokes both parts in sequence
 */

/* planner with wisdom: look up the codelet suggested by the wisdom */
static fftw_plan rplanner_wisdom(fftw_plan *table, int n,
				 fftw_direction dir, int flags,
				 int vector_size,
				 fftw_real *in, int istride,
				 fftw_real *out, int ostride)
{
     fftw_plan best = (fftw_plan) 0;
     fftw_plan_node *node;
     int have_wisdom;
     enum fftw_node_type wisdom_type;
     int wisdom_signature;
     fftw_recurse_kind wisdom_recurse_kind;

     /* see if we remember any wisdom for this case */
     have_wisdom = fftw_wisdom_lookup(n, flags, dir, RFFTW_WISDOM,
				      istride, ostride,
				      &wisdom_type, &wisdom_signature,
				      &wisdom_recurse_kind, 0);

     if (!have_wisdom)
	  return best;

     if (wisdom_type == FFTW_REAL2HC || wisdom_type == FFTW_HC2REAL) {
	  FOR_ALL_RCODELETS(p) {
	       if (p->dir == dir && p->type == wisdom_type) {
		    /* see if wisdom applies */
		    if (wisdom_signature == p->signature &&
			p->size == n) {
			 if (wisdom_type == FFTW_REAL2HC)
			      node = fftw_make_node_real2hc(n, p);
			 else
			      node = fftw_make_node_hc2real(n, p);
			 best = fftw_make_plan(n, dir, node, flags,
					       p->type, p->signature,
					       FFTW_NORMAL_RECURSE,
					       vector_size);
			 fftw_use_plan(best);
			 run_plan_hooks(best);
			 return best;
		    }
	       }
	  }
     }
     if (wisdom_type == FFTW_HC2HC) {
	  FOR_ALL_RCODELETS(p) {
	       if (p->dir == dir && p->type == wisdom_type) {

		    /* see if wisdom applies */
		    if (wisdom_signature == p->signature &&
			p->size > 1 &&
			(n % p->size) == 0) {
			 fftw_plan r = rplanner(table, n / p->size, dir,
						flags | FFTW_NO_VECTOR_RECURSE,
						wisdom_recurse_kind ==
						FFTW_VECTOR_RECURSE ?
						p->size : vector_size,
						in, istride, out,
						ostride);
			 if (!r)
			      continue;
			 node = fftw_make_node_hc2hc(n, dir, p,
						     r->root, flags);
			 best = fftw_make_plan(n, dir, node, flags,
					       p->type, p->signature,
					       wisdom_recurse_kind,
					       vector_size);
			 fftw_use_plan(best);
			 run_plan_hooks(best);
			 fftw_destroy_plan_internal(r);
			 return best;
		    }
	       }
	  }
     }
     /* 
      * BUG (or: TODO)  Can we have generic wisdom? This is probably
      * an academic question
      */

     return best;
}

/*
 * planner with no wisdom: try all combinations and pick
 * the best
 */

static fftw_plan rplanner_normal(fftw_plan *table, int n, fftw_direction dir,
				 int flags, int vector_size,
				 fftw_real *in, int istride,
				 fftw_real *out, int ostride)
{
     fftw_plan best = (fftw_plan) 0;
     fftw_plan newplan;
     fftw_plan_node *node;

     /* see if we have any codelet that solves the problem */
     {
	  FOR_ALL_RCODELETS(p) {
	       if (p->dir == dir &&
		   (p->type == FFTW_REAL2HC || p->type == FFTW_HC2REAL)) {
		    if (p->size == n) {
			 if (p->type == FFTW_REAL2HC)
			      node = fftw_make_node_real2hc(n, p);
			 else
			      node = fftw_make_node_hc2real(n, p);
			 newplan = fftw_make_plan(n, dir, node, flags,
						  p->type, p->signature,
						  FFTW_NORMAL_RECURSE,
                                                  vector_size);
			 fftw_use_plan(newplan);
			 run_plan_hooks(newplan);
			 rcompute_cost(newplan, in, istride, out, ostride);
			 best = fftw_pick_better(newplan, best);
		    }
	       }
	  }
     }

     /* Then, try all available twiddle codelets */
     {
	  FOR_ALL_RCODELETS(p) {
	       if (p->dir == dir && p->type == FFTW_HC2HC) {
		    if ((n % p->size) == 0 &&
			p->size > 1 &&
			(!best || n != p->size)) {
			 fftw_plan r = rplanner(table, n / p->size, dir,
						flags | FFTW_NO_VECTOR_RECURSE,
						vector_size,
						in, istride, out, ostride);
			 if (!r)
			      continue;
			 node = fftw_make_node_hc2hc(n, dir, p,
						     r->root, flags);
			 newplan = fftw_make_plan(n, dir, node, flags,
						  p->type, p->signature,
						  FFTW_NORMAL_RECURSE,
                                                  vector_size);
			 fftw_use_plan(newplan);
			 run_plan_hooks(newplan);
			 fftw_destroy_plan_internal(r);
			 rcompute_cost(newplan, in, istride, out, ostride);
			 best = fftw_pick_better(newplan, best);
		    }
	       }
	  }
     }

     /* try vector recursion unless prohibited by the flags: */
     if (! (flags & FFTW_NO_VECTOR_RECURSE)) {
	  FOR_ALL_RCODELETS(p) {
	       if (p->dir == dir && p->type == FFTW_HC2HC) {
		    if ((n % p->size) == 0 &&
			p->size > 1 &&
			(!best || n != p->size)) {
			 fftw_plan r = rplanner(table, n / p->size, dir, 
						flags | FFTW_NO_VECTOR_RECURSE,
						p->size,
						in, istride, out, ostride);
			 if (!r)
			      continue;
			 node = fftw_make_node_hc2hc(n, dir, p,
						     r->root, flags);
			 newplan = fftw_make_plan(n, dir, node, flags,
						  p->type, p->signature,
						  FFTW_VECTOR_RECURSE,
                                                  vector_size);
			 fftw_use_plan(newplan);
			 run_plan_hooks(newplan);
			 fftw_destroy_plan_internal(r);
			 rcompute_cost(newplan, in, istride, out, ostride);
			 best = fftw_pick_better(newplan, best);
		    }
	       }
	  }
     }

     /* 
      * Resort to generic codelets for unknown factors, but only if
      * n is odd--the rgeneric codelets can't handle even n's.
      */
     if (n % 2 != 0) {
	  fftw_rgeneric_codelet *codelet = (dir == FFTW_FORWARD ?
					    fftw_hc2hc_forward_generic :
					    fftw_hc2hc_backward_generic);
	  int size, prev_size = 0, remaining_factors = n;
	  fftw_plan r;

	  while (remaining_factors > 1) {
	       size = fftw_factor(remaining_factors);
	       remaining_factors /= size;

	       /* don't try the same factor more than once */
	       if (size == prev_size)
		    continue;
	       prev_size = size;

	       /* Look for codelets corresponding to this factor. */
	       {
		    FOR_ALL_RCODELETS(p) {
			 if (p->dir == dir && p->type == FFTW_HC2HC
			     && p->size == size) {
			      size = 0;
			      break;
			 }
		    }
	       }

	       /* 
	        * only try a generic/rader codelet if there were no 
	        * twiddle codelets for this factor 
	        */
	       if (!size)
		    continue;

	       r = rplanner(table, n / size, dir,
			    flags | FFTW_NO_VECTOR_RECURSE,
			    vector_size,
			    in, istride, out, ostride);

	       node = fftw_make_node_rgeneric(n, size, dir, codelet,
					      r->root, flags);
	       newplan = fftw_make_plan(n, dir, node, flags, FFTW_RGENERIC, 0,
					FFTW_NORMAL_RECURSE, vector_size);
	       fftw_use_plan(newplan);
	       run_plan_hooks(newplan);
	       fftw_destroy_plan_internal(r);
	       rcompute_cost(newplan, in, istride, out, ostride);
	       best = fftw_pick_better(newplan, best);
	  }
     }
     return best;
}

static fftw_plan rplanner(fftw_plan *table, int n, fftw_direction dir,
			  int flags, int vector_size,
			  fftw_real *in, int istride,
			  fftw_real *out, int ostride)
{
     fftw_plan best = (fftw_plan) 0;

     if (vector_size > 1)
          flags |= FFTW_NO_VECTOR_RECURSE;

     /* see if plan has already been computed */
     best = fftw_lookup(table, n, flags, vector_size);
     if (best) {
	  fftw_use_plan(best);
	  return best;
     }
     /* try a wise plan */
     best = rplanner_wisdom(table, n, dir, flags, vector_size,
			    in, istride, out, ostride);

     if (!best) {
	  /* No wisdom.  Plan normally. */
	  best = rplanner_normal(table, n, dir, flags,
				 vector_size,
				 in, istride, out, ostride);
     }
     if (best) {
	  fftw_insert(table, best);

	  /* remember the wisdom */
	  fftw_wisdom_add(n, flags, dir, RFFTW_WISDOM,
			  istride, ostride,
			  best->wisdom_type,
			  best->wisdom_signature,
			  best->recurse_kind);
     }
     return best;
}

fftw_plan rfftw_create_plan_specific(int n, fftw_direction dir, int flags,
				     fftw_real *in, int istride,
				     fftw_real *out, int ostride)
{
     fftw_plan table;
     fftw_plan p1;

     /* validate parameters */
     if (n <= 0)
	  return (fftw_plan) 0;

#ifndef FFTW_ENABLE_VECTOR_RECURSE
     /* TEMPORARY: disable vector recursion until it is more tested. */
     flags |= FFTW_NO_VECTOR_RECURSE;
#endif

     if ((dir != FFTW_FORWARD) && (dir != FFTW_BACKWARD))
	  return (fftw_plan) 0;

     fftw_make_empty_table(&table);
     p1 = rplanner(&table, n, dir, flags, 1,
		   in, istride, out, ostride);
     fftw_destroy_table(&table);

     if (p1)
	  fftw_complete_twiddle(p1->root, n);
     return p1;
}

fftw_plan rfftw_create_plan(int n, fftw_direction dir, int flags)
{
     fftw_real *tmp_in;
     fftw_real *tmp_out;
     fftw_plan p;

     if (flags & FFTW_MEASURE) {
	  tmp_in = (fftw_real *) fftw_malloc(2 * n * sizeof(fftw_real));
	  if (!tmp_in)
	       return 0;
	  tmp_out = tmp_in + n;

	  p = rfftw_create_plan_specific(n, dir, flags,
					 tmp_in, 1, tmp_out, 1);

	  fftw_free(tmp_in);
     } else
	  p = rfftw_create_plan_specific(n, dir, flags,
				 (fftw_real *) 0, 1, (fftw_real *) 0, 1);

     return p;
}

void rfftw_destroy_plan(fftw_plan plan)
{
     fftw_destroy_plan_internal(plan);
}

void rfftw_fprint_plan(FILE *f, fftw_plan p)
{
     fftw_fprint_plan(f, p);
}

void rfftw_print_plan(fftw_plan p)
{
     rfftw_fprint_plan(stdout, p);
}