WebSVN - shark - Blame - Rev 2 - /shark/branches/pj/ports/fftw/rfftw/rplanner.c

Rev	Author	Line No.	Line
2	pj	1	/*
		2	* Copyright (c) 1997-1999 Massachusetts Institute of Technology
		3	*
		4	* This program is free software; you can redistribute it and/or modify
		5	* it under the terms of the GNU General Public License as published by
		6	* the Free Software Foundation; either version 2 of the License, or
		7	* (at your option) any later version.
		8	*
		9	* This program is distributed in the hope that it will be useful,
		10	* but WITHOUT ANY WARRANTY; without even the implied warranty of
		11	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
		12	* GNU General Public License for more details.
		13	*
		14	* You should have received a copy of the GNU General Public License
		15	* along with this program; if not, write to the Free Software
		16	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
		17	*
		18	*/
		19
		20	/*
		21	* planner.c -- find the optimal plan
		22	*/
		23
		24	/* $Id: rplanner.c,v 1.1.1.1 2002-03-29 14:12:59 pj Exp $ */
		25	#ifdef FFTW_USING_CILK
		26	#include <cilk.h>
		27	#include <cilk-compat.h>
		28	#endif
		29
		30	#include <stdlib.h>
		31	//#include <stdio.h>
		32
		33	#include <ports/fftw-int.h>
		34	#include <ports/rfftw.h>
		35
		36	extern fftw_codelet_desc rfftw_config[]; / global from rconfig.c */
		37	extern fftw_rgeneric_codelet fftw_hc2hc_forward_generic;
		38	extern fftw_rgeneric_codelet fftw_hc2hc_backward_generic;
		39
		40	/* undocumented debugging hook */
		41	void (rfftw_plan_hook)(fftw_plan plan) = (void ()(fftw_plan)) 0;
		42
		43	/* timing rfftw plans: */
		44	static double rfftw_measure_runtime(fftw_plan plan,
		45	fftw_real *in, int istride,
		46	fftw_real *out, int ostride)
		47	{
		48	fftw_time begin, end, start;
		49	double t, tmin;
		50	int i, iter;
		51	int n;
		52	int repeat;
		53
		54	n = plan->n;
		55
		56	iter = 1;
		57
		58	for (;;) {
		59	tmin = 1.0E10;
		60	for (i = 0; i < n; ++i)
		61	in[istride * i] = 0.0;
		62
		63	start = fftw_get_time();
		64	/* repeat the measurement FFTW_TIME_REPEAT times */
		65	for (repeat = 0; repeat < FFTW_TIME_REPEAT; ++repeat) {
		66	begin = fftw_get_time();
		67	for (i = 0; i < iter; ++i)
		68	rfftw(plan, 1, in, istride, 0, out, ostride, 0);
		69	end = fftw_get_time();
		70
		71	t = fftw_time_to_sec(fftw_time_diff(end, begin));
		72	if (t < tmin)
		73	tmin = t;
		74
		75	/* do not run for too long */
		76	t = fftw_time_to_sec(fftw_time_diff(end, start));
		77	if (t > FFTW_TIME_LIMIT)
		78	break;
		79	}
		80
		81	if (tmin >= FFTW_TIME_MIN)
		82	break;
		83
		84	iter *= 2;
		85	}
		86
		87	tmin /= (double) iter;
		88	return tmin;
		89	}
		90
		91	/* auxiliary functions */
		92	static void rcompute_cost(fftw_plan plan,
		93	fftw_real *in, int istride,
		94	fftw_real *out, int ostride)
		95	{
		96	if (plan->flags & FFTW_MEASURE)
		97	plan->cost = rfftw_measure_runtime(plan, in, istride, out, ostride);
		98	else {
		99	double c;
		100	c = plan->n * fftw_estimate_node(plan->root);
		101	plan->cost = c;
		102	}
		103	}
		104
		105	static void run_plan_hooks(fftw_plan p)
		106	{
		107	if (rfftw_plan_hook && p) {
		108	fftw_complete_twiddle(p->root, p->n);
		109	rfftw_plan_hook(p);
		110	}
		111	}
		112
		113	/* macrology */
		114	#define FOR_ALL_RCODELETS(p) \
		115	fftw_codelet_desc *__q, p; \
		116	for (__q = &rfftw_config[0]; (p = (*__q)); ++__q)
		117
		118	/******************************************
		119	* Recursive planner *
		120	******************************************/
		121	static fftw_plan rplanner(fftw_plan *table, int n,
		122	fftw_direction dir, int flags,
		123	fftw_real , int, fftw_real , int);
		124
		125	/*
		126	* the planner consists of two parts: one that tries to
		127	* use accumulated wisdom, and one that does not.
		128	* A small driver invokes both parts in sequence
		129	*/
		130
		131	/* planner with wisdom: look up the codelet suggested by the wisdom */
		132	static fftw_plan rplanner_wisdom(fftw_plan *table, int n,
		133	fftw_direction dir, int flags,
		134	fftw_real *in, int istride,
		135	fftw_real *out, int ostride)
		136	{
		137	fftw_plan best = (fftw_plan) 0;
		138	fftw_plan_node *node;
		139	int have_wisdom;
		140	enum fftw_node_type wisdom_type;
		141	int wisdom_signature;
		142
		143	/* see if we remember any wisdom for this case */
		144	have_wisdom = fftw_wisdom_lookup(n, flags, dir, RFFTW_WISDOM,
		145	istride, ostride,
		146	&wisdom_type, &wisdom_signature, 0);
		147
		148	if (!have_wisdom)
		149	return best;
		150
		151	if (wisdom_type == FFTW_REAL2HC \|\| wisdom_type == FFTW_HC2REAL) {
		152	FOR_ALL_RCODELETS(p) {
		153	if (p->dir == dir && p->type == wisdom_type) {
		154	/* see if wisdom applies */
		155	if (wisdom_signature == p->signature &&
		156	p->size == n) {
		157	if (wisdom_type == FFTW_REAL2HC)
		158	node = fftw_make_node_real2hc(n, p);
		159	else
		160	node = fftw_make_node_hc2real(n, p);
		161	best = fftw_make_plan(n, dir, node, flags,
		162	p->type, p->signature);
		163	fftw_use_plan(best);
		164	run_plan_hooks(best);
		165	return best;
		166	}
		167	}
		168	}
		169	}
		170	if (wisdom_type == FFTW_HC2HC) {
		171	FOR_ALL_RCODELETS(p) {
		172	if (p->dir == dir && p->type == wisdom_type) {
		173
		174	/* see if wisdom applies */
		175	if (wisdom_signature == p->signature &&
		176	p->size > 1 &&
		177	(n % p->size) == 0) {
		178	fftw_plan r = rplanner(table, n / p->size, dir,
		179	flags, in, istride, out,
		180	ostride);
		181	if (!r)
		182	continue;
		183	node = fftw_make_node_hc2hc(n, dir, p,
		184	r->root, flags);
		185	best = fftw_make_plan(n, dir, node, flags,
		186	p->type, p->signature);
		187	fftw_use_plan(best);
		188	run_plan_hooks(best);
		189	fftw_destroy_plan_internal(r);
		190	return best;
		191	}
		192	}
		193	}
		194	}
		195	/*
		196	* BUG (or: TODO) Can we have generic wisdom? This is probably
		197	* an academic question
		198	*/
		199
		200	return best;
		201	}
		202
		203	/*
		204	* planner with no wisdom: try all combinations and pick
		205	* the best
		206	*/
		207
		208	static fftw_plan rplanner_normal(fftw_plan *table, int n, fftw_direction dir,
		209	int flags, fftw_real *in, int istride,
		210	fftw_real *out, int ostride)
		211	{
		212	fftw_plan best = (fftw_plan) 0;
		213	fftw_plan newplan;
		214	fftw_plan_node *node;
		215
		216	/* see if we have any codelet that solves the problem */
		217	{
		218	FOR_ALL_RCODELETS(p) {
		219	if (p->dir == dir &&
		220	(p->type == FFTW_REAL2HC \|\| p->type == FFTW_HC2REAL)) {
		221	if (p->size == n) {
		222	if (p->type == FFTW_REAL2HC)
		223	node = fftw_make_node_real2hc(n, p);
		224	else
		225	node = fftw_make_node_hc2real(n, p);
		226	newplan = fftw_make_plan(n, dir, node, flags,
		227	p->type, p->signature);
		228	fftw_use_plan(newplan);
		229	run_plan_hooks(newplan);
		230	rcompute_cost(newplan, in, istride, out, ostride);
		231	best = fftw_pick_better(newplan, best);
		232	}
		233	}
		234	}
		235	}
		236
		237	/* Then, try all available twiddle codelets */
		238	{
		239	FOR_ALL_RCODELETS(p) {
		240	if (p->dir == dir && p->type == FFTW_HC2HC) {
		241	if ((n % p->size) == 0 &&
		242	p->size > 1 &&
		243	(!best \|\| n != p->size)) {
		244	fftw_plan r = rplanner(table, n / p->size, dir, flags,
		245	in, istride, out, ostride);
		246	if (!r)
		247	continue;
		248	node = fftw_make_node_hc2hc(n, dir, p,
		249	r->root, flags);
		250	newplan = fftw_make_plan(n, dir, node, flags,
		251	p->type, p->signature);
		252	fftw_use_plan(newplan);
		253	run_plan_hooks(newplan);
		254	fftw_destroy_plan_internal(r);
		255	rcompute_cost(newplan, in, istride, out, ostride);
		256	best = fftw_pick_better(newplan, best);
		257	}
		258	}
		259	}
		260	}
		261
		262	/*
		263	* Resort to generic codelets for unknown factors, but only if
		264	* n is odd--the rgeneric codelets can't handle even n's.
		265	*/
		266	if (n % 2 != 0) {
		267	fftw_rgeneric_codelet *codelet = (dir == FFTW_FORWARD ?
		268	fftw_hc2hc_forward_generic :
		269	fftw_hc2hc_backward_generic);
		270	int size, prev_size = 0, remaining_factors = n;
		271	fftw_plan r;
		272
		273	while (remaining_factors > 1) {
		274	size = fftw_factor(remaining_factors);
		275	remaining_factors /= size;
		276
		277	/* don't try the same factor more than once */
		278	if (size == prev_size)
		279	continue;
		280	prev_size = size;
		281
		282	/* Look for codelets corresponding to this factor. */
		283	{
		284	FOR_ALL_RCODELETS(p) {
		285	if (p->dir == dir && p->type == FFTW_HC2HC
		286	&& p->size == size) {
		287	size = 0;
		288	break;
		289	}
		290	}
		291	}
		292
		293	/*
		294	* only try a generic/rader codelet if there were no
		295	* twiddle codelets for this factor
		296	*/
		297	if (!size)
		298	continue;
		299
		300	r = rplanner(table, n / size, dir, flags,
		301	in, istride, out, ostride);
		302
		303	node = fftw_make_node_rgeneric(n, size, dir, codelet,
		304	r->root, flags);
		305	newplan = fftw_make_plan(n, dir, node, flags, FFTW_RGENERIC, 0);
		306	fftw_use_plan(newplan);
		307	run_plan_hooks(newplan);
		308	fftw_destroy_plan_internal(r);
		309	rcompute_cost(newplan, in, istride, out, ostride);
		310	best = fftw_pick_better(newplan, best);
		311	}
		312	}
		313	return best;
		314	}
		315
		316	static fftw_plan rplanner(fftw_plan *table, int n, fftw_direction dir,
		317	int flags, fftw_real *in, int istride,
		318	fftw_real *out, int ostride)
		319	{
		320	fftw_plan best = (fftw_plan) 0;
		321
		322	/* see if plan has already been computed */
		323	best = fftw_lookup(table, n, flags);
		324	if (best) {
		325	fftw_use_plan(best);
		326	return best;
		327	}
		328	/* try a wise plan */
		329	best = rplanner_wisdom(table, n, dir, flags,
		330	in, istride, out, ostride);
		331
		332	if (!best) {
		333	/* No wisdom. Plan normally. */
		334	best = rplanner_normal(table, n, dir, flags,
		335	in, istride, out, ostride);
		336	}
		337	if (best) {
		338	fftw_insert(table, best, n);
		339
		340	/* remember the wisdom */
		341	fftw_wisdom_add(n, flags, dir, RFFTW_WISDOM,
		342	istride, ostride,
		343	best->wisdom_type,
		344	best->wisdom_signature);
		345	}
		346	return best;
		347	}
		348
		349	fftw_plan rfftw_create_plan_specific(int n, fftw_direction dir, int flags,
		350	fftw_real *in, int istride,
		351	fftw_real *out, int ostride)
		352	{
		353	fftw_plan table;
		354	fftw_plan p1;
		355
		356	/* validate parameters */
		357	if (n <= 0)
		358	return (fftw_plan) 0;
		359
		360	if ((dir != FFTW_FORWARD) && (dir != FFTW_BACKWARD))
		361	return (fftw_plan) 0;
		362
		363	fftw_make_empty_table(&table);
		364	p1 = rplanner(&table, n, dir, flags,
		365	in, istride, out, ostride);
		366	fftw_destroy_table(&table);
		367
		368	if (p1)
		369	fftw_complete_twiddle(p1->root, n);
		370	return p1;
		371	}
		372
		373	fftw_plan rfftw_create_plan(int n, fftw_direction dir, int flags)
		374	{
		375	fftw_real *tmp_in;
		376	fftw_real *tmp_out;
		377	fftw_plan p;
		378
		379	if (flags & FFTW_MEASURE) {
		380	tmp_in = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		381	tmp_out = (fftw_real ) fftw_malloc(n sizeof(fftw_real));
		382	if (!tmp_in \|\| !tmp_out)
		383	return 0;
		384
		385	p = rfftw_create_plan_specific(n, dir, flags,
		386	tmp_in, 1, tmp_out, 1);
		387
		388	fftw_free(tmp_in);
		389	fftw_free(tmp_out);
		390	} else
		391	p = rfftw_create_plan_specific(n, dir, flags,
		392	(fftw_real ) 0, 1, (fftw_real ) 0, 1);
		393
		394	return p;
		395	}
		396
		397	void rfftw_destroy_plan(fftw_plan plan)
		398	{
		399	fftw_destroy_plan_internal(plan);
		400	}
		401	/*
		402	void rfftw_fprint_plan(FILE *f, fftw_plan p)
		403	{
		404	fftw_fprint_plan(f, p);
		405	}
		406
		407	void rfftw_print_plan(fftw_plan p)
		408	{
		409	rfftw_fprint_plan(stdout, p);
		410	}
		411	*/

Subversion Repositories shark

(root)/shark/branches/pj/ports/fftw/rfftw/rplanner.c - Rev 2